Small field models with gravitational wave signature supported by CMB data

We study scale dependence of the cosmic microwave background (CMB) power spectrum in a class of small, single-field models of inflation which lead to a high value of the tensor to scalar ratio. The inflaton potentials that we consider are degree 5 polynomials, for which we precisely calculate the power spectrum, and extract the cosmological parameters: the scalar index ns, the running of the scalar index nrun and the tensor to scalar ratio r. We find that for non-vanishing nrun and for r as small as r = 0.001, the precisely calculated values of ns and nrun deviate significantly from what the standard analytic treatment predicts. We study in detail, and discuss the probable reasons for such deviations. As such, all previously considered models (of this kind) are based upon inaccurate assumptions. We scan the possible values of potential parameters for which the cosmological parameters are within the allowed range by observations. The 5 parameter class is able to reproduce all of the allowed values of ns and nrun for values of r that are as high as 0.001. Subsequently this study at once refutes previous such models built using the analytical Stewart-Lyth term, and revives the small field brand, by building models that do yield an appreciable r while conforming to known CMB observables.


Introduction
Recent years have shown an increase in cosmological observational data, largely due to the Planck mission [1], and the searches for primordial gravitational waves (GW) signal in the cosmic microwave background (CMB) by terrestrial experiments such as BICEP2 and the Keck Array [2,3]. Inflation [4][5][6][7] is widely accepted as a probable model for the origin of our universe, one of the hallmarks of which is the production of GW (for example [4,8]).
Sensitivity for detecting GW in the CMB have, over the years, improved constantly. Constraints on the tensor-to-scalar ratio r were tightened [1][2][3][9][10][11][12] and it is expected that a sensitivity level of r ≲ 0.03 be reached in the near future [13]. Furthermore one can optimistically expect the next decade to yield measurements of r ≲ 0.001 or better [14]. Constant headway is also made in the model building front, as some models become less probable, while others gain dominance.
We study a class of models that were proposed by Ben-Dayan & Brustein [15]. These models sport, along with the ability to conform to known observable quantities such as the a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 primordial power spectrum (PPS) scalar index (n s ), and its running (n run ), the generation of appreciable amplitude of GW signal. This type of models appear in many fundamental physics frameworks, such as effective field theory, supergravity and string theory. A discussion regarding small field models and the possibility of GW generation [16][17][18] soon followed. In these models, high values of r in the CMB are generally associated with a scale dependence of the scalar power spectrum. We study the models proposed by Ben-Dayan & Brustein using exact calculations. For each model, we solve the background eqautions and the Mukhanov-Sassaki (MS) equations [19][20][21] to obtain a primordial power spectrum. This process is applied to a large sample of models and allows us to study the dependence of cosmological parameters on the potential parameters with unprecedented accuracy.
Significant differences between analytical predictions of the commonly used Stewart-Lyth (SL) expressions [22,23] for CMB observables and the precise results were found. These discrepencies were already found in [24], however all previous discussions of such models [15][16][17][18] nevertheless heavily rely on the SL expression, thus the importance of this discrepancy is enhanced. These differences arise from several factors, chief among them is breaking of slowroll hierarchy. When the hierarchy is broken the time derivatives of the first and second slowroll parameters ( H , δ H ) cannot be neglected. Hence, rather than general arguments, these models require precise calculations in order to study their validity. This also means that, in some cases, it is not possible to use Hankel functions as an approximate solution of the MS equation. This was discussed in some length in [25]. In other cases the Hankel functions can still be used, but either require adjustments, or some additional requirements must be met as in [26].

The primordial power spectrum and the cosmological parameters
The primordial power spectrum (PPS) is traditionally characterized by its spectral index n s and the index running n run (sometimes also denoted as α), which are given by the first and second logarithmic derivatives of the logarithm of the PPS: where aH = k denotes the CMB scale.

A brief review
The process of relating slow-roll parameters to the power spectrum is documented extensively in [22] and described in broad strokes in [23]. Following is a brief review of the process. In principle, the process of deriving the PPS given an inflationary potential is straightforward. The background evolution equations are solved to construct the pump field: where a dot denotes a derivative with respect to cosmic time. The MS equations [19][20][21] are: in conformal time τ and in Fourier space with wave vector k, where ω(τ) is given by: Here a prime denotes a derivative with respect to conformal time. The eigenfunctions U k (τ) of these equations are recovered. Evaluating these at a time τ later than the latest freeze-out time yields the PPS generated by the inflationary potential V.
In [22], Stewart & Lyth derive an analytic expression for the spectral index of a wide array of inflationary scenarios. They first assume a slow-roll inflation, sufficiently slow, so that both slow roll parameters, can be approximated by constants. It is useful to rewrite the quantity Z 00 Z as: For strictly constant H , δ H , withC a constant. In this case, the background solution corresponds to power law inflation. The resulting MS equations becomes the Bessel equations which can be solved analytically. When the Bunch-Davies boundary conditions are imposed, the resulting solution is given by a Hankel function of the first kind: with the index ν given by: The resulting power spectrum is given by: Upon derivation of P R with respect to log(k) and evaluating the derivative at k = aH the scalar index is obtained, Inserting ν from Eq (11) one gets: with b = 2 − log(2) − γ, γ being the Euler number. The resulting scalar index running is given (for instance in [23]) by: We note that which in the slow-roll paradigm is usually taken to be small, appears in both Eqs (14) and (15). It might be tempting then, to drop these terms. However, this term was shown in [22], and later in [23] to be required for a better than *1% accurate prediction of the CMB observables. The authors of [22] then proceed to connect slow-roll parameters to the potential and its derivatives by a process of Taylor expanding with respect to cosmic time, and re-substituting the Friedman equations to 2 nd order. Thus they are able to obtain an analytical expression that connects the PPS observables directly to the potential and its derivatives to a high degree of accuracy. Following the same procedure for the running of the scalar index, yields (again, to 2 nd order): where the subscript 0 denotes evaluating the quantity at the CMB point and the subscript V denotes that these are potential derivatives: However, when H or δ H are not strictly constants, the analytic solution to the MS equation is not generally known. We show that the above analytic expressions are not accurate enough for certain models where H and δ H are time-dependent. Therefore, one has to use the precise calculations which takes into account the deviations of the MS equation solutions from the Hankel functions. Another observable, used to parametrize the amplitude of GW at the onset of inflation is the scalar-to-tensor ratio r, r = 16 H . In fact, should we ever detect a GW signal, we would be able to directly probe the energy scale of inflation [27].

Inflationary models
Small field models of inflation in which inflation occurs near a flat feature, a maximum, or a saddle point are studied (see [28] for a review). This class of models is interesting because they appear in many fundamental physics frameworks, effective field theory, supergravity [29] and string theory [30] in successive order of complexity. Our focus on such models is also motivated by the expected properties of the moduli potentials in string theory. More generally speaking these type of models can be viewed as a Taylor expansion approach to other models [31]. A different more observable-oriented classification of models can be found in [32], in which analysis our models fall into the toward-exit class.
In general, inflation will occur in a multi-dimensional space, however, the results for multifield inflation cannot usually be obtained in a simple way. In many known cases it is possible to identify a-posteriori a single degree of freedom along which inflation takes place. To gain some insight about the expected typical results effective single field potentials can be used.
Generic small field models predict a red spectrum of scalar perturbations, negligible spectral index running and non-gaussianity. They also predict a characteristic suppression of tensor perturbations [33]. Hence, they were not viewed as candidate models for high-r inflation. Large field models of inflation are thus the standard candidates for high-r inflation.
In [15], a new class of more complicated single small field models of inflation was considered (see also [16]) that can predict, contrary to popular wisdom [27,34], an observable GW signal in the CMB (see also [35]). The notion that observable signal GW precludes small field models partly stems from [34] and similar analyses that study monomial potential models as small field models. The spectral index, its running, the tensor to scalar ratio and the number of e-folds were claimed to cover all the parameter space currently allowed by cosmological observations. The main feature of these models is that the high value of r is accompanied by a relatively strong scale dependence of the resulting power spectrum. Another unique feature of models in this class is their ability to predict, again contrary to popular wisdom [36], a negative spectral index running. The single observable consequence that seems common to all single field models is the negligible amount of non-gaussianity. In [24] the inflationary potential was Taylor-expanded up to order 4. The approach applied in [24] is similar, however it seems only potentials that are monotonic in the entire CMB window were considered.
The current work yields corrected predictions of this class of models by a systematic highprecision analysis, thus providing a viable alternative to the large field-high r option. The analysis of [15] is extended, in preparation for a future detailed comparison of the models to data. This is done in order to simplify the parametrization of the potential and facilitate a comprehensive numerical study.

Inflaton potentials
The following class of polynomial inflationary potentials proposed in [15] is: The virtue of these models from a phenomenological point-of-view is the ability to separate the CMB region from the region of large e-fold production. Hence, these potentials can produce a very different spectrum early on, than in the later stages of inflation. Fig 1 illustrates this point, with separate CMB region and e-fold generation region. In the context of both classification systems mentioned, current observational data weakly support these [37,38]. However the small field model studied in [38] are monomial potential models of the form V / 1 − a p ϕ p , which are different from many of our models.
In many models ε V * 1/N 2 , η V * 1/N 2 , and the time derivative d Hdt can approximately be replaced with a factor of 1 N 2 [39]. In the above models this standard hierarchal dependence is broken, they have a more complicated dependence while obeying the slow-roll conditions H , δ H ( 1. In [15] it was shown that these models can be written as: Here r 0 , η 0 , α 0 are defined as respectively. The subscript 0 means that these are the values at the CMB point. Specifically for a potential of the form V / 1 þ P 5 p¼1 a p 0 p , the SL analytic expression for the scalar index and its running (Eqs (16) and (17)) is given by

Reduced parameter space
The potential in (22) is a small field candidate, which after some scaling and normalization, depends on four free parameters. One parameter is used for setting r 0 at the CMB point, and thus the predicted amplitude of the GW signal produced, while the other two parameters are used to parametrize the n s − n run -plane. The fourth parameter determines the number of efolds from the CMB point to the end of inflation. ϕ end is set to 1 to simplify the analysis. It follows that Suppose we want inflation to end at ϕ = α, we can rescale ϕ: In this formulation, whereã p ¼ a p a p . Since this is the exact same potential, it follows the exact same CMB observables are yielded. Thus, applying condition (25) can be viewed as a scaling scheme for the different terms in the potential which does not limit the generality of our results. Substituting the expression for the potential and its derivative at ϕ = 1 we get: a 4 is now given in terms of the other coefficients: Using the standard definition for the number of e-folds N ¼ , and the approximation Vð0Þ ¼ 1 þ P 5 p¼1 a p 0 p ' 1 yields a rough estimate for a 5 as a function of N, This estimate is then used as a starting point to refine a 5 by solving the background equations iteratively thereby obtaining the accurate coefficient a 5 that yields the correct N. Thus a 4-dimensional parameter space r 0 , a 2 , a 3 , N is defined. The parameters a 2 , a 3 are constrained by the requirement |a 2 |, |a 3 | ( 1, a 1 is constrained by the observable value of r and a 5 is determined by the other parameters and by the number of e-folds (taken to be in between 50 * 60). The PPS considered is in the range of the first log(2500)*8 e-folds of inflation.

Precise evaluation of the cosmological parameters
Using the analytic results in Eqs (16) and (17) it can be concluded that the above class of models can cover the part of the n s − n run plane of interest [15]. However, several approximations are made along the way. Significant deviations from analytic prediction are found, of the order of a percent or so in estimating n s and 50% or more in estimating n run . The unavoidable conclusion is that rather than a general argument, a precise calculation is necessary to extract the cosmological parameters these models yield.

From potentials to cosmological parameters
The process of calculating the cosmological parameters for a given potential is the following. A potential candidate is built by setting a parameter (for instance r 0 ), and randomly drawing the other parameters (in this example a 2 and a 3 ) from a uniform distribution function. The limits of this distribution function are set by hand and require a process of trial and error (guided by theoretical insights such as overall behaviour of precisely calculated n s and n run ). After the 3 first parameters are fixed, Eq (29) is used to relate a 4 to a 5 , and the value of a 5 is calculated for a the desired value of N. a 5 is found as explained above, with no approximations. The choice of which parameters to fix and which to randomly draw relies on the observables studied. For each potential the Friedmann equations and the inflaton scalar field equation are solved. The initial conditions are set such that integration starts 3.5 efolds before the CMB point with _ 0 ¼ 0. In that fashion we ensure that we are well within the slow roll regime, and on the attractor solution when the CMB point is reached. The solution is used to construct Z and ω k as described in (4) and (6). The eigenfunctions for the MS equations (5) are found and used to calculate the power spectrum. Finally we provide a fit for the power spectrum, from which we extract n s and n run .

Cosmological parameters of small field models
In this section we present the results of evaluating cosmological parameters for many small field models. In Fig 2 we show an example for which we calculate n s and n run for about 1100 models with a fixed scalar to tensor ratio r 0 = 0.001. The results are shown on a n s − n run joint probability graph with the 68%, 95% contours that are the probability estimators as yielded by a CosmoMC [40] ΛCDM +index running model run, with the most recent Bicep & Planck data (including WMAP 9-year mission) [12].
The reason for choosing the value of r 0 = 0.001 (and not a higher value, for example, r = 0.01) was the following. We discovered that as we increased the values of r, the inflaton potentials needed to be more complicated and additional parameters were required. Also, we encountered several technical difficulties which we were able to resolve for the lower values of r. Solving these difficulties and constructing a reliable framework for numerical calculations of the CMB observables is an essential step towards building models with higher values of r, which we intend to do in a future publication.
We allow the values of n s to vary quite substantially, rather than restrict them to the narrow range that is allowed by the data. Our idea is that when r and n run are free to vary, the constraints on n s are relaxed in a significant way. The reason is that there is some degeneracy among the parameters. This is validated in the preliminary analysis that we present in this paper. In addition, despite of the fact that some models have yielded an almost flat (and some even a blue) n s and therefore are in conflict with the data, we find their analysis useful because insight regarding the departure of precisely calculated results from what the analytic SL term ((16)) predicts (see below), is gained.

Evaluating cosmological parameters for fixed r 0
The n s − n run plane was covered with models which yield a fixed value of r 0 = 0.001. The cosmological parameters of some 3500 potentials were calculated. Fig 2 shows cosmological parameters of * 1100 models. A significant number of the models yield values of n s and n run within the 68% and 95% likelihood region. The most probable value for V 00 V ¼ À 0:0052 AE 0:0034. This The contour curves are the 68% and 95% confidence estimators, obtained from a CosmoMC ΛCDM + index running model run [40] using the Planck & Bicep joint data analysis [12].
The pivot scale used in the analysis is k pivot = 0.05 h Mpc −1 , which is the same scale as in [12]. https://doi.org/10.1371/journal.pone.0197735.g002 is within the 68% CL Planck results, with or without including high-l polarization data. The third coefficient values are given by V 000 V 0 V 2 ¼ 0:0138 AE 0:0065, which is in better agreement with the result without high-l data. However the 2015 Planck analysis [41] sets 4 0 which might bias the results slightly. In the 2013 analysis [42] this was not done, and our results agree with their analyses, including our values for V ð4Þ V 0 V 2 . Additional factors that contribute to the difference in analyses, are the approximate connection between Hubble flow functions i and the potential derivative quantities V ; Z V ; x 2 V . An interesting feature of these models is the departure of precisely calculated results from what the analytic SL expression (16) predicts, to be discussed later. It might be possible to cover the n s − n run allowed region with models with a higher scalar-totensor ratio. However the treatment of models which yield higher r is more complex, since by increasing r, one is forced to consider a larger Δϕ range CMB region. The CMB region (see Fig  1) is roughly 3 times larger in ϕ for models with r 0 = 0.01, thus it will typically result in a running of running of the power spectrum.

Evaluating cosmological parameters for fixed η 0
The effects of varying r 0 on the resulting power spectrum were studied. In order to do this η 0 was set to 0 for simplicity, and the n s − n run plane was covered with models of varying r 0 and α 0 . Fig 3 shows the results of this study.
Notice that the effect of varying both r 0 and α 0 on the changes in the value of n s is more pronounced than expected. Usually one expects n s − 1 to first order to be / À 3r 0 8 and thus Δn s /Δr 0 ' 10 −4 * 10 −5 . At second order, we expect n s −1 to be / a 0 15 and thus Δn s /Δα 0 ' 10 −3 , whereas in this case the change in n s is of the order of 10 −2 . A possible explanation to this phenomenon is a discrepancy between the analytic predictions made using (16) and the precise calculations (see below).

Comparison of calculated results and the Stewart-Lyth analytic predictions
An additional study of models with a larger value of r was conducted. This was done in order to confirm the ability of this class of models to produce significant GW signal, while yielding acceptable values of n s and n run . For r ≳ 10 −3 , a significant deviation from the analytical expressions in Eqs (16) and (17) was found. Potentials that by the standard analytic treatment should have yielded acceptable observables, were wide off the mark. On the other hand potentials which were supposed to be ruled out, yielded observables inside the n s − n run acceptable domain. Fig 4 elucidates this point, with a potential (Fig 4, upper panel), for which r 0 = 0.001. The resulting n s and n run are within the 68% probability allowed region, while the analytic expressions yield values outside the 95% probability allowed region. The example in Fig 4, lower panel shows the opposite also occurs. Table 1 contains as examples ten specific potentials that were chosen such that the precise results for n s and n run are within the accepted values. All of the models produce a tensor to scalar ratio r 0 = 0.001. The table also contains the analytic predictions made with Eqs (16 and 17). As can be seen from the table, the discrepancy between the analytic predictions and the precise calculations can be quite significant for n run . The spectrum is composed using 15 k − modes, and the error in the rightmost column is the cumulative error defined as   Small field models with gravitational wave signature supported by CMB data

Possible explanations of the source of deviation between precise results and analytical estimates
From the discussion in Appendix B, one can easily see that the definition of ν, is potentially the most significant discrepancy. The effect of this change in definition is an error of less than about 0.4%. Table 2 contains three examples of potentials. Two yield observables that are within acceptable limits, and a third shows an excluded precise result with an allowed analytic prediction. These examples are used to study the origin of discrepancy. The differences between the slowroll parameters defined via the potential vs. their definition in terms of time derivatives are also discussed in appendix B. We have found, that in the degree 5 polynomial potentials that were studied, small but significant departures from the relations in Eq (B:44) are detected. For instance δ H = −0.0016 and δ V = 0.001 at the time when n s is evaluated. Table 3 contains values   Table 1. Shown is a table of 10 potentials constructed such that r 0 = 0.001, and N = 60. The parameters a 2 and a 3 are constructed by randomly drawing from a uniform distribution as explained in Section 4. The discrepancy in n s is around 0.8%*1.25%, while the n run discrepancy is much more pronounced.   Table 4 contains the scalar index for the corresponding potentials (examples 1,2 and 3).
The overall effect of this discrepancy can sometimes amount to a 5 * 8% error towards higher values.
Finally there is also a significant difference in the derivatives of ν and ν SL , ν SL being ν in the SL formulation: where time dependency of the slow roll parameters is neglected. This difference is mainly due to neglecting the term d H 0 ⃛ H € 0 in the definition of Z 00 Z . This yields a difference in the derivative terms of the order of 0.02 * 0.04, which in turn is responsible for a difference in n s of the order of 4 * 8%. Using ν SL instead of the full term, tends to drive the resulting n s downwards.
The tendencies of the two aforementioned errors are opposite, and so they might sometimes cancel each other. This makes it possible to get an accurate result using the standard SL expression for a specific potential, but studying a collection of such potentials reveals the incomplete nature of this cancellation. Table 4 shows the different results using different methods of deriving the scalar index. We use three different analytical methods: (1) Eq (B:45)-The SL original method, extracting a term for the scalar index as a function of the potential and its derivatives, (2) Eq (14)-The SL original method, but not relating slow roll quantities to potential and derivatives, and (3) Using the same methods as the SL analysis, with the definition for ν as in Eqs (B: 36,40,42). From this analysis it seems the origin of the most significant error is the inaccurate relations Small field models with gravitational wave signature supported by CMB data between slow roll parameters and their potential and derivatives counterparts. Second in significance is the definition of ν with the full Z 00 Z t 2 expression, along with the proper derivation of @n @logðkÞ . The evaluation of −τaH(1 − H ) = 1 is off by * 0.04% and the difference between c 3 2 À Á and ψ(ν) yields a correction of the order of * 0.01%.
There might be additional factors that stem from the temporal dependence of ν in the MS equation, however, these mostly affect the running of the spectral index, and are harder to estimate accurately.
Taking these approximations into account, lowers the discrepancy to the order of 0.5%, in a consistent manner. Another possible explanation is that the time-dependence in (8), modifies the corresponding o 2 This could lead to modified solutions for the MS equation. An example of this phenomenon is given in [43], where the Hankel functions were replaced by the Whittaker functions (albeit these models are observationally excluded). It is worth mentioning that this avenue was studied analytically by Dodelson & Stewart [44,45]. They derived an expression for the scalar index in cases where the slow-roll hierarchy breaks down. However, this analysis was not checked numerically. Additional derivation attempts aiming at yielding better precision analytical expression for the scalar index n s were made in [26,46]. Specifically [26] supplies an analysis of the predicted level of accuracy as a function of the horizon flow functions 1 H and 2 2( H + δ h ), in Fig 6. The different approximation schemes were put to the numerical test in the context of our models. Fig 7 shows that all methods of approximation yield results varying in accuracy and precision levels, it also shows however that the SEG approximation is the best candidate to improve on, since on average they yield errors of less than 1%. Studying results where relative errors in n s are over 1%, for each expression and locating it on the 1 − | 2 | diagram in Fig 8 reveals that the analysis offered in [26] is not completely applicable to our models. Fig 9 shows that for the models studied, even though the conditions outlined in [26] are met, and 1 2 × ΔN < 10 −2*3 for ΔN = 60, the relative error between numerical result and SEG-CH expression can be above 1%.

Summary, conclusions and outlook
An interesting class of models that can produce a high tensor-to-scalar ratio while conforming to observable values of n s and n run was presented and studied. This work has shown that while the arguments for small field model validity presented in [ [15,16] generally apply, the method by which they choose favoured models is based on approximations that are not always accurate enough for the cases studied. While this work argued this possible weakness, it also supplied a remedy: The precise calculation method. Using precise calculations points to new candidates previously disregarded. Specifically, The predictions made using the standard SL analytic expressions were found to deviate by more than 1% from the actual results, for many models in this class. Other approximate expressions such as those suggested in [26,44,46] are, in general, better than the SL expressions, but still miss by more than 1% in some cases. We hope to extend this work to models that produce higher values of r, and determine the best candidate for small field inflationary models [47].

Power law inflation
The accuracy of the procedure is tested by using the benchmark case of power law inflation (a / t p ). This is the only case for which the analytic results are exact since H and δ H are Regions in the 1 -| 2 | parameter space where the spectral amplitudes could be calculated with an accuracy better than 1%, according to analysis presented in [26]. In the dark shaded region the Stewart-Lyth (SL) approximation [22], as well as all other approximations are supposedly sufficiently accurate. Second-order corrections, as calculated by Stewart and Gong (SG) [46], extend that region to the light shaded region. The constant horizon approximation at order n (chn), and the growing horizon approximation at order n (ghn), do well below the thick line. The rays indicate where the corresponding higher order corrections are necessary. The thick line itself is the condition 1 | 2 | < (A/100%)/ΔN, with ΔN = 10 and A = 1%. We study these approximations and others, and find that our models defy these analyses. Figure and caption adapted with author permission from [26]. https://doi.org/10.1371/journal.pone.0197735.g006 Small field models with gravitational wave signature supported by CMB data Fig 7. Around 50,000 of our models numerically simulated and compared to different analytical expressions reveals a varying level of accuracy in predicting the correct scalar index. The figure shows only a partial sample of *8000 restricted to 1 < 0.0275, | 2 | < 0.0275 and 0.96 < n s < 0.99. Each data point is a relative error between the numerical result of a model and an analytical expression from [44] (DS,green circles), [46] (SG,red diamonds), [26] (SEG-GH, growing horizon variant-blue triangle, and SEG-CH, constant horizon variant-inverted cyan triangle), and the usual SL [22] expression (purple squares).
https://doi.org/10.1371/journal.pone.0197735.g007 Small field models with gravitational wave signature supported by CMB data The cosmological parameters are given by: n run ¼ 0: In S1 Fig, results of the convergence of the precise calculations to analytic predictions for the case of power law inflation are shown. The overall shape of both precise calculations and analytic curves agree and a relative error in n s , estimated by: n precise s À n analytic s À Á = 1 2 n precise s þ n analytic s À Á ; is of the order of 10 −2 * 10 −3 %. The method for error estimate in n run is more subtle, since the correct value of n run is zero. In order to assess our error in n run the following diagnostic was therefore used: the difference between the precise and analytic n s is divided by the difference in log(k). The criterion for convergence is that the absolute value of n run is smaller than Dn s DlogðkÞ . S2 Fig displays the convergence of the precisely calculated results for n run , using the diagnostic Dn s DlogðkÞ . It is apparent that n run is always bounded from above by Dn s DlogðkÞ . Additionally the extracted n run is an order of magnitude or so below current observational bound.

Quadratic potentials
As we aim to study models that produce slow roll parameters which are time dependent, we need to check the precision of the numerical code against such models.
Consequently we tested the accuracy of our calculations for quadratic potentials of the type While satisfying the condition 1 | 2 | × ΔN < 10 −2*3 , for ΔN = 60, one finds a relative difference of well over 1% between analytical predictions and numerical results. This is in contrast to the analysis proposed in [26]. https://doi.org/10.1371/journal.pone.0197735.g009 In these cases the analytic expression for the scalar index is given by, Here N is the number of efolds and b is the same as in (14). S3 Fig presents the results of this study, as relative errors between precise calculations and the SL analytic expressions. These results are accurate to * 0.1%. However there is a systematic error that is traced back to the inaccuracy of the approximation: A shift N ! N − 0.8 is sufficient to reduce the systematic error such that the relative error is of the order of 10 −3 * 10 −4 %. Additional types of simple potentials, which yield time-dependent slow-roll parameters were also studied. In all cases the relative error between calculated results and the traditional SL expression (16) is bounded from above by *0.1%. Furthermore, a more careful analytical treatment leads to better accuracy, bounded from above by about 0.02% relative error. Additionally, we were able to recover the "Cosmic ring" phenomenon, that is the PPS response to a step function in the potential. This response feature in the PPS was first studied in [48].
We take all these results as a strong indication of sufficient accuracy of our calculations.
The embedded panel demonstrates the high accuracy level, as the initial typical error is of order 10 −2 * 10 −3 % and it decreases as a function of the power law index p.