1 Introduction

The natures on inflation [15] in the early universe have been revealed by the recent cosmological observations such as the Wilkinson Microwave Anisotropy Probe (WMAP) [6, 7], the Planck satellite [8, 9], and the BICEP2 experiment [10, 11] on the quite tiny anisotropy of the cosmic microwave background (CMB) radiation. Owing to the release of the recent observational data, in addition to seminal inflation with a single scalar field such as new inflation [3, 4], chaotic inflation [12], natural inflation [13], and power-law inflation with the exponential inflaton potential [14], novel models of single-field inflation have been proposed in Refs. [1590]Footnote 1 (for reviews on more various inflationary models, see, e.g., [124128]).

In addition to inflationary models driven by a scalar field (i.e., the inflaton field) described above, there have been considered the so-called Starobinsky inflation [5, 129] originating from the higher-order curvature term such as an \(R^2\) term,Footnote 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 theory including F(R) gravity to account for the late-time cosmic acceleration (for reviews on the dark energy problem and modified gravity theories, see, for instance, [130140]). Various inflationary models in modified gravity theories corresponding to extensions of Starobinsky inflation have been explored in Refs. [141153].

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.Footnote 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 the two scalar fields constant) and (ii) two dynamical scalar fields. As a consequence, we analyze the spectral index of the 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 the 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 [155159] including the so-called curvaton scenario [160164] have been constructed and the cosmological perturbations in these models have also been investigated [165177].

We use units of \(k_\mathrm {B} = c = \hbar = 1\) and express the gravitational constant \(8 \pi G_\mathrm {N}\) by \({\kappa }^2 \equiv 8\pi /{M_{\mathrm {Pl}}}^2\) with the Planck mass of \(M_{\mathrm {Pl}} = G_\mathrm {N}^{-1/2} = 1.2 \times 10^{19}\) GeV.

The organization of the paper is as follows. In Sect. 2, 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 Sect. 3, with the conformal transformation [178, 179], 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 the other scalar fields are set constant. In Sect. 4, 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 Sect. 5, 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 Sect. 3, 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 Sect. 6.

2 Two-scalar-field model with breaking the conformal invariance

2.1 Model action and its transformation into the canonical form

Our model action, which consists of two scalar fields \(\phi \) and u and an \(R^2\) term, is described by

$$\begin{aligned} S= & {} \int d^4 x \sqrt{-g}\left\{ \frac{\alpha }{2} R^2 + \frac{s}{2}\left[ \frac{(\phi ^2-u^2)}{6}R + (\nabla \phi )^2 - (\nabla u)^2 \right] \right. \nonumber \\&-\left. (\phi ^2-u^2)^2J(y) \right\} , \end{aligned}$$
(2.1)

where g is the determinant of the metric \(g_{\mu \nu }\), R is the scalar curvature, \(\alpha (\ne 0)\) is a non-zero constant, \(s=\pm 1\) is a model parameter, \(\nabla \) is the covariant derivative, and J (y) is a function of y defined as \(y \equiv u/\phi \).

The action in Eq. (2.1) without the \(R^2\) term has been proposed and studied in Refs. [180189].

Here and in the following, we take \(2\kappa ^2 = 1\). Regarding the property of the action in Eq. (2.1), we remark that if there does not exist an \(R^2\) term, this action is conformally invariant, while 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 [190].

By introducing an auxiliary field \(\Phi \), the action in Eq. (2.1) can be rewritten as

$$\begin{aligned} S&= \int d^4 x \sqrt{-g}\left\{ \left[ \Phi + \frac{s}{12}(\phi ^2-u^2) \right] R -\frac{\Phi ^2}{2\alpha }\right. \nonumber \\&\quad + \left. \frac{s}{2} \left[ (\nabla \phi )^2 - (\nabla u)^2 \right] - (\phi ^2-u^2)^2J(y) \right\} . \end{aligned}$$
(2.2)

The simplest way to transform the action in Eq. (2.2) into the canonical form is to take the following gauge:

$$\begin{aligned} \Phi + \frac{s}{12}(\phi ^2-u^2) = 1. \end{aligned}$$
(2.3)

With this gauge, the action in Eq. (2.2) is transformed into

$$\begin{aligned} S&= \int d^4 x \sqrt{-g}\left\{ R -\frac{1}{2\alpha }\left[ 1 - \frac{s}{12}(\phi ^2-u^2) \right] ^2\right. \nonumber \\&\quad \left. + \frac{s}{2}\left[ (\nabla \phi )^2 - (\nabla u)^2 \right] - (\phi ^2-u^2)^2J(y) \right\} . \end{aligned}$$
(2.4)

Here, we explicitly explain the properties of the actions described above. The representative feature of the action in Eq. (2.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. (2.2), we have used the gauge in Eq. (2.3). As a result, for the resultant action in Eq. (2.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. (2.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. (2.4), we obtain the gravitational equation:

$$\begin{aligned}&R_{\mu \nu }-\frac{1}{2}Rg_{\mu \nu }+\frac{s}{2}(\nabla _{\mu }\phi \nabla _{\nu }\phi -\nabla _{\mu } u \nabla _{\nu } u) \nonumber \\&\quad + \frac{1}{2}g_{\mu \nu } \left\{ \frac{1}{2\alpha }\left[ 1 - \frac{s}{12}(\phi ^2-u^2) \right] ^2\right. \nonumber \\&\quad \left. - \frac{s}{2}\left[ (\nabla \phi )^2-(\nabla u)^2\right] +(\phi ^2-u^2)^2 J(y) \right\} =0 , \end{aligned}$$
(2.5)

where \(R_{\mu \nu }\) is the Ricci tensor, and the equations of motion for the scalar fields \(\phi \) and u:

$$\begin{aligned}&s\Box \phi - \frac{s}{6\alpha }\phi \left[ 1 - \frac{s}{12}(\phi ^2-u^2) \right] \ \nonumber \\&\quad +\, 4\phi (\phi ^2-u^2)J(y) - \frac{u}{\phi ^2}(\phi ^2-u^2)^2 J'(y)=0 ,\end{aligned}$$
(2.6)
$$\begin{aligned}&s\Box u - \frac{s}{6\alpha }u\left[ 1 - \frac{s}{12}(\phi ^2-u^2) \right] \nonumber \\&\quad +\, 4u(\phi ^2-u^2)J(y) - \frac{1}{\phi }(\phi ^2-u^2)^2 J'(y)=0 . \end{aligned}$$
(2.7)

Here, \(\Box \equiv g^{\mu \nu } {\nabla }_{\mu } {\nabla }_{\nu }\) is the covariant d’Alembertian operator for scalar quantities, and the prime denotes the derivative with respect to y. Furthermore, the action in Eq. (2.4) is represented as

$$\begin{aligned} S\!=\! \int d^4 x \sqrt{-g}\left\{ R \!+ \frac{s}{2} \left[ (\nabla \phi )^2 \!-\! (\nabla u)^2 \right] \!-\! V(\phi , u, J) \right\} ,\nonumber \\ \end{aligned}$$
(2.8)

where, \(V(\phi , u, J)\) is a potential for u and \(\phi \), defined as

$$\begin{aligned} V(\phi , u, J) \!\equiv \! \frac{1}{2\alpha }\left[ 1 - \frac{s}{12}(\phi ^2-u^2) \right] ^2+(\phi ^2-u^2)^2 J(y).\nonumber \\ \end{aligned}$$
(2.9)

In this case, the gravitational field equation can be expressed by

$$\begin{aligned}&R_{\mu \nu }-\frac{1}{2}Rg_{\mu \nu }+\frac{s}{2}\left( \nabla _{\mu }\phi \nabla _{\nu }\phi -\nabla _{\mu } u \nabla _{\nu } u\right) \nonumber \\&\quad + \frac{1}{2}g_{\mu \nu } \left\{ V - \frac{s}{2}\left[ (\nabla \phi )^2-(\nabla u)^2\right] \right\} =0 . \end{aligned}$$
(2.10)

The equations of motion for u and \(\phi \) also become quite simple, namely

$$\begin{aligned} s\Box \phi + V_{\phi }= & {} 0,\end{aligned}$$
(2.11)
$$\begin{aligned} s\Box u - V_{u}= & {} 0, \end{aligned}$$
(2.12)

with \(V_{\phi } \equiv \partial V/\partial \phi \) and \(V_{u} \equiv \partial V/\partial u\).

The flat Friedmann–Lemaître–Robertson–Walker (FLRW) metric \(\hbox {d}s^2 = -\hbox {d}t^2 + a^2(t) \sum _{i=1,2,3}\left( dx^i\right) ^2\) is assumed, where a(t) is the scale factor. The Hubble parameter is given by \(H \equiv \dot{a}/a\), where the dot shows the time derivative. In this background space-time, the gravitational field equations read

$$\begin{aligned}&3H^2 + \frac{s}{4}(\dot{\phi }^2-\dot{u}^2)-\frac{1}{2}V=0 ,\end{aligned}$$
(2.13)
$$\begin{aligned}&2\dot{H} + 3H^2 - \frac{s}{4}(\dot{\phi }^2-\dot{u}^2)-\frac{1}{2}V=0. \end{aligned}$$
(2.14)

In addition, the equations of motion for \(\phi \) and u become

$$\begin{aligned} s(\ddot{\phi }+3H\dot{\phi }) - V_{\phi }= & {} 0, \end{aligned}$$
(2.15)
$$\begin{aligned} s(\ddot{u}+3H\dot{u}) + V_{u}= & {} 0. \end{aligned}$$
(2.16)

2.2 Single dynamical scalar field model

We consider that \(\phi \) is the inflaton field and the slow-roll inflation occurs. Namely, \(\left| \ddot{\phi }\right| \ll \left| 3H \dot{\phi } \right| \) in the equation of motion for \(\phi \) (2.15) and the kinetic energy \(\left( 1/2\right) \dot{\phi }^2\) of \(\phi \) 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 ({=}\mathrm {constant})\) at the inflationary stage. In this case, it follows from Eq. (2.12) that

$$\begin{aligned} V_{u}&=\frac{s}{6\alpha }u\left[ 1 - \frac{s}{12}(\phi ^2-u^2) \right] - 4u(\phi ^2-u^2)J(y)\nonumber \\&\quad + \frac{1}{\phi }(\phi ^2-u^2)^2J'(y) =0. \end{aligned}$$
(2.17)

This has to be satisfied for a value of \(\phi (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 separated from the equations during inflation. Consequently, the effective potential of the inflaton \(\phi \) is given by

$$\begin{aligned} V_\mathrm {eff}(\phi ) = \frac{1}{2\alpha } \left( 1 - \frac{s}{12}\phi ^2 \right) ^2+C\phi ^4 , \end{aligned}$$
(2.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

$$\begin{aligned} N_e(\phi )= & {} \int ^{\phi _\mathrm {f}}_{\phi }H(\hat{\phi })\frac{d\hat{\phi }}{\dot{\hat{\phi }}}=-\frac{s}{2}\int ^{\phi }_{\phi _\mathrm {f}}\frac{V(u, \hat{\phi }, J)}{V_\phi (u, \hat{\phi }, J)}d\hat{\phi } , \end{aligned}$$
(2.19)

where \(\phi _\mathrm {f}\) is the amplitude of \(\phi \) at the end of inflation. For the effective potential in Eq. (2.18) with \(\phi \gg \phi _\mathrm {f}\), we have

$$\begin{aligned}&N_e(\phi )\nonumber \\&\quad =-\frac{s}{2}\left[ \frac{\phi ^2}{8} + \frac{432\alpha C \ln (s^2\phi ^2+288\alpha C\phi ^2 -12s)}{s\left( s^2+288\alpha C\right) } -\frac{3\ln \phi }{s} \right] .\nonumber \\ \end{aligned}$$
(2.20)

Hence, it is seen that there exist second and third logarithmic correction terms in Eq. (2.20) in comparison with the number of e-folds for the quartic inflaton potential as \(\lambda \phi ^4\) with \(\lambda \) a constant.

The slow-roll parameters in the so-called kinematic approach are defined as (for reviews, see, for instance, [124, 125])

$$\begin{aligned} \epsilon\equiv & {} -\frac{\dot{H}}{H^2} , \end{aligned}$$
(2.21)
$$\begin{aligned} \eta\equiv & {} \epsilon - \frac{\ddot{H}}{2H\dot{H}} . \end{aligned}$$
(2.22)

For the slow-roll inflation, the spectral index \(n_\mathrm {s}\) of the curvature perturbations and the tensor-to-scalar ratio r of the density perturbations are described by [191, 192]

$$\begin{aligned} n_\mathrm {s}= & {} 1 -6\epsilon + 2\eta , \end{aligned}$$
(2.23)
$$\begin{aligned} r= & {} 16 \epsilon . \end{aligned}$$
(2.24)

In the case of the slow-roll inflation driven by the effective potential in Eq. (2.18), the slow-roll parameters are written as

$$\begin{aligned} \epsilon= & {} -\frac{1}{s}\left( \frac{1}{V_\mathrm {eff}(\phi )} \frac{dV_\mathrm {eff}(\phi )}{d\phi } \right) ^2\nonumber \\= & {} -\frac{16\phi ^2}{s}\left[ \frac{(s^2+288\alpha C)\phi ^2-12s}{288\alpha C\phi ^4+(12-s\phi ^2)^2} \right] ^2 ,\end{aligned}$$
(2.25)
$$\begin{aligned} \eta= & {} -\frac{2}{s} \frac{1}{V_\mathrm {eff}(\phi )} \frac{d^2 V_\mathrm {eff}(\phi )}{d \phi ^2}\nonumber \\= & {} -\frac{24}{s} \frac{(s^2+288\alpha C)\phi ^2-4s}{288\alpha C\phi ^4+(12-s\phi ^2)^2} , \end{aligned}$$
(2.26)

where we have used Eqs. (2.11), (2.13), and (2.14). With Eq. (2.25), we find that s has to be negative because \(\epsilon >0\). When \(\epsilon <0\), it follows from Eq. (2.21) that \(\dot{H} >0\). This means that not 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 enough, such as \(N_e=50\)–60, to solve the so-called horizon and flatness problems and the values of \(n_\mathrm {s}\) and r suggested by the observations, it is apparently regarded that there are three equations, Eqs. (2.20), (2.23), and (2.24), for the three independent variables (\(\phi \), \(\alpha \), C). By solving these three equations, we can estimate viable values of our model parameters. However, \(\alpha \) and C are incorporated into all the equations in the form of \(x \equiv \alpha C\). Hence, it is necessary to analyze a system of two equations, for example, Eqs. (2.23) and (2.20), whereas r should depend on \(n_\mathrm {s}\). For this reason, we may try to modify the initial action in Eq. (2.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 [127], although only for positive values of the parameter s.

According to the Planck 2015 results, \(n_{\mathrm {s}} = 0.968 \pm 0.006\, (68\,\%\,\mathrm {CL})\) [8, 9] and \(r < 0.11\, (95\,\%\,\mathrm {CL})\) [9]. These values are consistent with those obtained by the WMAP satellite [6, 7]. The BICEP2 experiment has suggested \(r=0.20^{+0.07}_{-0.05}\, (68\,\%\,\mathrm {CL})\) [10], but recently the B-mode polarization of the CMB radiation is attributed to an effect of the dust, and not of primordial gravitational waves [11].

In our model, \(n_\mathrm {s} \approx 0.96\) can be realized for \(50 \le N_e \le 60\) and a wide range of the 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_\mathrm {s}\) and r in a wide range of x and \(\phi \) for \(s=-1\). From this figure, we see that we can obtain \(r \approx 0.21\) for \(n_\mathrm {s} \approx 0.96\). In this case, the typical value of x is very large. For instance, for \(N_e = 55\) and \(n_\mathrm {s} = 0.9603\), we find \(\phi \approx 34.7\), \(x \approx 2.8 \times 10^9\), and eventually we obtain \(r \approx 0.212\).

Fig. 1
figure 1

\(n_\mathrm {s}\) and r as functions of x and \(\phi \) for \(s=-1\). The sheet whose height changes from \({\sim } 0\) to \({\sim } 1\) denotes \(n_\mathrm {s}\), while the other sheet whose height varies from \({\sim } 0\) to \({\sim } 5\) shows r

We explore the values of our initial parameters \(\alpha \) and C. For clarity, we take \(N_e = 55\). Using the definition of \(x \equiv \alpha C\), we find

$$\begin{aligned} V_\mathrm {eff} \sim \frac{1.0 \times 10^{15}}{\alpha } . \end{aligned}$$
(2.27)

This implies that the viable values are \(C < 3.0 \times 10^{-6}\) and hence \(\alpha > 1.0 \times 10^{15}\). We mention that for such a large value of \(\alpha \), the term \(1/\left( 4\alpha \right) \) in the Friedman equation becomes small, and therefore it could play the role of the effective cosmological constant \(\Lambda _\mathrm {eff}\). This term may lead to the late-time cosmic acceleration.

2.3 Equilibrium points

Next, we investigate equilibrium points in the system. We have a system consisting of two dynamical equations, Eqs. (2.11) and (2.12), with the constraint equation (2.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 an exploration of 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 the parameters.

For most of the possible shapes of the function J(y), there is no true minimum point. This is because the value of the potential in Eq. (2.9) on the lines \(u^2=\phi ^2\) is exactly equal to \(1/\left( 2\alpha \right) \). This fact is true for such kinds 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. (2.9) has a true minimum. For instance, we have

$$\begin{aligned} J(y) = \frac{C}{\left( y^2-1\right) ^2}. \end{aligned}$$
(2.28)

The typical behaviors of the potential V as a function of \(\phi \) and u in Eq. (2.9) are plotted in Figs. 2, 3, and 4.

Fig. 2
figure 2

V as a function of \(\phi \) and u in Eq. (2.9) near an extreme value for \(\alpha =1.0 \times 10^{16}\), \(C=3.0 \times 10^{-7}\), \(s=-1\), and \(J=C\)

Fig. 3
figure 3

V as a function of \(\phi \) and u in Eq. (2.9) around the minimum for \(\alpha =1.0 \times 10^{16}\), \(C=3.0 \times 10^{-7}\), \(s=-1\), and \(J=C(y^2-1)^{-2}\)

Fig. 4
figure 4

V as a function of \(\phi \) and u in Eq. (2.9) around \(\alpha =1.0 \times 10^{16}\), \(C=3.0 \times 10^{-7}\), \(s=-1\), and \(J=C(y^2-1)^{-2}\)

2.4 Re-collapse and bounce solutions

We explore the possible re-collapse and bounce solutions.Footnote 4 For such kinds of solutions, the relation \(H=0\) has to be satisfied at some time. At the time, from Eqs. (2.13) and (2.14), we have

$$\begin{aligned} \dot{H} =\frac{1}{2}V . \end{aligned}$$
(2.29)

It is clear that the condition \(\dot{H}>0\) is necessary for bouncing solutions, whereas the condition \(\dot{H}<0\) has to be met for re-collapsing solutions. Accordingly, it can be seen from the general form of the potential V in Eq. (2.9) that for \(J>0\), the possibility of bouncing solutions 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. (2.29) with \(\alpha >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 [180185, 198] including F(R) gravity [189] has been examined.

3 Inflationary cosmology

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

3.1 Conformal transformation

The action in Eq. (2.1) is written in the so-called Jordan frame. Instead of taking the gauge in Eq. (2.3), we first introduce an auxiliary field \(\Phi \) and then make the conformal transformation, i.e., the Weyl re-scaling, of the metric from the Jordan frame to the Einstein frame [178, 179]

$$\begin{aligned} g_{\mu \nu }=\Lambda \bar{g}_{\mu \nu }, \end{aligned}$$
(3.1)

where the bar shows the quantities in the Einstein frame. The Ricci scalar is transformed thus:

$$\begin{aligned} R=\Lambda ^{-1}\left[ \bar{R} -3\Lambda ^{-1} \bar{\Box }\Lambda +\frac{3}{2}\Lambda ^{-2}(\bar{\nabla }\Lambda )^2 \right] . \end{aligned}$$
(3.2)

By plugging this relation into Eq. (2.2) and setting

$$\begin{aligned} \Lambda \left[ \Phi +\frac{s}{12}\left( \phi ^2-u^2\right) \right] =1, \end{aligned}$$
(3.3)

we findFootnote 5

$$\begin{aligned} S= & {} \int d^4x\sqrt{-\bar{g}}\left\{ \bar{R} -\frac{3}{2}\Lambda ^{-2}(\bar{\nabla }\Lambda )^2 -\frac{1}{2\alpha }\right. \nonumber \\&+ \frac{s}{12\alpha }\Lambda (\phi ^2-u^2)-\frac{s^2}{288\alpha }\Lambda ^2(\phi ^2-u^2)^2 \nonumber \\&\left. {}+ \frac{s}{2}\Lambda \left[ (\bar{\nabla }\phi )^2 - (\bar{\nabla } u)^2 \right] - \Lambda ^2(\phi ^2-u^2)^2J(y) \right\} , \nonumber \\ \end{aligned}$$
(3.4)

where we have removed the auxiliary field \(\Phi \) by Eq. (3.3). If we define \(\Lambda \) as \(\Lambda \equiv \mathrm {e}^{\lambda }\) with \(\lambda \) a scalar field, Eq. (3.4) is represented as

$$\begin{aligned} S= & {} \int d^4x\sqrt{-\bar{g}}\left\{ \bar{R} \!-\!\frac{3}{2}(\bar{\nabla }\lambda )^2 \!-\!\frac{1}{2\alpha }\left[ 1\!-\!\frac{s}{12}\mathrm {e}^{\lambda }(\phi ^2-u^2)\right] ^2 \right. \nonumber \\&\left. + \frac{s}{2}\mathrm {e}^{\lambda }\left\{ (\bar{\nabla }\phi )^2 - (\bar{\nabla } u)^2 \right\} - \mathrm {e}^{2\lambda }(\phi ^2-u^2)^2J(y) \right\} . \nonumber \\ \end{aligned}$$
(3.5)

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:

$$\begin{aligned} V(\lambda ,\phi ,u,J)&=\frac{1}{2\alpha }\left[ 1-\frac{s}{12}\mathrm {e}^{\lambda }(\phi ^2-u^2)\right] ^2\nonumber \\&\quad + \mathrm {e}^{2\lambda }(\phi ^2-u^2)^2J(y). \end{aligned}$$
(3.6)

The equations of motion for the scalar fields \(\lambda \), \(\phi \), and u are given by

$$\begin{aligned}&3\Box \lambda +\frac{s}{2}\mathrm {e}^{\lambda }[(\nabla \phi )^2-(\nabla u)^2] - V_{\lambda }=0,\end{aligned}$$
(3.7)
$$\begin{aligned}&s\mathrm {e}^{\lambda }\Box \phi + s\mathrm {e}^{\lambda }\nabla ^{\mu }\lambda \nabla _{\mu }\phi + V_{\phi }=0,\end{aligned}$$
(3.8)
$$\begin{aligned}&s \mathrm {e}^{\lambda }\Box u + s\mathrm {e}^{\lambda }\nabla ^{\mu }\lambda \nabla _{\mu }u - V_{u}=0. \end{aligned}$$
(3.9)

Moreover, the gravitational field equation becomes

$$\begin{aligned}&R_{\mu \nu }-\frac{1}{2}Rg_{\mu \nu }+\frac{s}{2}\mathrm {e}^{\lambda }(\nabla _{\mu }\phi \nabla _{\nu }\phi -\nabla _{\mu } u \nabla _{\nu } u) -\frac{3}{2}\nabla _{\mu }\lambda \nabla _{\nu }\lambda \nonumber \\&\quad +\frac{1}{2}g_{\mu \nu } \left\{ V - \frac{s}{2}\mathrm {e}^{\lambda }\left[ (\nabla \phi )^2-(\nabla u)^2\right] +\frac{3}{2}(\nabla \lambda )^2 \right\} =0. \end{aligned}$$
(3.10)

In the FLRW background, from Eqs. (3.7)–(3.9), we obtain

$$\begin{aligned}&3\ddot{\lambda }+ 9H\dot{\lambda }+\frac{s}{2}\mathrm {e}^{\lambda }(\dot{\phi }^2-\dot{u}^2) + V_{\lambda }=0,\end{aligned}$$
(3.11)
$$\begin{aligned}&s\mathrm {e}^{\lambda }\ddot{\phi }+ 3s\mathrm {e}^{\lambda }H\dot{\phi }+ s\mathrm {e}^{\lambda }\dot{\lambda }\dot{\phi }- V_{\phi }=0,\end{aligned}$$
(3.12)
$$\begin{aligned}&s\mathrm {e}^{\lambda }\ddot{u} + 3s\mathrm {e}^{\lambda }H\dot{u} + s\mathrm {e}^{\lambda }\dot{\lambda }\dot{u} + V_{u}=0. \end{aligned}$$
(3.13)

Furthermore, from Eq. (3.10), we get

$$\begin{aligned}&3H^2 + \frac{s}{4}\mathrm {e}^{\lambda }(\dot{\phi }^2-\dot{u}^2) -\frac{3}{4}\dot{\lambda }^2-\frac{1}{2}V=0,\end{aligned}$$
(3.14)
$$\begin{aligned}&2\dot{H} + 3H^2 - \frac{s}{4}\mathrm {e}^{\lambda }(\dot{\phi }^2-\dot{u}^2) +\frac{3}{4}\dot{\lambda }^2-\frac{1}{2}V=0. \end{aligned}$$
(3.15)

3.2 Inflationary model

To calculate the observables for inflationary models including the spectral index of curvature perturbations \(n_\mathrm {s}\) and the tensor-to-scalar ratio r, it is necessary for two scalar fields to be made constant. If the scalar field \(\lambda \) is a constant and one of the scalar fields \(\phi \) or u plays the role of the inflaton field, we obtain a similar theory to that described in the previous section.

We suppose that the scalar field \(\lambda \) is the inflaton field and the other two scalar fields are set constant as \(\phi =\phi _0 ({=} \mathrm {constant})\) and \(u=u_0 ({=} \mathrm {constant})\). It should be cautioned that these two scalar fields, \(\phi \) and u, cannot be made zero, because in this case, \(V_\mathrm {eff}=1/(2\alpha )\), 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 the reason why the relations \(\phi =\phi _0 \ne 0\) and \(u=u_0 \ne 0\) with \(u_0 \ne \pm \phi _0\) are enforced. Therefore, with Eqs. (3.12) and (3.13), the conditions \(V_{\phi }=0\), \(V_u=0\) during inflation are expressed as

$$\begin{aligned} V_{\phi }= & {} -\frac{s}{6\alpha }\left[ 1-\frac{s}{12}\mathrm {e}^{\lambda }(\phi _0^2-u_0^2)\right] \mathrm {e}^{\lambda }\phi _0\nonumber \\&+ 4\mathrm {e}^{2\lambda }(\phi _0^2-u_0^2)\phi _0J(y_0)\nonumber \\&-\mathrm {e}^{2\lambda }(\phi _0^2-u_0^2)^2J'(y_0)\frac{u_0}{\phi _0^2}=0, \end{aligned}$$
(3.16)
$$\begin{aligned} V_{u}= & {} \frac{s}{6\alpha }\left[ 1-\frac{s}{12}\mathrm {e}^{\lambda }(\phi _0^2-u_0^2)\right] \mathrm {e}^{\lambda }u_0\nonumber \\&- 4\mathrm {e}^{2\lambda }(\phi _0^2-u_0^2)u_0J(y_0)\nonumber \\&+\mathrm {e}^{2\lambda }(\phi _0^2-u_0^2)^2J'(y_0)\frac{1}{\phi _0}=0. \end{aligned}$$
(3.17)

One has the following cases in which these equations are realized.

  • Case 1: The relation \(K\equiv \left( s/12\right) \mathrm {e}^{\lambda }(\phi _0^2-u_0^2)\gg 1\) is met. In this case, the first terms in the brackets \([\,\,]\) in Eqs. (3.16) and (3.17) may be neglected. Namely, the term with \(\mathrm {e}^{\lambda }\) is much smaller than any other terms. Accordingly, all the other terms have the multiplication factor \(\mathrm {e}^{2\lambda }\) incorporated in the same way, so that this overall factor can be removed from the other terms.

  • Case 2: The values of the first and second terms in the brackets \([\,\,]\) may be similar to each other, but the overall coefficient term \(s/\left( 6\alpha \right) \) is sufficiently small. As a consequence, the first terms in Eqs. (3.16) and (3.17) are suppressed so as to be much smaller than all the other terms in these equations.

Provided that (at least) one of Cases 1 and 2 is realized, namely, the first terms in Eqs. (3.16) and (3.17) are negligible, by multiplying Eqs. (3.16) and (3.17) by \(u_0\) and \(\phi _0\), respectively, and summing them, we obtainFootnote 6

$$\begin{aligned} J'(y_0)=0. \end{aligned}$$
(3.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. (3.16) and (3.17), we see that in Case 1, the relation \(J(y_0)=-s^2/(288\alpha )\) 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 these two conditions. In the following, we study these two cases in more detail.

3.2.1 Case 1

In Case 1, by combining the value of the function J with the expression of the potential in Eq. (3.6), we have the effective inflaton potential

$$\begin{aligned} V_\mathrm {eff} (\lambda ) =\frac{1}{2\alpha }\left[ 1-\frac{s}{6}\mathrm {e}^{\lambda }(\phi _0^2-u_0^2)\right] . \end{aligned}$$
(3.19)

This is a kind of potential occurring in the modified Higgs inflation model [90]. The number of e-folds in the slow-roll regime is expressed as

$$\begin{aligned} N_{e}(\lambda )=\int ^{\lambda _\mathrm {f}}_{\lambda }H(\lambda )\frac{d\lambda }{\dot{\lambda }}=\frac{3}{2}\int ^{\lambda }_{\lambda _\mathrm {f}}\frac{V}{V_{\lambda }}d\lambda , \end{aligned}$$
(3.20)

with \(\lambda _\mathrm {f}\) the value of \(\lambda \) at the end of inflation. If \(\lambda \gg \lambda _\mathrm {f}\), the representation of \(N_{e}(\lambda )\) for the effective potential in Eq. (3.19) becomes

$$\begin{aligned} N_{e}(\lambda )=\frac{3}{2}\left[ \lambda +\frac{6}{s(\phi _0^2-u_0^2)\mathrm {e}^{\lambda }} \right] , \end{aligned}$$
(3.21)

or by taking into account the assumption that \(K \equiv \left( s/12\right) \mathrm {e}^{\lambda }(\phi _0^2-u_0^2)\gg 1\), we may rewrite \(N_{e}(\lambda )\) as

$$\begin{aligned} N_{e}(\phi )\approx \frac{3}{2} \lambda . \end{aligned}$$
(3.22)

With Eqs. (2.21) and (2.22), the slow-roll parameters can be described by using the effective potential as

$$\begin{aligned} \epsilon= & {} \frac{1}{3} \left( \frac{1}{V_\mathrm {eff}(\lambda )} \frac{dV_\mathrm {eff}(\lambda )}{d\lambda } \right) ^2 , \end{aligned}$$
(3.23)
$$\begin{aligned} \eta= & {} \frac{2}{3} \left( \frac{1}{V_\mathrm {eff}(\lambda )} \frac{d^2 V_\mathrm {eff}(\lambda )}{d \lambda ^2} \right) . \end{aligned}$$
(3.24)

Here, the coefficients of 1 / 3 in Eq. (3.23) and 2 / 3 in Eq. (3.24) originate from the non-canonical definition of the scalar field \(\lambda \); namely, the coefficient of the kinetic term for \(\lambda \) in the action in Eq. (3.5) is 3 / 2, and not 1 / 2. Eventually, \(n_\mathrm {s}\) and r are described by

$$\begin{aligned} n_\mathrm {s}= & {} \frac{\zeta ^2 \mathrm {e}^{2\lambda }-60\zeta \mathrm {e}^{\lambda }+108}{3\left( \zeta \mathrm {e}^{\lambda }-6\right) ^2},\end{aligned}$$
(3.25)
$$\begin{aligned} r= & {} \frac{16 \zeta ^2 \mathrm {e}^{2\lambda }}{\left( \zeta \mathrm {e}^{\lambda }-6\right) ^2}, \end{aligned}$$
(3.26)

where we have defined \(\zeta \) as \(\zeta \equiv s(\phi _0^2-u_0^2)\). Consequently, we see that \(n_\mathrm {s}\) can be written as a function of r as \(n_\mathrm {s}=n_\mathrm {s}(r)\). However, the values of the parameter \(\zeta \) lead to those of \(n_\mathrm {s}\) and r, which are compatible with the observations. Indeed, for \(N_e=60\) and \(\zeta =0.01\), we obtain \(r=0.004\) and \(n_\mathrm {s}=0.9617\). In this case, we find \(K=0.013\) (\({\ll } 1\)). For \(N_e=50\) and \(\zeta =0.18\) (the value of \(\zeta \) for \(N_e=50\) becomes larger than that for \(N_e = 60\) by one order of magnitude), we obtain \(r=0.005\) and \(n_\mathrm {s}=0.9568\). In this case, we have \(K=0.015\) (\({\ll } 1\)). As a result, we see that this case is inconsistent, because the initial assumption is \(K \gg 1\), but the consequences suggest \(K \ll 1\). It means that this case cannot be realized, in contrast with Case 2, as is shown next.

3.2.2 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 a kind of function J is

$$\begin{aligned} J(y)=C(y-y_0)^q , \quad q \ge 2 , \end{aligned}$$
(3.27)

where q is a constant. The effective potential is given by

$$\begin{aligned} V_\mathrm {eff} (\lambda ) = \frac{1}{2\alpha }\left( 1-\frac{\zeta }{12} \mathrm {e}^{\lambda }\right) ^2 , \end{aligned}$$
(3.28)

where we have used \(\zeta \equiv s(\phi _0^2-u_0^2)\). In addition, the expression of \(N_e\) reads

$$\begin{aligned} N_{e}(\phi )=\frac{3}{4}\lambda +\frac{9}{\zeta \mathrm {e}^{\lambda }} . \end{aligned}$$
(3.29)

Accordingly, we find

$$\begin{aligned} n_\mathrm {s}= & {} \frac{432 - 5 \zeta ^2 \mathrm {e}^{2\lambda }-168\zeta \mathrm {e}^{\lambda }}{3(\zeta \mathrm {e}^{\lambda }-12)^2},\end{aligned}$$
(3.30)
$$\begin{aligned} r= & {} \frac{64 \zeta ^2 \mathrm {e}^{2\lambda }}{3(\zeta \mathrm {e}^{\lambda }-12)^2} . \end{aligned}$$
(3.31)

Thus, it is seen that this theory also has only one parameter, and therefore we obtain \(n_\mathrm {s}=n_\mathrm {s}(r)\). Through the numerical calculations, the value of \(n_\mathrm {s}\) may be set near the value of \(n_\mathrm {s}=0.96\), suggested by the observations. In fact, in the case that \(N_e=60\), for \(\zeta = 10^{-1}\), \(10^{-2}\), \(10^{-3}\), \(10^{-4}\), \(10^{-5}\), and \(10^{-6}\), we have \(n_\mathrm {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 \(\zeta = 10^{-1}\), \(10^{-2}\), \(10^{-3}\), \(10^{-4}\), \(10^{-5}\), and \(10^{-6}\), we obtain \(n_\mathrm {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_\mathrm {s} \approx 0.96\) and \(r < 0.11\). These results are compatible with the observations of the Planck 2015 data.

Furthermore, we examine another aspect of this theory. If \(N_e=60\) with \(\zeta = 10^{-1}\), \(10^{-2}\), \(10^{-3}\), \(10^{-4}\), \(10^{-5}\), and \(10^{-6}\), or if \(N_e=50\) with these values of \(\zeta \), namely, the same combinations of values of \(N_e\) and \(\zeta \) shown above, we obtain the estimation as \(\left( \zeta /12\right) \mathrm {e}^{\lambda _\mathrm {f}} \thickapprox 0.14\)–0.17. It means that for all the combinations of values of \(N_e = 50\)–60 (during inflation) and those of \(\zeta \) described above, the value of the effective potential \(V_\mathrm {eff}\) may be estimated as

$$\begin{aligned} V_\mathrm {eff} \thickapprox \frac{0.7}{2\alpha } . \end{aligned}$$
(3.32)

This relation implies that for sufficiently large values of \(\alpha \), the values of the potential fall below the Planck scale. This fact is in good agreement with our initial assumption that \(s/\left( 6\alpha \right) \ll 1\). As a result, it is considered that Case 2 is quite viable.

4 Dynamical two scalar field model

In this section, we explore the theory proposed in Sect. 2, whose action is given by Eq. (2.1), and we consider how to realize inflation by using the two dynamical two scalar fields \(\phi \) and u.Footnote 7

The definitions of the slow-roll parameters in Eqs. (2.21) and (2.22) are used to describe the slow-roll inflation in the present theory. In addition, the spectral index \(n_\mathrm {s}\) of the curvature perturbations and the tensor-to-scalar ratio r are supposed to be represented as Eqs. (2.24) and (2.23), respectively. This assumption has been justified in Ref. [200].

We consider the FLRW background. We start with the system of the Eqs. (2.11) and (2.12), namely, (2.13) and (2.14) in the FLRW space-time. By imposing the standard slow-rolling conditions \(\ddot{u} \ll \dot{u} H\), \(\ddot{\phi }\ll \dot{\phi }H\), and \(\dot{\phi }^2-\dot{u}^2\ll H^2\), the Friedmann equation and the equations of motion for \(\phi \) and u are reduced to the following simple forms:

$$\begin{aligned}&H^2 = \frac{1}{6} V ,\end{aligned}$$
(4.1)
$$\begin{aligned}&3sH \dot{\phi } = V_{\phi } ,\end{aligned}$$
(4.2)
$$\begin{aligned}&3sH \dot{u} = -V_u . \end{aligned}$$
(4.3)

First, we derive the number of e-folds by using the initial definition,

$$\begin{aligned} N_{e}\equiv \int ^{a_\mathrm {f}}_{a_\mathrm {i}}d\ln a= \int ^{t_\mathrm {f}}_{t_\mathrm {i}}H\,\hbox {d}t , \end{aligned}$$
(4.4)

where \(a_\mathrm {i}\) and \(a_\mathrm {f}\) are the values of the scale factor a(t) at the beginning \(t_\mathrm {i}\) and end \(t_\mathrm {f}\) of inflation, respectively. Moreover, we also have the following relation:

$$\begin{aligned} dV=V_udu+V_{\phi }d\phi =V_u\dot{u}\,\hbox {d}t+ V_{\phi }\dot{\phi }\,\hbox {d}t . \end{aligned}$$
(4.5)

By using Eqs. (4.1)–(4.3), we get

$$\begin{aligned} N_{e}=\frac{s}{2}\int _{\phi ,u}^{\phi _\mathrm {f}, u_\mathrm {f}}\frac{V \left( V_u \hbox {d}u+V_{\phi }\hbox {d}\phi \right) }{V_{\phi }^2-V_u^2} , \end{aligned}$$
(4.6)

with \(\phi _\mathrm {f}\) and \(u_\mathrm {f}\) the amplitudes of \(\phi \) and u at the end of inflation, respectively. It follows from the definition in Eq. (2.21) that \(\epsilon \) is described as

$$\begin{aligned} \epsilon =-\frac{\dot{V}}{2HV}= \frac{V_u^2 - V_{\phi }^2}{sV^2} . \end{aligned}$$
(4.7)

On the other hand, with the definition in Eq. (2.22) and Eq. (4.1), we find

$$\begin{aligned} \eta =-\frac{1}{4HV\dot{V}} \left( \dot{V}^2 + 2\ddot{V}V \right) =-\frac{2\left( V_{\phi }^2V_{\phi \phi }+V_u^2V_{uu}\right) }{sV\left( V_{\phi }^2-V_u^2 \right) },\nonumber \\ \end{aligned}$$
(4.8)

where in deriving the last equality, we have used Eqs. (4.2) and (4.3). The expressions (4.6)–(4.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, Eq. (4.6) leads only to a relation between the scalar fields \(\phi \) and u, and not both values of \(\phi \) and u simultaneously. Hence, another additional assumption to determine both values, of \(\phi \) and u (e.g., \(\phi \approx u\)), is necessary, or one of the values of \(\phi \) and u may be regarded as a free parameter.

In the simplest case \(J=C\), it is possible to obtain the values of the main parameters analytically.Footnote 8 In this case, from Eq. (4.6), we find

$$\begin{aligned} N_{e}&=\frac{s}{16}(u^2-\phi ^2)+\frac{3}{2}\ln (\phi ^2-u^2)\nonumber \\&\quad -\frac{432x}{288x+s^2}\ln \left[ \left( 288x+s^2\right) \left( \phi ^2-u^2\right) -12s\right] . \end{aligned}$$
(4.9)

The expressions for \(\epsilon \) and \(\eta \) take the form

$$\begin{aligned} \epsilon= & {} -\frac{16(\phi ^2-u^2)}{s}\left\{ \frac{(288x+s^2)(\phi ^2-u^2)-12s}{288x(\phi ^2-u^2)^2+\left[ 12-s(\phi ^2-u^2)\right] ^2} \right\} ^2 ,\nonumber \\ \end{aligned}$$
(4.10)
$$\begin{aligned} \eta= & {} - \left[ \frac{8}{s\left( \phi ^2-u^2\right) } \right] \nonumber \\&\times \frac{\left( 288x+s^2\right) \left( 3\phi ^4-2\phi ^2u^2+3u^4\right) -12s\left( \phi ^2-u^2\right) }{288x\left( \phi ^2-u^2\right) ^2+\left[ 12-s(\phi ^2-u^2)\right] ^2}.\nonumber \\ \end{aligned}$$
(4.11)

With the expression of \(N_e\) in (4.9), the value of \((\phi ^2-u^2)\) can be obtained. It determines the value of \(\epsilon \), but that of \(\eta \) cannot be found by using only that of \((\phi ^2-u^2)\). Since \(\epsilon (> 0)\) should be positive, from Eqs. (4.9) and (4.10), we have \(\phi ^2-u^2>0\) and \(s<0\). These relations and \(x = \alpha C >0\) lead to a positive value of \(\eta \).

We estimate numerical values. For simplicity, we take \(s=-1\). From Eqs. (4.9) and (4.10), we get \((\phi ^2-u^2)\) and x for any interesting values of \(N_e\) and \(r=16\epsilon \). 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_\mathrm {s} = n_\mathrm {s} (r, \mathcal {G}^4)\), where we have defined a new parameter \(\mathcal {G}^4 \equiv u^2\phi ^2\). This relation yields a value of \(n_\mathrm {s}\) consistent with the observations. Indeed, for \(N_e=50\) and \(x=10^{30}\), we have \(\phi ^2-u^2=2593\) and \(r\simeq 0.1\). In this case, \(n_\mathrm {s}\) is represented as

$$\begin{aligned} n_\mathrm {s}=3.7 \times 10^{-9} \mathcal {G}^4+0.9815 . \end{aligned}$$
(4.12)

Moreover, for \(N_e=60\) and \(x=10^{25}\), we obtain \(\phi ^2-u^2=2317\) and \(r\simeq 0.11\). In this case, \(n_\mathrm {s}\) is expressed as

$$\begin{aligned} n_\mathrm {s} = 5.0 \times 10^{-9} \mathcal {G}^4 + 0.9793 . \end{aligned}$$
(4.13)

The larger x is, the smaller r becomes, but the minimum value of \(n_\mathrm {s}\) in both Eqs. (4.12) and (4.13) increases. Equations (4.12) and (4.13) imply that it is impossible for the pair of \(n_\mathrm {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 \(\phi ^2-u^2=1267\), \(r\simeq 0.202\), and \(n_\mathrm {s}=3 \times 10^{-8}\mathcal {G}^4+0.9621\). Therefore, the case of \(\mathcal {G}^4 < 10^5\) is quite viable and consistent 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 a single scalar field. There is another possibility for considering the more viable types of the function J, e.g., Eq. (2.29). In such a case, it is quite difficult to analyze the equations analytically (for example, \(N_e\) is described as an integral equation).

5 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 Sect. 3.2.1, 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

$$\begin{aligned} H = H_\mathrm {inf} \left( 1 + \delta (t) \right) , \quad \left| \delta (t) \right| \ll 1 \end{aligned}$$
(5.1)

where \(H_\mathrm {inf} (>0)\) (\(=\)constant) is the Hubble parameter during inflation (whose value is positive), and \(\delta (t)\) means the perturbations from the de Sitter solution. In the present case, only the scalar field \(\lambda \) is dynamical, and \(\phi \) and u are constant. By taking the time derivative of the gravitational field equation (3.15) with \(V = V_\mathrm {eff}\), where \(V_\mathrm {eff}\) is given by Eq. (3.28), we have

$$\begin{aligned} 2\ddot{H} +6H\dot{H} +\frac{3}{2} \dot{\lambda } \ddot{\lambda } + \frac{\zeta }{24\alpha } \left( 1-\frac{\zeta }{12} \mathrm {e}^{\lambda } \right) \mathrm {e}^{\lambda } \dot{\lambda } = 0 . \end{aligned}$$
(5.2)

Moreover, from Eq. (3.11) with \(V = V_\mathrm {eff}\) [in Eq. (3.28)], the equation of motion for \(\lambda \) reads

$$\begin{aligned} 3\ddot{\lambda } + 9H\dot{\lambda } - \frac{\zeta }{12\alpha } \left( 1-\frac{\zeta }{12} \mathrm {e}^{\lambda } \right) \mathrm {e}^{\lambda } =0. \end{aligned}$$
(5.3)

Since we consider the solution at the inflationary stage in the early universe, we take the limit \(t \rightarrow 0\). In this limit, we could have an approximate solution \(\lambda \approx \ln \left( H t \right) \) for Eq. (5.3) with the quasi de Sitter solution \(H \approx 1/\left( 3t \right) \). By taking into account this solution, we see that in the limit \(t \rightarrow 0\), the third term, proportional to \(\ddot{\lambda }\) in the left-hand side of Eq. (5.2), is much larger than the fourth term, which is proportional only to \(\dot{\lambda }\), and therefore the fourth term could be neglected.

To examine the instability of the de Sitter solution, we represent \(\delta (t)\) as

$$\begin{aligned} \delta (t) = \mathrm {e}^{\beta t} , \end{aligned}$$
(5.4)

with \(\beta \) a constant. When \(\beta \) 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 \(\beta > 0\), the amplitude of \(\delta (t)\) increases in time. By substituting the expression of the perturbations in Eq. (5.1) with Eq. (5.4) into Eq. (5.2) and using approximate solutions \(\lambda \approx \ln \left( H t \right) \) and \(H \approx 1/\left( 3t \right) \), we find

$$\begin{aligned} 2H_\mathrm {inf} \beta ^2 + 6H_\mathrm {inf}^2 \beta + \frac{81}{2} H_\mathrm {inf}^3 = 0. \end{aligned}$$
(5.5)

The solutions for this equation are obtained as

$$\begin{aligned} \beta _\pm = \frac{3\left( -1\pm \sqrt{10}\right) H_\mathrm {inf}}{2}. \end{aligned}$$
(5.6)

As a consequence, we obtain the positive solution of \(\beta = \beta _+ > 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 may exist a solution representing 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 \(\lambda \) and its dynamics through the auxiliary field \(\Phi \) in the Einstein frame. This fact can be seen from the action in Eq. (2.2), and Eq. (3.3) with \(\Lambda = \mathrm {e}^{\lambda }\). 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.

6 Conclusions

In the present paper, we have studied inflationary cosmology in a theory where there exist two scalar fields non-minimally coupled 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_\mathrm {s}\) of the 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 a single dynamical scalar field in the Jordan frame, if the number of e-folds during inflation is \(50 \le N_e \le 60\), we have \(n_\mathrm {s} \approx 0.96\) and \(r = \mathcal {O} (0.1)\). On the other hand, in the Einstein frame, for the case of one dynamical scalar field, we can obtain \(n_\mathrm {s} \approx 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 \le N_e \le 60\), we obtain \(n_\mathrm {s} \approx 0.96\) and \(r = \mathcal {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 seem to be more in correspondence with the observational data obtained from the Planck satellite than single field inflation models.

It is remarked that the action in Eq. (2.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 come 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. [201], a theory without an \(R^2\) term and that with the general potential \((\phi ^2 - u^2 - \psi ^2)^2 J(u/\phi , u/\psi )\) have been explored. The action for the theory with the more specific potential \(\psi ^4 J(u/\phi , u/\psi )\) and without the \(R^2\) term) is given by

$$\begin{aligned} S&=\int d^4x\sqrt{-g}\left\{ \frac{s}{2}\left[ \frac{(\phi ^2-u^2 +\epsilon \psi ^2)}{6}R + (\nabla \phi )^2 - (\nabla u)^2\right. \right. \nonumber \\&\quad +\left. \left. \epsilon (\nabla \psi )^2\right] - \psi ^4J\left( \frac{u}{\phi },\frac{u}{\psi }\right) \right\} , \end{aligned}$$
(6.1)

where \(\epsilon \) is a constant and \(\psi \) is the additional third scalar field to two scalar fields \(\phi \) and u already introduced. If the \(R^2\) term is added to Eq. (6.1), the representation of the action reads

$$\begin{aligned} S&=\int d^4x\sqrt{-g}\left\{ \frac{\alpha }{2} R^2 + \frac{s}{2}\left[ \frac{(\phi ^2-u^2 +\epsilon \psi ^2)}{6}R\right. \right. \nonumber \\&\quad +\left. \left. (\nabla \phi )^2 - (\nabla u)^2 +\epsilon (\nabla \psi )^2 \right] - \psi ^4J\left( \frac{u}{\phi },\frac{u}{\psi }\right) \right\} . \end{aligned}$$
(6.2)

Moreover, when the general form of the potential is included and the signature of \(\epsilon \) is opposite to that in the action in Eq. (6.2), the action is represented as

$$\begin{aligned} S&=\int d^4x\sqrt{-g}\left\{ \frac{\alpha }{2} R^2 + \frac{s}{2}\left[ \frac{(\phi ^2-u^2 -\epsilon \psi ^2)}{6}R\right. \right. \nonumber \\&\quad +\left. \left. (\nabla \phi )^2 - (\nabla u)^2 -\epsilon (\nabla \psi )^2 \right] \!-\! (\phi ^2 - u^2 - \epsilon \psi ^2)^2J\left( \frac{u}{\phi },\frac{u}{\psi }\right) \right\} . \end{aligned}$$
(6.3)

The same analysis as in Sects. 24 may be adopted for the models in Eqs. (6.1)–(6.3), in which the additional third scalar field \(\psi \) may play the role of the other species of dark matter, or make an additional contribution to the dark energy dominated stage or inflation.