Neutrino masses and cosmological parameters from a Euclid-like survey: Markov Chain Monte Carlo forecasts including theoretical errors

We present forecasts for the accuracy of determining the parameters of a minimal cosmological model and the total neutrino mass based on combined mock data for a future Euclid-like galaxy survey and Planck. We consider two different galaxy surveys: a spectroscopic redshift survey and a cosmic shear survey. We make use of the Monte Carlo Markov Chains (MCMC) technique and assume two sets of theoretical errors. The first error is meant to account for uncertainties in the modelling of the effect of neutrinos on the non-linear galaxy power spectrum and we assume this error to be fully correlated in Fourier space. The second error is meant to parametrize the overall residual uncertainties in modelling the non-linear galaxy power spectrum at small scales, and is conservatively assumed to be uncorrelated and to increase with the ratio of a given scale to the scale of non-linearity. It hence increases with wavenumber and decreases with redshift. With these two assumptions for the errors and assuming further conservatively that the uncorrelated error rises above 2% at k = 0.4 h/Mpc and z = 0.5, we find that a future Euclid-like cosmic shear/galaxy survey achieves a 1-sigma error on Mnu close to 32 meV/25 meV, sufficient for detecting the total neutrino mass with good significance. If the residual uncorrelated errors indeed rises rapidly towards smaller scales in the non-linear regime as we have assumed here then the data on non-linear scales does not increase the sensitivity to the total neutrino mass. Assuming instead a ten times smaller theoretical error with the same scale dependence, the error on the total neutrino mass decreases moderately from sigma(Mnu) = 18 meV to 14 meV when mildly non-linear scales with 0.1 h/Mpc

Abstract: We present forecasts for the accuracy of determining the parameters of a minimal cosmological model and the total neutrino mass based on combined mock data for a future Euclid-like galaxy survey and Planck. We consider two different galaxy surveys: a spectroscopic redshift survey and a cosmic shear survey. We make use of the Monte Carlo Markov Chains (MCMC) technique and assume two sets of theoretical errors. The first error is meant to account for uncertainties in the modelling of the effect of neutrinos on the non-linear galaxy power spectrum and we assume this error to be fully correlated in Fourier space. The second error is meant to parametrize the overall residual uncertainties in modelling the non-linear galaxy power spectrum at small scales, and is conservatively assumed to be uncorrelated and to increase with the ratio of a given scale to the scale of non-linearity. It hence increases with wavenumber and decreases with redshift. With these two assumptions for the errors and assuming further conservatively that the uncorrelated error rises above 2% at k = 0.4 h/Mpc and z = 0.5, we find that a future Euclid-like cosmic shear/galaxy survey achieves a 1-σ error on M ν close to 32 meV/25 meV, sufficient for detecting the total neutrino mass with good significance. If the residual uncorrelated errors indeed rises rapidly towards smaller scales in the non-linear regime as we have assumed here then the data on non-linear scales does not increase the sensitivity to the total neutrino mass. Assuming instead a ten times smaller theoretical error with the same scale dependence, the error on the total neutrino mass decreases moderately from σ(M ν ) = 18 meV to 14 meV when mildly non-linear scales with 0.1 h/Mpc < k < 0.6 h/Mpc are included in the analysis of the galaxy survey data. notable exception of Lyman−α forest data, which gives an even lower bound of 0.17eV [10]. These constraints rely on a combination of data from Cosmic Microwave Background (CMB) experiments such as WMAP, Baryonic Acoustic Oscillations (BAOs), SuperNovae (SN) distance moduli, galaxy clustering and cosmic shear (especially from the SDSS 5 and CFHTLS 6 surveys). Data sets provided by Large Scale Structure (LSS) are particularly important, since they are able to probe scales and redshifts affected by neutrino free streaming both in the linear and non-linear regimes. Neutrino oscillation experiments provide a lower bound of 0.05eV on the total neutrino mass, meaning that the allowed range is now significantly squeezed by cosmological data, and well within reach of future planned surveys.
Several forecasts have already been published on the sensitivity of Euclid to cosmological parameters, with a focus on dark energy, modified gravity, the neutrino mass, or other extensions of the minimal ΛCDM model (see e.g. [11,12,13,14,15,16,17,18,19]). However reliable forecasts are difficult to obtain; interpreting Euclid data on small (nonlinear) scales will require a more accurate modeling of systematic effects than is currently achievable. This is true for both non-linear corrections to the matter power spectrum, and for effects specific to each survey. In the case of the galaxy redshift survey, for instance, redshift space distortions and scale-dependent bias. In the case of the cosmic shear survey, noise bias in shape measurements [20]. Some authors have pointed out that without considerable progress in modeling these effects, the sensitivity to cosmological parameters might degrade considerably (see e.g. [18]).
Current forecasts tend either to incorporate only linear scales and neglect these systematics, or to include a small range of mildly non-linear scales and model systematics by including nuisance parameters which are then marginalized over. Introducing such nuisance parameters (for instance, in order to describe redshift-space distortions) still assumes that we can predict the shape of these effects, and reduce them to a simple family of curves. Hence, this approach is not the most conservative.
On top of this, many forecasts are affected by a methodology issue: apart from two recent works [19,21], they are based on a Fisher matrix technique, whose results depend on the step chosen in the calculation of numerical derivatives of the spectrum with respect to the parameters (see e.g. [21,22]).
The present forecast has three objectives: • First, we wish to use a reliable forecast method for the sensitivity of a Euclid-like survey to ΛCDM parameters and to the total neutrino mass, based not on Fisher matrices, but on a parameter extraction from mock data with Markov Chain Monte Carlo (MCMC). This goal has also been achieved very recently by [19], although with a different approach for modeling the galaxy redshift survey. To our knowledge, the present analysis is the first MCMC forecast of a Euclid-like galaxy redshift survey using as an observable the power spectrum P (k) in wavenumber space.
• Second, we wish to incorporate non-linear corrections using the most accurate available fitting formula accounting for neutrino mass effects, namely the version of halofit [23] presented in Ref. [24]. This formula has been obtained by fitting to a suite of N-body simulations which incorporate neutrinos as free-streaming dark matter particles, using the code first presented in Ref. [25]. The error in this formula specific to the neutrino mass was estimated by Ref. [24] to be Gaussian, with squared variance where f ν = ω ν /ω m and k σ (z) is the non-linear wavenumber as defined and computed in halofit. We include this in the likelihood as a fully correlated error, as described in detail in Appendix A, associated to a unique nuisance parameter.
• In order to obtain conservative results while keeping the analysis simple, we will combine this correlated error with a second uncorrelated error. This second uncorrelated error is assumed to account for extra uncertainties in our approximate modeling of non-linear corrections, redshift space distortions, scale-dependent bias and other systematic effects. By assuming an uncorrelated error on each data point, we remain more conservative than if we marginalized over a small set of nuisance parameters representing several types of fully correlated errors. Throughout this work, we assumed for convenience that the relative theoretical error on the power spectrum was given by Eq. (1.1), with f ν replaced by a constant factor, by default 0.05. This error grows smoothly from zero on linear scales up to 5% on deeply non-linear scales. For a concordance cosmology and at redshift z = 0.5, it reaches 1% near k = 0.1 hMpc −1 and 2.3% around k = 0.6 hMpc −1 . We assume that ten years from now, this will provide a reasonable description of the total uncertainty coming from all systematic effects in each of the two surveys. Occasionally, we will consider the effect of dividing the magnitude of the error by two or ten, to evaluate the effect of better control of non-linear systematics. We emphasise that the exact form of the uncorrelated error is obviously just an educated guess and that a different k-dependence will e.g. influence the assessment of how useful pushing to smaller scales will be. Of course, introducing a fully uncorrelated error (or alternatively, form filling functions as in [26]) is very conservative in that it assumes that no modeling of systematics is accurate enough.
In several years from now, it might become realistic to model most systematics with several types of correlated errors, and to reduce the residual uncorrelated theoretical error to a smaller level than assumed in this work.

Galaxy redshift survey
Throughout this paper, our fiducial model is chosen to be a flat ΛCDM model with three degenerate massive neutrino species. The fiducial parameter values are taken to be ω b = 0.02258, ω c = 0.1109, A s = 2.43 × 10 −9 (pivot scale k * = 0.05 hMpc −1 ), n s = 0.963, h = 0.710, z reio = 10.3, m ν = 0.07 eV (so M ν = 0.21 eV). For the power spectrum of the mock data, we could take directly the fiducial power spectrum, or generate a random spectrum realization corresponding to the same model. As illustrated in [22], the two options lead to the same forecast errors, so for simplicity we assume an observed power spectrum equal to the theoretical power spectrum of the fiducial model.  The upper plots show a comparison between a model with massless neutrinos and our fiducial model (M ν = 3m ν = 0.21 eV). Solid lines are derived from the non-linear matter power spectrum using the updated halofit version of ref. [24], while dashed lines are derived from the linear power spectrum. The lower plots show the part of the relative error coming from observational or theoretical errors only (cosmic variance is included in the observational error). In these plots, the individual 1-σ error on each data point has been rescaled by the square root of the number of points, in such a way that the edges of the error bands correspond to a shift between theory and observation leading to ∆χ 2 = 1, when only the observational or theoretical error is incorporated in the likelihood expression. In these lower plots, we also show for comparison the ratio between a massless model and a model with the minimum total mass allowed by neutrino experiments, M ν = 0.05 eV.
We fit the mock and Euclid-like spectra using the MCMC code MontePython [27].
MontePython uses the Metropolis-Hastings algorithm like CosmoMC [28], but is interfaced with class [29,30] instead of camb [31], is written in python, and has extra functionality; it will soon be released publicly, including the Euclid-like likelihood codes used in this work. Technical details of the assumed likelihood and our analysis are presented in Appendix A. Let us summarize here the essential points. As in most of the recent Fishermatrix-based forecasts, we assume that the reduced data is described by a set of observable power spectra P obs (k ref , µ, z), related to the familiar non-linear matter power spectrum P N L (k, z) in a non-trivial way in order to take into account redshift space distortions, linear light-to-mass bias, spectroscopic redshift errors and the Alcock-Paczynsky effect (see A.1). Of course, this modeling is imperfect: for this reason we introduce a theoretical error. For instance, we do not take into account galactic feedback [32], assuming that this contamination can be predicted by simulations up to the level of our residual theoretical error function. The arguments k ref and µ of the observable power spectrum stand respectively for the observed wavenumber assuming the fiducial cosmology, and the cosine of the angle between the observed wavevector and the line of sight. We assume sixteen redshift bins with mean redshift ranging from 0.5 to 2, and bin widths of ∆z = 0.1. For a fixed theoretical model, each observed value of P obs in a bin centered on the point (k ref , µ, z) follows, to a good approximation, a Gaussian distribution with variance where dµ is the size of the bins in µ space, and [dk ref /k ref ] the size of the logarithmic bins in wavenumber space (see A.2). The characteristics of the survey are encoded in V survey , the survey volume, and n g , the comoving number density of galaxies accounting for shot noise (see A. 3). Hence, if for every observed data point the theory and the observation differed by this amount, the effective χ 2 would increase with respect to its minimum value by the number of data points, namely where B is the number of redshift bins.
To illustrate this error, in figure 1, we show the relative error bar on the observed spectrum in the first and last redshift bin, assuming no additional theoretical error. For the purpose of comparing with the theoretical error introduced below, we do not show as usual the error corresponding to a one-sigma deviation for each given data point; we divided each error by √ N , in such a way that the edge of the error band corresponds to a deviation between the observed and theoretical spectrum leading to ∆χ 2 = 1. Note that the displayed quantity ±∆P obs /(P obs √ N ) does not depend on the width of the bins in (k ref , µ, z) space, but only on P th , V survey and n g .
We incorporate the theoretical error in the likelihood in the way described in section A.4. In few words, this error is normalized in such a way that a shift between theory and observations by a relative amount α (the quantity defined in eq. (1.1)) leads to an increase of the χ 2 by one. This is achieved simply by adding a term N (αP th ) 2 to the total error variance. Figure 1 shows the relative theoretical error on the observed spectrum, normalized in such a way that the edge of the error band corresponds to a deviation between the observed and theoretical spectrum leading to ∆χ 2 = 1 when the observational error is switched off. These edges are directly given by ±α.
We see in this figure that our assumption for α leads to an error of 1% at k = 0.1hMpc −1 and 2.5% at k = 0.6hMpc −1 for the first redshift bin centered on z = 0.5. For the last redshift bin in the galaxy survey, centered on z = 2, non-linear corrections appear on smaller scales, and the error is only 1% at k = 0.6hMpc −1 .   Table 1: Marginalized 1-σ error for each model parameter, in a fit of Planck + Euclid-like galaxy survey data. The different lines correspond to different choices of k max , to the inclusion or not of the global uncorrelated theoretical error (un. err.), divided by ten (1/10), by two (1/2), or full (•), to that of the specific neutrino-related correlated error (co. err.), and to the use of the non-linear or linear power spectrum. The models with correlated error have one more nuisance parameter e ν not shown here, with unit 1-σ error.
We performed several forecasts for a combination of Planck data and a Euclid-like galaxy redshift survey data. It should be stressed that the characteristics of Euclid are not yet finalized. Our choice for V survey and n g (z), detailed in A.3, should be taken as indicative only. For Planck, we follow the method presented in [22] and do not include lensing extraction. For the experimental Planck sensitivity, we use the numbers presented in the Planck Bluebook 7 . This is a rather conservative model since the sensitivities are based on 14 months of observations instead of 30.
The differences between our forecasts reside in the maximum wavenumber, equal to k max = 0.1 or 0.6 hMpc −1 , and in various prescription for the theoretical error: no error at all, the uncorrelated error described above and in A.4 (divided by ten, by two, or full), or additionally the correlated error accounting for neutrino-mass-related effects (described in A.5). Since we are using an increasing theoretical error on non-linear scales, we expect the amount of information contained in the data to saturate above some value of k max : this is the reason we can consider such a high value as 0.6 hMpc −1 . We did not try even higher values, first because our result would not change, and second because our forecast would become unrealistic: deep in the non-linear regime, the Gaussian assumption for the likelihood breaks down.
Our results are presented in Table 1. Parameters like ω b and z reio are well determined by CMB data, and their forecast error depends very mildly on our different assumptions. For other parameters, the redshift survey plays a crucial role in removing parameter degeneracies. In that case, even with k max = 0.1 hMpc −1 , including the uncorrelated theoretical error makes a difference: the parameter sensitivity degrades by up to 70% for h. The 68% neutrino mass error bar degrades by 40%, from σ(M ν ) = 0.018 eV to σ(M ν ) = 0.025 eV.
Assuming only this uncorrelated error, the cases k max =0.1 hMpc −1 and k max =0.6 hMpc −1 give almost the same results. Hence, our assumption for the theoretical error magnitude is such that most of the information is contained on linear scales. Thanks to realistic (or at least conservative) assumptions for the theoretical error, the results of our forecast are nearly independent of the cut-off k max . Without a theoretical error, increasing k max to 0.6 hMpc −1 would lead to a spectacular (but totally unrealistic) decrease of the error bars, with σ(M ν ) = 0.0059 eV.
If we are more optimistic and half the uncorrelated error, the error bars decrease marginally, as can be seen in the Table (lines starting with "1/2"). The error on the neutrino mass only decrease by ∼ 10%. Assuming no error at all implies that the spectrum can be predicted up to the 0.1% level or better on small scales. In comparison, assuming a precision of one percent is not very different from assuming two percent. With the halved error, the sensitivity to the neutrino mass increases from σ(M ν ) = 0.023eV to σ(M ν ) = 0.022eV when including data in the range from 0.1 to 0.6 hMpc −1 .
Finally, in a very optimistic forecast with an error ten times smaller, we start to see how extra information can be extracted from non-linear scales; the error decreases from σ(M ν ) = 0.018 eV to σ(M ν ) = 0.014 eV when pushing k max from 0.1 to 0.6 hMpc −1 .
The inclusion of an additional correlated error accounting for neutrino-mass-related systematics has a negligible impact on our results. In our forecast, the uncorrelated and correlated part of the error have similar amplitudes and the same shape; however the uncorrelated error allows much more freedom and thus leads significantly more conservative results: this explains why the correlated error has a comparatively small effect. It should be stressed that our results depend not only on the assumed error amplitude at a given scale and redshift, but also on the wavenumber dependence of the error function α. Different assumptions, with a steeper or smoother step in the error function around the scale of nonlinearity, would lead to different forecasts. In particular, as already mentioned the actual benefit from pushing to smaller, non-linear scales depends on the assumed k-dependence of the residual uncorrelated theoretical error.
For the case with k max = 0.6 hMpc −1 and no neutrino-related correlated error, we show the one and two-dimensional posterior probability on cosmological parameters in figure 2. We see several pronounced parameter degeneracies. For instance, the neutrino mass is very correlated with ω c and h. This suggests that further progress could be made by including extra data sets, such as direct measurements of the Hubble parameter, the cluster mass function, supernovae luminosity, 21-cm anisotropies, and so on.
Our results are consistent with those of [14,17], although a direct comparison is   Table 1). difficult, since those authors include two extra parameters, w 0 and w a , in their forecast. The results in Table 2.1 of [33], based on the same cosmological model, match our prediction in the case with no non-linear scales and no theoretical error included. A similar sensitivity was found by [19] for a Euclid-like photometric redshift survey, referred to as "cg" in their Table 2. However, this reference presents other results based on even more conservative assumptions than ours. We assumed that the bias function for each redshift bin could be determined in advance (up to corrections on non-linear scales contained in our global theoretical error). This assumption has also been made in most recent forecasts, since both N-body simulations and higher-order statistics in the real data allow the prediction of the redshift-dependent bias of a given population of galaxies, at least on linear scales. Were this approach found to be unreliable, it would be necessary to marginalize over the linear bias in each redshift bin, b(z i ). Ref. [19] did such a marginalization in the runs called "cgb" and "cgbl", with no prior at all on each b(z i ). They found roughly the same error bar on ω c and h than in our forecast with theoretical error, but a much larger error on the neutrino mass. However, it seems unlikely that at the time when Euclid data will be analyzed, no information at all will be available on the linear bias of the observed population of galaxies.

Cosmic shear survey
For the case of a Euclid-like cosmic shear survey, we stick to the same fiducial model and methodology as in the previous section. The likelihood is now a function of the observed lensing power spectrum C obs ij l in harmonic space and for each pair ij of redshift bins, taking into account photometric redshift errors and shot noise (for details, see B.1 and B.2). We assume experimental sensitivities summarized in B.3, and cut the observations in five redshift bins covering the range 0 < z < 3.5 (although a negligible amount of galaxies contribute between 3 and 3.5). We do not take into account intrinsic alignment, assuming that this contamination can be removed up to the level of our residual theoretical error function [34,35].
As explained in detail in B.4, there is a small technical difference between the likelihood of the galaxy survey and the shear survey in the way we incorporate the uncorrelated theoretical error. For the galaxy survey, the theoretical error was encoded as an extra contribution to the total error variance. This can be justified mathematically by marginalizing over one nuisance parameter for each data point. The shape of the galaxy survey likelihood allows for an analytical minimization over each nuisance parameter, in such a way that nuisance parameters do not appear explicitly in the final likelihood. We found that no such scheme is accurate enough in the case of the (chi-square type) shear likelihood. Hence our likelihood routine performs an explicit minimization over one nuisance parameter per data point. For simplicity, we assume that the error is uncorrelated between different values of l, but not between different bins for a given l: this assumption could be relaxed, at the expense of increasing the computing time.
We fixed l max = 2000, since beyond this value both the shot noise term and the theoretical error are large, as shown in figure 3. This figure also shows the relative error on the observed spectrum in the first and last redshift bins, coming either from observational errors (including cosmic variance) or from the theoretical error, and using exactly the same conventions as in the previous section: the edges of each of the two error bands correspond to a shift between the theory and the observation leading to ∆χ 2 = 1 when either the observational or the theoretical error are included in the likelihood. The lowest redshift bin incorporates small non-linear scales: this explains why at l = 2000, the theoretical error reaches 3.5%.
Our results are presented in Table 2 for three cases: no theoretical error, uncorrelated error only (described in B.4), or additional neutrino-related correlated error (described in Solid lines are derived from the non-linear matter power spectrum using the recent update of halofit [24], while dotted lines are derived from the linear power spectrum. The lower plots show the part of the relative error coming from observational or theoretical errors only (cosmic variance is included in the observational error). In these plots, the individual 1-σ error on each data point has been rescaled by the square root of the number of points, in such a way that the edges of the error bands correspond to a shift between theory and observation leading to ∆χ 2 = 1, when only the observational or theoretical error is incorporated in the likelihood expression. In these lower plots, we also show for comparison the ratio between a massless model and a model with the minimum total mass allowed by neutrino experiments, M ν = 0.05 eV.
B.5). The impact of the uncorrelated error is again important, but not as pronounced as in the galaxy power spectrum case, because on small scales the precision of the shear survey is limited by a significant shot noise contribution. The neutrino mass error degrades only from σ(M ν ) = 0.026 eV to 0.028 eV. For the shear survey we did not perform runs with a twice or ten times smaller error: the result for σ(M ν ) would simply lie between those two numbers. The impact of the neutrino-related error is small but further degrades the sensitivity to σ(M ν ) = 0.032 eV. While in the absence of theoretical error the galaxy survey seems more sensitive to the neutrino mass, the performance of the two methods are roughly identical once the same theoretical error ansatz is included. un. co.    .), and of the specific neutrino-related correlated error (co. err.). Our preferred prediction is given on the last line, and is very close to that of the second line.
The triangle plot of figure 4 shows that the parameter degeneracies are very similar for the two cases of the galaxy survey and shear survey. Nevertheless, [19] showed that combining the two data sets (with a proper cross-correlation matrix) leads to sensitivity improvements. It would be interesting to test this conclusion in presence of theoretical errors.
Our result are consistent with those of [11], although a direct comparison is difficult, since these authors include several extra parameters (w 0 , w a , r, α s ) in their forecast. The predictions of [19] (case "cs" in their Table 2) lie between our results with and without theoretical errors. This is consistent since on the one hand, these authors use more optimistic survey characteristics (d, γ 2 rms , σ ph ), and on the other hand, we are including much larger values of l (which is legitimate if our theoretical error is realistic).

Conclusions
We have presented forecasts of cosmological parameters by using, in combination with Planck data, two Euclid-like mock future data sets: a galaxy spectroscopic redshift survey and a cosmic shear survey. We focused our attention on constraints that can be achieved on the total neutrino mass by using the data in the linear and non-linear regimes.
In order to do this conservatively we adopt the following improvements with respect to similar works performed recently in the literature: i) we make use of Markov Chain Monte Carlo rather than the Fisher Matrix, which results in more reliable error bars, as well as considering degeneracies between parameters. Ultimately, we found that the posterior probability is very close to a multivariate Gaussian for the model considered. However, a Fisher matrix approach could not have confirmed this, and would not have been explicitly independent of the stepsize in the numerical derivatives. ii) we rely on a modification of HALOFIT that accounts for massive neutrinos, and predicts the non-linear matter power spectrum to small scales, based on the results of N-body and hydro simulations. iii) we conservatively consider errors both on the non-linear observable power at small scales and on the neutrino induced suppression, and explictly show how to implement these errors in the likelihood calculation.
It is instructive to see that with the shape assumed for the uncorrelated theoretical error, and a conservative assumption on its amplitude (leading to a 2% error at k max = 0.4 h/Mpc and z = 0.5), the sensitivity to cosmological parameter is still satisfactory. The   Figure 4: Marginalized posteriors and two-dimensional probability contours in a fit of Planck + Euclid-like shear survey data, with a global uncorrelated error of 5% on non-linear scales (second model in Table 2). error bar on the total neutrino mass, of the order of 32 meV (cosmic shear) or 25 meV (redshift survey), would still allow for a two sigma detection of the total neutrino mass in the minimal normal hierarchy scenario. However, with this amplitude and k-dependence of the theoretical error, essentially all the information comes from linear scales. The next interesting question is to check how much the uncorrelated error should be controlled in order to start being sensitive to mildly non-linear scales. Assuming a twice smaller error does not change the parameter sensitivity by a significant amount. Extracting significant information from non-linear scales requires an error ten times smaller, at the level of 0.2%.
Here the error on the neutrino mass decreased from σ(M ν ) = 18 meV to 14 meV when adding scales with 0.1 < k < 0.6 h/Mpc to the analysis. This shows that it would be extremely useful to be able to predict the observable power spectrum of a given cosmological model up to a residual uncorrelated error of the order of 0.1% (resp. 0.2%) at k ∼0.1 h/Mpc (resp. k ∼0.4 h/Mpc) and z = 0.5. This will be a major challenge for theoretical and numerical cosmology in the next decade.

A.1 Observed spectrum
Let P obs be the observed/mock/fiducial power spectrum, and P th the spectrum that one would expect to see given the theoretical model. Each of these quantities relates to the galaxy spectrum P g and finally to the total non-linear matter spectrum P NL by taking into account redshift distortion effects, spectroscopic redshift errors and light-to-mass bias. A good approximation of such a relation is given by (see e.g. [16,15]): with the definitions Here b(z) is the bias, assumed to be scale-independent in the range of scales of interest, a is the scale factor, H(z) is the Hubble parameter, D A (z) the angular diameter distance, and β(z, k) accounts approximately for redshift space distortions. So we can treat k as a function of the arguments (k ref , µ, z) and write

A.2 Likelihood
For a narrow redshift bin b centered onz, the likelihood reads with an effective survey volume given by (A.10) Later, we will specify the sensitivity of the survey, parameterized by V survey , n g , σ r , k min and k max . We skip here the derivation of the Fisher matrix, obtained by differentiating the above formula twice with respect to the cosmological parameters on which P th depends, and evaluating this derivative at the maximum likelihood point. We checked that this calculation gives exactly the formula commonly used in the literature (see e.g. [16,15]). For the purpose of the discussion in the next section (and also of the numerical implementation), we wish to write explicitly the discrete limit of the integrals. We discretize µ in a set of equally spaced values µ i , and l ≡ ln k in a set of equally spaced values l j = ln k refj . The step sizes are denoted ∆µ and ∆l respectively. We then expand the integral as a sum, and for simplicity we omit the factors 1/2 that should weight the boundary terms of each of the two integrals. We introduce the short-cut notations: and we get This expression is easy to understand from first principles. Let us consider a single variable δ obeying a Gaussian distribution centered on zero and with variance δ 2 = P . If we observe N independent realization δ n of the variable δ, we can build an estimator of the variance P of δ, 14) The variance of this estimator can be computed by noticing that each δ 2 n follows a χ 2 distribution of order one, for which the mean is P and the variance 2P 2 . So the sum n δ 2 n has a variance 2N P 2 . Finally E has a variance (2N P 2 )/N 2 = 2P 2 /N . Moreover, E is nearly Gaussian if N is large, as a consequence of the central limit theorem. So the probability of the data E given the theory P is a Gaussian of mean P and of variance 2P 2 /N . In other words, The previous likelihood follows this form for each discrete term. Indeed each term corresponds to the likelihood of the estimator of the power spectrum in a thin shell in Fourier space. The number of independent measurements, i.e. of independent wavenumbers in each shell, is given by N ij . The role of E and P is played respectively by P obs ij and P th ij . Such a likelihood was first derived in pioneering papers like [36,37].

A.3 Survey specifications
We computed this likelihood for values of V survey (z), n g (z), σ r (z) inspired from currently plausible Euclid specifications, which are likely to change over the next years. We divide the observations into sixteen redshift bins of width ∆z = 0.1, ranging fromz = 0.5 tō z = 2.0. For each bin, we assumed: • a volume per bin V survey (z) = 4πf sky [r(z)] 2 (1+z) −3 ∂r(z) ∂z ∆z, where r(z) is the comoving distance up to a comoving object with redshift z, with the explicit assumption that a 0 = 1: .
We assume a sky coverage f sky = 0.375.
• a galaxy number density per comoving volume n g (z), related to the number of galaxies per square degree d g (z) through For d g (z), we start from the number presented in Table 2 of [38] for the case of a limiting flux of 3 × 10 −16 erg s −1 cm −2 . Following the recommendation of that paper, we divide these numbers by 1.37 in order to get conservative predictions. Finally, we multiply them by an efficiency factor = 0.25 (standing for the redshift success rate). For instance, for the first redshift bin, this gives d g (z) = 9376/1.37 × 0.25 = 1710 deg −2 .
• a scale-independent linear bias b(z). The choice of b(z) values affects the final result less crucially than that of d g (z). We could adopt the predictions of [39] inferred from N-body simulations, but for simplicity, our forecast is performed under the approximation b(z) = √ 1 +z. So, we assume in this forecast that the linear bias will be accurately measured or predicted for each bin, and that deviations from this prediction (coming from the non-linear evolution) will be known up to the level described by the theoretical error function.
• k min can be chosen arbitrarily close to zero without changing the results.

A.4 Accounting for a global uncorrelated theoretical error
To present a realistic forecast, one should model all the systematic effects not accounted for by the previous likelihood formula, such as: theoretical errors in the calculation of the linear and non-linear power spectrum, scale-dependence of the bias on small scales, residual shot noise in galaxy counts beyond the contribution already included in the definition of V eff , residual errors in the modeling of redshift space distortion beyond the above scheme. On top of these corrections, one may have to take into account the fact that the likelihood is not Gaussian on strongly non-linear scales. In this paper, we limit ourselves to mildly non-linear scales k ≤ k max = 0.6hMpc −1 , and assume that non-Gaussianity effects are sub-dominant to the previously mentioned systematics. We also neglect to marginalize over residual shot noise in each redshift bin, because Ref. [16,15] found that this has a negligible impact.
Understanding these various systematics is a major challenge for the future, which should be addressed with better simulations and analytical modeling. Here we want to keep the analysis simple, and model these systematic errors in a simple way, by adding to the spectrum an uncorrelated theoretical error function. By uncorrelated we mean that the errors made at different scales are independent from each other, which is the most conservative possible assumption. In this case, we can introduce an independent Gaussiandistributed nuisance parameter ij for each data point, and marginalize over it -or rather, to a very good approximation, minimize over it: where R ij is the theoretical error variance for a bin in (µ, k ref ) space centered on (µ i , k refj ). As long as the theoretical error is assumed to be small, it is also a valid approximation to neglect the ij -dependence of the denominator, in order to find a simple analytic solution for ij , which, injected back in eq. (A.18), gives In other words, the theoretical error variance simply adds up to the noise variance. Note that we explicitly checked that it is legitimate to neglect the ij -dependence of the likelihood denominator when minimizing over ij . We also coded the full likelihood with explicit minimization over each ij , and found the same results up to very good accuracy.
We choose a numerical value of R ij motivated mainly by the current level of precision of the halofit algorithm. We assume a relative error on the non-linear power spectrum of the form where k σ (z) is the scale of non-linearity computed by halofit. This function increases from zero to 5% around the scale of non-linearity. Using the function k(k ref , µ,z), this error can easily be propagated to the theoretical observable spectrum In terms of the discretized observable spectrum, the error reads The error variance R ij should be proportional to the power spectrum variance (α ij P th ij ) 2 . We also assume that the error makes a constant contribution to each logarithmic interval in the space where observations are performed, i.e. is of the form We normalize the error variance R ij in such a way that a one-sigma theoretical error in each data point results in increasing the effective χ 2 by one unit, namely, where B is the number of bins. The role of the normalization factor between squared brackets will become clear below. The likelihood becomes (using eq. (A. 19) and going back to the continuous limit) where we omitted the argument (k ref , µ,z b ) of the functions P th , P obs , V eff and α. If one assumes that the observed and theoretical spectra differ by αP th for each (k, µ, z), and that in the denominator the theoretical error dominates over the observational one (V eff = ∞), then which corresponds to a shift by ∆χ 2 eff = 1 with respect to the maximum likelihood L = Π b N b .
If we had assumed the error to be fully correlated, instead of increasing the denominator of the likelihood, we would have replaced P th by P th (1 + α), multiplied the likelihood by 1/2π exp[− 2 /2], and marginalized/minimized over . Then, the assumption P obs = P th (1 + α) would correspond to an optimal choice = 1 in the large V eff limit, and would also lead to a shift in ∆χ 2 eff by one unit with respect to the assumption P obs = P th . In our case, we obtain the same shifting while assuming statistically independent errors for each data point.
Finally, the likelihood can be simplified to where we omitted the argumentz b in the functions V survey , D A , H and n g . This is exactly the relation implemented in our code.

A.5 Accounting for an extra neutrino-related error
The impact of massive neutrinos on non-linear corrections to the power spectrum has been investigated in [24]. By comparing with N-body simulations including neutrino particles, the authors of [24] re-calibrated halofit, with a new neutrino mass dependent correction. This fitting procedure is of course not perfect and adds a systematic error growing with the neutrino mass. It was found that the leading error can be described with a correction with f ν ≡ ω ν /ω m , and e ν is an unknown correction of unit variance, that we will treat as a Gaussian nuisance parameter. Hence our final definition of the likelihood accounting for both types of error reads where we omitted the argument (k ref , µ,z) of the functions P obs , P th , σ ν , α ν and V eff . Note that the correction proportional to e ν should not be added to P obs since we are assuming for simplicity that the fiducial value of e ν in the mock data is zero.

B.1 Observed spectrum
As in e.g. [16,15], we define the likelihood of the shear auto or cross-correlation power spectrum in bins i and j: Here, W i (z) is the window function of the i'th bin. It can be evaluated as a function of the radial distribution of galaxies in each redshift bin, D i (z), obtained by convolving the full radial distribution D(z) with the photometric redshift uncertainty function P(z, z ph ), multiplied the top-hat window function of each bin: The radial distribution D(z) can be arbitrarily normalized, since n i (z) is anyway normalized to one. We will assume that the photometric redshift uncertainty function is normalized to ∞ 0 P(z, z ph )dz ph = 1, but a different normalization would not impact the final result for the same reason as for D(z). The noise spectrum contaminating the measurement of C ij l is given by the diagonal matrix in ij space: where γ 2 rms 1/2 is the root mean square intrinsic shear (like in the forecasts of the Euclid Red Book [1], we assume that this quantity is equal to 0.30), and n i is the number of galaxies per steradian in the i'th bin, given by where d is the full number of galaxies per square arcminute in all bins, andn i is the fraction of galaxies in the i'th bin, given by:n where a obs lm = a obs i lm is the N-dimensional vector of observed multipoles in each bin,C th l is the theoretical covariance matrix of element C th ij l and N is a normalisation factor. After some simple algebra 8  We assume five bins, with the first bin starting at z min 1 = 0, the last one ending at z max N = 3.5, and bin edges z min i = z max i−1 chosen such that each bin contains the same number of galaxies, i.e.n i = 1/N .

B.5 Accounting for an extra neutrino-related error
Finally, we account for the correlated error modelling neutrino-related uncertainties by multiplying the theoretical power spectrum P th (k, z) by a factor (1+e ν σ ν (k, z)), as in equation (A.28), as well as adding e 2 ν to ∆χ 2 eff . The nuisance parameter e ν is then marginalized over. Note that the factor (1 + e ν σ ν (k, z)) should not multiply the observed/fiducial spectrum, as long as we assume a fiducial value of e ν equal to zero.