Modified gravity and massive neutrinos: constraints from the full shape analysis of BOSS galaxies and forecasts for Stage IV surveys

We constrain the growth index γ by performing a full-shape analysis of the power spectrum multipoles measured from the BOSS DR12 data. We adopt a theoretical model based on the Effective Field theory of the Large Scale Structure (EFTofLSS) and focus on two different cosmologies: γCDM and γνCDM, where we also vary the total neutrino mass. We explore different choices for the priors on the primordial amplitude As and spectral index ns , finding that informative priors are necessary to alleviate degeneracies between the parameters and avoid strong projection effects in the posterior distributions. Our tightest constraints are obtained with 3σ Planck priors on As and ns : we obtain γ = 0.647 ± 0.085 for γCDM and γ = 0.612+0.075 -0.090, Mν < 0.30 for γνCDM at 68% c.l., in both cases ∼ 1σ consistent with the ΛCDM prediction γ ≃ 0.55. Additionally, we produce forecasts for a Stage-IV spectroscopic galaxy survey, focusing on a DESI-like sample. We fit synthetic data-vectors for three different galaxy samples generated at three different redshift bins, both individually and jointly. Focusing on the constraining power of the Large Scale Structure alone, we find that forthcoming data can give an improvement of up to ∼ 85% in the measurement of γ with respect to the BOSS dataset when no CMB priors are imposed. On the other hand, we find the neutrino mass constraints to be only marginally better than the current ones, with future data able to put an upper limit of Mν < 0.27 eV. This result can be improved with the inclusion of Planck priors on the primordial parameters, which yield Mν < 0.18 eV.


Introduction
The current concordance model of cosmology, the ΛCDM model, has been confirmed over the last few decades by increasingly precise observations, spanning from the Cosmic Microwave Background (CMB) [1] to the Large Scale Structure (LSS) [2][3][4].In the standard framework, the gravitational interaction is described by General Relativity (GR) and the energy content of the Universe is composed of matter, which includes cold dark matter (CDM) and baryons, and dark energy in the form of a cosmological constant (Λ), introduced to explain the observed accelerated expansion of the Universe [5,6].Despite the many successes of the model at fitting observations, major fundamental questions concerning the nature of the dark components remain unanswered.Furthermore, the increasing accuracy of recent observations has brought to light tensions in the cosmological parameters when measured from early versus late time probes [7][8][9].This context offers the perfect breeding ground for alternative theories, which can explain the accelerated expansion of the Universe without the need for a dark energy component, or mitigate the tensions in cosmological parameter measurements.
Confirming or disproving the standard picture is the primary goal of ongoing and forthcoming Stage-IV galaxy redshift surveys, such as DESI [10], the Euclid satellite mission [11], Rubin's Legacy Survey of Space and Time [12] and the Nancy Grace Roman space telescope [13].These will map the 3D galaxy distribution with unprecedented accuracy over extremely large volumes, delivering high-precision measurements of the cosmological observables.The latter will allow to measure the cosmological parameters to sub-percent precision and disentangle between different gravity models [14,15].A crucial element to achieve this goal is the availability of an accurate and reliable theoretical model for the cosmic observables, able to describe the nonlinear regime of structure formation.
The standard approach to describe the impact of nonlinear evolution on clustering observables is based on Perturbation Theory (see [16] for a review), which features contributions to the power spectrum in the form of convolution integrals.If computed with standard integration techniques, such integrals are too computationally expensive to be evaluated over a large number of points in parameter space; for this reason, previous analyses have relied on a template-fitting approach.However, recent advances on the theoretical and numerical fronts have allowed for a full-shape analysis of the summary statistics.On the one hand, the development of the Effective Field Theory of the Large Scale Structure (EFTofLSS, [17][18][19]) has allowed to include the impact of unknown small-scale physics on the intermediate scales, thus extending the validity range of the model.On the other hand, the development of FFTlog-based algorithms [20][21][22] has sped up the computation of the convolution integrals by orders of magnitude, enabling sampling of the full parameter space.
Constraints on the parameters of the ΛCDM model obtained from full-shape fits of the clustering measurements of the Baryon Oscillation Spectroscopic Survey (BOSS) have been presented in a number of recent papers, see e.g.[23][24][25].In terms of constraints on beyond-ΛCDM scenarios from the same dataset, the BOSS collaboration adopted a templatefitting approach [26][27][28].Additionally, a full-shape analysis has been performed in [29][30][31][32][33][34] in the context of specific beyond-ΛCDM models.In this work, we focus instead on a model independent way of checking deviations from base ΛCDM, the so-called growth index parameterisation [35,36].This is a one-parameter extension of the standard model which allows the growth of structures to deviate from the prediction of ΛCDM by means of the growth index, γ.Despite the model being purely phenomenological, a measurement of γ inconsistent with its ΛCDM prediction would provide a clear detection of deviations from the standard cosmological model.On the other hand, it has been shown that the growth index parameterisation can provide a good fit for specific modified gravity models, such as the DGP model [37], with γ ≃ 0.67 [38,39], and a class of Horndeski theories [40], although in that case a time dependent growth index is required.For this reason a precise measurement of the growth index is among the objectives of forthcoming Stage-IV surveys [10,14,41].In addition to testing the theory of gravity, a key science goal of forthcoming galaxy surveys is a measurement of the neutrino mass.It is well known that massive neutrinos can suppress structure formation below a specific 'free streaming' scale, thus leaving a distinct signature on cosmological observables (see [42] for a review).However, current LSS data alone are not sufficiently accurate to constrain the neutrino mass, or even provide competitive upper bounds with respect to CMB measurements.This picture is expected to change with Stage-IV measurements of the LSS, that promise to achieve sufficient precision to allow for a detection of the neutrino mass.
In this work, we perform the full-shape analysis of the galaxy power spectrum as measured from BOSS Data Release 12, providing joint constraints for the growth index and the total neutrino mass.The paper is structured as follows: in Sec. 2 we give an overview of the theoretical model for the power spectrum, and describe the dataset and our analysis setup in Sec. 3. Our main results are presented in Sec. 4, where we also explore the impact on the constraints of different prior choices for the primordial parameters A s and n s .Additionally, we present forecasts on γ and M ν for Stage-IV spectroscopic surveys, focusing on a DESI-like galaxy sample.Finally, we summarise our conclusions in Sec. 5.

Theoretical model
We model the nonlinear galaxy power spectrum in redshift space using the Effective Field Theory of the Large Scale Structure (EFTofLSS, [17,18,43,44]), which allows to model the effects of unknown small-scale physics on mildly nonlinear scales with the introduction of scale-dependent terms in the theoretical model.In particular, we follow a prescription analogous to the one of [45], but employ an independent implementation [46][47][48], whose redshift space modelling for the power spectrum has been validated on N-body simulations in [49,50].The model features several ingredients, including loop corrections, EFTofLSS counterterms, a bias prescription and an infrared resummation routine; we outline these below, but refer to [32] for a more detailed description.

Nonlinear power spectrum
The bias expansion we adopt is based on a perturbative expansion of the galaxy density field δ g [51][52][53][54], which includes contributions from the underlying matter density field and large-scale tidal fields.Here we only consider terms up to third order in perturbations of δ: where b 1 and b 2 are respectively the linear and quadratic local bias parameters, G 2 and Γ 3 are non-local functions of the density field and ϵ is a stochastic contribution.The one loop anisotropic galaxy power spectrum can then be written as: where P L is the linear power spectrum 1 and ) are the redshiftspace kernels [57], whose full expressions for the bias basis of Eq. 2.1 can be found in Appendix A of [24].The integrals in Eq. 2.2 are computed under the assumption that time and scale evolution are fully separable, i.e. we adopt the so-called Einstein-de Sitter approximation.The latter was shown to be better than 1% accurate in ΛCDM at the redshift of BOSS [58], and allows us to compute the integrals only once per set of cosmological parameters, and then re-scale them using the linear growth factor to get the observables at the desired redshift.The power spectrum of Eq. 2.2 includes two additional contributions: the EFTofLSS counterterms, and noise.The EFTofLSS counterterms can be written as: where µ is the cosine of the angle between the wavevector k and the line of sight.Following [45], we re-define the k2 -counterterm parameters in order to have separate contributions to each multipole.We model the noise as: where N is a constant that includes deviations from pure Poisson shot noise, and we have two additional scale-dependent terms.We account for the nonlinear evolution of the BAO peak [59][60][61] by means of an infrared resummation routine, applied by splitting the linear power spectrum in a smooth broadband component (computed using the fit of [62]) and a wiggle part.Full details of the approach we use can be found in Sec.2.2.2 of [32].In total we have a set of 11 nuisance parameters per redshift bin: In our baseline analysis we keep the total neutrino mass fixed to its fiducial value M ν = 0.06 eV.We model the scale-dependent suppression induced by massive neutrinos by computing the cdm+baryon linear power spectrum and using it as input for the theoretical model (instead of the total matter one, which also includes the neutrino contribution).Additionally, we consider a cosmology with free neutrino mass, and we modify the theoretical model as described in more detail in Sec.2.3.
In order to include the impact of the fiducial cosmology assumed when converting redshifts to distances in the data, we correct wavenumbers and angles by applying Alcock-Paczynski distortions [63]: as well as re-scaling the power spectrum by a factor where H 0 is the Hubble factor today, D A is the angular diameter distance, and fid refers to quantities evaluated in the fiducial cosmology.Finally, we project the anisotropic power spectrum to multipoles in the usual way: where P ℓ (µ) is the Legendre polynomial of order ℓ.
Our implementation allows for bypassing the camb2 Boltzmann solver [64] in order to use linear power spectrum emulators, namely bacco3 [65] and CosmoPower4 [66].We note that this speeds up the model evaluations by roughly two orders of magnitude.

The growth index γ parameterisation
In order to check for deviations from ΛCDM, we adopt the phenomenological parameterisation proposed in [35,36].Specifically, we compute the growth rate f = d ln D/d ln a as: where Ω m (a) = Ω m,0 a −3 /E 2 (a) is the time dependent matter density parameter, with E(a) the dimensionless Hubble factor.Eq. 2.10 has been shown to be 0.2% accurate for ΛCDM scenarios, with γ = 0.55 [36].We compute the linear growth factor D by numerically integrating Eq. 2.10: We normalise the growth rate to its ΛCDM value at high redshift, i.e. during matter domination.The growth functions from Eq. 2.10 and 2.11 are then used to re-scale the linear power spectrum computed at redshift z = 0 and construct the redshift space multipoles as described in Sec.2.1.An analysis of the behaviour of γ over cosmic history can be found in [67,68], while [40] propose an alternative parameterisation, that extends γ to be a function of redshift.

Massive neutrinos
The scale-dependent suppression in growth induced by massive neutrinos free streaming (see [42] for a review) is modelled following the prescription of [55], which consists in using the cold dark matter + baryon power spectrum P cb instead of the total matter P m to compute the galaxy power spectrum multipoles.While not exact, this approach was shown to be sufficiently accurate for survey volumes similar to that of BOSS in [69], especially for small neutrino masses.For our baseline analysis with fixed M ν = 0.06 eV we modify the model described in Sec.2.1 by substituting δ with δ cb in Eq. 2.1, and thus P L with P cb in what follows.
For the case where the neutrino mass is varied as a free parameter, given the bigger impact that neutrino masses larger than the minimal one have on the power spectrum, we perform the following modifications to the model: • we compute the linear P cb twice, once for each effective redshift of the galaxy sample, in order to properly model the scale-dependent suppression introduced by massive neutrinos; • we modify our infrared resummation routine, since the Eisenstein-Hu prescription we adopt for the broadband power spectrum does not include a dependence on massive neutrinos.We use instead a discrete sine transform approach to perform the wiggle-no wiggle decomposition [70,71]; • we compute the growth rate of Eq. 2.10 using only the cdm+baryon density Ω cb [72]; • the resulting growth factor is normalised to its ΛCDM counterpart at the effective redshift, so that differences in the amplitude only come from the γ parameter.
In principle, a more accurate modelling of the impact of massive neutrinos would require a modification of the perturbation theory kernels and the use of a scale dependent growth rate (see e.g.[73][74][75][76][77]).However, we do not expect any relevant effect for the case of BOSS measurements given the size of the error bars [69].

Data and analysis setup
Our main analysis setup and dataset follow the ones described in [32].We summarise here the main points, but refer to that work for more details.

Dataset
We use the galaxy power spectrum multipoles of BOSS DR12 [78][79][80], which is split in two galaxy samples (CMASS and LOWZ) and covers two different sky cuts (NGC and SGC).The measurements are performed by splitting the sample into two redshift bins, with effective redshift z 1 = 0.38 and z 3 = 0.61 and widths ∆z = 0.3 and 0.25, respectively.In particular, we use the windowless measurements provided in [81], based on the windowless estimator of [82,83].The power spectrum measurements are complemented with BAO measurements of the Alcock-Paczynski parameters α ⊥ , α ∥ obtained from the same BOSS data release and also provided by [81], as well as pre-reconstruction BAO measurements at low redshift (6DF survey [84] and SDSS DR7 MGS [85]) and high redshift measurements of the Hubble factor and angular diameter distance from the Ly-α forest auto and cross-correlation with quasars from eBOSS DR12 [86].We use a numerical covariance matrix computed from 2048 'MultiDark-Patchy' mock catalogs [87,88], also provided in its windowless form by [81].

Analysis setup
We fit all three multipoles of the galaxy power spectrum up to k max = 0.2 h Mpc −1 .The nonlinear power spectrum is modelled with the EFTofLSS as described in Sec.2.1, with a custom implementation which takes advantage of the FAST-PT5 algorithm [20,21] to compute the convolution integrals of Eq. 2.2 in O(10 −2 ) seconds.The theoretical prediction is then combined with an independently developed likelihood pipeline, where we sample the parameter space by means of the affine-invariant sampler implemented in the emcee6 package [89], using the integrated auto-correlation time 7 as convergence diagnostic for our MCMC chains [90].For the Stage-IV forecasts presented in Sec.4.3 we use the pre-conditioned MonteCarlo method implemented in pocomc8 [91,92], which allows for a speed up in the sampling of parameters.We check the two samplers give consistent results.
We marginalise analytically over the nuisance parameters that enter the model linearly [23,30], namely the EFT counterterm parameters, the noise parameters, and b Γ3 .The parameter space is thus restricted to 18 free parameters (19 when we also vary the total neutrino mass): where the bias parameters assume different values in each redshift bin and sky cut (e.g. ).
Concerning priors 9 , we impose a Gaussian BBN prior on ω b10 , and broad, flat priors on the cosmological parameters, matching the ranges of parameters allowed by the bacco linear emulator: Ω cb ∈ U(0.06, 0.7) , h ∈ U(0.5, 0.9) , M ν ∈ U(0, 1) , where and Ω cb = Ω c + Ω b .We impose no prior on the other cosmological parameters in our baseline analysis, but we also explore the impact of including Planck priors on the primordial parameters A s and n s .In particular, for each cosmology explored (γCDM, γνCDM), we have three options: • one with no priors on A s , n s (labelled baseline); • one with a 3σ Gaussian Planck prior on A s : ln(10 10 A s ) ∈ N (3.044, 0.042) (labelled prior A s ); • and one with a 3σ Gaussian Planck prior on both A s and n s : ln(10 10 A s ) ∈ N (3.044, 0.042) and n s ∈ N (0.9649, 0.0126) (labelled prior A s ,n s ).
For both parameters, the Gaussian priors are centered on the best-fit values from [1] and have width corresponding to the 3σ error from the same work.For the nuisance parameters, we adopt those of [81]:

BOSS analysis
Here we present the main results of the analysis of the galaxy power spectrum multipoles measured from the BOSS DR12; we cover the two extensions of the base ΛCDM scenario discussed in the previous sections: γCDM and γνCDM.For each cosmology, we explore the three options for priors described in Sec.3.2.We show the posterior distributions for the γCDM case with fixed neutrino mass M ν = 0.06 eV in Fig. 1 and summarise the 68% c.l., mean and best-fit values for the cosmological parameters in Table 1.Plots for the full parameter space can be found in Appendix A. We obtain γ = 0.007 +0.170  −0.229 (0.149) for the baseline priors, shown by green contours, γ = 0.617 +0.098 −0.110 (0.614) (prior on A s , orange contours) and γ = 0.647 ± 0.085 (0.655) (prior on A s and n s , purple contours), where values in parentheses refer to the best-fit values, obtained as the maximum of the analytically marginalised posterior.
We first note that our baseline analysis is affected by strong projection effects (also known as prior volume effects), in particular involving the parameters regulating the amplitude of the power spectrum, A s and γ: when no CMB priors are imposed, the two-dimensional marginalised posteriors are shifted towards extremely low values of A s and γ, along the degeneracy between the two.In fact, lowering the primordial amplitude A s or introducing a low value for γ, which effectively increases the growth factor, result in two opposite effects which can balance out and give the same power spectrum.An analogous behaviour for extensions of the standard model was highlighted in [32], especially for the case of exotic (interacting) dark energy cosmology (see their Fig. 7).Other works focusing on the full-shape analysis of the same dataset also highlighted issues when trying to constrain extended parameter models [34,96].In particular in [34] the authors provide constraints on nDGP modified gravity [37] using a similar EFTofLSS-based model, and discuss the presence of projections due to the degeneracy between the primordial amplitude A s and the nDGP Ω rc parameter.We study the presence of prior volume effects in more detail in Sec.4.2 by generating and fitting synthetic data, adopting the same covariance matrix used for our main BOSS analysis, and by profiling the posterior.
Imposing a prior on A s breaks these degeneracies and shifts the two-dimensional marginalised posterior closer to the Planck best-fit and ΛCDM values, as can be seen in the orange contours of Fig. 1.In fact, the other cosmological parameters are mostly unaffected, but the error on γ is reduced by ∼ 48%.Moreover, adopting a prior on n s yields an additional ∼ 18% improvement on γ.Both the mean and best-fit values obtained in these cases are slightly larger than the ΛCDM prediction of γ = 0.55, though still ∼ 1σ consistent.This slight deviation is likely equivalent to the finding in previous full-shape BOSS analyses of a low value of A s11 [23,24,32] (ln(10 10 A s ) ≃ 2.8 as opposed to the Planck value ln(10 10 A s ) = 3.044): when we impose a prior on A s the model tries to compensate by lowering the amplitude with a larger value for γ.We also notice that the degeneracy with A s was masking a degeneracy between γ, ω c and n s , which emerges when we impose a prior on A s and is then broken when we apply a Planck prior on n s .We then open up the parameter space and allow the neutrino mass to vary as a free parameter.It is worth noting that we do not expect to obtain tight constraints on the neutrino mass, as previous analyses already highlighted the inability of BOSS data alone to do so [23,31,98]; in fact, as discussed in [23] the precision of BOSS data does not allow to detect the step in the power spectrum which marks the free streaming length of massive neutrinos, restricting detectable signatures only to changes to the overall amplitude.Nevertheless, our goal is to study possible degeneracies between γ and M ν , and the impact of massive neutrinos on the constraints on γ.We modify the theoretical model as described in Sec.2.3, and explore the same three options for priors as for the γCDM case.The marginalised posteriors for γνCDM are shown in Fig. 2, and the 68% c.l., mean and best-fit values for the cosmological parameters are listed in Table 2.
We first note that the baseline shows a similar shift in the two-dimensional marginalised posteriors for γ and A s as the one of Fig. 1.This is again due to the presence of strong degeneracies between the parameters that control the amplitude of the power spectrum.However, the error on γ is not affected by the introduction of M ν , and in fact there seems to be no degeneracy between the two.This is actually due to the fact that such degeneracy is mild compared to the others and to the size of the error bars on BOSS data (see e.g.Fig. 12 where we perform forecasts for Stage-IV surveys, an anti-correlation between γ and M ν can clearly be seen).
When we impose Planck priors on the primordial parameters on the other hand we can see a mild degeneracy emerge between γ and M ν , which leads to marginally lower values for γ with respect to γCDM.Concerning neutrinos, as expected the data are not very constraining, and we only provide (somewhat large) upper limits.It is worth commenting on the case with a CMB prior on A s , where the model seems to be picking up a large value for M ν : we argue that this is to be ascribed to the model trying to lower the amplitude by exploiting the degeneracy between M ν and n s .More specifically, in the baseline case the amplitude is mostly controlled by A s ; when that is fixed by the prior, other parameters change to compensate: mostly γ, but also the neutrino mass, resulting in a looser constraint on M ν in that case.However, this is only possible if the change in the shape of the power spectrum due to massive neutrinos is compensated by a change in n s .When n s is also fixed by the Planck prior, there is less freedom in the shape and thus the constraint on M ν tightens again.In Fig. 3 we plot the one-dimensional marginalised posterior distributions for γ for all cases described in this section, to allow for an easier comparison between the constraints we get for this parameter.We mark the best-fit values with dashed lines, following the same color scheme of Fig. 1 and 2, and the ΛCDM prediction with solid grey lines.As we can see, adopting CMB priors on the primordial parameters shifts the peak of the posterior from (projection affected) negative values to values consistent with the ΛCDM prediction, albeit slightly larger.Moreover, the difference between the peak of the posterior and the best-fit for the case without CMB-based priors is a hint of the projection which affects the marginalised posteriors.
Our results are broadly consistent with previous analyses of the BOSS data, although an exact comparison is difficult due to different choices in the analysis setup and datasets included.The main novelty of this work is the use of a full-shape approach, as opposed to template fitting adopted in the official BOSS analysis, and the exploration of the impact of priors on the inferred parameters.Moreover, most previous works take advantage of a joint analysis with CMB data, which results in generally tighter constraints.Nevertheless, we report here some of the recent measurements of γ that used BOSS galaxies as the main dataset for comparison, highlighting the main differences with respect to this work: [28] found γ = 0.566 ± 0.058 combining BOSS DR12 galaxies with Planck CMB data and SNIa data.Additionally, [27] provides joint constraints on γ and the neutrino mass M ν , although they used a previous data release (DR11) and again performed a joint fit with CMB data: they found γ = 0.67 ± 0.14, M ν = 0.25 +0.13 −0.22 .A joint analysis with the galaxy bispectrum was performed in [99], finding γ = 0.733 +0.068 −0.069 .All these works are official analyses from the BOSS collaboration, and adopt a fixed template for the power spectrum multipoles.More recently, [100] performed a similar analysis combining data from BOSS, DES and Planck, finding γ = 0.633 +0.025 −0.024 .Also in this case the power spectrum was fitted using a fixed template, as opposed to the full-shape fit we perform in this work.
Overall, our strongest constraints (those where we impose CMB-based priors on the primordial parameters) show a marginal preference for a high value of γ, consistent for example with [100].However, in our case the deviation from the ΛCDM prediction is ∼ 1σ, as opposed to the ∼ 4σ discrepancy found in that work.In general, our constraints are weaker than those obtained in previous works.This can be attributed to the fact that we only rely on LSS data and include CMB information merely in terms of priors on the primordial parameters: given the degeneracies between the additional parameters of γCDM and γνCDM and the other parameters, some of which are tightly constrained by the CMB, it is no surprise that performing a joint analysis with Planck data can yield tighter constraints.This can also be approximated by including Gaussian priors derived from Planck on all cosmological parameters, see Appendix F of [101] for a comparison.
Another possible way of gaining in constraining power without resorting to CMB information would be to use higher-order statistics such as the galaxy bispectrum.Indeed, the tree-level bispectrum does not depend on b Γ 3 , and its inclusion in the analysis can therefore break the degeneracy between b G 2 and b Γ 3 , leading to improved constraints on the linear bias b 1 (see [23,81] for a full-shape joint analysis of BOSS data).While the impact on cosmological parameters is somewhat marginal in the context of ΛCDM and Stage-III surveys, it can become significant for extended models where there are strong degeneracies between parameters controlling the amplitude: a tighter constraint on b 1 implies a more precise mea-surement of A s and therefore of γ.Some work towards quantifying the impact of including the bispectrum for beyond-ΛCDM models has already been done in [50] in the context of simulations, where a ∼ 30% improvement was found for a power spectrum and bispectrum joint analysis for an interacting dark energy model.We leave a proper exploration of the impact of the bispectrum in the context of γCDM and γνCDM to a future work.

Projection effects
In this section we study in more details the prior volume effects found for γCDM and γνCDM when no CMB priors are applied (green contours in Fig. 1 and 2, in particular the A s -γ plane).The presence of projection effects is already highlighted by the shift between the peak of the one-dimensional marginalised posterior and the best-fit point, see Fig. 3. To better investigate this we adopt two approaches: first we perform a consistency check using synthetic data, then we carry out profiling of the posterior.
For the consistency check we perform the same analysis described in Sec.4.1 on a synthetic datavector generated with a fiducial cosmology and set of nuisance parameters.We note that the ability of our pipeline to recover the cosmological parameters correctly has been validated with N-body simulations, and we are using the same theoretical model to generate the datavector and fit it, which should lead to perfect agreement between the input and inferred parameters.We use the same covariance computed from the 'MultiDark-Patchy' mock catalogues adopted in the BOSS analysis, and include (synthetic) BAO information.As fiducial parameters we adopt the best-fit values of Planck [1], while the nuisance parameters are determined by maximising the BOSS likelihood with cosmology fixed to the Planck values.Our results are plotted in Fig. 4 (purple contours and lines), where we also include the baseline result of Fig. 1 for comparison (green contours and lines).
The two-dimensional marginalised posteriors for the synthetic data exhibit a very similar behaviour as the one we observed in the BOSS data, with a strong shift in the posterior towards low values of A s and γ and almost perfect overlap with the posterior obtained from BOSS data.As discussed above, this effect can be ascribed to the presence of extreme degeneracies in the extended parameter space, which lead to a highly non-Gaussian posterior.As a consequence, the marginalisation process can result in apparent biases in the two-dimensional confidence regions, as discussed thoroughly in [102]. 12In that work, the author proposes to use profile distributions to assess the impact of marginalisation.As opposed to marginalisation, this profile likelihood (PL) approach allows one to obtain a one-dimensional distribution by evaluating the posterior over a scan of a parameter of interest, θ i , while maximising the posterior with respect to all other parameters θ k , k ̸ = i.This profile distribution is not biased by projection effects, since it is centred on the maximum of the full posterior.For non-Gaussian posteriors, as is the case here, the resulting distribution will be different than the usual marginalised distribution obtained via integration over the marginalised parameters.Depending on the type of non-Gaussianity, this distribution can be broader or narrower over the parameter of interest.This difference between integration and maximisation of the posterior is approximately given by the so-called Laplace term, which accounts for volume effects as explained in Ref. [103].
Ref. [102] suggests performing this profiling step directly on the MCMC samples.However, in addition to resulting in noisy distributions, this is not ideal in our case, since we sample over the analytically marginalised posterior and we therefore expect additional projection due to that pre-marginalisation step.We profile instead the full posterior using the MINOS algorithm implemented in the iminuit package13 [104].Our parameter of interest is the growth index γ, for which we plot the PL in the left panel of Fig. 5, and compare to the one-dimensional marginalised posterior from the MCMC in the right panel.Additionally, we compare the confidence intervals obtained from the MCMC to those obtained from the profiling in Fig. 6.Fig. 5 shows a clear shift of the PL towards the ΛCDM prediction with respect to the marginalised posterior.Moreover, the confidence intervals are larger in the case of the PL, making the results consistent with ΛCDM.We also note how the best-fit points from the MCMC (orange crosses in Fig. 6) are close to the PL results, though still somewhat lower in A s .This demonstrates the presence of additional projection effects due to the analytic marginalisation of linear nuisance parameters in the posterior sampled in the MCMC.Overall, the use of the profile likelihood shows a more conservative result, which appears, at least in this case, to reflect better the effect of the large degeneracy between the amplitude parameters on the final one-dimensional constraint on γ.Note that should the posterior be closer to the Gaussian case, the difference between marginalising and profiling would disappear, so we expect projection effects to disappear with more constraining data.
As mentioned already in Sec 4, the presence of volume (or projection) effects was already highlighted in a number of recent works performing the full-shape analysis of BOSS data in the context of beyond-ΛCDM models [32,34,96,105].It is expected that the higher precision of forthcoming data will mitigate this effect, as also suggested by the forecasts we present in Sec.4.3.However, projections can also be alleviated by the combination with additional probes or the inclusion of higher order correlation functions.In general, we recommend extreme care in the choice of priors, especially when constraining beyond-ΛCDM models.

Forecasts for Stage IV surveys
In this section we present forecasts for a Stage-IV spectroscopic galaxy survey, focusing in particular on DESI-like galaxies.Similarly to what described in Sec.4.2, we generate and fit synthetic datavectors with the same EFTofLSS theoretical model.The aim is twofold: on one hand, we wish to provide forecasted constraints on γ and jointly on γ and M ν , on the other hand, we can assess the ability of the higher precision measurement of reducing the projection effects.Despite a number of approximations we make in our forecasting strategy, which we highlight below, we expect the results to give a more realistic picture of the outcome of future surveys, as opposed to the standard Fisher matrix approach.The latter, while extremely useful due to its low computational cost, has been shown to be unable to capture complex features in the likelihood, which can arise in the context of strong degeneracies between parameters and highly non-Gaussian posteriors (see e.g.[106][107][108][109]).
We generate three sets of synthetic datavectors, assuming Planck values for the cosmological parameters and simulating the three galaxy samples that are the target of DESI: the bright galaxy sample (BGS), luminous red galaxies (LRGs), and star-forming emission line galaxies (ELG).We follow [10] and compute expected values for the linear bias and number density of the samples, which we list in Table 3 together   Table 3. Parameters used to generate the synthetic DESI-like datavectors for the forecasts.We follow [10] in computing the effective redshift z eff , the volume V s and number density n, and determine the linear bias b 1 assuming constant, sample-specific values for the combination b 1 (z)D(z) as specified in the text.

Galaxy sample z
All redshift bins have size ∆z = 0.4 and are centered on the z eff values listed in the table.We compute an analytic covariance matrix, assuming no cross correlation between the redshift bins / samples.Specifically, we use the following [110]: where µ i is the cosine between k i and the line of sight, L ℓ is the Legendre polynomial of order ℓ, n is the number density and N k = 4πk 2 ∆k V s /(2π) 3 is the number of Fourier modes within a given k-bin of size ∆k, V s being the survey volume.For the power spectrum P (k) in Eq. 4.1 we use the Kaiser approximation [111]: As for priors, we use the same ones we adopted for the baseline BOSS analysis of Sec.3.2, specifically we impose no prior on the primordial parameters A s and n s , both to showcase the constraining power of LSS alone and to assess potential projection effects in the context of Stage-IV surveys.The resulting power spectrum multipoles for the three samples, with corresponding error bars, are shown in Fig. 7.We fit the three galaxy samples separately and then combine them, and explore two options for the scale cuts: a pessimistic one with k max = 0.15 h Mpc −1 and an optimistic one with k max = 0.25 h Mpc −1 .Our results for γCDM are shown in Fig. 9 (8) for the optimistic (pessimistic) case, and summarised in Table 4.Moreover, we perform the same analysis for γνCDM and show the results in Fig. 11 (10) for the optimistic (pessimistic) case, while the constraints are listed in Table 5.
For the joint analysis in γCDM we find σ(γ) = 0.058 (0.072) for the optimistic (pessimistic) case, which represents an improvement of ∼ 85% with respect to the BOSS results when no CMB priors are imposed.For γνCDM we find the same values, despite the introduction of M ν as a free parameter.As for massive neutrinos, we get M ν < 0.27 (< 0.314) at 68% c.l. for the optimistic (pessimistic) case, with a very modest improvement with respect to the constraints obtained from BOSS data.Such a dramatic improvement in the constraints on γ is due to the combination of the different samples: each sample features a slightly different orientation of the degeneracy between A s and γ, which is then broken when we perform the joint fit.
Focusing on the projection effects which heavily affected the BOSS analysis, we can clearly see how they are strongly mitigated by the increased precision of the mock data and the combination of multiple samples at different redshifts.In fact, for the joint analysis of the three galaxy samples the deviation between the input and best-fit values for γ is ∼ 0.1σ, with the exception of the γνCDM cosmology in the optimistic case where the deviation is ∼ 0.4σ.On the other hand, some projection is still present in the fits for the single samples, although we can always recover the input value at the 1σ level.It is also worth noting how the inclusion of more nonlinear modes can further alleviate the projection.For example, from the BGS sample we obtain γ = 0.43 +0.25  −0.22 (0.37) in the pessimistic case and γ = 0.68 +0.22 −0.19 (0.521) in the optimistic case, with the shift between the best-fit point and the input value reduced from 0.37σ to 0.06σ.
A proper comparison to other forecasts in the literature is complicated, because the results depend strongly on the forecasting strategy.Nonetheless, we try to summarise previous results and highlight the differences with this work.Overall we find weaker constraints on γ, which can mainly be attributed to the following reasons: • most works rely on the Fisher matrix formalism, known to be unable to properly describe the likelihood in presence of strong degeneracies in the parameter space; • most works use a linear model for the power spectrum, with a significantly smaller number of nuisance parameters; • most works perform a joint analysis with CMB information, which, as discussed above, can tightly constrain some of the parameters which are most degenerate with γ.Forecasts for the γ parameterisation for Stage-IV surveys can be found in [112].The authors use a Fisher matrix approach and a linear model for the power spectrum, limiting their analysis at k max < 0.1 h Mpc −1 , and find an overly optimistic σ(γ) = 0.0096 for the case with fixed neutrino mass and dark energy parameter.This result is mainly driven by the differences in the dimensionality of the parameter space considered: their modelling only involves six parameters (nine when they also consider massive neutrinos and evolving dark energy), with only one (linear) bias, and no parameter controlling the overall amplitude of the power spectrum, whose inclusion can undermine the ability to measure γ.
A more comprehensive collection of forecasts for future experiments can be found in [113], that however still relies on a Fisher matrix approach.All forecasts presented there include Planck CMB and are limited to a linear model for the power spectrum, and can therefore be expected to be more constraining than our analysis.Nonetheless, the authors find σ(γ) = 0.025, which is closer to our forecasted σ(γ) = 0.058 than the results of [112].We note that lifting some of the simplifications we make when generating the synthetic dataset could improve our constraints.For example, considering several redshift bins per sample can yield more precise measurements of the cosmological parameters, including γ [114,115], but it also leads to an increased number of nuisance parameters and increased shot-noise per redshift bin.Moreover, some of the effects we are not accounting for are expected to result in reduced constraining power, for example the assumption of a Gaussian covariance, neglecting the convolution with the survey window function and in general additional observational systematics which can increase the size of the error bars on the measurements [116,117].We finally comment on our forecasts on massive neutrinos.A measurement of the neutrino mass is one of the main goals of Stage-IV surveys, with predicted errors well below the detection threshold (e.g.[101,113,118]).In this sense, our forecasts might seem too pessimistic, given our best case scenario gives M ν < 0.27 eV at 68% c.l. Once again, most previous forecasts take advantage of a combined analysis with CMB data from Planck, able to put tight constraints on the primordial parameters A s and n s .The latter are degenerate with the neutrino mass M ν , as can be seen in Fig. 11 (or equivalently and more clearly in the green contours of Fig. 12): we therefore expect a joint analysis to be extremely beneficial.We do not perform such an analysis here, but we can get a sense of the impact by imposing a 3σ Planck prior on the primordial parameters, as done in Sec.4.1 for the BOSS analysis.We only focus on the combination of the three redshift bins for the optimistic case with   k max = 0.25h Mpc −1 .Our results are shown in Fig. 12 and Table 6.As expected, the CMB prior on A s and n s breaks the degeneracies and brings the upper limit on the neutrino mass down to M ν < 0.175 eV at 68% c.l., with a 35% improvement over the case with no priors.Despite being quite large in the context of a detection of the neutrino mass, our results are generally compatible with [69]; there, the authors performed a validation against Nbody simulations of a similar 1-loop model for the power spectrum as the one adopted here.For a 25 Gpc 3 /h 3 survey, similar in volume to our combined DESI-like mocks, they found M ν < 0.49 eV at 68% c.l. and conclude that 2-point summary statistics from spectroscopic clustering alone might not be sufficient detect the neutrino mass.Beside the combination with CMB data, we expect the combination with higher order statistics such as the bispectrum to bring significant improvements [119,120].We leave the study of the impact of the bispectrum on the determination of the neutrino mass scale to a future work.Further improvement should come from using photometric probes, namely cosmic shear and galaxy-galaxy lensing [41].

Conclusions
We presented the full-shape analysis of the power spectrum multipoles measured from BOSS DR12 galaxies with the inclusion of post-reconstruction BAO data.We used the windowless measurements of [81] combined with BAO measurements obtained from a variety of data-sets [81,[84][85][86].We provided constraints on a single parameter extension of ΛCDM which allows for deviations from the standard scenario in the growth functions, the so-called 'growth index' (or γ) parameterisation.We also explored the case with the total neutrino mass as a free parameter, and provided joint constraints for γ and M ν .Our theoretical model for the power spectrum is based on the EFTofLSS and takes advantage of the bacco emulator for the linear power spectrum and the FAST-PT algorithm for a fast evaluation of the likelihood.We explored different options for the priors on the parameters that determine the primordial power spectrum, A s and n s , finding that strong degeneracies in the extended parameter space result in large projection effects, especially concerning the parameters which regulate the amplitude of the power spectrum.For this reason, we applied a 3σ Planck prior on A s and on both A s and n s , and obtained γ = 0.647 ± 0.085 and γ = 0.612 +0.075 −0.090 , respectively, ∼ 1σ consistent with the ΛCDM prediction γ = 0.55.For neutrinos we found M ν < 0.478 eV (Planck prior on A s ) and M ν < 0.298 eV (Planck prior on A s , n s ), consistent with previous EFTofLSS-based studies of the BOSS dataset that did not perform a full joint analysis with CMB data [23,98].
To assess the presence of projection effects in the case when no priors are imposed on the primordial parameters we generated synthetic datavectors with a fiducial set of cosmological and nuisance parameters, and then fitted them using the same numerical covariance used for the BOSS analysis.We found a similar shift of the posterior in the A s -γ plane as the one found in our baseline BOSS analysis.Additionally, we performed a profile likelihood analysis, finding the maximum of the PL to be closer to the ΛCDM prediction than the peak of the marginalised posterior, and the confidence intervals derived from the PL to be larger and consistent with γ = 0.55 at 68% c.l.. Overall, we find the EFTofLSS model complemented with Planck priors to provide constraints on γ that are only ∼ 30% larger than the official BOSS analysis [28], despite the fact that we did not perform a full joint analysis with CMB data.Additionally, our theoretical model features a significantly larger number of parameters, which is expected to degrade the constraints to some extent.On the other hand, the full-shape approach allows to provide direct constraints on the cosmological parameters.
In the second part of this work we presented forecasts for a Stage-IV spectroscopic survey, focusing on a DESI-like galaxy sample.We generated three synthetic datavectors at three redshifts using different values for the nuisance parameters, to match the expected BGS, LRGs and ELGs samples that are the target of DESI.We computed Gaussian covariance matrices neglecting the cross-correlation between samples, and adopted two different scale-cuts: k max = 0.15 h Mpc −1 (pessimistic) and k max = 0.25 h Mpc −1 (optimistic).We performed separate fits for each sample and then combined them in a joint analysis, finding the combination to provide significantly tighter constraints, yielding σ(γ) = 0.058 (0.072) in the optimistic (pessimistic) case, with a ∼ 85% improvement with respect to our baseline BOSS analysis without CMB-based priors.Concerning neutrinos, we found that the improvement with respect to Stage-III constraints is only marginal, M ν < 0.27 eV (M ν < 0.314 eV) in the optimistic (pessimistic) case at 68% c.l..
In order to reduce the error bars on γ and M ν , but also to keep the projection effects under control, we advocate the combination with additional observables.For example, a joint analysis with Planck data (or the inclusion of CMB-based priors, at least on the primordial parameters), can have a significant impact on the constraints.We explored this by performing a fit of the synthetic DESI-like data where we adopted 3σ Planck priors on A s and n s .We obtained a ∼ 20% improvement in the measurement of γ, and ∼ 35% for M ν with respect to the case with no priors.Similarly, a combination with weak lensing can alleviate the degeneracies and yield tighter constraints [41].However, an optimal exploitation of galaxy clustering data alone can already go in this direction: as shown in [50], the ability of the bispectrum to better determine the bias parameters can prove crucial in the context of extended models.We leave the study of the impact of the bispectrum for γCDM and γνCDM to a future work.

A Full contourplots
We show here the two-dimensional marginalised posteriors for the full parameter space sampled in the BOSS analysis presented in Sec.4.1.The results for γCDM are shown in Fig. 13, corresponding to the constraints on cosmological parameters shown in Fig. 1, while γνCDM .Two-dimensional marginalised posterior for all parameters for the BOSS analysis for γCDM and k max = 0.2 h Mpc −1 (see Fig. 1).We use the same color scheme as the main text to show the three prior choices.Dashed grey lines mark the Planck best-fit values (ΛCDM prediction for γ).
results are plotted in Fig. 14, corresponding to the constraints on cosmolofical parameters shown in Fig. 2. The additional model parameters not shown are the ones we analytically marginalise over, namely b Γ 3 , the EFT counterterm parameters and the noise parameters.We follow the same color scheme as the main text for the three prior choices, and mark Planck best-fit values with dashed grey lines.

B Massive neutrinos with fixed γ
We perform a complementary analysis of the BOSS data with free neutrino mass but with the γ parameter fixed to its ΛCDM value, γ = 0.55.The two-dimensional marginalised posteriors are shown in Fig. 15, and the best-fits are summarised in Table 7.We impose a wide, flat prior on the neutrino mass U(0, 1) eV, covering the parameter space allowed by the emulator for the linear power spectrum.We explore again the three options for the priors on the primordial parameters, and obtain the following constraints for the total neutrino mass: M ν < 0.319, best-fit 0.137 (baseline), M ν = 0.34 +0.13 −0.28 , best-fit 0.141 (prior on A s ), M ν = 0.29 +0.12 −0.24 , best-fit 0.111 (prior on A s and n s ), where values are expressed in eV.Additionally, for the baseline analysis we find a higher value of the amplitude: log(10 10 A s ) = 2.90 +0.17 −0.21 , as opposed to log(10 10 A s ) = 2.67 ± 0.17 found in the ΛCDM analysis of [32] (that used the same pipeline and priors adopted here).Comparing to previous constraints on the total neutrino mass from the full-shape analysis of BOSS data, we find our results to be in excellent agreement with [98].There, the authors perform a similar EFT-based analysis of ΛCDM + massive neutrinos, also finding that having the neutrino mass as a free parameter allows for a slightly larger value of A s , closer to the Planck best-fit value, which is then balanced by a peak in the posterior at M ν = 0.29 +0.31  −0.18 eV.The same EFTofLSS model adopted in this work was also used in [121] to constrain the total neutrino mass from BOSS data.However, in that work the fit is performed jointly with CMB data from Planck, which results in significantly tighter constraints: M ν < 0.16 eV at 95% c.l..  7. Best-fit and 68% c.l. values for the analysis with free neutrino mass and γ = 0.55 for the three prior choices.We also derive best-fit values for σ 8 .

Figure 1 .
Figure1.Marginalised posterior distribution for the cosmological parameters for the γCDM cosmology and the three prior choices, as detailed in the legend.We fit all three multipoles and use k max = 0.2 h Mpc −1 .Grey dashed lines mark the Planck best-fit values (ΛCDM prediction for γ).

Figure 2 .
Figure 2. Marginalised posterior distribution for the cosmological parameters for γνCDM for the three prior choices.We use k max = 0.2 h Mpc −1 and a flat prior M ν ∼ U(0, 1) eV for the neutrino mass.Grey dashed lines mark the Planck best-fit values (with the ΛCDM prediction for γ and minimal neutrino mass M ν = 0.06 eV).

Figure 3 .
Figure 3. Marginalised one-dimensional posterior distribution for the γ parameter for γCDM (left) and γνCDM (right), for the three prior choices as stated in the legend.Dashed vertical lines mark the best-fit values corresponding to the maximum of the marginalised posterior distribution; the solid grey line corresponds to the ΛCDM prediction.

Figure 4 .
Figure 4. Posterior distribution for the cosmological parameters for the BOSS data (green) and synthetic data generated with the Planck fiducial cosmology (purple).We fix the total neutrino mass to M ν = 0.06 eV and use a scale-cut of k max = 0.2 h Mpc −1 .Grey dashed lines mark the fiducial cosmology used to generate the synthetic data.

Figure 5 .Figure 6 .
Figure 5. Left: ∆χ 2 min for a set of values for γ obtained from profiling the likelihood with MINOS.The horizontal dotted line marks ∆χ 2 min = 1, which is used to determine confidence intervals.Right: comparison between the marginalised posterior obtained from the mcmc (green curve) and the PL (orange curve and datapoints).The dashed vertical line marks the ΛCDM prediciton.
with the respective effective volume.Specifically, we assume constant values for the combination b 1 (z)D(z), where D(z) is the linear growth factor evaluated at the central redshift of each bin, and use b LRG (z)D(z) = 1.7, b ELG (z)D(z) = 0.84, and b BGS (z)D(z) = 1.34.

4 Figure 7 .
Figure 7. Power spectrum multipoles for the three galaxy samples in our synthetic DESI-like dataset: BGS (green, triangle markers), LRGs (orange, round markers) and ELGs (blue, crosses).The error bars are the square roots of the diagonal elements of the analytic covariance matrix.

Figure 8 .
Figure 8. Marginalised posterior for γCDM for the pessimistic case with k max = 0.15 h Mpc −1 .We fit the three samples separately and jointly, as detailed in the legend.Dashed grey lines mark the fiducial values.

Figure 9 .
Figure 9. Marginalised posterior for γCDM for the optimistic case with k max = 0.25 h Mpc −1 .We fit the three samples separately and jointly, as detailed in the legend.Dashed grey lines mark the fiducial values.

Figure 10 .
Figure 10.Marginalised posterior for γνCDM for the pessimistic case with k max = 0.15 h Mpc −1 .We fit the three samples separately and jointly, as detailed in the legend.Dashed grey lines mark the fiducial values.

Figure 11 .
Figure 11.Marginalised posterior for γνCDM for the optimistic case with k max = 0.25 h Mpc −1 .We fit the three samples separately and jointly, as detailed in the legend.Dashed grey lines mark the fiducial values.

Figure 12 .
Figure 12.Comparison between the forecast when no CMB prior is imposed (green lines and contours) and the case with 3σ Planck priors on A s and n s (purple lines and contours).We use k max = 0.25 h Mpc −1 .

Figure 13
Figure 13.Two-dimensional marginalised posterior for all parameters for the BOSS analysis for γCDM and k max = 0.2 h Mpc −1 (see Fig.1).We use the same color scheme as the main text to show the three prior choices.Dashed grey lines mark the Planck best-fit values (ΛCDM prediction for γ).

Figure 14 .
Figure 14.Two-dimensional marginalised posterior for all parameters for the BOSS analysis with γνCDM and k max = 0.2 h Mpc −1 (see Fig.2).We use the same color scheme as the main text to show the three prior choices.Dashed grey lines mark the Planck best-fit values (with the ΛCDM prediction for γ and minimal neutrino mass M ν = 0.06 eV).

Figure 15 .
Figure 15.Marginalised posterior distribution for the cosmological parameters for ΛCDM with massive neutrinos and the three prior choices, as detailed in the legend.We fit all three multipoles and use k max = 0.2 h Mpc −1 .Grey dashed lines mark the Planck best-fit values.

Table 1 .
Mean values and 68% c.l. values for γCDM with fixed neutrino mass M ν = 0.06 eV for the three prior choices.We show the best-fit values in parentheses, and include derived constraints on σ 8 .

Table 2 .
Mean values and 68% c.l. values for γνCDM for the three prior choices.We show the best-fit values in parentheses, and include derived constraints on σ 8 .

Table 4 .
Mean values and 68% c.l. values for γCDM with fixed neutrino mass M ν = 0.06 eV for the three DESI-like galaxy samples fitted separately and jointly.We show the best-fit values in parentheses, and include derived constaints for σ 8 .

Table 5 .
Mean values and 68% c.l. values for γνCDM for the three DESI-like galaxy samples fitted separately and jointly.We show the best-fit values in parentheses, and include derived constraints for σ 8 .

Table 6 .
Best-fit and 68% c.l. values for the DESI forecasts with 3σ Planck priors on A s ,n s .We use k max = 0.25 h Mpc −1 .