Likelihood analysis of small field polynomial models of inflation yielding a high Tensor-to-Scalar ratio

Inflationary potentials, with Planckian field excursions, described by a 6th degree polynomial are studied. We solve the Mukhanov-Sasaki equations exactly and employ a probabilistic approach as well as multinomial fitting to analyse the results. We identify the most likely models which yield a tensor-to-scalar ratio r = 0.01 in addition to currently allowed Cosmic Microwave Background (CMB) spectrum and observables. Additionally, we find a significant inter-dependence of CMB observables in these models. This might be an important effect for future analyses, since the different moments of the primordial power spectrum are taken to be independent in the usual Markov chain Monte Carlo methods.


Introduction
In a previous article [1], a class of small field inflationary models which are able to reproduce the currently measured Cosmologic Microwave Background (CMB) observables, while also generating an appreciable primordial Gravitational Wave (GW) signal were studied. The existence of such small field models provides a viable alternative to the large field models that generate a high Tensor-to-Scalar ratio. Our exact analysis was shown to give accurate results [1]. Models which yield Tensor to Scalar ratio (r), of less than r ≲ 0.003 were previously studied in [1]. The initial study demonstrated a significant difference between analytical Stewart-Lyth [2,3] estimates and the exact results. This result should be confronted with analyses such as in [4] where the Stewart-Lyth expression is relied upon, and [5] in which the authors use a Green's function approach and perturbation theory, but assume the log of the input is well behaved. Our method extends and improves the method of the model building technique employed in [6,7]. Previous analytical work [8][9][10] has shown that a fourth-order polynomial potential is sufficient to generate a high tensor-to-scalar ratio, even up to r ≳ 0.1. However, it was hard to realize this numerically. It was discovered in [1], that a fifth-order polynomial potential was required for generating 0.001 ≲ r ≲ 0.003. Furthermore a sixth-order polynomial seems to be required for a tensor-to-scalar ratio greater than r ≳ 0.003. A simple explanation is offered by observing (see Fig 1) that increasing r by factor � 10, causes the e-folds per field excursion generated at the CMB window to decrease by a factor of � 3. This means widening the CMB window, and losing the decoupling between the CMB window and the e-fold generating peak. a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 Adding the 6th coefficient pushes the peak from ϕ 2 [0.4, 0.5], to higher values of ϕ and decouples these regions.

Conventions
In this article we follow the conventions of [11]. The Primordial Power Spectrum (PPS) is given by: Our conventions are: The PPS is expanded about the pivot scale k 0 . The pivot scale is usually set a-posteriori, at the scale in which the parameters {n s − 1, α s , β s } are minimally dependent [12][13][14], in Planck + BICEP2 data analyses it is usually set at k 0 = 0.05 hMpc −1 .

Fig 1. A graph depicting À 1=
ffi ffi ffi ffi ffi 2� p � V V � as a function of the inflaton ϕ for two fifth-order polynomial models, and a sixth-order polynomial model. For a fifth-order polynomial model with r = 0.001 (Blue line) the CMB window width is �8 e-folds, while the field changes by about Δϕ � 0.1. Most of the e-folds are generated when ϕ reaches �0.4. When r is increased the CMB window widens, and approaches the e-fold generating peak (Red dots). While marginally affecting the CMB window width, the introduction of an additional coefficient, a 6 , allows to shift the peak to higher values of ϕ, thereby decoupling the CMB window and the e-fold generating peak (Green dash).

Inflationary models
The small field models previously studied in [1] yielded results that are consistent with observable data up to values of r ' 0.003. While these values agree with the current limits on r set by Planck [15,16], we are interested in studying models with higher r. The study of small field models is motivated by their appearance in many fundamental physics frameworks, effective field theory, supergravity [17] and string theory [18] in successive order of complexity. For models with r ≳ 0.003, significant running of running is found. This means that while three free parameters (corresponding to n s , α s , N) were previously needed, we now need an additional free parameter. Therefore we turn to a model of a degree six polynomial potential. Obviously considering higher degree models complicates the analysis by adding other tunable parameters. The potential is given by the following polynomial: It has been shown [6] that the potential can be written as: However, for simplicity, we express the potential as follows: with the subscript 0 denoting the value at the CMB point. By setting ϕ 0 = 0; ϕ end = 1 we limit ourselves to small field models in which Δϕ = 1 in Planck units, with little effect on CMB observables. According to the Lyth bound [19,20], given a Tensor-to-Scalar ratio of r ' 0.01, the lower bound on the field excursion is approximately given by Δϕ 4 ≳ 0.03 m pl . Here Δϕ 4 is the field excursion while the first �4 efolds are generated. Our models satisfy this strict bound, as the first 4 efolds or so typically result in Δϕ 4 � 0.15 which is well above 0.03. The Lyth bound was further extrapolated [21] to cover the entire inflationary period. Applying this approach to models with r � 0.01 yields Δϕ ' 2 m pl . However, in [7], it was shown that in models such as the ones we study, the value of Δϕ can be smaller because � H is non-monotonic. In this case, Δϕ = 1 m pl from the CMB point to the end of inflation is consistent with the Lyth bound.
When the coefficients {r 0 , a 2 , a 3 , a 4 } are fixed, the remaining coefficients are related by: The procedure of finding f 1 and f 2 was explained in detail for the degree 5 polynomial models in [1], and here we follow a similar procedure for the degree 6 models. So, ultimately, the model is parametrized by 5 parameters: the two physical parameters r 0 and N and the three other parameters (a 2 , a 3 , a 4 ) that are used to parametrize the n s , α s , β s parameter space. It should be pointed out that N is not an observable, rather N � 50 − 60 is a 'soft' constraint. Strictly speaking, N depends on the reheating temperature and only its maximum value can be determined. However for simplicity we treat N as an observable, in order to facilitate the study of a large sample of models.

Coefficient extraction methods
In this section we explain the two methods for calculating the most likely coefficients {a 2 , a 3 , a 4 }, given a large number of simulated models and the likelihood data for the CMB observables. This data is available through MCMC analysis of CMB data, such as the Planck data.

Likelihood assignment method-Gaussian extraction
To each potential, after calculating the observables n s , α s , β s , we assign a likelihood. For each observable we calculate the likelihood according to the CosmoMC [22] likelihood analysis of the data sets used. We then assign the product of the likelihoods L ðn s Þ � L ða s Þ � L ðb s Þ to the potential. A concrete example is the following: suppose we extract the trio (n s , α s , β s ) = (0.96, 0.011, 0.024), we look up the likelihoods: L ðn s ¼0:96Þ � L ða s ¼0:011Þ � L ðb s ¼0:024Þ . We now multiply them, and so the likelihood attached to that specific model which yielded these observables is given by L potential ¼ L ðn s ¼0:96Þ � L ða s ¼0:011Þ � L ðb s ¼0:024Þ . We proceed to extract likelihoods for the different coefficients by a process of marginalization. The expectation is that this method will yield a (roughly) Gaussian distribution for each of the values of a 2 , a 3 , a 4 . The advantage of this method is in yielding not only the most likely value, but also the width of the Gaussian. This width can then be used as an indication for the level of tuning that is needed in these models.

Possible pitfalls.
This method of likelihood assignment is vulnerable in two ways: a). To be valid, this method requires a uniform cover of the relevant parameter space, by the potential parameters. If the cover significantly deviates from uniform, the results might be skewed by overweighting areas of negligible weight, or underweighting areas of significant weight. Fig 2 shows a mostly uniform cover. b). Since L n s ;a s ;b s ' L ðn s Þ � L ða s Þ � L ðb s Þ only if the paired covariance is small, we must make sure that this is the case. In our underlying CosmoMC [22] analysis this is indeed the case. The covariance terms are, in general, one to two orders of magnitude smaller than the likelihoods at the tails of the Gaussian. c). We also run the risk of false results if the fit we apply to the data points produced by the numerical analysis yields a large fitting error. However the fitting error of the polynomial function to the log PPS − log k data, is usually of the order of 10 −6 . This fitting is done over 30 data points generated by the MS equation numerical evaluation, for each potential. The error is calculated as D ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi P 30 1 ððlog PPSÞ i À fitððlog kÞ i ÞÞ 2 q , thus roughly speaking the error per data point is of the order of 10 −6�7 . We conclude that the log PPS function is well fitted.

Multinomial fit
Another method for calculating the most likely coefficients is by fitting the simulated data with a multinomial function of the CMB observables. We aim to find a set of functions F i such that, for example, a 2 = F 2 (n s , α s , β s ). We assume that this function is smooth and thus can be expanded in the vicinity of the most likely CMB observables. Hence, we can find a set of multinomials (F 2 , F 3 , F 4 ), such that: We have found that a quadratic multinomial is sufficiently accurate and that using a higher degree multinomial does not improve the accuracy significantly. Thus we may represent these by a symmetric bilinear form plus a linear term, as follows: where O = (n s , α s , β s ), B i is the bilinear matrix, and the linear coefficient vector is A i .

Pivot scale
So far, we discussed matching potentials and their resulting PPS around the CMB point. However, in order to correctly compare results of the PPS to observables, one has to take into account the pivot scale at which the CMB observables are defined. Since, in this case, the pivot scale is given by k 0 = 0.05 hMpc −1 , and the CMB point is at k � 10 −4 hMpc −1 , the observables in the CMB point and k 0 should be related in a simple way only if the spectrum varies slowly with k. This is not true for the case at hand. Two potentials can yield very different power spectra near the CMB point, and nevertheless yield the same observables at the pivot scale. These degeneracies, stem from our limited knowledge of the power spectra on small scales, and at the CMB point. For concreteness take two PPS functions, one that is well approximated by a cubic fit near the pivot scale, and the other that is well approximated only when we consider a quartic fit. Suppose, additionally, that these two PPS functions have the exact same first three coefficients, it follows that they yield the exact same observables {n s , α s , β s }. However if we go to sufficiently small scales, or sufficiently high k values, these functions will diverge. This is also true at the large scale end, where the CMB point is set. Hence the degeneracy. A possible solution to this problem, is classifying the resulting power spectra by the level of minimal good fit. We define a good fit as one in which the cumulative relative error D ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi P k ðlogðPPSðlogðkÞÞÞ À fitðlogðkÞÞÞ 2 q , is less than 10 −7 . Given a single power spectrum, we fit our result with a polynomial fit, increasing in order until the accumulated relative error is sufficiently small. The minimal degree polynomial fit that approximates the log(PPS) function to the aforementioned accuracy is called the minimal good fit. We then study separately power spectra that are well fitted by cubic polynomials, quartic polynomials etc. In this way we make sure that we compare non-degenerate cases.

Monte Carlo analysis of Cosmic Microwave Background with running of running
In [11] it was shown that the inclusion of additional parameters, i.e., the running of the spectral index (α s ), and the running of the running (β s ) resolves much of the tension between different data sets. In this section we briefly discuss the effect of considering non-vanishing α s and β s on the most likely shape of the PPS. First we find n s when it is the only free parameter. We then use n s , and α s as the free parameters, and finally we conduct an analysis with n s , α s , and β s as the free parameters. The shape of the power spectrum changes significantly when running of running is considered.
The data sets that were used are the latest BICEP2+Planck baseline [15], along with the low l's [23], low TEB and lensing likelihoods. The results of these analyses are given in Table 1, as well as in Fig 3. As expected the resulting power spectra converge at the pivot scale k 0 = 0.05 hMpc −1 . However for lower k's, the resulting spectra diverge considerably, consistent with  Likelihood analysis of small field polynomial models of inflation yielding a high Tensor-to-Scalar ratio cosmic variance. Notably, the spectra also diverge at higher k's. This indicates the inability of current observational data to constrain the models in this range of k's. This inability is also demonstrated in Fig 4 where, for l > 1500, the most restrictive data cannot rule out models with significant running, or running of running. Fig 4 also shows that the three models are virtually indistinguishable in terms of the observed C l 's. The conclusion is that we will need additional accurate data from smaller cosmic scales to be able to differentiate between the three scenarios. These extra e-folds might come from future missions such as Euclid [24], or from μ-type distortion data [25,26].

Results
We apply the methods discussed in Section 3 to the degree six polynomial inflationary potentials that yield r = 0.01. We calculate the most likely coefficients and extract the resulting most likely polynomial inflationary potential, the results are given in Tables 2 and 3. The PPS resulting from this inflationary potential is then calculated in order to confirm that the most likely coefficients reconstruct the most likely observables.

Results for degree six polynomials that yield r = 0.01
In Fig 2 we showed a cover for the joint likelihood map of n s − α s , of about 2000 potentials with r = 0.01. The cover is approximately uniform, thus we were able to assign likelihoods to every potential we study, as previously discussed.  Likelihood analysis of small field polynomial models of inflation yielding a high Tensor-to-Scalar ratio By a process of marginalization, as discussed in Section 3, we extract the most likely coefficients, which yield the likeliest observables. This process is represented graphically in Fig 3. The results are shown in Table 3 and the PPS reconstruction is shown in Fig 5. The advantage of this method is that it also determines the deviation from the average value. This can be used as an indicator for the level of tuning that is required to construct the most likely small field model. A discussion of tuning in field theoretic models can be found in [27], as well as in [28] and [29]. In most cases the tuning level can be viewed as simply Dx i x i , which in this case are given by (0.375, 0.27, 5.5) for a 2 , a 3 , a 4 respectively. The width of the Gaussian fits for {a 2 , a 3 , a 4 } are {0.015, 0.041, 0.112} respectively This is shown in Fig 6. These widths represent the effective area in parameter space that yields observables within the 68% CL. Which is another measure of the tuning required in the sixth-order polynomial models.
Recalculating the CMB observables that this most likely model yields, we find n s = 0.9687, α s = 0.0089, β s = 0.0176. These values are very close to the most likely values found from the previously discussed MCMC analysis of the BICEP2+Planck data. The resulting scalar index fits the most likely value in Table 1 exactly, while α s and β s deviate from these values by no more than 12.5%. We found that this is a relic of the binning method. Adding more models to  Likelihood analysis of small field polynomial models of inflation yielding a high Tensor-to-Scalar ratio the simulated data and refining the binning process results in even better proximity to the desired values.
Using the method of multinomial evaluation (3.2), we found the multinomial coefficients for each of the model degrees of freedom. For instance for a 2 Since n s � 1, and α s and β s are of the order of 10 −2 , the above result suggests that a 2 is primarily dominated by n s . Similarly, we have found that a 3 is dominated by a linear combination of α s , and β s , and a 4 is primarily dominated by β s . This method yields the most likely CMB observables with comparable accuracy to the previous method upon recalculation: n s = 0.9684, α s = 0.0077, β s = 0.020. Table 2 contains a comparison between the results of the two methods.

Most likely potentials
Since n s is better constrained, we opt for the analysis that yields a more precise value of n s . The leading 6th degree polynomial which yields r = 0.01 at the proper pivot scale, is thus given by: Upon initial investigation, it seems these models produce a relatively flat tensor power spectrum. This might motivate future research of the tensor power spectrum predictions and constraints.
this is evident in Fig 7. It should be stressed that this is a phenomenon associated with the models and not with the observational data. This is supported by the small n s , α s , β s pairedcovariance found in [11], implying weak dependence among observables in the data itself. One might think that the previous findings in [16, Figure 23] indicate that n s and α s are dependent. However, the graph shows a dependence between α s and n s,0.002 which is the scalar index evaluated at k = 0.002 hMpc −1 . Taking some initial n s evaluated at k 0 , it follows that n s evaluated at some other scale, depends on the index running α s . Indeed, when one examines the color coding in [16, Figure 23], which represents n s at the pivot scale, it is clear that n s and α s are independent.

Summary and outlook
A large sample of potentials that yield r = 0.01 and conform to the allowed observable values was successfully generated. The sample provided a uniform cover of the allowed region of parameters which enabled us to assign likelihoods to each of the potentials and extract each coefficient's likelihood (3.1). Another approach was implemented, representing each coefficient as a multinomial function of the observables (3.2), which yielded similar results. A most likely small field potential giving rise to r = 0.01 and (n s ' 0.9694, α s ' 0.009, β s ' 0.0175), was identified, and its power spectrum simulated. An interesting inter-dependence of (n s , α s , β s ) was found in these models, which may have some bearing on future MCMC analysis. We hope to perform such an MCMC analysis, either with priors that include this dependence or with a numerical scheme that reflects it.
The Planck collaboration may soon release additional analysis products, and BICEP3 is also expected to release results in the near future. Thus, it might be possible to check our prediction for the tensor-to-scalar ratio. However, ruling out models of the class discussed in this paper might be a more difficult task due to the lack of constraining power of current observations in the range l > 1500. We expect that either S4 cosmology or μ-distortion data will be able to resolve this in the foreseeable future by adding observational data on smaller scales.

Author Contributions
Formal analysis: Ira Wolfson. Likelihood analysis of small field polynomial models of inflation yielding a high Tensor-to-Scalar ratio