Second-order corrections to Starobinsky inflation

Higher-order theories of gravity are extensions to general relativity (GR) motivated mainly by high-energy physics searching for GR ultraviolet completeness. They are characterized by the inclusion of correction terms in the Einstein-Hilbert action that leads to higher-order field equations. In this paper, we propose investigating inflation due to the GR extension built with all correction terms up to the second-order involving only the scalar curvature $R$, namely, $R^{2}$, $R^{3}$, $R\square R$. We investigate inflation within the Friedmann cosmological background, where we study the phase space of the model, as well as explore inflation in slow-roll leading-order. Furthermore, we describe the evolution of scalar perturbations and properly establish the curvature perturbation. Finally, we confront the proposed model with recent observations from Planck, BICEP3/Keck, and BAO data.


I. INTRODUCTION
Despite the immense predictive power of general relativity (GR), extensions to it have been motivated by several areas.In high-energy physics, which aims for the ultraviolet completeness of GR, quantum gravity and inflation models are included.On the other hand, models involving low-energy physics include, among others, the phenomenology of the dark sector of the universe and spherically symmetric solutions in a weak-field regime.
According to Lovelock's theorem, fundamentally, GR is constructed based on some hypotheses: it is a 4-dimensional Riemannian metric gravity theory, containing the metric g µν as the only fundamental field, invariant by diffeomorphism and with second-order field equations.In this sense, extensions to GR are achieved by violating any of these hypotheses [1].By violating the first hypothesis, we can allow a higherdimensional spacetime or even consider a gravitational action constructed with curvature and torsion invariants due to a Riemann-Cartan spacetime [2].If we violate the hypothesis that the theory of gravity has the metric as the only fundamental field, we can obtain, for example, the Horndeski theories.These, in turn, are the more general 4-dimensional theories of gravity whose action, constructed with the metric and a scalar field, leads to second-order field equations [3].On the other hand, by allowing field equations above the secondorder and preserving all other assumptions, we find the higherorder gravities.
Higher-order theories of gravity are characterized by the inclusion of correction terms in the Einstein-Hilbert (EH) action that lead to higher-order field equations.Such corrections can be conveniently classified according to their mass (energy) scale.In this scenario, EH plus the cosmological constant represents the usual zero-order term.First-order corrections to EH are fourth mass terms constructed from the 4 possible invariants1 R2 , R µν R µν , R µναβ R µναβ and □R.
In turn, the second-order corrections to EH are sixth mass terms, built with the invariants 2   R□R, R µν □R µν , And so on, we will have more higher-order correction terms as we increase the energy scales.Models involving higher-order gravities have been explored in various contexts.There are papers in the literature whose purpose is to show the equivalence between different classes of gravity theories, in particular, between f (R) or f R, □ k R and scalar-tensor theories [4][5][6][7][8][9][10].In some contexts, it becomes more convenient to pass from the original frame to the Jordan or Einstein frames, through a conformal transformation, in order to handle equations for scalar fields rather than higher-order equations for the metric.Another topic of great interest is the investigation of spherically symmetric and static solutions in higher-order gravities, with Stelle's paper [11] being one of those responsible for shedding light on this line of research.In particular, the study of the possibility of non-Schwarzschild black hole solutions through different approaches is addressed in Refs.[12][13][14][15][16][17][18][19], whereas researches involving weak-field regime solutions are covered in Refs.[19][20][21][22].There are also models that study the generation and properties of gravitational waves [23][24][25][26][27][28][29][30][31][32][33].The latter is a topic of great current appeal due to the direct detections [34][35][36] that allow the rising of gravitational wave astrophysics.
In this paper, we propose to investigate the extension to the Starobinsky model due to the inclusion of all correction terms up to the second-order involving only the scalar curvature R. In this sense, we have the following gravitational action where κ 0 has squared mass unit and parameters α 0 and β 0 are dimensionless quantities.Furthermore, M P l is the reduced Planck mass, such that M 2 P l ≡ (8πG) −1 and □ ≡ ∇ σ ∇ σ represents the covariant d'Alembertian operator.In this scenario, where we only address the scalar sector of corrections, R 2 represents the first-order correction, while the last two terms correspond to the second-order corrections to EH.Note that the parameter κ 0 is responsible for establishing the energy scale of inflation, while the parameters α 0 and β 0 give us a measure of the Starobinsky deviation.Since R3 and R□R are both second-order correction terms on energy scales, they must contribute similarly to inflation, so there is a joint effect that must be considered.In that regard, it is worth noting that our paper goes a step further in recent researches developed in [44] and [58,59], which address the models Starobinsky+R 3 and Starobinsky+R□R, respectively.In this paper, the multifield treatment associated with the R□R term is different from that used in Ref. [58].While in that paper, inflation is described by a scalar and a vector field, here, inflation is driven through the dynamics of two scalar fields.Furthermore, by properly constructing the curvature perturbation, we can obtain observational constraints different from those obtained in Ref. [58] for the tensor-to-scalar ratio.In turn, by assuming the R□R sixth-derivative term as a small perturbation to Starobinsky inflation, Ref. [59] uses a somewhat different approach, being able to map the model into a one-scalar theory.
It is important to comment that the discussed model ( 1) is not seen as a fundamental theory of gravity.On the other hand, it is seen as a classical model of gravity in a context of effec-tive theory.One could legitimately worry about the ghost-type instabilities introduced with the R□R sixth-derivative term 3 .Nevertheles, as previously pointed out by Refs.[65,66], the complications of the growing up explosive behaviour of the ghost-type perturbations will not take place only until the initial seeds of such perturbations do not have sufficiently high frequencies.Usually, as long as the energy scales involved are close to the Planck order of magnitude, cosmological solutions are stable.
The paper is structured as follows.In section II, we start from the original frame for the action (1) and rewrite it in the scalar-tensor representation in the Einstein frame, where the theory is described through a metric and two auxiliary scalar fields, only one of which is associated with a canonical kinetic term.Then we write the field equations for each of the fields.Section III is responsible for making the full description of inflation in the cosmological background.In section III A, we study the critical points and the 4-dimensional phase space of the model.Next, we explore inflation in the slow-roll leading order regime by defining the slow-roll factor, and thus, we obtain the slow-roll parameters and the number of e-folds.In section IV, we give a complete description of the evolution of scalar perturbations.In addition to writing the perturbed field equations in the slow-roll leading order regime, we define the adiabatic and isocurvature perturbations by separating of the background phase space trajectories in the tangent (adiabatic perturbation) and orthogonal (isocurvature perturbation) directions.This allows us to properly establish the curvature perturbation, which is essential to connect our model with the observations.In section V, we confront the proposed model with the recent observations of Ref. [40], where by using the constraint for the number of inflation e-folds found in [44], we build the usual n s ×r 0.002 plane and the Plot for the parameter space α 0 × β 0 .In section VI, we make some final comments.

II. FIELD EQUATIONS
The first step is to rewrite the action (1) in the Einstein frame.Performing this calculation, we get with the potential associated with the model.The quantities with bar are defined from the metric as ḡµν = e χ g µν and the di-mensionless fields χ and λ are defined as where in the Einstein frame, □λ = e χ □λ − ∂µ λ ∂µ χ .
Addendum By recovering the usual notation and the dimensions of the scalar fields, and the potential, we must take This way, we can rewrite the action (2) as By starting from the action (2), we obtain three field equations: one for ḡµν and another two for each of the scalar fields χ and λ.Taking the variation concerning the metric ḡµν , we find where we define an effective energy-momentum tensor as The variation concerning the χ and λ fields results in where V χ = ∂ χ V and V λ = ∂ λ V represent derivatives concerning the fields χ and λ, respectively.

III. INFLATION IN FRIEDMANN COSMOLOGICAL BACKGROUND
On large scales (≳ 100 Mpc), we can consider the universe to be homogeneous and isotropic.Furthermore, for a spatially flat universe, the line element that describes the evolution of a comoving frame of reference is given by where a (t) is the scale factor.By obtaining the field equations in Friedmann background is to write the field equations ( 7), ( 9) and (10) for the metric (11).From the field equation for the metric, we get two independent ones, namely the Friedmann equations where H = ȧ/a.In addition to these equations, we also have the equations for the χ and λ fields.Since, for a scalar field Φ, for the equation of χ given in (9), we have In turn, for the equation of λ given in (10), we have A. Phase space In this section, we will analyze the phase space of the model.Therefore, it becomes convenient to rewrite the field equations in a dimensionless way: we define the dimensionless time derivative and the dimensionless potential V as With that, it is possible to rewrite the equations of cosmological dynamics ( 12), ( 13), ( 14) and ( 15) as follows: and We already know the inflationary dynamics of Starobinsky+R 3 model, which in the scalar-tensor approach in the Einstein frame is characterized by its specific potential V (χ) [44], as well as the dynamics inflation of Starobinsky+R□R model, explored in Ref. [58] via a scalar-vector approach.A first step in order to understand the dynamics of our current case is through the study of its phase space, having as reference the known particular cases mentioned above.In this first part, we will investigate the existence of an attracting inflationary regime in some region of the phase space.
Since the dimensionless equations governing the dynamics of the χ and λ fields are written as in (18) and (19), that is, two autonomous second-order differential equations concerning time, we can rewrite them as a system of four first-order differential equations.Taking χ t = ψ and λ t = ϕ, we have where which is associated with a physically consistent system when its root argument is positive 4 .
From that point on, we will study the approximate behavior of the solutions of the system at critical points.Critical points are equilibrium points of the system, and it is our interest to investigate their stability, which is directly related to the necessary conditions for the occurrence of a physical inflationary regime 5 .The analysis of the previous system allows us to conclude that there are two critical points: The study on the stability of these critical points is done through the linearization of the 4-dimensional autonomous system (χ, λ, ψ, ϕ).Linearizing the system given by Eqs. ( 20), ( 21), ( 22) and ( 23) around P 0 , we verify that the Lyapunov exponents r 0 , associated with the stability of the critical point, satisfy the fourth-order characteristic equation whose solution is A center or spiral point occurs when we obtain pure imaginary roots.Looking at the previous expression, we see that this occurs whenever the condition is satisfied.Any value of β 0 outside this range contains at least one Lyapunov exponent with positive real part.That is, outside the range (27) the point P 0 is unstable. 6A numerical analysis of the system (21) shows that within the interval (27) the point P 0 is an attracting spiral point and therefore stable (see figure 1).This behavior is essential for the existence of a graceful exit.In fact, the spiral dynamics around P 0 constitute the period of coherent oscillations consistent with the initial phases of reheating.It is also worth noting that Eq. ( 26) is independent of α 0 , and therefore the term R 3 plays no role at the end of the inflationary period.
In turn, linearizing the system (21) around P c , we verify that the Lyapunov exponents r c satisfy the characteristic fourth-order equation where .
A numerical study of this characteristic equation, considering α 0 > 0 and β 0 > 0, shows that at least two of the four roots of Eq. ( 28) are real and have opposite signs.This shows that P c is a saddle point and therefore unstable.This conclusion also remains valid for β 0 = 0 and α 0 > 0, in which case we have only two roots. 7See Ref. [44] for details.
To better understand the dynamics of the χ and λ fields, we will numerically study the 4-dimensional phase space.In this study, we will analyze two 2-dimensional slices of this space given by χ t × χ and λ t × λ.For that, we manipulate Eqs. ( 18) and ( 19) writing them as where h is given by ( 16).
Numerical analysis of the equation dχ t /dχ is more easily performed if we write λ = λ (χ, χ t , λ t , λ tt , α 0 , β 0 ).For that, it is necessary to work with the equations ( 19) and ( 16).Solving the quadratic equation for λ in Eq. ( 19), we get where we choose the positive sign to guarantee the Starobinsky limit.In principle, we can substitute (16) in the previous expression, obtain a third-degree algebraic equation for λ and solve it to obtain λ = λ (χ, χ t , λ t , λ tt , α 0 , β 0 ).However, we will see in Sec.V that the values of interest for α 0 and β 0 are such that α 0 < 10 −3 and β 0 < 3 × 10 −2 . 8In this case, it is licit to consider only linearized corrections of α 0 and disregard terms of the type α 0 β 0 .Performing these approximations, we obtain the functional forms where with In figures 1 and 2, we show direction fields associated with equations ( 29) and (30).
The first (and most relevant) point that can be seen in figure 1 is that there is an attractor line close to χ t ≃ 0. The existence of this region is consistent with any value of α 0 < 10 −3 and β 0 < 3 × 10 −2 and for any interval of λ t and λ tt that yields real results in the region of interest χ ∈ [0, 8]. 9 At the same time that the χ field tends to the attracting line (χ t ≃ 0), figure 2 indicates that λ tends to a finite value and λ t → 0. This finite value of λ essentially depends on the value of χ with variations on a smaller scale due to changes in the parameter α 0 .The other fixed parameters χ t and β 0 in figure 2 change how λ approaches the accumulation point but does not change its value.We will see in the Sec.III B that this attractor region in the 4-dimensional phase space where (χ, λ, χ t ,λ t ) ≃ (χ, λ (χ) , 0, 0) corresponds to a slow-roll inflationary regime.
Once the attractor region is reached, we must ask ourselves if inflation occurs enough, i.e., if it generates a sufficient number of e-folds and if it ends in a reheating phase.The answer to this question essentially depends on the position where the χ field hits the attractor line in figure 1.If the χ field is to the 8 See also Refs.[44] and [58]. 9These ranges are typically between −10 e 10. left of the critical point P c (black dots in the graphs of figure 1), inflation proceeds normally and ends in a phase of coherent oscillations associated with the beginning of reheating.On the other hand, if χ is to the right of P c the value of χ increases indefinitely, and inflation never ends (see Ref. [44] for details).Thus, a physical inflationary regime, i.e., consistent with a graceful exit, only occurs if χ < χ c ⇒ α 0 < 3 (e χ − 4) −2 which for sufficiently large χ corresponds to α 0 < 3e −2χ .
In the next section, we will see how to describe the dynamics of the χ and λ fields during the slow-roll inflationary phase.

B. Inflation in the slow-roll leading order regime
This section aims to describe the dynamics of χ, λ, and their derivatives during the inflationary regime considering the slow-roll approximation.In the region associated with physical inflation, the parameter χ is a monotonic decreasing function of time, so we can parameterize the various quantities in terms of χ.For the case of Starobinsky model, we know that in the slow-roll leading order regime χ t ∼ δ and χ tt ∼ δ 2 , where δ is the slow-roll factor defined as δ ≡ e −χ .And for Starobinsky plus R 3 model (i.e., β 0 = 0 and α 0 ̸ = 0), we have [44] Note that since α 0 < 3δ 2 , the second term of the previous expression is of the same order or less than δ.The χt × χ graphs considering phase space cuts (χ, λ, χt, λt) fixing λt = λtt = 0 and β0 = 0.001 with (λ, α0) = (173, 0.0001) (top graph) and (λ, α0) = (94, 0.00034) (bottom graph).The red and black points correspond to the critical points P0 and Pc, respectively.For α0 = 0.0001, we have Pc = (5.18,173, 0, 0) and for α0 = 0.00034, we have Pc = (4.58,94, 0, 0).The red (cyan) trajectories represent trajectories that, when reaching the attractor line close to χ = 0, approach (depart) from the origin.Details on the interpretation of the graphics are presented in the body of the text.
The previous discussion allows us to associate the factor δ as a parameter that controls the slow-roll approximation order, i.e., a quantity f ∼ δ n will be an nth-order slow-roll quantity.In this case, χ t present in (33) is first-order in slow-roll, since both δ and α 0 δ −1 are first-order.To apply this reasoning in our model, it is also necessary to establish what is the maximum slow-roll order of the parameter β 0 .A more detailed analysis of the field equations in the attractor region shows us that, for slow-roll inflation, β 0 ≲ δ, i.e., β 0 is a (at most) first-order slow-roll parameter (for details see Ref. [58]).
By following similar reasoning, we propose the following ansatz for λ: In this case, the first term is of order O (−1) in slow-roll, and the others are zero-order terms.The O (−1) order term is necessary as we know that in the case of Starobinsky λ = δ −1 − 1 (see Eq. ( 19) with α 0 = β 0 = 0).Analogously to χ t , each derivative of λ with respect to t increases by one the slow-roll order.
The next step is substituting these two ansatzes and their derivatives into Eqs.( 16), (18), and (19) taking into account only the slow-roll leading order.In this situation, we get where By explicitly substituting Eqs. ( 34) and ( 35) in these last three expressions, we get after a long calculation For details see appendix A. It is worth noting that the previous expressions are well defined only for β 0 δ −1 < 3. Note in Eq. ( 39) the existence of two terms that, in the slow-roll leading order, are first-order terms.In Eq. ( 40), we have the O (−1) order term in addition to the zero-order corrections.Finally, h 2 , related to the Hubble parameter, is given by the zero-order slow-roll leading term plus first-order corrections (independent of β 0 ).

Calculation of slow-roll parameters and number of e-folds
The characterization of the inflationary regime is done through the slow-roll parameters By substituting (39) and ( 40) in (17), we get Thus, in the slow-roll leading order, we have The next step is calculating η.Differentiating ϵ and using this result together with ϵ itself in Eq. ( 43), we get Note that by construction α 0 δ −2 < 3 and β 0 δ −1 < 3.
In order to have robust inflation, i.e., with enough number of e-folds, we must have ϵ ≪ 1 and η ≪ 1.Thus, from the equations ( 44) and ( 45), we see that this occurs for δ ≪ 1 (typically χ ≳ 4).However, unlike the Starobinsky case, we also have lower bounds for δ.In fact, the slow-roll inflationary regime only occurs if The first condition does not represent a real difficulty for the existence of slow-roll inflation, because even if at some point we have δ < β 0 ,10 the dynamics of the phase space guarantees that χ decreases monotonically so that at some point δ becomes greater than β 0 (see figure 1).The second condition represents a real constraint for carrying out a physical inflation (see discussion at the end of Sec.III A).For a discussion of the implications of this second constraint and the initial conditions of inflation see Ref. [44].Next we will calculate the number of e-folds N in the slowroll leading order.By the definition of N , we have where the index e corresponds to the end of inflation.To integrate this expression, it is convenient to perform the following change of variable: In this case, we get whose solution is By considering only leading terms and taking into account that x e ≪ x, we finally get By construction, physical inflation occurs in the interval 0 ≤ x < 1.In fact, when x → 1 we have δ → δ m which corresponds approximately to χ → χ c (see eq. ( 25)).However, the expression (48) has an extra restriction due to the presence of the β 0 term.For β 2 0 > 3α 0 ⇒ γ > 1, we have that when x → 1, the value of N diverges to −∞, and this is clearly not physical.What happens is that for γ > 1, the function N presents a maximum point within the interval 0 ≤ x < 1. Differentiating N with respect to time, we have So, for γ > 1, we have On the other hand, for values of x such that x max ≤ x < 1, we get which violates the first condition of Eq. (46).Therefore, based on the previous analysis, we conclude that Eqs.(44), (45) and (48) referring to the quantities ϵ, η and N are valid in the following intervals: In the next section, we will study the inflationary regime from the perturbative point of view.

IV. INFLATION VIA COSMOLOGICAL PERTURBATION THEORY
In this section, we investigate inflation of the model (2) via cosmological perturbations.Recall that its background dynamic equations are the Friedmann ones, given by Eqs. ( 12) and ( 13), and the equations of motion for the scalar fields χ and λ, given by Eqs. ( 14) and (15).
Before proceeding with our developments, it is worth commenting on scalar perturbations.In addition to the perturbations of the two scalar fields, which we will denote by δχ and δλ, we have the scalar perturbations of the metric.The line element in the perturbed Friedmann-Lemaître-Robertson-Walker (FLRW) metric is given by with A, B, ψ and E being the scalar perturbations of the metric [67,68].In order to obtain the perturbative field equations through a perturbation directly in the action, we need to write it up to the second order in the perturbations.In this case, we must consider second-order terms for perturbations involving only scalar field perturbations (e.g., δχ 2 ), second-order terms involving only metric scalar perturbations (e.g., A 2 ) and cross terms, that is, a product of first-order terms (e.g., Aδχ).In the following subsection, by following a perturbative procedure directly in the action, along the lines of that found in Refs.
[ [68][69][70], and assuming the spatially flat gauge, 11 we obtain and discuss the equations of motion for the perturbations.

A. Equations for scalar perturbations
The first step in order to perturbate the action is defining the perturbations of the scalar fields.For an inhomogeneous distribution of matter, we write In turn, the metric in Eq. ( 51) is written as By writing Eq. ( 2) up to the second order in the perturbations and taking their variations with respect to each one of the perturbations, we are able to obtain the following equations of motion for the perturbations δχ and δλ and as well as the Einstein equations and The double subscript in potential V represents second-order differentiation with respect to the corresponding scalar fields.

B. Equations in the slow-roll leading order regime
Once we obtain Eqs. ( 55), ( 56), ( 57) and ( 58), which completely describe the evolution of scalar perturbations, the next step is to write them in the slow-roll leading order regime.This is not a trivial task and therefore, we will initially present the particular case of the Starobinsky model (α 0 = β 0 = 0).In this case, Eqs. ( 55), ( 56), ( 57) and ( 58) reduce to and with In Sec.III B, we saw, in the context of the background, the behavior of the scalar fields, their derivatives and the relationships they keep between them.Once the slow-roll factor δ was established, we recall that in the slow-roll leading order regime, we obtain and that with each differentiation with respect to time in the scalar fields, an order of slow-roll is increased, that is, χ ∼ δ 2 and λ ∼ δ 0 .By making the constructions in this section, some assumption is necessary, namely, to find the slow-roll orders of the perturbations, we need to establish the slow-roll order of one of them.In this sense, we take the perturbation δχ as a zero-order slow-roll quantity.Also, we keep in mind that derivatives do not change the slow-roll order of perturbations.We are now able to write the equations of motion in the slow-roll leading order.Analyzing Eq. ( 62), note that since χδχ ∼ δ, perturbation A must be at most first-order in slowroll.Regarding Eq. ( 61), in its right member, we have the first and third terms, which are of first-order, the second term, which is subdominant 3 χ2 A of third-order in slow-roll, and the last one is null, since Vλ = 0. Furthermore, since in its left member we have 3H 2 A ∼ δ, we conclude that perturbation B is at most first-order in slow-roll.In turn, as Vχλ ∼ δ and Vλλ ∼ δ 2 , we see that Eq. ( 60) establish the perturbation leading order of δλ, namely, δλ ∼ δ −1 .With Eq. ( 60), we can still write δλ in terms of δχ and substitute in Eq. ( 59).Thus, on the left side of Eq. ( 59), the first three terms are of zero order in slow-roll, while all other terms of the equation give us subdominant contributions.So we can write When developing the previous analysis now for the case of the complete equations, we find that the perturbations evolve, in the slow-roll leading order regime, with the same orders obtained previously.In short, the scalar perturbations evolve in the form A ∼ B ∼ δ and δλ ∼ δ −1 .It is interesting to note that the perturbation δλ goes in leading order regime with δ −1 , and that if it were otherwise, it would seriously compromise the slow-roll dynamics.
By applying all the discussion raised above, we find, in the slow-roll leading order regime, the following equations of motion for the perturbations of the scalar fields and These results are in agreement with those obtained in Ref. [58], where the Starobinsky+R□R model is explored.

C. Adiabatic and isocurvature perturbations
In this subsection, we define adiabatic and isocurvature perturbations, obtain expressions that describe their dynamics, and study their solutions.
The action (2) can be rewritten, along the lines of Ref. [71], compactly as where the scalars Φ I (x) are seen as local coordinates of the scalar field space with metric G IJ (Φ) and V represents the potential of the model, Eq. ( 3).In a twofield scalar model, the field space is 2-dimensional and characterized by G IJ (Φ).To conveniently describe the evolution of perturbations, we can define a basis having a tangent direction, which we will denote by σI , and another orthogonal, ŝI , to the background trajectories.Tangent directions to background trajectories are associated with adiabatic perturbation, while orthogonal directions are associated with isocurvature perturbation.In this sense, we build the basis through the definitions, respectively, of the module of the velocity vector, the unit velocity vector in the tangent direction and the normalization rule and for the orthogonal direction, the normalization 12 and orthogonality rules For our case, we have for the velocity module σ σ = 3 χ2 − β 0 e −χ λ2 .
Note that it is directly related to the Friedmann equation ( 13).
In turn, for the unit vectors in the orthogonal direction to the background trajectories, we have ŝχ = β 0 e −χ 3 λ σ and ŝλ = 3 Continuing our study on the evolution of scalar perturbations, we point out that the quantity δΦ I g given by is gauge invariant.It turns out that ψ = 0 when working on a spatially flat gauge, so that δΦ I g = δΦ I .That said, by projecting δΦ I in the σ and ŝ directions, we construct the adiabatic Q σ and isocurvature Q s perturbations, respectively.In that sense, we have and It is worth noting that from the point of view of the slowroll approximation both Eq. ( 75) and Eq. ( 76) are zero-order.By having written the expressions for the adiabatic and isocurvature perturbations, the next step is to invert the relations in order to obtain δΦ I = δΦ I (Q).We obtain these relations by solving the linear system given by Eqs. ( 75) and (76).Thus, we find the expressions and Once the field perturbations in terms of the adiabatic and isocurvature perturbations were obtained, we can write the second order perturbed action for Q σ and Q s .Such an action is fundamentally constituted by quadratic terms involving Q σ and Q s (e.g., Q 2 σ ), cross terms involving a Q σ or Q s and a metric perturbation (e.g., AQ σ ), and quadratic terms for metric perturbations (e.g., A 2 ).On the other hand, we can express the cross term only in terms of the perturbations Q σ and Q s by making use of the constraints from the Einstein equations, Eqs. ( 57) and (58).By taking these considerations into account, we find the following structure for the part of the action that depends only on Q σ and Q s where the coefficient of the cross kinetic term ∂ κ Q σ ∂ κ Q s is zero, and the others can be found in Appendix B. It is interesting to analyze the behavior of kinetic terms and their possible contribution to the emergence of ghost-type instabilities.Note that the kinetic terms are canonical, equal and with reversed signs.This characteristic irremediably indicates that the existence of ghost-type instabilities is something intrinsic to the model and that it is essential to take this into account when performing the perturbation quantization process. 13Furthermore, the fact of the non-existence of the cross kinetic term is something expected and is directly related to our approach of making a consistent decomposition of the perturbations in the tangent and orthogonal directions to the trajectories of the background phase space.
When writing Eq. ( 79) considering a slow-roll leading order regime, we obtain We now turn our attention to the task of writing the equations for the evolution of adiabatic and isocurvature perturbations.This is done by substituting the expressions (77) and (78) in the dynamic equations ( 65) and (66).In this case, taking the first derivatives of Eqs. ( 77) and ( 78), remembering that the background quantities can be considered constant, we are able to write and where Note that relations (81) and ( 82) indicate that the adiabatic Q σ and isocurvature Q s perturbations are decoupled, that is, they evolve independently in our model.That is an interesting result since such a decoupling usually does not occur.Generally, the isocurvature perturbation enter as source of the adiabatic one.[71].
From this point on, it becomes convenient to treat the field equations for the perturbations in a Mukhanov-Sasaki form.By making a redefinition of the perturbations and assuming conformal time and Fourier space, we get with the prime representing derivative with respect to conformal time and where In a de Sitter background (slow-roll zero order), we have Furthermore, in the slow-roll zero order regime, we have For consistency with several previous results, we have β 0 e χ < 3 so that we can define a quantity and in this way we write the expressions Next, we will explore the solutions of Eqs. ( 90) and (91).

D. Solutions to the perturbations
Once the equations for the dynamics of adiabatic and isocurvature perturbations have been established in the appropriate form, given by Eqs. ( 90) and (91), we can write and analyze their solutions.
In a subhorizon regime, kη ≫ 1, equations to the perturbations are approximated by The quantization process in de Sitter takes the following initial conditions [72,73] or Note that due to the ghost-type behavior of isocurvature perturbation, the quantization of the Q s field was performed in the same way as in Refs.[48,74].In principle, this behavior can raise questions about the unitarity of the theory [75] (see also the discussions in Refs.[65,66]).However, as we will see below, the Q s field decays rapidly after crossing the horizon, suppressing any observable effects associated with isocurvature perturbation. 14he exact solution of Eqs. ( 90) and (91) can be written using a combination of Hankel's functions as [73,76] where for the adiabatic perturbation δφ σ , we have ν σ = 3/2, and for the isocurvature perturbation δφ s , we have To determine the constants, we compare the general solution with the initial conditions in Eq. (92).For the adiabatic case (ν σ = 3/2) in the subhorizon regime, we find so that In turn, for the case of isocurvature perturbation, where ν s is given by Eq. ( 95), for a subhorizon regime, we get15 Thus, the solution to the isocurvature perturbation is written as Since the solutions (97) and (99) for the perturbations have been established, we can analyze their behavior in the superhorizon regime (kη ≪ 1).Taking into account Eq. (87), which gives us the relation between η and a in a de Sitter background, for adiabatic perturbation, we find while for isocurvature perturbation, we have The result in Eq. (100) tells us that the adiabatic perturbation is constant in the superhorizon limit, whereas Eq. (102) reveals a decaying behavior for the isocurvature one.Remembering that M > 0, we have the following situations: • If the quantity ν s is real, we have representing a solution to Q s that decays with a.
• If the quantity ν s is imaginary, also providing a decaying solution, since we have the product of an oscillatory term and the term decaying with a − 3 2 .
We can provide a quantitative measure of isocurvature perturbation by considering scales of interest during inflation (measured in CMB anisotropies).These ones are within the range 10 −3 M pc −1 < k < 10 4 M pc −1 , where the pivot scale is k * = 0.002 with 50 < N * < 60.The point is that during the inflationary regime, a given scale k crosses the horizon at a specific value of the number of e-folds N .The smaller/larger the scale k is, the smaller/larger the number of e-folds N it experiences after crossing the horizon.Taking the smallest scale 16 k sm = 10 4 M pc −1 , we find N sm = N * − 15.4.In this sense, for N * = 50 and ν s = 1/2, we get which shows us that isocurvature perturbation is negligible after inflation.Since, in addition to this, they do not enter as a source of adiabatic perturbation, we can consider them negligible.All the previous analysis was carried out considering the slow-roll approximation.
In the next section, we will analyze the connection of our model with the observations.

V. OBSERVATIONAL CONSTRAINTS
At the end of the previous section, we show why isocurvature perturbation is negligible after inflation.Furthermore, since adiabatic perturbation has the same behavior as in the case of a single scalar field, we easily recognize the power spectrum and its connection with observational parameters.To make the connection with the observations, specifically to write the power spectrum, it is interesting to recover the mass units of the fields.In this sense, equations such as ( 12), ( 13), ( 69) and (75) need to be written in terms of massive fields given in Eq. ( 5).Since the curvature perturbation is given by 16 Since all others decay further. 17the power spectrum of adiabatic perturbation is written as where we evaluate it to k = Ha at the instant when k crosses the horizon.The result in Eq. ( 106) is identical to the power spectrum for single-field inflationary models.Thus, in the slow-roll leading order regime, the scalar spectral index n s and the tensor-to-scalar ratio r are, respectively, where ϵ and η are the slow-roll parameters of the model given by the Eqs.( 44) and ( 45).These equations depend on α 0 , β 0 , and the number of e-folds N through Eq. ( 48) which carries the dependency between N and δ.
In our paper, there are two types of Plots where we compare our model with observational data [40], built from Eq. (107) taking the three independent parameters α 0 , β 0 , and N : the usual n s × r 0.002 plane and the parameter space α 0 × β 0 .The Plots are constructed by setting one of the parameters and varying the others.We use the range 52 ≤ N ≤ 59 for the number of inflation e-folds N based on a reheating modeling.For details, see appendix C.
The figure 3 shows the n s × r 0.002 plane containing the observational constraints (in blue) obtained from Ref. [40] and the theoretical evolution of the model in two different situations.
In the top graph of figure 3, we fixed the parameter β 0 and varied the others.In it, the light red region represents Starobinsky+R 3 model, which starts at the light red points.In turn, the light yellow region represents the complete model with β 0 = 1.5 × 10 −2 , starting at the yellow points.As we increase the values of the parameter α 0 , the region predicted by the model shifts to the left and slightly downwards, until it crosses the region of 95% C.L..This behavior can also be seen in Ref. [44], and it is consistent with the results obtained in Ref. [42], where β 0 = 0.These constraints establish, in the most conservative way, a maximum value for α 0 ∼ 10 −4 .
In the bottom graph of figure 3, on the other hand, we fixed the parameter α 0 and varied the others.The light red region represents the Starobinsky+R□R model, which starts at the light red points.In turn, the light green region represents the complete model with α 0 = 10 −5 , starting at the light green points.As the values of β 0 increase, the region predicted by the model moves to the right and slightly upwards, until it crosses the region of 95% C.L..These constraints establish a maximum value for β 0 ∼ 10 −2 .Similar results were obtained in Refs.[58,59] for the Starobinsky+R□R case.However, a considerable difference between our results and those in Ref. [58] is checked for the constraint on the tensor-to-scalar ratio r 0.002 .There, r 0.002 can assume larger values, so that the growth of the region predicted by the model is more accentuated.This difference is due to the fact that in Ref. [58], the definition for the curvature perturbation was not established properly by not making a separation of the background phase space trajectories in the tangent (adiabatic perturbation) and orthogonal (isocurvature perturbation) directions.On the other hand, our results are closer to those in Ref. [59], indicating that the approach of treating the R□R term as a small perturbation is relevant and consistent.Another plot developed, figure 4, is the parameter space α 0 ×β 0 allowed by the observations.In the top graph of figure 4, we have the Plot for N = 52, while in the bottom graph of figure 4, we have it for N = 59.The blue regions represent the allowed regions for the α 0 and β 0 parameters in 68% and 95% C.L.. Note that the Plot for N = 52 gives us a smaller region for the parameters if we compare it to the Plot for N = 59.In addition, we can see two regions on each of the Plots.One is an approximated rectangular region completely within the 95% C.L. (for N = 52, the sides correspond to α 0 = 2.8× 10 −5 and β 0 = 1.7×10 −2 , and for N = 59, α 0 = 5.3×10 −5 and β 0 = 1.5 × 10 −2 ), where the parameters α 0 and β 0 do not keep a dependency between them, being able to assume any values independently.In this region of independence between the model parameters, we reproduced the results obtained in Refs.[42,59] for each model separately.The other region is the asymptotic one for large values of α 0 and β 0 , whose occurrence suggests a dependence β 0 = β 0 (α 0 ) between the parameters.

VI. FINAL COMMENTS
The Starobinsky model is one of the most competitive candidates for describing physical inflation.In addition to having a well-grounded theoretical motivation, it better fits the recent observations [39,40].Motivated by the success of such a model, we propose to investigate inflation based on the higherorder gravitational action characterized by the inclusion of all terms up to the second-order correction involving only the scalar curvature, namely, the terms R 2 , R 3 , and R□R.In this sense, our proposed model has two additional dimensionless parameters, α 0 , and β 0 , whose values represent deviations from Starobinsky.
Unlike Ref. [58], whose multi-field treatment used to address the term R□R gives us an inflation described by a scalar and a vector field, here, when passing from the original frame to the representation in the Einstein frame, the model is described through the dynamics of two scalar fields χ and λ, where only one of them is associated with a canonical kinetic term, and whose potential is V (χ, λ) given in Eq. ( 3).The study of inflation in a Friedmann background, through the analysis of the critical points and phase space of the model, is essential to verify the existence of an attractor region associated with the occurrence of an inflationary regime and to know if such a regime has a graceful exit.We took as a basis the study of particular cases developed in Refs.[44,58], which deal with the Starobinsky+R□R and Starobinsky+R 3 extensions.We saw that there is an attractor line near χ t ≃ 0, corresponding to the slow-roll inflation, for any value of α 0 < 10 −3 and β 0 < 3 × 10 −2 .Furthermore, the occurrence of such a physical inflation regime essentially depends on the initial conditions for the χ field.If they are such that the χ field is to the right of the critical point P c , the value of χ increases indefinitely, and inflation never ends.On the other hand, the occurrence of a consistent physical inflationary regime that has a graceful exit essentially requires that the initial conditions be such that χ < χ c , i.e., that it is to the left of the critical point P c .Finally, we conclude the background analysis with the study of inflation considering the slow-roll approximation.By defining the slow-roll factor δ, which in our analysis is responsible for controlling the slow-roll approximation order, we obtain all relevant quantities, such as ε and η, in the slow-roll leading order.
There is considerable literature about multi-field inflation models, which we took into account to develop the analysis at the perturbative level [68,70].The equations of motion for the scalar perturbations were obtained using the spatially flat gauge.By writing the equations in the slow-roll leading order approximation, we saw that the scalar perturbations of the metric are sub-dominant concerning the perturbations δχ and δλ.At this point, we performed a correct decomposition of the perturbations in the tangent (adiabatic perturbations) and orthogonal (isocurvature perturbations) directions to the phase space background trajectories.This way, adiabatic Q σ and isocurvature Q s perturbations are completely separated.Such a decomposition allows us to consistently establish the curvature perturbation, which led us to obtain observational constraints different from those obtained in Ref. [58].The action written in terms of Q σ and Q s makes it clear that there are irremediably ghost-type instabilities in the model since the kinetic terms have opposite signs.Next, we write the equations of motion in a Mukhanov-Sasaki form in order to study their solutions.We obtained the exact solutions for the perturbations through a linear combination of the Hankel functions.Their analysis leads us to conclude that the isocurvature perturbation associated with the ghost field is negligible after inflation and that the adiabatic one has the same behavior as in the case of a single-field inflation.All previous results were obtained considering the slow-roll approximation.Thus, a question that remains is whether the suppression of isocurvature perturbation holds beyond the slow-roll regime.This issue will be addressed in a further work.
Finally, we confront our model with recent observations from the Planck satellite, BICEP3/Keck and BAO [39,40], making use of a constraint on the number of e-folds N of inflation (52 ≤ N ≤ 59) based on reheating modeling [44].For that, we made two types of Plots, namely, the usual n s ×r 0.002 plane and the parameter space α 0 × β 0 .In this analysis, we have three parameters: α 0 , β 0 , and N .Thus, to build the Plots, we set one of the parameters and vary the others.Fixing the parameter β 0 , we observe that the region predicted by the model in the n s × r 0.002 plane shifts to the left and slightly downwards.On the other hand, fixing the parameter α 0 , we notice that the predicted region shifts to the right and slightly upwards.By setting α 0 = 0, we get the Starobinsky+R□R model.In this context, we saw that inconsistency in establishing the curvature perturbation in Ref. [58] led them to obtain values higher than ours for the tensor-to-scalar ratio.In turn, by fixing the number of e-folds N , we construct the parameter space α 0 × β 0 constrained by the observations.In general, the model predictions are more in agreement with the observations for a number of e-folds N = 59.Our analysis, conservatively, restrict the parameters to maximum values of α 0 ∼ 10 −4 and β 0 ∼ 10 −2 .It is also worth pointing out the behavior of the α 0 × β 0 parameter space.The R 3 and R□R terms are second-order correction terms on energy scales and, therefore, should contribute similarly to inflation.In this sense, the joint effect of such terms is reflected in the plot α 0 × β 0 .In fact, there is a considerable region in which the parameters do not depend on each other and which we can associate with the models separately discussed in Refs.[42,44,58,59].However, there is an asymptotic region for large values of α 0 and β 0 , where such parameters seem to keep a constraint.In this particular region, a change in one of the parameters necessarily implies a change in the other, so that the possibility of a dependence β 0 = β 0 (α 0 ) is something to be investigated.This is a topic that the authors will address in a future research.In Ref. [44], we saw that in the particular case of β 0 = 0, the uncertainty in the number of e-folds N k for the reference scale k = 0.002 M pc −1 , defines the interval 52 ≤ N k ≤ 59. (C1) This result was obtained through a very general modeling of the reheating phase considering that at least the fields of the standard model of particles are present during this phase.
The basic equations of the performed modeling are [44] N re = 4 3 w a − 1 where N re is the number of e-folds of the reheating, T re is the temperature of the reheating, w a is the average of the effective equation of state during the reheating, ρ e is the energy density at the end of inflation, T 0 is the CMB temperature in the present day and g re and g 0 are the relativistic degrees of freedom in reheating and in the present day, respectively.The range (C1) is obtained by imposing the bounds N re ≥ 0 and T re ≥ T (min) re , where T (min) re is determined from the decay of the inflaton field in the matter fields [44].
Based on the previous equations and results, it is relatively simple to conclude that the range obtained in (C1) also applies to β 0 ̸ = 0.The main point is that at the end and after inflation where χ < 1 ⇒ δ ∼ 1, the proposed model behaves essentially like the Starobinsky model. 18Thus, the only term present in Eqs.(C2) and (C3), which may have some relevance when β 0 ̸ = 0 is the term H k .However, by (41), we see that H k weakly depends on β 0 even taking into account slow-roll first-order corrections.With this, we conclude that for our model, it is licit to consider the range for the number of e-folds N as given by Eq. (C1).

Figure 4 .
Figure 4.The regions in blue represent the allowed regions for the parameters α0 and β0 in 68% and 95% C.L., due to observational data from Planck plus BICEP3/Keck plus BAO [40].In the top graph we have the Plot for N = 52, while in the bottom graph we have it for N = 59.Note that the constraints for N = 59 allow a larger region for the parameters α0 and β0 in line with what we saw in the figure 3, whose predictions for N = 59 are more within the region of 68% C.L. Note that for large values of α0 and β0, around α0 = 1.5 × 10 −4 and β0 = 2.5 × 10 −2 (N = 52) and α0 = 2.2 × 10 −4 and β0 = 1.2 × 10 −2 (N = 59), the predicted regions for the parameters converge to an asymptotic region.In this region, the values of α0 and β0 suggest to keep a constraint.