A focused review of quintom cosmology: from quintom dark energy to quintom bounce

Figures(13)

Get Citation
Tao-tao Qiu, Yifu Cai, Yang Liu, Si-Yu Li, Jarah Evslin and Xinmin Zhang. A focused review of quintom cosmology: from quintom dark energy to quintom bounce[J]. Chinese Physics C, 2026, 50(1): 012001. doi: 10.1088/1674-1137/ae30e9
Tao-tao Qiu, Yifu Cai, Yang Liu, Si-Yu Li, Jarah Evslin and Xinmin Zhang. A focused review of quintom cosmology: from quintom dark energy to quintom bounce[J]. Chinese Physics C, 2026, 50(1): 012001.  doi: 10.1088/1674-1137/ae30e9 shu
Milestone
Received: 2025-12-01
Article Metric

Article Views(12)
PDF Downloads(0)
Cited by(0)
Policy on re-use
To reuse of subscription content published by CPC, the users need to request permission from CPC, unless the content was published under an Open Access license which automatically permits that type of reuse.
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Email This Article

Title:
Email:

A focused review of quintom cosmology: from quintom dark energy to quintom bounce

  • 1. School of Physics, Huazhong University of Science and Technology, Wuhan 430074, China
  • 2. Department of Astronomy, School of Physical Sciences, University of Science and Technology of China, Hefei 230026, China
  • 3. CAS Key Laboratory for Research in Galaxies and Cosmology, School of Astronomy and Space Science, University of Science and Technology of China, Hefei 230026, China
  • 4. Key Laboratory of Particle Astrophysics, Institute of High Energy Physics, Chinese Academy of Sciences, Beijing 100049, China
  • 5. Institute of Modern Physics, Lanzhou 730000, China
  • 6. University of the Chinese Academy of Sciences, Beijing 100049, China
  • 7. Theoretical Physics Division, Institute of High Energy Physics, Chinese Academy of Sciences, Beijing 100049, China
  • 8. School of Nuclear Science and Technology, University of Chinese Academy of Sciences, Beijing 101408, China

Abstract: The recently released data of DESI DR2 favors a dynamical dark energy theory, with the equation of state crossing the cosmological constant boundary w = −1. In this paper, we briefly review quintom cosmology, especially the quintom bounce. We will give three examples of a quintom bounce and one example of a cyclic universe with quintom matter.

    HTML

  • The accelerated expansion of the Universe was first discovered in 1998 through observations of high-redshift Type Ia supernovae (SNe) [1, 2], and has since been robustly confirmed by cosmic microwave background (CMB) anisotropies, baryon acoustic oscillations (BAO), and large-scale structure data. Within the standard cosmological model, this acceleration is attributed to a cosmological constant Λ, yielding the highly successful ΛCDM paradigm for past 20 years.

    Recent results from Dark Energy Spectroscopic Instrument (DESI) Data Release 2, when combined with CMB and supernova data, provide compelling evidence that dark energy may be dynamical in nature [3, 4]. The joint analysis reveals a statistically significant time evolution of the dark energy equation of state (EoS), w, with the DESI DR2+CMB+DESY5 dataset excluding the ΛCDM model at the $ 4.2\sigma $ level. Notably, the reconstructed w crosses the cosmological constant boundary at $ w = -1 $, a distinctive signature of the Quintom scenario [5].

    To quantify the statistical support for the Quintom scenario, we perform an independent cosmological analysis using the same combined DESI DR2+CMB+DESY5 dataset. Here, CMB includes the full Planck temperature (TT), polarization (EE), and cross (TE) power spectra (using $ \mathtt{simall} $, $ \mathtt{Commander} $ for $ \ell < 30 $, and $ \mathtt{CamSpec} $ for $ \ell \geqslant 30 $), supplemented by the Planck+ACT DR6 CMB lensing likelihood. Using the $ \mathtt{Cobaya} $ inference framework with Metropolis–Hastings MCMC sampling, interfaced with $ \mathtt{CAMB} $ for theory predictions, we adopt the CPL parametrization $ (w_0, w_a) $ [6, 7] for dark energy, and the standard six-parameter base ($ \omega_b $, $ \omega_c $, $ 100\theta_{\rm{MC}} $, $ \ln(10^{10}A_s) $, $ n_s $, τ). Chains are deemed converged when the Gelman–Rubin statistic satisfies $ R - 1 < 0.01 $, and posterior summaries are generated using the $ \mathtt{GetDist} $ package.

    Our result yields $ \Delta \chi_{{\rm{MAP}}}^2 = -21.2 $ relative to the ΛCDM model, corresponding to a $ 4.22\,\sigma $ deviation (In a recent reanalysis by the DES group combined with CMB data, this confidence level has been reduced, but a $ 3.2\,\sigma $ preference for evolving dark energy remains [8]). Moreover, we find a posterior probability of 99.997439% for the Quintom-B region defined by $ w_0 > -1 $ and $ w_0 + w_a < -1 $, equivalent to a significance of approximately $ 4.05\,\sigma $. As shown in Fig. 1, our results are fully consistent with the DESI findings [3] and very strongly support the Quintom scenario. 1

    Figure 1.  (color online) Results for the posterior distributions of $ w_0 $ and $ w_a $, from fits of the $ w_0w_aCDM $ model to DESI in combination with CMB and DESY5 datasets as labelled. Quintom-A: w evolves from Quintessence to Phantom; Quintom-B: from Phantom to Quintessence.

    The Quintom scenario was first proposed in April of 2004 [5]. The Quintom theory of dark energy differs from Quintessence and the Phantom, and predicts the EoS of dark energy crossing 2 the cosmological constant boundary at $ w = -1 $. In general, the basic single-field models or single perfect fluid models can not exhibit Quintom behavior due to the No-Go theorem [1113]. For a detailed description of the Quintom model, please refer to Refs. [1419].

    In 2007, we considered the applications of Quintom scenarios on the very early universe and demonstrated that non-singular bounce models can emerge naturally with quintom matter [20]. As we know, non-singular bounce cosmology encompasses scenarios in which the Universe initially undergoes a phase of contraction, reaches a finite minimum scale, then subsequently enters an expanding phase. This behavior can be understood by examining the dynamics of the scale factor $ a(t) $ and the Hubble parameter $ H(t) $. It is well-known that in GR, the Friedmann equations that govern the universe are written as:

    $ H^2=\frac{8\pi G}{3}\rho\; ,\; \; \; \dot H=-4\pi G(\rho+p)\; , $

    (1)

    where $ H={\rm d}\ln a/{\rm d}t $. When the Universe goes from the contracting phase to the expanding phase, the Hubble parameter goes from $ H<0 $ to $ H>0 $, with $ H=\rho=0 $ at the bounce point. This requires that the time derivative of H must be larger than zero, $ \dot H>0 $, which gives $ \rho+p<0 $. Therefore, at this point, one has $ w\equiv p/\rho\rightarrow -\infty $, namely below the $ w=-1 $ line. On the other hand, after the universe enters into the normal expanding phase, it should experience radiation-dominant era ($ w=1/3 $), matter-dominant era ($ w=0 $), etc., as required by the standard Big-Bang Theory. This means that in the whole bounce process, the equation of state of our universe will cross -1 from below. Such a crossing behavior is the characteristic property of the quintom model.

    After the bounce, the Universe will enter into an expanding phase, and according to its behavior, in the future it will evolve towards various fates [17]. Among these fates, an interesting one is that the Universe will return to a contracting phase and then bounce again, thus performing a cyclic behavior forever [2127]. Like the bounce, a 4D non-singular cyclic universe also needs its EoS to cross -1, and thus it can be realized by quintom matter [13, 28, 29].

    In the following, we will give three models as examples of a quintom bounce, and one example of a cyclic universe with quintom matter.

  • The accelerated expansion of the Universe was first discovered in 1998 through observations of high-redshift Type Ia supernovae (SNe) [1, 2], and has since been robustly confirmed by cosmic microwave background (CMB) anisotropies, baryon acoustic oscillations (BAO), and large-scale structure data. Within the standard cosmological model, this acceleration is attributed to a cosmological constant Λ, yielding the highly successful ΛCDM paradigm for past 20 years.

    Recent results from Dark Energy Spectroscopic Instrument (DESI) Data Release 2, when combined with CMB and supernova data, provide compelling evidence that dark energy may be dynamical in nature [3, 4]. The joint analysis reveals a statistically significant time evolution of the dark energy equation of state (EoS), w, with the DESI DR2+CMB+DESY5 dataset excluding the ΛCDM model at the $ 4.2\sigma $ level. Notably, the reconstructed w crosses the cosmological constant boundary at $ w = -1 $, a distinctive signature of the Quintom scenario [5].

    To quantify the statistical support for the Quintom scenario, we perform an independent cosmological analysis using the same combined DESI DR2+CMB+DESY5 dataset. Here, CMB includes the full Planck temperature (TT), polarization (EE), and cross (TE) power spectra (using $ \mathtt{simall} $, $ \mathtt{Commander} $ for $ \ell < 30 $, and $ \mathtt{CamSpec} $ for $ \ell \geqslant 30 $), supplemented by the Planck+ACT DR6 CMB lensing likelihood. Using the $ \mathtt{Cobaya} $ inference framework with Metropolis–Hastings MCMC sampling, interfaced with $ \mathtt{CAMB} $ for theory predictions, we adopt the CPL parametrization $ (w_0, w_a) $ [6, 7] for dark energy, and the standard six-parameter base ($ \omega_b $, $ \omega_c $, $ 100\theta_{\rm{MC}} $, $ \ln(10^{10}A_s) $, $ n_s $, τ). Chains are deemed converged when the Gelman–Rubin statistic satisfies $ R - 1 < 0.01 $, and posterior summaries are generated using the $ \mathtt{GetDist} $ package.

    Our result yields $ \Delta \chi_{{\rm{MAP}}}^2 = -21.2 $ relative to the ΛCDM model, corresponding to a $ 4.22\,\sigma $ deviation (In a recent reanalysis by the DES group combined with CMB data, this confidence level has been reduced, but a $ 3.2\,\sigma $ preference for evolving dark energy remains [8]). Moreover, we find a posterior probability of 99.997439% for the Quintom-B region defined by $ w_0 > -1 $ and $ w_0 + w_a < -1 $, equivalent to a significance of approximately $ 4.05\,\sigma $. As shown in Fig. 1, our results are fully consistent with the DESI findings [3] and very strongly support the Quintom scenario. 1

    Figure 1.  (color online) Results for the posterior distributions of $ w_0 $ and $ w_a $, from fits of the $ w_0w_aCDM $ model to DESI in combination with CMB and DESY5 datasets as labelled. Quintom-A: w evolves from Quintessence to Phantom; Quintom-B: from Phantom to Quintessence.

    The Quintom scenario was first proposed in April of 2004 [5]. The Quintom theory of dark energy differs from Quintessence and the Phantom, and predicts the EoS of dark energy crossing 2 the cosmological constant boundary at $ w = -1 $. In general, the basic single-field models or single perfect fluid models can not exhibit Quintom behavior due to the No-Go theorem [1113]. For a detailed description of the Quintom model, please refer to Refs. [1419].

    In 2007, we considered the applications of Quintom scenarios on the very early universe and demonstrated that non-singular bounce models can emerge naturally with quintom matter [20]. As we know, non-singular bounce cosmology encompasses scenarios in which the Universe initially undergoes a phase of contraction, reaches a finite minimum scale, then subsequently enters an expanding phase. This behavior can be understood by examining the dynamics of the scale factor $ a(t) $ and the Hubble parameter $ H(t) $. It is well-known that in GR, the Friedmann equations that govern the universe are written as:

    $ H^2=\frac{8\pi G}{3}\rho\; ,\; \; \; \dot H=-4\pi G(\rho+p)\; , $

    (1)

    where $ H={\rm d}\ln a/{\rm d}t $. When the Universe goes from the contracting phase to the expanding phase, the Hubble parameter goes from $ H<0 $ to $ H>0 $, with $ H=\rho=0 $ at the bounce point. This requires that the time derivative of H must be larger than zero, $ \dot H>0 $, which gives $ \rho+p<0 $. Therefore, at this point, one has $ w\equiv p/\rho\rightarrow -\infty $, namely below the $ w=-1 $ line. On the other hand, after the universe enters into the normal expanding phase, it should experience radiation-dominant era ($ w=1/3 $), matter-dominant era ($ w=0 $), etc., as required by the standard Big-Bang Theory. This means that in the whole bounce process, the equation of state of our universe will cross -1 from below. Such a crossing behavior is the characteristic property of the quintom model.

    After the bounce, the Universe will enter into an expanding phase, and according to its behavior, in the future it will evolve towards various fates [17]. Among these fates, an interesting one is that the Universe will return to a contracting phase and then bounce again, thus performing a cyclic behavior forever [2127]. Like the bounce, a 4D non-singular cyclic universe also needs its EoS to cross -1, and thus it can be realized by quintom matter [13, 28, 29].

    In the following, we will give three models as examples of a quintom bounce, and one example of a cyclic universe with quintom matter.

  • The accelerated expansion of the Universe was first discovered in 1998 through observations of high-redshift Type Ia supernovae (SNe) [1, 2], and has since been robustly confirmed by cosmic microwave background (CMB) anisotropies, baryon acoustic oscillations (BAO), and large-scale structure data. Within the standard cosmological model, this acceleration is attributed to a cosmological constant Λ, yielding the highly successful ΛCDM paradigm for past 20 years.

    Recent results from Dark Energy Spectroscopic Instrument (DESI) Data Release 2, when combined with CMB and supernova data, provide compelling evidence that dark energy may be dynamical in nature [3, 4]. The joint analysis reveals a statistically significant time evolution of the dark energy equation of state (EoS), w, with the DESI DR2+CMB+DESY5 dataset excluding the ΛCDM model at the $ 4.2\sigma $ level. Notably, the reconstructed w crosses the cosmological constant boundary at $ w = -1 $, a distinctive signature of the Quintom scenario [5].

    To quantify the statistical support for the Quintom scenario, we perform an independent cosmological analysis using the same combined DESI DR2+CMB+DESY5 dataset. Here, CMB includes the full Planck temperature (TT), polarization (EE), and cross (TE) power spectra (using $ \mathtt{simall} $, $ \mathtt{Commander} $ for $ \ell < 30 $, and $ \mathtt{CamSpec} $ for $ \ell \geqslant 30 $), supplemented by the Planck+ACT DR6 CMB lensing likelihood. Using the $ \mathtt{Cobaya} $ inference framework with Metropolis–Hastings MCMC sampling, interfaced with $ \mathtt{CAMB} $ for theory predictions, we adopt the CPL parametrization $ (w_0, w_a) $ [6, 7] for dark energy, and the standard six-parameter base ($ \omega_b $, $ \omega_c $, $ 100\theta_{\rm{MC}} $, $ \ln(10^{10}A_s) $, $ n_s $, τ). Chains are deemed converged when the Gelman–Rubin statistic satisfies $ R - 1 < 0.01 $, and posterior summaries are generated using the $ \mathtt{GetDist} $ package.

    Our result yields $ \Delta \chi_{{\rm{MAP}}}^2 = -21.2 $ relative to the ΛCDM model, corresponding to a $ 4.22\,\sigma $ deviation (In a recent reanalysis by the DES group combined with CMB data, this confidence level has been reduced, but a $ 3.2\,\sigma $ preference for evolving dark energy remains [8]). Moreover, we find a posterior probability of 99.997439% for the Quintom-B region defined by $ w_0 > -1 $ and $ w_0 + w_a < -1 $, equivalent to a significance of approximately $ 4.05\,\sigma $. As shown in Fig. 1, our results are fully consistent with the DESI findings [3] and very strongly support the Quintom scenario. 1

    Figure 1.  (color online) Results for the posterior distributions of $ w_0 $ and $ w_a $, from fits of the $ w_0w_aCDM $ model to DESI in combination with CMB and DESY5 datasets as labelled. Quintom-A: w evolves from Quintessence to Phantom; Quintom-B: from Phantom to Quintessence.

    The Quintom scenario was first proposed in April of 2004 [5]. The Quintom theory of dark energy differs from Quintessence and the Phantom, and predicts the EoS of dark energy crossing 2 the cosmological constant boundary at $ w = -1 $. In general, the basic single-field models or single perfect fluid models can not exhibit Quintom behavior due to the No-Go theorem [1113]. For a detailed description of the Quintom model, please refer to Refs. [1419].

    In 2007, we considered the applications of Quintom scenarios on the very early universe and demonstrated that non-singular bounce models can emerge naturally with quintom matter [20]. As we know, non-singular bounce cosmology encompasses scenarios in which the Universe initially undergoes a phase of contraction, reaches a finite minimum scale, then subsequently enters an expanding phase. This behavior can be understood by examining the dynamics of the scale factor $ a(t) $ and the Hubble parameter $ H(t) $. It is well-known that in GR, the Friedmann equations that govern the universe are written as:

    $ H^2=\frac{8\pi G}{3}\rho\; ,\; \; \; \dot H=-4\pi G(\rho+p)\; , $

    (1)

    where $ H={\rm d}\ln a/{\rm d}t $. When the Universe goes from the contracting phase to the expanding phase, the Hubble parameter goes from $ H<0 $ to $ H>0 $, with $ H=\rho=0 $ at the bounce point. This requires that the time derivative of H must be larger than zero, $ \dot H>0 $, which gives $ \rho+p<0 $. Therefore, at this point, one has $ w\equiv p/\rho\rightarrow -\infty $, namely below the $ w=-1 $ line. On the other hand, after the universe enters into the normal expanding phase, it should experience radiation-dominant era ($ w=1/3 $), matter-dominant era ($ w=0 $), etc., as required by the standard Big-Bang Theory. This means that in the whole bounce process, the equation of state of our universe will cross -1 from below. Such a crossing behavior is the characteristic property of the quintom model.

    After the bounce, the Universe will enter into an expanding phase, and according to its behavior, in the future it will evolve towards various fates [17]. Among these fates, an interesting one is that the Universe will return to a contracting phase and then bounce again, thus performing a cyclic behavior forever [2127]. Like the bounce, a 4D non-singular cyclic universe also needs its EoS to cross -1, and thus it can be realized by quintom matter [13, 28, 29].

    In the following, we will give three models as examples of a quintom bounce, and one example of a cyclic universe with quintom matter.

  • I.   BOUNCE WITH TWO SCALAR FIELDS
    • The simplest Quintom model contains a quintessence-like scalar field and a phantom-like scalar field, with the Lagrangian [5, 14, 20] 3:

      $ S=\int {\rm d}^4x\sqrt{-g}\left[-\frac{1}{2}\nabla_\mu\phi\nabla^\mu\phi-V(\phi)+\frac{1}{2}\nabla_\mu\sigma\nabla^\mu\sigma-V(\sigma)\right]\; , $

      (2)

      where ϕ and σ stand for quintessence and phantom fields, respectively. The Friedmann equations are:

      $ \begin{aligned} H^2=&\ \frac{8\pi G}{3}\left[\frac{1}{2}\dot\phi^2-\frac{1}{2}\dot\sigma^2+V(\phi)+V(\sigma)\right]\; , \end{aligned} $

      (3)

      $ \begin{aligned} \dot H=&\ -4\pi G(\dot\phi^2-\dot\sigma^2)\; , \end{aligned} $

      (4)

      while the equations of motion of the two fields are:

      $ \begin{aligned} \ddot\phi+3H\dot\phi+V_{,\phi}= 0\; , \end{aligned} $

      (5)

      $ \begin{aligned} \ddot\sigma+3H\dot\sigma+V_{,\sigma}=&\ 0\; . \end{aligned} $

      (6)

      From the Friedmann equation one can see that, provided that the potentials of both fields are not less than zero, the only possibility for making $ H=0 $, as is required by bounce process, is to have the kinetic term of the σ-field become large and compensate the rest of the energy density. On the other hand, this term should not be dominant in order not to give rise to a negative energy density. This can be achieved by setting $ V(\sigma)=0 $, so that the equation of motion for the σ-field becomes $ \ddot\sigma+3H\dot\sigma=0 $ with the solution $ \dot\sigma^2\sim a^{-6} $. Therefore at the bounce point where the scale factor $ a(t) $ reaches its minimum, $ \dot\sigma^2 $ gets its maximum value, while whether the universe goes backwards or forwards, $ a(t) $ gets enlarged and $ \dot\sigma^2 $ will decay. Thus the evolution of the universe away from the bounce depends on the potential of the normal field ϕ only.

      To obtain the observational signals, it is important to analyze the perturbations generated by the model. The longitudinal metric perturbations are given by:

      $ {\rm d}s^2=a^2(\eta)[(1+\Phi){\rm d}\eta^2-(1-2\Psi){\rm d}x^i{\rm d}x^j]\; , $

      (7)

      where η is the conformal time with $ {\rm d}\eta=a^{-1}(t){\rm d}t $, while Φ and Ψ are gauge-invariant metric perturbations. The perturbed Einstein equation gives rise to the equations of Φ and Ψ [3134]:

      $ \begin{split} \nabla^2\Psi-3{\cal{H}}(\Psi'+{\cal{H}}\Phi)=\;&\ 4\pi G[\phi'(\delta\phi'-\phi'\Phi)\\ &+a^2V_{,\phi}\delta\phi-\sigma'(\delta\sigma'-\sigma'\Phi)]\; , \end{split} $

      (8)

      $ \begin{aligned} \Psi'+{\cal{H}}\Phi=&\ -4\pi G [-\phi'\delta\phi+\sigma'\delta\sigma]\; , \end{aligned} $

      (9)

      $ \begin{split} \Phi''+3{\cal{H}}\Phi'+(2{\cal{H}}'+{\cal{H}}^2)\Phi=\;&\ 4\pi G [\phi'(\delta\phi'-\phi'\Phi)\\ &-a^2V_{,\phi}\delta\phi-\sigma'(\delta\sigma'-\sigma'\Phi)]\; , \end{split} $

      (10)

      which can be combined to give the main equation of the Newtonian potential Φ:

      $\begin{split} &\Phi''+2\left({\cal{H}}-\frac{\phi''}{\phi'}\right)\Phi'+2\left({\cal{H}}'-{\cal{H}}\frac{\phi''}{\phi'}\right)\Phi-\nabla^2\Phi\\ &\quad =8\pi G\left(2{\cal{H}}+\frac{\phi''}{\phi'}\right)\sigma'\delta\sigma\; ,\end{split} $

      (11)

      where prime denotes derivative with respect to the conformal time η, and $ {\cal{H}} $ is the conformal Hubble parameter, $ {\cal{H}}\equiv a'/a $. Moreover, the curvature fluctuation ζ is related with Φ via:

      $ \zeta=\Phi+\frac{\cal{H}}{{\cal{H}}^2-{\cal{H}}'}(\Phi'+{\cal{H}}\Phi)\; , $

      (12)

      which is conserved on super-Hubble scales:

      $ (1+w)\zeta'=\frac{2k^2(\Phi'+{\cal{H}}\Phi)}{9{\cal{H}}^2}\; . $

      (13)

      In the following, we will see that different forms of potential $ V(\phi) $ will give rise to different backgrounds as well as different evolutions of the perturbations. In order to understand their main influence, we consider two broad categories of potentials, namely large field potentials and small field potentials.

    • I.   BOUNCE WITH TWO SCALAR FIELDS
      • The simplest Quintom model contains a quintessence-like scalar field and a phantom-like scalar field, with the Lagrangian [5, 14, 20] 3:

        $ S=\int {\rm d}^4x\sqrt{-g}\left[-\frac{1}{2}\nabla_\mu\phi\nabla^\mu\phi-V(\phi)+\frac{1}{2}\nabla_\mu\sigma\nabla^\mu\sigma-V(\sigma)\right]\; , $

        (2)

        where ϕ and σ stand for quintessence and phantom fields, respectively. The Friedmann equations are:

        $ \begin{aligned} H^2=&\ \frac{8\pi G}{3}\left[\frac{1}{2}\dot\phi^2-\frac{1}{2}\dot\sigma^2+V(\phi)+V(\sigma)\right]\; , \end{aligned} $

        (3)

        $ \begin{aligned} \dot H=&\ -4\pi G(\dot\phi^2-\dot\sigma^2)\; , \end{aligned} $

        (4)

        while the equations of motion of the two fields are:

        $ \begin{aligned} \ddot\phi+3H\dot\phi+V_{,\phi}= 0\; , \end{aligned} $

        (5)

        $ \begin{aligned} \ddot\sigma+3H\dot\sigma+V_{,\sigma}=&\ 0\; . \end{aligned} $

        (6)

        From the Friedmann equation one can see that, provided that the potentials of both fields are not less than zero, the only possibility for making $ H=0 $, as is required by bounce process, is to have the kinetic term of the σ-field become large and compensate the rest of the energy density. On the other hand, this term should not be dominant in order not to give rise to a negative energy density. This can be achieved by setting $ V(\sigma)=0 $, so that the equation of motion for the σ-field becomes $ \ddot\sigma+3H\dot\sigma=0 $ with the solution $ \dot\sigma^2\sim a^{-6} $. Therefore at the bounce point where the scale factor $ a(t) $ reaches its minimum, $ \dot\sigma^2 $ gets its maximum value, while whether the universe goes backwards or forwards, $ a(t) $ gets enlarged and $ \dot\sigma^2 $ will decay. Thus the evolution of the universe away from the bounce depends on the potential of the normal field ϕ only.

        To obtain the observational signals, it is important to analyze the perturbations generated by the model. The longitudinal metric perturbations are given by:

        $ {\rm d}s^2=a^2(\eta)[(1+\Phi){\rm d}\eta^2-(1-2\Psi){\rm d}x^i{\rm d}x^j]\; , $

        (7)

        where η is the conformal time with $ {\rm d}\eta=a^{-1}(t){\rm d}t $, while Φ and Ψ are gauge-invariant metric perturbations. The perturbed Einstein equation gives rise to the equations of Φ and Ψ [3134]:

        $ \begin{split} \nabla^2\Psi-3{\cal{H}}(\Psi'+{\cal{H}}\Phi)=\;&\ 4\pi G[\phi'(\delta\phi'-\phi'\Phi)\\ &+a^2V_{,\phi}\delta\phi-\sigma'(\delta\sigma'-\sigma'\Phi)]\; , \end{split} $

        (8)

        $ \begin{aligned} \Psi'+{\cal{H}}\Phi=&\ -4\pi G [-\phi'\delta\phi+\sigma'\delta\sigma]\; , \end{aligned} $

        (9)

        $ \begin{split} \Phi''+3{\cal{H}}\Phi'+(2{\cal{H}}'+{\cal{H}}^2)\Phi=\;&\ 4\pi G [\phi'(\delta\phi'-\phi'\Phi)\\ &-a^2V_{,\phi}\delta\phi-\sigma'(\delta\sigma'-\sigma'\Phi)]\; , \end{split} $

        (10)

        which can be combined to give the main equation of the Newtonian potential Φ:

        $\begin{split} &\Phi''+2\left({\cal{H}}-\frac{\phi''}{\phi'}\right)\Phi'+2\left({\cal{H}}'-{\cal{H}}\frac{\phi''}{\phi'}\right)\Phi-\nabla^2\Phi\\ &\quad =8\pi G\left(2{\cal{H}}+\frac{\phi''}{\phi'}\right)\sigma'\delta\sigma\; ,\end{split} $

        (11)

        where prime denotes derivative with respect to the conformal time η, and $ {\cal{H}} $ is the conformal Hubble parameter, $ {\cal{H}}\equiv a'/a $. Moreover, the curvature fluctuation ζ is related with Φ via:

        $ \zeta=\Phi+\frac{\cal{H}}{{\cal{H}}^2-{\cal{H}}'}(\Phi'+{\cal{H}}\Phi)\; , $

        (12)

        which is conserved on super-Hubble scales:

        $ (1+w)\zeta'=\frac{2k^2(\Phi'+{\cal{H}}\Phi)}{9{\cal{H}}^2}\; . $

        (13)

        In the following, we will see that different forms of potential $ V(\phi) $ will give rise to different backgrounds as well as different evolutions of the perturbations. In order to understand their main influence, we consider two broad categories of potentials, namely large field potentials and small field potentials.

      • I.   BOUNCE WITH TWO SCALAR FIELDS
        • The simplest Quintom model contains a quintessence-like scalar field and a phantom-like scalar field, with the Lagrangian [5, 14, 20] 3:

          $ S=\int {\rm d}^4x\sqrt{-g}\left[-\frac{1}{2}\nabla_\mu\phi\nabla^\mu\phi-V(\phi)+\frac{1}{2}\nabla_\mu\sigma\nabla^\mu\sigma-V(\sigma)\right]\; , $

          (2)

          where ϕ and σ stand for quintessence and phantom fields, respectively. The Friedmann equations are:

          $ \begin{aligned} H^2=&\ \frac{8\pi G}{3}\left[\frac{1}{2}\dot\phi^2-\frac{1}{2}\dot\sigma^2+V(\phi)+V(\sigma)\right]\; , \end{aligned} $

          (3)

          $ \begin{aligned} \dot H=&\ -4\pi G(\dot\phi^2-\dot\sigma^2)\; , \end{aligned} $

          (4)

          while the equations of motion of the two fields are:

          $ \begin{aligned} \ddot\phi+3H\dot\phi+V_{,\phi}= 0\; , \end{aligned} $

          (5)

          $ \begin{aligned} \ddot\sigma+3H\dot\sigma+V_{,\sigma}=&\ 0\; . \end{aligned} $

          (6)

          From the Friedmann equation one can see that, provided that the potentials of both fields are not less than zero, the only possibility for making $ H=0 $, as is required by bounce process, is to have the kinetic term of the σ-field become large and compensate the rest of the energy density. On the other hand, this term should not be dominant in order not to give rise to a negative energy density. This can be achieved by setting $ V(\sigma)=0 $, so that the equation of motion for the σ-field becomes $ \ddot\sigma+3H\dot\sigma=0 $ with the solution $ \dot\sigma^2\sim a^{-6} $. Therefore at the bounce point where the scale factor $ a(t) $ reaches its minimum, $ \dot\sigma^2 $ gets its maximum value, while whether the universe goes backwards or forwards, $ a(t) $ gets enlarged and $ \dot\sigma^2 $ will decay. Thus the evolution of the universe away from the bounce depends on the potential of the normal field ϕ only.

          To obtain the observational signals, it is important to analyze the perturbations generated by the model. The longitudinal metric perturbations are given by:

          $ {\rm d}s^2=a^2(\eta)[(1+\Phi){\rm d}\eta^2-(1-2\Psi){\rm d}x^i{\rm d}x^j]\; , $

          (7)

          where η is the conformal time with $ {\rm d}\eta=a^{-1}(t){\rm d}t $, while Φ and Ψ are gauge-invariant metric perturbations. The perturbed Einstein equation gives rise to the equations of Φ and Ψ [3134]:

          $ \begin{split} \nabla^2\Psi-3{\cal{H}}(\Psi'+{\cal{H}}\Phi)=\;&\ 4\pi G[\phi'(\delta\phi'-\phi'\Phi)\\ &+a^2V_{,\phi}\delta\phi-\sigma'(\delta\sigma'-\sigma'\Phi)]\; , \end{split} $

          (8)

          $ \begin{aligned} \Psi'+{\cal{H}}\Phi=&\ -4\pi G [-\phi'\delta\phi+\sigma'\delta\sigma]\; , \end{aligned} $

          (9)

          $ \begin{split} \Phi''+3{\cal{H}}\Phi'+(2{\cal{H}}'+{\cal{H}}^2)\Phi=\;&\ 4\pi G [\phi'(\delta\phi'-\phi'\Phi)\\ &-a^2V_{,\phi}\delta\phi-\sigma'(\delta\sigma'-\sigma'\Phi)]\; , \end{split} $

          (10)

          which can be combined to give the main equation of the Newtonian potential Φ:

          $\begin{split} &\Phi''+2\left({\cal{H}}-\frac{\phi''}{\phi'}\right)\Phi'+2\left({\cal{H}}'-{\cal{H}}\frac{\phi''}{\phi'}\right)\Phi-\nabla^2\Phi\\ &\quad =8\pi G\left(2{\cal{H}}+\frac{\phi''}{\phi'}\right)\sigma'\delta\sigma\; ,\end{split} $

          (11)

          where prime denotes derivative with respect to the conformal time η, and $ {\cal{H}} $ is the conformal Hubble parameter, $ {\cal{H}}\equiv a'/a $. Moreover, the curvature fluctuation ζ is related with Φ via:

          $ \zeta=\Phi+\frac{\cal{H}}{{\cal{H}}^2-{\cal{H}}'}(\Phi'+{\cal{H}}\Phi)\; , $

          (12)

          which is conserved on super-Hubble scales:

          $ (1+w)\zeta'=\frac{2k^2(\Phi'+{\cal{H}}\Phi)}{9{\cal{H}}^2}\; . $

          (13)

          In the following, we will see that different forms of potential $ V(\phi) $ will give rise to different backgrounds as well as different evolutions of the perturbations. In order to understand their main influence, we consider two broad categories of potentials, namely large field potentials and small field potentials.

        • A.   Large field potentials

        • The typical large field potential is the mass squared potential, $ V(\phi)=m^2\phi^2/2 $ [33]. This potential is symmetric with respect to its minimum at $ \phi=0 $, therefore it can give rise to symmetric evolution before and after the bounce.

          At the very beginning, the ϕ-field is set near the minimum and oscillates with increasing amplitude, due to the contraction of the Universe and the friction term. This will cause an oscillating behavior of the EoS, with the average value of $ \langle w\rangle=0 $. This phase is called the "heating phase". When the friction becomes less important, the Universe enters into a "slow-climbing" phase, where ϕ evolves slowly along the potential upwards, with the EoS $ w\simeq-1 $. Meanwhile, the kinetic term of the σ-field becomes more and more important. When it reaches the value of the ϕ-field’s energy density, the two parts cancel out and the bounce takes place. At the bounce point, the EoS reaches $ -\infty $, as has been demonstrated previously. After the bounce, the σ-field decays again, and the ϕ-field also rolls slowly along the potential to its minimum, and the EoS approaches -1 again. Finally, when the ϕ-field reaches the bottom, it oscillates around its minimum again with decreasing amplitude, with $ \langle w\rangle=0 $. The evolution behavior of w through the whole process is shown in Fig. 2.

          Figure 2.  (color online) Plot of evolution of EoS w in a double-field quintom bounce with large field potential. The figure is taken from Ref. [33].

          The evolution of the perturbation is given by the main equation (11). Before going into the details, we draw the evolution of perturbation wavelength scales vs. Hubble radius in Fig. 3. From this plot we can see that, the evolution is symmetric before and after the bounce, thus like the expanding phase, in the contracting phase, the Hubble radius will also shrink and expand (to infinity at the bounce point), and the perturbations will also exit and reenter the horizon. This will change the evolution of the perturbations and will affect the final observational signals in the CMB. Moreover, for large k (small scale) modes, the reentrance into the horizon takes place in the slow-climbing region, for small k (large scale) modes it happens in the bounce phase. Therefore the evolution of the perturbations is highly nontrivial.

          Figure 3.  (color online) The evolution of perturbation wavelengths with different comoving wavenumbers k as well as the Hubble radius $ |H^{-1}| $. The figure is taken from Ref. [33].

          In the heating phase, since the EoS is effectively zero, one has

          $ a\propto \eta^2\; ,\; \; \; {\cal{H}}=\frac{2}{\eta}\; ,\; \; \; \phi'=\frac{1}{\eta}\; . $

          (14)

          Moreover, since in this phase the kinetic term of σ-field is negligible, we can omit the right hand side of Eq. (11) for simplicity. Then the equation becomes (in momentum space)

          $ \Phi''_k+\frac{6}{\eta}\Phi'_k+k^2\Phi_k=0\; , $

          (15)

          with the solution

          $ \Phi_k=\eta^{-5/2}[k^{-5/2}A_kJ_{5/2}(k\eta)+k^{5/2}B_kJ_{-5/2}(k\eta)]\; . $

          (16)

          The coefficients can be determined by matching the sub-horizon approximation ($ |k\eta|\gg 1 $) with the initial condition of $ \Phi_k $, which we take to be the Bunch-Davies vacuum solution:

          $ \Phi_k=A\eta^{-3}\frac{{\rm e}^{-{\rm i}k\eta}}{\sqrt{2k^3}}\; ,\; \; \; A\equiv4\pi G\sqrt{\rho_0}\eta_0^3\; , $

          (17)

          from which we get $ A_k=({\rm i}\sqrt{\pi}/2)Ak^{3/2} $, $ B_k= -(\sqrt{\pi}/2) \times Ak^{-7/2} $. On the other hand, in the super-horizon approximation with $ |k\eta|\ll 1 $, the approximate solution becomes:

          $ \Phi_k\simeq \frac{{\rm i}Ak^{3/2}}{15\sqrt{2}}-\frac{3A}{\sqrt{2}}k^{-7/2}\eta^{-5}\; , $

          (18)

          which has one constant mode and one growing mode.

          We assume that the Universe enters into the slow-climbing phase at the time $ \eta_c $. After that, Eq. (11) becomes

          $ \Phi''_k+(2\epsilon_H-\delta_H){\cal{H}}\Phi'_k-\delta_H{\cal{H}}^2\Phi_k+k^2\Phi_k=0\; , $

          (19)

          where we defined the slow-climb parameters: $ \epsilon_H\equiv -\dot H/H^2 $, $ \delta_H\equiv\dot\epsilon_H/(H\epsilon_H) $, satisfying $ |\epsilon_H|,|\delta_H|\ll 1 $. The solution is

          $ \Phi_k=(\eta-\tilde{\eta}_c)^\alpha[k^{-\nu}C_kJ_\nu(k(\eta-\tilde{\eta}_c))+k^\nu D_kJ_{-\nu}(k(\eta-\tilde{\eta}_c))]\; , $

          (20)

          where $ \alpha\simeq 1/2 $, $ \nu\simeq 1/2 $, and $ \tilde{\eta}_c\equiv \eta_c+1/{\cal{H}}_c $. The super-horizon approximation, which is to be connected with the super-horizon approximation in the heating phase (18), turns out to be:

          $ \Phi_k\simeq \sqrt{\frac{2}{\pi}}[C_k(\eta-\tilde{\eta}_c)+D_k]\; . $

          (21)

          By matching with (18) using matching conditions [35, 36] (see also [37]), one gets the coefficients:

          $ \begin{aligned} C_k=&\ -{\cal{H}}_c\left[\frac{1}{15}\left(1-\frac{2}{3}\epsilon_H\right)A_k+3(1+\epsilon_H)B_k\eta_c^{-5}\right]\; , \end{aligned} $

          (22)

          $ \begin{aligned} D_k=&\ \epsilon_H\left(\frac{2}{45}A_k-3B_k\eta_c^{-5}\right)\simeq 0\; . \end{aligned} $

          (23)

          Thus the main contribution comes from the growing mode. On the other hand, the sub-horizon approximation of (20) is

          $ \Phi_k\simeq\sqrt{\frac{2}{\pi}}\left[\frac{C_k}{k}\sin(k(\eta-\tilde{\eta}_c))+D_k\cos(k(\eta-\tilde{\eta}_c))\right]\; . $

          (24)

          When the universe enters into a bounce phase, the EoS goes down towards $ -\infty $ and goes up again to above -1, so it will be complicated to analyze the dynamics of the perturbations. Nevertheless, it is convenient to parametrize the Hubble parameter as a linear function which crosses zero at the bounce point, namely

          $ {\cal{H}}\simeq \frac{y}{2} (\eta-\eta_B)\; ,\; \; \; a\simeq \frac{a_B}{1-y(\eta-\eta_B)^2/4}\; , $

          (25)

          where $ \eta_B $ is the time when bounce occurs, and $ a_B $ is the scale factor at $ \eta_B $. Moreover, during the bounce phase $ \dot\sigma $ becomes important. However, from the background equation (4) one approximately has $ \phi''\simeq-2{\cal{H}}\phi' $, so the right hand side of Eq. (11) can still be neglected. This can be easily seen considering calculations around Eq. (4). Differentiating Eq. (4) with respect to t we can get

          $ \ddot H=-8\pi G(\dot\phi\ddot\phi-\dot\sigma\ddot\sigma)\; . $

          (26)

          As the σ field has no potential, we simply get $ \ddot\sigma+3H\dot\sigma=0 $. Combining Eqs. (4) and (26) and after some calculations one finds

          $ \frac{\ddot\phi}{\dot\phi}=-3H\times\frac{3H\dot\sigma^2+\ddot H/(8\pi G)}{3H\dot\sigma^2-3H\dot H/(4\pi G)}\simeq -3H\times\frac{2\sigma'^2-y/(6\pi G)}{2\sigma'^2-y/(4\pi G)}\; , $

          (27)

          where in the last step we made use of the parametrization (25). It implies that as long as $ \sigma'^2\gg y/G $ is satisfied, one has $ \ddot\phi/\dot\phi\simeq -3H $, which is equivalent to $ \phi''+2{\cal{H}}\phi'\simeq 0 $. We verify this result by numerical calculation shown in Fig. 4. Therefore, the perturbation equation during the bounce phase becomes:

          Figure 4.  The evolution of the factor $ 3H+\ddot\phi/\dot\phi $ with respect to cosmic time t in the bouncing phase. The figure is taken from Ref. [33].

          $ \Phi''_k+3y(\eta-\eta_B)\Phi'_k+(k^2+y)\Phi_k\simeq 0\; , $

          (28)

          and the solution is

          $ \Phi_k\simeq {\rm e}^{-\frac{3}{4}y(\eta-\eta_B)^2}\left[E_k\sin[k(\eta-\eta_B)]+F_k\cos[k(\eta-\eta_B)]\right]\; . $

          (29)

          The coefficients in solution (29), namely $ E_k $ and $ F_k $, are obtained by matching (29) to the slow-climbing solution (20) with the same matching conditions. From Fig. 3, one can see that there are two cases, namely matching in the sub-horizon region (for large k) and matching in the super-horizon region (for small k). In the first case (29) matches (24), while in the second case it matches (21). This will give two sets of $ E_k $ and $ F_k $, and the detailed expressions are given in [33]. From the solution we can also see that, while in the first case $ \Phi_k $ in the bounce phase is dominated by the growing mode in slow-climbing phase, in the second case it is dominated by the constant mode.

          After the bounce, the Universe enters into an expanding phase at the time $ \eta_{B+} $, where the ϕ field is slow-rolling. In this phase, the perturbation equation is the same as Eq. (19), with the solution:

          $ \Phi_k=(\eta-\tilde{\eta}_{B+})^\alpha[k^{-\nu}G_kJ_\nu(k(\eta-\tilde{\eta}_{B+}))+k^\nu H_kJ_{-\nu}(k(\eta-\tilde{\eta}_{B+}))]\; , $

          (30)

          where $ \tilde{\eta}_{B+}\equiv\eta_{B+}+1/{\cal{H}}_{B+} $. The sub-horizon approximation is:

          $ \Phi_k\simeq\sqrt{\frac{2}{\pi}}\left[\frac{G_k}{k}\sin(k(\eta-\tilde{\eta}_{B+}))+H_k\cos(k(\eta-\tilde{\eta}_{B+}))\right]\; , $

          (31)

          while the super-horizon approximation is:

          $ \Phi_k\simeq \sqrt{\frac{2}{\pi}}[G_k(\eta-\tilde{\eta}_{B+})+H_k]\; , $

          (32)

          where the first and second terms correspond to the decaying and constant mode, respectively. The coefficients $ G_k $ and $ H_k $ will also be obtained by matching the solution (30) to the solution (29) in the bounce phase. Similarly, there are two cases of matching at sub-horizon and super horizon regions, where (29) matches to (31) and (32), respectively. The detailed expressions are given in [33]. In the first case where both $ G_k $ and $ H_k $ are important, the final perturbation $ \Phi_k $ is dominated by its constant mode. While in the second case where $ H_k\simeq 0 $ is obtained, $ \Phi_k $ is dominated by its decaying mode [22, 3850].

          CMB observations suggest a nearly scale-invariant power spectrum of primordial curvature perturbations, which could be realized by either de-Sitter like expansion (inflation) or matter-like contraction. However, due to the addition of the slow-climbing phase, the k-dependence of the perturbations become more complicated, and in general cannot give rise to a scale-invariant power spectrum (see numerical results in [33]).

        • A.   Large field potentials

        • The typical large field potential is the mass squared potential, $ V(\phi)=m^2\phi^2/2 $ [33]. This potential is symmetric with respect to its minimum at $ \phi=0 $, therefore it can give rise to symmetric evolution before and after the bounce.

          At the very beginning, the ϕ-field is set near the minimum and oscillates with increasing amplitude, due to the contraction of the Universe and the friction term. This will cause an oscillating behavior of the EoS, with the average value of $ \langle w\rangle=0 $. This phase is called the "heating phase". When the friction becomes less important, the Universe enters into a "slow-climbing" phase, where ϕ evolves slowly along the potential upwards, with the EoS $ w\simeq-1 $. Meanwhile, the kinetic term of the σ-field becomes more and more important. When it reaches the value of the ϕ-field’s energy density, the two parts cancel out and the bounce takes place. At the bounce point, the EoS reaches $ -\infty $, as has been demonstrated previously. After the bounce, the σ-field decays again, and the ϕ-field also rolls slowly along the potential to its minimum, and the EoS approaches -1 again. Finally, when the ϕ-field reaches the bottom, it oscillates around its minimum again with decreasing amplitude, with $ \langle w\rangle=0 $. The evolution behavior of w through the whole process is shown in Fig. 2.

          Figure 2.  (color online) Plot of evolution of EoS w in a double-field quintom bounce with large field potential. The figure is taken from Ref. [33].

          The evolution of the perturbation is given by the main equation (11). Before going into the details, we draw the evolution of perturbation wavelength scales vs. Hubble radius in Fig. 3. From this plot we can see that, the evolution is symmetric before and after the bounce, thus like the expanding phase, in the contracting phase, the Hubble radius will also shrink and expand (to infinity at the bounce point), and the perturbations will also exit and reenter the horizon. This will change the evolution of the perturbations and will affect the final observational signals in the CMB. Moreover, for large k (small scale) modes, the reentrance into the horizon takes place in the slow-climbing region, for small k (large scale) modes it happens in the bounce phase. Therefore the evolution of the perturbations is highly nontrivial.

          Figure 3.  (color online) The evolution of perturbation wavelengths with different comoving wavenumbers k as well as the Hubble radius $ |H^{-1}| $. The figure is taken from Ref. [33].

          In the heating phase, since the EoS is effectively zero, one has

          $ a\propto \eta^2\; ,\; \; \; {\cal{H}}=\frac{2}{\eta}\; ,\; \; \; \phi'=\frac{1}{\eta}\; . $

          (14)

          Moreover, since in this phase the kinetic term of σ-field is negligible, we can omit the right hand side of Eq. (11) for simplicity. Then the equation becomes (in momentum space)

          $ \Phi''_k+\frac{6}{\eta}\Phi'_k+k^2\Phi_k=0\; , $

          (15)

          with the solution

          $ \Phi_k=\eta^{-5/2}[k^{-5/2}A_kJ_{5/2}(k\eta)+k^{5/2}B_kJ_{-5/2}(k\eta)]\; . $

          (16)

          The coefficients can be determined by matching the sub-horizon approximation ($ |k\eta|\gg 1 $) with the initial condition of $ \Phi_k $, which we take to be the Bunch-Davies vacuum solution:

          $ \Phi_k=A\eta^{-3}\frac{{\rm e}^{-{\rm i}k\eta}}{\sqrt{2k^3}}\; ,\; \; \; A\equiv4\pi G\sqrt{\rho_0}\eta_0^3\; , $

          (17)

          from which we get $ A_k=({\rm i}\sqrt{\pi}/2)Ak^{3/2} $, $ B_k= -(\sqrt{\pi}/2) \times Ak^{-7/2} $. On the other hand, in the super-horizon approximation with $ |k\eta|\ll 1 $, the approximate solution becomes:

          $ \Phi_k\simeq \frac{{\rm i}Ak^{3/2}}{15\sqrt{2}}-\frac{3A}{\sqrt{2}}k^{-7/2}\eta^{-5}\; , $

          (18)

          which has one constant mode and one growing mode.

          We assume that the Universe enters into the slow-climbing phase at the time $ \eta_c $. After that, Eq. (11) becomes

          $ \Phi''_k+(2\epsilon_H-\delta_H){\cal{H}}\Phi'_k-\delta_H{\cal{H}}^2\Phi_k+k^2\Phi_k=0\; , $

          (19)

          where we defined the slow-climb parameters: $ \epsilon_H\equiv -\dot H/H^2 $, $ \delta_H\equiv\dot\epsilon_H/(H\epsilon_H) $, satisfying $ |\epsilon_H|,|\delta_H|\ll 1 $. The solution is

          $ \Phi_k=(\eta-\tilde{\eta}_c)^\alpha[k^{-\nu}C_kJ_\nu(k(\eta-\tilde{\eta}_c))+k^\nu D_kJ_{-\nu}(k(\eta-\tilde{\eta}_c))]\; , $

          (20)

          where $ \alpha\simeq 1/2 $, $ \nu\simeq 1/2 $, and $ \tilde{\eta}_c\equiv \eta_c+1/{\cal{H}}_c $. The super-horizon approximation, which is to be connected with the super-horizon approximation in the heating phase (18), turns out to be:

          $ \Phi_k\simeq \sqrt{\frac{2}{\pi}}[C_k(\eta-\tilde{\eta}_c)+D_k]\; . $

          (21)

          By matching with (18) using matching conditions [35, 36] (see also [37]), one gets the coefficients:

          $ \begin{aligned} C_k=&\ -{\cal{H}}_c\left[\frac{1}{15}\left(1-\frac{2}{3}\epsilon_H\right)A_k+3(1+\epsilon_H)B_k\eta_c^{-5}\right]\; , \end{aligned} $

          (22)

          $ \begin{aligned} D_k=&\ \epsilon_H\left(\frac{2}{45}A_k-3B_k\eta_c^{-5}\right)\simeq 0\; . \end{aligned} $

          (23)

          Thus the main contribution comes from the growing mode. On the other hand, the sub-horizon approximation of (20) is

          $ \Phi_k\simeq\sqrt{\frac{2}{\pi}}\left[\frac{C_k}{k}\sin(k(\eta-\tilde{\eta}_c))+D_k\cos(k(\eta-\tilde{\eta}_c))\right]\; . $

          (24)

          When the universe enters into a bounce phase, the EoS goes down towards $ -\infty $ and goes up again to above -1, so it will be complicated to analyze the dynamics of the perturbations. Nevertheless, it is convenient to parametrize the Hubble parameter as a linear function which crosses zero at the bounce point, namely

          $ {\cal{H}}\simeq \frac{y}{2} (\eta-\eta_B)\; ,\; \; \; a\simeq \frac{a_B}{1-y(\eta-\eta_B)^2/4}\; , $

          (25)

          where $ \eta_B $ is the time when bounce occurs, and $ a_B $ is the scale factor at $ \eta_B $. Moreover, during the bounce phase $ \dot\sigma $ becomes important. However, from the background equation (4) one approximately has $ \phi''\simeq-2{\cal{H}}\phi' $, so the right hand side of Eq. (11) can still be neglected. This can be easily seen considering calculations around Eq. (4). Differentiating Eq. (4) with respect to t we can get

          $ \ddot H=-8\pi G(\dot\phi\ddot\phi-\dot\sigma\ddot\sigma)\; . $

          (26)

          As the σ field has no potential, we simply get $ \ddot\sigma+3H\dot\sigma=0 $. Combining Eqs. (4) and (26) and after some calculations one finds

          $ \frac{\ddot\phi}{\dot\phi}=-3H\times\frac{3H\dot\sigma^2+\ddot H/(8\pi G)}{3H\dot\sigma^2-3H\dot H/(4\pi G)}\simeq -3H\times\frac{2\sigma'^2-y/(6\pi G)}{2\sigma'^2-y/(4\pi G)}\; , $

          (27)

          where in the last step we made use of the parametrization (25). It implies that as long as $ \sigma'^2\gg y/G $ is satisfied, one has $ \ddot\phi/\dot\phi\simeq -3H $, which is equivalent to $ \phi''+2{\cal{H}}\phi'\simeq 0 $. We verify this result by numerical calculation shown in Fig. 4. Therefore, the perturbation equation during the bounce phase becomes:

          Figure 4.  The evolution of the factor $ 3H+\ddot\phi/\dot\phi $ with respect to cosmic time t in the bouncing phase. The figure is taken from Ref. [33].

          $ \Phi''_k+3y(\eta-\eta_B)\Phi'_k+(k^2+y)\Phi_k\simeq 0\; , $

          (28)

          and the solution is

          $ \Phi_k\simeq {\rm e}^{-\frac{3}{4}y(\eta-\eta_B)^2}\left[E_k\sin[k(\eta-\eta_B)]+F_k\cos[k(\eta-\eta_B)]\right]\; . $

          (29)

          The coefficients in solution (29), namely $ E_k $ and $ F_k $, are obtained by matching (29) to the slow-climbing solution (20) with the same matching conditions. From Fig. 3, one can see that there are two cases, namely matching in the sub-horizon region (for large k) and matching in the super-horizon region (for small k). In the first case (29) matches (24), while in the second case it matches (21). This will give two sets of $ E_k $ and $ F_k $, and the detailed expressions are given in [33]. From the solution we can also see that, while in the first case $ \Phi_k $ in the bounce phase is dominated by the growing mode in slow-climbing phase, in the second case it is dominated by the constant mode.

          After the bounce, the Universe enters into an expanding phase at the time $ \eta_{B+} $, where the ϕ field is slow-rolling. In this phase, the perturbation equation is the same as Eq. (19), with the solution:

          $ \Phi_k=(\eta-\tilde{\eta}_{B+})^\alpha[k^{-\nu}G_kJ_\nu(k(\eta-\tilde{\eta}_{B+}))+k^\nu H_kJ_{-\nu}(k(\eta-\tilde{\eta}_{B+}))]\; , $

          (30)

          where $ \tilde{\eta}_{B+}\equiv\eta_{B+}+1/{\cal{H}}_{B+} $. The sub-horizon approximation is:

          $ \Phi_k\simeq\sqrt{\frac{2}{\pi}}\left[\frac{G_k}{k}\sin(k(\eta-\tilde{\eta}_{B+}))+H_k\cos(k(\eta-\tilde{\eta}_{B+}))\right]\; , $

          (31)

          while the super-horizon approximation is:

          $ \Phi_k\simeq \sqrt{\frac{2}{\pi}}[G_k(\eta-\tilde{\eta}_{B+})+H_k]\; , $

          (32)

          where the first and second terms correspond to the decaying and constant mode, respectively. The coefficients $ G_k $ and $ H_k $ will also be obtained by matching the solution (30) to the solution (29) in the bounce phase. Similarly, there are two cases of matching at sub-horizon and super horizon regions, where (29) matches to (31) and (32), respectively. The detailed expressions are given in [33]. In the first case where both $ G_k $ and $ H_k $ are important, the final perturbation $ \Phi_k $ is dominated by its constant mode. While in the second case where $ H_k\simeq 0 $ is obtained, $ \Phi_k $ is dominated by its decaying mode [22, 3850].

          CMB observations suggest a nearly scale-invariant power spectrum of primordial curvature perturbations, which could be realized by either de-Sitter like expansion (inflation) or matter-like contraction. However, due to the addition of the slow-climbing phase, the k-dependence of the perturbations become more complicated, and in general cannot give rise to a scale-invariant power spectrum (see numerical results in [33]).

        • A.   Large field potentials

        • The typical large field potential is the mass squared potential, $ V(\phi)=m^2\phi^2/2 $ [33]. This potential is symmetric with respect to its minimum at $ \phi=0 $, therefore it can give rise to symmetric evolution before and after the bounce.

          At the very beginning, the ϕ-field is set near the minimum and oscillates with increasing amplitude, due to the contraction of the Universe and the friction term. This will cause an oscillating behavior of the EoS, with the average value of $ \langle w\rangle=0 $. This phase is called the "heating phase". When the friction becomes less important, the Universe enters into a "slow-climbing" phase, where ϕ evolves slowly along the potential upwards, with the EoS $ w\simeq-1 $. Meanwhile, the kinetic term of the σ-field becomes more and more important. When it reaches the value of the ϕ-field’s energy density, the two parts cancel out and the bounce takes place. At the bounce point, the EoS reaches $ -\infty $, as has been demonstrated previously. After the bounce, the σ-field decays again, and the ϕ-field also rolls slowly along the potential to its minimum, and the EoS approaches -1 again. Finally, when the ϕ-field reaches the bottom, it oscillates around its minimum again with decreasing amplitude, with $ \langle w\rangle=0 $. The evolution behavior of w through the whole process is shown in Fig. 2.

          Figure 2.  (color online) Plot of evolution of EoS w in a double-field quintom bounce with large field potential. The figure is taken from Ref. [33].

          The evolution of the perturbation is given by the main equation (11). Before going into the details, we draw the evolution of perturbation wavelength scales vs. Hubble radius in Fig. 3. From this plot we can see that, the evolution is symmetric before and after the bounce, thus like the expanding phase, in the contracting phase, the Hubble radius will also shrink and expand (to infinity at the bounce point), and the perturbations will also exit and reenter the horizon. This will change the evolution of the perturbations and will affect the final observational signals in the CMB. Moreover, for large k (small scale) modes, the reentrance into the horizon takes place in the slow-climbing region, for small k (large scale) modes it happens in the bounce phase. Therefore the evolution of the perturbations is highly nontrivial.

          Figure 3.  (color online) The evolution of perturbation wavelengths with different comoving wavenumbers k as well as the Hubble radius $ |H^{-1}| $. The figure is taken from Ref. [33].

          In the heating phase, since the EoS is effectively zero, one has

          $ a\propto \eta^2\; ,\; \; \; {\cal{H}}=\frac{2}{\eta}\; ,\; \; \; \phi'=\frac{1}{\eta}\; . $

          (14)

          Moreover, since in this phase the kinetic term of σ-field is negligible, we can omit the right hand side of Eq. (11) for simplicity. Then the equation becomes (in momentum space)

          $ \Phi''_k+\frac{6}{\eta}\Phi'_k+k^2\Phi_k=0\; , $

          (15)

          with the solution

          $ \Phi_k=\eta^{-5/2}[k^{-5/2}A_kJ_{5/2}(k\eta)+k^{5/2}B_kJ_{-5/2}(k\eta)]\; . $

          (16)

          The coefficients can be determined by matching the sub-horizon approximation ($ |k\eta|\gg 1 $) with the initial condition of $ \Phi_k $, which we take to be the Bunch-Davies vacuum solution:

          $ \Phi_k=A\eta^{-3}\frac{{\rm e}^{-{\rm i}k\eta}}{\sqrt{2k^3}}\; ,\; \; \; A\equiv4\pi G\sqrt{\rho_0}\eta_0^3\; , $

          (17)

          from which we get $ A_k=({\rm i}\sqrt{\pi}/2)Ak^{3/2} $, $ B_k= -(\sqrt{\pi}/2) \times Ak^{-7/2} $. On the other hand, in the super-horizon approximation with $ |k\eta|\ll 1 $, the approximate solution becomes:

          $ \Phi_k\simeq \frac{{\rm i}Ak^{3/2}}{15\sqrt{2}}-\frac{3A}{\sqrt{2}}k^{-7/2}\eta^{-5}\; , $

          (18)

          which has one constant mode and one growing mode.

          We assume that the Universe enters into the slow-climbing phase at the time $ \eta_c $. After that, Eq. (11) becomes

          $ \Phi''_k+(2\epsilon_H-\delta_H){\cal{H}}\Phi'_k-\delta_H{\cal{H}}^2\Phi_k+k^2\Phi_k=0\; , $

          (19)

          where we defined the slow-climb parameters: $ \epsilon_H\equiv -\dot H/H^2 $, $ \delta_H\equiv\dot\epsilon_H/(H\epsilon_H) $, satisfying $ |\epsilon_H|,|\delta_H|\ll 1 $. The solution is

          $ \Phi_k=(\eta-\tilde{\eta}_c)^\alpha[k^{-\nu}C_kJ_\nu(k(\eta-\tilde{\eta}_c))+k^\nu D_kJ_{-\nu}(k(\eta-\tilde{\eta}_c))]\; , $

          (20)

          where $ \alpha\simeq 1/2 $, $ \nu\simeq 1/2 $, and $ \tilde{\eta}_c\equiv \eta_c+1/{\cal{H}}_c $. The super-horizon approximation, which is to be connected with the super-horizon approximation in the heating phase (18), turns out to be:

          $ \Phi_k\simeq \sqrt{\frac{2}{\pi}}[C_k(\eta-\tilde{\eta}_c)+D_k]\; . $

          (21)

          By matching with (18) using matching conditions [35, 36] (see also [37]), one gets the coefficients:

          $ \begin{aligned} C_k=&\ -{\cal{H}}_c\left[\frac{1}{15}\left(1-\frac{2}{3}\epsilon_H\right)A_k+3(1+\epsilon_H)B_k\eta_c^{-5}\right]\; , \end{aligned} $

          (22)

          $ \begin{aligned} D_k=&\ \epsilon_H\left(\frac{2}{45}A_k-3B_k\eta_c^{-5}\right)\simeq 0\; . \end{aligned} $

          (23)

          Thus the main contribution comes from the growing mode. On the other hand, the sub-horizon approximation of (20) is

          $ \Phi_k\simeq\sqrt{\frac{2}{\pi}}\left[\frac{C_k}{k}\sin(k(\eta-\tilde{\eta}_c))+D_k\cos(k(\eta-\tilde{\eta}_c))\right]\; . $

          (24)

          When the universe enters into a bounce phase, the EoS goes down towards $ -\infty $ and goes up again to above -1, so it will be complicated to analyze the dynamics of the perturbations. Nevertheless, it is convenient to parametrize the Hubble parameter as a linear function which crosses zero at the bounce point, namely

          $ {\cal{H}}\simeq \frac{y}{2} (\eta-\eta_B)\; ,\; \; \; a\simeq \frac{a_B}{1-y(\eta-\eta_B)^2/4}\; , $

          (25)

          where $ \eta_B $ is the time when bounce occurs, and $ a_B $ is the scale factor at $ \eta_B $. Moreover, during the bounce phase $ \dot\sigma $ becomes important. However, from the background equation (4) one approximately has $ \phi''\simeq-2{\cal{H}}\phi' $, so the right hand side of Eq. (11) can still be neglected. This can be easily seen considering calculations around Eq. (4). Differentiating Eq. (4) with respect to t we can get

          $ \ddot H=-8\pi G(\dot\phi\ddot\phi-\dot\sigma\ddot\sigma)\; . $

          (26)

          As the σ field has no potential, we simply get $ \ddot\sigma+3H\dot\sigma=0 $. Combining Eqs. (4) and (26) and after some calculations one finds

          $ \frac{\ddot\phi}{\dot\phi}=-3H\times\frac{3H\dot\sigma^2+\ddot H/(8\pi G)}{3H\dot\sigma^2-3H\dot H/(4\pi G)}\simeq -3H\times\frac{2\sigma'^2-y/(6\pi G)}{2\sigma'^2-y/(4\pi G)}\; , $

          (27)

          where in the last step we made use of the parametrization (25). It implies that as long as $ \sigma'^2\gg y/G $ is satisfied, one has $ \ddot\phi/\dot\phi\simeq -3H $, which is equivalent to $ \phi''+2{\cal{H}}\phi'\simeq 0 $. We verify this result by numerical calculation shown in Fig. 4. Therefore, the perturbation equation during the bounce phase becomes:

          Figure 4.  The evolution of the factor $ 3H+\ddot\phi/\dot\phi $ with respect to cosmic time t in the bouncing phase. The figure is taken from Ref. [33].

          $ \Phi''_k+3y(\eta-\eta_B)\Phi'_k+(k^2+y)\Phi_k\simeq 0\; , $

          (28)

          and the solution is

          $ \Phi_k\simeq {\rm e}^{-\frac{3}{4}y(\eta-\eta_B)^2}\left[E_k\sin[k(\eta-\eta_B)]+F_k\cos[k(\eta-\eta_B)]\right]\; . $

          (29)

          The coefficients in solution (29), namely $ E_k $ and $ F_k $, are obtained by matching (29) to the slow-climbing solution (20) with the same matching conditions. From Fig. 3, one can see that there are two cases, namely matching in the sub-horizon region (for large k) and matching in the super-horizon region (for small k). In the first case (29) matches (24), while in the second case it matches (21). This will give two sets of $ E_k $ and $ F_k $, and the detailed expressions are given in [33]. From the solution we can also see that, while in the first case $ \Phi_k $ in the bounce phase is dominated by the growing mode in slow-climbing phase, in the second case it is dominated by the constant mode.

          After the bounce, the Universe enters into an expanding phase at the time $ \eta_{B+} $, where the ϕ field is slow-rolling. In this phase, the perturbation equation is the same as Eq. (19), with the solution:

          $ \Phi_k=(\eta-\tilde{\eta}_{B+})^\alpha[k^{-\nu}G_kJ_\nu(k(\eta-\tilde{\eta}_{B+}))+k^\nu H_kJ_{-\nu}(k(\eta-\tilde{\eta}_{B+}))]\; , $

          (30)

          where $ \tilde{\eta}_{B+}\equiv\eta_{B+}+1/{\cal{H}}_{B+} $. The sub-horizon approximation is:

          $ \Phi_k\simeq\sqrt{\frac{2}{\pi}}\left[\frac{G_k}{k}\sin(k(\eta-\tilde{\eta}_{B+}))+H_k\cos(k(\eta-\tilde{\eta}_{B+}))\right]\; , $

          (31)

          while the super-horizon approximation is:

          $ \Phi_k\simeq \sqrt{\frac{2}{\pi}}[G_k(\eta-\tilde{\eta}_{B+})+H_k]\; , $

          (32)

          where the first and second terms correspond to the decaying and constant mode, respectively. The coefficients $ G_k $ and $ H_k $ will also be obtained by matching the solution (30) to the solution (29) in the bounce phase. Similarly, there are two cases of matching at sub-horizon and super horizon regions, where (29) matches to (31) and (32), respectively. The detailed expressions are given in [33]. In the first case where both $ G_k $ and $ H_k $ are important, the final perturbation $ \Phi_k $ is dominated by its constant mode. While in the second case where $ H_k\simeq 0 $ is obtained, $ \Phi_k $ is dominated by its decaying mode [22, 3850].

          CMB observations suggest a nearly scale-invariant power spectrum of primordial curvature perturbations, which could be realized by either de-Sitter like expansion (inflation) or matter-like contraction. However, due to the addition of the slow-climbing phase, the k-dependence of the perturbations become more complicated, and in general cannot give rise to a scale-invariant power spectrum (see numerical results in [33]).

        • B.   Small field potential

        • The typical small field potential is the Coleman-Weinberg potential, $ V(\phi)=(\lambda\phi^4/4)\ln(|\phi|/v)-\lambda(\phi^4- v^4)/16 $ [34, 51]. As a result of this potential, the vacuum is shifted to $ \phi=\pm v $. Therefore, the symmetry is spontaneously broken in the vacuum. It will also give rise to an asymmetric bounce behavior.

          Figure 5 shows the evolution behavior of w over the process in this model. Initially, the ϕ field is still set around the minimum value, and it can oscillate with its amplitude increasing and $ \langle w\rangle=0 $. When the field reaches the plateau, the σ field catches up, and the bounce takes place when w goes to minus infinity. After the bounce, when the ϕ-field rolls along the plateau, the Universe enters into a slow-roll phase where $ w\simeq-1 $, very much like inflation. Finally when the field drops into another minimum and oscillates, the inflation ends and the field becomes oscillating again, with $ \langle w\rangle=0 $.

          Figure 5.  (color online) Plot of evolution of EoS w in a double-field quintom bounce with small field potential. The figure is taken from Ref. [34].

          We also draw the evolution plot of perturbation wavelength scales v.s. Hubble radius in Fig. 6. From the plot we can see that, since now the symmetry before and after the bounce is broken, the evolution of the perturbations will also be different. While the small k modes of perturbations still need to exit and reenter the horizon in the contracting phase, the large k modes needn't, and will stay inside the horizon until the expanding phase. In the following, we will focus on these modes.

          Figure 6.  (color online) The evolution of perturbation wavelengths with different comoving wavenumbers k as well as Hubble radius $ |H^{-1}| $. The figure is taken from Ref. [34].

          In the contracting phase where the EoS of the Universe oscillates continuously with $ \langle w\rangle=0 $, the background and perturbation evolutions are still described by Eqs. (14) and (15). While the solution is given in Eq. (16), since we are only interested in the subhorizon solution, we only have:

          $ \Phi_k=4\pi G\frac{\sqrt{\rho_0}\eta_0^3}{\eta^3}\frac{{\rm e}^{-{\rm i}k\eta}}{\sqrt{2k^3}}\; . $

          (33)

          Since in this case the evolution is not symmetric and for large k modes there is no horizon-exit behavior before bounce, we don't have a "slow-climbing" phase. In bounce phase, the equation is the same as (28), and the solution is given by (29). After the bounce, the universe will enter a slow-rolling expanding phase, in which the perturbation equation is (19) and the solution is the same as (30).

          The coefficients $ E_k $, $ F_k $, $ G_k $ and $ H_k $ are obtained using the same matching method as above. However, since here the matchings are only performed at sub-horizon regions, the final perturbations are determined by the constant mode of $ H_k $. From the tedious calculation done in [34], one obtains

          $ H_k=\sqrt{\frac{\pi}{2}}\frac{4\pi G}{\sqrt{2k^3}}|\dot\phi|{\rm e}^{-{\rm i}k\tilde{\eta}_{B+}}\left\{1+\frac{3{\rm e}^{-{\rm i}k(\eta_{B-}-\tilde{\eta}_{B+})}}{k\eta_{B-}}\sin[k(\eta_{B-}-\tilde{\eta}_{B+})]\right\}\; . $

          (34)

          where $ \eta_{B-} $ is the time when the Universe enters the bounce phase. Moreover, in the super-horizon region with nearly de-Sitter spacetime we have $ \zeta\simeq \Phi/\epsilon_H $, Thus one obtains the power spectrum of the curvature perturbations

          $ P_\zeta\equiv\frac{k^3}{2\pi^2}|\zeta|^2=\frac{8G^2\rho}{3\epsilon_H}\left\{1-\frac{3{\cal{H}}_{B-}}{2k}\sin\frac{2k}{{\cal{H}}_{B+}}\right\}\; . $

          (35)

          This result indicates that for large k (small scale) region where the last term could be neglected, the power spectrum is both constant and scale-invariant, while for small k (large scale) region where the second term becomes important, the spectrum gets some wiggle-like corrections. This is due to the effect of perturbations during bounce region. Moreover, for even larger scales, as the perturbations exit the horizon before the bounce, the power spectrum will be blue-tilted and will thus get suppressed [5254]. There are already some hints of these signals in the observations [5558]. If they are furtherly confirmed, they will provide a smoking gun for a bouncing cosmology.

        • B.   Small field potential

        • The typical small field potential is the Coleman-Weinberg potential, $ V(\phi)=(\lambda\phi^4/4)\ln(|\phi|/v)-\lambda(\phi^4- v^4)/16 $ [34, 51]. As a result of this potential, the vacuum is shifted to $ \phi=\pm v $. Therefore, the symmetry is spontaneously broken in the vacuum. It will also give rise to an asymmetric bounce behavior.

          Figure 5 shows the evolution behavior of w over the process in this model. Initially, the ϕ field is still set around the minimum value, and it can oscillate with its amplitude increasing and $ \langle w\rangle=0 $. When the field reaches the plateau, the σ field catches up, and the bounce takes place when w goes to minus infinity. After the bounce, when the ϕ-field rolls along the plateau, the Universe enters into a slow-roll phase where $ w\simeq-1 $, very much like inflation. Finally when the field drops into another minimum and oscillates, the inflation ends and the field becomes oscillating again, with $ \langle w\rangle=0 $.

          Figure 5.  (color online) Plot of evolution of EoS w in a double-field quintom bounce with small field potential. The figure is taken from Ref. [34].

          We also draw the evolution plot of perturbation wavelength scales v.s. Hubble radius in Fig. 6. From the plot we can see that, since now the symmetry before and after the bounce is broken, the evolution of the perturbations will also be different. While the small k modes of perturbations still need to exit and reenter the horizon in the contracting phase, the large k modes needn't, and will stay inside the horizon until the expanding phase. In the following, we will focus on these modes.

          Figure 6.  (color online) The evolution of perturbation wavelengths with different comoving wavenumbers k as well as Hubble radius $ |H^{-1}| $. The figure is taken from Ref. [34].

          In the contracting phase where the EoS of the Universe oscillates continuously with $ \langle w\rangle=0 $, the background and perturbation evolutions are still described by Eqs. (14) and (15). While the solution is given in Eq. (16), since we are only interested in the subhorizon solution, we only have:

          $ \Phi_k=4\pi G\frac{\sqrt{\rho_0}\eta_0^3}{\eta^3}\frac{{\rm e}^{-{\rm i}k\eta}}{\sqrt{2k^3}}\; . $

          (33)

          Since in this case the evolution is not symmetric and for large k modes there is no horizon-exit behavior before bounce, we don't have a "slow-climbing" phase. In bounce phase, the equation is the same as (28), and the solution is given by (29). After the bounce, the universe will enter a slow-rolling expanding phase, in which the perturbation equation is (19) and the solution is the same as (30).

          The coefficients $ E_k $, $ F_k $, $ G_k $ and $ H_k $ are obtained using the same matching method as above. However, since here the matchings are only performed at sub-horizon regions, the final perturbations are determined by the constant mode of $ H_k $. From the tedious calculation done in [34], one obtains

          $ H_k=\sqrt{\frac{\pi}{2}}\frac{4\pi G}{\sqrt{2k^3}}|\dot\phi|{\rm e}^{-{\rm i}k\tilde{\eta}_{B+}}\left\{1+\frac{3{\rm e}^{-{\rm i}k(\eta_{B-}-\tilde{\eta}_{B+})}}{k\eta_{B-}}\sin[k(\eta_{B-}-\tilde{\eta}_{B+})]\right\}\; . $

          (34)

          where $ \eta_{B-} $ is the time when the Universe enters the bounce phase. Moreover, in the super-horizon region with nearly de-Sitter spacetime we have $ \zeta\simeq \Phi/\epsilon_H $, Thus one obtains the power spectrum of the curvature perturbations

          $ P_\zeta\equiv\frac{k^3}{2\pi^2}|\zeta|^2=\frac{8G^2\rho}{3\epsilon_H}\left\{1-\frac{3{\cal{H}}_{B-}}{2k}\sin\frac{2k}{{\cal{H}}_{B+}}\right\}\; . $

          (35)

          This result indicates that for large k (small scale) region where the last term could be neglected, the power spectrum is both constant and scale-invariant, while for small k (large scale) region where the second term becomes important, the spectrum gets some wiggle-like corrections. This is due to the effect of perturbations during bounce region. Moreover, for even larger scales, as the perturbations exit the horizon before the bounce, the power spectrum will be blue-tilted and will thus get suppressed [5254]. There are already some hints of these signals in the observations [5558]. If they are furtherly confirmed, they will provide a smoking gun for a bouncing cosmology.

        • B.   Small field potential

        • The typical small field potential is the Coleman-Weinberg potential, $ V(\phi)=(\lambda\phi^4/4)\ln(|\phi|/v)-\lambda(\phi^4- v^4)/16 $ [34, 51]. As a result of this potential, the vacuum is shifted to $ \phi=\pm v $. Therefore, the symmetry is spontaneously broken in the vacuum. It will also give rise to an asymmetric bounce behavior.

          Figure 5 shows the evolution behavior of w over the process in this model. Initially, the ϕ field is still set around the minimum value, and it can oscillate with its amplitude increasing and $ \langle w\rangle=0 $. When the field reaches the plateau, the σ field catches up, and the bounce takes place when w goes to minus infinity. After the bounce, when the ϕ-field rolls along the plateau, the Universe enters into a slow-roll phase where $ w\simeq-1 $, very much like inflation. Finally when the field drops into another minimum and oscillates, the inflation ends and the field becomes oscillating again, with $ \langle w\rangle=0 $.

          Figure 5.  (color online) Plot of evolution of EoS w in a double-field quintom bounce with small field potential. The figure is taken from Ref. [34].

          We also draw the evolution plot of perturbation wavelength scales v.s. Hubble radius in Fig. 6. From the plot we can see that, since now the symmetry before and after the bounce is broken, the evolution of the perturbations will also be different. While the small k modes of perturbations still need to exit and reenter the horizon in the contracting phase, the large k modes needn't, and will stay inside the horizon until the expanding phase. In the following, we will focus on these modes.

          Figure 6.  (color online) The evolution of perturbation wavelengths with different comoving wavenumbers k as well as Hubble radius $ |H^{-1}| $. The figure is taken from Ref. [34].

          In the contracting phase where the EoS of the Universe oscillates continuously with $ \langle w\rangle=0 $, the background and perturbation evolutions are still described by Eqs. (14) and (15). While the solution is given in Eq. (16), since we are only interested in the subhorizon solution, we only have:

          $ \Phi_k=4\pi G\frac{\sqrt{\rho_0}\eta_0^3}{\eta^3}\frac{{\rm e}^{-{\rm i}k\eta}}{\sqrt{2k^3}}\; . $

          (33)

          Since in this case the evolution is not symmetric and for large k modes there is no horizon-exit behavior before bounce, we don't have a "slow-climbing" phase. In bounce phase, the equation is the same as (28), and the solution is given by (29). After the bounce, the universe will enter a slow-rolling expanding phase, in which the perturbation equation is (19) and the solution is the same as (30).

          The coefficients $ E_k $, $ F_k $, $ G_k $ and $ H_k $ are obtained using the same matching method as above. However, since here the matchings are only performed at sub-horizon regions, the final perturbations are determined by the constant mode of $ H_k $. From the tedious calculation done in [34], one obtains

          $ H_k=\sqrt{\frac{\pi}{2}}\frac{4\pi G}{\sqrt{2k^3}}|\dot\phi|{\rm e}^{-{\rm i}k\tilde{\eta}_{B+}}\left\{1+\frac{3{\rm e}^{-{\rm i}k(\eta_{B-}-\tilde{\eta}_{B+})}}{k\eta_{B-}}\sin[k(\eta_{B-}-\tilde{\eta}_{B+})]\right\}\; . $

          (34)

          where $ \eta_{B-} $ is the time when the Universe enters the bounce phase. Moreover, in the super-horizon region with nearly de-Sitter spacetime we have $ \zeta\simeq \Phi/\epsilon_H $, Thus one obtains the power spectrum of the curvature perturbations

          $ P_\zeta\equiv\frac{k^3}{2\pi^2}|\zeta|^2=\frac{8G^2\rho}{3\epsilon_H}\left\{1-\frac{3{\cal{H}}_{B-}}{2k}\sin\frac{2k}{{\cal{H}}_{B+}}\right\}\; . $

          (35)

          This result indicates that for large k (small scale) region where the last term could be neglected, the power spectrum is both constant and scale-invariant, while for small k (large scale) region where the second term becomes important, the spectrum gets some wiggle-like corrections. This is due to the effect of perturbations during bounce region. Moreover, for even larger scales, as the perturbations exit the horizon before the bounce, the power spectrum will be blue-tilted and will thus get suppressed [5254]. There are already some hints of these signals in the observations [5558]. If they are furtherly confirmed, they will provide a smoking gun for a bouncing cosmology.

        II.   BOUNCE WITH A HIGHER-DERIVATIVE FIELD
        • The single field Quintom model with higher-derivative interactions was proposed in [59]. With higher-derivative terms, the single field gains an additional degree of freedom that can make the EoS cross −1. When applied to the early Universe, it can also give rise to a bounce solution.

          A detailed analysis of a bounce model with higher-derivative interactions was performed in Ref. [60]. In this work, the Lagrangian of the single field $ \hat{\phi} $ is written as:

          $ L=\frac{1}{2}\nabla_\mu\hat{\phi}\nabla^\mu\hat{\phi}-\frac{1}{2M^2}(\Box\hat{\phi})^2-\frac{1}{2}m^2\hat{\phi}^2\; . $

          (36)

          In particle physics, this Lagrangian is called the ''Lee-Wick model'', in which the Higgs mass is stabilized against radiative corrections to solve the hierarchy problem [6163]. Note that this model can also be transformed into a two-scalar-field model. By introducing an auxiliary field $ \tilde{\phi} $ and redefining a scalar field as:

          $ \phi=\hat{\phi}+\tilde{\phi}\; , $

          (37)

          the Lagrangian (36) then becomes

          $ L=\frac{1}{2}\nabla_\mu\phi\nabla^\mu\phi-\frac{1}{2}\nabla_\mu\tilde{\phi}\nabla^\mu\tilde{\phi}-\frac{1}{2}m^2\phi^2+\frac{1}{2}M^2\tilde{\phi}^2\; , $

          (38)

          where the equations of motion (EoM) for the ϕ and $ \tilde{\phi} $ fields are

          $ \ddot\phi+3H\dot\phi+m^2\phi=0\; ,\; \; \; \ddot{\tilde{\phi}}+3H\dot{\tilde{\phi}}+M^2\tilde{\phi}=0\; , $

          (39)

          respectively. From the EoM one can see that, since both fields now get masses, in the contracting phase they will both oscillate around their vacuum with increasing amplitude,

          $ |\phi(t)|\sim |\tilde{\phi}(t)|\sim a^{-3/2}(t)\; , $

          (40)

          and will freeze out when their amplitude reaches $ (12\pi)^{-1/2}m_{pl} $. At the initial time let’s assume the energy density of $ \phi(t) $ is larger than $ \tilde{\phi}(t) $, so that $ \phi(t) $ will dominate the Universe. When $ \phi(t) $ is oscillating, the Universe will behave like matter with its average value of EoS $ \langle w\rangle=0 $. Similarly to the large field bounce story, when the ϕ-field reaches $ (12\pi)^{-1/2}m_{pl} $ and then freezes out, it will slow-climb along the potential and tend to drive a period of “deflation”. But what is different is that, now the second field $ \tilde{\phi} $ also has mass and will also oscillate, which will destroy the deflation behavior. If the $ \tilde{\phi} $ field oscillates until its energy density catches up with that of ϕ, the total energy density will vanish and the bounce will occur. If the bounce happens before $ \tilde{\phi} $ freezes out, there will be no deflation era and we will get a pure matter bounce [46, 64, 65]. Moreover, since both fields have symmetric potentials, the evolution of the Universe will also be symmetric before and after the bounce, namely the Universe will directly enter into a matter-dominated era without inflation. The evolutions of both fields and the EoS of the Universe are plotted in Fig. 7.

          Figure 7.  (color online) Plot of evolutions of fields ϕ, $ \tilde{\phi} $ and EoS w in the Lee-Wick type quintom bounce. The figure is taken from Ref. [60].

          The perturbation equation is the same as in case of the two-field models, namely Eq. (11). However, as is described above, the cosmological periods that the Universe experiences are different. Here in Fig. 8 we plot the perturbation wavelength scales vs. Hubble radius (in conformal time). The process is just from matter dominance to the bounce, then to matter dominance again. As is demonstrated above, the initial condition is given by the Bunch-Davies vacuum solution, namely Eq. (33).

          Figure 8.  (color online) The comoving wavelength $ \lambda =1/k $ as well as the comoving Hubble radius $ |{\cal{H}}^{-1}| $. The figure is taken from Ref. [60].

          In the previous section, we have shown that if the matter contracting phase is followed by an inflationary phase after the bounce, where the perturbations did not exit the horizon before the bounce, a scale-invariant power spectrum is acquired. This is because the evolution inside of the horizon actually does not affect the averaged k-dependence of the perturbations and it still behaves like the Bunch-Davies vacuum. However, in this model, in the absence of inflation era, one may not be able to obtain a scale-invariant power spectrum. As we will see below, this issue can be overcome by making the perturbations exit the horizon before the bounce. In this case, the growing mode of Φ, which dominates in the contracting phase, gets deeply red-tilted, and when crossing the bounce point, it will be compensated by the same amount of blue-tilting.

          The super-horizon approximation of the contracting phase was obtained previously

          $ \Phi=D_-+S_-(\eta-\tilde{\eta}_{B-})^{-5}\; , $

          (41)

          where $ D_-\sim k^{3/2} $, $ S_-\sim k^{-7/2} $, $ \tilde{\eta}_{B-}\equiv \eta_{B-}-2/{\cal{H}}_{B-} $. While the first mode is constant, the second mode is growing.

          During the bounce phase, the equation of motion is given by Eq. (28). However, since we are considering the fluctuation modes which exit the horizon in the heating phase, unlike the case of the previous solution (29), now we consider the solution with much smaller k, namely $ k^2\ll y $. Thus the y-term will dominate $ k^2 $ term and the solution will present approximately no k-dependence. This is a very important property of the large-scale solution which will affect the scale-dependence as follows. The solution is given by

          $ \Phi_k=F_k+E_k\sqrt{y}(\eta-\eta_B)-\frac{3}{5}F_ky(\eta-\eta_B)^2+{\cal O}(y^{3/2}(\eta-\eta_B)^3)\; , $

          (42)

          where we assume that the bounce period is very short and it is enough to expand up to $ {\cal O}(\eta-\eta_B)^3 $.

          In the expanding phase, the solution is similar to the contracting phase since both have an EoS of matter

          $ \Phi=D_++S_+(\eta-\tilde{\eta}_{B+})^{-5}\; , $

          (43)

          where $ \tilde{\eta}_{B+}\equiv \eta_{B+}-2/{\cal{H}}_{B+} $. Now the second mode becomes decaying, so the constant mode will be dominant. Matching the solutions (41), (42), and (43) one obtains the final result for $ D_+ $:

          $ \begin{split} D_+=\;& D_-+\left[-\frac{4}{5}D_-+\frac{3}{5}\frac{{\cal{H}}_{B-}^5}{2^4}S_-\right]\frac{2k^2}{9{\cal{H}}_{B-}^2}\; ,\\ = \;& -\frac{\sqrt{\rho_{B-}}}{10\sqrt{2}}k^{-3/2}\; , \end{split} $

          (44)

          where in the last line we made use of $ D_-\sim k^{3/2} $, $ S_-\sim k^{-7/2} $. Here $ \rho_{B-} $ is the energy density at $ \eta_{B-} $. Therefore the power spectrum turns out to be

          $ P_\Phi\equiv\frac{k^3}{2\pi^2}|D_+|^2=\frac{\rho_{B-}}{(20\pi)^2}\; , $

          (45)

          which is scale-invariant. From Eq. (44) one can see that, this is due to the fact that in the small k region, while the constant mode is suppressed, the growing mode, which is deep red, will contribute to $ D+ $, and will be blue-tilted by the right amount. In Fig. 9 we plot $ P_\Phi $ and the spectral index $ n_s $ with respect to k. We can see from the plot that, while in the large k region the power spectrum is blue-tilted, in the small k region it will be scale-invariant.

          Figure 9.  (color online) Left: numerical results for the power spectrum of curvature perturbations $ P_\Phi $ with respect to time (t for cosmic time and η for comoving time). Right: the power spectrum $ P_\Phi $ as well as the spectral index $ n_s $ with respect to comoving wavenumber k. The figures are taken from Ref. [60].

          For the matter bounce, initial fluctuations could also be generated classically by the thermal matter perturbations inside the Hubble radius [60, 6669]. From the perturbed Einstein equations, the Fourier space correlation function of the Newtonian potential Φ is related to density fluctuations $ \langle|\delta T^0_0(k)|\rangle^2 $ (or mass fluctuations in position space $ \delta M(R) $),

          $\begin{split} \langle|\Phi(k)|\rangle^2=\;&16\pi^2G^2k^{-4}a^4\langle|\delta T^0_0(k)|\rangle^2\\ = \;& 16\pi^2G^2k^{-4}a^4R^{-3}\delta M(R)^2\; ,\end{split} $

          (46)

          where $ R=a/k $ is the physical radius of the region corresponding to the wavenumber k. On the other hand, $ \delta M(R)^2\sim T^5 R^3 $. Therefore one has

          $ P_\Phi(k)\sim k^3\langle|\Phi(k)|\rangle^2\sim k^{-1}T^5\sim k^{-1}a^{-5}\; , $

          (47)

          where we used the fact $ T\sim a^{-1} $ for thermal matter. At horizon crossing time $ t_c $ where $ a(t_c)H(t_c)\sim k $, one can get $ a\sim k^{-(1+3w)/2} $. For matter contraction $ \langle|w|\rangle=0 $, this leads to $ P_\Phi\sim k^9 $ (constant mode) and $ k^{-1} $ (growing mode) [60].

          One can also discuss the tensor perturbations in the bounce model, which contribute to the primordial gravitational waves and will have signals in the B-mode polarizations of the CMB map. Since the tensor perturbation is also generated in the contracting phase and propagates during the contracting and bounce phases, it is expected to be different from that generated in the normal inflationary scenario. To see this, one can start with the FRW metric with tensorial perturbation

          $ {\rm d}s^2=a^2(\eta)[{\rm d}\eta^2-(\delta_{ij}+h_{ij}){\rm d}x^i{\rm d}x^j]\; , $

          (48)

          with the properties $ h_{ij}=h_{ji} $, $ h_{ii}=0 $ and $ h_{ij,j}=0 $. It is well known that, in Einstein’s gravity, it satisfies the equation of motion

          $ h_{ij}''+2\frac{a'}{a}h_{ij}'-\nabla^2h_{ij}=0\; . $

          (49)

          Considering the Fourier modes of $ h_{ij} $

          $ h_{ij}=\sum\limits_{\lambda=1,2}\int \frac{{\rm d}^3k}{(2\pi)^{3/2}}h^{(\lambda)}_k(\eta,{\bf k}){\rm e}^{(\lambda)}_{ij}{\rm e}^{{\rm i}{\bf k}{\bf x}}\; , $

          (50)

          and defining $ v_k=ah_k/\sqrt{2} $ (omitting superscript "$ (\lambda) $"), one gets the Mukhanov-Sasaki equation for $ v_k $:

          $ v_k''+\left(k^2-\frac{a''}{a}\right)v_k=0\; . $

          (51)

          The procedure for dealing with this equation is the same as that for scalar perturbations. The solution is divided into three stages, namely the contracting and expanding phases with $ w=0 $, and the bounce phase. For the former, the behavior of $ a(\eta) $ is given by Eq. (14), so one has $ a''/a\propto \eta^{-2} $. For the latter one has

          $ \frac{a''}{a}\simeq \frac{4}{\pi}\alpha a_B^2=\frac{y}{3}\; . $

          (52)

          One can then solve Eq. (51) for each stage and connect with matching conditions. The final solution for $ v_k $ will take an asymptotic form,

          $ v_k\simeq {\rm i} \frac{\sqrt{2}}{24}\frac{{\cal{H}}^3_{B+}}{k^{3/2}}(\eta-\tilde{\eta}_{B+})^2\; . $

          (53)

          Then the tensor power spectrum is given by

          $ P_T(k)\equiv 2\frac{k^3}{\pi^2}\left|\frac{v_k}{a}\right|^2=\frac{2\rho_{B+}}{27\pi^2}\; , $

          (54)

          indicating a scale-invariant and constant tensor power spectrum [60, 70].

          In 2008, a class of single scalar models called the Galileon has been proposed [7173]. These models contain higher derivative terms in their Lagrangians, however, due to the delicate design of the Lagrangians (for example, the coupling of the higher derivative term to the kinetic term), their equations of motion can be made 2nd order, therefore avoiding a dynamical ghost mode. These models can be extended to include modified gravity or nonminimal coupling, which is a realization of an idea pioneered by Horndeski in 1974 [74, 75]. Due to this interesting property, these models have been applied to many aspects in cosmology, see Ref. [7685].

          A Galileon model can also realize quintom behavior and drive a bounce [8697]. Here we give a concrete model:

          $ \begin{split} S=\;&\int {\rm d}^4x\sqrt{-g}\left[\frac{R}{16\pi G}+F^2 {\rm e}^{2\phi}(\partial\phi)^2+\frac{F^3}{M^3}(\partial\phi)^2\Box\phi\right.\\ &\left.+\frac{F^3}{2M^3}(\partial\phi)^4\right]\; , \end{split}$

          (55)

          which is called the conformal Galileon model [98]. The coefficient F is the mass scale of the Galileon field, while M is the scale of higher derivative operators. Thus we can get the energy density and pressure as

          $ \begin{aligned} \rho=&\ F^2\left[-{\rm e}^{2\phi}\dot\phi^2+\left(\frac{3F}{2M^3}\right)(\dot\phi^4+4H\dot\phi^3)\right]\; , \end{aligned} $

          (56)

          $ \begin{aligned} p=&\ F^2\left[-{\rm e}^{2\phi}\dot\phi^2+\left(\frac{F}{2M^3}\right)(\dot\phi^4-4\dot\phi^2\ddot\phi)\right]\; . \end{aligned} $

          (57)

          In the absence of potential terms in ρ and p, the Friedmann equation can be easily solved to get an analytical solution in the contracting phase ($ H<0 $):

          $ \dot\phi\sim \frac{1}{\sqrt{2(t_0-t)}}\; ,\; \; \; H=\frac{1}{2(t-t_0)}\; . $

          (58)

          This indicates that this model will give a bounce solution with $ w=1/3 $ (radiation-like) in the contracting phase. In Fig. 10 we plot the scale factor from the contracting phase to the expanding phase.

          Figure 10.  (color online) The evolution of the scale factor a with respect to cosmic time t in a Galileon type quintom bounce. The figure is taken from Ref. [86].

          To see how the ghost mode is eliminated in this kind of model, it is useful to apply the ADM formalism [99] to write down the second order perturbed action. In this formalism, the perturbed metric is

          $ {\rm d}s^2=-N^2{\rm d}t^2+a^2{\rm e}^{2\zeta}\delta_{ij}({\rm d}x^i+N^i{\rm d}t)({\rm d}x^j+N^j{\rm d}t)\; , $

          (59)

          where $ N(x) $ and $ N_i(x) $ are the lapse function and shift vector, while ζ is the curvature perturbation. Moreover, in uniform-ϕ gauge we set $ \delta\phi=0 $. Then the 2nd order perturbed equation for the curvature perturbation is [100]:

          $ S_2=3\int {\rm d}t{\rm d}^3x a^3DM_p^2\left[\dot\zeta^2-\frac{c_s^2}{a^2}(\partial_i\zeta)^2\right]\; , $

          (60)

          where

          $ \begin{aligned} D=&\ \frac{2M_p^4H^2+2(F/M)^6\dot\phi^6+(F/M)^3M_p^2\dot\phi^4}{2[M_p^2H-(F/M)^3\dot\phi^3]^2}\; , \end{aligned} $

          (61)

          $ \begin{aligned} c_s^2=&\ \frac{-2M_p^4\dot H+2(F/M)^3M_p^2H\dot\phi^3-2(F/M)^6\dot\phi^6+6(F/M)^3M_p^2\dot\phi^2\ddot\phi}{3[2M_p^4H^2+(F/M)^3M_p^2\dot\phi^4+2(F/M)^6\dot\phi^6]}\; . \end{aligned} $

          (62)

          From this result, we can see that the factor in front of the kinetic term of the perturbation ζ is positive definite, therefore there is no ghost mode. We also plot D and $ c_s^2 $ in the bounce solution in Fig. 11 and Fig. 12. We can see that, the parameter D stays above 0 throughout the process, indicating that there is no ghost instability. The sound speed squared $ c_s^2 $ is also positive in most of the region from the contraction to the expansion phase, especially until the bounce phase, but will go below 0 after the bounce, which would break down the Galileon effective theory. Although the exact time point is model-dependent and it may not necessarily ruin the bounce, the existence of the negative-$ c_s^2 $ period is proved to be general in cubic Galileon theories [101103]. However, in [93, 94, 104107] it has been shown that, when some modified gravity effects beyond cubic action are included, $ c_s^2 $ can be made positive at all times without affecting the background.

          Figure 11.  (color online) The evolution of the parameter D with respect to cosmic time t in a Galileon type quintom bounce. The figure is taken from Ref. [86].

          Figure 12.  (color online) The evolution of the parameter $ c_s^2 $ with respect to cosmic time t in a Galileon type quintom bounce. Left: $ c_s^2 $ in the contracting and bouncing phases. Right: $ c_s^2 $ after the bounce. The figure is taken from Ref. [86].

          Moreover, as a bounce solution with a radiation-like contracting phase, it cannot give rise to the scale-invariant power spectrum required by the observational data. To remedy this defect, we introduce another curvaton field which has nonminimal coupling to the Galileon field to generate the power spectrum. A usual curvaton is a light field in an inflationary phase, which is used to give rise to a scale-invariant spectrum instead of the background field [108110]. But in a radiation-like contracting background, the curvaton needs to couple to the background field to make itself “feel” that it was in an inflationary background. A possible action of the coupled curvaton is [86, 89, 96, 98]

          $ S_\sigma=\int {\rm d}^4x\sqrt{-g}[-\phi^q(\partial\sigma)^2]\; , $

          (63)

          where σ is the curvaton field, while q is a parameter. The equation of motion for the curvaton is

          $ \sigma''+\left(2{\cal{H}}+\frac{q\phi'}{\phi}\right)\sigma'-\nabla^2\sigma=0\; . $

          (64)

          Since, in a radiation-dominant phase, from solution (58) one has $ \phi\sim \eta $, $ {\cal{H}}\sim \eta $, then the solution of Eq. (64) (in momentum space) is

          $ \sigma_k\sim c_1 k^{-(1+q)/2}|\eta|^{-(1+q)}+c_2k^{(1+q)/2} $

          (65)

          in the super-horizon approximation. The power spectrum can be scale-invariant if:

          1) $ q=2 $:

          $ \sigma_k\sim c_1 k^{-3/2}|\eta|^{-3}+c_2k^{3/2}\; ,\; \; \; P_\sigma\equiv\frac{k^3}{2\pi^2}|\sigma|^2\simeq \frac{c_1^2}{2\pi^2}|\eta|^{-6}\; , $

          (66)

          or 2) $ q=-4 $:

          $ \sigma_k\sim c_1 k^{3/2}|\eta|^3+c_2k^{-3/2}\; ,\; \; \; P_\sigma\equiv\frac{k^3}{2\pi^2}|\sigma|^2\simeq \frac{c_2^2}{2\pi^2}\; . $

          (67)

          However, the equation of motion (64) for the background gives the energy density for σ field: $ \rho_\sigma \sim (t_0-t)^{-3-q/2} $, while the radiation-like background scales as $ \rho_{bg}\sim (t_0-t)^{-2} $. To obtain a stable evolution where the energy density of the curvaton field does not overwhelm the background and ruin the bounce, we require $ q\geqslant -2 $. Therefore the solution with $ q=2 $ is chosen.

          Meanwhile, one can also calculate the tensor power spectrum in this model. Since the gravity part is not altered from Einstein’s Gravity, the equation of motion for the tensor power spectrum is the same as the previous one, namely Eq. (51). For the radiation-dominated background $ a\sim \eta $, one has

          $ P^T(k)\equiv \frac{k^3}{2\pi^2}\left|\frac{v_k}{a}\right|^2\sim k^2\; ,\; \; \; n_T\equiv \frac{d\ln P^T(k)}{d\ln k}=2\; . $

          (68)

          Therefore, a blue-tilted tensor power spectrum is obtained. However, it obviously depends on the contracting background as well as the gravity sector of the model. When the model is different, it is possible to obtain scale-invariant tensor spectrum as well.

        II.   BOUNCE WITH A HIGHER-DERIVATIVE FIELD
        • The single field Quintom model with higher-derivative interactions was proposed in [59]. With higher-derivative terms, the single field gains an additional degree of freedom that can make the EoS cross −1. When applied to the early Universe, it can also give rise to a bounce solution.

          A detailed analysis of a bounce model with higher-derivative interactions was performed in Ref. [60]. In this work, the Lagrangian of the single field $ \hat{\phi} $ is written as:

          $ L=\frac{1}{2}\nabla_\mu\hat{\phi}\nabla^\mu\hat{\phi}-\frac{1}{2M^2}(\Box\hat{\phi})^2-\frac{1}{2}m^2\hat{\phi}^2\; . $

          (36)

          In particle physics, this Lagrangian is called the ''Lee-Wick model'', in which the Higgs mass is stabilized against radiative corrections to solve the hierarchy problem [6163]. Note that this model can also be transformed into a two-scalar-field model. By introducing an auxiliary field $ \tilde{\phi} $ and redefining a scalar field as:

          $ \phi=\hat{\phi}+\tilde{\phi}\; , $

          (37)

          the Lagrangian (36) then becomes

          $ L=\frac{1}{2}\nabla_\mu\phi\nabla^\mu\phi-\frac{1}{2}\nabla_\mu\tilde{\phi}\nabla^\mu\tilde{\phi}-\frac{1}{2}m^2\phi^2+\frac{1}{2}M^2\tilde{\phi}^2\; , $

          (38)

          where the equations of motion (EoM) for the ϕ and $ \tilde{\phi} $ fields are

          $ \ddot\phi+3H\dot\phi+m^2\phi=0\; ,\; \; \; \ddot{\tilde{\phi}}+3H\dot{\tilde{\phi}}+M^2\tilde{\phi}=0\; , $

          (39)

          respectively. From the EoM one can see that, since both fields now get masses, in the contracting phase they will both oscillate around their vacuum with increasing amplitude,

          $ |\phi(t)|\sim |\tilde{\phi}(t)|\sim a^{-3/2}(t)\; , $

          (40)

          and will freeze out when their amplitude reaches $ (12\pi)^{-1/2}m_{pl} $. At the initial time let’s assume the energy density of $ \phi(t) $ is larger than $ \tilde{\phi}(t) $, so that $ \phi(t) $ will dominate the Universe. When $ \phi(t) $ is oscillating, the Universe will behave like matter with its average value of EoS $ \langle w\rangle=0 $. Similarly to the large field bounce story, when the ϕ-field reaches $ (12\pi)^{-1/2}m_{pl} $ and then freezes out, it will slow-climb along the potential and tend to drive a period of “deflation”. But what is different is that, now the second field $ \tilde{\phi} $ also has mass and will also oscillate, which will destroy the deflation behavior. If the $ \tilde{\phi} $ field oscillates until its energy density catches up with that of ϕ, the total energy density will vanish and the bounce will occur. If the bounce happens before $ \tilde{\phi} $ freezes out, there will be no deflation era and we will get a pure matter bounce [46, 64, 65]. Moreover, since both fields have symmetric potentials, the evolution of the Universe will also be symmetric before and after the bounce, namely the Universe will directly enter into a matter-dominated era without inflation. The evolutions of both fields and the EoS of the Universe are plotted in Fig. 7.

          Figure 7.  (color online) Plot of evolutions of fields ϕ, $ \tilde{\phi} $ and EoS w in the Lee-Wick type quintom bounce. The figure is taken from Ref. [60].

          The perturbation equation is the same as in case of the two-field models, namely Eq. (11). However, as is described above, the cosmological periods that the Universe experiences are different. Here in Fig. 8 we plot the perturbation wavelength scales vs. Hubble radius (in conformal time). The process is just from matter dominance to the bounce, then to matter dominance again. As is demonstrated above, the initial condition is given by the Bunch-Davies vacuum solution, namely Eq. (33).

          Figure 8.  (color online) The comoving wavelength $ \lambda =1/k $ as well as the comoving Hubble radius $ |{\cal{H}}^{-1}| $. The figure is taken from Ref. [60].

          In the previous section, we have shown that if the matter contracting phase is followed by an inflationary phase after the bounce, where the perturbations did not exit the horizon before the bounce, a scale-invariant power spectrum is acquired. This is because the evolution inside of the horizon actually does not affect the averaged k-dependence of the perturbations and it still behaves like the Bunch-Davies vacuum. However, in this model, in the absence of inflation era, one may not be able to obtain a scale-invariant power spectrum. As we will see below, this issue can be overcome by making the perturbations exit the horizon before the bounce. In this case, the growing mode of Φ, which dominates in the contracting phase, gets deeply red-tilted, and when crossing the bounce point, it will be compensated by the same amount of blue-tilting.

          The super-horizon approximation of the contracting phase was obtained previously

          $ \Phi=D_-+S_-(\eta-\tilde{\eta}_{B-})^{-5}\; , $

          (41)

          where $ D_-\sim k^{3/2} $, $ S_-\sim k^{-7/2} $, $ \tilde{\eta}_{B-}\equiv \eta_{B-}-2/{\cal{H}}_{B-} $. While the first mode is constant, the second mode is growing.

          During the bounce phase, the equation of motion is given by Eq. (28). However, since we are considering the fluctuation modes which exit the horizon in the heating phase, unlike the case of the previous solution (29), now we consider the solution with much smaller k, namely $ k^2\ll y $. Thus the y-term will dominate $ k^2 $ term and the solution will present approximately no k-dependence. This is a very important property of the large-scale solution which will affect the scale-dependence as follows. The solution is given by

          $ \Phi_k=F_k+E_k\sqrt{y}(\eta-\eta_B)-\frac{3}{5}F_ky(\eta-\eta_B)^2+{\cal O}(y^{3/2}(\eta-\eta_B)^3)\; , $

          (42)

          where we assume that the bounce period is very short and it is enough to expand up to $ {\cal O}(\eta-\eta_B)^3 $.

          In the expanding phase, the solution is similar to the contracting phase since both have an EoS of matter

          $ \Phi=D_++S_+(\eta-\tilde{\eta}_{B+})^{-5}\; , $

          (43)

          where $ \tilde{\eta}_{B+}\equiv \eta_{B+}-2/{\cal{H}}_{B+} $. Now the second mode becomes decaying, so the constant mode will be dominant. Matching the solutions (41), (42), and (43) one obtains the final result for $ D_+ $:

          $ \begin{split} D_+=\;& D_-+\left[-\frac{4}{5}D_-+\frac{3}{5}\frac{{\cal{H}}_{B-}^5}{2^4}S_-\right]\frac{2k^2}{9{\cal{H}}_{B-}^2}\; ,\\ = \;& -\frac{\sqrt{\rho_{B-}}}{10\sqrt{2}}k^{-3/2}\; , \end{split} $

          (44)

          where in the last line we made use of $ D_-\sim k^{3/2} $, $ S_-\sim k^{-7/2} $. Here $ \rho_{B-} $ is the energy density at $ \eta_{B-} $. Therefore the power spectrum turns out to be

          $ P_\Phi\equiv\frac{k^3}{2\pi^2}|D_+|^2=\frac{\rho_{B-}}{(20\pi)^2}\; , $

          (45)

          which is scale-invariant. From Eq. (44) one can see that, this is due to the fact that in the small k region, while the constant mode is suppressed, the growing mode, which is deep red, will contribute to $ D+ $, and will be blue-tilted by the right amount. In Fig. 9 we plot $ P_\Phi $ and the spectral index $ n_s $ with respect to k. We can see from the plot that, while in the large k region the power spectrum is blue-tilted, in the small k region it will be scale-invariant.

          Figure 9.  (color online) Left: numerical results for the power spectrum of curvature perturbations $ P_\Phi $ with respect to time (t for cosmic time and η for comoving time). Right: the power spectrum $ P_\Phi $ as well as the spectral index $ n_s $ with respect to comoving wavenumber k. The figures are taken from Ref. [60].

          For the matter bounce, initial fluctuations could also be generated classically by the thermal matter perturbations inside the Hubble radius [60, 6669]. From the perturbed Einstein equations, the Fourier space correlation function of the Newtonian potential Φ is related to density fluctuations $ \langle|\delta T^0_0(k)|\rangle^2 $ (or mass fluctuations in position space $ \delta M(R) $),

          $\begin{split} \langle|\Phi(k)|\rangle^2=\;&16\pi^2G^2k^{-4}a^4\langle|\delta T^0_0(k)|\rangle^2\\ = \;& 16\pi^2G^2k^{-4}a^4R^{-3}\delta M(R)^2\; ,\end{split} $

          (46)

          where $ R=a/k $ is the physical radius of the region corresponding to the wavenumber k. On the other hand, $ \delta M(R)^2\sim T^5 R^3 $. Therefore one has

          $ P_\Phi(k)\sim k^3\langle|\Phi(k)|\rangle^2\sim k^{-1}T^5\sim k^{-1}a^{-5}\; , $

          (47)

          where we used the fact $ T\sim a^{-1} $ for thermal matter. At horizon crossing time $ t_c $ where $ a(t_c)H(t_c)\sim k $, one can get $ a\sim k^{-(1+3w)/2} $. For matter contraction $ \langle|w|\rangle=0 $, this leads to $ P_\Phi\sim k^9 $ (constant mode) and $ k^{-1} $ (growing mode) [60].

          One can also discuss the tensor perturbations in the bounce model, which contribute to the primordial gravitational waves and will have signals in the B-mode polarizations of the CMB map. Since the tensor perturbation is also generated in the contracting phase and propagates during the contracting and bounce phases, it is expected to be different from that generated in the normal inflationary scenario. To see this, one can start with the FRW metric with tensorial perturbation

          $ {\rm d}s^2=a^2(\eta)[{\rm d}\eta^2-(\delta_{ij}+h_{ij}){\rm d}x^i{\rm d}x^j]\; , $

          (48)

          with the properties $ h_{ij}=h_{ji} $, $ h_{ii}=0 $ and $ h_{ij,j}=0 $. It is well known that, in Einstein’s gravity, it satisfies the equation of motion

          $ h_{ij}''+2\frac{a'}{a}h_{ij}'-\nabla^2h_{ij}=0\; . $

          (49)

          Considering the Fourier modes of $ h_{ij} $

          $ h_{ij}=\sum\limits_{\lambda=1,2}\int \frac{{\rm d}^3k}{(2\pi)^{3/2}}h^{(\lambda)}_k(\eta,{\bf k}){\rm e}^{(\lambda)}_{ij}{\rm e}^{{\rm i}{\bf k}{\bf x}}\; , $

          (50)

          and defining $ v_k=ah_k/\sqrt{2} $ (omitting superscript "$ (\lambda) $"), one gets the Mukhanov-Sasaki equation for $ v_k $:

          $ v_k''+\left(k^2-\frac{a''}{a}\right)v_k=0\; . $

          (51)

          The procedure for dealing with this equation is the same as that for scalar perturbations. The solution is divided into three stages, namely the contracting and expanding phases with $ w=0 $, and the bounce phase. For the former, the behavior of $ a(\eta) $ is given by Eq. (14), so one has $ a''/a\propto \eta^{-2} $. For the latter one has

          $ \frac{a''}{a}\simeq \frac{4}{\pi}\alpha a_B^2=\frac{y}{3}\; . $

          (52)

          One can then solve Eq. (51) for each stage and connect with matching conditions. The final solution for $ v_k $ will take an asymptotic form,

          $ v_k\simeq {\rm i} \frac{\sqrt{2}}{24}\frac{{\cal{H}}^3_{B+}}{k^{3/2}}(\eta-\tilde{\eta}_{B+})^2\; . $

          (53)

          Then the tensor power spectrum is given by

          $ P_T(k)\equiv 2\frac{k^3}{\pi^2}\left|\frac{v_k}{a}\right|^2=\frac{2\rho_{B+}}{27\pi^2}\; , $

          (54)

          indicating a scale-invariant and constant tensor power spectrum [60, 70].

          In 2008, a class of single scalar models called the Galileon has been proposed [7173]. These models contain higher derivative terms in their Lagrangians, however, due to the delicate design of the Lagrangians (for example, the coupling of the higher derivative term to the kinetic term), their equations of motion can be made 2nd order, therefore avoiding a dynamical ghost mode. These models can be extended to include modified gravity or nonminimal coupling, which is a realization of an idea pioneered by Horndeski in 1974 [74, 75]. Due to this interesting property, these models have been applied to many aspects in cosmology, see Ref. [7685].

          A Galileon model can also realize quintom behavior and drive a bounce [8697]. Here we give a concrete model:

          $ \begin{split} S=\;&\int {\rm d}^4x\sqrt{-g}\left[\frac{R}{16\pi G}+F^2 {\rm e}^{2\phi}(\partial\phi)^2+\frac{F^3}{M^3}(\partial\phi)^2\Box\phi\right.\\ &\left.+\frac{F^3}{2M^3}(\partial\phi)^4\right]\; , \end{split}$

          (55)

          which is called the conformal Galileon model [98]. The coefficient F is the mass scale of the Galileon field, while M is the scale of higher derivative operators. Thus we can get the energy density and pressure as

          $ \begin{aligned} \rho=&\ F^2\left[-{\rm e}^{2\phi}\dot\phi^2+\left(\frac{3F}{2M^3}\right)(\dot\phi^4+4H\dot\phi^3)\right]\; , \end{aligned} $

          (56)

          $ \begin{aligned} p=&\ F^2\left[-{\rm e}^{2\phi}\dot\phi^2+\left(\frac{F}{2M^3}\right)(\dot\phi^4-4\dot\phi^2\ddot\phi)\right]\; . \end{aligned} $

          (57)

          In the absence of potential terms in ρ and p, the Friedmann equation can be easily solved to get an analytical solution in the contracting phase ($ H<0 $):

          $ \dot\phi\sim \frac{1}{\sqrt{2(t_0-t)}}\; ,\; \; \; H=\frac{1}{2(t-t_0)}\; . $

          (58)

          This indicates that this model will give a bounce solution with $ w=1/3 $ (radiation-like) in the contracting phase. In Fig. 10 we plot the scale factor from the contracting phase to the expanding phase.

          Figure 10.  (color online) The evolution of the scale factor a with respect to cosmic time t in a Galileon type quintom bounce. The figure is taken from Ref. [86].

          To see how the ghost mode is eliminated in this kind of model, it is useful to apply the ADM formalism [99] to write down the second order perturbed action. In this formalism, the perturbed metric is

          $ {\rm d}s^2=-N^2{\rm d}t^2+a^2{\rm e}^{2\zeta}\delta_{ij}({\rm d}x^i+N^i{\rm d}t)({\rm d}x^j+N^j{\rm d}t)\; , $

          (59)

          where $ N(x) $ and $ N_i(x) $ are the lapse function and shift vector, while ζ is the curvature perturbation. Moreover, in uniform-ϕ gauge we set $ \delta\phi=0 $. Then the 2nd order perturbed equation for the curvature perturbation is [100]:

          $ S_2=3\int {\rm d}t{\rm d}^3x a^3DM_p^2\left[\dot\zeta^2-\frac{c_s^2}{a^2}(\partial_i\zeta)^2\right]\; , $

          (60)

          where

          $ \begin{aligned} D=&\ \frac{2M_p^4H^2+2(F/M)^6\dot\phi^6+(F/M)^3M_p^2\dot\phi^4}{2[M_p^2H-(F/M)^3\dot\phi^3]^2}\; , \end{aligned} $

          (61)

          $ \begin{aligned} c_s^2=&\ \frac{-2M_p^4\dot H+2(F/M)^3M_p^2H\dot\phi^3-2(F/M)^6\dot\phi^6+6(F/M)^3M_p^2\dot\phi^2\ddot\phi}{3[2M_p^4H^2+(F/M)^3M_p^2\dot\phi^4+2(F/M)^6\dot\phi^6]}\; . \end{aligned} $

          (62)

          From this result, we can see that the factor in front of the kinetic term of the perturbation ζ is positive definite, therefore there is no ghost mode. We also plot D and $ c_s^2 $ in the bounce solution in Fig. 11 and Fig. 12. We can see that, the parameter D stays above 0 throughout the process, indicating that there is no ghost instability. The sound speed squared $ c_s^2 $ is also positive in most of the region from the contraction to the expansion phase, especially until the bounce phase, but will go below 0 after the bounce, which would break down the Galileon effective theory. Although the exact time point is model-dependent and it may not necessarily ruin the bounce, the existence of the negative-$ c_s^2 $ period is proved to be general in cubic Galileon theories [101103]. However, in [93, 94, 104107] it has been shown that, when some modified gravity effects beyond cubic action are included, $ c_s^2 $ can be made positive at all times without affecting the background.

          Figure 11.  (color online) The evolution of the parameter D with respect to cosmic time t in a Galileon type quintom bounce. The figure is taken from Ref. [86].

          Figure 12.  (color online) The evolution of the parameter $ c_s^2 $ with respect to cosmic time t in a Galileon type quintom bounce. Left: $ c_s^2 $ in the contracting and bouncing phases. Right: $ c_s^2 $ after the bounce. The figure is taken from Ref. [86].

          Moreover, as a bounce solution with a radiation-like contracting phase, it cannot give rise to the scale-invariant power spectrum required by the observational data. To remedy this defect, we introduce another curvaton field which has nonminimal coupling to the Galileon field to generate the power spectrum. A usual curvaton is a light field in an inflationary phase, which is used to give rise to a scale-invariant spectrum instead of the background field [108110]. But in a radiation-like contracting background, the curvaton needs to couple to the background field to make itself “feel” that it was in an inflationary background. A possible action of the coupled curvaton is [86, 89, 96, 98]

          $ S_\sigma=\int {\rm d}^4x\sqrt{-g}[-\phi^q(\partial\sigma)^2]\; , $

          (63)

          where σ is the curvaton field, while q is a parameter. The equation of motion for the curvaton is

          $ \sigma''+\left(2{\cal{H}}+\frac{q\phi'}{\phi}\right)\sigma'-\nabla^2\sigma=0\; . $

          (64)

          Since, in a radiation-dominant phase, from solution (58) one has $ \phi\sim \eta $, $ {\cal{H}}\sim \eta $, then the solution of Eq. (64) (in momentum space) is

          $ \sigma_k\sim c_1 k^{-(1+q)/2}|\eta|^{-(1+q)}+c_2k^{(1+q)/2} $

          (65)

          in the super-horizon approximation. The power spectrum can be scale-invariant if:

          1) $ q=2 $:

          $ \sigma_k\sim c_1 k^{-3/2}|\eta|^{-3}+c_2k^{3/2}\; ,\; \; \; P_\sigma\equiv\frac{k^3}{2\pi^2}|\sigma|^2\simeq \frac{c_1^2}{2\pi^2}|\eta|^{-6}\; , $

          (66)

          or 2) $ q=-4 $:

          $ \sigma_k\sim c_1 k^{3/2}|\eta|^3+c_2k^{-3/2}\; ,\; \; \; P_\sigma\equiv\frac{k^3}{2\pi^2}|\sigma|^2\simeq \frac{c_2^2}{2\pi^2}\; . $

          (67)

          However, the equation of motion (64) for the background gives the energy density for σ field: $ \rho_\sigma \sim (t_0-t)^{-3-q/2} $, while the radiation-like background scales as $ \rho_{bg}\sim (t_0-t)^{-2} $. To obtain a stable evolution where the energy density of the curvaton field does not overwhelm the background and ruin the bounce, we require $ q\geqslant -2 $. Therefore the solution with $ q=2 $ is chosen.

          Meanwhile, one can also calculate the tensor power spectrum in this model. Since the gravity part is not altered from Einstein’s Gravity, the equation of motion for the tensor power spectrum is the same as the previous one, namely Eq. (51). For the radiation-dominated background $ a\sim \eta $, one has

          $ P^T(k)\equiv \frac{k^3}{2\pi^2}\left|\frac{v_k}{a}\right|^2\sim k^2\; ,\; \; \; n_T\equiv \frac{d\ln P^T(k)}{d\ln k}=2\; . $

          (68)

          Therefore, a blue-tilted tensor power spectrum is obtained. However, it obviously depends on the contracting background as well as the gravity sector of the model. When the model is different, it is possible to obtain scale-invariant tensor spectrum as well.

        II.   BOUNCE WITH A HIGHER-DERIVATIVE FIELD
        • The single field Quintom model with higher-derivative interactions was proposed in [59]. With higher-derivative terms, the single field gains an additional degree of freedom that can make the EoS cross −1. When applied to the early Universe, it can also give rise to a bounce solution.

          A detailed analysis of a bounce model with higher-derivative interactions was performed in Ref. [60]. In this work, the Lagrangian of the single field $ \hat{\phi} $ is written as:

          $ L=\frac{1}{2}\nabla_\mu\hat{\phi}\nabla^\mu\hat{\phi}-\frac{1}{2M^2}(\Box\hat{\phi})^2-\frac{1}{2}m^2\hat{\phi}^2\; . $

          (36)

          In particle physics, this Lagrangian is called the ''Lee-Wick model'', in which the Higgs mass is stabilized against radiative corrections to solve the hierarchy problem [6163]. Note that this model can also be transformed into a two-scalar-field model. By introducing an auxiliary field $ \tilde{\phi} $ and redefining a scalar field as:

          $ \phi=\hat{\phi}+\tilde{\phi}\; , $

          (37)

          the Lagrangian (36) then becomes

          $ L=\frac{1}{2}\nabla_\mu\phi\nabla^\mu\phi-\frac{1}{2}\nabla_\mu\tilde{\phi}\nabla^\mu\tilde{\phi}-\frac{1}{2}m^2\phi^2+\frac{1}{2}M^2\tilde{\phi}^2\; , $

          (38)

          where the equations of motion (EoM) for the ϕ and $ \tilde{\phi} $ fields are

          $ \ddot\phi+3H\dot\phi+m^2\phi=0\; ,\; \; \; \ddot{\tilde{\phi}}+3H\dot{\tilde{\phi}}+M^2\tilde{\phi}=0\; , $

          (39)

          respectively. From the EoM one can see that, since both fields now get masses, in the contracting phase they will both oscillate around their vacuum with increasing amplitude,

          $ |\phi(t)|\sim |\tilde{\phi}(t)|\sim a^{-3/2}(t)\; , $

          (40)

          and will freeze out when their amplitude reaches $ (12\pi)^{-1/2}m_{pl} $. At the initial time let’s assume the energy density of $ \phi(t) $ is larger than $ \tilde{\phi}(t) $, so that $ \phi(t) $ will dominate the Universe. When $ \phi(t) $ is oscillating, the Universe will behave like matter with its average value of EoS $ \langle w\rangle=0 $. Similarly to the large field bounce story, when the ϕ-field reaches $ (12\pi)^{-1/2}m_{pl} $ and then freezes out, it will slow-climb along the potential and tend to drive a period of “deflation”. But what is different is that, now the second field $ \tilde{\phi} $ also has mass and will also oscillate, which will destroy the deflation behavior. If the $ \tilde{\phi} $ field oscillates until its energy density catches up with that of ϕ, the total energy density will vanish and the bounce will occur. If the bounce happens before $ \tilde{\phi} $ freezes out, there will be no deflation era and we will get a pure matter bounce [46, 64, 65]. Moreover, since both fields have symmetric potentials, the evolution of the Universe will also be symmetric before and after the bounce, namely the Universe will directly enter into a matter-dominated era without inflation. The evolutions of both fields and the EoS of the Universe are plotted in Fig. 7.

          Figure 7.  (color online) Plot of evolutions of fields ϕ, $ \tilde{\phi} $ and EoS w in the Lee-Wick type quintom bounce. The figure is taken from Ref. [60].

          The perturbation equation is the same as in case of the two-field models, namely Eq. (11). However, as is described above, the cosmological periods that the Universe experiences are different. Here in Fig. 8 we plot the perturbation wavelength scales vs. Hubble radius (in conformal time). The process is just from matter dominance to the bounce, then to matter dominance again. As is demonstrated above, the initial condition is given by the Bunch-Davies vacuum solution, namely Eq. (33).

          Figure 8.  (color online) The comoving wavelength $ \lambda =1/k $ as well as the comoving Hubble radius $ |{\cal{H}}^{-1}| $. The figure is taken from Ref. [60].

          In the previous section, we have shown that if the matter contracting phase is followed by an inflationary phase after the bounce, where the perturbations did not exit the horizon before the bounce, a scale-invariant power spectrum is acquired. This is because the evolution inside of the horizon actually does not affect the averaged k-dependence of the perturbations and it still behaves like the Bunch-Davies vacuum. However, in this model, in the absence of inflation era, one may not be able to obtain a scale-invariant power spectrum. As we will see below, this issue can be overcome by making the perturbations exit the horizon before the bounce. In this case, the growing mode of Φ, which dominates in the contracting phase, gets deeply red-tilted, and when crossing the bounce point, it will be compensated by the same amount of blue-tilting.

          The super-horizon approximation of the contracting phase was obtained previously

          $ \Phi=D_-+S_-(\eta-\tilde{\eta}_{B-})^{-5}\; , $

          (41)

          where $ D_-\sim k^{3/2} $, $ S_-\sim k^{-7/2} $, $ \tilde{\eta}_{B-}\equiv \eta_{B-}-2/{\cal{H}}_{B-} $. While the first mode is constant, the second mode is growing.

          During the bounce phase, the equation of motion is given by Eq. (28). However, since we are considering the fluctuation modes which exit the horizon in the heating phase, unlike the case of the previous solution (29), now we consider the solution with much smaller k, namely $ k^2\ll y $. Thus the y-term will dominate $ k^2 $ term and the solution will present approximately no k-dependence. This is a very important property of the large-scale solution which will affect the scale-dependence as follows. The solution is given by

          $ \Phi_k=F_k+E_k\sqrt{y}(\eta-\eta_B)-\frac{3}{5}F_ky(\eta-\eta_B)^2+{\cal O}(y^{3/2}(\eta-\eta_B)^3)\; , $

          (42)

          where we assume that the bounce period is very short and it is enough to expand up to $ {\cal O}(\eta-\eta_B)^3 $.

          In the expanding phase, the solution is similar to the contracting phase since both have an EoS of matter

          $ \Phi=D_++S_+(\eta-\tilde{\eta}_{B+})^{-5}\; , $

          (43)

          where $ \tilde{\eta}_{B+}\equiv \eta_{B+}-2/{\cal{H}}_{B+} $. Now the second mode becomes decaying, so the constant mode will be dominant. Matching the solutions (41), (42), and (43) one obtains the final result for $ D_+ $:

          $ \begin{split} D_+=\;& D_-+\left[-\frac{4}{5}D_-+\frac{3}{5}\frac{{\cal{H}}_{B-}^5}{2^4}S_-\right]\frac{2k^2}{9{\cal{H}}_{B-}^2}\; ,\\ = \;& -\frac{\sqrt{\rho_{B-}}}{10\sqrt{2}}k^{-3/2}\; , \end{split} $

          (44)

          where in the last line we made use of $ D_-\sim k^{3/2} $, $ S_-\sim k^{-7/2} $. Here $ \rho_{B-} $ is the energy density at $ \eta_{B-} $. Therefore the power spectrum turns out to be

          $ P_\Phi\equiv\frac{k^3}{2\pi^2}|D_+|^2=\frac{\rho_{B-}}{(20\pi)^2}\; , $

          (45)

          which is scale-invariant. From Eq. (44) one can see that, this is due to the fact that in the small k region, while the constant mode is suppressed, the growing mode, which is deep red, will contribute to $ D+ $, and will be blue-tilted by the right amount. In Fig. 9 we plot $ P_\Phi $ and the spectral index $ n_s $ with respect to k. We can see from the plot that, while in the large k region the power spectrum is blue-tilted, in the small k region it will be scale-invariant.

          Figure 9.  (color online) Left: numerical results for the power spectrum of curvature perturbations $ P_\Phi $ with respect to time (t for cosmic time and η for comoving time). Right: the power spectrum $ P_\Phi $ as well as the spectral index $ n_s $ with respect to comoving wavenumber k. The figures are taken from Ref. [60].

          For the matter bounce, initial fluctuations could also be generated classically by the thermal matter perturbations inside the Hubble radius [60, 6669]. From the perturbed Einstein equations, the Fourier space correlation function of the Newtonian potential Φ is related to density fluctuations $ \langle|\delta T^0_0(k)|\rangle^2 $ (or mass fluctuations in position space $ \delta M(R) $),

          $\begin{split} \langle|\Phi(k)|\rangle^2=\;&16\pi^2G^2k^{-4}a^4\langle|\delta T^0_0(k)|\rangle^2\\ = \;& 16\pi^2G^2k^{-4}a^4R^{-3}\delta M(R)^2\; ,\end{split} $

          (46)

          where $ R=a/k $ is the physical radius of the region corresponding to the wavenumber k. On the other hand, $ \delta M(R)^2\sim T^5 R^3 $. Therefore one has

          $ P_\Phi(k)\sim k^3\langle|\Phi(k)|\rangle^2\sim k^{-1}T^5\sim k^{-1}a^{-5}\; , $

          (47)

          where we used the fact $ T\sim a^{-1} $ for thermal matter. At horizon crossing time $ t_c $ where $ a(t_c)H(t_c)\sim k $, one can get $ a\sim k^{-(1+3w)/2} $. For matter contraction $ \langle|w|\rangle=0 $, this leads to $ P_\Phi\sim k^9 $ (constant mode) and $ k^{-1} $ (growing mode) [60].

          One can also discuss the tensor perturbations in the bounce model, which contribute to the primordial gravitational waves and will have signals in the B-mode polarizations of the CMB map. Since the tensor perturbation is also generated in the contracting phase and propagates during the contracting and bounce phases, it is expected to be different from that generated in the normal inflationary scenario. To see this, one can start with the FRW metric with tensorial perturbation

          $ {\rm d}s^2=a^2(\eta)[{\rm d}\eta^2-(\delta_{ij}+h_{ij}){\rm d}x^i{\rm d}x^j]\; , $

          (48)

          with the properties $ h_{ij}=h_{ji} $, $ h_{ii}=0 $ and $ h_{ij,j}=0 $. It is well known that, in Einstein’s gravity, it satisfies the equation of motion

          $ h_{ij}''+2\frac{a'}{a}h_{ij}'-\nabla^2h_{ij}=0\; . $

          (49)

          Considering the Fourier modes of $ h_{ij} $

          $ h_{ij}=\sum\limits_{\lambda=1,2}\int \frac{{\rm d}^3k}{(2\pi)^{3/2}}h^{(\lambda)}_k(\eta,{\bf k}){\rm e}^{(\lambda)}_{ij}{\rm e}^{{\rm i}{\bf k}{\bf x}}\; , $

          (50)

          and defining $ v_k=ah_k/\sqrt{2} $ (omitting superscript "$ (\lambda) $"), one gets the Mukhanov-Sasaki equation for $ v_k $:

          $ v_k''+\left(k^2-\frac{a''}{a}\right)v_k=0\; . $

          (51)

          The procedure for dealing with this equation is the same as that for scalar perturbations. The solution is divided into three stages, namely the contracting and expanding phases with $ w=0 $, and the bounce phase. For the former, the behavior of $ a(\eta) $ is given by Eq. (14), so one has $ a''/a\propto \eta^{-2} $. For the latter one has

          $ \frac{a''}{a}\simeq \frac{4}{\pi}\alpha a_B^2=\frac{y}{3}\; . $

          (52)

          One can then solve Eq. (51) for each stage and connect with matching conditions. The final solution for $ v_k $ will take an asymptotic form,

          $ v_k\simeq {\rm i} \frac{\sqrt{2}}{24}\frac{{\cal{H}}^3_{B+}}{k^{3/2}}(\eta-\tilde{\eta}_{B+})^2\; . $

          (53)

          Then the tensor power spectrum is given by

          $ P_T(k)\equiv 2\frac{k^3}{\pi^2}\left|\frac{v_k}{a}\right|^2=\frac{2\rho_{B+}}{27\pi^2}\; , $

          (54)

          indicating a scale-invariant and constant tensor power spectrum [60, 70].

          In 2008, a class of single scalar models called the Galileon has been proposed [7173]. These models contain higher derivative terms in their Lagrangians, however, due to the delicate design of the Lagrangians (for example, the coupling of the higher derivative term to the kinetic term), their equations of motion can be made 2nd order, therefore avoiding a dynamical ghost mode. These models can be extended to include modified gravity or nonminimal coupling, which is a realization of an idea pioneered by Horndeski in 1974 [74, 75]. Due to this interesting property, these models have been applied to many aspects in cosmology, see Ref. [7685].

          A Galileon model can also realize quintom behavior and drive a bounce [8697]. Here we give a concrete model:

          $ \begin{split} S=\;&\int {\rm d}^4x\sqrt{-g}\left[\frac{R}{16\pi G}+F^2 {\rm e}^{2\phi}(\partial\phi)^2+\frac{F^3}{M^3}(\partial\phi)^2\Box\phi\right.\\ &\left.+\frac{F^3}{2M^3}(\partial\phi)^4\right]\; , \end{split}$

          (55)

          which is called the conformal Galileon model [98]. The coefficient F is the mass scale of the Galileon field, while M is the scale of higher derivative operators. Thus we can get the energy density and pressure as

          $ \begin{aligned} \rho=&\ F^2\left[-{\rm e}^{2\phi}\dot\phi^2+\left(\frac{3F}{2M^3}\right)(\dot\phi^4+4H\dot\phi^3)\right]\; , \end{aligned} $

          (56)

          $ \begin{aligned} p=&\ F^2\left[-{\rm e}^{2\phi}\dot\phi^2+\left(\frac{F}{2M^3}\right)(\dot\phi^4-4\dot\phi^2\ddot\phi)\right]\; . \end{aligned} $

          (57)

          In the absence of potential terms in ρ and p, the Friedmann equation can be easily solved to get an analytical solution in the contracting phase ($ H<0 $):

          $ \dot\phi\sim \frac{1}{\sqrt{2(t_0-t)}}\; ,\; \; \; H=\frac{1}{2(t-t_0)}\; . $

          (58)

          This indicates that this model will give a bounce solution with $ w=1/3 $ (radiation-like) in the contracting phase. In Fig. 10 we plot the scale factor from the contracting phase to the expanding phase.

          Figure 10.  (color online) The evolution of the scale factor a with respect to cosmic time t in a Galileon type quintom bounce. The figure is taken from Ref. [86].

          To see how the ghost mode is eliminated in this kind of model, it is useful to apply the ADM formalism [99] to write down the second order perturbed action. In this formalism, the perturbed metric is

          $ {\rm d}s^2=-N^2{\rm d}t^2+a^2{\rm e}^{2\zeta}\delta_{ij}({\rm d}x^i+N^i{\rm d}t)({\rm d}x^j+N^j{\rm d}t)\; , $

          (59)

          where $ N(x) $ and $ N_i(x) $ are the lapse function and shift vector, while ζ is the curvature perturbation. Moreover, in uniform-ϕ gauge we set $ \delta\phi=0 $. Then the 2nd order perturbed equation for the curvature perturbation is [100]:

          $ S_2=3\int {\rm d}t{\rm d}^3x a^3DM_p^2\left[\dot\zeta^2-\frac{c_s^2}{a^2}(\partial_i\zeta)^2\right]\; , $

          (60)

          where

          $ \begin{aligned} D=&\ \frac{2M_p^4H^2+2(F/M)^6\dot\phi^6+(F/M)^3M_p^2\dot\phi^4}{2[M_p^2H-(F/M)^3\dot\phi^3]^2}\; , \end{aligned} $

          (61)

          $ \begin{aligned} c_s^2=&\ \frac{-2M_p^4\dot H+2(F/M)^3M_p^2H\dot\phi^3-2(F/M)^6\dot\phi^6+6(F/M)^3M_p^2\dot\phi^2\ddot\phi}{3[2M_p^4H^2+(F/M)^3M_p^2\dot\phi^4+2(F/M)^6\dot\phi^6]}\; . \end{aligned} $

          (62)

          From this result, we can see that the factor in front of the kinetic term of the perturbation ζ is positive definite, therefore there is no ghost mode. We also plot D and $ c_s^2 $ in the bounce solution in Fig. 11 and Fig. 12. We can see that, the parameter D stays above 0 throughout the process, indicating that there is no ghost instability. The sound speed squared $ c_s^2 $ is also positive in most of the region from the contraction to the expansion phase, especially until the bounce phase, but will go below 0 after the bounce, which would break down the Galileon effective theory. Although the exact time point is model-dependent and it may not necessarily ruin the bounce, the existence of the negative-$ c_s^2 $ period is proved to be general in cubic Galileon theories [101103]. However, in [93, 94, 104107] it has been shown that, when some modified gravity effects beyond cubic action are included, $ c_s^2 $ can be made positive at all times without affecting the background.

          Figure 11.  (color online) The evolution of the parameter D with respect to cosmic time t in a Galileon type quintom bounce. The figure is taken from Ref. [86].

          Figure 12.  (color online) The evolution of the parameter $ c_s^2 $ with respect to cosmic time t in a Galileon type quintom bounce. Left: $ c_s^2 $ in the contracting and bouncing phases. Right: $ c_s^2 $ after the bounce. The figure is taken from Ref. [86].

          Moreover, as a bounce solution with a radiation-like contracting phase, it cannot give rise to the scale-invariant power spectrum required by the observational data. To remedy this defect, we introduce another curvaton field which has nonminimal coupling to the Galileon field to generate the power spectrum. A usual curvaton is a light field in an inflationary phase, which is used to give rise to a scale-invariant spectrum instead of the background field [108110]. But in a radiation-like contracting background, the curvaton needs to couple to the background field to make itself “feel” that it was in an inflationary background. A possible action of the coupled curvaton is [86, 89, 96, 98]

          $ S_\sigma=\int {\rm d}^4x\sqrt{-g}[-\phi^q(\partial\sigma)^2]\; , $

          (63)

          where σ is the curvaton field, while q is a parameter. The equation of motion for the curvaton is

          $ \sigma''+\left(2{\cal{H}}+\frac{q\phi'}{\phi}\right)\sigma'-\nabla^2\sigma=0\; . $

          (64)

          Since, in a radiation-dominant phase, from solution (58) one has $ \phi\sim \eta $, $ {\cal{H}}\sim \eta $, then the solution of Eq. (64) (in momentum space) is

          $ \sigma_k\sim c_1 k^{-(1+q)/2}|\eta|^{-(1+q)}+c_2k^{(1+q)/2} $

          (65)

          in the super-horizon approximation. The power spectrum can be scale-invariant if:

          1) $ q=2 $:

          $ \sigma_k\sim c_1 k^{-3/2}|\eta|^{-3}+c_2k^{3/2}\; ,\; \; \; P_\sigma\equiv\frac{k^3}{2\pi^2}|\sigma|^2\simeq \frac{c_1^2}{2\pi^2}|\eta|^{-6}\; , $

          (66)

          or 2) $ q=-4 $:

          $ \sigma_k\sim c_1 k^{3/2}|\eta|^3+c_2k^{-3/2}\; ,\; \; \; P_\sigma\equiv\frac{k^3}{2\pi^2}|\sigma|^2\simeq \frac{c_2^2}{2\pi^2}\; . $

          (67)

          However, the equation of motion (64) for the background gives the energy density for σ field: $ \rho_\sigma \sim (t_0-t)^{-3-q/2} $, while the radiation-like background scales as $ \rho_{bg}\sim (t_0-t)^{-2} $. To obtain a stable evolution where the energy density of the curvaton field does not overwhelm the background and ruin the bounce, we require $ q\geqslant -2 $. Therefore the solution with $ q=2 $ is chosen.

          Meanwhile, one can also calculate the tensor power spectrum in this model. Since the gravity part is not altered from Einstein’s Gravity, the equation of motion for the tensor power spectrum is the same as the previous one, namely Eq. (51). For the radiation-dominated background $ a\sim \eta $, one has

          $ P^T(k)\equiv \frac{k^3}{2\pi^2}\left|\frac{v_k}{a}\right|^2\sim k^2\; ,\; \; \; n_T\equiv \frac{d\ln P^T(k)}{d\ln k}=2\; . $

          (68)

          Therefore, a blue-tilted tensor power spectrum is obtained. However, it obviously depends on the contracting background as well as the gravity sector of the model. When the model is different, it is possible to obtain scale-invariant tensor spectrum as well.

        III.   BOUNCE WITH MODIFIED GRAVITY
        • Bouncing cosmologies can also be realized in modified gravity models [111]. In modified gravity, the Friedmann equations are no longer of the form (1), while the modified gravity will act as effective ρ and p. Thus even in the absence of exotic matter, they can still effectively realize a quintom scenario and bounce the Universe. A large class of modified gravity models is called metric affine gravity, which contains other geometric quantities such as the torsion tensor [112118]:

          $ T^\alpha_{\mu\nu}=\Gamma^\alpha _{\nu\mu}-\Gamma^\alpha_{\mu\nu}\; , $

          (69)

          whose existence means that the connection is no longer symmetric, and the nonmetricity tensor [119121]:

          $ Q_{\alpha\mu\nu}\equiv\nabla_\alpha g_{\mu\nu}=\frac{\partial g_{\mu\nu}}{\partial x^\alpha}-g_{\nu\sigma}\Gamma^{\sigma}_{\mu\alpha}-g_{\sigma\mu}\Gamma^{\sigma}_{\nu\alpha}\; , $

          (70)

          whose existence means that the covariant derivative of metric is no longer vanishing. Both of these tensors can be contracted into scalars:

          $ \begin{split} T\equiv\;& -\frac{1}{4}T^\rho_{\,\mu\nu}T^{\mu\nu}_{\,\,\rho}+\frac{1}{4}T^\rho_{\,\mu\nu}T^{\nu\mu}_{\,\,\rho}+\frac{1}{4}T^\rho_{\,\mu\nu}T^{\,\mu\nu}_{\rho}+\frac{1}{2}T^\mu_{\,\mu\nu}T^{\alpha\nu}_{\,\,\alpha}\\ &-\frac{1}{2}T^\nu_{\,\mu\nu}T^{\alpha\mu}_{\,\,\alpha}\; , \end{split} $

          (71)

          $ \begin{aligned} Q\equiv&\ \frac{1}{2}Q_{\mu\nu\lambda}Q^{\lambda\mu\nu}-\frac{1}{4}Q_{\mu\nu\lambda}Q^{\mu\nu\lambda}+\frac{1}{4}Q^{\,\,\,\,\mu}_{\lambda\mu}Q^{\lambda\,\,\mu}_{\,\,\mu}-\frac{1}{2}Q^{\,\,\,\,\mu}_{\lambda\mu}Q^{\mu}_{\,\,\mu\lambda}\; , \end{aligned} $

          (72)

          which are called the torsion scalar and the nonmetricity scalar, respectively. Moreover, people interestingly found that these two scalars have a simple relation to the Ricci scalar R:

          $ R=T-2\hat{\nabla}^\nu T^\mu_{\,\mu\nu}\; ,\; \; \; R=Q-\hat{\nabla}_\alpha(Q^{\alpha\,\,\mu}_{\,\,\mu}-Q^{\mu\,\,\alpha}_{\,\,\mu})\; , $

          (73)

          where $ \hat{\nabla} $ denotes the covariant derivative with respect to the Christoffel symbol. For this reason, when these scalars act as the gravity sector of the Universe’s action, in the minimal case they will be mathematically equivalent to each other [122125]:

          $ \int {\rm d}^4x\sqrt{-g}R=\int {\rm d}^4x\sqrt{-g}T=\int {\rm d}^4x\sqrt{-g}Q\; , $

          (74)

          and all of the physical laws originating from this action are the same. Thus the gravity theories of T and Q are equivalent to GR. However, when the actions are extended to arbitrary functions of R, T and Q, they will be different from each other since the total derivative cannot be removed

          $ \int {\rm d}^4x\sqrt{-g}f(R)\neq\int {\rm d}^4x\sqrt{-g}f(T)\neq\int {\rm d}^4x\sqrt{-g}f(Q)\; . $

          (75)

          Then we will have different actions with different physical consequences [126128]. For example, when applied to FRW spacetime where the metric is $ {\rm d}s^2=-{\rm d}t^2+a^2(t)\delta_{ij}{\rm d}x^i{\rm d}x^j $, the Friedmann equations for the cases of T and Q should be:

          $ \begin{aligned} 6f_TH^2+\frac{1}{2}f=&\ \rho_m\; , \end{aligned} $

          (76)

          $ \begin{aligned} (f_T-12H^2f_{TT})\dot H=\ -\frac{1}{2}(\rho_m+p_m)\; ,\; \; \; (or\; T\rightarrow Q ) \end{aligned} $

          (77)

          where $ T=Q=-6H^2 $ 4. Moreover, the Friedmann equations for $ f(R) $ will be [129]:

          $ \begin{aligned} 3f_RH^2+3H\dot f_R-\frac{1}{2}f_RR+\frac{1}{2}f=&\ \rho_m\; , \end{aligned} $

          (78)

          $ \begin{aligned} f_R\dot H+\frac{1}{2}\ddot f_R-\frac{1}{2}H\dot f_R=&\ -\frac{1}{2}(\rho_m+p_m)\; . \end{aligned} $

          (79)

          Various works on constructing bounce models in framework of modified gravity theories such as $ f(R) $, $ f(T) $ and $ f(Q) $ has been done [130142]. For such theories, it is convenient to construct a bounce scenario by finding an Ansatz bounce solution first, then reconstructing the forms of the functionals. We consider a simple Ansatz

          $ a(t)=a_B(1+\alpha t^2)^{\tfrac{1}{3(\gamma+1)}}\; , $

          (80)

          therefore the equation of state far from the bounce is $ w=\gamma $. According to Eq. (80), the Hubble parameter is given by

          $ H(t)=\frac{\dot a}{a}=\frac{2\alpha t}{3(1+\gamma)(1+\alpha t^2)}\; . $

          (81)

          In principle, it is straightforward to substitute (81) into the Friedmann equations (76), (77) to obtain the functional of f. However, the calculation process is difficult and we usually do not have an analytical solution over the whole period, but only have a numerical solution, or analytical solutions in certain periods, for which we refer the readers to e. g. Ref. [131, 142].

          More detailed calculations and results can be found in Ref. [131, 142], where realistic bounce models were obtained, and the evolutions of background and perturbations are discussed extensively. Although the action is different, at the leading order in perturbation theory, one can get similar equation of motion as in field theory models, and a nearly scale-invariant power spectrum can be obtained. See also [143148] for bounce models constructed from other modified gravity theories.

        III.   BOUNCE WITH MODIFIED GRAVITY
        • Bouncing cosmologies can also be realized in modified gravity models [111]. In modified gravity, the Friedmann equations are no longer of the form (1), while the modified gravity will act as effective ρ and p. Thus even in the absence of exotic matter, they can still effectively realize a quintom scenario and bounce the Universe. A large class of modified gravity models is called metric affine gravity, which contains other geometric quantities such as the torsion tensor [112118]:

          $ T^\alpha_{\mu\nu}=\Gamma^\alpha _{\nu\mu}-\Gamma^\alpha_{\mu\nu}\; , $

          (69)

          whose existence means that the connection is no longer symmetric, and the nonmetricity tensor [119121]:

          $ Q_{\alpha\mu\nu}\equiv\nabla_\alpha g_{\mu\nu}=\frac{\partial g_{\mu\nu}}{\partial x^\alpha}-g_{\nu\sigma}\Gamma^{\sigma}_{\mu\alpha}-g_{\sigma\mu}\Gamma^{\sigma}_{\nu\alpha}\; , $

          (70)

          whose existence means that the covariant derivative of metric is no longer vanishing. Both of these tensors can be contracted into scalars:

          $ \begin{split} T\equiv\;& -\frac{1}{4}T^\rho_{\,\mu\nu}T^{\mu\nu}_{\,\,\rho}+\frac{1}{4}T^\rho_{\,\mu\nu}T^{\nu\mu}_{\,\,\rho}+\frac{1}{4}T^\rho_{\,\mu\nu}T^{\,\mu\nu}_{\rho}+\frac{1}{2}T^\mu_{\,\mu\nu}T^{\alpha\nu}_{\,\,\alpha}\\ &-\frac{1}{2}T^\nu_{\,\mu\nu}T^{\alpha\mu}_{\,\,\alpha}\; , \end{split} $

          (71)

          $ \begin{aligned} Q\equiv&\ \frac{1}{2}Q_{\mu\nu\lambda}Q^{\lambda\mu\nu}-\frac{1}{4}Q_{\mu\nu\lambda}Q^{\mu\nu\lambda}+\frac{1}{4}Q^{\,\,\,\,\mu}_{\lambda\mu}Q^{\lambda\,\,\mu}_{\,\,\mu}-\frac{1}{2}Q^{\,\,\,\,\mu}_{\lambda\mu}Q^{\mu}_{\,\,\mu\lambda}\; , \end{aligned} $

          (72)

          which are called the torsion scalar and the nonmetricity scalar, respectively. Moreover, people interestingly found that these two scalars have a simple relation to the Ricci scalar R:

          $ R=T-2\hat{\nabla}^\nu T^\mu_{\,\mu\nu}\; ,\; \; \; R=Q-\hat{\nabla}_\alpha(Q^{\alpha\,\,\mu}_{\,\,\mu}-Q^{\mu\,\,\alpha}_{\,\,\mu})\; , $

          (73)

          where $ \hat{\nabla} $ denotes the covariant derivative with respect to the Christoffel symbol. For this reason, when these scalars act as the gravity sector of the Universe’s action, in the minimal case they will be mathematically equivalent to each other [122125]:

          $ \int {\rm d}^4x\sqrt{-g}R=\int {\rm d}^4x\sqrt{-g}T=\int {\rm d}^4x\sqrt{-g}Q\; , $

          (74)

          and all of the physical laws originating from this action are the same. Thus the gravity theories of T and Q are equivalent to GR. However, when the actions are extended to arbitrary functions of R, T and Q, they will be different from each other since the total derivative cannot be removed

          $ \int {\rm d}^4x\sqrt{-g}f(R)\neq\int {\rm d}^4x\sqrt{-g}f(T)\neq\int {\rm d}^4x\sqrt{-g}f(Q)\; . $

          (75)

          Then we will have different actions with different physical consequences [126128]. For example, when applied to FRW spacetime where the metric is $ {\rm d}s^2=-{\rm d}t^2+a^2(t)\delta_{ij}{\rm d}x^i{\rm d}x^j $, the Friedmann equations for the cases of T and Q should be:

          $ \begin{aligned} 6f_TH^2+\frac{1}{2}f=&\ \rho_m\; , \end{aligned} $

          (76)

          $ \begin{aligned} (f_T-12H^2f_{TT})\dot H=\ -\frac{1}{2}(\rho_m+p_m)\; ,\; \; \; (or\; T\rightarrow Q ) \end{aligned} $

          (77)

          where $ T=Q=-6H^2 $ 4. Moreover, the Friedmann equations for $ f(R) $ will be [129]:

          $ \begin{aligned} 3f_RH^2+3H\dot f_R-\frac{1}{2}f_RR+\frac{1}{2}f=&\ \rho_m\; , \end{aligned} $

          (78)

          $ \begin{aligned} f_R\dot H+\frac{1}{2}\ddot f_R-\frac{1}{2}H\dot f_R=&\ -\frac{1}{2}(\rho_m+p_m)\; . \end{aligned} $

          (79)

          Various works on constructing bounce models in framework of modified gravity theories such as $ f(R) $, $ f(T) $ and $ f(Q) $ has been done [130142]. For such theories, it is convenient to construct a bounce scenario by finding an Ansatz bounce solution first, then reconstructing the forms of the functionals. We consider a simple Ansatz

          $ a(t)=a_B(1+\alpha t^2)^{\tfrac{1}{3(\gamma+1)}}\; , $

          (80)

          therefore the equation of state far from the bounce is $ w=\gamma $. According to Eq. (80), the Hubble parameter is given by

          $ H(t)=\frac{\dot a}{a}=\frac{2\alpha t}{3(1+\gamma)(1+\alpha t^2)}\; . $

          (81)

          In principle, it is straightforward to substitute (81) into the Friedmann equations (76), (77) to obtain the functional of f. However, the calculation process is difficult and we usually do not have an analytical solution over the whole period, but only have a numerical solution, or analytical solutions in certain periods, for which we refer the readers to e. g. Ref. [131, 142].

          More detailed calculations and results can be found in Ref. [131, 142], where realistic bounce models were obtained, and the evolutions of background and perturbations are discussed extensively. Although the action is different, at the leading order in perturbation theory, one can get similar equation of motion as in field theory models, and a nearly scale-invariant power spectrum can be obtained. See also [143148] for bounce models constructed from other modified gravity theories.

        III.   BOUNCE WITH MODIFIED GRAVITY
        • Bouncing cosmologies can also be realized in modified gravity models [111]. In modified gravity, the Friedmann equations are no longer of the form (1), while the modified gravity will act as effective ρ and p. Thus even in the absence of exotic matter, they can still effectively realize a quintom scenario and bounce the Universe. A large class of modified gravity models is called metric affine gravity, which contains other geometric quantities such as the torsion tensor [112118]:

          $ T^\alpha_{\mu\nu}=\Gamma^\alpha _{\nu\mu}-\Gamma^\alpha_{\mu\nu}\; , $

          (69)

          whose existence means that the connection is no longer symmetric, and the nonmetricity tensor [119121]:

          $ Q_{\alpha\mu\nu}\equiv\nabla_\alpha g_{\mu\nu}=\frac{\partial g_{\mu\nu}}{\partial x^\alpha}-g_{\nu\sigma}\Gamma^{\sigma}_{\mu\alpha}-g_{\sigma\mu}\Gamma^{\sigma}_{\nu\alpha}\; , $

          (70)

          whose existence means that the covariant derivative of metric is no longer vanishing. Both of these tensors can be contracted into scalars:

          $ \begin{split} T\equiv\;& -\frac{1}{4}T^\rho_{\,\mu\nu}T^{\mu\nu}_{\,\,\rho}+\frac{1}{4}T^\rho_{\,\mu\nu}T^{\nu\mu}_{\,\,\rho}+\frac{1}{4}T^\rho_{\,\mu\nu}T^{\,\mu\nu}_{\rho}+\frac{1}{2}T^\mu_{\,\mu\nu}T^{\alpha\nu}_{\,\,\alpha}\\ &-\frac{1}{2}T^\nu_{\,\mu\nu}T^{\alpha\mu}_{\,\,\alpha}\; , \end{split} $

          (71)

          $ \begin{aligned} Q\equiv&\ \frac{1}{2}Q_{\mu\nu\lambda}Q^{\lambda\mu\nu}-\frac{1}{4}Q_{\mu\nu\lambda}Q^{\mu\nu\lambda}+\frac{1}{4}Q^{\,\,\,\,\mu}_{\lambda\mu}Q^{\lambda\,\,\mu}_{\,\,\mu}-\frac{1}{2}Q^{\,\,\,\,\mu}_{\lambda\mu}Q^{\mu}_{\,\,\mu\lambda}\; , \end{aligned} $

          (72)

          which are called the torsion scalar and the nonmetricity scalar, respectively. Moreover, people interestingly found that these two scalars have a simple relation to the Ricci scalar R:

          $ R=T-2\hat{\nabla}^\nu T^\mu_{\,\mu\nu}\; ,\; \; \; R=Q-\hat{\nabla}_\alpha(Q^{\alpha\,\,\mu}_{\,\,\mu}-Q^{\mu\,\,\alpha}_{\,\,\mu})\; , $

          (73)

          where $ \hat{\nabla} $ denotes the covariant derivative with respect to the Christoffel symbol. For this reason, when these scalars act as the gravity sector of the Universe’s action, in the minimal case they will be mathematically equivalent to each other [122125]:

          $ \int {\rm d}^4x\sqrt{-g}R=\int {\rm d}^4x\sqrt{-g}T=\int {\rm d}^4x\sqrt{-g}Q\; , $

          (74)

          and all of the physical laws originating from this action are the same. Thus the gravity theories of T and Q are equivalent to GR. However, when the actions are extended to arbitrary functions of R, T and Q, they will be different from each other since the total derivative cannot be removed

          $ \int {\rm d}^4x\sqrt{-g}f(R)\neq\int {\rm d}^4x\sqrt{-g}f(T)\neq\int {\rm d}^4x\sqrt{-g}f(Q)\; . $

          (75)

          Then we will have different actions with different physical consequences [126128]. For example, when applied to FRW spacetime where the metric is $ {\rm d}s^2=-{\rm d}t^2+a^2(t)\delta_{ij}{\rm d}x^i{\rm d}x^j $, the Friedmann equations for the cases of T and Q should be:

          $ \begin{aligned} 6f_TH^2+\frac{1}{2}f=&\ \rho_m\; , \end{aligned} $

          (76)

          $ \begin{aligned} (f_T-12H^2f_{TT})\dot H=\ -\frac{1}{2}(\rho_m+p_m)\; ,\; \; \; (or\; T\rightarrow Q ) \end{aligned} $

          (77)

          where $ T=Q=-6H^2 $ 4. Moreover, the Friedmann equations for $ f(R) $ will be [129]:

          $ \begin{aligned} 3f_RH^2+3H\dot f_R-\frac{1}{2}f_RR+\frac{1}{2}f=&\ \rho_m\; , \end{aligned} $

          (78)

          $ \begin{aligned} f_R\dot H+\frac{1}{2}\ddot f_R-\frac{1}{2}H\dot f_R=&\ -\frac{1}{2}(\rho_m+p_m)\; . \end{aligned} $

          (79)

          Various works on constructing bounce models in framework of modified gravity theories such as $ f(R) $, $ f(T) $ and $ f(Q) $ has been done [130142]. For such theories, it is convenient to construct a bounce scenario by finding an Ansatz bounce solution first, then reconstructing the forms of the functionals. We consider a simple Ansatz

          $ a(t)=a_B(1+\alpha t^2)^{\tfrac{1}{3(\gamma+1)}}\; , $

          (80)

          therefore the equation of state far from the bounce is $ w=\gamma $. According to Eq. (80), the Hubble parameter is given by

          $ H(t)=\frac{\dot a}{a}=\frac{2\alpha t}{3(1+\gamma)(1+\alpha t^2)}\; . $

          (81)

          In principle, it is straightforward to substitute (81) into the Friedmann equations (76), (77) to obtain the functional of f. However, the calculation process is difficult and we usually do not have an analytical solution over the whole period, but only have a numerical solution, or analytical solutions in certain periods, for which we refer the readers to e. g. Ref. [131, 142].

          More detailed calculations and results can be found in Ref. [131, 142], where realistic bounce models were obtained, and the evolutions of background and perturbations are discussed extensively. Although the action is different, at the leading order in perturbation theory, one can get similar equation of motion as in field theory models, and a nearly scale-invariant power spectrum can be obtained. See also [143148] for bounce models constructed from other modified gravity theories.

        D.   Beyond a bounce: Cyclic scenarios
        • A simple and precise example of a quintom cyclic model is given by the two-field action [29]:

          $ S=\int {\rm d}^4x\sqrt{-g}\left[\frac{1}{2}\partial_\mu\phi\partial^\mu\phi-\frac{1}{2}\partial_\mu\psi\partial^\mu\psi-V(\phi,\psi)\right]\; , $

          (82)

          Thus the energy density and pressure of the universe are

          $ \rho=\frac{1}{2}\dot\phi^2-\frac{1}{2}\dot\psi^2+V(\phi,\psi)\; ,\; \; \; p=\frac{1}{2}\dot\phi^2-\frac{1}{2}\dot\psi^2-V(\phi,\psi)\; , $

          (83)

          while the equations of motion for the ϕ and ψ fields are

          $ \ddot\phi+3H\dot\phi+V_{,\phi}=0\; ,\; \; \; \ddot\psi+3H\dot\psi-V_{,\psi}=0\; . $

          (84)

          To realize a cyclic scenario, we consider a coupled potential of the form

          $ V(\phi,\psi)=(\Lambda_0+\lambda\phi\psi)^2+\frac{1}{2}m^2\phi^2-\frac{1}{2}m^2\psi^2\; . $

          (85)

          Due to the symmetry of the potential under the transformation: $ \phi\rightarrow i\psi $, $ \psi\rightarrow -i\phi $, the two scalar fields will dominate alternately. Due to the couplings between two fields, the potential will be bounded from below. For Eq. (84) (with Hubble parameter satisfying the Friedmann equation and the energy density comes from Eq. (83)), we can obtain an Ansatz solution

          $ \phi=\sqrt{A_0}\cos mt\; ,\; \; \; \psi=\sqrt{A_0}\sin mt\; , $

          (86)

          where $ A_0 $ is the oscillation amplitude of the fields. Moreover, the Hubble parameter turns out to be:

          $ H=\frac{\sqrt{3}}{3M_p}(\Lambda_0+\Lambda_1\sin 2mt)\; , $

          (87)

          where we define $ \Lambda_1=\sqrt{3}mA_0/(4M_p) $.

          From Eq. (87) one can see that the Hubble parameter shows an oscillating behavior, and since $ H=d\ln a/{\rm d}t $, the scale factor can be oscillating as well. Moreover, different choices of the parameters $ \Lambda_0 $ and $ \Lambda_1 $ can give rise to different oscillating behaviors. In Fig. 13 we provide 4 possible kinds of oscillating behaviors of $ \ln a $, H and EoS w given by Eq. (87). The first case is for $ \Lambda_0=0 $ where both H and $ \ln a $ are oscillating with an equal amplitude. This is a standard oscillating Universe, but it suffers from the entropy problem. The second case is for $ 0<\Lambda_0<\Lambda_1 $. In this case, the center value of the Hubble parameter is positive, so the period of $ H<0 $ in each cycle will be shorter than that of $ H>0 $, therefore, the $ \ln a $ is oscillating with the total amplitude increasing. The third case is for $ \Lambda_0\geqslant \Lambda_1 $. In this case the Hubble parameter is always larger than zero, so the Universe is always expanding with increasing $ \ln a $. However, the expansion will be accelerating or decelerating alternately. Since there is no transition process from contraction to expansion, the EoS will also be regular all the time. The last case is for $ -|\Lambda_1|<\Lambda_0<0 $ where the center value of Hubble parameter is negative. In this case, period of $ H<0 $ in each cycle will be longer than that of $ H>0 $, therefore, the $ \ln a $ will be decreasing in total, but have short periods of increase in each cycle. Therefore, it still has room to explain our current expanding universe.

          Figure 13.  (color online) The evolution of $ \ln a $, H, ρ and w with respect to cosmic time t, corresponding to blue, purple, red and green lines, respectively. In the numerical calculation, the parameters are taken to be $ m=3\times 10^{-3} $ and $ \Lambda_1\simeq 0.39 $, and $ \Lambda_0 $ is taken as 0 (left-top), $ 0.10 $ (right-top), $ 0.402 $ (left-bottom) and $ -0.10 $ (right-bottom) respectively, corresponding to the cases of $ \Lambda_0=0 $, $ 0<\Lambda_0<\Lambda_1 $, $ \Lambda_0\geq\Lambda_1 $ and $ -|\Lambda_1|<\Lambda_0<0 $. The figures are taken from Ref. [29].

          Apart from bouncing and cyclic scenarios, one should note that a non-singular Universe can also be realized via emergent behavior, when the Universe approaches a static model in the infinite past [149151]. In this case, $ H\rightarrow 0 $ and $ \rho\rightarrow 0 $, while in order to enter into a realistic expanding phase, H still needs to gradually increase, with $ \dot H>0 $. Thus the EoS of the Universe will go from $ -\infty $ to above -1 as well. Therefore, in the emergent scenario, the Universe also realizes quintom-like behavior [152, 153].

          In conclusion, in this paper we briefly reviewed the Quintom model of dark energy and its applications to constructing bouncing cosmologies. Specifically, we have considered three classes of quintom scenarios: a two-field model, a single scalar with higher derivatives and modified gravity.

          Recently, other non-singular bounce scenarios making use of matter with EoS crossing -1 have been proposed. In Ref. [154] Bardeen, Bars, Hanson and Peccei (BBHP) introduced the folded string, which is a classical solution of a limit of the Nambu model. This was generalized by Itzhaki [155], who found that in a classical (1+1)-dimensional string model with an increasing dilaton field, there are Instant Folded String (IFS) solutions. These are essentially the BBHP folded strings, with the roles of time and space reversed so that they spontaneously nucleate 5. The IFS violates the null energy condition (NEC) [158], a fact which was used to create a non-singular bouncing cosmology in Ref. [159]. In Ref. [160], Itzhaki, Peleg and Steinhardt showed that by introducing a particular potential for the dilaton, together with an extra field that becomes massless at an enhanced symmetry point [161], the bounce of Ref. [159] can be incorporated into a cyclic bouncing cosmology, complete with reheating performed by the extra field.

        D.   Beyond a bounce: Cyclic scenarios
        • A simple and precise example of a quintom cyclic model is given by the two-field action [29]:

          $ S=\int {\rm d}^4x\sqrt{-g}\left[\frac{1}{2}\partial_\mu\phi\partial^\mu\phi-\frac{1}{2}\partial_\mu\psi\partial^\mu\psi-V(\phi,\psi)\right]\; , $

          (82)

          Thus the energy density and pressure of the universe are

          $ \rho=\frac{1}{2}\dot\phi^2-\frac{1}{2}\dot\psi^2+V(\phi,\psi)\; ,\; \; \; p=\frac{1}{2}\dot\phi^2-\frac{1}{2}\dot\psi^2-V(\phi,\psi)\; , $

          (83)

          while the equations of motion for the ϕ and ψ fields are

          $ \ddot\phi+3H\dot\phi+V_{,\phi}=0\; ,\; \; \; \ddot\psi+3H\dot\psi-V_{,\psi}=0\; . $

          (84)

          To realize a cyclic scenario, we consider a coupled potential of the form

          $ V(\phi,\psi)=(\Lambda_0+\lambda\phi\psi)^2+\frac{1}{2}m^2\phi^2-\frac{1}{2}m^2\psi^2\; . $

          (85)

          Due to the symmetry of the potential under the transformation: $ \phi\rightarrow i\psi $, $ \psi\rightarrow -i\phi $, the two scalar fields will dominate alternately. Due to the couplings between two fields, the potential will be bounded from below. For Eq. (84) (with Hubble parameter satisfying the Friedmann equation and the energy density comes from Eq. (83)), we can obtain an Ansatz solution

          $ \phi=\sqrt{A_0}\cos mt\; ,\; \; \; \psi=\sqrt{A_0}\sin mt\; , $

          (86)

          where $ A_0 $ is the oscillation amplitude of the fields. Moreover, the Hubble parameter turns out to be:

          $ H=\frac{\sqrt{3}}{3M_p}(\Lambda_0+\Lambda_1\sin 2mt)\; , $

          (87)

          where we define $ \Lambda_1=\sqrt{3}mA_0/(4M_p) $.

          From Eq. (87) one can see that the Hubble parameter shows an oscillating behavior, and since $ H=d\ln a/{\rm d}t $, the scale factor can be oscillating as well. Moreover, different choices of the parameters $ \Lambda_0 $ and $ \Lambda_1 $ can give rise to different oscillating behaviors. In Fig. 13 we provide 4 possible kinds of oscillating behaviors of $ \ln a $, H and EoS w given by Eq. (87). The first case is for $ \Lambda_0=0 $ where both H and $ \ln a $ are oscillating with an equal amplitude. This is a standard oscillating Universe, but it suffers from the entropy problem. The second case is for $ 0<\Lambda_0<\Lambda_1 $. In this case, the center value of the Hubble parameter is positive, so the period of $ H<0 $ in each cycle will be shorter than that of $ H>0 $, therefore, the $ \ln a $ is oscillating with the total amplitude increasing. The third case is for $ \Lambda_0\geqslant \Lambda_1 $. In this case the Hubble parameter is always larger than zero, so the Universe is always expanding with increasing $ \ln a $. However, the expansion will be accelerating or decelerating alternately. Since there is no transition process from contraction to expansion, the EoS will also be regular all the time. The last case is for $ -|\Lambda_1|<\Lambda_0<0 $ where the center value of Hubble parameter is negative. In this case, period of $ H<0 $ in each cycle will be longer than that of $ H>0 $, therefore, the $ \ln a $ will be decreasing in total, but have short periods of increase in each cycle. Therefore, it still has room to explain our current expanding universe.

          Figure 13.  (color online) The evolution of $ \ln a $, H, ρ and w with respect to cosmic time t, corresponding to blue, purple, red and green lines, respectively. In the numerical calculation, the parameters are taken to be $ m=3\times 10^{-3} $ and $ \Lambda_1\simeq 0.39 $, and $ \Lambda_0 $ is taken as 0 (left-top), $ 0.10 $ (right-top), $ 0.402 $ (left-bottom) and $ -0.10 $ (right-bottom) respectively, corresponding to the cases of $ \Lambda_0=0 $, $ 0<\Lambda_0<\Lambda_1 $, $ \Lambda_0\geq\Lambda_1 $ and $ -|\Lambda_1|<\Lambda_0<0 $. The figures are taken from Ref. [29].

          Apart from bouncing and cyclic scenarios, one should note that a non-singular Universe can also be realized via emergent behavior, when the Universe approaches a static model in the infinite past [149151]. In this case, $ H\rightarrow 0 $ and $ \rho\rightarrow 0 $, while in order to enter into a realistic expanding phase, H still needs to gradually increase, with $ \dot H>0 $. Thus the EoS of the Universe will go from $ -\infty $ to above -1 as well. Therefore, in the emergent scenario, the Universe also realizes quintom-like behavior [152, 153].

          In conclusion, in this paper we briefly reviewed the Quintom model of dark energy and its applications to constructing bouncing cosmologies. Specifically, we have considered three classes of quintom scenarios: a two-field model, a single scalar with higher derivatives and modified gravity.

          Recently, other non-singular bounce scenarios making use of matter with EoS crossing -1 have been proposed. In Ref. [154] Bardeen, Bars, Hanson and Peccei (BBHP) introduced the folded string, which is a classical solution of a limit of the Nambu model. This was generalized by Itzhaki [155], who found that in a classical (1+1)-dimensional string model with an increasing dilaton field, there are Instant Folded String (IFS) solutions. These are essentially the BBHP folded strings, with the roles of time and space reversed so that they spontaneously nucleate 5. The IFS violates the null energy condition (NEC) [158], a fact which was used to create a non-singular bouncing cosmology in Ref. [159]. In Ref. [160], Itzhaki, Peleg and Steinhardt showed that by introducing a particular potential for the dilaton, together with an extra field that becomes massless at an enhanced symmetry point [161], the bounce of Ref. [159] can be incorporated into a cyclic bouncing cosmology, complete with reheating performed by the extra field.

        D.   Beyond a bounce: Cyclic scenarios
        • A simple and precise example of a quintom cyclic model is given by the two-field action [29]:

          $ S=\int {\rm d}^4x\sqrt{-g}\left[\frac{1}{2}\partial_\mu\phi\partial^\mu\phi-\frac{1}{2}\partial_\mu\psi\partial^\mu\psi-V(\phi,\psi)\right]\; , $

          (82)

          Thus the energy density and pressure of the universe are

          $ \rho=\frac{1}{2}\dot\phi^2-\frac{1}{2}\dot\psi^2+V(\phi,\psi)\; ,\; \; \; p=\frac{1}{2}\dot\phi^2-\frac{1}{2}\dot\psi^2-V(\phi,\psi)\; , $

          (83)

          while the equations of motion for the ϕ and ψ fields are

          $ \ddot\phi+3H\dot\phi+V_{,\phi}=0\; ,\; \; \; \ddot\psi+3H\dot\psi-V_{,\psi}=0\; . $

          (84)

          To realize a cyclic scenario, we consider a coupled potential of the form

          $ V(\phi,\psi)=(\Lambda_0+\lambda\phi\psi)^2+\frac{1}{2}m^2\phi^2-\frac{1}{2}m^2\psi^2\; . $

          (85)

          Due to the symmetry of the potential under the transformation: $ \phi\rightarrow i\psi $, $ \psi\rightarrow -i\phi $, the two scalar fields will dominate alternately. Due to the couplings between two fields, the potential will be bounded from below. For Eq. (84) (with Hubble parameter satisfying the Friedmann equation and the energy density comes from Eq. (83)), we can obtain an Ansatz solution

          $ \phi=\sqrt{A_0}\cos mt\; ,\; \; \; \psi=\sqrt{A_0}\sin mt\; , $

          (86)

          where $ A_0 $ is the oscillation amplitude of the fields. Moreover, the Hubble parameter turns out to be:

          $ H=\frac{\sqrt{3}}{3M_p}(\Lambda_0+\Lambda_1\sin 2mt)\; , $

          (87)

          where we define $ \Lambda_1=\sqrt{3}mA_0/(4M_p) $.

          From Eq. (87) one can see that the Hubble parameter shows an oscillating behavior, and since $ H=d\ln a/{\rm d}t $, the scale factor can be oscillating as well. Moreover, different choices of the parameters $ \Lambda_0 $ and $ \Lambda_1 $ can give rise to different oscillating behaviors. In Fig. 13 we provide 4 possible kinds of oscillating behaviors of $ \ln a $, H and EoS w given by Eq. (87). The first case is for $ \Lambda_0=0 $ where both H and $ \ln a $ are oscillating with an equal amplitude. This is a standard oscillating Universe, but it suffers from the entropy problem. The second case is for $ 0<\Lambda_0<\Lambda_1 $. In this case, the center value of the Hubble parameter is positive, so the period of $ H<0 $ in each cycle will be shorter than that of $ H>0 $, therefore, the $ \ln a $ is oscillating with the total amplitude increasing. The third case is for $ \Lambda_0\geqslant \Lambda_1 $. In this case the Hubble parameter is always larger than zero, so the Universe is always expanding with increasing $ \ln a $. However, the expansion will be accelerating or decelerating alternately. Since there is no transition process from contraction to expansion, the EoS will also be regular all the time. The last case is for $ -|\Lambda_1|<\Lambda_0<0 $ where the center value of Hubble parameter is negative. In this case, period of $ H<0 $ in each cycle will be longer than that of $ H>0 $, therefore, the $ \ln a $ will be decreasing in total, but have short periods of increase in each cycle. Therefore, it still has room to explain our current expanding universe.

          Figure 13.  (color online) The evolution of $ \ln a $, H, ρ and w with respect to cosmic time t, corresponding to blue, purple, red and green lines, respectively. In the numerical calculation, the parameters are taken to be $ m=3\times 10^{-3} $ and $ \Lambda_1\simeq 0.39 $, and $ \Lambda_0 $ is taken as 0 (left-top), $ 0.10 $ (right-top), $ 0.402 $ (left-bottom) and $ -0.10 $ (right-bottom) respectively, corresponding to the cases of $ \Lambda_0=0 $, $ 0<\Lambda_0<\Lambda_1 $, $ \Lambda_0\geq\Lambda_1 $ and $ -|\Lambda_1|<\Lambda_0<0 $. The figures are taken from Ref. [29].

          Apart from bouncing and cyclic scenarios, one should note that a non-singular Universe can also be realized via emergent behavior, when the Universe approaches a static model in the infinite past [149151]. In this case, $ H\rightarrow 0 $ and $ \rho\rightarrow 0 $, while in order to enter into a realistic expanding phase, H still needs to gradually increase, with $ \dot H>0 $. Thus the EoS of the Universe will go from $ -\infty $ to above -1 as well. Therefore, in the emergent scenario, the Universe also realizes quintom-like behavior [152, 153].

          In conclusion, in this paper we briefly reviewed the Quintom model of dark energy and its applications to constructing bouncing cosmologies. Specifically, we have considered three classes of quintom scenarios: a two-field model, a single scalar with higher derivatives and modified gravity.

          Recently, other non-singular bounce scenarios making use of matter with EoS crossing -1 have been proposed. In Ref. [154] Bardeen, Bars, Hanson and Peccei (BBHP) introduced the folded string, which is a classical solution of a limit of the Nambu model. This was generalized by Itzhaki [155], who found that in a classical (1+1)-dimensional string model with an increasing dilaton field, there are Instant Folded String (IFS) solutions. These are essentially the BBHP folded strings, with the roles of time and space reversed so that they spontaneously nucleate 5. The IFS violates the null energy condition (NEC) [158], a fact which was used to create a non-singular bouncing cosmology in Ref. [159]. In Ref. [160], Itzhaki, Peleg and Steinhardt showed that by introducing a particular potential for the dilaton, together with an extra field that becomes massless at an enhanced symmetry point [161], the bounce of Ref. [159] can be incorporated into a cyclic bouncing cosmology, complete with reheating performed by the extra field.

      Reference (161)

目录

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return