Inflation in a conformally-invariant two-scalar-field theory with an extra $R^2$ term

We explore inflationary cosmology in a theory where there are two scalar fields which non-minimally couple to the Ricci scalar and an additional $R^2$ term, which breaks the conformal invariance. Particularly, we investigate the slow-roll inflation in the case of one dynamical scalar field and that of two dynamical scalar fields. It is explicitly demonstrated that the spectral index of scalar mode of the density perturbations and the tensor-to-scalar ratio can be consistent with the observations acquired by the recent Planck satellite. The graceful exit from the inflationary stage is achieved as in convenient $R^2$ gravity. We also propose the generalization of the model under discussion with three scalar fields.

In addition to inflationary models driven by the scalar filed (i.e., the inflaton field) described above, there have been considered the so-called Starobinsky inflation [3,22] originating from the higher-order curvature term such as R 2 term 2 , where R is the Ricci scalar. This model is observationally supported by the Planck results. Such a theory can be interpreted as a kind of modified gravity theories including F (R) gravity to account for the late-time cosmic acceleration (for reviews on dark energy problem and modified gravity theories, see, for instance, [23][24][25][26][27][28][29][30][31][32][33]). Various inflationary models in modified gravity theories corresponding to extensions of the Starobinsky inflation have been explored in Refs. [34,35].
In this paper, we investigate inflation in a theory consisting of two scalar fields which non-minimally couple to the Ricci scalar and an additional R 2 term 3 . We consider the conformally-invariant two-scalar-field theory in which the conformal invariance is broken by adding an R 2 term. In particular, we explore the slow-roll inflation in the cases of (i) one dynamical scalar field (namely, we set one of two scalar fields a constant) and (ii) two dynamical scalar fields. As a consequence, we analyze the spectral index of scalar mode of the density perturbations and the tensor-to-scalar ratio and compare the theoretical results with the observational data obtained by the recent Planck satellite and the BICEP2 experiment. It is clearly shown that the spectral index and tensor-to-scalar ratio can be compatible with the recent Planck results.
The motivation to propose our theory is to unify inflation in the early universe originating from the R 2 term and the late-time cosmic acceleration, i.e., the dark energy dominated stage with dark matter. The R 2 term is interpreted as the contribution from modified gravity, dark energy is described by one of the scalar fields, and dark matter is represented by the other scalar field. Furthermore, it seems that multiple field inflation models can fit the Planck data better than single field inflation models. Inflationary models with multi-scalar fields [37] including the so-called curvaton scenario [38] have been constructed and the cosmological perturbations in these models have also been investigated [39]. We use units of k B = c = = 1 and express the gravitational constant 8πG N by κ 2 ≡ 8π/M Pl 2 with the Planck mass of M Pl = G −1/2 N = 1.2 × 10 19 GeV. The organization of the paper is the following. In Sec. II, we explain our model action and derive the gravitational field equation and the equations of motions for the scalar fields. We also examine the slow-roll inflation in the case of one dynamical scalar field and study the dynamics of the system including the equilibrium points in detail. In Sec. III, with the conformal transformation [40], we explore inflationary cosmology in the Einstein frame. Especially, we study inflationary models in the case that the conformal scalar is the dynamical inflaton field and other scalar fields are set to be constants. In Sec. IV, we investigate the slow-roll inflation in the case of two dynamical scalar fields. Particularly, we consider the resultant spectral index and the tensor-to-scalar ratio in the Jordan frame (i.e., the original conformal frame). In Sec. V, we explore the graceful exit from inflation, namely, the instability of the de Sitter solution in the present theory. As a demonstration, we concentrate on the case of one dynamical scalar field in the Einstein frame, because, as is described in Sec. III, in this case the spectral index and the tensor-to-scalar ratio can be consistent with the recent Planck results. Finally, conclusions are described in Sec. VI.

A. Model action and its transformation into the canonical form
Our model action, which consists of two scalar fields φ and u and an R 2 term, is described as where g is the determinant of the metric g µν , R is the scalar curvature, α( = 0) is a non-zero constant, s = ±1 is a model parameter, ∇ is the covariant derivative, and J(y) is a function of y defined as y ≡ u/φ. The action in Eq. (II.1) without the R 2 term has been proposed and studied in Refs. [41][42][43][44][45][46][47][48]. Here and the following, we have taken 2κ 2 = 1. Regarding the property of the action in Eq. (II.1), we remark that if there does not exist the R 2 term, this action is conformally invariant, while that when the R 2 term is added, the conformal invariance of this action is effectively broken. We also note that our action may be considered to be invariant yet under the restricted conformal invariance [49]. By introducing an auxiliary field Φ, the action in Eq. (II.1) can be rewritten as The simplest way to transform the action in Eq. (II.2) into the canonical form is to take the following gauge: With this gauge, the action in Eq. (II.2) is transformed into Here, we explicitly explain the properties of the actions described above. The representative feature of the action in Eq. (II.1) is that this action is not conformally invariant, while it is scale invariant. To eliminate the non-minimal coupling between the scalar fields and the scalar curvature in the action in Eq. (II.2), we have used the gauge in Eq. (II.3). As a result, for the resultant action in Eq. (II.4), the scale invariance is broken, although the canonical form, i.e., the Einstein-Hilbert term, is recovered. Therefore, it is considered that the form of the gauge in Eq. (II.2) has two aspects. One is to make the action canonical, but the other is to apparently break its scale invariance of the action.
From the action in Eq. (II.4), we obtain the gravitational equation: where R µν is the Ricci tensor, and the equations of motion for the scalar fields φ and u: where ≡ g µν ∇ µ ∇ ν is the covariant d'Alembertian operator for scalar quantities, and the prime denotes the derivative with respect to y. Furthermore, the action in Eq. (II.4) is represented as Here, V (φ, u, J) is a potential for u and φ, defined as In this case, the gravitational field equation can be expressed by The equations of motion for u and φ also become quite simple as is the scale factor. The Hubble parameter is given by H ≡ȧ/a, where the dot shows the time derivative. In this background space-time, the gravitational field equations read (II.14) In addition, the equation of motions for φ and u become (II.16)

B. Single dynamical scalar field model
We consider that φ is the inflaton field and the slow-roll inflation occurs. Namely, φ ≪ 3Hφ in the equation of motion for φ (II.15) and the kinetic energy (1/2)φ 2 of φ is much smaller than the potential energy V . We also suppose that the mass of the other scalar field u is much smaller than the Hubble parameter at the inflationary stage and hence the amplitude of u does not vary during inflation. Accordingly, we set u = u 0 (= constant) at the inflationary stage. In this case, it follows from Eq. (II.12) that This has to be satisfied for a value of φ(t) and u = u 0 , although it cannot be true for an arbitrary function J(y). In general, we can take u 0 = 0 and J = J(y 2 ). The simplest case is u 0 = 0 and J = 1. For such a case, the scalar field u is totally decomposed from equations during inflation. Consequently, the effective potential of the inflaton φ is given by 18) with C a constant. This expression is appropriate for any function J = J(y 2 ). The number of e-folds during inflation is represented as where φ f is the amplitude of φ at the end of inflation. For the effective potential in Eq. (II. 18) Hence, it is seen that there exist the second and third logarithmic correction terms in Eq. (II.20) in comparison with the number of e-folds for the quartic inflaton potential as λφ 4 with λ a constant. The slow-roll parameters in the so-called kinematic approach are defined as (for reviews, see, for instance, [17,18]) For the slow-roll inflation, the spectral index n s of the curvature perturbations and the tensor-to-scalar ratio r of the density perturbations are described by [50,51] For the slow-roll inflation driven by the effective potential in Eq. (II.18), the slow-roll parameters are written as where we have used Eqs. (II.11), (II.13), and (II.14). With Eq. (II.25), we find that s has to be negative because ǫ > 0. When ǫ < 0, it follows from Eq. (II.21) thatḢ > 0. This means that not the slow-roll inflation but the so-called super-inflation can occur. If we take the number of e-folds, whose value has to be large such as N e = 50-60 enough to solve the so-called horizon and flatness problem and the values of n s and r suggested by observations, it is apparently regarded that there are three equations (II.20), (II.23) and (II.24) for three independent variables (φ, α, C). By solving these three equations, we can estimate viable values of our model parameters. However, α and C are incorporated into all the equations in the form of x ≡ αC. Hence, it is necessary to analyze a system of two equations, for example, Eqs. (II.23) and (II.20), whereas r should depend on n s . From this reason, we may try to modify the initial action in Eq. (II.1) in the following way. Let us suppose that s is an arbitrary numerical parameter. In this case, we have a system of three equations for three variables. We also remark that our potential is well known as a potential of "Spontaneous symmetry breaking inflation (SSBI)" and a viable inflationary model can be constructed [20], although only for positive values of the parameter s.
In our model, n s ≈ 0.96 can be realized for 50 ≤ N e ≤ 60 and a wide range of parameter s, but r < 0.11 cannot be satisfied for any values of s. Therefore, we may put s = −1, although it is difficult to produce the value of r compatible with the Planck/WMAP data in our model. In Fig. 1, we show the values of n s and r in the wide range of x and φ for s = −1. From this figure, we see that we can acquire r ≈ 0.21 for n s ≈ 0.96. In this case, a typical value of x is very large. For instance, for N e = 55 and n s = 0.9603, we find φ ≈ 34.7, x ≈ 2.8 × 10 9 , and eventually obtain r ≈ 0.212.
We explore the values of our initial parameters α and C. For clarity, we take N e = 55. Using the definition of x ≡ αC, we find The sheet whose height changes from ∼ 0 to ∼ 1 denotes ns, while the other sheet whose height varies from ∼ 0 to ∼ 5 shows r.
This implies that the viable values are C < 3.0 × 10 −6 and hence α > 1.0 × 10 15 . We mention that for such a large value of α, the term 1/ (4α) in the Friedman equation becomes small, and therefore it could play a role of the effective cosmological constant Λ eff . This term may lead to the late-time cosmic acceleration.

C. Equilibrium points
Next, we investigate equilibrium points in the system. We have the system consisting of two dynamical equations (II.11) and (II.12) with a constraint equation (II.13). To examine equilibrium points in this system, we need to rewrite the system of two second order differential equations as that of four first order differential equations. It is clear that this task is equivalent to explore the shape of the potential, namely, to find its extreme values and study their natures. We execute the numerical analysis by using the graphics of the potential for several values of parameters.
For most of possible shapes of the function J(y), there is no true minimum point. This is because the value of the potential in Eq. (II.9) on the lines u 2 = φ 2 is exactly equal to 1/ (2α). This fact is true for such kind of functions as J = C, J = C(y 2 − 1) m with m a constant, J = C cos 2 (y 2 ), J = C exp(y 2 ), and so on. The typical behavior of a potential for such kind of functions is drawn in Fig. 2. Nevertheless, it is possible to construct functions J for which the potential in Eq. (II.9) has a true minimum. For instance, we have The typical behaviors of the potential V as a function of φ and u in Eq. (II.9) are plotted in Figs. 2-4.

D. Re-collapse and bounce solutions
We explore the possible re-collapse and bounce solutions 4 . For such kinds of solutions, the relation H = 0 has to be satisfied at some time. At the time, from Eqs. (II.13) and (II.14), we havė It is clear that the conditionḢ > 0 is necessary for bouncing solutions, whereas the conditionḢ < 0 has to be met for re-collapsing solutions. Accordingly, it can be seen from the general form of the potential V in Eq. (II.9) that for J > 0, the possibility for bouncing solutions to exist becomes higher because the value of V is defined to be positive.
On the other hand, if the potential can have negative values, the existence of re-collapsing solutions is possible as well. However, in the case that the value of the potential is always positive, it is impossible for re-collapsing solutions to exist for, e.g., J given by Eq. (II.29) with α > 0 and C > 0. Related to the bouncing solutions, we mention that the anti-gravity regime in the extended gravity theories with the Weyl invariance [41][42][43][44]53] including F (R) gravity [48] has been examined.

III. INFLATIONARY COSMOLOGY
In this section, we reconsider the theory whose action is described by Eq. (II.1) and build an inflationary model in another way.

A. Conformal transformation
The action in Eq. (II.1) is written in the so-called Jordan frame. Instead of taking the gauge in Eq. (II.3), we first introduce an auxiliary field Φ and then make the conformal transformation, i.e., the Weyl re-scaling, of the metric from the Jordan frame to the Einstein frame [40] g µν = Λḡ µν , (III.1) where the bar shows the quantities in the Einstein frame. The Ricci scalar is transform as By plugging this relation into Eq. (II.2) and setting where we have removed the auxiliary field Φ with the relation (III.3). If we define Λ as Λ ≡ e λ with λ a scalar field, Eq. (III.4) is represented as For simplicity, from this point, we will not write the bar over the quantities and operators in the Einstein frame. We introduce the following form of the potential The equations of motion for the scalar fields λ, φ, and u are given by Moreover, the gravitational field equation becomes In the FLRW background, from Eqs. (III.7)-(III.9), we acquire Furthermore, from Eq. (III.10), we get

B. Inflationary model
To calculate the observables for inflationary models including the spectral index of curvature perturbations n s and the tensor-to-scalar ratio r, it is necessary for two scalar fields to be made constants. If the scalar field λ is a constant and one of the scalar fields φ or u plays a role of the inflaton field, we obtain the similar theory to that described in the previous section.
We suppose that the scalar field λ is the inflaton field and other two scalar fields are set to be constants as φ = φ 0 (= constant) and u = u 0 (= constant). It should be cautioned that both these scalar fields φ and u cannot be made zero, because in this case, V eff = 1/(2α) and hence the spectrum of the curvature perturbations is flat. Nevertheless, it is possible for one of these scalar fields (e.g., u) to be taken as zero, while the other has a non-zero value. In fact, however, it is impossible to construct any inflationary model with at least (even) one viable parameter in this way. This is a reason why the relations φ = φ 0 = 0 and u = u 0 = 0 with u 0 = ±φ 0 are forced to be set. Therefore, with Eqs. (III.12) and (III.13), the conditions V φ = 0, V u = 0 during inflation are expressed as There are the following cases in which these equations are realized: • Case 1: The relation K ≡ (s/12) e λ (φ 2 0 − u 2 0 ) ≫ 1 is met. In this case, the first terms in the brackets [ ] in Eqs. (III.16) and (III.17) may be neglected. Namely, the term with e λ is much smaller than any other terms. Accordingly, all the other terms have the multiplier factor e 2λ incorporated in the same way, so that this overall factor can removed from the other terms.
• Case 2: The values of the first and second terms in the brackets [ ] may be similar with each other, but the overall coefficient term s/ (6α) is sufficiently small. As a consequence, the first terms in Eqs. (III.16) and (III.17) are suppressed to be much smaller than all the other terms in these equations.
(III.18) Therefore, the function J has to reach its extreme value when the argument becomes y = y 0 . Moreover, by substituting this value into Eqs. (III.16) and (III.17), we see that in Case 1, the relation J(y 0 ) = −s 2 /(288α) has to be met, whereas in Case 2, the relation J(y 0 ) = 0 has to be realized. Note also that a very wide class of functions may satisfy both these conditions. In the following, we study these two cases in more detail.

Case 1
In Case 1, by combining the value of the function J with the expression of the potential in Eq. (III.6), we have the effective inflaton potential This is some kind of the potential in the modified Higgs inflation model [14]. The number of e-folds in the slow-roll regime is expressed as 20) with λ f the value of λ at the end of inflation. If λ ≫ λ f , the representation of N e (λ) for the effective potential in Eq. (III.19) becomes (III.21) or by taking into account the assumption that K ≡ (s/12) e λ (φ 2 0 − u 2 0 ) ≫ 1, we may rewrite N e (λ) as With Eqs. (II.21) and (II.22), the slow-roll parameters can be described by using the effective potential as Here, the coefficients of 1/3 in Eq. (III.23) and 2/3 in Eq. (III.24) originate from the non-canonical definition of the scalar field λ, namely, the coefficient of the kinetic term for λ in the action in Eq. (III.5) is 3/2, and not 1/2. Eventually, n s and r are described as where we have defined ζ as ζ ≡ s(φ 2 0 −u 2 0 ). Consequently, we see that n s can be written as a function of r as n s = n s (r). However, the values of the parameter ζ lead to those of n s and r, which are compatible with the observations. Indeed, for N e = 60 and ζ = 0.01, we obtain r = 0.004 and n s = 0.9617. In this case, we find K = 0.013 (≪ 1). For N e = 50 and ζ = 0.18 (the value of ζ for N e = 50 becomes larger than that for N e = 60 by one order of magnitude), we acquire r = 0.005 and n s = 0.9568. In this case, we have K = 0.015 (≪ 1). As a result, we see that this case is inconsistent, because the initial assumption is K ≫ 1, but the consequences suggest K ≪ 1. It means that this case cannot be realized, in contrast with Case 2, as is shown next.

Case 2
In Case 2, as already mentioned above, we have the following relations for the function J: J(y 0 ) = 0 and J ′ (y 0 ) = 0. A simple example of such kind of the function J is where q is a constant. The effective potential is given by (III. 28) where we have used ζ ≡ s(φ 2 0 − u 2 0 ). In addition, the expression of N e reads N e (φ) = 3 4 λ + 9 ζe λ . (III.29) Accordingly, we find (III.31) Thus, it is seen that this theory also has only one parameter, and therefore we obtain n s = n s (r). Through the numerical calculations, the value of n s may be set near the value of n s = 0.96 suggested by the observations. In fact, in the case that N e = 60, for ζ = 10 −1 , 10 −2 , 10 −3 , 10 −4 , 10 −5 , and 10 −6 , we have n s = 0.9652, 0.9641, 0.9629, 0.9617, 0.9604, and 0.9589, respectively. For all of these cases, the value of r is about r = 0.004. Moreover, similarly to the above results, in the case that N e = 50, for ζ = 10 −1 , 10 −2 , 10 −3 , 10 −4 , 10 −5 , and 10 −6 , we acquire n s = 0.9577, 0.9561, 0.9543, 0.9524, 0.9504, and 0.9481, respectively. For all of these cases, the value of r is in the range of r = 0.005-0.008. Consequently, in this case, we can obtain n s ≈ 0.96 and r < 0.11. These results are compatible with the observations of the Planck 2015. Furthermore, we examine another aspect of this theory. If N e = 60 with ζ = 10 −1 , 10 −2 , 10 −3 , 10 −4 , 10 −5 , and 10 −6 , or if N e = 50 with these values of ζ, namely, the same combinations of values of N e and ζ shown above, we obtain the estimation as (ζ/12) e λ f ≈ 0.14-0.17. It means that for all the combinations of values of N e = 50-60 (during inflation) and those of ζ described above, the value of the effective potential V eff may be estimated as This relation implies that for sufficiently large values of α, the values of the potential become under the Planck scale. This fact is in good agreement with our initial assumption that s/ (6α) ≪ 1. As a result, it is considered that Case 2 is quite viable.

IV. DYNAMICAL TWO SCALAR FIELD MODEL
In this section, we explore the theory proposed in Sec. II, whose action is given by Eq. (II.1), and consider how to realize inflation by using both dynamical two scalar fields φ and u 7 .
The definitions of the slow-roll parameters in Eqs. (II.21) and (II.22) are used to describe the slow-roll inflation in the present theory. In addition, the spectral index n s of the curvature perturbations and the tensor-to-scalar ratio r are supposed to be represented as Eqs. (II.24) and (II.23), respectively. This assumption has been justified in Ref. [55].
We consider the FLRW background. We start with the system of the equations (II.11) and (II.12), namely, (II.13) and (II.14) in the FLRW space-time. By imposing the standard slow-rolling conditionsü ≪uH,φ ≪φH, anḋ φ 2 −u 2 ≪ H 2 , the Friedmann equation and the equations of motion for φ and u are reduced to be the following simple forms: First, we derive the number of e-folds by using the initial definition where a i and a f are the values of the scale factor a(t) at the beginning t i and end t f of inflation, respectively. Moreover, we also have the following relation By using Eqs. (IV.1)-(IV.3), we get with φ f and u f the amplitudes of φ and u at the end of inflation, respectively. It follows from the definition in Eq. (II.21) that ǫ is described as On the other hand, with the definition in Eq. (II.22) and Eq. (IV.1), we find where in deriving the last equality, we have used Eqs. (IV.2) and (IV.3). The expressions (IV.6)-(IV.8) can be used for any choice of the function J, but in the most general case, the calculations would be too cumbersome. In addition, the expression (IV.6) leads to only a relation between the scalar fields φ and u, and not both the values of φ and u simultaneously. Hence, another additional assumption to determine both the values of φ and u (e.g., φ ≈ u) is necessary, or one of the values of φ and u may be regarded as a free parameter.
In the simplest case J = C, it is possible to acquire the values of the main parameters analytically 8 . In this case, from Eq. (IV.6), we find (IV.9) Expressions for ǫ and η takes the form With the expression of N e in (IV.9), the value of (φ 2 − u 2 ) can be obtained. It determines the value of ǫ, but that of η cannot be found by using only that of (φ 2 − u 2 ). Since ǫ(> 0) should be positive, from Eqs. (IV.9) and (IV.10), we have φ 2 − u 2 > 0 and s < 0. These relations and x = αC > 0 lead to the positive value of η. We estimate numerical values. For simplicity, we take s = −1. From Eqs. (IV.9) and (IV.10), we get (φ 2 − u 2 ) and x for any interesting values of N e and r = 16ǫ. Thanks to the structure of the equations, we obtain the value of r as r < 0.11. On the other hand, it is seen that n s = n s (r, G 4 ), where we have defined a new parameter G 4 ≡ u 2 φ 2 . This relation yields the value of n s consistent with the observations. Indeed, for N e = 50 and x = 10 30 , we have φ 2 − u 2 = 2593 and r ≃ 0.1. In this case, n s is represented as n s = 3.7 × 10 −9 G 4 + 0.9815 . (IV.12) Moreover, for N e = 60 and x = 10 25 , we acquire φ 2 − u 2 = 2317 and r ≃ 0.11. In this case, n s is expressed as If the larger x is, the smaller r becomes, but the minimum value of n s in both Eqs. (IV.12) and (IV.13) increases. Equations (IV.12) and (IV.13) imply that it is impossible to for the pair of n s and r to take their values compatible with the observations. The consequence is the same for the different values of s and N e . Thus, for the simplest case of J = C, it is impossible to realize the Planck analysis with two dynamical scalar fields. However, even in this simplest case, the BICEP2 results may be realized. For N e = 50 and x = 10 6 , we obtain φ 2 − u 2 = 1267, r ≃ 0.202, and n s = 3 × 10 −8 G 4 + 0.9621. Therefore, the case of G 4 < 10 5 is quite viable and consist with the BICEP2 experiment. As a result, in this simplest case (of function J), the Planck results cannot be realized. Nevertheless, the model with two dynamical scalar fields leads to significantly different results in comparison with the inflationary model with the single scalar field. There is another possibility to consider the more viable types of the function J, e.g., Eq. (II.29). In such a case, it is quite difficult to analyze the equations analytically (for example, N e is described as an integral equation).

V. GRACEFUL EXIT FROM INFLATION
In this section, we investigate the graceful exit from inflation, namely, the instability of the de Sitter solution at the inflationary stage for the present theory. Especially, we demonstrate the instability of the de Sitter solution in the case of one dynamical scalar field in the Einstein frame. In this case, as shown in Sec. III B 2, the spectral index and the tensor-to-scalar ratio can be compatible with the observations by the Planck satellite. We also examine the contribution of an R 2 term in the Jordan frame to the instability of the de Sitter solution during inflation.
We set the perturbations of the Hubble parameter at the inflationary stage as where H inf (> 0) (= constant) is the Hubble parameter during inflation (whose value is positive), and δ(t) means the perturbations from the de Sitter solution. In the present case, only the scalar field λ is dynamical, and φ and u are constants. By taking the time derivative of the gravitational field equation (III.15) with V = V eff , where V eff is given by Eq. (III.28), we have Moreover, from Eq. (III.11) with V = V eff (in Eq. (III.28)), the equation of motion for λ reads Since we consider the solution at the inflationary stage in the early universe, we take the limit t → 0. In this limit, we could have an approximate solution λ ≈ ln (Ht) for Eq. (V.3) with the quasi de Sitter solution H ≈ 1/ (3t). By taking into account this solution, we see that in the limit t → 0, the third term proportional toλ in the left-hand side of Eq. (V.2) is much larger than the fourth term, which is proportional only toλ, and therefore the fourth term could be neglected.
To examine the instability of the de Sitter solution, we represent δ(t) as with β a constant. When β is positive, the de Sitter solution during inflation becomes unstable, and hence the universe can exit from the inflationary stage and then enter the reheating stage. This is because for β > 0, the amplitude of δ(t) increases in time. By substituting the expression of the perturbations in Eq. (V.1) with Eq. (V.4) into Eq. (V.2) and using approximate solutions λ ≈ ln (Ht) and H ≈ 1/ (3t), we find The solutions for this equation are obtained as As a consequence, we acquire the positive solution of β = β + > 0, and thus the universe can exit from the inflationary stage. We note that even if the other scalar field becomes dynamical, the procedure to examine the instability of the de Sitter solution is basically the same as the one demonstrated above. We examine the perturbations of the Hubble parameter by using the gravitational field equation with solutions for the equation of motions in terms of two dynamical scalar fields. Qualitatively, the form of the solution for the perturbations will be changed, but in principle there can exist a solution to represent the property that the de Sitter solution is unstable.
The contribution of the R 2 term in the Jordan frame to the instability of the de Sitter solution is included in the scalar field λ and its dynamics through the auxiliary field Φ in the Einstein frame. This fact can be seen from the action in Eq. (II.2), and the relation (III.3) with Λ = e λ . Accordingly, it is considered that the R 2 term can be related to the graceful exit from inflation, i.e., the instability of the de Sitter solution to describe the slow-roll inflation.

VI. CONCLUSIONS
In the present paper, we have studied inflationary cosmology in a theory where there exist two scalar fields nonminimally coupling to the Ricci scalar and an additional R 2 term. We have investigated the slow-roll inflation in the case of one dynamical scalar field and that of two dynamical scalar fields. We have analyzed the spectral index n s of scalar mode of the density perturbations and the tensor-to-scalar ratio r in comparison with the observations of the recent Planck and BICEP2 results.
For the case of single dynamical scalar field in the Jordan frame, if the number of e-folds during inflation is 50 ≤ N e ≤ 60, we have n s ≈ 0.96 and r = O(0.1). On the other hand, in the Einstein frame, for the case of one dynamical scalar field, we can acquire n s ≈ 0.96 and r < 0.11. These are consistent with the Planck 2015 results. Furthermore, for the case of two dynamical scalar fields in the Jordan frame, when 50 ≤ N e ≤ 60, we obtain n s ≈ 0.96 and r = O(0.1). As a result, we have found that in the present theory, the spectral index and the tensor-to-scalar ratio can be compatible with the recent Planck analysis. We have also shown that in the present theory, the de Sitter solution representing the inflationary stage is unstable, and therefore the universe can successfully exit from inflation. The R 2 term in the Jordan frame is considered to be related to the instability of the de Sitter solution, namely, the graceful exit from the inflationary stage.
As the further developments on the cosmological consistent scenario in the present theory, it is possible to unify inflation in the early universe realized by the R 2 term and the late-time cosmic acceleration by the dark energy component of one of the scalar fields. In this theory, dark matter can also be explained by the other scalar field. To construct such a unified scenario is the significant purpose in this work. In addition, the consequences in multiple field inflation models seems to be more suitable for the observational data obtained form the Planck satellite than single field inflation models.
It is remarked that the action in Eq. (II.1) may be generalized so that the conformal invariance can be broken by an arbitrary function of the F (R) term. Thanks to the additional term in the gravity sector, the novel contributions to cosmology yield from the breaking of conformal invariance of the theory. Thus, we have more possibilities to realize the unified scenario of inflation, dark energy, and dark matter mentioned above.
In fact, we may further extend the investigations in this work. We start from the theory with the conformally invariance and three scalar fields, and then add an R 2 term. We explicitly show the expressions of these actions below. In Ref. [56], a theory without an R 2 term and that with the general potential (φ 2 − u 2 − ψ 2 ) 2 J(u/φ, u/ψ) have been explored. The action for the theory with the more specific potential ψ 4 J(u/φ, u/ψ) and without the R 2 term) is given by where ǫ is a constant and ψ is the additional third scalar field to two scalar fields φ and u already introduced. If the R 2 term is added to Eq. (VI.1), the representation of the action reads Moreover, when the general form of the potential is included and the signature of ǫ is opposite to that in the action in Eq. (VI.2), the action is represented as The same analysis as in Secs. II-IV may be adopted to the models in Eqs. (VI.1)-(VI.3) in which the additional third scalar field ψ may play a role of the other specie of dark matter, or the additional contribution to the dark energy dominated stage or inflation.