---
abstract: 'We study the evolution of cosmological perturbations in a non-singular bouncing cosmology with a bounce phase which has superimposed oscillations of the scale factor. We identify length scales for which the final spectrum of fluctuations obtains imprints of the non-trivial bounce dynamics. These imprints in the spectrum are manifested in the form of damped oscillation features at scales smaller than a characteristic value and an increased reddening of the spectrum at all the scales as the number of small bounces increases.'
author:
- Robert Brandenberger
- Qiuyue Liang
- 'Rudnei O. Ramos'
- Siyi Zhou
title: Fluctuations through a Vibrating Bounce
---

Introduction
============

Some bouncing cosmologies provide an alternative to cosmological inflation as a way to obtain primordial cosmological fluctuations (see, e.g., Ref. [@BP] for a recent review). Specifically, in a model which contains a matter-dominated phase of contraction, initial vacuum fluctuations in the far past which exit the Hubble radius during the matter-dominated contracting phase evolve into a scale-invariant spectrum of cosmological perturbations [@FB; @Wands]. For example, in the case of an Ekpyrotic contracting universe [@Ekp] entropy fields can source scale-invariant curvature fluctuations [@newEkp]. In all bouncing cosmologies, new physics is required in order to obtain a non-singular cosmological bounce. Such new physics could come from the matter sector (see, e.g., Refs. [@matterBounce-1; @matterBounce-2]), from modifications of the classical gravitational action (as for example in Horava-Lifshitz gravity [@HLbounce] or in the non-local gravity construction of Ref. [@Biswas]), or from quantum gravity effects. Examples of the latter are the bounce in loop quantum cosmology (see, e.g., Ref. [@LQCbounce] for reviews), in deformed AdS/CFT cosmology [@Elisa], the S-brane bounce of Ref. [@Costas] and the temperature bounce in String Gas Cosmology [@BV].

Concerning the robustness of the computations of the spectrum of cosmological fluctuations, an advantage of bouncing cosmologies (without an inflationary phase after the bounce) is that the physical length of modes which are probed in current observations remain in the far infrared throughout the cosmological evolution as long as the energy density at the bounce point is smaller than the Planck density. Hence, the computations can be done in the realm where effective field theory is well justified. This is in contrast to the situation in inflationary cosmology [@MB] where the physical wavelengths of even the largest scales which are currently observed are smaller than the Planck length at the beginning of inflation (provided that the inflationary phase lasts slightly longer than the minimal period it has to last in order to solve the horizon and flatness problems of Standard Big Bang Cosmology).

A key question is to whether the predictions for cosmological perturbations at late times in the expanding phase are sensitive to the details of the bounce phase. For simple parametrizations of the bounce phase, detailed studies have shown that the spectral shape does not change during the bounce phase provided that the duration of the bounce phase is shorter than the length scale of the fluctuations at the bounce point (see, e.g., Ref. [@matterBounce-2] in the case of matter-driven bounces, Ref. [@HLflucts] in the case of the Horava-Lifshitz bounce, Ref. [@Elisa] in the case of the AdS/CFT bounce, and Ref. [@Subodh] for the S-brane bounce). On the other hand, there are examples where the bounce phase yields dramatic changes in the spectrum [@BXue]. The reason why such dramatic changes are possible is that the Hubble radius at the bounce point is infinite, and we cannot invoke the freezing of cosmological perturbations on super-Hubble scales to argue for a constancy of the spectrum[^1].

To further analyze the sensitivity of the spectrum of cosmological fluctuations on the details of the bounce phase, we here consider a toy model where the scale factor undergoes small amplitude oscillations during the bounce phase. Such a behavior may emerge from certain models motivated by ideas from loop quantum gravity [@Alesci]. Heuristically, one would argue that those small bounces will not influence the large scale modes provided that the wavelengths of these modes are so large scale that they would not feel the small scale fluctuations of the scale factor. On the other hand, smaller scale modes whose wavelength is comparable or smaller to the total duration of the bounce phase should be sensitive to the details of the dynamics during the bounce. In this work, we would like to give a careful treatment to see if this is really the case.

In the case of a cyclic cosmology, when the time interval between cycles is larger than the wavelength of the modes being considered, it is generally sufficient [@RHBcyclic] to consider only the dominant modes in each phase (except during the bounce phase and when mode matching conditions are applied [@Durrer]). In our case however, the time scale between the small bounces is very small compared to the length scales of interest, and hence we cannot just focus on the dominant modes because for those small time durations the subdominant modes can also have an effect on the primordial power spectra. We need to keep all the contributions and give a comprehensive analysis.

This paper is organized as follows. In Sec. \[sec2\], we specify our setup for the intermediate small bounce feature and discuss about the relevant scales involved. In Sec. \[sec3\], we present the calculation of power spectrum. A detailed presentation of the required matching conditions is given and the specific results are given for the two specific small inter bounce features we considered. An analysis of these results for the power spectrum is then given in Sec. \[sec4\]. In Sec. \[sec5\], we give a generalization for the case of a large number of small bounces. In Sec. \[sec6\] we present our conclusions. An appendix is included to discuss some of the technical details.

Setup {#sec2}
=====

We will consider a spatially flat Friedmann-Lemaitre-Robertson-Walker space-time in which the metric is given by $$ds^2 \, = \, dt^2 - a(t)^2 d{\bf x}^2 \, ,$$ where $t$ is physical time, ${\bf x}$ are the comoving spatial coordinates, and $a(t)$ is the cosmological scale factor. It will be convenient to use conformal time $\tau$ related to the physical time via $dt = a(t) d\tau$. We will consider only linear cosmological perturbations (see, e.g., Ref. [@MFB] for a comprehensive review). In this case, fluctuations evolve independently in Fourier space. We will label the fluctuation modes in terms of their comoving wave number $k$.

We consider a non-singular symmetric bouncing cosmology in which the cosmological scale factor $a(\tau)$ has the form shown in Fig. \[bounceplot\], i.e., for which $a(\tau)$ has one “oscillation” between the onset and end of the bounce phase. We consider the two forms depicted in Fig. \[bounceplot\], with one peak between the time $- \tau_B$ when $a(\tau)$ reaches its first minimum, and the time $+ \tau_B$ when the second minimum of $a(\tau)$ is taken on. Specifically, we consider two specific models: Model 1, with a flat plateau about $\tau = 0$; and the Model 2, with a kink of $a(\tau)$ at $\tau = 0$, which is the limiting case of the first model when the duration of the flat plateau equals to zero. The time interval of the plateau region is $-\Delta < \tau < \Delta$ given by some conformal time $\Delta$ with $\Delta < \tau_B$. The forms shown in Fig. \[bounceplot\] are simpler enough such as to allow an analytical study, but already of sufficient complexity such as to provide the main relevant features in the power spectrum that we might also observe in some more complex setup for the inter-bounce features. In particular, our results can be easily generalized to the case of many oscillations, as we will later discuss in Sec. \[sec5\].

As seen in Fig. \[bounceplot\], for the two forms shown the time interval can be divided into five intervals. The first is the initial contracting phase (Phase I) $\tau < - \tau_B$. The second is the intermediate expanding phase (Phase II). The third is Phase III with static scale factor, the fourth (Phase IV) is the intermediate contracting phase and Phase V is the final expanding phase. Fluctuation modes exit the Hubble radius in the initial contracting phase. The top panels in Fig. \[bounceplot\] give a sketch of the scale factor, the lower panels show the corresponding time evolution of the comoving Hubble radius. The vertical axis of the lower panels can also be viewed as a label for comoving wavelength. In this way, it is easy to read off when various modes enter and exit the Hubble radius.

All scales re-enter and re-exit the Hubble radius several times since at the extrema of $a(t)$ the Hubble radius is infinite. We will treat the transitions at $-\tau_B$ and $\tau_B$ as instantaneous[^2]. More specifically, we will cut out a time interval $$-\tau_B - \epsilon < \tau < - \tau_B + \epsilon,$$ (with $\epsilon \ll \tau_B$) and correspondingly another time interval of the same length about $\tau_B$ and we will match the solutions between the neighboring phases making use of the matching conditions given in Refs. [@HV; @DM], which are the cosmological version of the Israel [@Israel] ones[^3]. Hence, the only Hubble re-entry which is important to us is the one which occurs between $-
\tau_B$ and $+ \tau_B$.

In the first model, given by the plots on the left shown in Fig. \[bounceplot\], there are two characteristic comoving length scales. The first is $k_*^{-1}$ which is defined as the length which re-enters the Hubble radius at time $- \Delta$. The second one, $k_{\Delta}^{-1}$, is the mode which undergoes one oscillation between $\tau = - \Delta$ and $\tau = + \Delta$ (we are assuming $k_{\Delta}^{-1} >  k_*^{-1}$ if this is not satisfied then we recover the results for the second model, given by the plots on the right shown in Fig. \[bounceplot\]). Modes with wavelength smaller than $k_*^{-1}$ enter the Hubble radius during Phase II and exit again during Phase IV. Modes with $k_*^{-1} < k^{-1}
< k_{\Delta}^{-1}$ are inside the Hubble radius only during Phase III. For these modes the matching occurs at times $- \Delta$ and $+
\Delta$. This is also true for modes with $k^{-1} > k_{\Delta}^{-1}$. These modes, however, undergo a negligible amount of oscillations in Phase III. In the case that $k_\Delta^{-1}>k_*^{-1}$, we have three different behaviors of the power spectrum. For the very large scale modes $k^{-1}\gg k_\Delta^{-1}$, the power spectrum does not feel the influence of the small bump of the scale factor. For the modes $k_*^{-1}<k^{-1}<k_\Delta^{-1}$, there is a complicated change of the power spectrum induced by the flat plateau. For the modes $k^{-1}<k_*^{-1}$, the change of the power spectrum approaches the well known result for cyclic cosmologies [@RHBcyclic], as we will explicitly verify later on below. In the case that $k_\Delta^{-1}<k_*^{-1}$, there is only one characteristic scale $k_*^{-1}$. The mode with $k^{-1}>k_*^{-1}$ will not feel the influence of the bump, while the mode with $k^{-1}<k_*^{-1}$ will be changed by the bump according to the well known result for cyclic cosmologies. A special situation belonging to this case is the limiting case $\Delta \rightarrow 0$.

We divide the evolution of fluctuation modes into five regions as shown on Fig. \[bounceplot\] [^4]. The five regions are denoted by Region I, Region II, Region III, Region IV and Region V, respectively. Region I and Region II are separated by the time $-\tau_B$. Region IV and Region V are separated by $\tau_B$. The separation between Regions II and III, and between Regions III and IV are more complicated. Because of the existence of the flat plateau (or of the local maximum of the comoving Hubble radius in the case of Model 2), we can see that there is a clear distinction between the large scale and small scale modes separated by a characteristic scale $k_*^{-1}$. Small scale modes (i.e. $k^{-1}<k_*^{-1}$) enter the Hubble radius at time $-\tau_H(k)
\leq - \Delta$ and exit the Hubble radius at $\tau_H(k) > \Delta$, and the separation between Regions II and III and between Regions III and IV are given by $-\tau_H(k)$ and $\tau_H(k)$, respectively. For large scale modes, the separations between Regions II and III, and between Regions III and IV are given by the times $-\Delta$ and $\Delta$, respectively, because the modes enter and exit the Hubble radius at these two times. The situation in the case of Model 2 is simpler. The evolution of small scale modes is the same as the case with a flat plateau, while the large scale modes have only four regions which we denote by Regions I, II, IV and V, respectively. The separation between Regions II and IV, in this case, is the time $\tau= 0$.

The Computation of the Power Spectrum {#sec3}
=====================================

We are interested in the power spectrum of the primordial curvature perturbation $\zeta$ (see, e.g., Ref. [@MFB] for a review of the theory of cosmological perturbations). We quantize the linear fluctuations and write them in terms of the more convenient Mukhanov-Sasaki variable $v$. In the case of a constant equation of state, the relation between $\zeta$ and $v$ is $$\begin{aligned}
v = C a \zeta \, ,\end{aligned}$$ where $C$ is a constant. Thus, the equation of motion for the mode function $v$, in momentum space, is given by $$\begin{aligned}
v'' + \bigg(k^2 - \frac{a''}{a} \bigg) v = 0 \, .\end{aligned}$$ For the scale factor $a\sim \tau^q$, the solution of the Mukhanov-Sasaki equation is given by $$\begin{aligned}
v (\tau) = c_1(k) \sqrt{\tau} J_{\alpha} (k\tau) +  c_2(k) \sqrt{\tau}
Y_{\alpha} (k\tau),\quad \alpha\equiv\ q-\frac{1}{2}  ,\end{aligned}$$ where $J(x)$ and $Y(x)$ are Bessel functions of the first and second kind, respectively. On sub-Hubble scales, the solutions are oscillatory, on super-Hubble scales they can be approximated by a power law. To see this, we note that the expansion of the Bessel function solutions for small argument, $x\ll 1$, is $$\begin{aligned}
J_\alpha (x)  &= \sum_{m = 0}^\infty
\frac{(-1)^m}{m!\Gamma(m+\alpha+1)} \left( \frac{x}{2}
\right)^{2m+\alpha}, \\ Y_\alpha (x)  &=
\frac{\cos(\alpha\pi)}{\sin(\alpha\pi)} J_\alpha(x) -
\frac{1}{\sin(\alpha\pi)}J_{-\alpha}(x)  \nonumber \\ &=
\frac{\cos(\alpha \pi)}{\sin(\alpha \pi)} \sum_{m = 0}^\infty
\frac{(-1)^m}{m!\Gamma(m+\alpha+1)} \left( \frac{x}{2}
\right)^{2m+\alpha} - \frac{1}{\sin(\alpha\pi)}\sum_{m = 0}^\infty
\frac{(-1)^m}{m!\Gamma(m-\alpha+1)} \left( \frac{x}{2}
\right)^{2m-\alpha},\end{aligned}$$ we can express the mode function in terms of a series expansion $$\begin{aligned}
v(\tau) = \sum_{m=0}^{\infty} d_{1m} (k) \tau^{q+2m} +
\sum_{m=0}^{\infty} d_{2m} (k) \tau^{1-q+2m}  ,
\label{vtaum}\end{aligned}$$ where $d_{1m}(k)$ and $d_{2m}(k)$ are given, respectively, by $$\begin{aligned}
d_{1m}(k) & = \left[ c_1(k) + \frac{\cos(\alpha\pi)}{\sin(\alpha\pi)}
  c_2(k) \right]  \frac{(-1)^m}{m!\Gamma(m+\alpha+1)} \left(
\frac{k}{2} \right)^{2m+\alpha}, \\ d_{2m}(k) & =
-\frac{c_2(k)}{\sin(\alpha\pi)} \frac{(-1)^m}{m!\Gamma(m-\alpha+1)}
\left( \frac{k}{2} \right)^{2m-\alpha}  .\end{aligned}$$ Since we are interested in those modes that went classical (crossed the Hubble radius), such that $k \tau \ll 1$, the higher order terms of $v(\tau)$ given by $m> 0$ are subleading. Thus, in the following, we can just focus on the $m=0$ terms in Eq. (\[vtaum\]).

The scale factors and the solutions to the mode functions of the five regions can be obtained by shifting the time coordinate. Thus, they are given, respectively, by

-   Region I (contracting), where $\tau<-\tau_B$, we have that

    $$\begin{aligned}
    a \sim  (-\tau-\tau_B)^{q_2} \sim(-t-t_B)^{p_2}, \quad v_{1} = c_{11}
    (-\tau_B-\tau)^{1-q_2} + c_{12} (-\tau_B-\tau)^{q_2},\end{aligned}$$

-   Region II (expanding), where $ -\tau_B<\tau<-\Delta\,\,\, {\rm
      for}\,\,\, k^{-1}>k_*^{-1},   \,\,\, {\rm
      and}\,\,\,-\tau_B<\tau<-\tau_H \,\,\, {\rm for}\,\,\,
    k^{-1}<k_*^{-1}$, we have that

    $$\begin{aligned}
    a\sim (\tau+\tau_B)^{q_1}\sim(t+t_B)^{p_1},\quad v_2 =
    c_{21}(\tau+\tau_B)^{1-q_1} + c_{22} (\tau+\tau_B)^{q_1},\end{aligned}$$

-   Region III (intermediate), where $-\Delta<\tau<\Delta\,\,\, {\rm
      for}\,\,\, k^{-1}>k_*^{-1}, \,\,\, {\rm and}\,\,\,
    -\tau_H<\tau<\tau_H \,\,\, {\rm for}\,\,\, k^{-1}<k_*^{-1}$, we have that

    $$\begin{aligned}
    a\sim {\rm constant\,\,\,in\,\,\,-\Delta<\tau<\Delta} , \quad v_3 =
    c_{31} e^{i k \tau} + c_{32} e^{-i k \tau},\end{aligned}$$

-   Region IV (contracting), where $\Delta<\tau<\tau_B\,\,\, {\rm
      for}\,\,\, k^{-1}>k_*^{-1},   \,\,\, {\rm and}\,\,\,
    \tau_H<\tau<\tau_B \,\,\, {\rm for}\,\,\, k^{-1}<k_*^{-1}$, we have that

    $$\begin{aligned}
    a \sim  (-\tau+\tau_B)^{q_1} \sim (-t+t_B)^{p_1},\quad v_4 = c_{41}
    (-\tau+\tau_B)^{1-q_1} +c_{42} (-\tau+\tau_B)^{q_1},\end{aligned}$$

-   Region V (expanding), where $\tau>\tau_B$, we have that

    $$\begin{aligned}
    a\sim  (\tau-\tau_B)^{q_2} \sim (t-t_B)^{p_2},\quad v_5 = c_{51} (\tau
    - \tau_B)^{1-q_2} + c_{52} (\tau-\tau_B)^{q_2}  .\end{aligned}$$

Note that in the above, the time $\tau_H$ depends on $k$, $\tau_H\equiv \tau_H(k)$. To simplify the notation we do not write the k-dependence explicitly. Note also that $q_2$ (the power of the scale factor in conformal time) and $p_2$ (the power of the scale factor in cosmological time) are the indices of the scale factor during the initial contracting and final expanding phase, and that $q_1$ and $p_1$ are the indices in the intervening periods. In particular, the index 1 in these quantities should not be confused with the index in Region I. Note also that the $q_i$ are related to the $p_i$ through the relation

$$q_i = \frac{p_i}{1-p_i},\;\;\;i=1,2.
\label{qipi}$$

Matching conditions
-------------------

Let us now discuss the process of matching between the different regions discussed above. The matching conditions for the metric across a space-like hypersurface were derived in Hwang-Vishniac [@HV] and Deruelle-Mukhanov [@DM] and are the generalization of the Israel matching conditions [@Israel]. For cosmological fluctuations, the matching conditions say the solutions in different regions can be connected by enforcing two conditions, namely the continuity of both $v$ and its derivative across the boundary surface.

As mentioned earlier, for the two non-singular bouncing points $-
\tau_B$ and $\tau_B$ we match the solutions at times $\mp\tau_B-\epsilon$ and $\mp\tau_B+\epsilon$, neglecting any evolution in the intervening time period. This is similar to what was done in Refs. [@Elisa; @Subodh]. A second justification of this method (in addition to the one given earlier) is that for modes we are interested in, the time interval $2\epsilon$ is so small, thus the mode functions do not have enough time to oscillate inside the Hubble radius. On the other hand, in our first model (with a flat plateau for $a(t)$, model 1) we consider the case that the interval $2\Delta$ is sufficiently long such that some of the modes we are interested in have time to oscillate while the mode is inside the Hubble radius. Very large scale modes, on the other hand, still do not oscillate inside the Hubble radius.

### Between region I and region II

The matching conditions are $$\begin{aligned}
v_1(-\tau_B-\epsilon) = v_2 (-\tau_B+\epsilon), \quad
v'_1(-\tau_B-\epsilon) = v'_2 (-\tau_B+\epsilon) \, .\end{aligned}$$ We write down the equations explicitly in terms of coefficients $c_{ij}$ of the fundamental solutions of the equation of motion. The index $i$ stands for the region, the index $j$ (either $1$ or $2$) running over the two different modes:

$$\begin{aligned}
\begin{pmatrix}
\epsilon^{1-q_2} & \epsilon^{q_2} \\ -(1-q_2)  \epsilon^{-q_2} &
-q_2\epsilon^{q_2-1}
\end{pmatrix} \begin{pmatrix}
c_{11} \\ c_{12}
\end{pmatrix} = \begin{pmatrix}
\epsilon^{1-q_1} & \epsilon^{q_1} \\ (1-q_1) \epsilon^{-q_1} & q_1
\epsilon^{q_1-1}
\end{pmatrix} \begin{pmatrix}
c_{21} \\ c_{22}
\end{pmatrix} .\end{aligned}$$

### Between region II and region III

To be explicit, we focus in this case on the large scale modes for which the time duration of Region III is $-\Delta<\tau<\Delta$. For the small scale modes which enter the Hubble radius before $- \Delta$, we just make the substitution $\Delta\rightarrow\tau_H(k)$. Apart from that the discussion is the same. For the next subsection the convention will be the same. The matching conditions in this case are $$\begin{aligned}
v_2(-\Delta) = v_3 (-\Delta),\quad  v'_2(-\Delta) = v'_3 (-\Delta),\end{aligned}$$ which in matrix form can be expressed as

$$\begin{aligned}
\begin{pmatrix}
\delta^{1-q_1} & \delta^{q_1} \\ (1-q_1) \delta^{-q_1} & q_1
\delta^{q_1-1} 
\end{pmatrix} \begin{pmatrix}
c_{21} \\ c_{22}
\end{pmatrix} = \begin{pmatrix}
e^{-ik \Delta} & e^{i k \Delta} \\ i k e^{-ik\Delta} & - i k e^{i k
  \Delta}
\end{pmatrix} \begin{pmatrix}
c_{31} \\ c_{32}
\end{pmatrix},\end{aligned}$$

and where in the above equation we have defined $\delta$ as

$$\delta = \tau_B - \Delta .
\label{delta}$$

Note that for the small scale modes the definition of $\delta$ should be changed to $\delta\to \delta(k) \equiv \tau_B-\tau_H(k)$.

### Between region III and region IV

The matching conditions are $$\begin{aligned}
v_3(\Delta) = v_4 (\Delta), \quad v'_3(\Delta) = v'_4 (\Delta) \, .\end{aligned}$$ In matrix form this yields $$\begin{aligned}
\begin{pmatrix} 
e^{ik\Delta} & e^{-ik\Delta} \\ i k e^{ik\Delta} & -i k e^{-ik\Delta}  
\end{pmatrix} \begin{pmatrix}
c_{31} \\ c_{32}
\end{pmatrix} = \begin{pmatrix}
\delta^{1-q_1} & \delta^{q_1} \\ -(1-q_1) \delta^{-q_1} &
-q_1\delta^{q_1-1}
\end{pmatrix} \begin{pmatrix}
c_{41} \\ c_{42}
\end{pmatrix}.\end{aligned}$$

### Between region IV and region V

The matching conditions in this case are $$\begin{aligned}
v_4(\tau_B-\epsilon) = v_5 (\tau_B+\epsilon),\quad
v'_4(\tau_B-\epsilon) = v'_5 (\tau_B+\epsilon) \, .\end{aligned}$$ In matrix form this yields $$\begin{aligned}
\begin{pmatrix}
\epsilon^{1-q_1} & \epsilon^{q_1} \\ -(1-q_1) \epsilon^{-q_1} & -q_1
\epsilon^{q_1-1}
\end{pmatrix} \begin{pmatrix}
c_{41} \\ c_{42}
\end{pmatrix} = \begin{pmatrix}
\epsilon^{1-q_2} & \epsilon^{q_2} \\ (1-q_2) \epsilon^{-q_2} & q_2
\epsilon^{q_2-1}
\end{pmatrix} \begin{pmatrix}
c_{51} \\ c_{52}
\end{pmatrix}.\end{aligned}$$

Analysis and Results for the Power Spectrum {#sec4}
===========================================

Combining the results of the previous section we find that the final mode coefficients can be written in terms of the initial ones via $$\begin{aligned}
\nonumber
& \mathcal C_5 = 
\begin{pmatrix}
c_{51} \\ c_{52}
\end{pmatrix},\end{aligned}$$ where[^5]

$$\begin{aligned}
c_{51}&=&\frac{c_{11}}{(1-2q_2)(1-2q_1)2k(1-2q_1)} \left[ a_{11} (1 -
  2 (q_1-1) q_1 - 2 (q_2 -1 ) q_2) \right.  \nonumber \\ &-&
  \left. a_{12} ( (q_1 - q_2 - 1) (q_1 + q_2 - 2) \epsilon^{1 - 2 q_1}
  ) + a_{21} (q_1 - q_2 + 1) (q_1 + q_2) \epsilon^{2 q_1 - 1} \right],
\label{c51}
\\ c_{52}&=&\frac{c_{11}}{(1-2q_2)(1-2q_1)2k(1-2q_1)} \left[ 2 a_{11}
  (1 + q_1 - q_2) (-2 + q_1 + q_2) \epsilon^{1 - 2 q_2} \right.
  \nonumber \\ &+& \left. a_{12} (-2 + q_1 + q_2)^2 \epsilon^{-2 (-1 +
    q_1 + q_2)} + a_{21} (1 + q_1 - q_2)^2 \epsilon^{2 (q_1 -
    q_2)}\right],
\label{c52}\end{aligned}$$

and $$\begin{aligned}
\label{a11}
a_{11} & = \frac{2\left[k^2\delta^2+(q_1-1)q_1\right]\sin
  (2k\Delta)}{\delta} - 2 k \cos (2 k \Delta), \\
\label{a12}
a_{12} & = \delta^{2q_1-2}
\left[-4kq_1\delta\cos(2k\Delta)-2(q_1-k\delta)(q_1+k\delta)
  \sin(2k\Delta)\right],
\\
\label{a21}
a_{21} & =
\delta^{-2q_1}\left[2(q_1-k\delta-1)(q_1+k\delta-1)\sin(2k\Delta)-
  4k(q_1-1)\delta\cos(2k\Delta)\right]
.\end{aligned}$$ Note that these coefficients oscillate as a function of $k$. These oscillations are important, however, only for small wavelength fluctuations. For these we will obtain oscillations in the power spectrum. The final general result for the power spectrum is given by $$\begin{aligned}
P_{\zeta} & = \zeta^2 k^{3}  \sim \left( \frac{v}{z} \right)^2 k^{3}
\nonumber  \\ & = \left[
  \frac{c_{51}(\tau-\tau_B)^{1-q_2}+c_{52}(\tau-\tau_B)^{q_2}}{(\tau-\tau_B)^{q_2}}
  \right]^2 k^{3}    \nonumber \\  &= \left[ c_{51}
  (\tau-\tau_B)^{1-2q_2} + c_{52} \right]^2 k^{3} \, .
\label{finalPzeta}\end{aligned}$$ Below we will analyze some of the specific cases given by our two models when applying the result given by Eq. (\[finalPzeta\]).

Limiting case of Instantaneous matching
---------------------------------------

We first consider the limit as the duration of the plateau region of $a(t)$ goes to zero, corresponding to what we have denoted by Model 2 in Sec. \[sec2\]. This is the limit $\Delta\rightarrow 0$. In this case, large scale modes $k^{-1} > k_*^{-1}$ do not enter the Hubble radius in the region near $t = 0$, and we can set $\Delta = 0$ in the matching condition equations, i.e.,

$$\begin{aligned}
\sin(2k\Delta) \rightarrow 0, \quad \cos(2 k \Delta) \rightarrow 1 \,
.\end{aligned}$$

On the other hand, small scale modes $k^{-1} < k_*^{-1}$ will enter the Hubble radius at a time given by $-\tau_H(k)$, and in the matching condition equations we must replace $\Delta$ by $\tau_H(k)$.

### Large scale modes $k^{-1} > k_*^{-1}$

Let us first consider the case for large scale modes $k^{-1} >
k_*^{-1}$. In this case we have that $$\begin{aligned}
a_{11}\rightarrow -2k, \quad a_{12}\rightarrow -4k q_1\delta^{2q_1-1},
\quad a_{21}\rightarrow -4k (q_1-1)\delta^{-2q_1+1}.\end{aligned}$$ Because we are interested in the parameter region $1/3<p<1$, then, written in terms of $q$, we have $q>1/2$. So the $c_{51}$ mode in the expression for the power spectrum Eq. (\[finalPzeta\]) is a decaying solution. Hence, we can focus on the constant mode $c_{52}$, and thus the power spectrum in this case becomes $$\begin{aligned}
 \nonumber
&P_{\zeta} \sim c_{52}^2 k^{3} =\left\{
  \frac{c_{11}}{(1-2q_2)(1-2q_1)2(1-2q_1)}  \right. \\ \nonumber &
  \left.  \times \left[ -4 (1 + q_1 - q_2) (-2 + q_1 + q_2)
    \epsilon^{1 - 2 q_2}-4q_1\delta^{2q_1-1}  (-2 + q_1 + q_2)^2
    \epsilon^{-2 (-1 + q_1 + q_2)} \right. \right. \nonumber \\   &
    \left. \left. -4(q_1-1)\delta^{-2q_1+1} (1 + q_1 - q_2)^2
    \epsilon^{2 (q_1 - q_2) } \right] \right\}^2 k^{3}  .\end{aligned}$$

The initial power spectrum is $$\begin{aligned}
P_i = P_{\zeta} (-\tau_B-\epsilon) = \zeta^2 k^{3}  = c_{11}^2
\epsilon^{2-4q_2}  k^{3}\end{aligned}$$ and, thus, we can relate the final to the initial power spectrum as $$\begin{aligned}
\label{powerspectrum}
P_{\zeta} = \left(A_1  + A_2   \delta^{2q_1-1} \epsilon^{1-2 q_1 } +
A_3\delta^{-2q_1+1} \epsilon^{2q_1-1} \right)^2 P_i \, ,\end{aligned}$$ where $A_1$, $A_2$ and $A_3$ are constants that do not depend on $k$. Their explicit forms are $$\begin{aligned}
&A_1 = \frac{-2(1+q_1-q_2)(-2+q_1+q_2)}{(1-2q_2)(1-2q_1)^2},\\  &A_2 =
  \frac{-2q_1(-2+q_1+q_2)^2}{(1-2q_2)(1-2q_1)^2}, \\ &A_3 =
  \frac{-2(q_1-1)(1+q_1-q_2)^2}{(1-2q_2)(1-2q_1)^2}.\end{aligned}$$

For very large scale modes $k^{-1} \gg  k_*^{-1}$, $\delta \to
\tau_B$ and $\delta$ can be approximated as a constant time interval. Thus, the power spectrum in this case becomes

$$\begin{aligned}
\label{powerspectrum1}
(A_1  + A_2  \tau_B^{2q_1-1} \epsilon^{1-2 q_1 } + A_3
  \tau_B^{-2q_1+1} \epsilon^{2q_1-1} )^2 P_i .\end{aligned}$$

The first conclusion we draw from this result is that the shape of the spectrum for large scale modes does not change during the bounce. This agrees with the conclusions of previous work on simple bounce models [@matterBounce-2]. The amplitude, on the other hand, is amplified. Recall that $2q_i - 1 > 0$, and that $\epsilon \ll \tau_B$. Hence, it is the second term in Eq. (\[powerspectrum1\]) which dominates, and we conclude that the amplitude of the spectrum is amplified by a factor of $${\cal A} \, = \, A_2^2 \bigg(\frac{\tau_B}{\epsilon}\bigg)^{4q_1 - 2} .
\label{Amplitude}$$ This result can also be understood easily: Fluctuations grow both in the contracting and in the expanding phase. In fact, the fluctuations diverge in the limit when the scale factor becomes zero. Hence, without an effective cutoff $\epsilon$ we would get a divergence in the spectrum. With a cutoff, the enhancement factor of the amplitude of the power spectrum will be determined by the dimensionless ratio between $\tau_B$ and $\epsilon$ to a power which depends on the growth rate of the fluctuations on super-Hubble scales, i.e., on $q_1$ (see the discussion of these issues in a more general context in the review article Ref. [@RHBbounceReview]).

### Small scale modes $k^{-1} < k_*^{-1}$

For small scale modes $k^{-1}<k_{*}^{-1}$, we set $\delta$ equal to the Hubble crossing time. Thus, we can use the Hubble crossing condition $a H= k$, which from $a\sim \tau^q$ gives $$\begin{aligned}
\delta =q_1  k^{-1} .\end{aligned}$$ But we need to have $k \tau_H = k\tau_B - k\delta= k\tau_B - q_1 $. As a consequence of the oscillations in the coefficients $a_{ij}$, Eqs. (\[a11\]), (\[a12\]) and (\[a21\]), the final power spectrum of fluctuations will oscillate for small wavelengths. This is explicitly manifested when we show a numerical example for the power spectrum in Fig. \[bounceplot1\], where we chose an initial pre-bounce spectrum which is scale-invariant. We see that the scale-invariance of the spectrum is maintained on large scales, but that on small length scales there is both a change in the slope of the spectrum, and superimposed oscillations.

In the following we discuss in what range we can reproduce the results of Ref. [@RHBcyclic], which hold for a cyclic cosmology. In that work, it was found that for modes which re-enter the Hubble radius during the bounce phase, the index of the spectrum of cosmological perturbations changes during each cycle. For a matter-dominated contracting phase the change in the index $n_s$ of the power spectrum was determined to be $\Delta n_s = - 2$.

The results of Ref. [@RHBcyclic] are applicable when $k^{-1}<k_*^{-1}$, but for quite large scales such that $k^{-1}\rightarrow k_*^{-1}$. In this range we have $k\tau_H\rightarrow 0$. Using this in Eq. (\[powerspectrum\]), we obtain that $$\begin{aligned}
\label{powerspectrum_cyclic}
P_f \, = \,  \left[ A_1  + A_2 (k \epsilon/q_1)^{1-2q_1} +  A_3
  (k\epsilon/q_1)^{2q_1-1} \right]^2 P_i.\end{aligned}$$ Since $2q_1 - 1 > 0$ it is the second term in Eq. (\[powerspectrum\_cyclic\]) which dominates. Hence, we conclude that there is a change in the index of the power spectrum by $$\Delta n_s = 2-4 q_1= - 2 \frac{3 p_1 - 1}{1 - p_1} ,
\label{tiltns}$$ which coincides with the results of Ref. [@RHBcyclic]. This is as expected because the case studied in Ref. [@RHBcyclic] corresponds to a big bounce where $\delta$ is (cosmologically) large.

Case with a flat plateau
------------------------

In the case with a flat plateau and when $\Delta$ is very small, we have just one characteristic comoving mometum. However, when $\Delta$ is big, we have two key comoving momenta which are characterized by the mode which cross the Hubble radius at $\Delta$ and $\tau_B-
\Delta$, respectively. In this subsection, we would like to analyze in detail these two cases.

First we would like to calculate the critical comoving momentum $k_*^{-1}$. We start by analyzing the Hubble parameter $H$. The corresponding comoving Hubble parameter in region II is

$$\begin{aligned}
aH =q_1 (\tau+\tau_B)^{-1}.\end{aligned}$$

The critical scale $k_{*}$, which is obtained by $k_*=aH(\tau=-\Delta)$, is therefore $$\begin{aligned}
k_{*} =  q_1 (\tau_B - \Delta)^{-1}.\end{aligned}$$ The analysis here is similar to the instantaneous matching case of the previous subsection and we can obtain the power spectrum as $$\begin{aligned}
&P_{\zeta} \sim c_{52}^2 k^{3} = \left\{ \frac{ 1
  }{(1-2q_2)(1-2q_1)2k(1-2q_1)}  \right.  \nonumber \\ & \left. \times
  \left[  2 a_{11} (1 + q_1 - q_2) (-2 + q_1 + q_2)   + a_{12} (-2 +
    q_1 + q_2)^2 \epsilon^{1 -2 q_1 } \right. \right.  \nonumber \\ &
    \left. \left. + a_{21} (1 + q_1 - q_2)^2 \epsilon^{2 q_1 -1}
    \right] \right\}^2 P_i ,
\label{powerspectrumforinstantaneous}\end{aligned}$$ where $a_{11}$, $a_{12}$ and $a_{21}$ were already defined by Eqs. , and , respectively.

Model with no Region III
------------------------

Let us here consider the model with no Region III. We expect that the result we obtained in the previous Subsection will approach the result derived here in the limit when $\Delta\rightarrow 0$. The matching condition of Region I and II, Region IV and V are completely the same as in the flat plateau case, so here we only write down the matching condition between Region II and IV:

$$\begin{aligned}
v_2(0) = v_4 (0),\quad v'_2(0) = v'_4 (0) ,\end{aligned}$$

which can be written in terms of the more convenient matrix form $$\begin{aligned}
\begin{pmatrix}
\tau_B^{1-q_1} & \tau_B^{q_1} \\ (1-q_1) \tau_B^{-q_1} & q_1
\tau_B^{q_1-1} 
\end{pmatrix} \begin{pmatrix}
c_{21} \\ c_{22}
\end{pmatrix} = \begin{pmatrix} 
\tau_B^{1-q_1} & \tau_B^{q_1} \\ -(1-q_1) \tau_B^{-q_1} & -q_1
\tau_B^{q_1-1}
\end{pmatrix} \begin{pmatrix}
c_{41} \\ c_{42}
\end{pmatrix}.\end{aligned}$$ Combining these matching results we obtain the power spectrum completely the same as that of the instantaneous matching of the previous section.

Numerical Examples
------------------

In Fig. \[bounceplot1\] we show the form of the spectrum of cosmological perturbations and its tilt as a function of the comoving wavenumber $k$ for the two models we have considered in the previous sections.

We can clearly identify in Fig. \[bounceplot1\](a) the characteristic scales for each of the two models we have defined in Sec. \[sec2\]. For the flat plateau model (Model 1), there are two relevant comoving scales,

$$\begin{aligned}
&&k_{*,{\rm flat}}= \frac{q_1}{\tau_B-\Delta},
\label{kflat}
\\ &&k_{*,{\rm osc}}= \frac{q_1}{\tau_B}.
\label{kosc}\end{aligned}$$

In model 1 the spectrum is always evolving. On large scales $k <
k_{*,{\rm osc}}$ there is both an amplification of the spectrum and a damping evolution. On small scales $k > k_{*,{\rm osc}}$ the power spectrum shows superimposed damped oscillations.

In the instantaneous case $\Delta=0$ (Model 2), the characteristic comoving scale is $k_{*,{\rm inst}} \equiv k_{*,{\rm osc}}$, the same as Eq. (\[kosc\]). In the model 2, on large scales $k < k_{*,{\rm
    inst}}$ the spectral shape is unchanged during the bounce and only the amplitude increases, as identified in Eq. (\[Amplitude\]). On smaller scales $k > k_{*,{\rm inst}}$ there is a change in the spectral index and the power spectrum, as in the case of model 1, shows superimposed damped oscillations.

The results Fig. \[bounceplot1\](b) show that for Model 1 (black dashed line) the spectral tilt always decreases with the momentum. The discontinuity at $k_{*,{\rm flat}}$ (denoted by the black dashed vertical line) is an unphysical feature that appears as a consequence of the shape we have considered and should not appear in realistic smooth shapes. The same is true for the Model 2 case (red solid line), where the discontinuity happens at the characteristic scale $k_{*,{\rm
    inst}}$ in this case and comes from the kink like shape considered in this model. Other than that, the spectral index is unchanged (and null) for large scales modes $k< k_{*,{\rm inst}}$ and then decreases for small scale modes $k> k_{*,{\rm inst}}$ and agrees with that of model 1 from this point on, where the index of the power spectrum for both models acquires a large red tilt, and there are superimposed oscillations.

![\[bounceplot2\] The power spectrum for different parameters, namely $q_1= q_2 = 2/3$, $\epsilon = 0.01 \tau_B$, for the cases of $\Delta=0.7 \tau_B$ (blue dotted line), and $\Delta=0.9
    \tau_B$ (red solid line) for the case of the flat plateau (Model 1) and for the instantaneous case (Model 2), where $\Delta=0$ (black dashed line). Again, the vertical lines denote the positions of the characteristic scales $k_*$ for each model.](bouncefig3.pdf){width="45.00000%"}

Figure \[bounceplot2\] shows the results for different parameter values, parameters for which the two scales $k_*$ for the models with and without a plateau for $a(t)$ are more widely separated than they are for the earlier parameter values. For both models there are oscillations of the power spectrum for $k$ values between the two critical $k_*$ values. These results in particular show that as $\Delta \to \tau_B$, the scale $k_{*,{\rm flat}}$ can occur deeper in the oscillating regime $k>k_{*,{\rm osc}}$ for the spectrum.

A study of the oscillating regime for small scales $k>k_{*,{\rm
    osc}}$, and which is common for both models considered here, is given in the Appendix \[appA\]. In particular, it is shown that the envelope function of the power spectrum for the small scale modes keeps the spectral tilt $n_s = -4q_1+2$, as also seen in the previous Eq. (\[tiltns\]).

Generalization to $n$ small bounces {#sec5}
===================================

In this section, we would like to analyze the case where there are $n$ small bounces see Fig \[nbounceplot\]. Since the transfer matrix of the flat plateau case is quite involved, we would like to first consider Model 2 (no plateau interval) for illustrative purposes.

![\[nbounceplot\] An illustration of the generalization to n-vibrations (or small bounces). ](bouncefig4.pdf){width="45.00000%"}

Based on our previous calculations, we can easily write down the transfer matrices for the coefficient vector. We define the following useful matrices $$\begin{aligned}
M_1(\tau,q)  = \begin{pmatrix} \tau^{1-q} & \tau^q \\ -(1-q) \tau^{-q}
  & -q \tau^{q-1}
\end{pmatrix},\quad
M_2(\tau,q)  = \begin{pmatrix} \tau^{1-q} & \tau^q \\ (1-q) \tau^{-q}
  & q \tau^{q-1}
\end{pmatrix} .\end{aligned}$$

For large scale modes, we define the combination of matrices $$\begin{aligned}
N =
M^{-1}_2(\tau_B,q_1)M_1(\tau_B,q_1)M^{-1}_2(\epsilon,q_1)M_1(\epsilon,q_1)
,\end{aligned}$$ which becomes $$\begin{aligned}
N = \frac{1}{(1-2q_1)^2} \begin{pmatrix} 4(q_1-1)q_1
  \tau_B^{2q_1-1}\epsilon^{1-2q_1} + 1 &
  2q_1(\epsilon^{2q_1-1}-\tau_B^{2q_1-1})
  \\ 2(q_1-1)(\tau_B^{1-2q_1}-\epsilon^{1-21_1}) &
  4(q_1-1)q_1\epsilon^{2q_1-1}\tau_B^{1-2q_1}+1 
\end{pmatrix}.\end{aligned}$$

Taking the two bump model as an example, we obtain the final coefficient vector to be $$\begin{aligned}
\mathcal C_F = M_2^{-1} (\epsilon ,q_2) M_1(\epsilon ,q_1) N M_2^{-1}
(\tau_B, q_1) M_1(\tau_B, q_1) M_2^{-1} (\epsilon,q_1) M_1(\epsilon,
q_2) \mathcal C_I  .\end{aligned}$$ We set the initial coefficient matrix $\mathcal C_I$ to be $$\begin{aligned}
\mathcal C_I = \begin{pmatrix} c_{11} \\ 0
\end{pmatrix},\end{aligned}$$ and then we get $$\begin{aligned}
\nonumber
c_{52} =&\frac{1}{(1-2q_1)(1-2q_2^3)} \left[ 8  c_{11}
  \left(q_1-1\right){}^2 q_1   \left(q_1-q_2+1\right){}^2   \epsilon
  ^{4 q_1-2 q_2-1}   \tau _B^{2-4 q_1} \right.  \\ \nonumber &
  \left. -4 c_{11}  \left(q_1-1\right)   \left(q_1-q_2+1\right)
  \left(2 q_1^2+2 q_2 q_1-5   q_1+q_2-1\right) \epsilon   ^{2 q_1-2
    q_2} \tau _B^{1-2   q_1} \right.  \\ \nonumber & \left. -4 c_{11}
  q_1   \left(q_1+q_2-2\right)   \left(2 q_1^2-2 q_2   q_1+q_1+3
  q_2-4\right)   \epsilon ^{-2 q_1-2 q_2+2}   \tau _B^{2 q_1-1}
  \right.  \\ \nonumber & \left. -8 c_{11}   \left(q_1-1\right) q_1^2
  \left(q_1+q_2-2\right){}^2   \epsilon ^{-4 q_1-2 q_2+3}    \tau
  _B^{4 q_1-2} \right.  \\ \nonumber & \left. +2 c_{11}   \left(4
  q_1^4-8 q_1^3-4   q_2^2 q_1^2+16 q_2 q_1^2-10  q_1^2+4 q_2^2 q_1-16
  q_2   q_1+14 q_1+2 q_2^2-5   q_2+3\right) \epsilon ^{1-2   q_2}
  \right] .\end{aligned}$$ On very large scale there is no change in the spectral slope, as expected.

Now we want to deal with the small scale case. We need to define two more matrices $$\begin{aligned}
L_1 = \begin{pmatrix} e^{-ik\tau_H} & e^{i k \tau_H} \\ ik e^{-i k
    \tau_H}  & -i k e^{i k \tau_H}
\end{pmatrix}, \quad
L_2 = \begin{pmatrix} e^{ik\tau_H} & e^{-i k \tau_H} \\ ik e^{i k
    \tau_H}  & -i k e^{-i k \tau_H}
\end{pmatrix}.\end{aligned}$$ Then we have $$\begin{aligned}
\nonumber
\mathcal C_F = & M_2^{-1} (\epsilon,q_2) M_1(\epsilon,q_1) M_1^{-1}
(\delta,q_1) L_2 L_1^{-1} M_2(\delta,q_1) M_2^{-1} (\epsilon,q_1)
\\ &\times M_1(\epsilon,q_1) M_1^{-1} (\delta,q_1) L_2 L_1^{-1}
M_2(\delta,q_1) M_2^{-1} (\epsilon,q_1) M_1(\epsilon,q_2) \mathcal
C_I,\end{aligned}$$ and the general result has the form $$\begin{aligned}
P_{\zeta} = \left[\# + \# (\epsilon k)^{2-4q_1} + \# (\epsilon
  k)^{4q_1-2} + \# (\epsilon k)^{1-2q_1}  + \# (\epsilon
  k)^{2q_1-1}\right]^2 P_i \, .\end{aligned}$$ Since we are interested in modes which exit the Hubble radius before the time $-(\tau_B + \epsilon)$, we consider values of $k$ for which $\epsilon k \ll 1$. Hence, in this range of $k$ values it is the second term above which dominates and we find the scaling $$P_{\zeta} \, \sim \, (\epsilon k)^{4 - 8q_1} P_i  ~.$$ Thus, small scale modes acquire a red tilt compared to the initial spectrum. If the initial spectrum is scale invariant, then the resulting spectral index for small scale modes is $$n_s - 1  =  4 - 8q_1  ~.$$ Similarly, we can obtain the spectral index change for $n$ small bounces, which is $$n_s - 1  =  (2 - 4q_1)n  ~.$$

Summary and Conclusions {#sec6}
=======================

In this paper we have analyzed in detail the power spectrum of curvature fluctuations in a bouncing cosmology in which the bounce phase has small vibrations, i.e., small bounces. To be specific we have mostly considered the case of one small bounce with characteristic time scales $\tau_B$ and $\Delta < \tau_B$ which are much smaller than cosmological times. We have given a detailed study of the necessary matching conditions required to obtain the complete form for the power spectrum. The matchings connect at least five different phases for a given momentum scale which need to be treated with care.

In our study, we have adopted two simplified models for the shape of the vibrations, allowing a complete analytical study. Despite the apparent simplicity of these models, they are already of sufficient complexity to allow to extract similar features that can emerge in more realistic models. In particular, similar structures that we have considered here can appear in bounce models coming from quantum gravity, as those recently proposed in Ref. [@Alesci], which makes this study of particular importance. Our results for the power spectrum shows that there is an amplification of its amplitude and it also tends to get redder at large scales as the number of vibrations increase. At small scales the power spectrum features superimposed damped oscillations.

The reddening of the spectrum for scales which enter the small bounce agrees with the results found in Ref. [@RHBcyclic]. The oscillations in the power spectrum which are seen on small scales are reminiscent of oscillations which are obtained in some other approaches to the [*Trans-Planckian problem*]{} for cosmological fluctuations. For example, if initial conditions are set on a time-like [*new physics hypersurface*]{} [@newphysics] such that modes $k$ are initiated when the physical wavelength associated with $k$ equals a fixed physical length (e.g. the Planck length), and they are initiated in the same state (e.g. the state which locally looks like the Bunch-Davies vacuum [@BD]), then oscillations in the spectrum result.

Both the qualitative and quantitative changes in the power spectrum that we have obtained can produce observed effects in spectrum of cosmological perturbations accessible through the measurements of the cosmic microwave background radiation. These effects can manifest themselves both in pure bouncing cosmologies (no subsequent inflationary period) and in scenarios where there is a post-bounce inflationary phase. For instance, those bounce vibrations can induce particle production, changing the vacuum state such as to be different from the usual Bunch-Davis one, similar to recent pre-inflationary studies in Loop Quantum Cosmology [@Zhu:2017jew]. The modifications we have obtained in this work could then be used to put constraints on these possible features that can appear in these bounce models and which deserve further study. The results we have presented here provides then an important first step in understanding these effects and which we hope to address elsewhere.

Envelope of the Power Spectrum for Small Scale Modes {#appA}
====================================================

In this section, we would like to calculate the envelope of the power spectrum for small scale modes. Since the model without plateau is a special limit of the model with a non-vanishing flat plateau, we just focus on the latter. We can simply set $\Delta \rightarrow 0$ to get the answer for the model without a plateau.

By collecting the relevant terms in the power spectrum, we can write it in the form $$\begin{aligned}
P_{\zeta } = \left[C_1\sin(2k\Delta) + C_2 \cos(2k\Delta)\right]^2 P_i
\, .\end{aligned}$$ The envelope of the power spectrum is thus $$\begin{aligned}
P_{\zeta {\rm (env)}} = (C_1^2 + C_2^2 )P_i \, ,\end{aligned}$$ where the coefficients $C_1$ and $C_2$ are given by $$\begin{aligned}
\nonumber
C_1&  = \frac{k^{-1}}{(1-2q_2)(1-2q_1)^2} \left\{ 2 \delta^{-1}
\left[k^2\delta^2+(q_1-1)q_1 \right](1+q_1-q_2)(-2+q_1+q_2)
\right. \\  & \left. - \delta^{2q_1-2} (q_1-k\delta) (q_1+k\delta)
(-2+q_1+q_2)^2 \epsilon^{1-2q_1} + \delta^{-2q_1}
(q_1-k\delta-1)(q_1+k\delta-1)(1+q_1-q_2)^2\epsilon^{2q_1-1} \right\}, 
\label{C1env}
\\\nonumber C_2 & =  \frac{1}{(1-2q_2)(1-2q_1)^2} \left[
  -2(1+q_1-q_2)(-2+q_1+q_2)  \right. \\ & \left. -\delta^{2q_1-1} 2
  q_1 (-2+q_1+q_2)^2 \epsilon^{1-2q_1} - \delta^{-2q_1+1} 2 (q_1-1)
  (1+q_1-q_2)^2 \epsilon^{2q_1-1} \right] \, .
\label{C2env}\end{aligned}$$

We are interested in the parameter region $\epsilon / \delta \ll1$ (recall that the time scale $\epsilon$ is expected to be of the order of the Planck scale, whereas $\delta$ will be parametrically larger since it is associated with the time scale of the bounce). We are also interested in the range of values $1/3<p<1$, or equivalently, $1/2<q<+\infty$. We can then determine which are the dominant terms in $C_1$ and $C_2$, which from Eqs. (\[C1env\]) and (\[C1env\]), they are given by $$\begin{aligned}
C_1&  \simeq \frac{k^{-1}}{(1-2q_2)(1-2q_1)^2} \left[  -
  \delta^{2q_1-2} (q_1-k\delta) (q_1+k\delta)  (-2+q_1+q_2)^2
  \epsilon^{1-2q_1} \right], \\ C_2 & \simeq
\frac{1}{(1-2q_2)(1-2q_1)^2} \left[ -\delta^{2q_1-1} 2 q_1
  (-2+q_1+q_2)^2 \epsilon^{1-2q_1} \right] \, .\end{aligned}$$ When $k$ is close to the $k_*$, then $\delta \rightarrow
\tau_B-\Delta$ which is a constant. In this case, taking the square of $C_1$, we get terms with different with spectral indices, but the dominant contribution is the term with the lowest power of $k$, which the gives that the slope of the envelope (for an initial spectrum which is scale-invariant) will be $$\begin{aligned}
n_s -1 = -2 ,\end{aligned}$$ because we have $k\delta<q_1$ in this range. This can be seen from the numerical results shown in Fig. \[bounceplot3\] for the two models we have considered. The change in the spectral slope is due to the matching conditions. Each time, we can get factors of $1/k$ or $k$ when we match the solution across the boundaries of Regions II and III, and of Regions III and IV.

Note that in a generic case when we have a smooth evolution of the scale factor, we expect that there will be no discontinuities in the power spectrum. Thus, in a generic case, we do not expect that we always get an interval of wavenumber with a spectrum of slope $n_s=-2$. What we expect in the case of a smoothly evolving scale factor is that on very large scales, we get a scale invariant spectrum (the actual spectrum, not just the envelope), and then it will smoothly transit to a spectrum with tilt $n_s=-4q_1+2$ when we look at the envelope only. We see oscillations with amplitude given by the envelope function on intermediate and small scales.

The coefficient $C_2^2$ gives a scale invariant power spectrum $$\begin{aligned}
n_s -1 = 0,\end{aligned}$$ but its amplitude is suppressed by $k\delta$ compared to the amplitude of $C_1$. To be a bit more precise (still in the case of constant $\delta$), we can write $$\begin{aligned}
C_1  = A_1 k^{-1} + A_2 k  \, ,\end{aligned}$$ where the constants $A_1$ and $A_2$ are $$\begin{aligned}
\nonumber
A_1 & = \frac{1}{(1-2q_2)(1-2q_1)^2} \left[ 2 \delta^{-1} (q_1-1)q_1
  (1+q_1-q_2) (-2+q_1+q_2)  - \delta^{2q_1-2} q_1^2 (-2+q_1+q_2)^2
  \epsilon^{1-2q_1} \right. \\ &\left. -\delta^{-2q_1} (q_1-1)^2
  (1+q_1-q_2)^2 \epsilon^{2q_1-1}  \right],   \\ \nonumber A_2 & =
\frac{1}{(1-2q_2)(1-2q_1)^2} \left[ 2\delta(1+q_1-q_2)(-2+q_1+q_2)  +
  \delta^{2q_1} (-2+q_1+q_2)^2 \epsilon^{1-2q_1} \right.  \\ &- \left.
  \delta^{-2q_1+2} (1+q_1-q_2)^2 \epsilon^{2q_1-1}  \right] .\end{aligned}$$ The spectral index is computed as $$\begin{aligned}
n_s - 1 = \frac{d\ln P_{\zeta {\rm (env)}}}{d\ln k} = \frac{2
  P_i}{P_{\zeta {\rm (env)}}} (A_2^2 k^2 -A_1^2 k^{-2}).\end{aligned}$$ The power spectrum is hence comprised of several terms with different spectral tilts $n_s$ $$\begin{aligned}
n_s-1 = 2, 1, 0, -1, -2,\end{aligned}$$ More generally (for larger values of $k$ when $\delta$ is not constant), we have

$$\begin{aligned}
P_{\zeta } = \left\{D_1\sin[2(k\tau_B-q_1)] + D_2
\cos[2k(k\tau_B-q_1)]\right\}^2 P_i \, .\end{aligned}$$

The envelope of the power spectrum is thus $$\begin{aligned}
P_{\zeta {\rm (env)}} = (D_1^2 + D_2^2 )P_i \, ,\end{aligned}$$ where the coefficients $D_1$ and $D_2$ are given by $$\begin{aligned}
D_1 & = \frac{1}{(1-2q_2)(1-2q_1)} \left[ -2 (1+q_1-q_2) (-2+q_1+q_2)
  +  k^{2q_1-1} q_1^{-2q_1} (1+q_1-q_2)^2 \epsilon^{2q_1-1} \right] ,
\\ D_2 & = \frac{1}{(1-2q_2)(1-2q_1)^2} \left[
  -2(1+q_1-q_2)(-2+q_1+q_2) +  k^{1-2q_1} q_1^{2q_1-1}
  (-2q_1)(-2+q_1+q_2)^2\epsilon^{1-2q_1}  \right. \nonumber \\ &
  \left. -2 k^{2q_1-1} q_1^{-2q_1+1}  (q_1-1) (1+q_1-q_2)^2
  \epsilon^{2q_1-1} \right] \, .\end{aligned}$$

We now can see that this envelope function reproduces the result of Ref. [@RHBcyclic]. We have $$\begin{aligned}
D_1  &= B_1 + B_2 k^{2q_1-1},\\ D_2  &= E_1 + E_2 k^{1-2q_1} + E_3
k^{2q_1-1} \, ,\end{aligned}$$ where the constants $B_1$, $B_2$, $E_1$, $E_2$ and $E_3$ are given by $$\begin{aligned}
B_1  &= \frac{-2 (1+q_1-q_2) (-2+q_1+q_2)}{(1-2q_2)(1-2q_1)},\\ B_2
&= \frac{ q_1^{-2q_1} (1+q_1-q_2)^2 \epsilon^{2q_1-1}
}{(1-2q_2)(1-2q_1)},\\  E_1  &=
\frac{-2(1+q_1-q_2)(-2+q_1+q_2)}{(1-2q_2)(1-2q_1)^2},\\ E_2  &=
\frac{q_1^{2q_1-1}
  (-2q_1)(-2+q_1+q_2)^2\epsilon^{1-2q_1}}{(1-2q_2)(1-2q_1)^2}, \\ E_3
&= \frac{-2q_1^{-2q_1+1}  (q_1-1) (1+q_1-q_2)^2
  \epsilon^{2q_1-1}}{(1-2q_2)(1-2q_1)^2} .\end{aligned}$$

The spectral tilt is then given by $$\begin{aligned}
n_s - 1  & = \frac{P_i}{P_{\zeta {\rm (env)}}}  \left\{ 2 D_1 B_2
(2q_1-1) k^{2q_1-1} +  2 D_2 \left[E_2(1-2q_1) k^{-2q_1+1} + E_3
  (2q_1-1) k^{2q_1-1} \right]  \right\} \, .
\label{nsenv}\end{aligned}$$ The expression (\[nsenv\]) is comprised of several terms with spectral tilts $n_s$ given by $$\begin{aligned}
n_s-1 = 4q_1-2, 2q_1-1, 0, -4q_1+2, -2q_1+1 \, .\end{aligned}$$ Since are interested in modes with $k\epsilon<1$ and parameter values $1/3<p<1$ (or, equivalently, $1/2<q<+\infty$) we can determine the dominant terms in $D_1$ and $D_2$ and find them to be $$\begin{aligned}
D_1  \rightarrow 0,\quad D_2  \rightarrow E_2 k^{1-2q_1} \, .\end{aligned}$$ Thus, the dominant contribution to the power spectrum is $$\begin{aligned}
P_{\zeta {\rm (env)} } = E_2^2 k^{2-4q_1} P_i \, ,\end{aligned}$$ which corresponds to a spectral tilt of $$\begin{aligned}
n_s = -4q_1+2 \, .\end{aligned}$$

Acknowledgments {#acknowledgments .unnumbered}
===============

One of us (R.B.) is grateful to Emanuele Alesci and Stefano Liberati for discussions about the model of [@Alesci] which led to this project. He also thanks Stefano Liberati and the other organizers of the [*Probing the Spacetime Fabric: from Concepts to Phenomenology*]{} workshop help in July 2017 at SISSA for inviting him to participate and speak. The research at McGill was supported in part by an NSERC Discovery grant and by the Canada Research Chair program. Q.L acknowledge financial support from the University of Science and Technology of China, and from the CAST Young Elite Scientists Sponsorship Program (2016QNRC001), and by the NSFC (grant Nos. 11421303, 11653002). SZ is supported by the Hong Kong PhD Fellowship Scheme (HKPFS) issued by the Research Grants Council (RGC) of Hong Kong. R.O.R is partially supported by research grants from Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), grant No. 303377/2013-5 and Fundação Carlos Chagas Filho de Amparo á Pesquisa do Estado do Rio de Janeiro (FAPERJ), grant No. E - 26/201.424/2014. Q.L., R.O.R. and S.Z. are grateful for the hospitality of the Physics Department at McGill University during research visits when this work was initiated.

[999]{}

R. Brandenberger and P. Peter, “Bouncing Cosmologies: Progress and Problems,” Found. Phys.  [**47**]{}, no. 6, 797 (2017). F. Finelli and R. Brandenberger, “On the generation of a scale invariant spectrum of adiabatic fluctuations in cosmological models with a contracting phase,” Phys. Rev. D [**65**]{}, 103522 (2002). D. Wands, “Duality invariance of cosmological perturbation spectra,” Phys. Rev.  D [**60**]{}, 023507 (1999). J. Khoury, B. A. Ovrut, P. J. Steinhardt and N. Turok, “The Ekpyrotic universe: Colliding branes and the origin of the hot big bang,” Phys. Rev. D [**64**]{}, 123522 (2001). A. Notari and A. Riotto, “Isocurvature perturbations in the ekpyrotic universe,” Nucl. Phys. B [**644**]{}, 371 (2002);\
F. Finelli, “Assisted contraction,” Phys. Lett. B [**545**]{}, 1 (2002);\
F. Di Marco, F. Finelli and R. Brandenberger, “Adiabatic and isocurvature perturbations for multifield generalized Einstein models,” Phys. Rev. D [**67**]{}, 063512 (2003);\
J. L. Lehners, P. McFadden, N. Turok and P. J. Steinhardt, “Generating ekpyrotic curvature perturbations before the big bang,” Phys. Rev. D [**76**]{}, 103501 (2007);\
E. I. Buchbinder, J. Khoury and B. A. Ovrut, “New Ekpyrotic cosmology,” Phys. Rev. D [**76**]{}, 123503 (2007);\
P. Creminelli and L. Senatore, “A Smooth bouncing cosmology with scale invariant spectrum,” JCAP [**0711**]{}, 010 (2007). Y. F. Cai, T. Qiu, Y. S. Piao, M. Li and X. Zhang, “Bouncing universe with quintom matter,” JHEP [**0710**]{}, 071 (2007);\
C. Lin, R. H. Brandenberger and L. Perreault Levasseur, “A Matter Bounce By Means of Ghost Condensation,” JCAP [**1104**]{}, 019 (2011);\
T. Qiu, J. Evslin, Y. F. Cai, M. Li and X. Zhang, “Bouncing Galileon Cosmologies,” JCAP [**1110**]{}, 036 (2011);\
D. A. Easson, I. Sawicki and A. Vikman, “G-Bounce,” JCAP [**1111**]{}, 021 (2011);\
Y. F. Cai, D. A. Easson and R. Brandenberger, “Towards a Nonsingular Bouncing Cosmology,” JCAP [**1208**]{}, 020 (2012);\
Y. F. Cai, E. McDonough, F. Duplessis and R. H. Brandenberger, “Two Field Matter Bounce Cosmology,” JCAP [**1310**]{}, 024 (2013). Y. F. Cai, T. Qiu, R. Brandenberger, Y. S. Piao and X. Zhang, “On Perturbations of Quintom Bounce,” JCAP [**0803**]{}, 013 (2008);\
Y. F. Cai and X. Zhang, “Evolution of Metric Perturbations in Quintom Bounce model,” JCAP [**0906**]{}, 003 (2009);\
Y. F. Cai, T. t. Qiu, R. Brandenberger and X. m. Zhang, “A Nonsingular Cosmology with a Scale-Invariant Spectrum of Cosmological Perturbations from Lee-Wick Theory,” Phys. Rev. D [**80**]{}, 023511 (2009). R. Brandenberger, “Matter Bounce in Horava-Lifshitz Cosmology,” Phys. Rev. D [**80**]{}, 043516 (2009). T. Biswas, A. Mazumdar and W. Siegel, “Bouncing universes in string-inspired gravity,” JCAP [**0603**]{}, 009 (2006);\
T. Biswas, R. Brandenberger, A. Mazumdar and W. Siegel, “Non-perturbative Gravity, Hagedorn Bounce & CMB,” JCAP [**0712**]{}, 011 (2007);\
A. S. Koshelev, “Stable analytic bounce in non-local Einstein-Gauss-Bonnet cosmology,” Class. Quant. Grav.  [**30**]{}, 155001 (2013). M. Bojowald, “Loop quantum cosmology,” Living Rev. Rel.  [**11**]{}, 4 (2008);\
A. Ashtekar, “Singularity Resolution in Loop Quantum Cosmology: A Brief Overview,” J. Phys. Conf. Ser.  [**189**]{}, 012003 (2009);\
A. Ashtekar and A. Barrau, “Loop quantum cosmology: From pre-inflationary dynamics to observations,” Class. Quant. Grav.  [**32**]{}, no. 23, 234001 (2015). N. Turok, B. Craps and T. Hertog, “From big crunch to big bang with AdS/CFT,” arXiv:0711.1824 \[hep-th\];\
B. Craps, T. Hertog and N. Turok, “On the Quantum Resolution of Cosmological Singularities using AdS/CFT,” Phys. Rev. D [**86**]{}, 043513 (2012);\
C. S. Chu and P. M. Ho, “Time-dependent AdS/CFT duality and null singularity,” JHEP [**0604**]{}, 013 (2006);\
C. S. Chu and P. M. Ho, “Time-dependent AdS/CFT duality. II. Holographic reconstruction of bulk metric and possible resolution of singularity,” JHEP [**0802**]{}, 058 (2008);\
R. H. Brandenberger, E. G. M. Ferreira, I. A. Morrison, Y. F. Cai, S. R. Das and Y. Wang, “Fluctuations in a cosmology with a spacelike singularity and their gauge theory dual description,” Phys. Rev. D [**94**]{}, no. 8, 083508 (2016);\
E. G. M. Ferreira and R. Brandenberger, “Holographic Curvature Perturbations in a Cosmology with a Space-Like Singularity,” JCAP [**1607**]{}, no. 07, 030 (2016). C. Kounnas, H. Partouche and N. Toumbas, “S-brane to thermal non-singular string cosmology,” Class. Quant. Grav.  [**29**]{}, 095014 (2012). R. H. Brandenberger and C. Vafa, “Superstrings in the Early Universe,” Nucl. Phys. B [**316**]{}, 391 (1989). J. Martin and R. H. Brandenberger, “The TransPlanckian problem of inflationary cosmology,” Phys. Rev. D [**63**]{}, 123501 (2001);\
R. H. Brandenberger and J. Martin, “The Robustness of inflation to changes in superPlanck scale physics,” Mod. Phys. Lett. A [**16**]{}, 999 (2001). Y. S. Piao, “Primordial Perturbation in Horava-Lifshitz Cosmology,” Phys. Lett. B [**681**]{}, 1 (2009);\
X. Gao, Y. Wang, R. Brandenberger and A. Riotto, “Cosmological Perturbations in Horava-Lifshitz Gravity,” Phys. Rev. D [**81**]{}, 083508 (2010);\
X. Gao, Y. Wang, W. Xue and R. Brandenberger, “Fluctuations in a Horava-Lifshitz Bouncing Cosmology,” JCAP [**1002**]{}, 020 (2010);\
A. Cerioni and R. H. Brandenberger, “Cosmological Perturbations in the Projectable Version of Horava-Lifshitz Gravity,” JCAP [**1108**]{}, 015 (2011). R. H. Brandenberger, C. Kounnas, H. Partouche, S. P. Patil and N. Toumbas, “Cosmological Perturbations Across an S-brane,” JCAP [**1403**]{}, 015 (2014). B. Xue and P. J. Steinhardt, “Unstable growth of curvature perturbation in non-singular bouncing cosmologies,” Phys. Rev. Lett.  [**105**]{}, 261301 (2010);\
B. Xue and P. J. Steinhardt, “Evolution of curvature and anisotropy near a nonsingular bounce,” Phys. Rev. D [**84**]{}, 083520 (2011). J. Quintin, Z. Sherkatghanad, Y. F. Cai and R. H. Brandenberger, “Evolution of cosmological perturbations and the production of non-Gaussianities through a nonsingular bounce: Indications for a no-go theorem in single field matter bounce cosmologies,” Phys. Rev. D [**92**]{}, no. 6, 063532 (2015);\
Y. B. Li, J. Quintin, D. G. Wang and Y. F. Cai, “Matter bounce cosmology with a generalized single field: non-Gaussianity and an extended no-go theorem,” JCAP [**1703**]{}, no. 03, 031 (2017). E. Alesci, G. Botta, F. Cianfrani and S. Liberati, “Cosmological singularity resolution from quantum gravity: the emergent-bouncing universe,” Phys. Rev. D [**96**]{}, no. 4, 046008 (2017). R. H. Brandenberger, “Processing of Cosmological Perturbations in a Cyclic Cosmology,” Phys. Rev. D [**80**]{}, 023535 (2009). R. Durrer and F. Vernizzi, “Adiabatic perturbations in pre - big bang models: Matching conditions and scale invariance,” Phys. Rev. D [**66**]{}, 083503 (2002). V. F. Mukhanov, H. A. Feldman and R. H. Brandenberger, “Theory of cosmological perturbations. Part 1. Classical perturbations. Part 2. Quantum theory of perturbations. Part 3. Extensions,” Phys. Rept.  [**215**]{}, 203 (1992). J. c. Hwang and E. T. Vishniac, “Gauge-invariant joining conditions for cosmological perturbations,” Astrophys. J.  [**382**]{}, 363 (1991). N. Deruelle and V. F. Mukhanov, “On matching conditions for cosmological perturbations,” Phys. Rev. D [**52**]{}, 5549 (1995). W. Israel, “Singular hypersurfaces and thin shells in general relativity,” Nuovo Cim. B [**44S10**]{}, 1 (1966) \[Nuovo Cim. B [**48**]{}, 463 (1967)\] \[Nuovo Cim. B [**44**]{}, 1 (1966)\]. R. H. Brandenberger, “The Matter Bounce Alternative to Inflationary Cosmology,” arXiv:1206.4196 \[astro-ph.CO\]. R. Easther, B. R. Greene, W. H. Kinney and G. Shiu, “A Generic estimate of transPlanckian modifications to the primordial power spectrum in inflation,” Phys. Rev. D [**66**]{}, 023518 (2002);\
V. Bozza, M. Giovannini and G. Veneziano, “Cosmological perturbations from a new physics hypersurface,” JCAP [**0305**]{}, 001 (2003). T. S. Bunch and P. C. W. Davies, “Quantum Field Theory in de Sitter Space: Renormalization by Point Splitting,” Proc. Roy. Soc. Lond. A [**360**]{}, 117 (1978). T. Zhu, A. Wang, G. Cleaver, K. Kirsten and Q. Sheng, “Pre-inflationary universe in loop quantum cosmology,” Phys. Rev. D [**96**]{}, 083520 (2017).

[^1]: Note that there are models in which the spectrum of scalar fluctuations is boosted by a factor independent of wavelength on large scales during the bounce phase [@Jerome].

[^2]: From the point of view of string theory, one may view the time interval we are cutting out as the string time scale, the time scale where the effective field theory description will break down.

[^3]: Note that applying these matching conditions directly between a contracting phase and an expanding phase may be a bit suspect since the background does not satisfy the matching conditions (see Ref. [@Durrer] for a detailed discussion of this point). However, as long as the matching surface is unambiguously determined, the matching conditions for the fluctuations can indeed be applied.

[^4]: Note that the “regions” defined here are not the same as the “phases” defined above. The “phases” refer to particular behaviors of the scale factor, the “regions” to particular behaviors of the fluctuation modes. Phases I and V are equal to Regions I and V, but for the others there is a difference.

[^5]: Where we are neglecting the coefficient $c_{12}$ of the decaying mode in the initial phase.
