1 Introduction↩︎

Non-singular bounce cosmologies [1]–[3] serve as a competitive early universe scenario alternative to inflation. By introducing a contracting phase prior to the standard Big Bang expansion, non-singular bounce cosmologies resolve the initial singularity [4], [5] and the trans-Planckian problem [6], [7] that puzzles inflation. Furthermore, bounce cosmologies can equally explain the formation of large-scale structure and the observed cosmic microwave backgrounds [8]–[11], thus consistent with astrophysical observations [12], [13].

Given the appealing feature of non-singular bounce cosmologies, it is natural to ask if it can be experimentally distinguished from inflation. In the phase of contraction, the energy density of anisotropic stress grows rapidly since it scales as \(a^{-6}\) where \(a\) is the scale factor of the universe. To evade the over-production of anisotropies, there must be a contracting phase with effective Equation-of-state (EoS) parameter \(w_c > 1/3\) so that the background energy density \(\rho_{\rm bg} \propto a^{-3(1+w_c)}\) grows faster than that of anisotropies [14]–[18]. Primordial tensor fluctuations that cross the Hubble horizon in this phase has a strongly blue-tilted spectrum, i.e., the spectra index satisfying \(2 < n_T < 3\). The existence of the blue spectrum at certain frequency regime is regarded as a crucial feature of bounce cosmologies, and have wide phenomenological implications, e.g., the interpretations of power deficit in the CMB TT spectrum [19]–[21], the recent ACT observation [22], and pulsar timing array (PTA) signals [23]–[27]. Nonetheless, inflation can as well yield a blue tensor spectrum [28]–[30] at certain scales, see [31] for a review. Hence, a blue-tilted primordial GWs spectrum can hardly be a deterministic feature of bounce scenarios. We are thus in a hot pursuit to search for distinctive features of bounce cosmologies.

In this Letter, we report a unique pattern of primordial GWs in non-singular bounce cosmologies that can serve as a distinctive signature to experimentally test it in the GW astronomy. The production of primordial GWs can be understood as the scattering process of relic gravitons by the effective potential \(M_p^2V(\tau) \equiv M_p^2 a^{\prime \prime}/a\). The effective potential has one peak in generic inflation models, and two peaks in non-singular bounce cosmologies due to the existence of contracting phase. This double-peak structure in non-singular bounce cosmologies results in the inteference of GWs between the two peaks in the effective potential, which is absent in inflation. Therefore, the GWs spectrum in non-singular bounce cosmologies has an distinctive oscillatory pattern on high frequency regime that is not seen in generic inflation models. The amplitude of GWs spectrum on such scales is high enough to be detected by current and forthcoming GWs instruments. Our findings then provide a unique avenue to test bounce cosmologies and the underline new physics in GWs astronomy.

2 Primordial Gravitational Waves↩︎

In a spatially flat Friedmann Lemaître Robertson Walker (FLRW) universe, primordial tensor fluctuations \(h_k\) evolve as \[\begin{align} \label{eq:vkdynamical} v_k^{\prime \prime} + \left( k^2 - V(\tau) \right) v_k = 0 ~, \end{align}\tag{1}\] with \(v_k \equiv a h_k/2\) the canonically normalized tensor fluctuations, \(V(\tau) \equiv a^{\prime \prime}/a\), and a prime denotes differentiation with respect to conformal time \(\tau\). In the far past, primordial fluctuations are generated on sub-horizon scales \(k^2\gg V(\tau)\). The vacuum initial condition \(v_k(\tau) = \frac{e^{-i k \tau}}{\sqrt{2k}}\) represents a plane wave propogate forward in time. In the present era, the GWs we can measure are only those with wavelengths shorter than the cosmological horizon \(k^2\gg V(\tau)\). For those modes \(v_k\) also evolves as a plane wave \(v_k = \alpha_k \frac{e^{-i k \tau}}{\sqrt{2k}}+\beta_k \frac{e^{i k \tau}}{\sqrt{2k}}\). Thus, the evolution of primordial GWs can be understood as the scattering of a plane wave by a localized effective potential \(M_p^2V(\tau)\), see Fig. 1 for illustration.

Eq. 1 is reminiscent to the famous one-dimentional potential well problem in quantum mechanics, where \(V(\tau)\) plays the role of potential barrier. It is well-known that in quantum mechanics, resonant tunneling effect takes place when the potential barrier contains multiple peaks, which is absent in the case of single-peak potential barrier. The inteference of waves between the multiple peaks leads to oscillations in the reflected and transmitted waves.

Likewise, the oscillatory pattern emergies in the spectrum of GWs when the effective potential develops multiple-peak structure. The double-peak structure of effective potential \(V(\tau)\) acts as a temporal scattering cavity, and the primordial tensor mode is partially reflected by each peak in \(V(\tau)\). Tensor modes are therefore partially reflected by both peaks, and the two reflected components acquire a relative \(2k|(\tau_2-\tau_1)|\), where \(\tau_1\) and \(\tau_2\) denote the conformal-time positions of the two peaks, where \(\tau_2\) and \(\tau_1\) are the temporal location of the peaks in \(V(\tau)\). Their superposition produces an interference term in \(|\beta_k|^2\), leading to oscillatory features in \(\Omega_{\rm GW}(f)\). Accordingly, the oscillation frequency in \(\Omega_{\rm GW}(f)\) is \(\Delta f \simeq [2a_0|\tau_2-\tau_1|]^{-1}\) where \(a_0\) is the scale factor at today.

For illustrative purposes, we parametrize the double-peak structure as \[\label{eq:VLorentz} V(\tau) = A_2^2 L(\tau; \tau_2,\Delta \tau) - A_1^2 L(\tau;\tau_1,\Delta \tau) ~,\tag{2}\] with \(L(\tau; \tau',\Delta\tau)=\frac{1}{\pi}\frac{\Delta\tau}{(\tau-\tau')^2+\Delta\tau^2}\) the Lorentzian peak. On high frequency regime the Born approximation tells \[\begin{align} \label{eq:betakUV} \beta_k \simeq \frac{i e^{-2k \Delta \tau}}{2k} \left( A_2^2 e^{2ik\tau_2} - A_1^2 e^{2ik\tau_1} \right) ~, \end{align}\tag{3}\] so \(|\beta_k|^2\) develops an oscillatory pattern in momentum space with period \(\Delta k = \pi/|\tau_2 - \tau_1|\), or equivalently \(\Delta f = [2a_0|\tau_2-\tau_1|]^{-1}\) in the observed GWs spectrum, as expected. We provide the technical details of 3 and related discussion in Appendix 7.

In general, the effective potential contains one/two peak(s) in inflationary/bounce cosmologies, respectively (left panel of Fig. 1). As a result, on high-frequency regime, \(|\beta_k|^2\) in bounce cosmologies obtain an additional oscillatory features. Notably, the observed energy density spectrum of relic gravitons is associated to \(|\beta_k|^2\) via \(\Omega_{\rm GW} = \frac{k^4 |\beta_k|^2}{3\pi^2M_p^2H^2a^4}\), with \(\rho_c \equiv 3H_0^2 M_p^2\) the critial energy density and \(H_0\) the Hubble parameter today. Therefore, bounce cosmologies predict an oscillatory pattern on high-frequency GWs that is absent in inflation (right panel of Fig. 1). This feature can be distinguishable and serve as a smoking-gun signature of bounce cosmologies.

Now we analyze the structure of effective potential in inflationary and bounce cosmologies. In inflationary epoch, \(V(\tau) \simeq 2/\tau^2\), so the effective potential is initially positive and monotonously increases as the universe inflates from \(\tau \to -\infty\) to \(\tau \to 0\). After that, effective potential monotonously decreases in reheating phase and finally shrinks to 0 in the standard radiation domonated (RD) epoch. Thus, the effective potential generically contains only one peak in inflationary cosmologies.

To analyze the structure of effective potential in bounce cosmologies, we come to the function \(V^{\prime} (\tau) = 2 \mathcal{H} \mathcal{H}^{\prime} + \mathcal{H}^{\prime \prime}\). The number of peaks in the effective potential equals to the number of zero-points in \(V^{\prime} (\tau)\). In non-singular bounce cosmologies, there must be a bouncing phase where the universe transits from contraction (\(\mathcal{H}<0\)) to expansion (\(\mathcal{H}>0\)). In the bouncing phase, there must be two transition points satisfying \(\mathcal{H}^{\prime} = 0\), which we label as \(\tau_I\) and \(\tau_{II}\), as shown in Fig. 2. It’s straightforward to see from Fig. 2 that \(\mathcal{H}^{\prime \prime}(\tau_I) > 0\) and \(\mathcal{H}^{\prime \prime}(\tau_{II}) < 0\), so \[V^{\prime}(\tau_I) > 0 ~,~ V^{\prime}(\tau_{II}) < 0 ~.\] Namely, there is a zero-point of \(V^{\prime}(\tau)\) in the regime \(\tau_I < \tau < \tau_{II}\) due to the mean value theorem. As a result, a peak of effective potential emerges at that place.

The emergence of the second peak is a direct relic of the dynamics governing the contracting phase. According to the previous discussion, there should be a contracting phase with \(w_c > 1/3\) to bypass the anisotropic problem. In this phase \[V^{\prime}(\tau) = \frac{4(1-3w_c)}{(1+3w_c)^2 (-\tau)^3} < 0 ~.\] This contracting phase takes place before the first transition event \(\tau = \tau_I\), where \(V^{\prime}(\tau_I) > 0\). Therefore, there must be a second peak in the regime \(\tau < \tau_I\) due to the mean value theorem.

We conclude that the effective potential always has at least two peaks in non-singular bounce models that is free from anisotropic stress problem. This double-peak structure of effective potential in those models leads to unique oscillatory features in energy density spectrum of GWs on high-frequency regime absent in generical inflation models. Consequently, this oscillatory signature provides a key observational discriminant between early universe scenarios.

3 Features of Gravitational Waves in Bouncing Cosmologies↩︎

The dynamics of a bouncing universe introduce three pivot scales that shape the GWs spectrum. The infrared and intermediate scales, \(k_{\rm IR}\) and \(k_{\rm IM}\) , are those who crosses the Hubble horizons at the end of the reheating and contracting phase. The ultraviolet (UV) scale, \(k_{\rm UV} \equiv \sqrt{\max |V(\tau)|}\), corresponds to the smallest perturbations, which cross the Hubble horizon exactly once. Of specific phenomenological interest is the regime \(k_{\rm IR} < k_{\rm IM}\), wherein the reheating epoch contributes meaningfully to the primordial gravitational wave background. This is the scenario we shall consider throughout the paper. The spectrum of GWs exhibits generical properties on different scales:

• Primordial fluctuations with \(k > k_{\rm IM}\) are sub-horizon during the contracting phase. They can hardly feel the contraction of the universe, so \(|\beta_k|\) is relevant to the peak-structure of effective potential only. The oscillatory feature takes place when \(k \geq k_{\rm UV}\), where the Born approximation holds and 3 tells \[\Omega_{\rm GW} \propto k^4 \left[ 1 + \kappa^2 \sin (\omega k/k_{\rm UV}) \right] e^{-\mu (k/k_{\rm UV} )} ~.\] Notably, the pole structure of effective potential gives an additional \(k^2\) factor when performing contour integrations, so \(\Omega_{\rm GW}\) scales as \(k^4\) instead of \(k^2\) [32]. In the regime \(k_{\rm IM} < k \ll k_{\rm UV}\), the solution of Eq. 1 is the Jost function, whose leading-order amplitude is \(|\beta_k|^2 \simeq 1 - 4k^2 a_s^2 + \mathcal{O}(k^3)\) where \(a_s\) reflects the structure of peaks.

• Primordial fluctuations with \(k < k_{\rm IM}\) cross the Hubble horizon during the contracting phase. Since these fluctuations remain super-horizon at the peak locations, their subsequent evolution is left essentially untouched by the double-peak structure. For fluctuations re-enter the Hubble horizon during the RD epoch, namely \(k < k_{\rm IR}\), the corresponding spectra index of primoridal GWs is \(n_{\rm IR} \equiv 3 - \frac{3|w_c - 1|}{1+3w_c}\) [33]. We see the spectra index satisfies \(2 < n_{\rm IR} < 3\) for \(w_c > 1/3\), which is in agreement with the previous analysis. On the other hand, primordial fluctuations with \(k_{\rm IR} < k < k_{\rm IM}\) re-enter the horizon in reheating phase instead. The corresponding tensor spectrum is modified by a factor \(|\chi_k|^2 \propto k^{2\frac{3w_{\rm rh} - 1}{3w_{\rm rh} + 1}}\) by assuming a constant effective EoS parameter \(w_{\rm rh}\); see Appendix 8 for details. Therefore, the spectra index in the regime \(k_{\rm IR} \lesssim k\lesssim k_{\rm IM}\) is \(n_{\rm IM} \equiv n_{\rm IR} + 2 \frac{3w_{\rm rh} - 1}{3w_{\rm rh} + 1}\).

We numerically verify the asymptotic behaviors discussed above in Appendix 6. At low-frequency regime, the spectrum exhibits a broken-power-law behavior, with spectra index \(n_{\rm IR}\) on \(f \leq f_{\rm IR}\) and \(n_{\rm IM}\) on \(f_{\rm IR} < f \leq f_{\rm IM}\). At high-frequency regimes, the spectrum slowly decreases on \(f_{\rm IM} \leq f \leq f_{\rm UV}\), and then decays as \(k^4e^{-\mu k}\) with a characteristic osclllations on \(f \geq f_{\rm UV}\).

In the end, the amplitude of GWs spectrum can be estimated at the intermediate scale \(f = f_{\rm IM}\) in Appendix 9 \[\begin{align} \label{eq:amplitude} (\Omega_{\rm GW}h^2)(f_{\rm IM}) & \nonumber \sim 3.7 \times 10^{-33} \mathcal{A}^2 \left( \frac{H_c}{H_{\rm RD}} \right)^2 \\ & \times \left( \frac{f_{\rm IM}}{1 {\rm Hz}} \right)^4 \left( \frac{f_{\rm IR}}{f_{\rm IM}} \right)^{\frac{1 + 3w_{\rm rh}}{3(1+w_{\rm rh})}} ~. \end{align}\tag{4}\] where \(H_c\), \(H_{\rm RD}\) the Hubble parameter at the end of contraction/reheating phase, respectively, \(\mathcal{A}\) represents the amplitude factor of GWs from bouncing phase due to the tachyonic amplifications.

Notably, the amplitude 4 can be large enough and detectable by future observations with a physically motivated model parameters for two reasons. First, the blue-tilted GWs spectrum from the contraction phase with \(n_{\rm IR}\) can lead to significant growth of energy density spectrum of GWs on high-frequency regime. For instance, even a vanishingly small tensor-to-scalar ratio at Cosmic Microwave Background (CMB) regime can lead to a sizable primordial GWs signatures testable by Pulsar Timing Array (PTA) experiments at Nano-Hertz regime [34]–[36]. This scale dependence has been translated into \(H_c/H_{\rm RD}\) in 4 . A possible huge hierarchy \(H_c/H_{\rm RD}\) can arise from the energy scale difference between the Null-Energy-Condition violation physics and the reheating phase, or equivalently a relatively long reheating phase, see Appendix 9 for more details. Second, the bouncing phase may further amplify the primordial GWs due to tachyonic instabilities [9], [37].

Armed with the amplitude at \(f = f_{\rm IM}\) and the scale-dependent behavior discussed above, we can translate the model parameters (summarized in Table. 1) into a predicted spectrum for primordial GWs.

4 Forecast for High-frequency Gravitational-wave Experiments↩︎

We present the resulting GWs spectrum and compare it with the experimental sensitivity curves in Fig. 3. It is worth noting that the infrared tail of the spectrum lies beyond measurable reach. This is unsurprising given its origin: the IR pivot frequency \(f_{\rm IR}\) is governed by the Hubble parameter during radiation domination via \(\left( \frac{f_{\rm IR}}{1 {\rm Hz}} \right)^2 = 1.1 \times 10^{-27} \left( \frac{H_{\rm RD}}{1 {\rm GeV}} \right)\) in Appendix 9. The condition \(H_{\rm RD} \ll M_p\) then forces \(f_{\rm IR} \ll 10^{-5} {\rm Hz}\), relegating it to frequencies far lower than those we can probe.

On the other hand, high-frequency primordial GWs \(f >f_{\rm IR}\) are experimentally testable. The blue-tilted GWs spectrum on intermediate scales \(f_{\rm IR} \ll f < f_{\rm IM}\) falls within the sensitivity windows of TianQin [40], Taiji [41], and LISA [42], as illustrated in the left panel of Fig. 3. Meanwhile, the oscillatory features emerging at higher frequencies \(f \gtrsim f_{\rm UV}\) are accessible to next-generation detectors such as BBO [43] and DECIGO [44]. A combined observation of the primordial GWs spectrum spanning the LISA and BBO frequency bands would therefore constitute a distinctive fingerprint of bounce cosmologies, as the specific pattern of signals across these different frequency regimes is both unique and characteristic.

Remarkably, the distinctive signal persists into higher frequency regimes. The right panel of Fig. 3 demonstrates that around \(10^2\) Hz, the GWs spectrum enters the detection windows of upcoming terrestrial facilities, including Cosmic Explorer (CE) [45] and the Einstein Telescope (ET) [46]. A detection in this band, revealing the characteristic blue spectrum followed by oscillatory patterns, would offer independent and corroborating evidence for the bounce scenario.

As suggested by Eq. 4 , the amplitude of the gravitational wave spectrum grows with \(f_{\rm IM}\). Consequently, the signal becomes prominent for large \(f_{\rm IM}\), highlighting the high-frequency gravitational waves (HFGWs) as a promising observational target. While direct detection above MHz frequencies remains technically challenging and largely unexplored [47], [48], recent years have witnessed a growing interest in this field, motivated by its unique significance to early universe cosmology and new physics beyond the Standard Model. Considerable attention has been devoted to not only the HFGWs detection strategies such as the Gertsenshtein effect [49], the electromagnetic phenomenon [50], [51], and the atomic quantum sensors [52], [53], but also experimental initiatives including BAW [54], FLASH [55], [56], ADMX [57]–[59], and ALPHA [60], [61]. Those proposals may soon overcome the Big Bang Nucleosynthesis bound on stochastic backgrounds, offering a compelling opportunity to test our results.

5 Conclusions↩︎

In this Letter, we report a distinctive signature of primordial GWs in nonsingular bounce cosmologies. For bounce models that are free from the anisotropic stress problem, the effective potential inevitably develops a double-peak structure, resulting in the oscillatory patterns of GWs spectrum at high frequency regime due to the interference of GWs between the two peaks. This oscillatory pattern is unique and can serve as a distinctive signature of bounce cosmologies. Our results reveal that the energy spectrum of primordial GWs can be at reach for current and forthcoming terrestrial and space-based GWs experiments. Thus, our finding offers an appealing approach to experimentally test bounce scenarios.

Notably, our main conclusion of Eq. 4 show explicitly that the amplitude of primordial GWs at high frequency after the nonsingular bounce is related to \(H_{\rm RD}\), which is identified as the Hubble parameter at the end of the reheating phase. The detection of GWs signature would therefore enable an independent constraint of \(H_{\rm RD}\), complementing other cosmological probes. In this regard, our findings thereby build up a promisingly novel connection between the GW astronomy and new physics potentially operating during reheating.

Finally, the bouncing phase that connects the contracting and expanding phase may lead to richer structures than the simple two-peak scenario considered in the current Letter. Namely, multiple bouncing phases [62], [63], cyclic cosmologies [64]–[66] and time crystal cosmologies [67], [68] can yields multiple peaks in the effective potential for the propagation of primordial GWs. In such cases, oscillatory signatures still persist, though their detailed forms shall depend on the number and structure of the peaks, analogous to the resonant tunneling effect. A systematic investigation of the GWs spectrum in these scenarios therefore presents a compelling direction for future work.

6 Plotting of Figure 1 in the main text↩︎

In this section, we provide details of the models that we consider in the main text.

The bouncing scenarios we consider is the Ekpyrotic bounce, in which the contraction phase is governed by a stiff matter with effective equation-of-state \(w_c = -1 + 2/(3q) > 1\) [69]. As a result, the conformal Hubble parameter \(\mathcal{H}\) scales as \[\label{eq:Hekp} \mathcal{H} = \frac{q}{(1-q)\tau} < 0 ~;~ \tau < 0 ~,~ 0< q <1/3 ~.\tag{5}\] In the end, the universe enters into the standard radiation-dominated epoch, where the conformal Hubble parameter scales as \[\label{eq:Hrd} \mathcal{H} = \frac{1}{\tau} ~,~ \tau > 0 ~.\tag{6}\] We shall also take into account the possible reheating phase. Conventionally, one defines an effective EoS parameter \(w_{\rm rh}\) during the period of reheating, evaluated by the average of the instantaneous EoS parameter [70]. When drawing the figures, we adpot the minimal setup with \(w_{\rm re} = 0\), a choice wildly applied in inflationary cosmology which can be realized by a scalar field with canonical kinetic term and a mass term around its local minimum. Therefore, the conformal Hubble parameter scales as \(\mathcal{H} \simeq 2/\tau\) in reheating phase. Notice that we select \(\tau = 0\) to be the bouncing point \(\mathcal{H} = 0\), so the reheating phase takes place when \(\tau > 0\) and \(\mathcal{H}_{\rm rh} > 0\) accordingly.

In light of this fact, we adopt the following parametrization \[\begin{align} \label{eq:calHpara} \mathcal{H} \nonumber = \frac{1}{2} \frac{\tau}{\tau^2 + \tau_0^2} & \left\{ \frac{q}{1-q} \left( 1 + \tanh [\omega_c (\tau_c - \tau)] \right) \right. \\ & \nonumber \left. + 2\left( \tanh [\omega_r (\tau_r - \tau)] - \tanh [\omega_c (\tau_c - \tau)] \right) \right. \\ & \left. + \left( \tanh [\omega_r (\tau - \tau_r)] + 1 \right) \right\} ~;~ \tau_r \geq \tau_c ~,~ \tau_c < 0 ~, \end{align}\tag{7}\] where the cosmic evolution reduces to 5 when \(\tau \ll \tau_c\) and to 6 when \(\tau \gg \tau_r\). In addition, a straightforward Taylor expansion around \(\tau = 0\) yields \(\mathcal{H} \simeq 2\tau/\tau_0^2\), so this parametrization 7 also describe a sharp bouncing phase with \(\mathcal{H}^{\prime} \simeq 2\tau_0^{-2}\). When \(\tau_c = \tau_r\), the solution 7 depicts a scenario without reheating phase, and the value \(\tau_r - \tau_c\) can indicate the duration of reheating phase.

We use the following parameter sets to draw Fig. 1 in the main text: \[w_c = 1.2 ~,~ w_{\rm rh} = 0 ~;~ \omega_c = 5 ~,~ \omega_r = 0.5 ~;~ \tau_0 = 1/2 ~,~ \tau_c = -1 ~,~ \tau_r = 8 ~.\] The same parameters are used to numerically evaluate \(\Omega_{\rm GW}\), and we show the result in Fig. 4. The asymptotic behaviors in Fig. 4 are plotted via \[\frac{\Omega_{\rm GW}}{\max |\Omega_{\rm GW}|} = 40 \left( \frac{f}{f_{\rm IM}} \right)^{\frac{66}{23}} ~,~ \text{ Red dashed line in left channel} ~,\] \[\frac{\Omega_{\rm GW}}{\max |\Omega_{\rm GW}|} = \left( \frac{f}{f_{\rm IM}} \right)^{\frac{20}{23}} ~,~ \text{ Blue dashed line in left channel} ~,\] \[\frac{\Omega_{\rm GW}}{\max |\Omega_{\rm GW}|} = 1-0.9^2 \left( \frac{f}{f_{\rm IM}} -1 \right)^{2} ~,~ \text{ Red dashed line in right channel} ~,\] \[\frac{\Omega_{\rm GW}}{\max |\Omega_{\rm GW}|} = 0.1^2 \left( \frac{f}{f_{\rm IM}} \right)^{4} e^{-\frac{5}{3} \frac{f}{f_{\rm IM}}} ~,~ \text{ Blue dashed line in right channel} ~.\] We see that the numerical result indeed satisfies the asymptotic behaviors discussed in the main text.

Finally, the reference scales in Fig. 4 are chosen as \(\mathcal{H}_{\ast} = |\mathcal{H}|_{\rm max}\) and \(\tau_{\ast} = \mathcal{H}_{\ast}^{-1}\).

7 Physical origin of the oscillations in GWs energy density spectrum↩︎

In this section, we explain the physical mechanism producing the oscillations in the GWs spectrum. The propagation of primordial tensor fluctuations is governed by \[\begin{align} \label{eq:app95vkdynamical} v_k^{\prime \prime} + \left( k^2 - V(\tau) \right) v_k = 0 ~, \end{align}\tag{8}\] where \[V(\tau) \equiv a^{\prime \prime}/a ~,\] is the effective potential. On high-frequency regime \[k^2 \gg |V(\tau)| ~,\] the equation 8 can be resolved perturbatively as a series of \(|V|/k^2\). We write 8 as \[v_k''+k^2v_k=V(\tau)v_k ~.\] When the right-hand side is treated perturbatively, the solution is \[v_k(\tau) = \frac{e^{-ik\tau}}{\sqrt{2k}} + \int_{-\infty}^{\tau}d\tau'\, G_k(\tau,\tau')V(\tau')v_k(\tau') ~,\] with the retarded Green’s function of the free operator, \[G_k(\tau,\tau') = \frac{\sin[k(\tau-\tau')]}{k} \Theta(\tau-\tau') ~.\] The leading-order result is acquired by taking \(v_k(\tau')\) inside the integral as the plane wave solution \(e^{-ik\tau'}/\sqrt{2k}\), known as the Born approximation. In the late-time limit, we obtain \[\lim_{\tau \to +\infty}v_k \simeq \frac{1}{\sqrt{2k}} \left[ \left(1 + \alpha_k\right) e^{-ik\tau} + \beta_k e^{ik\tau} \right],\] with \[\beta_k \simeq -\frac{i}{2k} \int_{-\infty}^{+\infty}d\tau\, V(\tau)e^{-2ik\tau} ~. \label{eq:Born95beta}\tag{9}\] Namely, in the high-frequency regime, the production of relic gravitons is controlled by the Fourier component of the effective potential at frequency \(2k\).

Notably, the Bogoliubov coefficient \(\beta_k\) corresponds to the particle production of gravitons, so the energy density spectrum of produced gravitons are \[\Omega_{\rm GW} = \frac{k^4 |\beta_k|^2}{3\pi^2M_p^2H^2a^4} ~,\] with \(\rho_c \equiv 3H_0^2 M_p^2\) the critical energy density and \(H_0\) the Hubble parameter today. Therefore, the interference pattern in the Bogoliubov coefficient \(\beta_k\) leads to the oscillations in the GW energy density spectrum, with a frequency inversely proportional to the separation of the two-peaks in the effective potential \(V(\tau)\).

For illustrative purpose, we parametrize the double-peak in the effective potential as two localized Lotentzian fucntions: \[\label{eq:app95VLorentz} V(\tau) = A_2^2 L(\tau; \tau_2,\Delta \tau) - A_1^2 L(\tau;\tau_1,\Delta \tau) ~,~ L(\tau; \tau',\Delta\tau)=\frac{1}{\pi}\frac{\Delta\tau}{(\tau-\tau')^2+\Delta\tau^2} ~,\tag{10}\] where \(A_1\), \(A_2\), \(\tau_1\), \(\tau_2\) label the amplitude and location of the two peaks. The resulting Bogoliubov coefficient is \[|\beta_k|^2 \simeq \frac{e^{-4k\Delta\tau}}{4k^2}\left[ A_2^4+A_1^4-2A_1^2A_2^2 \cos\left(2k(\tau_2-\tau_1)\right) \right] ~.\] Thus, the oscillatory frequency in the GWs spectrum is \(f_{\rm osci} = (2a_0 |\tau_2 - \tau_1|)^{-1}\), inversely proportional to the distance of the two peaks, as expected.

We conclude that the double-peak structure of effective potential in non-singular bouncing cosmology leads to unique oscillatory features in energy density spectrum of GWs with a frequency inversely proportional to the distance of the two peaks on high-frequency regime.

In the end, we comment on the location of the oscillations. The analyze relies on the high-frequency assumption \(k^2 \gg |V(\tau)|\). In the main text, we have defined the characteristic scale \(k_{\rm UV} \equiv \sqrt{\max |V(\tau)|}\), so the oscillatory feature takes place when \(k \geq k_{\rm UV}\). In terms of frequency ,the oscillations in the GWs spectrum start at \(f \simeq f_{\rm UV} = k_{\rm UV}/(2\pi a_0)\).

8 Evolutions of primordial GWs in reheating phase↩︎

In this section we explain the behavior of primordial GWs in reheating phase. The readers are referred to Ref. [71], [72] for more technical details. The modofication of reheating phase is encoded in the tensor transfer function \(\chi_k(\tau)\), defined as \[h_k(\tau) = h_k^{\rm RH} \chi_k(\tau) ~,\] with \(h_k^{\rm RH}\) the tensor fluctuation at the beginning of reheating. This quantity obeys the same differential equation as \(h_k\): \[\chi_k^{\prime \prime} + 2 \frac{a^{\prime}}{a} \chi_k^{\prime} + k^2 \chi_k = 0 ~.\] Assuming the reheating phase is governed by an effective equation-of-state (EoS) parameter \(w_{\rm rh}\), then the Hubble parameter scales as \(H \propto a^{-3(1+ w_{\rm rh})/2}\), and the general solution of transfer function is the Bessel function \[\label{eq:chik} \chi_k \propto a^{-\frac{3}{4}(1-w_{\rm rh})} \left[ C_k J_{\nu_{\rm rh}} \left( \frac{k}{k_{\rm IR}} a^{\gamma} \right) + D_k J_{-\nu_{\rm rh}} \left( \frac{k}{\gamma k_{\rm IR}} a^{\gamma} \right) \right] ~,\tag{11}\] \[\nu_{\rm rh} \equiv \frac{3(1-w_{\rm rh})}{2(1+3w_{\rm rh})} ~,~ \gamma \equiv \frac{1+3w_{\rm rh}}{2} ~.\] This solution holds only if \(-1/3 < w_{\rm rh} < 1\), a quite generic condition as indicated by [70]. The parameter \(k_{\rm IR}\) represents the scale re-enter the horizon at the end of reheating. Notice that, small-scale primordial fluctuations that never cross the horizon shall be correctly regularized to eliminate the vacuum contributions. We introduce another characteristic scale \(k > k_{\rm UV}\) to denote the regime where the effect of regularization becomes important. So, 11 is valid only for \(k < k_{\rm UV}\), and primordial fluctuations on those scales receive negligible contributions from the subsequent RD epoch, as they re-enter the horizon during reheating. As a result, the primordial energy spectrum receives an correction factor \(\left| \chi_k \right|^2\) on \(k < k_{\rm UV}\).

In IR regime \(k \ll k_{\rm IR}\), the modifications shall be negligible as those fluctuations become sub-horizon before reheating. This condition fixes the coefficients in 11 via \[\lim_{k/k_{\rm IR} \to 0} \chi_k = 1 ~.\] resulting in \[\chi_k \simeq \Gamma(\nu_{\rm rh} + 1) \left( \frac{k}{2\gamma k_{\rm IR}} \right)^{-\nu_{\rm rh}} J_{\nu_{\rm rh}} \left( \frac{k}{\gamma k_{\rm IR}} \right) ~,\] at today where we set the \(a_0 \equiv a_{\rm today} = 1\).Notably, in the intermediate regime \(k_{\rm IR} < k < k_{\rm UV}\), the Bessel function asymptotics to \[\label{eq:Jasym} |J_{\alpha}(z)| \simeq \sqrt{\frac{2}{\pi z}} \cos \left( z - (2\alpha + 1) \frac{\pi}{4} \right) ~,\tag{12}\] so the mean value of spectrum develops an approximate power-law behavior \[\label{eq:nTIM} n_{T,IM} = n_{T,{\rm IR}} + 2 \frac{3w_{\rm rh} - 1}{3w_{\rm rh} + 1} ~.\tag{13}\] Notably, the oscillating term shall contribute a constant value after we averaging the Bessel function over time.

9 Amplitude of GW spectrum↩︎

In this section we briefly explain how we estimate the amplitude of \(\Omega_{\rm GW}\). We label the end of contraction phase as \(\tau_{\rm c}\). Notice that \(\tau_{\rm c}\) labels the point where the constant \(w_c\) assumption is about to break, and one shall not confuse it with \(\tau_I\), the first transition point. Also, we are interested in the case \(k_{\rm IM} > k_{\rm IR}\), otherwise, the reheating phase doesn’t introduce interesting observational signatures.

We can analytically track \(P_{T}(k_{\rm IM})\) at \(\tau = \tau_{\rm IM}\) via [33] \[P_{T}(k_{\rm IM},\tau_{c}) = \frac{2^{2|\nu|-1} \Gamma^2(|\nu|)}{(2\nu-1)^2\Gamma^2(3/2)} \left( \frac{H(\tau_{c})}{2\pi M_p^2} \right)^2 (-k_{\rm IM}\tau_{\rm c})^{3-2|\nu|} ~,~ \nu \equiv \frac{3(w_c - 1)}{2(1+3w_c)} ~.\]

Since we asssumed \(k_{\rm IM} > k_{\rm IR}\), \(k_{\rm IM}\) re-enters the horizon at reheating phase at a time \(\tau_{\rm IM}\), and the corresponding tensor fluctuations receive negligible modifications in the bouncing phase, then we can estimate \[P_{T,{\rm today}}(k_{\rm IM}) = P_{T}(k_{\rm IM},\tau_{\rm IM}) \left( \frac{a_0}{a(\tau_{\rm IM})} \right)^{-2} \simeq P_{T}(k_{\rm IM},\tau_c) \left( \frac{a_0}{a(\tau_{\rm IM})} \right)^{-2} ~.\] In addition, the energy spectrum at today is related to \(P_T\) as \[\label{eq:OGWtoday} \Omega_{\rm GW}(k_{\rm IM}) \simeq \frac{1}{12} \left( \frac{k_{\rm IM}}{a_0 H_0} \right)^2 P_{T,{\rm today}}(k_{\rm IM}) ~,\tag{14}\] where \(a_0\), \(H_0\) the scale factor and Hubble parameter today, so \[\Omega_{\rm GW}(k_{\rm IM}) \simeq \frac{1}{12} \frac{(k_{\rm IM}/a_0)^4}{H_0^2M_p^2} \frac{2^{2|\nu|-3} \Gamma^2(|\nu|)}{(2\nu-1)^2\Gamma^2(3/2) \pi^2} \left( \frac{H(\tau_c)}{H(\tau_{\rm IM})} \right)^2 \left( \frac{2}{1+3w_c} \right)^{3-2|\nu|} ~,\] where we’ve used the horizon crossing condition \(-k_{\rm IM}\tau_c = 2/(1+3w_c)\).

We can simplify \(H(\tau_c)/H(\tau_{\rm IM})\) by noticing that \[\frac{H(\tau_c)}{H(\tau_{\rm IM})} = \frac{H(\tau_c)}{H(\tau_{\rm IR})} \frac{H(\tau_{\rm IR})}{H(\tau_{\rm IM})} = \frac{H(\tau_c)}{H(\tau_{\rm IR})} \left( \frac{k_{\rm IR}}{k_{\rm IM}} \right)^{\frac{1 + 3w_{\rm rh}}{3(1+w_{\rm rh})}} ~.\] Notice that, \(|H(\tau_c)|\) also labels the highest energy scale in contraction phase, and \(H(\tau_{\rm IR})\) is the Hubble parameter at the beginning of RD epoch. Let us introduce \(|H(\tau_c)|\) as \(H_c\) and \(H(\tau_{\rm IR})\) as \(H_{\rm RD}\), then the \(|H_c|/H_{\rm RD}\) characterizes the ratio of energy scales between contraction and RD epochs.

We estimate \(2^{3-2|\nu|}/(1-2\nu)^{1-2\nu}\) as 3, since it ranges from 1 to 8.3 when \(|\nu| \leq 1/2\). We also notice that \(\Gamma(|\nu|)\) and \(2/(1+3w_c)\) is usually order unity, unless we assign some extreme value \(|\nu| \ll 1\) or \(w_c \gg 1\). This leads to \[\label{eq:OGWsimple} (\Omega_{\rm GW}h^2)(f_{\rm IM}) \sim 3.7 \times 10^{-33} \left( \frac{H_c}{H_{\rm RD}} \right)^2 \left( \frac{f_{\rm IM}}{1 {\rm Hz}} \right)^4 \left( \frac{f_{\rm IR}}{f_{\rm IM}} \right)^{\frac{1 + 3w_{\rm rh}}{3(1+w_{\rm rh})}} ~.\tag{15}\]

Unfortunately, the estimation may ignore an important contributions from bouncing phase. Especially, primordial fluctuations may experience an exponentially super-horizon growth due to tachyonic instabilities during bouncing phase [9]. In this case, 15 shall acquire an large amplification factor. For example, the amplification factor for scalar is designed to be around \(\mathcal{O}(10^2)\) larger than tensor to suppress the tensor-to-scalar ratio [37]. We introduce a new parameter \(\mathcal{A}\) to characterize the amplification of tensor fluctuations during bouncing phase, and the final result shall be \[(\Omega_{\rm GW}h^2)(f_{\rm IM}) \sim 3.7 \times 10^{-33} \mathcal{A}^2 \left( \frac{H_c}{H_{\rm RD}} \right)^2 \left( \frac{f_{\rm IM}}{1 {\rm Hz}} \right)^4 \left( \frac{f_{\rm IR}}{f_{\rm IM}} \right)^{\frac{1 + 3w_{\rm rh}}{3(1+w_{\rm rh})}} ~.\] Since \(\mathcal{A}\) represents the influence on primordial GWs from the bouncing phase, the measurement of \(\mathcal{A}\) from GW experiments can probe the microscopic physics of bouncing phase and the underlying new physics of NEC violation.

In the end, we comment that if we assume the standard \(\Lambda\)CDM scenario, the relativistic degree of freedom in RD epoch is \(g_{\ast} = 106.75\) and \(H_{\rm RD}\) is related to \(f_{\rm IR}\) via \[\left( \frac{f_{\rm IR}}{1 {\rm Hz}} \right)^2 = 1.1 \times 10^{-27} \left( \frac{H_{\rm RD}}{1 {\rm GeV}} \right) ~.\]

10 A generic parametrization of primordial GWs and plotting details of Fig. 3↩︎

In light of the properties of primordial GWs discussed in the main text, we capture the essense of GWs spectrum as \[\begin{align} \label{eq:OGWtemplete} (\Omega_{\rm GW}h^2) & \nonumber = 3.7 \times 10^{-33} \mathcal{A}^2 \left( \frac{H_c}{H_{\rm RD}} \right)^2 \left( \frac{f_{\rm IM}}{1 {\rm Hz}} \right)^4 \left( \frac{f_{\rm IR}}{f_{\rm IM}} \right)^{\frac{1 + 3w_{\rm rh}}{3(1+w_{\rm rh})}} \\ & \nonumber \times \left\{ \left( \frac{f}{f_{\rm IR}} \right)^{n_{\rm IR}} \left( \frac{f_{\rm IR}}{f_{\rm IM}} \right)^{n_{\rm IM}} [1 - \mathcal{D}(f;f_{\rm IR})] \right. \\ & \nonumber \left. + \left( \frac{f}{f_{\rm IM}} \right)^{n_{\rm IM}} [ \mathcal{D}(f;f_{\rm IR}) - \mathcal{D}(f;f_{\rm IM})] \right. \\ & \nonumber \left. + \frac{1 - 4a_s^2 \frac{f^2}{f_{\rm UV}^2}}{1 - 4a_s^2 \frac{f_{\rm IM}^2}{f_{\rm UV}^2}} [ \mathcal{D}(f;f_{\rm IM}) - \mathcal{D}(f;f_{\rm UV})] \right. \\ & + \left. \frac{1 - 4a_s^2}{1 - 4a_s^2 \frac{f_{\rm IM}^2}{f_{\rm UV}^2}} e^{-\mu \frac{f-f_{\rm UV}}{f_{\rm UV}} } \mathcal{D}(f;f_{\rm UV}) \left[ 1 + \kappa^2 \sin \left(\omega \frac{f - f_{\rm UV}}{f_{\rm UV}} \right) \right] \right\} ~, \end{align}\tag{16}\] where \[\mathcal{D}(f;f_0) \equiv \frac{1 + \tanh [(f-f_0)/f_{\#}]}{2} ~,\] are introduced to smoothly connect each pieces.

We plot the GWs signal in Fig.3 using 16 with the following parameters \[\mu = 0.5 ~,~ \kappa = 0.8 ~,~ \omega = 2 ~,\] \[w_c = 1.2 ~,~ w_{\rm rh} = 0 ~,~ H_c = 10^{16} {\rm GeV} ~,~ a_s = 0.4 ~.\]

The rest parameters are summarized in Tab. 2.

References↩︎