Light Sterile Neutrinos and Inflationary Freedom
We perform a cosmological analysis in which we allow the primordial power spectrum of scalar perturbations to assume a shape that is different from the usual power-law predicted by the simplest models of cosmological inflation. We parameterize the free primordial power spectrum with a “piecewise cubic Hermite interpolating polynomial” (PCHIP). We consider a 3+1 neutrino mixing model with a sterile neutrino having a mass at the eV scale, which can explain the anomalies observed in short-baseline neutrino oscillation experiments. We find that the freedom of the primordial power spectrum allows to reconcile the cosmological data with a fully thermalized sterile neutrino in the early Universe. Moreover, the cosmological analysis gives us some information on the shape of the primordial power spectrum, which presents a feature around the wavenumber .
[a]Department of Physics, University of Torino, Via P. Giuria 1, I–10125 Torino, Italy \affiliation[b]INFN, Sezione di Torino, Via P. Giuria 1, I–10125 Torino, Italy \affiliation[c]Dipartimento di Fisica e Astronomia “G. Galilei”, Università di Padova, and INFN, Sezione di Padova, Via F. Marzolo 8, I–35131 Padova, Italy
firstname.lastname@example.org, email@example.com, firstname.lastname@example.org
In typical analyses of cosmological data one of the main assumptions about the early Universe is the form of the primordial power-spectrum (PPS) of scalar fluctuations. The PPS is usually assumed to be a power-law (PL), as predicted by the simplest models of inflation (see Refs. [1, 2, 3]). However, if inflation is generated by a more complicated mechanism, the PPS can assume a different shape or present various features with respect to the power-law form (see Refs. [4, 5] and references therein). Since we cannot test directly the physics at the scale of cosmological inflation in order to check the correctness of the simplest inflationary models, any cosmological analysis performed assuming a power-law PPS can suffer of a biased constraint.
The cosmological observable that we can access is the late-time power spectrum of scalar perturbations, which is a convolution of the PPS and the transfer function, that can be calculated numerically as a function of a small number of cosmological parameters. Several experiments are designed to measure the late-time power spectrum at different redshifts (see Refs. [6, 7, 8]).
The physics of the transfer function is well understood and the experiments that measure the Cosmic Microwave Background (CMB) radiation give very strong constraints on the values of the cosmological parameters which determine the transfer function. The current most precise measurements of the angular power spectrum of the CMB are those of the Planck experiment  for the unpolarized data and those of the WMAP experiment  for the polarization spectra. However, the next Planck data release is expected to improve the current sensitivity on the unpolarized spectra and to include the new polarized spectra obtained by Planck.
On the other hand, since the inflationary scale cannot be directly tested, we can only try to reconstruct indirectly the PPS. In the literature there are several approaches for reconstructing a completely unknown PPS given the available experimental data. Among them we can list the “cosmic inversion” methods [11, 12, 13, 14, 15], maximum entropy deconvolution  and regularization methods like Richardson-Lucy iteration [17, 18, 19, 20], truncated singular value decomposition  and Tikhonov regularization .
The effects on cosmological parameter estimation of considering a PPS which is different from a power-law has been studied in several works: for example, the power-law PPS has been simply modified with the introduction of a running in the tilt of the power-law [23, 24, 25, 26], a running of the running , or a sharp cut-off in the power-law . Our main goal is to study how the freedom of the form of the PPS can affect the existing bounds on the presence in the early Universe of additional sterile neutrinos. In particular, we want to explore the impact of a light sterile neutrino with mass which has been thermalized by neutrino oscillations before neutrino decoupling at a temperature of the order of 1 MeV [28, 29]. Previous analyses of the cosmological data with a standard power-law PPS have found that the case of a fully thermalized sterile neutrino is quite disfavored [30, 31, 32, 33, 34, 35]. This result motivated the study of mechanisms which can suppress the thermalization of sterile neutrinos in the early Universe due to active-sterile oscillations before neutrino decoupling [28, 29]. Examples are a large lepton asymmetry [36, 37, 38, 39, 40], an enhanced background potential due to new interactions in the sterile sector [41, 42, 43, 44, 45, 46, 47], a larger cosmic expansion rate at the time of sterile neutrino production , and MeV dark matter annihilation .
Besides our main objective, which is to find out how the constraints on the presence in the early Universe of additional sterile neutrinos change if the PPS is not forced to be a power-law, we are also interested in obtaining information on the form of the PPS. With these aims, we considered a general form of the PPS that allows the presence of features without forcing a particular shape. In the literature several model-independent parameterizations have been used: for example, a free PPS can be described with wavelets [50, 51, 52, 53], principal components , top-hat bins without interpolation , power-law bins [56, 57], linear interpolation [58, 59, 60, 61, 62, 63, 64, 25], broken power-law [57, 65], and interpolating spline functions [66, 67, 68, 69, 70, 71, 26, 72]. We decided to follow part of the prescriptions of the interpolating spline form presented in Refs. [69, 71], improving the parametrization by using a “piecewise cubic Hermite interpolating polynomial” (PCHIP), which is described in Appendix A. This method allows us to avoid the spurious oscillating behavior that can appear between the nodes of interpolating splines.
This article is structured as follows: in Sec. 2 we introduce the neutrino mixing scheme, the cosmological model and the cosmological data used in the paper, in Sec. 3 and in Sec. 4 we discuss the results obtained from the analysis for the cosmological parameters and for the PPS respectively, and in Sec. 5 we present our conclusions.
2 Neutrino mixing scheme, cosmological model and data
In this Section we introduce the models and the datasets used in this paper. In Subsection 2.1 we present the neutrino mixing scheme, in Subsection 2.2 we introduce the cosmological model, and in Subsection 2.3 we present the cosmological data used in our analysis.
2.1 Neutrino mixing scheme
In this paper we consider the 3+1 neutrino mixing scheme, which is motivated by indications in favor of short-baseline neutrino oscillations found in the LSND experiment , in Gallium experiments [74, 75, 76, 77, 78] and in reactor experiments [79, 80, 81]. In this scheme, besides the three standard active flavor neutrinos , , , there is a sterile neutrino which does not interact through standard weak interactions. This sterile neutrino is a new particle beyond the Standard Model which cannot be detected directly in current experiments (see [82, 83, 84]).
The four flavor neutrinos , , , are superpositions of four massive neutrinos, , , , with respective masses , , , . The squared mass differences and (where ) generate the neutrino oscillations which have been observed in many solar, atmospheric and long-baseline experiments (see [85, 86, 87, 88]). An additional much larger squared mass difference is required in order to explain the indications in favor of short-baseline oscillations [89, 90, 91, 92, 93, 94, 95]. In the 3+1 scheme the three standard active flavor neutrinos , , are mainly composed of the three massive neutrinos , , , but they have a small component of in order to generate the observed short-baseline oscillations through . On the other hand, the sterile neutrino is mainly composed of the massive neutrino and in the following we use the common notation .
Since the case of a very light and almost degenerate , , at the eV scale is strongly disfavored by cosmological data (see Ref. ) we consider the case of and much lighter , , .
The combined analysis of cosmological data and short-baseline neutrino oscillation data is performed by using the posterior distribution of obtained from the analysis of SBL data  as a prior in the CosmoMC analysis of cosmological data [97, 31, 33, 34, 98]. As shown in Tab. 3 of Ref. , the best-fit value of obtained from short-baseline neutrino oscillation data is 1.27 eV and its 95.45% probability range () is between 0.97 and 1.42 eV.
2.2 Cosmological model
We used an extended flat CDM model to accommodate the presence of a sterile neutrino and inflationary freedom in the production of the primordial power spectra.
In the analysis with a power-law PPS we consider a flat CDM+ cosmological model with a total of eight parameters:
where and are the present-day physical CDM and baryon densities, the angular the sound horizon, the optical depth to reionisation, and and denote respectively the amplitude and spectral index of the initial scalar fluctuations at the pivot scale of 0.002 Mpc. , where  is the effective number of relativistic degrees of freedom before photon decoupling in the Standard Model with three massless neutrinos (see [100, 101]).
In contrast with previous analyses [33, 34, 98], we limit the allowed range of in the interval , assuming that the additional sterile neutrino cannot contribute to the relativistic energy density more than a standard active neutrino. This happens if sterile neutrinos are produced in the early Universe by neutrino oscillations before neutrino decoupling [28, 29].
We assume a flat prior for all the parameters in Eq. (1), except , for which we use a flat prior for only in the analyses which do not take into account the constraints from short-baseline neutrino oscillation data. In the analyses which take into account these constraints we use as prior for the posterior obtained from the analysis of SBL data presented in Ref. , as explained at the end of Subsection 2.1. We neglect the masses of the three light neutrinos , , , which are assumed to be much smaller than 1 eV.
In order to parameterize a free PPS we follow partially the prescriptions described in [69, 71], but instead of the cubic spline function we interpolate with a “piecewise cubic Hermite interpolating polynomial” (PCHIP) [102, 103], that is described only by the values of the PPS in a discrete number of nodes, as discussed in Appendix A. We used 12 nodes which span a wide range of values:
In the range , that is well constrained from the data , we choose equally spaced nodes in the logarithmic scale. The nodes and are used to parameterize a non-constant behaviour of the PPS outside this range and their position is chosen in order to have all the CosmoMC PPS evaluations inside the interval covered by our parametrization. The PCHIP PPS is described by
In the PCHIP PPS analysis we consider a flat CDM+ cosmological model with a total of 18 parameters:
where are the same as those in the set (1). We assume a flat prior on in the range .
2.3 Cosmological data sets
CMB (Cosmic Microwave Radiation): the temperature fluctuations power spectra provided by the Planck satellite  up to , by Atacama Cosmology Telescope (ACT)  and South Pole Telescope (SPT)  whose likelihoods cover the high multipole range, and , respectively. Concerning polarization we include the data of the Wilkinson Microwave Anisotropy Probe (WMAP) nine year data release .
LSS (Large Scale Structure): the matter power spectrum at four different redshifts , , and from the WiggleZ Dark Energy Survey .
: the Hubble parameter as obtained with the Hubble Space Telescope (HST) , which acts as a prior on the derived cosmological parameter .
PSZ: The Planck Sunayev Zel’Dovich catalogue  contains 189 galaxy clusters identified through the Sunayev Zel’Dovich effect. The number counts allows to compute the cluster mass function, which is related to a combination of and : . This result contributes as an additional in our analysis.
CFHTLenS: the 2D cosmic shear correlation function as determined by the Canada-France Hawaii Telescope Lensing Survey (CFHTLenS) [113, 114] through the measurements of redshifts and shapes of 4.2 million galaxies spanning the range . The weak gravitational lensing signal extracted from these measurements constrains a combination of the total matter density and the standard deviation of the amplitude of the matter density fluctuations on a sphere of radius : . This result is incorporated in our analysis following the same prescription used for PSZ.
In the following we denote the analyses of all these cosmological data alone as “COSMO” and those which include also the short-baseline neutrino oscillation prior as “COSMO+SBL”.
|; No limit; No limit||; ;|
||; ; No limit|
|; No limit; No limit||; ; No limit|
||; No limit|
|; ; No limit||; ; No limit|
||; ;||; ;|
3 Cosmological Parameters and Sterile Neutrinos
The results of our COSMO and COSMO+SBL analyses are presented in Tab. 1 for the standard case of a power-law PPS and in Tab. 2 for the free PPS with the PCHIP parameterization. In the upper part of the tables we listed the common parameters of the CDM model, in the central part we listed the neutrino parameters and , while the lower part concerns the parameters used to parameterize the PPS: and for the power-law PPS and for the PCHIP PPS. The constraints on the PPS parameters are discussed in the next section. In this section we discuss first the results relative to the parameters in the upper part of the tables, , , and , and then the results relative to the parameters in the central part of the tables, and .
The bounds on the parameters of the CDM model change slightly when more freedom is admitted for the PPS. Comparing Tabs. 1 and 2, one can see that the limits on the parameters of the CDM model are slightly weakened in the PCHIP PPS case and for some parameters there is also a small shift in the marginalized best-fit value. In all the cases in which this happens, the marginalized best-fit values move inside the uncertainties. The freedom of the form of the PPS affects the COSMO results more than the COSMO+SBL results: in the former case the and best values change by about , while a smaller shift is obtained for . On the other hand, in the COSMO+SBL analysis all the shifts are much smaller than the uncertainties.
Figure 1 shows the marginalized , and allowed intervals for and that we obtained in the COSMO(PL) and COSMO(PCHIP) analyses, without the SBL prior. Figure 3 shows the corresponding marginalized , and allowed regions in the – plane. We can notice some major changes in the allowed values of both and in the PCHIP PPS case with respect to the power-law PPS case. With a power-law PPS the best-fit value of is around 0.4, whereas with the PCHIP PPS it is at , that is the upper limit for assumed in the analysis. The reason of this behavior is that the effects of the presence of additional relativistic energy in the primordial universe can be compensated by an increase of the PCHIP PPS at large . As a result, the marginalized posterior for is increased in the region towards , in correspondence with higher values in the PCHIP PPS for .
Without the SBL constraint on , the different preferences for the value of in the power-law and PCHIP PPS analyses correspond to different allowed intervals for . As shown in Fig. 1, although in both cases the best-fit value of is near 0.6 eV, the intermediate preferred region for in the power-law PPS analysis gives for an upper limit of about 2 eV at , whereas the large preferred values for in the PCHIP PPS analysis gives a tighter upper limit of about 1.5 eV at .
The SBL prior on the sterile neutrino mass puts a constraint so strong that in practice the value of this parameter does not depend on the inclusion or not of the freedom of the PPS. In fact, the limits in Tabs. 1 and 2 are similar in the power-law PPS and PCHIP PPS analyses. This can be seen also from the marginalized allowed intervals of in Fig. 1, comparing the COSMO+SBL(PL) and COSMO+SBL(PCHIP) allowed intervals.
A major difference occurs, instead, in the limits for , because the effects of the presence of additional relativistic energy in the primordial universe can be compensated by an increase in the PCHIP PPS at large , as in the case without the SBL constraint on . As shown in Fig. 1, the best-fit and upper limits on in the COSMO+SBL(PL) and COSMO+SBL(PCHIP) are different. In particular, in the COSMO+SBL(PCHIP) the upper limit on allows the presence of a fully thermalized sterile neutrino compatible with the SBL constraint on .
Figure 3 shows the contour plots of the marginalized , and regions in the – plane that we obtained in the COSMO+SBL(PL) and COSMO+SBL(PCHIP) analyses. The allowed regions in the left panel are similar111 The only difference is that the analysis in Ref.  took into account also the BICEP2 data on the B-mode polarization of the CMB . to those obtained in Ref.  with a standard power-law PPS. One can see that in this case a fully thermalized sterile neutrino is quite disfavored. On the other hand, from the right panel one can see that in the PCHIP PPS analysis a fully thermalized sterile neutrino with a mass just below and with is even inside the region. This means that a fully thermalized sterile neutrino can be accommodated in the cosmological model if the PPS is not forced to be described by a power-law.
4 Best-fitting Primordial Power Spectrum
The results of our PCHIP PPS analyses without and with the SBL prior on give interesting information on the shape of the PPS.
The marginalized posterior limits for the values in Eq. (3) are listed in Tab. 2. One can see that the least constrained nodes are the first and last, in and , for which there are only upper limits on the corresponding . This was expected, because there are no data at the edges of the wide interval of that we have considered. For these two extreme nodes the marginalized posterior is peaked on the lowest value that we allowed in the fit (0.01).
On the contrary, the nodes from 5 to 10 are well constrained, at the level of a few percent at . Considering the nodes from 2 to 4, one can see that the second node has preferred values higher than 1 within , the third node value is more than below 1 (around 0.6), the fourth node value is again higher than 1 at more than . This implies that the PPS that we obtained from the MCMC has a significant dip around and a less significant bump around .
To help the reader to visualize this feature, we present in Fig. 4 a comparison of the best-fitting PPS 222 We consider as the best-fitting PPS that which corresponds to the lower value of , where is the marginalized posterior probability in the space of the parameters . However, one must take into account that in a parameter space with a large number of dimensions the MCMC is not expected to explore well the region near the true global best-fit corresponding to . In fact, the points are distributed mainly in a region where . Therefore, the PPS that we consider as best-fitting can be different from the true best-fitting PPS in the intervals of which are not well constrained by the data. in the power-law parametrization and in the PCHIP parametrization, without and with the SBL constraint. One can see that the best-fitting PCHIP curves with and without the SBL prior are significantly different only for . The dip around and the bump around are clearly seen in the PCHIP parametrization.
From Fig. 4 one can also see that the PCHIP parametrization has an approximate power-law behavior between about and .
Another helpful way to visualize the behaviour of the PPS obtained in the analyses without and with the SBL prior with the PCHIP parametrization is presented in Fig. 5, which shows the , and bands obtained by marginalizing the posterior distribution for each value of the wavenumber in a fine grid. The two plots in Fig. 5 show a well collimated band corresponding to the region in which the power-law gives a good approximation of the PCHIP PPS, between about and . Moreover, the dip at is well visible, as well as the bump at . On the other hand, the widths of the bands diverge for low and high values of , where there are large uncertainties.
The major features that we have noticed in the reconstructed PPS are in agreement with those found in Ref. , in which the scalar PPS has been reconstructed with a totally different technique, the Richardson-Lucy iteration algorithm, using the transfer function corresponding to the Planck 2013 best-fit for the CDM model. Apart for the suppression that they found around and the features at higher , the main differences with respect to the power-law PPS are the same that we found in our analysis. According to the authors of Ref. , these major features are related to the low- spectrum of the temperature perturbations measured by the Planck experiment, that obtained a dip in the power around and a slight excess around .
Although the parametrization with a natural cubic spline is noisy due to spurious oscillations between the nodes, also in Fig. 8 of Ref.  it is possible to guess the presence of a dip around , especially in the curves with more than 20 nodes. However, our parametrization is much cleaner and permits a better visualization of these features.
In this work we analyzed the effects of a free form of the primordial power-spectrum of scalar fluctuations, which is not constrained to the usual power-law form that is predicted by the simplest models of inflation (see Refs. [1, 2, 3]). This freedom in the PPS could arise from a more complicate inflationary mechanism (see Refs. [4, 5] and references therein).
We parameterized the PPS with a “piecewise cubic Hermite interpolating polynomial” (PCHIP) described in details in Appendix A. Our PCHIP parameterization of the PPS depends from the values of the PPS in twelve nodes (given in Eq. (2)) which cover a wide range of values of the wavenumber . We choose the PCHIP method in order to avoid spurious oscillations of the interpolated function between the nodes that can be obtained with spline interpolations (see Refs. [69, 71]).
We performed an analysis of cosmological data in which only the primordial spectrum of scalar perturbations is considered, neglecting the controversial [96, 115] data on the B-mode polarization of the CMB  which would require to take into account also the primordial spectrum of tensor perturbations. We used the most precise CMB measurements together with low-redshift measurements of the Hubble parameter, the galaxy distribution and the matter distribution in the Universe (see Section 2.3).
We found that the freedom of the form of the PPS does not affect significantly the fitted values of the parameters in the CDM model, while the results concerning the existence of a sterile neutrino in the early Universe can change drastically. If we do not impose any prior on the sterile neutrino mass from the results of short-baseline oscillation experiments (see Section 2.1), a larger value for the sterile neutrino contribution to the effective number of relativistic degrees of freedom before photon decoupling is preferred in the PCHIP PPS parameterization with respect to the standard power-law parameterization. The marginalized best fit of is moved towards one, which corresponds to a fully thermalized sterile neutrino. This shift corresponds to a tightening of the cosmological preferred values for .
In the analysis with a prior on obtained from the fit of short-baseline oscillation experiments , the freedom of the PCHIP PPS affects only the bound on , because the allowed range of is strongly constrained by the SBL prior. We found that a free form of the PPS allows the existence in the early Universe of a fully thermalized sterile neutrino with a mass of about 1 eV [28, 29]. This possibility is quite disfavored by the analysis of cosmological data with a power-law PPS [30, 31, 32, 33, 34, 35]. Hence, the freedom of the PPS allows us to reconcile the cosmological data with short-baseline neutrino oscillations without the need of an additional mechanism which suppresses the thermalization of the sterile neutrino [36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49].
We obtained also some interesting information on the form of the PPS, whose behavior is well constrained by the cosmological analysis for . In particular, we have shown that in the range the PPS can be approximated with a power-law and the values of the PPS in the nodes of the PCHIP parameterization lying in this range of have only a few-percent uncertainty. The PPS in the range presents a clear dip at , with a statistical significance of more than 2, and a small bump at , with a statistical significance of about . These features of the PPS are in agreement with those found in Ref.  with a completely different method.
In the future the analysis presented in this work could be repeated with the inclusion of a parametrization for the primordial spectrum of tensor perturbations when improved data on the B-mode polarization of the CMB will be available. This will allow us to study with more precision the few relics of cosmological inflation that we can access.
We would like to thank M. Archidiacono, E. di Valentino, N. Fornengo, S. Hannestad, A. Melchiorri, Y.F Li and H.W. Long for stimulating discussions and fruitful collaboration in previous works. This work is supported by the research grant Theoretical Astroparticle Physics number 2012CPPYP7 under the program PRIN 2012 funded by the Ministero dell’Istruzione, Università e della Ricerca (MIUR).
Appendix A Pchip Parametrization of the Primordial Power Spectrum
In this work we parameterized the PPS with a “piecewise cubic Hermite interpolating polynomial” (PCHIP) [102, 103]. We decided to adopt this interpolating function in order to avoid spurious oscillations of the interpolating function between the nodes which is often obtained in spline interpolations. This problem occurs because a natural cubic spline requires the values of the function, the first and the second derivatives to be continuous in the nodes .
The PCHIP function, instead, is constructed in order to preserve the shape of the set of points to be interpolated. This is achieved with a modification of the “monotone piecewise cubic interpolation”  which can accommodate non-monotone functions and preserves the local monotonicity.
Let us consider a function with known values in N nodes , with . A piecewise cubic interpolation is performed with cubic functions between the nodes. The determination of these cubic functions requires the determination of coefficients. Besides the constraints obtained by requiring that the initial and final point of each cubic function match the known values of the original function in the corresponding nodes, one needs a prescription for the other necessary constraints. In the case of a natural cubic spline interpolation one gets constraints by requiring the continuity of the first and second derivatives in the nodes and the remaining two constraints are obtained by requiring that the second derivatives in the first and last nodes vanish. The drawback of this method is that the interpolating curve is determined by a set of linear equations without any local control. In fact, all the interpolating curve is affected by the change of a single point.
Local control of the interpolating curve can be achieved by relaxing the requirement of continuity of the second derivatives in the nodes and using the resulting freedom to adjust the first derivatives with a local prescription. In order to see how it can be done, it is convenient to write the cubic interpolating polynomial between the nodes and in the Hermite form
where and . Here and are the values of the derivatives in the two nodes. In the PCHIP method the derivatives are chosen in order to preserve the local monotonicity of the interpolated points. This is done by considering the relative differences
The PCHIP prescription is:
If and have opposite signs, then is a discrete local minimum or maximum and .
If and have the same sign, then is determined by the weighted harmonic mean
with and .
The derivatives in the first and last nodes are determined by a shape-preserving prescription based on a quadratic fit of three points. For we consider the three points , , . The derivative in of the parabola which passes through these three points is given by
The shape-preserving prescription for is:
If the signs of and are different, then .
If the signs of and are different and , then .
For one must replace and .
We fit the power spectrum with Eq. (3), in which the function is calculated with the PCHIP prescription in the logarithmic scale of :
A comparison between the natural cubic spline and the PCHIP interpolations of the PPS is presented in Fig. 6. We choose the same nodes positions that we used for the PPS parametrization in our cosmological analysis and we choose the values of the function in the nodes in order to show the difference between the natural cubic spline and the PCHIP interpolations. One can see that the PCHIP interpolation can reproduce the shape of the points without adding the spurious features between the points that are clearly visible in the natural cubic spline interpolation.
-  D. H. Lyth and A. Riotto, “Particle physics models of inflation and the cosmological density perturbation,” Phys.Rept. 314 (1999) 1–146, hep-ph/9807278.
-  B. A. Bassett, S. Tsujikawa, and D. Wands, “Inflation dynamics and reheating,” Rev.Mod.Phys. 78 (2006) 537–589, astro-ph/0507632.
-  D. Baumann and H. V. Peiris, “Cosmological Inflation: Theory and Observations,” Adv.Sci.Lett. 2 (2009) 105–120, 0810.3022 [astro-ph].
-  J. Martin, C. Ringeval, and V. Vennin, “Encyclopædia Inflationaris,” Phys.Dark Univ. (2014) , 1303.3787 [astro-ph.CO].
-  N. Kitazawa and A. Sagnotti, “Pre-inflationary clues from String Theory?,” JCAP 1404 (2014) 017, 1402.1418 [hep-th].
-  SDSS Collaboration, M. Tegmark et al., “Cosmological parameters from SDSS and WMAP,” Phys.Rev. D69 (2004) 103501, astro-ph/0310723.
-  D. Parkinson, S. Riemer-Sorensen, C. Blake, G. B. Poole, T. M. Davis, et al., “The WiggleZ Dark Energy Survey: Final data release and cosmological results,” Phys.Rev. D86 (2012) 103518, 1210.2130 [astro-ph.CO].
-  SDSS Collaboration Collaboration, C. P. Ahn et al., “The Ninth Data Release of the Sloan Digital Sky Survey: First Spectroscopic Data from the SDSS-III Baryon Oscillation Spectroscopic Survey,” Astrophys.J.Suppl. 203 (2012) 21, 1207.7137 [astro-ph.IM].
-  Planck Collaboration, P. Ade et al., “Planck 2013 results. I. Overview of products and scientific results,” Astron.Astrophys. 571 (2014) A1, 1303.5062 [astro-ph].
-  WMAP Collaboration, C. Bennett et al., “Nine-Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Final Maps and Results,” Astrophys.J.Suppl. 208 (2013) 20, 1212.5225 [astro-ph.CO].
-  M. Matsumiya, M. Sasaki, and J. Yokoyama, “Cosmic inversion: Reconstructing primordial spectrum from CMB anisotropy,” Phys.Rev. D65 (2002) 083007, astro-ph/0111549 [astro-ph].
-  M. Matsumiya, M. Sasaki, and J. Yokoyama, “Cosmic inversion. 2. An iterative method for reproducing the primordial spectrum from the CMB data,” JCAP 0302 (2003) 003, astro-ph/0210365 [astro-ph].
-  N. Kogo, M. Matsumiya, M. Sasaki, and J. Yokoyama, “Reconstructing the primordial spectrum from WMAP data by the cosmic inversion method,” Astrophys.J. 607 (2004) 32–39, astro-ph/0309662.
-  N. Kogo, M. Sasaki, and J. Yokoyama, “Constraining cosmological parameters by the cosmic inversion method,” Prog.Theor.Phys. 114 (2005) 555–572, astro-ph/0504471.
-  R. Nagata and J. Yokoyama, “Reconstruction of the primordial fluctuation spectrum from the five-year WMAP data by the cosmic inversion method with band-power decorrelation analysis,” Phys.Rev. D78 (2008) 123002, 0809.4537 [astro-ph].
-  G. Goswami and J. Prasad, “Maximum Entropy deconvolution of Primordial Power Spectrum,” Phys.Rev. D88 no. 2, (2013) 023522, 1303.4747 [astro-ph.CO].
-  A. Shafieloo and T. Souradeep, “Primordial power spectrum from WMAP,” Phys.Rev. D70 (2004) 043523, astro-ph/0312174.
-  G. Nicholson and C. R. Contaldi, “Reconstruction of the Primordial Power Spectrum using Temperature and Polarisation Data from Multiple Experiments,” JCAP 0907 (2009) 011, 0903.1106 [astro-ph.CO].
-  D. K. Hazra, A. Shafieloo, and T. Souradeep, “Cosmological parameter estimation with free-form primordial power spectrum,” Phys.Rev. D87 no. 12, (2013) 123528, 1303.5336 [astro-ph.CO].
-  D. K. Hazra, A. Shafieloo, and T. Souradeep, “Primordial power spectrum from Planck,” JCAP 1411 (2014) 011, arXiv:1406.4827 [astro-ph.CO].
-  G. Nicholson, C. R. Contaldi, and P. Paykari, “Reconstruction of the Primordial Power Spectrum by Direct Inversion,” JCAP 1001 (2010) 016, 0909.5092 [astro-ph.CO].
-  P. Hunt and S. Sarkar, “Reconstruction of the primordial power spectrum of curvature perturbations using multiple data sets,” JCAP 1401 (2014) 025, 1308.2317 [astro-ph.CO].
-  M. Archidiacono, E. Giusarma, A. Melchiorri, and O. Mena, “Dark Radiation in extended cosmological scenarios,” Phys.Rev. D86 (2012) 043509, 1206.0109 [astro-ph.CO].
-  G.-B. Zhao, S. Saito, W. J. Percival, A. J. Ross, F. Montesano, et al., “The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: weighing the neutrino mass using the galaxy power spectrum of the CMASS sample,” Mon.Not.Roy.Astron.Soc. 436 (2013) 2038–2053, 1211.3741 [astro-ph.CO].
-  J. A. Vazquez, M. Bridges, Y.-Z. Ma, and M. Hobson, “Constraints on the tensor-to-scalar ratio for non-power-law models,” JCAP 1308 (2013) 001, 1303.4014 [astro-ph.CO].
-  K. N. Abazajian, G. Aslanyan, R. Easther, and L. C. Price, “The Knotted Sky II: Does BICEP2 require a nontrivial primordial power spectrum?,” JCAP 1408 (2014) 053, 1403.5922 [astro-ph].
-  C. Cheng, Q.-G. Huang, and W. Zhao, “Constraints on the extensions to the base CDM model from BICEP2, Planck and WMAP,” Sci.China Phys.Mech.Astron. 57 (2014) 1460–1465, 1404.3467 [astro-ph.CO].
-  A. Dolgov and F. Villante, “BBN bounds on active sterile neutrino mixing,” Nucl.Phys. B679 (2004) 261–298, hep-ph/0308083.
-  M. Cirelli, G. Marandella, A. Strumia, and F. Vissani, “Probing oscillations into sterile neutrinos with cosmology, astrophysics and experiments,” Nucl.Phys. B708 (2005) 215–267, hep-ph/0403158.
-  E. Di Valentino, A. Melchiorri, and O. Mena, “Dark radiation sterile neutrino candidates after Planck data,” JCAP 1311 (2013) 018, 1304.5981 [astro-ph].
-  M. Archidiacono, N. Fornengo, C. Giunti, S. Hannestad, and A. Melchiorri, “Sterile neutrinos: Cosmology versus short-baseline experiments,” Phys.Rev. D87 (2013) 125034, 1302.6720 [astro-ph].
-  A. Mirizzi, G. Mangano, N. Saviano, E. Borriello, C. Giunti, et al., “The strongest bounds on active-sterile neutrino mixing after Planck data,” Phys.Lett. B726 (2013) 8–14, 1303.5368 [astro-ph].
-  S. Gariazzo, C. Giunti, and M. Laveder, “Light Sterile Neutrinos in Cosmology and Short-Baseline Oscillation Experiments,” JHEP 1311 (2013) 211, 1309.3192 [hep-ph].
-  M. Archidiacono, N. Fornengo, S. Gariazzo, C. Giunti, S. Hannestad, et al., “Light sterile neutrinos after BICEP-2,” JCAP 1406 (2014) 031, 1404.1794 [astro-ph.CO].
-  J. Bergström, M. Gonzalez-Garcia, V. Niro, and J. Salvado, “Statistical tests of sterile neutrinos using cosmology and short-baseline data,” JHEP 1410 (2014) 104, 1407.3806 [hep-ph].
-  Y.-Z. Chu and M. Cirelli, “Sterile neutrinos, lepton asymmetries, primordial elements: How much of each?,” Phys.Rev. D74 (2006) 085015, astro-ph/0608206.
-  S. Hannestad, I. Tamborra, and T. Tram, “Thermalisation of light sterile neutrinos in the early universe,” JCAP 1207 (2012) 025, 1204.5861 [astro-ph].
-  A. Mirizzi, N. Saviano, G. Miele, and P. D. Serpico, “Light sterile neutrino production in the early universe with dynamical neutrino asymmetries,” Phys.Rev. D86 (2012) 053009, 1206.1046 [hep-ph].
-  N. Saviano, A. Mirizzi, O. Pisanti, P. D. Serpico, G. Mangano, et al., “Multi-momentum and multi-flavour active-sterile neutrino oscillations in the early universe: role of neutrino asymmetries and effects on nucleosynthesis,” Phys.Rev. D87 (2013) 073006, 1302.1200 [astro-ph].
-  S. Hannestad, R. S. Hansen, and T. Tram, “Can active-sterile neutrino oscillations lead to chaotic behavior of the cosmological lepton asymmetry?,” JCAP 1304 (2013) 032, 1302.7279 [astro-ph].
-  S. Hannestad, R. S. Hansen, and T. Tram, “How Self-Interactions can Reconcile Sterile Neutrinos with Cosmology,” Phys.Rev.Lett. 112 (2014) 031802, 1310.5926 [astro-ph].
-  B. Dasgupta and J. Kopp, “Cosmologically Safe eV-Scale Sterile Neutrinos and Improved Dark Matter Structure,” Phys.Rev.Lett. 112 (2014) 031803, 1310.6337 [hep-ph].
-  T. Bringmann, J. Hasenkamp, and J. Kersten, “Tight bonds between sterile neutrinos and dark matter,” JCAP 1407 (2014) 042, 1312.4947 [hep-ph].
-  P. Ko and Y. Tang, “MDM: A Model for Sterile Neutrino and Dark Matter Reconciles Cosmological and Neutrino Oscillation Data after BICEP2,” Phys.Lett. B739 (2014) 62–67, 1404.0236 [hep-ph].
-  M. Archidiacono, S. Hannestad, R. S. Hansen, and T. Tram, “Cosmology with self-interacting sterile neutrinos and dark matter - A pseudoscalar model,” 1404.5915 [astro-ph.CO].
-  N. Saviano, O. Pisanti, G. Mangano, and A. Mirizzi, “Unveiling secret interactions among sterile neutrinos with big-bang nucleosynthesis,” Phys.Rev. D90 (2014) 113009, 1409.1680 [astro-ph].
-  A. Mirizzi, G. Mangano, O. Pisanti, and N. Saviano, “Collisional production of sterile neutrinos via secret interactions and cosmological implications,” Phys.Rev. D91 (2015) 025019, 1410.1385 [hep-ph].
-  T. Rehagen and G. B. Gelmini, “Effects of kination and scalar-tensor cosmologies on sterile neutrinos,” JCAP 1406 (2014) 044, 1402.0607 [hep-ph].
-  C. M. Ho and R. J. Scherrer, “Sterile Neutrinos and Light Dark Matter Save Each Other,” Phys.Rev. D87 (2013) 065016, 1212.1689 [hep-ph].
-  P. Mukherjee and Y. Wang, “Wavelet band powers of the primordial power spectrum from CMB data,” Astrophys.J. 593 (2003) 38, astro-ph/0301058 [astro-ph].
-  P. Mukherjee and Y. Wang, “Direct wavelet expansion of the primordial power spectrum: Results from pre-MAP CMB data,” Astrophys.J. 598 (2003) 779–784, astro-ph/0301562.
-  P. Mukherjee and Y. Wang, “Model-independent reconstruction of the primordial power spectrum from WMAP data,” Astrophys.J. 599 (2003) 1–6, astro-ph/0303211.
-  P. Mukherjee and Y. Wang, “Primordial power spectrum reconstruction,” JCAP 0512 (2005) 007, astro-ph/0502136.
-  S. M. Leach, “Measuring the primordial power spectrum: Principal component analysis of the cosmic microwave background,” Mon.Not.Roy.Astron.Soc. 372 (2006) 646–654, astro-ph/0506390.
-  Y. Wang, D. N. Spergel, and M. A. Strauss, “Cosmology in the next millennium: Combining MAP and SDSS data to constrain inflationary models,” Astrophys.J. 510 (1999) 20, astro-ph/9802231 [astro-ph].
-  S. Hannestad, “Reconstructing the inflationary power spectrum from CMBR data,” Phys.Rev. D63 (2001) 043009, astro-ph/0009296 [astro-ph].
-  D. K. Hazra, A. Shafieloo, and G. F. Smoot, “Reconstruction of broad features in the primordial spectrum and inflaton potential from Planck,” JCAP 1312 (2013) 035, 1310.3038 [astro-ph].
-  Y. Wang and G. Mathews, “A Measurement of the primordial power spectrum from Maxima and Boomerang data,” Astrophys.J. 573 (2002) 1, astro-ph/0011351 [astro-ph].
-  S. Bridle, A. Lewis, J. Weller, and G. Efstathiou, “Reconstructing the primordial power spectrum,” Mon.Not.Roy.Astron.Soc. 342 (2003) L72, astro-ph/0302306.
-  S. Hannestad, “Reconstructing the primordial power spectrum - A New algorithm,” JCAP 0404 (2004) 002, astro-ph/0311491.
-  M. Bridges, A. Lasenby, and M. Hobson, “A bayesian analysis of the primordial power spectrum,” Mon.Not.Roy.Astron.Soc. 369 (2006) 1123–1130, astro-ph/0511573.
-  WMAP Collaboration, D. Spergel et al., “Wilkinson Microwave Anisotropy Probe (WMAP) three year results: implications for cosmology,” Astrophys.J.Suppl. 170 (2007) 377, astro-ph/0603449.
-  M. Bridges, A. Lasenby, and M. Hobson, “WMAP 3-year primordial power spectrum,” Mon.Not.Roy.Astron.Soc. 381 (2007) 68–74, astro-ph/0607404.
-  M. Bridges, F. Feroz, M. Hobson, and A. Lasenby, “Bayesian optimal reconstruction of the primordial power spectrum,” Mon.Not.Roy.Astron.Soc. 400 (2009) 1075–1084, 0812.3541 [astro-ph].
-  D. K. Hazra, A. Shafieloo, G. F. Smoot, and A. A. Starobinsky, “Ruling out the power-law form of the scalar primordial spectrum,” JCAP 1406 (2014) 061, 1403.7786 [astro-ph].
-  C. Sealfon, L. Verde, and R. Jimenez, “Smoothing spline primordial power spectrum reconstruction,” Phys.Rev. D72 (2005) 103520, astro-ph/0506707 [astro-ph].
-  L. Verde and H. V. Peiris, “On Minimally-Parametric Primordial Power Spectrum Reconstruction and the Evidence for a Red Tilt,” JCAP 0807 (2008) 009, 0802.1219 [astro-ph].
-  H. V. Peiris and L. Verde, “The Shape of the Primordial Power Spectrum: A Last Stand Before Planck,” Phys.Rev. D81 (2010) 021302, 0912.0268 [astro-ph.CO].
-  R. Hlozek, J. Dunkley, G. Addison, J. W. Appel, J. R. Bond, et al., “The Atacama Cosmology Telescope: a measurement of the primordial power spectrum,” Astrophys.J. 749 (2012) 90, 1105.4887 [astro-ph.CO].
-  C. Gauthier and M. Bucher, “Reconstructing the primordial power spectrum from the CMB,” JCAP 1210 (2012) 050, 1209.2147 [astro-ph.CO].
-  R. de Putter, E. V. Linder, and A. Mishra, “Inflationary Freedom and Cosmological Neutrino Constraints,” Phys.Rev. D89 (2014) 103502, 1401.7022 [astro-ph].
-  B. Hu, J.-W. Hu, Z.-K. Guo, and R.-G. Cai, “Reconstruction of the primordial power spectra with Planck and BICEP2 data,” Phys.Rev. D90 (2014) 023544, 1404.3690 [astro-ph].
-  LSND Collaboration, A. Aguilar-Arevalo et al., “Evidence for neutrino oscillations from the observation of anti-neutrino(electron) appearance in a anti-neutrino(muon) beam,” Phys.Rev. D64 (2001) 112007, hep-ex/0104049.
-  SAGE Collaboration, J. Abdurashitov, V. Gavrin, S. Girin, V. Gorbachev, P. Gurkina, et al., “Measurement of the response of a Ga solar neutrino experiment to neutrinos from an Ar-37 source,” Phys.Rev. C73 (2006) 045805, nucl-ex/0512041.
-  M. Laveder, “Unbound neutrino roadmaps,” Nucl.Phys.Proc.Suppl. 168 (2007) 344–346.
-  C. Giunti and M. Laveder, “Short-Baseline Active-Sterile Neutrino Oscillations?,” Mod.Phys.Lett. A22 (2007) 2499–2509, hep-ph/0610352.
-  M. A. Acero, C. Giunti, and M. Laveder, “Limits on nu(e) and anti-nu(e) disappearance from Gallium and reactor experiments,” Phys.Rev. D78 (2008) 073009, 0711.4222 [hep-ph].
-  C. Giunti and M. Laveder, “Statistical Significance of the Gallium Anomaly,” Phys.Rev. C83 (2011) 065504, 1006.3244 [hep-ph].
-  T. Mueller, D. Lhuillier, M. Fallot, A. Letourneau, S. Cormon, et al., “Improved Predictions of Reactor Antineutrino Spectra,” Phys.Rev. C83 (2011) 054615, 1101.2663 [hep-ex].
-  G. Mention, M. Fechner, T. Lasserre, T. Mueller, D. Lhuillier, et al., “The Reactor Antineutrino Anomaly,” Phys.Rev. D83 (2011) 073006, 1101.2755 [hep-ex].
-  P. Huber, “On the determination of anti-neutrino spectra from nuclear reactors,” Phys.Rev. C84 (2011) 024617, 1106.0687 [hep-ph].
-  S. M. Bilenky, C. Giunti, and W. Grimus, “Phenomenology of neutrino oscillations,” Prog.Part.Nucl.Phys. 43 (1999) 1–86, hep-ph/9812360.
-  M. Gonzalez-Garcia and M. Maltoni, “Phenomenology with Massive Neutrinos,” Phys.Rept. 460 (2008) 1–129, 0704.1800 [hep-ph].
-  K. Abazajian, M. Acero, S. Agarwalla, A. Aguilar-Arevalo, C. Albright, et al., “Light Sterile Neutrinos: A White Paper,” 1204.5379 [hep-ph].
-  M. Gonzalez-Garcia, M. Maltoni, J. Salvado, and T. Schwetz, “Global fit to three neutrino mixing: critical look at present precision,” JHEP 1212 (2012) 123, 1209.3023 [hep-ph].
-  G. Bellini, L. Ludhova, G. Ranucci, and F. Villante, “Neutrino oscillations,” Adv.High Energy Phys. 2014 (2014) 191960, 1310.7858 [hep-ph].
-  F. Capozzi, G. Fogli, E. Lisi, A. Marrone, D. Montanino, et al., “Status of three-neutrino oscillation parameters, circa 2013,” Phys.Rev. D89 (2014) 093018, 1312.2878 [hep-ph].
-  M. Gonzalez-Garcia, M. Maltoni, and T. Schwetz. NuFIT, http://www.nu-fit.org/.
-  J. Kopp, M. Maltoni, and T. Schwetz, “Are there sterile neutrinos at the eV scale?,” Phys.Rev.Lett. 107 (2011) 091801, 1103.4570 [hep-ph].
-  C. Giunti and M. Laveder, “3+1 and 3+2 Sterile Neutrino Fits,” Phys.Rev. D84 (2011) 073008, 1107.1452 [hep-ph].
-  C. Giunti and M. Laveder, “Status of 3+1 Neutrino Mixing,” Phys.Rev. D84 (2011) 093006, 1109.4033 [hep-ph].
-  C. Giunti and M. Laveder, “Implications of 3+1 Short-Baseline Neutrino Oscillations,” Phys.Lett. B706 (2011) 200–207, 1111.1069 [hep-ph].
-  J. Conrad, C. Ignarra, G. Karagiorgi, M. Shaevitz, and J. Spitz, “Sterile Neutrino Fits to Short Baseline Neutrino Oscillation Measurements,” Adv.High Energy Phys. 2013 (2013) 163897, 1207.4765 [hep-ex].
-  J. Kopp, P. A. N. Machado, M. Maltoni, and T. Schwetz, “Sterile Neutrino Oscillations: The Global Picture,” JHEP 1305 (2013) 050, 1303.3011 [hep-ph].
-  C. Giunti, M. Laveder, Y. Li, and H. Long, “Pragmatic View of Short-Baseline Neutrino Oscillations,” Phys.Rev. D88 (2013) 073008, 1308.5288 [hep-ph].
-  Planck Collaboration Collaboration, P. Ade et al., “Planck 2013 results. XVI. Cosmological parameters,” Astron.Astrophys. 571 (2014) A16, 1303.5076 [astro-ph.CO].
-  M. Archidiacono, N. Fornengo, C. Giunti, and A. Melchiorri, “Testing 3+1 and 3+2 neutrino mass models with cosmology and short baseline experiments,” Phys.Rev. D86 (2012) 065028, 1207.6515 [astro-ph].
-  S. Gariazzo, C. Giunti, and M. Laveder, “Cosmological Invisible Decay of Light Sterile Neutrinos,” 1404.6160 [astro-ph].
-  G. Mangano, G. Miele, S. Pastor, T. Pinto, O. Pisanti, et al., “Relic neutrino decoupling including flavor oscillations,” Nucl.Phys. B729 (2005) 221–234, hep-ph/0506164.
-  M. Archidiacono, E. Giusarma, S. Hannestad, and O. Mena, “Cosmic dark radiation and neutrinos,” Adv.High Energy Phys. 2013 (2013) 191047, 1307.0637 [astro-ph].
-  J. Lesgourgues and S. Pastor, “Neutrino cosmology and Planck,” New J.Phys. 16 (2014) 065002, 1404.1740 [hep-ph].
-  F. Fritsch and R. Carlson, “Monotone Piecewise Cubic Interpolation,” SIAM Journal on Numerical Analysis 17 no. 2, (1980) 238.
-  J. B. Fred Fritsch, “A Method for Constructing Local Monotone Piecewise Cubic Interpolants,” SIAM Journal on Scientific and Statistical Computing 5 no. 2, (1984) 300.
-  D. Larson, J. Dunkley, G. Hinshaw, E. Komatsu, M. Nolta, et al., “Seven-Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Power Spectra and WMAP-Derived Parameters,” Astrophys.J.Suppl. 192 (2011) 16, 1001.4635 [astro-ph.CO].
-  A. Lewis and S. Bridle, “Cosmological parameters from CMB and other data: A Monte Carlo approach,” Phys.Rev. D66 (2002) 103511, astro-ph/0205436.
-  A. Lewis, A. Challinor, and A. Lasenby, “Efficient computation of CMB anisotropies in closed FRW models,” Astrophys.J. 538 (2000) 473–476, astro-ph/9911177 [astro-ph].
-  BICEP2 Collaboration, P. Ade et al., “Detection of -Mode Polarization at Degree Angular Scales by BICEP2,” Phys.Rev.Lett. 112 (2014) 241101, 1403.3985 [astro-ph.CO].
-  Planck Collaboration, P. Ade et al., “Planck 2013 results. XV. CMB power spectra and likelihood,” Astron.Astrophys. 571 (2014) A15, 1303.5075 [astro-ph].
-  J. Dunkley, E. Calabrese, J. Sievers, G. Addison, N. Battaglia, et al., “The Atacama Cosmology Telescope: likelihood for small-scale CMB data,” JCAP 1307 (2013) 025, 1301.0776 [astro-ph.CO].
-  K. Story, C. Reichardt, Z. Hou, R. Keisler, K. Aird, et al., “A Measurement of the Cosmic Microwave Background Damping Tail from the 2500-square-degree SPT-SZ survey,” Astrophys.J. 779 (2013) 86, 1210.7231 [astro-ph.CO].
-  A. G. Riess, L. Macri, S. Casertano, H. Lampeitl, H. C. Ferguson, et al., “A 3Telescope and Wide Field Camera 3,” Astrophys.J. 730 (2011) 119, 1103.2976 [astro-ph.CO].
-  Planck Collaboration, P. A. R. Ade et al., “Planck 2013 results. XX. Cosmology from Sunyaev-Zeldovich cluster counts,” Astron.Astrophys. 571 (2014) A20, arXiv:1303.5080 [astro-ph].
-  M. Kilbinger, L. Fu, C. Heymans, F. Simpson, J. Benjamin, et al., “CFHTLenS: Combined probe cosmological model comparison using 2D weak gravitational lensing,” Monthly Notices of the Royal Astronomical Society 430 no. 3, (2013) 2200–2220, 1212.3338 [astro-ph.CO].
-  C. Heymans, E. Grocutt, A. Heavens, M. Kilbinger, T. D. Kitching, et al., “CFHTLenS tomographic weak lensing cosmological parameter constraints: Mitigating the impact of intrinsic galaxy alignments,” Mon.Not.Roy.Astron.Soc. 432 (2013) 2433, 1303.1808 [astro-ph.CO].
-  Planck Collaboration Collaboration, R. Adam et al., “Planck intermediate results. XXX. The angular power spectrum of polarized dust emission at intermediate and high Galactic latitudes,” 1409.5738 [astro-ph.CO].
-  W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes 3rd Edition: The Art of Scientific Computing. Cambridge University Press, New York, NY, USA, 3 ed., 2007.