Red, Straight, no bends: primordial power spectrum reconstruction from CMB and large-scale structure

We present a minimally parametric, model independent reconstruction of the shape of the primordial power spectrum. Our smoothing spline technique is well-suited to search for smooth features such as deviations from scale invariance, and deviations from a power law such as running of the spectral index or small-scale power suppression. We use a comprehensive set of the state-of the art cosmological data: {\it Planck} observations of the temperature and polarisation anisotropies of the cosmic microwave background, WiggleZ and Sloan Digital Sky Survey Data Release 7 galaxy power spectra and the Canada-France-Hawaii Lensing Survey correlation function. This reconstruction strongly supports the evidence for a power law primordial power spectrum with a red tilt and disfavours deviations from a power law power spectrum including small-scale power suppression such as that induced by significantly massive neutrinos. This offers a powerful confirmation of the inflationary paradigm, justifying the adoption of the inflationary prior in cosmological analyses.

power law, using a set of theoretically motivated shapes for the PPS, and a minimally parametric analysis to reconstruct the PPS. In all cases there is no strong evidence for deviations from a power law.
Our analysis differs from that of the Planck collaboration and from others existing in the literature as we analyse jointly a comprehensive set of state-of-the-art experiments probing the matter power spectrum and the latest Planck measurements.
Because we assume standard late-time evolution of density perturbations and consider both early-time observables (CMB) and late-time ones (i.e., large-scale structure), our reconstruction is also sensitive to late-time effects on structure formation. In particular a non-negligible neutrino mass would suppress the growth of structures below the neutrino free-streaming scale, inducing an "effective" loss of small scale power in our reconstructed PPS. Reconstructing in a model-independent way a possible neutrino signature on the shape of the matter power spectrum is of particular importance as [11][12][13][14][15][16] claims that relatively large neutrino masses (Σ ν 0.4 eV) could solve the tension between CMB and local measurements, whilst other studies [17][18][19][20][21][22][23][24] rule out this possibility.
The rest of the paper is organised as follows: in section 2 we briefly summarise the methodology adopted, the data chosen and how they are analysed. In section 3 we present our findings; we discuss and present the conclusions in section 4.

Spline reconstruction
We perform a minimally-parametric reconstruction of the primordial power spectrum based on the method presented in [7] and further refined in [8,9]. Here we only briefly summarise the approach; it is based on the cubic smoothing spline technique (for details see [25]). In this approach to recover a smooth function f (x), given its value f i only on a set of n points x i , hereafter knots, one fits the pairs (x i , f i ) with a cubic spline s(x). The spline, its first, and second derivatives are continuous on the knots by definition. The full function is then uniquely defined by the values at the knots and two boundary conditions. We choose to require that the jump in the third derivative across the first and last knots is forced to zero.
In our application the resulting spline function is the reconstructed primordial power spectrum. The f i are free parameters we wish to determine and we place the knots equally spaced in log k as it is the most conservative choice to recover deviations from a power law. The whole s(k) is used as the PPS to compute the observables and evaluate the likelihood of the parameters f i . Including the roughness penalty, the effective likelihood becomes − log(L) = − log(L exp ) + α p where s denotes the second derivative of s with respect to ln k, k i and k f are respectively the position of the first and of the last knots, α p is a weight that controls the penalty, and L exp is the likelihood given by the experiments. The roughness penalty effectively reduces the degrees of freedom, disfavouring jagged functions that "fit the noise". As α p goes to infinity, one effectively implements linear regression; as α p goes to zero one is interpolating. The use of cubic spline -instead of other possible interpolating functions -is motivated by the fact that such a function minimises the roughness penalty for a given set of knots (f i , x i ).
In generic applications of smoothing splines, cross-validation is a rigorous statistical technique for choosing the optimal roughness penalty [25]. Cross-validation (CV) quantifies the notion that if the PPS has been correctly recovered, we should be able to accurately predict new independent data. To make the problem computationally manageable, we adopt the following. We split the data set in two halves A and B. A Markov chain Monte Carlo (MCMC) parameter estimation analysis (for a given roughness penalty) is carried out on A, finding the best fit model. Then the − log likelihood of B given the best fit model for A, CV AB , is computed and stored. This is repeated by switching the roles of the two halves, obtaining CV BA . The sum CV AB +CV BA , gives the "CV score" for that penalty weight. With this construction, the smoothing parameter that best describes the entire data set is the one that minimises the CV score. The cross validation data sets are described below (see table  1).
We choose to use 5 knots equally spaced in log k between k = 10 −5 Mpc −1 and k = 1 Mpc −1 , i.e., k 1 = 10 −5 Mpc −1 , figure 1 bottom panel for knots placement visualisation). The number and position of the knots is held fixed throughout the analysis. As discussed in reference [8], beyond a minimum number of knots, there is a trade-off between the number of knots and the penalty, and the form of the reconstructed function does not depend significantly on the number of knots beyond this minimum number. As the main goal of this work is to explore, in a minimally parametric way, smooth deviations from a power law, a few (> 3) knots are sufficient.
The basic cosmological parameters, ω b = Ω b h 2 , ω c = Ω c h 2 , h, and τ reio -physical baryonic matter density parameter, physical cold dark matter parameter, dimensionless Hubble parameter and optical depth to last scattering surface -are varied in the MCMC alongside the values f i of the reconstruction at the knots. A flat geometry is assumed so that Ω m + Ω Λ = 1.
The prediction for cosmological observables, the calculation of the likelihood and the MCMC parameter inference are implemented using the standard Boltzmann code CLASS [26] and its Monte Carlo code, Monte Python (MP) [27], suitably modified. 12 Even though we reconstruct the primordial power spectrum, we are sensitive to late-time cosmological effects. Our main focus is on massive neutrinos: the presence of non-negligibly massive neutrinos would distort our reconstruction in a way that is predictable due to the linearity of the growth functions [28] (see Appendix). Thus in the analysis we will assume massless neutrinos.
Of course neutrino masses do not actually affect the physical PPS. But assuming standard gravity, standard growth of structure, and massless neutrinos in the analysis, would yield a reconstructed PPS with an artificial distortion, if neutrino masses were not negligible. In fact a detectable signature of massive neutrinos in the real data would appear as a power suppression feature in the reconstructed PPS. Of course a detection of power suppression cannot be univocally interpreted as signature of neutrino masses; other particles beyond the standard model could easily share the same properties of neutrinos when it comes to damping perturbations or it could be a real feature in the PPS.

Datasets
We use a comprehensive set of power spectra obtained from observations of CMB and of large scale structure (including both weak gravitational lensing and galaxies redshift surveys) as follows: • Planck power spectra of temperature and polarisation of the CMB. The Planck collaboration released in 2013 the temperature data from the first half of the mission [29]. We complement the Planck 2013 data with the WMAP polarisation. We refer to this as PlanckCMB13. In 2015 the results of the full analysis has been released [30]. Temperature and E-mode polarisation power spectra (and their cross-correlation) data and likelihoods come in two sets: a low from = 2 to = 49, and the high angular power spectrum. We use the temperature and polarisation data up to = 2500 and we refer to this as PlanckCMB15.
• Beside the CMB power spectrum, Planck reconstructed the CMB lensing potential [31], which contains information on the amplitude of large scale structure integrated from recombination to present time. We will refer to it as PlanckLens.
• The Canada-France-Hawaii Lensing Survey (CFHTLenS) [32] provides the two point correlation function of the tomographic weak lensing signal. • The WiggleZ Dark Energy Survey (WiggleZ), through the measurement of position and redshift of 238,000 galaxies, mapped a volume of one cubic gigaparsec over seven regions of the sky up to a redshift z 1. The corresponding galaxy power spectrum is presented in [33].
• The Sloan Digital Sky Survey collaboration, in Data release 7 (SDSS DR7), used a sample of luminous red galaxies to reconstruct the halo density field and its power spectrum roughly between k = 0.02 h/Mpc and k = 0.2 h/Mpc [34].
In figure 1 we show the scales probed by each experiment along with the location of the knots.

Runs set-up
We now describe the cross validation set up. In order to constrain both the shape of the PPS and the cosmological parameters, we have to consider CMB primary data in all CV runs. Because of time constraints PlanckCMB2013 is used in the set up CV runs but PlanckCMB2015 is used in the final run. This choice is conservative, favouring slightly more freedom (lower penalty) to the reconstructed PPS. Besides these, we have 4 other experiments: 2 measuring weak lensing and 2 using galaxy catalogues. We perform 3 CV runs in a pyramidal scheme as summarised in table 1. We start performing in parallel two different cross-validation analysis on two pairs of experiments where each pair is formed by a weak lensing experiment and by a galaxy catalogue. The dependence of the CV score on α p was mapped by sampling several α p values. The results of these preliminary runs show no unexpected behaviour or tension, i.e., the reconstructed PPS shows no significant deviation from a power-law, and the shape of the CV score is the same for both run 1.1 and run 1.2. Knowing this, we then combine the large scale structure data to have one weak lensing and one galaxy survey in each CV set. The best roughness penalty found from this CV is used in the final run which includes all experiments (this is called "Rec." run in the table). The penalty parameter value to use in the reconstruction is determined by the CV score of run 2 alone: its dependence on α p is illustrated in figure 2. The fact that the shape of the three CV scores -from run 1.1, 1.2, and 2 -shown in figure 2 is very similar, indicate robustness and that there are no significant tensions between the datasets.
The CV score has a fairly well defined "wall" for high penalties , but is quite constant under a certain threshold at α p ∼ 10. For high α p the penalty starts being the dominant contribution to the likelihood, so the behaviour in the limit of high α p is expected. On the other hand, if small values of the penalty were to lead to overfitting, the CV score should increase as α p decreases. This is not what we see and can be understood as follows. CMB angular power spectra are always included in the analysis and in this limit, it is the statistical power of these data (not the penalty) that drives the smoothness of the reconstruction and therefore the CV score. In other words, for low values of the penalty below α p ∼ 10, all datasets are well consistent with the Planck-inferred PPS reconstruction: the CMB data alone disfavour unnecessarily wiggly shapes, even when there is a low penalty.
Since there is not a well defined minimum for the CV score, we opt for presenting two different cases. One is more conservative, in the sense that it has a stronger penalty that allows only small deviations from the concordance power-law model. For this one we choose α p = 1.
The other leaves more freedom to the data, as we choose a more relaxed penalty α p = 0.01. A reconstruction with α p 0.01 is pretty much uninformative. In fact recall that the free parameters in our MCMC runs are the physical baryon density ω b , the physical cold dark matter density ω cdm , the rescaled Hubble parameter h, the optical depth to reionization τ reio , and the value of the five knots of the spline that we used to parametrize the shape of the PPS. At such low penalty values the reconstruction transfers in part the features of the radiation transfer function and the effect of the optical depth to reionization into the PPS opening up degeneracies in parameter space.

Results
Here we present the results with the latest Planck likelihood (2015 release) and all the large scale structure power spectrum data (Planck Lensing 2015, WiggleZ, CFHTLenS, and SDSS DR7), with the two different roughness penalties (α p = 1 and α p = 0.01) justified above.
As discussed in refs. [17,[35][36][37][38] there is a tension between the inferred matter power spectrum amplitude from CMB and from CFHTLenS, which may arise from possible systematic errors in the photometric redshifts of CFHTLens. For this reason we present results first without and then with CFHTLens.

Reconstruction without CFHTLens
In figure 3a and 3b we show the reconstructed PPS for α p = 1 and α p = 0.01 respectively. The colour-bars on the upper side show the scales probed by each experiment as in figure 1,   green for PlanckLens, red for WiggleZ, gold for SDSS DR7. PlanckCMB15 covers the whole plot. The best fit reconstruction is shown in yellow and errors are shown by plotting in dark blue (light blue) a random sample of 400 reconstructions chosen among the 68.27% most likely points (points in the range 68.27% -95.45%) in the MCMC. The 95.5% confidence regions appear to coincide with the 68.3%: this is because the reconstructed spectra are simply more wiggly and are not allowed to deviate more, and consistently across scales, from the best fit. In the figure the red and pale red regions show the 68 and 95% confidence intervals for the standard power law ΛCDM Planck 2015 TT, TE, EE + Low P analysis [2].
Note that for the more conservative choice of the penalty, errors of the reconstructed PPS are comparable with errors from Planck parametric fit at all scales. For the less conservative penalty this is also true on scales corresponding to > 30. This did not happen with the previous generation of cosmological data (see [9]) where the reconstructed PPS was significantly less constrained than with a power law fit.
The additional freedom in the PPS allowed by the lower penalty α p = 0.01 is used on scales corresponding to low CMB multipoles < 30. These scales are dominated by cosmic variance and are known to be lower than the standard ΛCDM prediction e.g., [29,[39][40][41] and refs therein.
In figure 4a and 4b we also show the reconstructed n(k) ≡ d ln P (k)/d ln k (α p = 1 and α p = 0.01) for ease of comparison with the standard power law results. 3 We find no evidence that any scale dependence of the power spectrum spectral slope is necessary, which is in agreement with previous analyses. However with this new data set we find that n = 1 is highly disfavoured by the data, in particular for α p = 1 the significance of the departure from scale invariance is comparable with that obtained when adopting the "inflation-motivated" power-law prior. Even for the more flexible reconstruction, not even one point of the more than 4 ×    paradigm, justify the adoption of the inflationary prior in cosmological analyses. Finally in figs. 5a and 5b we show the ratio of the reconstructed PPS to the best fit Planck 2015 (temperature, polarisation, and lensing) power law model.
The reconstruction is fully compatible with the parametric fit. The figure also shows the expected effect of small scale power suppression due to massive neutrino free-streaming for two representative values of neutrino masses Σm ν = 0 eV and 0.2 eV. The two models are the conditional (i.e., keeping Σm ν fixed at the required value) best fit to the data (Planck 2015 TT, TE, EE + Low P + Lensing + BAO + JLA + H 0 data). Clearly models with Σm ν > 0.2 eV are highly disfavoured by the data even with this minimally parametric reconstruction: not a single step of a 4 × 10 5 size MCMC goes near the Σm ν = 0.2 eV line. This of course does not exclude the -admittedly contrived -case with a arbitrarily large neutrino mass    inducing a small scale power suppression which is cancelled by a compensating boost of the PPS on the same scales. Occam's razor disfavours this scenario.

Reconstruction with CFHTLens
The reconstructed P (k), n(k) and P (k) relative to the power law best fit are shown in figures 6, 7, and 8 using the same conventions as in figures 3, 4, and 5. Comparison with the results of section 3.1 (in figures 3,4,5) shows that qualitatively the reconstructions are very similar, there is no strong evidence for deviations from the power law behaviour and scale invariance is still excluded. However quantitatively some differences may be appreciated. Adding the CFHTLenS datasets has the effect of lowering the overall PPS normalisation (clearly visible by comparison with figure 3).
The ratio with Planck power law best-fit in figure 8 highlights how, independently from our choice of datasets, high neutrino masses are disfavoured. Quantitatively the Σm ν > 0.2 eV bound is excluded at more than 95% confidence if we assume a power law PPS, as discussed in section 3.1.
For completeness we also report the recovered values and errors for all the model parameters in table 2 and table 3 for the two penalties α p = 1 and α p = 0.01 respectively.    The degeneracies among the parameters for the PPS value at the knots can be appreciated in the triangle plots of figure 9a for α p = 1 and figure 9b for α p = 0.01. Correlations with and among the cosmological parameters not shown are negligible. As expected, higher penalty induce correlations among the knots which are stronger between neighbouring ones.
Interestingly the only cosmological parameter that correlates with the knots is τ reio ,  Table 3. Best fit, mean and confidence intervals for the MCMC parameters in the reconstruction with α p = 0.01 which show degeneracy with the knots at higher k ( figure 9c and 9d). This behaviour is however not unexpected. The τ reio parameter only affects the CMB and in particular its main effect is to suppress the temperature power spectrum at multipoles 80. With our choice for the location of the knots, the most affected knots are therefore K 4 and K 5 . Improved polarisation data at low should reduce this degeneracy. The figure excluding the CFHTLenS dataset is qualitatively very similar and thus is not shown here.

Discussion and conclusions
The analysis of the latest cosmological data [2] indicates a highly significant deviation from scale invariance of the primordial power spectrum (PPS) when parameterized by a power law or by a spectral index and a "running". This offers a powerful tool to discriminate among theories for the origin of perturbations and among inflationary models. In fact, the deviation from scale invariance of the PPS is a critical prediction of inflation and is the only signature that is generic to all inflationary models. It is therefore a vital test of the inflationary paradigm.
One may wonder if a strong theory prior on the form of the power spectrum, such as the power law prescription, can lead to artificially tight constraints or even a spurious detection of a deviation from scale invariance, if the adopted model were not a good fit to the data.
Here we have built on the work of [7,8] to reconstruct the PPS with a minimally parametric approach, using the cross-validation technique as the smoothness criterion. We (d) α p = 0.01. consider a comprehensive set of state-of-the art cosmological data including probes of the Cosmic Microwave Background, and of large scale structure via gravitational lensing and galaxy redshift surveys. While the spline reconstruction used here is best suited for smooth features in the PPS, it is also sensitive to sharp features if they have high enough signal-tonoise. We find that there is no evidence for deviations from a power law PPS, and that errors of the reconstructed PPS are comparable with errors obtained with a power law fit. These results should be compared with those presented in [9], to appreciate the increase in statistical power brought about by the latest generation of experiments. In fact with current data a scale-invariant power spectrum is highly disfavoured even with this minimally parametric reconstruction. In particular for our conservative choice of smoothness penalty parameter values the significance of the departure from scale invariance is comparable with that obtained when adopting the "inflation-motivated" power-law prior. Constraints no longer relax significantly when generic forms of the PPS are allowed.
Because of its flexibility, our reconstruction would be able to detect the tell-tale signature of small scale power suppression induced by free streaming of neutrino if they are sufficiently massive. Of course in reality the suppression happens in the late-time power spectrum, not in the primordial one. But as we do not include the effect of neutrino masses in the matter transfer function, the reconstruction would recover an "effective" small scale damping. Our reconstruction detects no such signature, ruling out a model with a power law PPS and sum of neutrino masses of 0.2 eV or larger.
Our results, which recover in a model independent way a power law power spectrum with a small but highly significant red tilt, offer a powerful confirmation of the inflationary paradigm, justifying adoption of the inflationary prior in cosmological analyses.
(A.2) Then the evolution takes into account the presence of massive neutrino. (B. 2) The initial condition is a power-law multiplied by the neutrino power suppression. Then the evolution equation with massless neutrino is used.  Appendix: Reconstruction sensitivity to non-primordial effects Figure 10A.1-10R.2 visualises the concept -exploited here -that in the reconstructed power spectrum the effect of a non-zero neutrino mass is degenerate with a power suppression. This is a good approximation especially on scales where the evolution is linear or mildly nonlinear, i.e., k < 0.2 h/Mpc. Consider a ΛCDM universe with massive neutrinos, where all the cosmological parameters are known and with a power law PPS at the end of inflation (figure 10A.1). From these initial conditions we evolve the perturbations assuming massive neutrino (different values for the total mass are shown). On small physical scales neutrino free streaming [42] suppresses power 10A.2) yielding a resulting power spectrum shown in figure 10R.1. Now we can think of an alternative method: implement the neutrino power suppression (figure 10B.1) directly on the initial PPS as a deviation from a power law as shown in figure 10B.2. This initial power spectrum is then evolved assuming massless neutrinos. The linearity of the perturbation evolution equations guarantees that the generated matter and CMB power spectra would be the same as in the first case ( figure 10R.1). In figure 10R.2 we can appreciate the fact that discrepancies in the prediction made in the two cases come from non-linearities. For example, when considering neutrino with Σm ν = 0.4 eV we expect a small scale suppression in the linear power spectrum of 15% ( figure 10A.2). The differences due to non linearities exceed 1% only above k = 0.1 h/Mpc, and are never more than 4% at the scales of interest. This means that in principle this approach we should be able to distinguish the effect for Σm ν ≥ 0.06 eV, even though it would prevent us to obtain an unbiased measure in case of detection. Another source of error that might contribute is given by the use of spline with a limited number of knots. If the number of knots, or their position, is not suitably chosen, one could be unable to reconstruct a given signal. It is not our case, with our choice of knots we have verified that we can reconstruct any neutrino power suppression with a 10 −3 accuracy.