{'fraction_non_alphanumeric': 0.03544923714447165, 'fraction_numerical': 0.012657750988886796, 'mean_word_length': 4.638487680543755, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Features in the primordial power spectrum? A frequentist analysis

Features in the primordial power spectrum have been suggested as an explanation for glitches in the angular power spectrum of temperature anisotropies measured by the WMAP satellite. However, these glitches might just as well be artifacts of noise or cosmic variance. Using the effective Delta chi^2 between the best-fit power-law spectrum and a deconvolved primordial spectrum as a measure of"featureness"of the data, we perform a full Monte-Carlo analysis to address the question of how significant the recovered features are. We find that in 26% of the simulated data sets the reconstructed spectrum yields a greater improvement in the likelihood than for the actually observed data. While features cannot be categorically ruled out by this analysis, and the possibility remains that simple theoretical models which predict some of the observed features might stand up to rigorous statistical testing, our results suggest that WMAP data are consistent with the assumption of a featureless power-law primordial spectrum.


Introduction
The advent of a multitude of cosmological precision data in the past decade has led to the emergence of the so-called concordance, or "vanilla" model of cosmology. A key part of the vanilla model is the assumption that the spectrum of primordial curvature perturbations is smooth and featureless, and can be described by a simple power-law -consistent with their origin from an earlier period of slow-roll inflation [1][2][3][4].
Nonetheless, the inflationary mechanism also allows for more complex shapes of the spectrum, caused for instance by non-standard initial conditions [5][6][7][8], or unusual dynamics of the inflaton field due to, e.g., a phase transition [9,10], nonsmoothness of the inflaton potential [11][12][13] or particle production [14,15] during inflation (see also [16][17][18][19][20][21][22]). In any case, any detection of features in the spectrum, i.e., a deviation from the standard power-law behaviour, would yield invaluable clues on the physics of the early Universe. Presently, the most powerful single source of information about the primordial state of perturbations are the observations of the temperature anisotropies of the Cosmic Microwave Background (CMB) by the WMAP satellite [23][24][25]. Since the first data release, a number of efforts have been undertaken to ascertain the compatibility of the data with the power-law paradigm.
Generally, these results have shown that by introducing suitable features in the primordial spectrum, the fit to the data can be improved by ∆χ 2 eff ∼ O(10) over the power-law fit. Taken by itself, however, ∆χ 2 eff does not really help answer the crucial question whether the better fit indicates a real feature in the primordial spectrum or whether it just stems from over-fitting the scatter in the data from noise and cosmic variance. Indeed, attempts at interpreting results of parameter inference in terms of 'ruling out the power-law model in favour of a feature model' are likely to be prone to grossly overestimating the 'evidence' for features. A simplistic goodness-of-fit χ 2 /dof analysis would not be very enlightening either, due to the large number of degrees of freedom involved. One might thus be tempted to turn to the formalism of Bayesian model selection [54,55]. Though aside of the potential technical problems of evaluating Bayes factors for feature models with a large number of free parameters, the interpretation of results from a Bayesian model selection analysis may be problematic [56], owing to the absence of a unique well-motivated choice of priors for these mostly empirical parameterisations.
We believe that for this particular problem the toolbox of frequentist statistics provides a useful implement, and propose to perform an analysis based on the tech-nique of hypothesis testing. Given the null hypothesis of an underlying power-law spectrum, this involves the evaluation of a suitably chosen statistic on a large number of Monte Carlo simulated mock CMB temperature anisotropy data sets. From the resulting frequency distribution of this statistic one can derive a p-value, i.e., the probability that, given the null hypothesis, the value of the statistic is larger than the one observed. We note that in defining the statistic one has to beware of a posteriori interpretations of the data; a particular feature observed in the real data may be very unlikely (and lead to a low p-value; see e.g., [57,58]), but the probability of observing some feature may be quite large.
In Section 2 (and in the Appendix) we describe in detail the methods used in our analysis, the results of which are presented in Section 3. We discuss and interpret our findings in Section 4.

Null hypothesis
We assume as our null hypothesis that the temperature angular power spectrum inferred from WMAP can be explained by the present cosmological standard ("vanilla") model. In particular, the underlying primordial spectrum of curvature perturbations is taken to have a smooth power-law form determined by the two parameters n S , the scalar spectral index, and A S , the amplitude of fluctuations at the pivot scale k * = 0.05 Mpc −1 . The remaining four free parameters of the vanilla model are the baryon and cold dark matter densities ω b and ω c , the ratio of sound horizon to angular diameter distance at decoupling θ, and the reionisation optical depth τ .
The fiducial spectrum C fid used to generate random realisations of WMAP5 data should be chosen as the maximum likelihood spectrum. However, since temperature data alone are not very sensitive to τ , using temperature data alone leads to a value of τ which is at odds with results of a combined temperature/polarisation analysis. We therefore fix τ to the best-fit result for the full WMAP5 likelihood function, and then determine the values of the other five parameters by fitting the thus reduced fiveparameter vanilla model to the temperature data only (using the likelihood function described in Appendix A.1).

Alternative hypothesis
The null hypothesis is to be tested against an alternative hypothesis. Our alternative hypothesis is that the primordial spectrum is not given by a power-law, but possesses features of some sort. It now remains to define a suitable test statistic that can be used to assess the validity of the null hypothesis.

Statistic
By relaxing the assumption of a precise power-law form of the primordial spectrum one can generally achieve a better fit to the data. A commonly used measure of this improvement is the effective delta-chi-squared where L V max is the maximum likelihood of a fit to the vanilla model, and L X max the maximum likelihood of a fit to a model with an alternative shape of the spectrum.
In the present work we wish to allow a very general form of the primordial spectrum, so we focus on a reconstruction of the primordial spectrum from a binned version of the temperature angular power spectrum with the aid of a modified Richardson-Lucy (RL) deconvolution algorithm as performed in Refs. [45,48] (see Appendix A.2). Given a transfer function T (k) and an observed or simulated angular power spectrum C , this algorithm results in an "optimised" primordial spectrum P R RL (k), and is a very powerful tool to find potential features in the spectrum. We take L X max to be the likelihood of the data for the primordial spectrum P R RL (k) combined with the best-fit transfer function of the vanilla model. This method does have two limitations though: as long as the transfer function is held fixed, one is only sensitive to features that are not degenerate with the cosmological parameters that determine the transfer function (e.g., typically sharp features such as spikes or oscillations, but not broad distortions). It would be possible to circumvent this problem by also optimising the transfer function, but such a procedure would become computationally prohibitively expensive in the context of our analysis. Also, since we deconvolve binned angular power spectrum data, we are not sensitive to extremely high-frequency modulations of the primordial spectrum as reported in Ref. [39]. Nevertheless, the ∆χ 2 eff obtained with our method provides a good measure of the potential deviation from smoothness for realistic primordial spectra.

Numerical implementation
Applied to real data, the calculation of ∆χ 2 eff involves a maximisation of the likelihood using the reduced (i.e., τ = 0.089) vanilla model, and the application of the RLdeconvolution algorithm with the best-fit parameter transfer function.
Naturally, the simulated data sets should be treated in exactly the same way as the real data. Schematically, the procedure to be performed for each random realisation of WMAP data should hence read as follows: (i) Generate a random realisation (i) C sim of simulated WMAP temperature data.
(ii) Find minimal χ 2 eff,V for a fit of the vanilla model to (i) C sim .
(iii) Calculate the transfer function (i) T (k) for the best-fit parameters found in step (ii).
(iv) Apply the modified Richardson-Lucy deconvolution algorithm to (i) C sim , using the transfer function (i) T (k).
(v) Determine χ 2 eff,RL from the deconvolved primordial spectrum and calculate the improvement ∆χ 2 eff over the power-law fit.
To keep the task of simulating mock WMAP data at the level of generating angular power spectra (i) C sim and avoid having to generate random maps, we adopt a slightly simplified form of the WMAP likelihood function (described in Appendix A.1) throughout this paper. For details of simulating mock WMAP data (step (i)) we refer the reader to Appendix A.3. The minimisation of step (ii) turns out to be the computationally most expensive part of the algorithm. Deterministic, simplex-like minimisation routines like Numerical Recipes' amoeba [59] are not very reliable for our purpose since they run the risk of getting stuck in local minima. Random-walk based routines, such as simulated annealing are more suitable here. We choose a combination of the amebsa [59] routine to find a good starting point for a subsequent low-temperature run of the Metropolis-Hastings algorithm (based on the CosmoMC code [60]), which results in an estimate of the minimal χ 2 eff,V that is accurate to about 0.1. Nonetheless, this algorithm requires several thousand evaluations of CMB angular power spectra and the corresponding likelihoods, which would take a considerable time if one were to use a conventional Boltzmann code, such as CAMB [61]. We resort instead to the interpolation code PICO [62,63] for calculating the C s, which can considerably speed up this task, returning a reliable estimate for the minimal χ 2 eff,V within a few minutes on a regular desktop CPU.
Having found a good estimate of the maximum likelihood point, we then use CAMB to calculate the transfer function for the best-fit cosmological parameters in step (iii). Having applied the RL-deconvolution (step (iv)), giving us the optimised spectrum P R RL (k), we can finally compute the improvement ∆χ 2 eff the optimised spectrum would yield over the best power-law one (step (v)).

Real data
We commence by applying the RL-deconvolution method to the observed WMAP5 temperature data. The resulting best fit "optimised" primordial spectrum is shown in Fig. 1; its dominant features are qualitatively similar to those found in Ref. [48] for the 3-year WMAP data: a cutoff-like suppression at the largest scales accompanied by two prominent wiggles around k ∼ 0.002 Mpc −1 and k ∼ 0.07 Mpc −1 . The optimised spectrum gives χ 2 eff,RL = 1029.45, which marks an improvement of ∆χ 2 = 24.41 over the power-law best fit of χ 2 eff,V = 1053.86.

Simulated data
We have generated 2000 random realisations of simulated WMAP5 temperature data for multipoles 2 ≤ ≤ 1000 and applied the procedure of Section 2.3 to determine the ∆χ 2 eff for each of them. A histogram of the distribution of ∆χ 2 eff values is shown in Figure 2. Since the RL-algorithm is applied to binned data whereas ∆χ 2 eff is defined using un-binned data, for some simulated spectra the resulting ∆χ 2 eff can be negative, even after wavelet smoothing. In other words, in these cases the RLoptimised spectra yield no improvement over a power-law, and we set ∆χ 2 eff = 0. This has occurred in about 5% of our simulated spectra and is marked by the red part of the lowest-∆χ 2 eff -bin in the histogram. The vertical blue line in Figure 2 marks the observed value ∆χ 2 eff = 24.41. Of the 2000 mock data sets, 525 have ∆χ 2 eff > 24.41, corresponding to a p-value of ∼ 0.26. To ensure this number is not affected by our simplified likelihood-approximation, we calculated the observed ∆χ 2 eff using the full WMAP5 likelihood code (version v3p2), and found a difference of ∼ 1.
Additionally, we considered statistics which are less susceptible to low multipoles where the likelihood approximation is most inaccurate: ∆χ 2 eff evaluated between min and max = 1000 with min = 10, 20, and 30. The resulting p-values are 0.33, 0.14, and 0.31, respectively. This behaviour can be understood from the features observed in the real data: removing < 10 gets rid of a feature, leading to an increase in the p-value; taking away 10 ≤ < 20 removes a part of the data which does not show features, causing the p-value to drop, only to go up again when one removes the feature at 20 ≤ < 30.
In summary, the improvement in the fit from using an "optimised" primordial spectrum over a power-law spectrum is larger than the observed one in roughly 30% of our simulated data sets which assume an underlying smooth spectrum. We thus conclude that the null hypothesis cannot be rejected at high significance within the limitations of our chosen statistic, and the data are compatible with a smooth powerlaw primordial spectrum. We illustrate this qualitatively in Fig. 3, which shows that the reconstructed spectrum from the observed data does not have an unusual amount of "featureness" compared to the simulated data.

Conclusions
We have addressed the question whether the 5-year WMAP temperature anisotropy data are compatible with the assumption of an underlying smooth power-law primordial spectrum of curvature perturbations, or whether they show any indication for features or other unaccounted-for systematic effects. Assuming the true underlying primordial spectrum to be of power-law form, we generated 2000 simulated WMAP angular temperature power spectra, and estimated the amount of features by looking at the improvement to the likelihood gained from fitting an "optimised" primordial spectrum, obtained by deconvolving the data, instead of the usual power- For comparison we show the primordial spectrum reconstructed from WMAP5 data (red line) and our fiducial power-law spectrum (black line). It is evident that the spectrum reconstructed from real data does not have an unusual amount of features. The apparent feature at 0.05 Mpc −1 < k < 0.07 Mpc −1 is caused by the noise term becoming dominant at the corresponding multipoles in the binned WMAP data. For the last bin, 951 ≤ ≤ 1000, the averaged transfer function peaks at k ∼ 0.07 Mpc −1 . Beyond that the reconstructed primordial spectrum is no longer dominated by the peak of one data bin, but rather by the tails of the transfer functions of the last few bins; more structure here does not further improve the fit.
law. We found that 26% of all simulated data sets show a greater improvement than the one observed, which leads us to conclude that the features seen in the WMAP 5-year temperature data are not at all incompatible with the assumption of a smooth primordial spectrum. 1 However, we emphasise that our analysis does not disprove the existence of features -they are merely not strictly required by present data. Specifically, 'simple' theoretical models that predict similar features could conceivably remain favoured over a power-law. Additionally, since the spectrum reconstruction method used here is not sensitive to features below the binning scale of the WMAP data, we cannot rule out the possibility of extremely high-frequency features. We have demonstrated the feasibility of a full-scale frequentist hypothesis testing analysis in the search for deviations from smooth spectra though, and our method can just as well be applied to other statistics which allow for even more general shapes of the primordial spectrum.
Future analyses and new measurements may well reveal evidence for the existence of features, and a search for them will certainly remain a worthwhile endeavour, considering that any detection would profoundly impact our understanding of the physics of inflation and also precise estimation of the cosmological parameters [64]. Besides the imminent improvements to temperature anisotropy data from the Planck satellite [66], better measurements of the E-polarisation of the CMB will greatly enhance sensitivity to primordial features [52,67] -a dedicated mission like CMBPol [68] may prove particularly helpful here.

Acknowledgments
We thank Pedro Ferreira, Subir Sarkar and Yvonne Wong for interesting discussions and comments. JH thanks the Oxford Physics group for their hospitality during the initial stages of this work.

A. Appendix: Technical details A.1 Likelihood function
For reasons of simplicity we adopt here a form of the likelihood function similar to the one used by the WMAP team in their first data release [69]. For an input theoretical angular power spectrum C th , the likelihood of the data reads and a log-normal part where z th = ln(C th + N ) andẑ = ln(Ĉ + N ), withĈ the angular power spectrum estimated from observation. The curvature matrix Q is given by The off-diagonal terms induced by the sky cut r , the effective noise power spectrum N and the effective sky fraction f sky are supplied by the WMAP team [23] and available for download on the lambda web site 2 .
A.2 RL-algorithm for reconstructing the primordial spectrum The Richardson-Lucy (RL) algorithm was developed and is widely used in the context of image reconstruction in astronomy [70,71]. However, the method has also been successfully used in cosmology, to deproject the 3-D correlation function and power spectrum from the measured 2-D angular correlation and 2-D power spectrum [72,73]. The angular power spectrum, C , is a convolution of the initial power spectrum P R (k) generated in the early universe, with a radiative transfer kernel T (k) (the transfer function), that is determined by the values of the cosmological parameters. In our application, we solve the inverse problem of determining the primordial power spectrum, P R (k), from the measured angular power spectrum, C , using the relation In the above equation, the target measured angular power spectrum, C ≡ C obs , is the data given by observations, and the radiative transport kernel, encodes the response of the present multipoles of the CMB perturbed photon distribution function ∆ T (k i , η 0 ) to unit of power per logarithm interval of wavenumber, k, in the primordial perturbation spectrum. The kernel T (k) is completely fixed by the cosmological parameters of the base cosmological model. The kernel T (k) also includes the effect of geometrical projection from the three-dimensional wavenumber, k, to the harmonic multipole, l on the two-dimensional sphere.
Obtaining P R (k) from the measured C for a given T (k) is clearly a deconvolution problem. An important feature of the problem is that C obs , T (k) and P R (k) are all positive definite. We employ an improved RL method to solve the inverse problem for P R (k) in Eq. (A.6). The advantage of the RL method is that positivity of the recovered P R (k) is automatically ensured, given positive definite T (k) and C s. The RL method, readily derived from elementary probability theory of distributions [70], is an iterative method that can be neatly encoded into a simple recurrence relation. The power spectrum P R (i+1) (k) recovered after iteration (i + 1) is given by is the normalised kernel (on the space for all k wavenumbers),C obs is the normalised measured data (target) and C (i) is the angular power spectrum at the i th iteration obtained from C (i) = T (k) P R (i) (k) using the recovered power spectrum P R (i) (k). It is important to keep in mind that, due to the formulation in terms of conditional probability distributions, the RL method requires the kernel T (k), data C , and the target vector P R (k), all to be normalised at the beginning, whereP R (1) (k) is the normalised initial guess model of the primordial spectrum and the normalisation factor ζ k is defined by ζ k = T (k) . Note that Eq.(A.8) represents a modified form of the Richardson-Lucy algorithm specifically tailored for our purpose in this problem [45,47,48].
Following Refs. [45,47,48], we apply the RL-algorithm to a binned version of the (real and simulated) WMAP5 data. The resulting raw deconvolved spectrum has spurious oscillations and features on scales smaller than the bin size, arising largely due to the k space sampling and binning in space, which adversely impact the (unbinned) full likelihood of the data given the reconstructed spectrum. As shown in Ref. [47], a subsequent smoothing of the raw reconstructed spectrum with a discrete wavelet transform can lead to a significant improvement in the likelihood.

A.2.1 Discrete Wavelet Transform
Wavelet transforms provide a powerful tool for the analysis of transient and nonstationary data and is particularly useful in picking out characteristic variations at different resolutions or scales. This linear transform separates a data set in the form of low-pass or average coefficients, which reflect the average behaviour of the data, and wavelet or high-pass coefficients at different levels, which capture the variations at corresponding scales. As compared to Fourier or window Fourier transform, wavelets allow optimal "time-frequency" localisation in the real as well as in the Fourier domain. The vocabulary of DWT stems from applications in one dimensional timestream signal trains, but has found wide application in signal in other domains and dimensions. Specifically in our case, the signal being transformed is the primordial power spectrum, P R (k), a one dimensional function of wavenumber, k.
Wavelets are an orthonormal basis of small waves, with their variations primarily concentrated in a finite region, which make them ideal for analysing localised transient signals. Wavelets can be continuous or discrete. In the latter case, the basis elements are strictly finite in size, enabling them to achieve localisation, while disentangling characteristic variations at different frequencies [74]. For more details about DWT and its theoretical basis, see [47].
In this paper, we use the prescription of Ref. [48], and smooth the raw deconvolved spectrum with the following method: which can be split up into an original signal C obs and a contribution from experimental noise N obs . Theory, on the other hand, can only predict the average of this quantity over an ensemble of independent observations, C th = C obs . In the following, we shall describe how to generate independent realisations of simulated mock data (which replace theĈ in the associated mock likelihood function) from a theoretical input spectrum C fid . Since theĈ are actually defined as the difference between the observed spectrum and the expected noise power spectrum, C = C obs − N th , the corresponding simulated quantity is C sim − N th .
For the sake of clarity we commence with the case of an ideal (i.e., noise-free, full-sky) observation.

A.3.1 No experimental noise, full sky-coverage, no correlations
Even in the case of an ideal measurement of the CMB temperature anisotropies (N = 0), the observed angular power spectrum would still be subject to cosmic variance. Under the assumption of Gaussian fluctuations, the coefficients a m can be considered Gaussian random variables for the purpose of simulating data sets C sim . Isotropy dictates a m a * m = C fid δ δ mm , so different multipoles will be independent of each other. Hence, for each the simulated C sim can be constructed by drawing random numbers x from a χ 2 2 +1 -distribution:

A.3.2 Adding experimental noise
Realistic observations are of course subject to instrumental noise.
Assuming the noise to be isotropic and Gaussian, its effect on the distribution the a m are drawn from can be considered as a convolution with a Gaussian of width N th , and we have

A.3.3 Partial sky coverage and correlations
Even for nominally full-sky observations, some parts of the sky cannot be used to construct spectra and need to be cut out, since the signal is obscured by point sources and galactic foregrounds. This has two consequences on the recoveredĈ . First, using only part of the sky leads to a loss of information, resulting in a larger scatter ofĈ around the true underlying spectrum C fid . Second, since the sky cut breaks isotropy, the (pseudo)-Ĉ constructed from an incomplete map will no longer be uncorrelated [75]. Since the curvature matrix (A.4) appearing in the likelihood function scales like f 2 sky , the standard deviation of C sim − C fid should scale like f −1 sky . This increase in scatter around the fiducial model can be modelled by drawing the random variables x from a χ 2 f 2 sky (2 +1) -distribution. Since it is numerically much easier to generate χ 2 variates with an integer number of degrees of freedom we round f 2 sky (2 + 1) to the nearest integer and eventually rescale the resulting y ≡ x / x − 1 accordingly, to recover the correct standard deviation.
The so-generated y are uncorrelated random variables. For the simulated spectrum to have the same correlation properties as the real WMAP data, we take instead the vector (ỹ ) ≡ Chol(K )(y ), (A.14) where the lower triangular matrix Chol(K ) is the Cholesky decomposition of the correlation matrix K = δ + r . Finally, the simulated spectra are given by C sim − N th = C fid + C fid + N th ỹ . (A.15) The above procedure can rarely result in negative values of C sim + N sim at low multipoles, due to the rescaling of the y or due to the transformation of Eq. (A.15). Such data is unphysical and not compatible with the log-normal piece of the likelihood function and must therefore be discarded. As a result, the simulated data can be slightly biased; however, the bias is negligible and remains smaller than O(1%) of cosmic variance for all .