Gravitational-wave background in bouncing models from semi-classical, quantum and string gravity

Abstract

We study the primordial spectra and the gravitational-wave background (GWB) of three models of semi-classical, quantum or string gravity where the big bang is replaced by a bounce and the tensor spectrum is blue-tilted: ekpyrotic universe with fast-rolling Galileons, string-gas cosmology with Atick–Witten conjecture and pre-big-bang cosmology. We find that the ekpyrotic scenario does not produce a GWB amplitude detectable by present or third-generation interferometers, while the string-gas model is ruled out for producing too large a signal. In contrast, the GWB of the pre-big-bang scenario falls within the sensitivity window of both LISA and Einstein Telescope, where it takes the form of a single or a broken power law depending on the choice of parameters. The latter will be tightly constrained by both detectors.

1 Introduction

Gravitational waves (GWs) are one of the newest and most exciting windows into the workings of gravity. Thanks to the observation of GWs from binary systems by the LIGO-Virgo-KAGRA (LVK) network of ground-based interferometers [1, 2], and of a gravitational-wave background (GWB) of possibly astrophysical (but yet uncertain) origin in Pulsar Timing Arrays (PTAs) [3, 4, 5, 6, 7, 8], we have acquired a better understanding of the physics of solar-mass black holes and neutron stars as well as of the behaviour of the gravitational force near these compact objects and of its propagation across cosmological distances. Both LVK and third-generation instruments such as the Laser Interferometer Space Antenna (LISA) [9, 10, 11, 12], Einstein Telescope (ET) [13, 14] and DECIGO [15, 16, 17] will be able to further probe Einstein’s theory and its extensions to modified gravity and quantum gravity [18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33].

The detection of a relic primordial GWB [34, 35, 36] would be a momentous opportunity to look into the early universe and the gravitational interaction in extreme curvature or energy regimes. At the level of fundamental physics it is not easy to construct robust cosmological models embedded in realistic scenarios of modified or quantum gravity, and even less so to obtain one such model predicting an observable signal without invoking an ad hoc matter field dynamics. With a first, rapid scan of the literature we may find five candidates attempting to fulfil these characteristics, heterogeneous in terms of robustness and predictive power [29]: nonlocal Starobinsky inflation [37, 38], Brandenberger–Ho non-commutative inflation [39, 40], the new ekpyrotic universe [41, 42, 43], string-gas cosmology [44, 45, 46, 47] and multi-fractional inflation [48]. However, only the last four are possibly able to produce a detectable signal and only crossing the DECIGO sensitivity curve in the most optimistic cases. A sixth model of quantum gravity, nonlocal and non-minimally coupled with radiation, appeared afterwards with similar characteristics [32]. Recently, however, a more detailed exploration of the landscape of quantum and string cosmologies [10] singled out three more scenarios with a signal potentially reaching LISA and the Einstein Telescope:

• A bouncing model with a slow ekpyrotic contraction phase sustained by fast rolling Galileons with a non-canonical kinetic term and where perturbations are sourced by a $U(1)$ gauge field [49, 50, 51, 52, 53]. This particle content is the main difference with respect to the new ekpyrotic scenario of [41, 42, 43]. While the original ekpyrotic scenario was inspired by string-theory concepts, the current model is based on effective quantum field theory, and is not necessarily tied to a specific quantum-gravity theory.

• A bouncing model where the contracting phase is dominated by a string gas behaving like a stiff fluid and evolving according to Einstein’s gravity [54]. The main difference with respect to the string-gas cosmology of [44, 45, 46, 47] is that, while the latter model is based on the behaviour of closed-string modes below the Hagedorn temperature, in [54] the string thermodynamics was studied above the Hagedorn temperature. This implies that the free energy of the strings grows with the temperature more slowly than for ordinary radiation.

• A pre-big-bang model evolving from the string perturbative vacuum, proposed long ago [55, 56, 57, 58] on the grounds of the string cosmological equations, which enjoy T-duality [59] and may be characterized by a non-singular bounce thanks to all-orders (higher-curvature) $\alpha^{\prime}$ corrections [60, 61].

The first two models are phenomenological because they are not fully embedded in any high-energy quantum theory of gravity. The third one is based on the modified gravitational dynamics uniquely fixed, in principle, by the string unification of all fundamental interactions, at all energy scales including their quantum limit; however, it also contains phenomenological aspects concerning the (presently unknown) details of the dilaton dynamics in the strong coupling regime. In any case, all three models above are among the very few in quantum gravity possibly able to produce a GWB, arising from the vacuum fluctuations of the metric tensor, with high enough amplitude in the frequency range of third-generation detectors. They have in common an initial phase of growing curvature scale (described by a contracting kinematics with the metric of the Einstein frame), preceding the final decelerated expansion and passing through a non-singular bounce of spacetime curvature. The presence of accelerated contraction ($\dot{a}<0$, $\ddot{a}<0$, $\dot{H}<0$) or of super-inflationary expansion ($\dot{a}>0$, $\ddot{a}>0$, $\dot{H}>0$) in different metric frames is responsible for a blue tilt in the primordial tensor spectrum.

The first two models have been mainly explored at the level of the primordial tensor and scalar spectra, while for the third an approximate GWB profile is known (first computed in [62, 63, 64, 65] and recently discussed in [58]). The GWB of none of them, however, has been studied systematically so far, and the question of whether their signal can reach LISA, ET or DECIGO remains open. Also, one may wonder whether the common characteristic of having a contracting phase would produce a unique type of signal. It is the purpose of this paper to give an answer to these questions and to complement the analysis of [29] done for other models of quantum gravity with a blue-tilted primordial tensor spectrum. We find that:

• The GWB amplitude of the ekpyrotic model with Galileons is highly suppressed and unobservable, contrary to the new ekpyrotic scenario of [41, 42, 43];

• The string-gas cosmological model following the Atick–Witten conjecture is ruled out because its signal is too high and violates current bounds, contrary to string-gas models not adopting this conjecture [44, 45, 46, 47];

• Within the theoretically allowed parameter space, the GWB of the pre-big-bang model reaches the LISA and ET observational window while respecting present bounds. In all plots, we use the ET-D sensitivity curve for a three-year observation run [66, 67].

To the best of our knowledge, the pre-big-bang model and those studied in [29] are the only ones motivated by quantum gravity and generating a detectable GWB directly from the primordial tensor sector. Recently, another bouncing model was proposed where curvature perturbations evolving through a bounce can trigger the formation of primordial black holes and also induce a GWB signal with high amplitude crossing also the LISA and ET windows [68]. However, we do not consider scalar-induced GWBs here.

The paper is organized as follows. Basic expressions connecting the primordial tensor spectrum and the GWB are recalled in section 1.1. The primordial spectra and the GWB of the three models above are studied, respectively, in sections 2, 3 and 4. Conclusions are in section 5. Technical material is relegated to appendices.

1.1 Basic formulæ

Primordial GWs can be described by a small set of quantities and observables, independently of the underlying model. We denote the primordial tensor and scalar spectra as, respectively, $\mathcal{P}_{\rm t}(k)$ and $\mathcal{P}_{\rm s}(k)$, where $k$ is the comoving wave-number. From here, one calculates the tensor-to-scalar ratio at any given pivot scale $k=k_{*}$,

$$r\coloneqq\frac{\mathcal{P}_{\rm t}}{\mathcal{P}_{\rm s}}\,,$$ (1.1)

as well as the tensor and scalar spectral indices

$$n_{\rm t}\coloneqq\frac{d\ln\mathcal{P}_{\rm t}(k)}{d\ln k}\,,\qquad n_{\rm s}-1\coloneqq\frac{d\ln\mathcal{P}_{\rm s}(k)}{d\ln k}\,.$$ (1.2)

The current estimate for the scalar amplitude is $\mathcal{P}_{\rm s}(k_{*})\approx 2.1\times 10^{-9}$ at the cosmic microwave background (CMB) scale $k_{*}=0.05\,{\rm Mpc}^{-1}$ [69], while $n_{\rm s}=0.9649\pm 0.0042$ at 68% confidence level at the same pivot scale, assuming $dn_{\rm s}/d\ln k=0$ [70, 71]. The upper bound on (1.1) is $r(k_{*})<0.036$ at the same scale at $95\%$ confidence level [72].

In general, the amplitude of the GWB is defined as

$$\Omega_{\textsc{gw}}(k,\tau_{0})=\frac{1}{\rho_{\rm crit}(\tau_{0})}\frac{d\rho_{k}(\tau_{0})}{d\ln k}\,,$$ (1.3)

where $\tau_{0}$ is the present value of the conformal time $\tau$, defined in terms of the cosmic time $t$ as $\tau\coloneqq\int dt/a(t)$, where $a(t)$ is the scale factor, $\rho_{\rm crit}=3M_{\rm Pl}^{2}H^{2}$ is the critical density, $M_{\rm Pl}^{2}=(8\pi G)^{-1}$ is the reduced Planck mass and $\rho_{k}(\tau_{0})$ is the energy density of the $k$-th Fourier mode of tensor perturbations amplified by the given model of the early universe and evaluated at the present time $\tau_{0}$. For GWs generated by tensor perturbations, the GWB spectral shape can be recast as

$$\Omega_{\textsc{gw}}(f)=\frac{k^{2}}{12a_{0}^{2}H_{0}^{2}}\mathcal{P}_{\rm t}(k)\,\mathcal{T}^{2}(k,\tau_{0})\Big{|}_{k=2\pi f}\,,$$ (1.4)

where $f=k/(2\pi)$ is the GW frequency measured in Hz, $a_{0}$ is the scale factor today ($a_{0}=1$), $H_{0}$ is the value of the Hubble parameter today and $\mathcal{T}(k,\tau_{0})$ is the transfer function encoding how the primordial spectrum evolved after horizon crossing until today [73, 74, 75]. The expressions below are valid for any model where observable perturbations have re-entered the horizon in the past, either because they were originally generated inside the horizon and where later expelled out (e.g., by inflation), or because they were generated outside the horizon in the first place. In the case of instantaneous reheating and a quick transition to a radiation-dominated phase, which applies to all the models discussed below, one has [76]

$$\mathcal{T}^{2}(k,\tau_{0})=\Omega_{{\rm m}0}^{2}\left[\frac{g_{\ast}(T_{\rm in})}{g_{\ast 0}}\right]\left[\frac{g_{\ast s0}}{g_{\ast s}(T_{\rm in})}\right]^{4/3}\left[\frac{3j_{1}(k\tau_{0})}{k\tau_{0}}\right]^{2}\mathcal{T}_{\rm eq}^{2}(k)\,,$$ (1.5)

where $j_{1}$ is the first spherical Bessel function, the fitting function $g_{\ast}(T_{\rm in})$ is [77]

$$g_{\ast}[T_{\rm in}(k)] = g_{\ast 0}\frac{A_{1}+\tanh\left\{-2.5\log_{10}\left[\frac{k/(2\pi)}{2.5\times 10^{-12}\leavevmode\nobreak\ {\rm Hz}}\right]\right\}}{A_{1}+1}$$ (1.7)

$$\times\frac{A_{2}+\tanh\left\{-2.0\log_{10}\left[\frac{k/(2\pi)}{6.0\times 10^{-9}\leavevmode\nobreak\ {\rm Hz}}\right]\right\}}{A_{2}+1}\,,$$

with $g_{\ast 0}=3.36$, $A_{1}=(-1-10.75/g_{\ast 0})/(-1+10.75/g_{\ast 0})$, $A_{2}=(-1-g_{\rm max}/10.75)/(-1+g_{\rm max}/10.75)$, $g_{\rm max}=106.75$ for a Standard-Model particle content and the function $g_{\ast s}(T_{\rm in})$ is the same as (1.7) upon replacing $g_{\ast 0}\to g_{\ast s0}=3.91$. The fitting function at the end of (1.5) is [73]

$$\mathcal{T}_{\rm eq}^{2}(k)=1+1.57\frac{k}{k_{\rm eq}}+3.42\left(\frac{k}{k_{\rm eq}}\right)^{2}\,,\qquad k_{\rm eq}=7.1\times 10^{-2}\,\Omega_{\rm m}h^{2}\leavevmode\nobreak\ {\rm Mpc}^{-1},$$ (1.8)

$k_{\rm eq}$ being the comoving wave-number at radiation-matter equality, so that $f_{\rm eq}\approx 9.9\times 10^{-17}$ $(\Omega_{\rm m}h^{2}/0.143)\,{\rm Hz}$, where $\Omega_{\rm m}=0.3153$ is the matter-energy density and $h$ is the dimensionless Hubble parameter, which we take at its CMB-scale value $h=0.6736$ [69].

2 Ekpyrotic universe with fast-rolling Galileons

In the article [53], a non-singular bounce model was proposed as an alternative model for structure formation in the early Universe. Non-singular bounce models often predict a blue-tilted scalar spectrum contrary to CMB observations. The difference in the spectral index was shown to be remedied by the inclusion of a gauge field [52, 51]; the scalar spectrum matches current observations and the model predicts a tensor-to-scalar ratio with values below the upper bound [53]

$$r\lesssim 10^{-2}\,.$$ (2.1)

The model involves Galileons with a non-canonical kinetic term and a coupling with a $U(1)$ gauge field. The action is

$$S=M_{\rm Pl}^{2}\int d^{4}x\sqrt{|g|}\left[\frac{R}{2}-K(\phi,X)+G(\phi,X)\Box\phi-I^{2}(\phi)\left(\frac{1}{4}F^{\mu\nu}F_{\mu\nu}-\frac{\delta}{4}\tilde{F}^{\mu\nu}F_{\mu\nu}\right)\right],$$ (2.2)

where $M_{\rm Pl}$ is the reduced Planck mass, $\phi$ is the Galileon field, $A_{\mu}$ is the $U(1)$ gauge vector, $F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$, $\tilde{F}^{\mu\nu}=\frac{1}{2}\epsilon^{\mu\nu\rho\sigma}F_{\rho\sigma}$ and $\delta>0$ is a coupling constant. The functions $K$ and $G$ in (2.2) are

$$K(\phi,X)=[1-g(\phi)]X+\beta X^{2}-V(\phi)\,,$$ (2.3)

$$G(\phi,X)=\gamma X\,,$$ (2.4)

$$g(\phi)=\frac{2g_{0}}{e^{-\sqrt{\frac{2}{p}}\phi}+e^{b_{g}\sqrt{\frac{2}{p}}\phi}}\,,\;\;\;\;\;V(\phi)=-\frac{2V_{0}}{e^{-\sqrt{\frac{2}{q}}\phi}+e^{b_{V}\sqrt{\frac{2}{q}}\phi}}\,,$$ (2.5)

where $X=-\frac{1}{2}\partial_{\mu}\phi\partial^{\mu}\phi$ and $\beta$, $\gamma$ and $g_{0}$ are parameters. The background dynamics of the Universe is determined by a single scalar field $\phi$ with a non-trivial kinetic term typical of Galileons. This model gives rise to a non-singular bounce and the potential $V$ is chosen in such a way as to obtain ekpyrotic contraction away from the bounce. The Universe starts at $\phi\ll-1$, far away from the bounce, with a slow ekpyrotic contraction. As $\phi$ accelerates towards $\phi$ = 0, the value of $g$ increases. Since $g(0)>1$ (which we require), at some point in time, $g$ exceeds the critical value $g=1$ and the sign of the kinetic term $X$ in (2.3) becomes negative, giving rise to a phase of ghost condensation coinciding with the bouncing phase. This is the region in field space where the null energy condition is violated, which in turn triggers the bounce at $\phi=0$. The Universe continues to roll to positive larger values of the field $\phi>\phi_{B+}$, after which the Lagrangian regains a canonical form and the Universe enters an era of kinetic energy domination. $I$, the coupling of the scalar field to the gauge field, is

$$I(\phi)=\frac{1}{1+e^{-a_{1}n(\phi-\phi_{B-})}}\,,$$ (2.6)

where $a_{1}=1/\sqrt{2q}$ and $\phi_{B-}$ is value of field $\phi$ at the beginning of the bounce. $I(\phi)$ is defined such that during the regime of ekpyrosis, i.e., for large and negative $\phi$,

$$I(\phi)\simeq e^{a_{1}n\phi}.$$ (2.7)

The values $n=2$ or $n=-1$ lead to scale-invariant sourced perturbations [52]. Introducing gauge fields can source second-order inhomogeneities that dominate over the blue-tilted first-order perturbations at CMB scales. This leads to sourced and unsourced perturbations that are linearly independent since the equation for perturbations is linear and creation/annihilation operators of sourced and unsourced fluctuations are uncorrelated. The total power spectrum is of the form

$$\mathcal{P}^{\rm tot}=\mathcal{P}^{v}+\mathcal{P}^{s}$$ (2.8)

where $\mathcal{P}^{v}$ is the unsourced or vacuum spectrum and $\mathcal{P}^{s}$ is the sourced spectrum. The vacuum scalar ($\mathcal{P}_{\rm s}^{v}$) and tensor ($\mathcal{P}_{\rm t}^{v}$) spectra as well as the sourced scalar ($\mathcal{P}_{\rm s}^{s}$) and tensor ($\mathcal{P}_{\rm t}^{s}$) spectra are [53]

$$\begin{split}\mathcal{P}^{v}_{\rm t}(k)&=\left(\frac{k}{H_{B^{-}}}\right)^{\frac{6(1+w_{\rm ekp})}{1+3w_{\rm ekp}}}\frac{\gamma_{E}^{2}H_{B^{-}}^{2}}{2\pi^{2}M_{\rm Pl}^{2}}\,,\\ \mathcal{P}^{v}_{\rm s}(k)&=\left(\frac{k}{H_{B^{-}}}\right)^{\frac{6(1+w_{ekp})}{1+3w_{ekp}}}\frac{\gamma_{E}^{2}H_{B^{-}}^{2}}{48\pi^{2}M_{\rm Pl}^{2}}F_{\zeta}^{2}\,,\\ \mathcal{P}^{s}_{\rm t}(k)&=2.8\frac{9e^{4\pi\xi}}{8\pi q^{4}\xi^{6}}\left(\frac{H_{B-}}{M_{\rm Pl}}\right)^{4}\left(\frac{k}{H_{B-}}\right)^{n_{\rm t}}\,,\\ \mathcal{P}^{s}_{\rm s}(k)&=2.8\frac{225e^{4\pi\xi}}{32\pi q^{4}\xi^{6}}\left(\frac{H_{B-}}{M_{\rm Pl}}\right)^{4}\left(\frac{k}{H_{B-}}\right)^{n_{\rm s}-1}F_{\zeta}^{2}\,,\\ \end{split}$$ (2.9)

where $\gamma_{E}$ is the Euler–Mascheroni constant, $q$ is related to the equation of state of the ekpyrotic phase by $w_{\rm ekp}=-1+\frac{2}{3q}$, $\xi$ is related to the strength of gauge coupling by $\xi=2\pi\delta$ and $H_{B-}$ is the Hubble parameter at the end of the bounce determined by the scale of the bounce. $n_{\rm s}$ and $n_{\rm t}$ are the spectral indices of scalar and tensor perturbations, respectively. By choosing $n=-2.01$, it is possible to obtain a scalar spectral index $n_{\rm s}=0.96$, a red tilted scalar power spectrum in agreement with observations. The tensor-to-scalar ratio $r$ is determined by the factor $F_{\zeta}$, governing amplification across the bounce. Here lies a main qualitative difference between the different spectra. The scalar spectra are amplified across the bounce while the tensor ones are unchanged. Since the sourced spectrum is significantly larger compared to the vacuum spectrum, $\mathcal{P}^{\rm tot}\simeq\mathcal{P}^{s}$ and

$$r\simeq\frac{\mathcal{P}_{\rm t}^{s}}{\mathcal{P}_{\rm s}^{s}}=\frac{4}{25}\left(\frac{1}{F_{\zeta}}\right)^{2}.$$ (2.10)

The sourced tensor spectrum has the same red-tilted spectral index as the scalar spectrum and is unobservable by current and third-generation interferometers . However, the vacuum perturbations are blue-tilted, and although tensor vacuum perturbations are small compared to the sourced spectrum at CMB scales, the blue-tilted spectrum may be observable at higher frequencies.

After the bounce, kinetic domination takes place. The gauge field in our model can be either a putative gauge field, or the actual $U(1)$ of electromagnetism. The different options will imply different reheating epochs and may modify the predicted tensor spectrum. We now examine the possibility of the gauge field sourcing perturbations as electromagnetic radiation. It is possible to calculate the scale factor at which the Universe ends a reheating phase with this assumption. The Universe becomes radiation dominated when the energy density of the scalar field responsible for the bounce is smaller than the energy density of the gauge field. The energy density of the scalar field $\phi$ at the end of the bounce phase is approximately

$$\rho_{\phi}(\tau_{B+})\simeq\frac{1}{2}\dot{\phi}^{2}(\tau_{B+})[1-g_{0}+3\beta\dot{\phi}^{2}(\tau_{B+})]\simeq\frac{g_{0}-1}{3\beta}e^{-2\frac{\tau_{B+}}{T}}\left(1-e^{-2\frac{\tau_{B+}}{T}}\right)\,,$$ (2.11)

where $\tau_{B+}$ and $\tau_{B-}$ are the conformal time at the beginning and at the end of the bounce phase, respectively, and $T$ is a quarter of the duration of bounce phase. During the phase of kinetic domination, the energy density decays as $a^{-6}$. Thus, the energy density of the scalar field during this phase is

$$\rho_{\phi}=\frac{g_{0}-1}{3\beta}e^{-2\frac{\tau_{B+}}{T}}\left(1-e^{-2\frac{\tau_{B+}}{T}}\right)\left(\frac{a_{B}}{a}\right)^{6}.$$ (2.12)

The gauge field energy density attains the maximum

$$\rho_{A,\tau_{B-}}=D(n)\frac{e^{2\pi\xi}}{\xi^{3}\tau_{B-}^{4}}\,,\qquad D(n)=\frac{1}{4\pi^{2}}\frac{(n-1)^{2}\Gamma(2n-1)}{2^{2n+1}(n-2)\pi}\,,$$ (2.13)

at the bounce, where the expression for $D(n)$ holds for $n$ close to 2. After the bounce, the energy density decays as $a^{-4}$ and the gauge field energy density is

$$\rho_{A,\tau}=D(n)\frac{e^{2\pi\xi}}{\xi^{3}\tau_{B-}^{4}}\left(\frac{a_{B}}{a}\right)^{4}\,.$$ (2.14)

From the expressions for energy densities, one can calculate the scale factor corresponding to which the Universe ends the reheating phase and find that reheating will last for only a few e-folds.

In the derivation of primordial spectra in (LABEL:e:finspec1), we have assumed the scale factor at the bounce to be one. However, the transfer function is derived with normalization $a_{0}=1$. For consistency, we reparametrize $a_{0}$ as

$$a_{0}=a_{B}\frac{a_{\rm r}}{a_{B+}}\frac{a_{B+}}{a_{B}}\frac{a_{0}}{a_{\rm r}}$$ (2.15)

where $a_{\rm r}$ is the scale factor at the onset of radiation domination. The reheating phase for this model is dominated by the kinetic term from $t_{B+}$ to the beginning of radiation domination, implying

$$\frac{a_{\rm r}}{a_{B+}}=\left(\frac{t_{\rm r}}{t_{B+}}\right)^{\frac{1}{3}}\simeq\left(\frac{H_{B+}}{H_{\rm r}}\right)^{\frac{1}{3}}.$$ (2.16)

From (2.15) and (2.16), we obtain

$$\frac{a_{0}}{a_{B}}=\frac{34\times 10^{13}}{87\,\text{GeV}}H_{B+}^{\frac{1}{3}}H_{\rm r}^{\frac{1}{6}}M_{\rm Pl}^{\frac{1}{2}}\,.$$ (2.17)

Instantaneous reheating implies $H_{\rm r}\simeq H_{B+}$. In our case, the value of $H_{B+}$ is determined by the theory of sourced perturbations and is dependent on the parameters of the theory. In our discussion, we set parameters for sourced perturbations to values optimal for obtaining a reasonable tensor-to-scalar ratio and a red spectral tilt. This value is around $H_{B+}\simeq 10^{-5}M_{\rm Pl}$. We notice that, for the case where the gauge field responsible for sourcing gravitational waves behaves as radiation, the reheating phase is short, $cN\sim 3-4$, and the spectrum of GWs is similar to that of the case where reheating is instantaneous. For models with instant reheating, the transfer function is described in (1.5) and present-day GWB is given by (1.4). We present our results in Fig. 1, where we used vacuum tensor spectrum from (LABEL:e:finspec1) assuming $w_{\rm ekp}\gg 1$ which gives

$$\mathcal{P}_{\rm t}^{v}(k)\simeq\frac{\gamma_{E}^{2}k^{2}}{2\pi^{2}M_{\rm Pl}^{2}}\,.$$ (2.18)

If, instead, we assume that the gauge field responsible for sourced perturbations is not radiation and is something more exotic such as $U(1)$ axions or dark radiation, then it is possible to have an extended reheating period where the Universe is dominated by a kinetic term. We have also examined this possibility. The transfer function in the presence of a reheating phase is slightly modified and is well known [76]. Figure 1 also shows the GWB in the presence of a prolonged reheating phase, where we assumed a large reheating period lasting till the onset of big-bang nucleosynthesis (BBN).

Figure 1: GWB generated by the primordial vacuum tensor spectrum of the ekpyrotic model with fast-rolling Galileons for the instantaneous reheating scenario with $H_{B+}=H_{\rm r}=10^{-5}\,M_{\rm Pl}$ (solid thin curve) and the prolonged reheating scenario with $H_{B+}=10^{-5}\,M_{\rm Pl}$ and $H_{\rm r}=10^{-43}\,M_{\rm Pl}$ (solid thick curve), compared with the sensitivity curves (dashed) of LVK, SKA, LISA, ET and DECIGO. Also shown are the CMB and BBN exclusion (where the CMB upper bound is given thickness to improve visibility in the plot) and, as a thick red segment, the signal detected by PTAs. The bottom plot captures the high-amplitude part of the GWB.

This analysis shows that the predicted GWB will not be observable with the current sensitivity of GW observatories, nor by third-generation interferometers. Also, the long reheating scenario is ruled out because it violates the BBN upper bound. This constraint, however, can be bypassed assuming fewer e-folds of contraction. Considering a contraction period of $cN=55$ instead of 60 e-folds, which is the minimum conservatively assumed to isotropize the universe, is enough to avoid violations of the BBN bound [78].

3 String-gas cosmology with Atick–Witten conjecture

One of the characteristic features of string theory is the existence of the Hagedorn phase at temperatures close to the string scale $M_{\rm s}$, where the energy is not dominated by the massless modes but rather by the most massive string states, leading to a pressureless fluid [79, 80, 81, 82]. In fact, a canonical description of the thermal phase indicates a limiting Hagedorn temperature [79, 80]. However, it was also argued that the limiting temperature only corresponds to the emergence of a thermal tachyonic mode, making the description of the system in terms of fundamental string excitations invalid [83, 84]. Atick and Witten [85] conjectured that, at temperatures larger than the Hagedorn temperature, the free energy grows much more slowly than in conventional field theories. Therefore, the system represents fewer degrees of freedom than expected from the zero-temperature string spectrum or from point-like particle field theories.

According to [85], the partition function only grows as $T^{2}$, hence the authors modeled the pressure with quadratic, linear and logarithmic terms in temperature. However, for convenience and transparency we may as well model the pressure in the form [54]

$$p(T)=M_{s}^{4}\left[\left(\frac{T}{M_{s}}\right)^{2}+c_{1}\left(\frac{T}{M_{s}}\right)^{1+\alpha}\right]\,,\qquad|\alpha|\ll 1\,.$$ (3.1)

where $\alpha$ is a real constant.

The possibility of thermal fluctuations being the origin of small inhomogeneities and anisotropies in the CMB was already proposed by Peebles [86]. The fluid fluctuations may arise naturally from two different sources. There might be fluctuations in the energy density and the associated temperature. And, even if there is a unique temperature in a given volume, there are fluctuations in energy within this volume due to the very statistical nature of thermal physics. This could also potentially seed primordial fluctuations; see, for instance, [87, 88, 89, 90] and references therein.

The statistical fluctuations in the energy inside a given volume $V=L^{3}$ is given by

$$\langle\Delta E\rangle_{L}^{2} \coloneqq \langle E^{2}\rangle-\langle E\rangle^{2}=\frac{\partial^{2}\ln Z}{\partial\beta^{2}}=T^{2}C_{L}$$

$$\Longrightarrow\quad\langle\delta\rho^{2}\rangle_{L} = \frac{T^{2}C_{V}}{L^{6}}=\frac{T^{2}}{L^{3}}\frac{\partial\rho}{\partial T}\,,$$ (3.2)

where $C_{L}$ is the heat capacity of the thermal system for the volume $L^{3}$. These are classical random fluctuations that exist in any finite-temperature system as long as the fluid is in local thermal equilibrium. Hence, as such there is no need for any seed quantum fluctuations here. Therefore, the power spectrum for the seed fluctuations could then be sourced by these sub-Hubble wavelengths. Once the modes become super-Hubble, thermal correlations over the relevant physical wavelengths can no longer be maintained, then the coupled metric and matter fluctuations can evolve according to the usual general-relativistic hydrodynamical differential equations.

A precise understanding of how these statistical fluctuations get encoded in the curvature perturbation $\zeta$ at the Hubble crossing, was achieved in [92] for a general extensive thermodynamic fluid whose pressure $p$ can be an arbitrary function of the temperature $T$. The derived curvature power spectrum reads

$${\cal P}_{\zeta}=k^{3}\zeta_{k}^{2}=8\sqrt{3\pi^{3}}A^{2}(T_{k})\frac{T_{k}^{2}\rho_{k}^{\prime}}{M_{\rm Pl}^{3}\sqrt{\rho_{k}}}\,,$$ (3.3)

where a prime denotes differentiation with respect to $T$ and

$$A(T)=\frac{3(1+w)\Omega+2(3+R)}{6(1+w)\Omega}\,,\qquad R=-\frac{3}{2}\left[1+\frac{(1+w)\rho(2\rho^{\prime}+T\rho^{\prime\prime})}{T{\rho^{\prime}}^{2}}\right].$$ (3.4)

The subscript $k$ (which we are going to subsequently drop) refers to the fact that all these quantities have to be evaluated at horizon crossing, $H_{k}=k/a$. Note, that all the above functions of temperature can be calculated if we know $p(T)$, as the energy density is related rather straightforwardly to pressure:

$$\rho(T)=T\frac{dp(T)}{dT}-p(T)\,.$$ (3.5)

The primordial tensor and scalar spectra as functions of the temperature $T$ are [54]

$$\mathcal{P}_{\rm t}(k) = \frac{2}{\pi^{2}}\left[\frac{H(k)}{M_{\rm Pl}}\right]^{2},$$ (3.6)

$$\mathcal{P}_{\rm s}(k) = \frac{\sqrt{3\pi^{3}}c_{1}^{2}}{4}\left(\frac{M_{\rm s}}{M_{\rm Pl}}\right)^{3}\left[\frac{T(k)}{M_{\rm s}}\right]^{2\alpha},$$ (3.7)

where $M_{\rm s}$ is the string mass scale and $c_{1}$ is the constant coefficient in the $O(T^{1+\alpha})$ term in (3.1). By tuning the parameters $c_{1}$ and $M_{\rm s}$ in (3.7), one can easily make the tensor-to-scalar ratio

$$r=\frac{8}{\pi^{3}\sqrt{3\pi}c_{1}^{2}}\left(\frac{H}{T}\right)^{2}\frac{M_{\rm Pl}}{M_{\rm s}}\left(\frac{T}{M_{\rm s}}\right)^{2(1-\alpha)},$$ (3.8)

as large as the upper bound $r=0.036$ at $k_{*}=0.05\,{\rm Mpc}^{-1}$. This, however, turns out to be ruled out by observations, as we shall see presently.

At high temperature, the string gas behaves like a stiff fluid $p\simeq\rho$ and the gravitational setting is Einstein’s gravity, so that the standard Friedmann equations hold and

$$H^{2}\propto\rho\propto T^{2}\propto a^{-6}\,.$$ (3.9)

Therefore, from the horizon crossing relation $k=aH$ we have $k\propto a^{-2}$ and

$$\mathcal{P}_{\rm t}\propto a^{-6}\,,\qquad\mathcal{P}_{\rm s}\propto a^{-6\alpha}\,,\qquad\frac{d}{d\ln k}=-\frac{1}{2}\frac{d}{d\ln a}\,,$$ (3.10)

so that the tensor and scalar spectral indices (1.2) are constant and given by

$$n_{\rm t}=3\,,\qquad n_{\rm s}-1=3\alpha\,.$$ (3.11)

The tensor spectrum (3.6) can thus be written as

$$\mathcal{P}_{\rm t}(k)=r(k_{*})\,\mathcal{P}_{\rm s}(k_{*})\left(\frac{k}{k_{*}}\right)^{n_{\rm t}},$$ (3.12)

which we plugged into (1.4) to get the spectral shape shown in Fig. 2.

Figure 2: GWB for the bouncing model of string-gas cosmology with Atick–Witten conjecture with $r=0.036,10^{-18},5\times 10^{-22}$ (increasing thickness) at the CMB pivot scale (where the upper bound is given thickness to improve visibility in the plot), compared with the sensitivity curves (dashed) of LVK, SKA, LISA, ET and DECIGO. Also shown are the BBN exclusion region and, as a thick red segment, the signal detected by PTAs. The bottom plot captures the high-amplitude part of the GWB.

As the reader can see, the signal predicted by this model is too strong and is excluded by all present constraints. The problem is the high tensor index in (3.11), two orders of magnitude higher than the minimum tilt required for detection in DECIGO, ET and LISA ($n_{\rm t}\gtrsim 0.06\!-\!0.34$) [33]. At $r=0.036$, the curve hits the Gaia bound around $\Omega_{\textsc{gw}}(f=5\times 10^{-12}\,{\rm Hz})\approx 7\times 10^{-12}$ [93]. The tensor-to-scalar ratio must be reduced down to $r\sim 10^{-18}$ to cross the signal detected by PTA, with which it is incompatible since $n_{\rm t}\approx 2/3$ there [7]. To shift past the PTA-SKA curve, it must be $r<5\times 10^{-22}$, at which point the parameter space of the model becomes too compressed to yield a realistic scenario. In the meanwhile, the BBN bound $\Omega_{\textsc{gw}}<5\times 10^{-6}$ truncates the top of all these signals before they can reach the cut-off frequency, which is much larger than the range depicted in the plot [29].

4 Pre-big-bang cosmology

A possible example of strongly blue-tilted spectrum of primordial tensor perturbations, peaked at high frequencies, is provided by the so-called Pre-Big-Bang (PBB) scenario [55, 56], based on the scale-factor duality of the string cosmology equations [59]. Such symmetry is a peculiar property of string theory, and is a crucial ingredient not only to fix the slope of the primordial spectrum [62, 63, 64] but also, as confirmed by recent results [60, 61], to implement a smooth transition from the initial growing curvature (PBB) regime to the standard regime of decreasing curvature evolution. According to this model, the evolution of our Universe, characterized by decelerated expansion (at intermediate times), decreasing temperature and curvature, weak gravitational coupling, should be preceded in time by an almost specularly symmetric phase of accelerated expansion, increasing curvature, increasing density and temperature and growing coupling. Such a dual counterpart of the present one describes a “pre-big-bang” evolution from a flat, cold, empty initial state with negligible interactions to a final high-curvature, high-energy, explosive bounce, marking the transition to the more standard cosmological regime (see, e.g., [94] for a recent non-technical introduction, and [95] for a more detailed and complete discussion).

Here we recall the derivation of the associated GWB in a self-contained way. We shall introduce the spectral energy density of the relic GW radiation present today inside our cosmic horizon, and produced by a simple model of PBB scenario which satisfies all present observational constraints and depends on four constant parameters (see also [58]). Two of these parameters control the inflationary growth of the scale factor and of the string coupling in the high-energy regime preceding the bounce; the other ones control the beginning and the end of the axion-dominated phase occurring after the curvature bounce. We recall that the presence of a dust-like phase dominated by the oscillations of the Kalb–Ramond axion is in general needed in the PBB scenario to obtain (via the curvaton mechanism) a flat spectrum of scalar metric perturbations [96, 97].

At present, the first two parameters are largely arbitrary, while the other two may vary in a rather small range of values. We stress, however, that the dependence of the amplitude of the GWB on the full set of the above four parameters is given and discussed, for the first time, in this paper. We have neglected a possible further parameter, the effective propagation speed of tensor perturbations during the high curvature string phase: it may be important for the production of primordial black holes [98] but it seems to have small effects on the energy density of the GW spectrum. We thus postpone this study to a future paper.

Let us now compute the spectral energy density in critical units of the GWB, eq. (1.3) with $k_{0}\leqslant k\leqslant k_{1}$, in such a way that all modes $k$ of eq. (1.3) satisfy the condition $k\tau_{0}\leavevmode\nobreak\ \raise 1.72218pt\hbox{$>$}\kern-6.99997pt\lower 2.6694pt\hbox{$\sim$}\leavevmode\nobreak\ 1$ (i.e., they are all inside our present Hubble horizon, $k\leavevmode\nobreak\ \raise 1.72218pt\hbox{$>$}\kern-6.99997pt\lower 2.6694pt\hbox{$\sim$}\leavevmode\nobreak\ k_{0}=\tau_{0}^{-1}$). The highest mode $k_{1}=\tau_{1}^{-1}$ is the maximum amplified frequency crossing the horizon just at the end of the phase of PBB inflation. Higher frequency modes ($k\gg k_{1}$) can, in principle, be included into Eq. (1.3), but their amplitude is exponentially suppressed [99, 100, 101] and their contribution to $\Omega_{\textsc{gw}}$ is negligible. For the explicit computation of $\Omega_{\textsc{gw}}$, two comments must be made.

The first is that, even if we are working in general with a higher-dimensional spacetime manifold (an unavoidable choice in the string theory context), we are mainly interested in the tensor perturbations of the four-dimensional metric, $\delta g_{\mu\nu}=h_{\mu\nu}$, assuming that the extra spatial dimensions are today compactified with frozen dynamics. However, this does not mean that we are neglecting the possible effects of the higher-dimensional geometry during its initial, non-trivial evolution: indeed, all such higher-dimensional contributions will be included into the canonical equation which controls the dynamics of $h_{\mu\nu}$ [56, 58, 95].

The second point is that, as usual, we are interested in the contributions to $\Omega_{\textsc{gw}}$ arising from the cosmological amplification of the quantum vacuum fluctuations of the metric tensor. This implies that we can describe the amplification of tensor perturbations as a quantum (or semi-classical) field-theory process of production of pairs of gravitons from the initial vacuum state (see, e.g., [95]) and we can write, for each mode $k$, the differential energy density of the amplified perturbations as follows:

$$d\rho_{k}(\tau_{0})=2k\,\langle n_{k}(\tau_{0})\rangle\,\frac{d^{3}\bm{k}}{(2\pi)^{3}}=\frac{k^{4}}{\pi^{2}}\,\langle n_{k}(\tau_{0})\rangle\,d\ln k\,,$$ (4.1)

where $2$ is the number of polarization states and $\langle n_{k}(\tau_{0})\rangle$ the number density of produced gravitons at the final epoch $\tau_{0}$. The last equality follows from assuming an isotropic final distribution.

To obtain $\langle n_{k}\rangle$ and then $\Omega_{\textsc{gw}}$, we need now to solve the evolution equation for the Fourier component of the (Mukhanov-Sasaki) canonical variable, $v_{k}(\tau)$, defined by putting in canonical form the quadratic action for the tensor fluctuation mode $h_{k}$ [102]:

$$v_{k}^{\prime\prime}+\left(k^{2}-\frac{\xi^{\prime\prime}}{\xi}\right)v_{k}=0\,.$$ (4.2)

Here a prime denotes differentiation concerning the conformal time, $v_{k}=\xi h_{k}$ and $\xi(\tau)$ is the so-called “pump field” which controls, according to the above equation, the dynamics of the perturbation modes in the given background.

In the model we are considering, the background may be approximated as a sequence of various cosmic phases, and in each of them the pump field $\xi$ is characterized by a simple power-law behaviour. In particular, in the initial PBB regime $-\infty<\tau\leqslant-\tau_{1}$, starting asymptotically from the string perturbative vacuum and approaching the curvature bounce at $\tau=-\tau_{1}$, we have two different phases with the following (canonically normalized) pump field behaviour [58]:

$$\xi\sim\frac{M_{\rm Pl}}{\sqrt{2}}(-\tau)^{1/2},\qquad\tau<-\tau_{s}\leavevmode\nobreak\ ;\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \xi\sim\frac{M_{\rm Pl}}{\sqrt{2}}(-\tau)^{\beta-1},\!\!\!\qquad-\tau_{s}<\tau<-\tau_{1}\,.$$ (4.3)

The parameter $\beta$ describes the high-energy growth of the dilaton and the dynamics of the internal dimensions, while the time scale $\tau_{s}$ is a free parameter which marks the transition from the low-energy initial phase to a possible late-time attractor, where spacetime curvature stays frozen at the value controlled by the fundamental string mass scale [103]. In both phases, the evolution of the pump field takes into account not only, as usual, the inflationary growth of the scale factor but also the additional string-theory effects [56, 95], such as the dynamics of the extra dimensions and the growth of the string coupling controlled by the scalar dilaton field (see also [103, 65]).

In the subsequent post-bouncing regime $-\tau_{1}<\tau\leqslant\tau_{0}$, the cosmic evolution is decelerated and we may have in principle three phases with the following pump-field behaviour [58]:

$$\xi\sim\frac{M_{\rm Pl}}{\sqrt{2}}\,\tau,\leavevmode\nobreak\ \leavevmode\nobreak\ -\tau_{1}<\tau<\tau_{\sigma};\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \xi\sim\frac{M_{\rm Pl}}{\sqrt{2}}\,\tau^{2},\leavevmode\nobreak\ \leavevmode\nobreak\ \tau_{\sigma}<\tau<\tau_{d};\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \xi\sim\frac{M_{\rm Pl}}{\sqrt{2}}\,\tau,\leavevmode\nobreak\ \leavevmode\nobreak\ \tau_{d}<\tau<\tau_{\rm eq}.$$ (4.4)

Here we are assuming that the extra-dimensional geometry and the string coupling (i.e., the dilaton) are frozen after the bounce, so that the pump field simply coincides with the scale factor. Here, again, we have two free parameters: the time scale $\tau_{\sigma}$, marking the beginning of the dust-like phase dominated by the axion oscillations, and the time scale $\tau_{d}$, marking the epoch of axion decay associated with the conventional reheating (source of the CMB radiation that we are presently observing) and corresponding to the beginning of the standard post-big-bang evolution.

As discussed in previous papers [58, 96, 97], instead of $\tau_{\sigma}$ and $\tau_{d}$ we can conveniently use as parameters the corresponding curvature scales $H_{\sigma}\coloneqq H(\tau_{\sigma})$ and $H_{d}\coloneqq H(\tau_{d})$ which can be expressed in terms of the (unknown) mass $m$ of the Kalb–Ramond axion and of its initial amplitude $\sigma_{i}$ after the bounce:

$$H_{\sigma}\simeq m\left(\frac{\sigma_{i}}{M_{\rm Pl}}\right)^{4},\qquad\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ H_{d}\simeq m\left(\frac{m}{M_{\rm Pl}}\right)^{2}.$$ (4.5)

For a consistent model, the allowed values of the parameters $m$ and $\sigma_{i}$ must satisfy the scale hierarchy $M_{\rm Pl}\leavevmode\nobreak\ \raise 1.72218pt\hbox{$>$}\kern-6.99997pt\lower 2.6694pt\hbox{$\sim$}\leavevmode\nobreak\ H_{1}\leavevmode\nobreak\ \raise 1.72218pt\hbox{$>$}\kern-6.99997pt\lower 2.6694pt\hbox{$\sim$}\leavevmode\nobreak\ H_{\sigma}\leavevmode\nobreak\ \raise 1.72218pt\hbox{$>$}\kern-6.99997pt\lower 2.6694pt\hbox{$\sim$}\leavevmode\nobreak\ H_{d}$. In addition, we have the obvious condition $H_{d}>H_{N}$, where $H_{N}\sim(1\,{\rm Mev})^{2}/M_{\rm Pl}$ is the nucleosynthesis scale of the standard cosmological scenario.

Given the full model of background evolution from the initial state at $\tau\rightarrow-\infty$ down to the present epoch $\tau_{0}$, and given the power-law behaviour of $\xi$ in the various phases, we can now work in the so-called “sudden approximation” [100] by imposing on the pump field to be continuous at the transitions scales, and solving, in each phase, the canonical equation (4.2). We recall that, in general, for a pump field given by $\xi=(M_{\rm Pl}/\sqrt{2})|\tau/\tau_{1}|^{\alpha}$, the exact solution for $h_{k}$ obtained from (4.2) can be written in terms of the first- and second-kind Hankel functions, $H_{\nu}^{(1)}$, $H_{\nu}^{(2)}$, as

$$h_{k}(\tau)=\frac{v_{k}}{\xi}=\left(\frac{2\tau_{1}}{Mp^{2}}\right)^{1/2}\left|\frac{\tau}{\tau_{1}}\right|^{\nu}\left[A_{+}(k)H_{\nu}^{(2)}(k\tau)+A_{-}(k)H_{\nu}^{(1)}(k\tau)\right],$$ (4.6)

where $\nu=-\alpha+1/2$. The complete solution for $h_{k}(\tau)$, describing its evolution from $-\infty$ to $\tau_{0}$, is then obtained by solving the canonical equation in the various phases and matching $h_{k}$ and $h_{k}^{\prime}$ at the transition scales.

In our model, in particular, we have four transitions (at $\tau_{s}$, $\tau_{1}$, $\tau_{\sigma}$, $\tau_{d}$), which means five different phases of background evolution (see Eqs. (4.3) and (4.4)), which implies five different solutions like (4.6), and thus ten different coefficients $A_{\pm}(k)$ to be determined at the various epochs. The continuity of $h_{k}$ and $h_{k}^{\prime}$ only gives eight conditions. The two remaining conditions are obtained by imposing on the canonical variable to initially describe a positive-frequency mode normalized to the quantum fluctuations of the Bunch–Davies vacuum, namely, by imposing $v_{k}=(1/\sqrt{2k})\exp(-ik\tau)$ for $\tau\rightarrow-\infty$: this implies (using the large-argument limit of the Hankel functions [95]) $A_{-}=0$ and $A_{+}=\sqrt{\pi/4}$ for the solution describing perturbations in the initial regime $\tau\rightarrow-\infty$.

With such a canonical normalization, the sought value of the number density $\langle n_{k}(\tau_{0})\rangle$ of produced gravitons is then automatically obtained from the coefficient $A_{-}(k)$ of the first-kind Hankel function describing the perturbation mode $h_{k}(\tau)$ in the final regime $\tau\rightarrow+\infty$ (actually, for our purpose, in the limit $\tau\rightarrow\tau_{0}$). More precisely, one finds (see, e.g., [95])

$$\langle n_{k}(\tau_{0})\rangle=\frac{4}{\pi}\left|A_{-}(k)\right|_{\tau=\tau_{0}}.$$ (4.7)

By performing the above computation and varying $k$ in the allowed frequency range we find that there are, in principle, four different branches of the energy density spectrum (1.3), (4.1), depending on the epochs of horizon crossing of the various modes. Noting that the axion-dominated phase is expected to occur soon after the bounce, in order to have a short duration with respect to the preceding high-curvature string phase (namely, $\tau_{d}/\tau_{\sigma}\ll\tau_{1}/\tau_{s}$) we may consistently assume that all modes re-entering the horizon before, during, or soon after the axion phase are amplified modes leaving the horizon during the high-curvature string phase. This means, in other words, that we can work with the following hierarchy of wave-number scales: $k_{1}\geqslant k_{\sigma}\geqslant k_{d}>k_{s}$, where $k_{i}\coloneqq\tau_{i}^{-1}$ is the limiting frequency of a mode crossing the horizon at the transition epoch $\tau_{i}$.

To obtain the explicit (parameter-dependent) form of the GWB, it is useful to first compute the frequency ratios of the four scales $k_{i}$. We note that in our model there are two phases of radiation-dominated evolution (i.e., $a\sim\tau\sim H^{-1/2}$) for $-\tau_{1}<\tau<\tau_{\sigma}$ and for $\tau>\tau_{d}$, one phase of matter-dominated evolution (i.e., $a\sim\tau^{2}\sim H^{-2/3}$) for $\tau_{\sigma}<\tau<\tau_{d}$, and one phase where the string frame scale factor undergoes a de Sitter-like evolution (i.e., $a\sim|\tau|^{-1}$) for $-\tau_{s}<\tau<-\tau_{1}$. By defining the convenient parameters

$$z_{s}\coloneqq\frac{\tau_{s}}{\tau_{1}}=\frac{k_{1}}{k_{s}},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ z_{\sigma}\coloneqq\frac{\tau_{\sigma}}{\tau_{1}}=\frac{k_{1}}{k_{\sigma}},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ z_{d}\coloneqq\frac{\tau_{d}}{\tau_{1}}=\frac{k_{1}}{k_{d}},$$ (4.8)

controlling the time extension of the pre-bouncung high curvature regime and of the two post-bouncing, non-standard phases, we find (using eq. (4.5))

$$\!\!\!\!\!\!\!\!z_{\sigma}=\frac{k_{1}}{k_{\sigma}}=\frac{H_{1}a_{1}}{H_{\sigma}a_{\sigma}}=\left(\frac{H_{1}}{H_{\sigma}}\right)^{\frac{1}{2}}\simeq\left(\frac{H_{1}}{M_{\rm Pl}}\right)^{\frac{1}{2}}\left(\frac{m}{M_{\rm Pl}}\right)^{-\frac{1}{2}}\left(\frac{\sigma_{i}}{M_{\rm Pl}}\right)^{-2},$$

$$\!\!\!\!\!\!\!\!\frac{z_{d}}{z_{\sigma}}=\frac{k_{\sigma}}{k_{d}}=\frac{H_{\sigma}a_{\sigma}}{H_{d}a_{d}}=\left(\frac{H_{\sigma}}{H_{d}}\right)^{\frac{1}{3}}\simeq\left(\frac{m}{M_{\rm Pl}}\right)^{-\frac{2}{3}}\left(\frac{\sigma_{i}}{M_{\rm Pl}}\right)^{\frac{4}{3}},$$

$$\!\!\!\!\!\!\!\!\frac{z_{s}}{z_{d}}=\frac{k_{d}}{k_{s}}=\frac{k_{d}}{k_{1}}\frac{k_{1}}{k_{s}}=z_{s}\frac{H_{d}a_{d}}{H_{1}a_{1}}=z_{s}\frac{H_{d}}{H_{1}}\frac{a_{d}}{a_{\sigma}}\frac{a_{\sigma}}{a_{1}}\simeq z_{s}\left(\frac{H_{1}}{M_{\rm Pl}}\right)^{-\frac{1}{2}}\left(\frac{m}{M_{\rm Pl}}\right)^{\frac{7}{6}}\left(\frac{\sigma_{i}}{M_{\rm Pl}}\right)^{\frac{2}{3}}.$$ (4.9)

By inverting the above relations we can also express $\sigma_{i}$ and $m$ in terms of $H_{1}$ and of the parameters $z_{s}$, $z_{\sigma}$ and $z_{d}$ as follows:

$$\frac{\sigma_{i}}{M_{\rm Pl}}\simeq\left(\frac{H_{1}}{M_{\rm Pl}}\right)^{\frac{1}{6}}z_{d}^{\frac{1}{4}}z_{\sigma}^{-\frac{7}{12}}\,,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \qquad\frac{m}{M_{\rm Pl}}\simeq\left(\frac{H_{1}}{M_{\rm Pl}}\right)^{\frac{1}{3}}z_{d}^{-1}z_{\sigma}^{\frac{1}{3}}.$$ (4.10)

Let us now give an example of full computation of the spectral distribution (1.3) for the highest frequency branch of the spectrum ($k_{\sigma}<k<k_{1}$) and for a generic epoch $\tau$ (late enough, however, to have all such modes inside the horizon, $k\tau\gg 1$). To this purpose, let us express the energy density of the perturbations in terms of their (time-dependent) proper frequency $\omega$ scaling in time like the inverse scale factor, $\omega(\tau)=k/a(\tau)$. From eqs. (1.3) and (4.1), we obtain

$$\Omega_{\textsc{gw}}(k,\tau)=\frac{\omega^{4}}{\pi^{2}\rho_{\rm crit}(\tau)}\,\langle n_{\omega}(\tau)\rangle\,,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \qquad k_{\sigma}<k<k_{1}.$$ (4.11)

The modes we are considering are amplified by the pump field (4.3) of the high-curvature string phase, and the canonical equations (4.6), (4.7) then give [95]

$$\langle n_{\omega}(\tau)\rangle\simeq\left(\frac{\omega}{\omega_{1}}\right)^{-1-|3-2\beta|}=\left(\frac{k}{k_{1}}\right)^{-1-|3-2\beta|}.$$ (4.12)

Also, it is convenient to express the critical density $\rho_{\rm crit}$ in terms of the critical fraction of radiation energy density, $\Omega_{\rm r}=\rho_{r}/\rho_{\rm crit}$, so that, referring to radiation produced at the axion-decay scale, we have

$$\rho_{\rm crit}(\tau)=\frac{\rho_{r}(\tau)}{\Omega_{\rm r}(\tau)}=\frac{3M_{\rm Pl}^{2}H_{d}}{\Omega_{\rm r}(\tau)}\left(\frac{a_{d}}{a}\right)^{4}.$$ (4.13)

Finally, for each mode of proper frequency $\omega(\tau)$ we can express its $\omega^{4}$ contribution to eq. (4.11) as

$$\omega^{4}=\left(\frac{k}{a}\right)^{4}=\left(\frac{k}{k_{1}}\right)^{4}\left(\frac{H_{1}a_{1}}{a}\right)^{4}.$$ (4.14)

Combining all these results and using eq. (4.5) we then find, for the considered band of frequency,

$$\Omega_{\textsc{gw}}(k,\tau)=\frac{\Omega_{\rm r}(\tau)}{\Omega_{{\rm r}0}}\,\Omega_{\rm PBB}\left(\frac{k}{k_{1}}\right)^{3-|3-2\beta|},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ k_{\sigma}<k<k_{1},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \tau>\tau_{d},$$ (4.15)

where $\Omega_{{\rm r}0}\equiv\Omega_{\rm r}(\tau_{0})\approx 4.15\times 10^{-5}h^{-2}$ is the present critical fraction of radiation energy density (including neutrinos) and we have defined the constant (parameter dependent) dimensionless amplitude

$$\Omega_{\rm PBB}\coloneqq\Omega_{{\rm r}0}\left(\frac{H_{1}}{M_{\rm Pl}}\right)^{2}\left(\frac{mM_{\rm Pl}}{\sigma_{i}^{2}}\right)^{4/3}=\Omega_{{\rm r}0}\left(\frac{H_{1}}{M_{\rm Pl}}\right)^{2}\left(\frac{z_{\sigma}}{z_{d}}\right)^{2}.$$ (4.16)

For simplicity we have absorbed all numerical factors of order of unity into the unknown scale $H_{1}$. Obviously, the result (4.15) is also valid if applied in particular to the present epoch $\tau=\tau_{0}$, with $\Omega_{\rm r}(\tau)=\Omega_{{\rm r}0}$.

By following the same procedure for the other (lower frequency) branches of the spectrum, and turning to the more conventional frequency variable $f=k/(2\pi)$, we obtain that the full GWB can be written synthetically as

$$\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\Omega_{\textsc{gw}}(f)=\left\{\begin{matrix}\Omega_{\rm PBB}\left(\dfrac{f}{f_{1}}\right)^{3-|3-2\beta|},\qquad\qquad\qquad\qquad\qquad f_{\sigma}\leavevmode\nobreak\ \raise 1.72218pt\hbox{$<$}\kern-8.00003pt\lower 2.6694pt\hbox{$\sim$}\leavevmode\nobreak\ f\leavevmode\nobreak\ \raise 1.72218pt\hbox{$<$}\kern-8.00003pt\lower 2.6694pt\hbox{$\sim$}\leavevmode\nobreak\ f_{1}\\ \\ \Omega_{\textsc{gw}}(f_{1})\left(\dfrac{f_{\sigma}}{f_{1}}\right)^{3-|3-2\beta|}\left(\dfrac{f}{f_{\sigma}}\right)^{1-|3-2\beta|},\qquad\quad f_{d}\leavevmode\nobreak\ \raise 1.72218pt\hbox{$<$}\kern-8.00003pt\lower 2.6694pt\hbox{$\sim$}\leavevmode\nobreak\ f\leavevmode\nobreak\ \raise 1.72218pt\hbox{$<$}\kern-8.00003pt\lower 2.6694pt\hbox{$\sim$}\leavevmode\nobreak\ f_{\sigma}\\ \\ \Omega_{\textsc{gw}}(f_{\sigma})\left(\dfrac{f_{d}}{f_{\sigma}}\right)^{1-|3-2\beta|}\left(\dfrac{f}{f_{d}}\right)^{3-|3-2\beta|},\qquad\,\,f_{s}\leavevmode\nobreak\ \raise 1.72218pt\hbox{$<$}\kern-8.00003pt\lower 2.6694pt\hbox{$\sim$}\leavevmode\nobreak\ f\leavevmode\nobreak\ \raise 1.72218pt\hbox{$<$}\kern-8.00003pt\lower 2.6694pt\hbox{$\sim$}\leavevmode\nobreak\ f_{d}\\ \\ \Omega_{\textsc{gw}}(f_{d})\left(\dfrac{f_{s}}{f_{d}}\right)^{3-|3-2\beta|}\left(\dfrac{f}{f_{s}}\right)^{3},\qquad\qquad\qquad\quad\,\,\,\,\,f\leavevmode\nobreak\ \raise 1.72218pt\hbox{$<$}\kern-8.00003pt\lower 2.6694pt\hbox{$\sim$}\leavevmode\nobreak\ f_{s}\end{matrix}\right.$$ (4.17)

where the three ratios of frequency scales $f_{\sigma}/f_{1}$, $f_{d}/f_{\sigma}$, $f_{s}/f_{d}$ can be expressed in terms of the parameters $z_{s}$, $z_{\sigma}$, $z_{d}$ (or $m$, $\sigma_{i}$, $H_{1}$) according to eq. (4.9). It is important to note that, for $f>f_{1}$, the spectrum is exponentially suppressed as $\Omega_{\textsc{gw}}(f)=\Omega_{\rm PBB}\exp{[-(f-f_{1})/f_{1}}]$, so that a smooth interpolation of all the branches can be given by the following expression proven in Appendix A:

$$\Omega_{\textsc{gw}}(f) = \,\Omega_{\rm PBB}\,f^{3}\left(f^{2}+f_{s}^{2}\right)^{-\frac{\left|3-2\beta\right|}{2}}\left(f^{2}+f_{d}^{2}\right)^{-1}\left(f^{2}+f_{\sigma}^{2}\right)\left(f^{2}+f_{1}^{2}\right)^{\frac{\left|3-2\beta\right|-3}{2}}$$ (4.18)

$$\times\exp{\left({-\frac{f}{f_{1}}+\arctan{\frac{f}{f_{1}}}}\right)}\,.$$

The smoothing of the piecewise profile (4.17) does not change the underlying scenario because the transition epochs from one phase to another are of very short, negligible duration compared with the time extension of such phases.

There are now two important points to be stressed. The first is that the overall GWB, and in particular the peak amplitude, is controlled not only by the bouncing curvature scale $H_{1}$ but also by the parameters $m$ and $\sigma_{i}$ and, thus, by the details of the post-bounce evolution. This implies, in particular, that the amplitude may result strongly suppressed with respect to the natural value fixed by the fundamental string mass scale $M_{\rm s}$, even for the highest frequency modes crossing the horizon at such scale.

The second point concerns the number of parameters controlling the shape of the spectrum. There are four time (or curvature) scales and one dimensionless number, the power $\beta$ of the pump field in the string phase; see (4.3). However, these five parameters must satisfy an important phenomenological condition. The PBB scenario we are considering, in fact, besides producing relic GW radiation must also produce a suitable large-scale background of scalar curvature perturbations with a nearly flat spectrum, in order to be compatible with CMB observations.

This is known to be possible via the curvaton mechanism [104, 105, 106] triggered by the Kalb–Ramond axion [96, 97], but this imposes constraints on the previous spectral parameters [58]. In particular, the primordial scalar amplitude $\mathcal{P}_{\rm s}(k_{*})$ and spectral index $n_{\rm s}$ must be in agreement with the observational results reported below (1.2). By imposing such conditions on the scalar perturbations produced by the PBB model we are considering, whose spectrum is controlled by the same set of parameters as the GWB of (4.17) or its smooth version (4.18), we can then eliminate one of the previous five parameters and fix, for instance, the transition scale $H_{1}$ as a function of $z_{s}$, $\beta$, $m$, $\sigma_{i}$ and of the two observables $\mathcal{P}_{\rm s}(k_{*})$ and $n_{\rm s}$. This is done in Appendix B.1, where we conclude that the parameter space of the model is

$$\{\beta,z_{s},z_{d},z_{\sigma}\}\,.$$ (4.19)

Without assuming any prior on these parameter, one can plot (4.18) and find some general trends. For instance, the larger $\beta$, the smaller the amplitude, so that, in practice, only values of $\beta$ near zero generate a detectable signal. Also, as the parameter $z_{d}$ increases, the spectral shape is squeezed and the peak becomes narrower but neither its frequency nor its amplitude vary appreciably. On the other hand, the parameter $z_{\sigma}$ changes the shape but not much the frequency peak or the amplitude, while the parameter $z_{s}$ which affects both the range and the amplitude of the intermediate plateau or slope.

Figure 3: The shaded area shows the allowed region in the $(\log_{10}z_{s},\log_{10}z_{d})$ plane for the limiting case $\beta=0$ and $\sigma_{i}=M_{\rm Pl}$, which maximizes the peak amplitude of the primordial GWB. The maximal allowed intensity is reached for parameter values lying on the dashed straight line marking the right border of the region, corresponding to $\log_{10}(H_{1}/M_{\rm Pl})\approx-3.29$, as explained in Appendix B.3.

However, the four parameters (4.19) do obey a non-trivial set of theoretical priors determined in Appendix B.2. Within this space, in Appendix B.3 we circumscribe the region in the parameter space where the peak amplitude of the GWB is maximized. This happens for $\beta=0$, $\sigma_{i}=M_{\rm Pl}$, $z_{\sigma}$ determined by eq. (B.13) and the values of $z_{s}$ and $z_{d}$ shown in Fig. 3.

Any given set of parameters $z_{s},z_{d},z_{\sigma}$ satisfying eq. (B.13) together with all other constraints, and implementing the additional limiting condition $\log_{10}(H_{1}/M_{\rm Pl})\approx-3.29$, produces a GWB with a peak of maximum intensity $\Omega_{\textsc{gw}}^{\rm max}\sim 10^{-11}\!-\!10^{-10}$ (eq. (B.15)). For phenomenological reasons, however, we are interested not only in the maximal amplitude but also in the maximal extension in frequency (in particular, towards small frequencies) of the allowed spectral region. This last property can be easily obtained by choosing, among all possible combinations of parameters producing the maximal amplitude, the combination selecting the maximal allowed value of $z_{s}$ (i.e., $z_{s}\approx 10^{22.3}$), together with the corresponding minimal value of $z_{d}$ (i.e., $z_{d}\approx 10^{2.19}$), together with the associated minimal value of $z_{\sigma}$ (i.e., $z_{\sigma}\approx 1$, according to eq. (B.13)).

Figure 4: The shaded area is the maximal allowed region for the GWB (4.18) extended in the low-frequency range, with all its parameters satisfying the self-consistency constraints (B.7)–(B.11). We have also plotted four possible spectra corresponding to different sets of parameters giving rise to different kinematic details of the phases of earlier, non-standard cosmic evolution, preceding and following the curvature bounce.

The GWB (4.18) with the above numerical values of the parameters (satisfying $\sigma_{i}=M_{\rm Pl}$) and with $\beta=0$ provides the border of the shaded area in Fig. 4 within which the primordial GWB of pre-big-bang cosmology with maximal peak amplitude falls. Such region respects known constraints obtained from the observations of millisecond pulsars [107], which imply $\Omega_{\textsc{gw}}\leavevmode\nobreak\ \raise 1.72218pt\hbox{$<$}\kern-8.00003pt\lower 2.6694pt\hbox{$\sim$}\leavevmode\nobreak\ 10^{-8}$ at a frequency scale $f\sim 10^{-8}$ Hz. It is also consistent with large-scale CMB constraints on the tensor spectrum. In the frequency range of PTA, the maximal GWB amplitude is smaller than the signal detected [3, 4, 5, 6, 7] and is therefore consistent with those observations.

For an illustrative purpose, in Fig. 4 we have also plotted a few examples of GWBs produced by different sets of parameters satisfying all required constraints. The plotted spectra maximize neither the peak amplitude nor the extension in frequency but they are well inside the border of the allowed region and they can still produce a detectable signal within the LISA and the ET frequency range.

5 Discussion

In this paper, we have continued the study of [29] on the tensor-produced primordial GWB in bouncing scenarios embedded in, or inspired by, semi-classical models, quantum gravity and string theory. A grand summary of our present results is the following.

• Cosmological models of string and quantum gravity such as inflation in flux compactification (see [109, 108] for reviews), nonlocal Starobinsky inflation [37, 38], Wheeler–DeWitt quantum cosmology [110, 111, 112, 113, 114, 115] and various incarnations of loop quantum cosmology [116, 117, 118, 119, 120] generate a red-tilted primordial tensor spectrum ($n_{\rm t}<0$) and the resulting GWB is not enhanced with respect to the standard inflationary one. It is thus unobservable by present and planned interferometers [29].

• Scenarios such as Brandenberger–Ho non-commutative inflation [39, 40], multi-fractional inflation [48] and the conformal early universe in nonlocal quantum gravity [32] also generate a blue-tilted tensor spectrum but the primordial GWB amplitude can only reach the DECIGO window [29, 32]. These models may allow for a strong signal at CMB scales already within observability reach but this is a fixed prediction only for the conformal early-universe scenario [32].

• The GWB of the new ekpyrotic universe [41, 42, 43] has a very low amplitude at CMB scales and it can reach the DECIGO window only if the tensor spectral index is running [29]. In some extreme but not very realistic cases, it can even reach the ET window. The ekpyrotic scenario with Galileons and a $U(1)$ gauge field studied here [49, 50, 51, 52, 53] predicts a sourced tensor amplitude which is observable only on CMB frequency range, while the vacuum GWB is unobservable in all other upcoming GWB observations.

• The GWB of string-gas cosmology is highly sensitive to the underlying assumptions. If the signal is produced below the Hagedorn temperature [44, 45, 46, 47], then it can reach at most the DECIGO sensitivity window assuming small or no running of the spectral indices [29]. If, as explored here, it is produced above the Hagedorn temperature [54], then it becomes too high and is observationally excluded. A possible way out of this would be to revise the assumptions underlying the model, in particular, to be already in Einstein’s gravity regime (3.9) when the primordial spectra were produced. Relaxing this condition and allowing generation of the spectra before recovering Einstein’s gravity might yield different results.

• The only model so far striking a balance between observability in the high-frequency range and observational consistency is the pre-big-bang scenario [55, 56, 57, 58]. We found that its primordial GWB has a convex shape with an intermediate flat plateau. The parameter space of the model is tightly constrained theoretically but it still gives enough phenomenological freedom, to the point that this GWB can comfortably fall within the sensitivity window of both LISA and ET. When the parameter $\beta$ is close to zero, in this range of frequencies the GWB of this model is a single power law

  $$\Omega_{\text{GW}}(f)\sim f^{n_{\textsc{gw}}}\,,\qquad n_{\textsc{gw}}=2+n_{\rm t}\simeq 3-|3-2\beta|\approx 0\,.$$ (5.1)

  When $z_{d}$ grows, the signal becomes a broken power law,

  $$\Omega_{\textsc{gw}}(f)\sim\left\{\begin{matrix}\,\,f^{n_{\textsc{gw},1}},\qquad\quad f\ll\bar{f}\\ f^{n_{\textsc{gw},2}},\qquad\quad f\gg\bar{f}\end{matrix}\right.,$$ (5.2)

  where the slopes depend on where the transition scale $\bar{f}$ lies in the model’s scales hierarchy and on whether the intermediate plateau is extended enough:

  $$\bar{f}=f_{d}\,:\hphantom{\simeq f_{d}}\quad n_{\textsc{gw},1}\simeq 3-|3-2\beta|\approx 0\,,\quad\,\,n_{\textsc{gw},2}\simeq 1-|3-2\beta|\approx-2\,,$$

  $$\bar{f}=f_{s}\,:\hphantom{\simeq f_{d}}\quad n_{\textsc{gw},1}\approx 3\,,\qquad\qquad\qquad\qquad n_{\textsc{gw},2}\simeq 3-|3-2\beta|\approx 0\,,$$ (5.3)

  $$\bar{f}=f_{s}\simeq f_{d}\,:\quad\!\!n_{\textsc{gw},1}\approx 3\,,\qquad\qquad\qquad\qquad n_{\textsc{gw},2}\simeq 1-|3-2\beta|\approx-2\,.$$

Overall, among all these early-universe scenarios only the pre-big-bang model is of interest for LISA and ET. Both the single- and the broken-power-law shapes (5.1) and (5.3) fit the two simplest templates commonly used for inflationary models [121] and they could be submitted to the same type of analysis performed in [122] as well as to a discussion on degeneracies and model selection when compared with alternative candidates with similar spectral shapes.

In general, the higher the signal-to-noise ratio, the smaller the error on the parameters of the spectral shape, as has been verified for a variety of cases [33, 121, 122]. Thus, for the PBB model one would be able to identify the corner of the parameter space such that the tilt or tilts of the GWB (and, thus, $\beta$) would be determined with an accuracy of, say, 1% by LISA and ET in their respective observability window.

In the case of the broken power law, the sharper the transition from one slope to the other with respect to the sensitivity curve the more difficult the determination of the parameters, due to the fact that the binning of the data introduces a certain level of coarse graining of short-range features [123]. We do not expect to have this problem for the PBB model, since the aforementioned transition is gentle enough (Fig. 4). These and other aspects of the PBB model will deserve to be explored in the future.

Appendix A Smooth interpolation of a broken power law

Suppose we have a piecewise continuous function $f(x):\mathbb{R}^{+}\rightarrow\mathbb{R}^{+}$ such that, in each interval, its behaviour is given by a power law $f(x)\sim x^{m}$. Then, we can define a sequence of exponents $m_{1},m_{2},\dots$ and interval extrema $x_{1,2},x_{2,3},\dots$ such that $x_{i,i+1}<x_{i+1,i+2}$. The derivative of the logarithm of the function with respect to the logarithm of its argument is simply given by

$$\frac{d\ln f(x)}{d\ln x}=m_{1}+(m_{2}-m_{1})\Theta(\ln x-\ln x_{1,2})+(m_{3}-m_{2})\Theta(\ln x-\ln x_{2,3})+\dots,$$ (A.1)

where $\Theta$ is the Heaviside step function: $\Theta(y)=1$ for $y\geqslant 0$ and $\Theta(y)=0$ otherwise. Let us approximate the Heaviside function with a logistic function:

$$\Theta(y-y_{0})\simeq\frac{1}{1+e^{-l(y-y_{0})}},$$ (A.2)

where $l$ is a positive constant and the greater $l$ the better the approximation. Substituting the latter into (A.1) and defining a different $l$ for each $\Theta$, we have in compact form

$$\frac{d\ln f(x)}{d\ln x}=m_{1}+\sum_{i=1}^{N-1}\frac{m_{i+1}-m_{i}}{1+(\frac{x}{x_{i,i+1}})^{-l_{i+1}}}\,.$$ (A.3)

We can moreover rewrite the last equality with respect to the $x$-derivative as

$$\frac{d\ln f(x)}{dx}=\frac{m_{1}}{x}+\sum_{i=1}^{N-1}\frac{m_{i+1}-m_{i}}{x}\frac{1}{1+(\frac{x}{x_{i,i+1}})^{-l_{i+1}}}\,.$$ (A.4)

Integrating the last expression, we obtain

$$\ln f(x)=\ln A+\ln(x^{m_{1}})+\sum_{i=1}^{N-1}\ln(x^{l_{i+1}}+x_{i,i+1}^{l_{i+1}})^{(m_{i+1}-m_{i})/l_{i+1}},$$

where $A$ is a constant, and finally

$$f(x) = A\,x^{m_{1}}\prod_{i=1}^{N-1}\left(x^{l_{i+1}}+x_{i,i+1}^{l_{i+1}}\right)^{(m_{i+1}-m_{i})/l_{i+1}}$$ (A.5)

$$= A\,x^{m_{1}}\left(x^{l_{2}}+x_{1,2}^{l_{2}}\right)^{(m_{2}-m_{1})/l_{2}}\left(x^{l_{3}}+x_{2,3}^{l_{3}}\right)^{(m_{3}-m_{2})/l_{3}}\times\dots.\,$$

If we wish to introduce an exponential cut-off for $x>x_{N,N+1}$ such that $f(x)\sim\exp[-(x-x_{N,N+1})/x_{N,N+1}]$, then we have to add to (A.1) an additional contribution given by $(-x/x_{N,N+1}-m_{N+1})\Theta(\ln x-\ln x_{N,N+1})$. Let us denote $x_{N,N+1}\equiv x_{M}$ and $l_{N+1}\equiv l$. Following the same steps as before, one finds upon integration and exponentiation that

$$f(x) = A\,x^{m_{1}}\prod_{i=1}^{N-1}\left(x^{l_{i+1}}+x_{i,i+1}^{l_{i+1}}\right)^{(m_{i+1}-m_{i})/l_{i+1}}\left(x^{l}+x_{M}^{l}\right)^{-m_{N+1}/l}$$ (A.6)

$$\times\exp\left[-\frac{(x/x_{M})^{1+l}\,_{2}F_{1}(1,\,1+1/l,\,2+1/l,\,-(x/x_{M})^{l})}{1+l}\right],$$

where ${}_{2}F_{1}(a,b,c,z)$ is the Gauss hypergeometric function. Translating into the general language we introduced before for the power spectrum of (4.17), we have $m_{1}=3$, $m_{2}=3-\left|3-2\beta\right|$, $m_{3}=1-\left|3-2\beta\right|$ , $m_{4}=3-\left|3-2\beta\right|$, $x_{1,2}=f_{s}$, $x_{2,3}=f_{d}$ and $x_{3,4}=f_{\sigma}$ and $x_{M}=f_{1}$, so that (setting $l_{i}=l$ for all $i$)

$$\Omega_{\text{GW}}^{\text{smooth}}(f) = A\,f^{3}\left(f^{l}+f_{s}^{l}\right)^{-\frac{\left|3-2\beta\right|}{l}}\left(f^{l}+f_{d}^{l}\right)^{-\frac{2}{l}}\left(f^{l}+f_{\sigma}^{l}\right)^{\frac{2}{l}}$$ (A.7)

$$\times\left(f^{l}+f_{1}^{l}\right)^{\frac{\left|3-2\beta\right|-3}{l}}\mathcal{F}(f,f_{1},l)\,,$$

where we defined the function $\mathcal{F}(f,f_{1},l)$ as the exponential cutoff appearing in the last line of (A.6). Setting $l=2$, we obtain (4.18):

$$\Omega_{\textsc{gw}}^{\text{smooth}}(f) = A\,f^{3}\left(f^{2}+f_{s}^{2}\right)^{-\frac{\left|3-2\beta\right|}{2}}\left(f^{2}+f_{d}^{2}\right)^{-1}\left(f^{2}+f_{\sigma}^{2}\right)\left(f^{2}+f_{1}^{2}\right)^{\frac{\left|3-2\beta\right|-3}{2}}$$ (A.8)

$$\times\exp{\left({-\frac{f}{f_{1}}+\arctan{\frac{f}{f_{1}}}}\right)}.$$

The integration constant $A$ is fixed to match the asymptotic behaviour in the limit $f\rightarrow 0$ with the explicit formula given in (4.17),

$$\Omega_{\textsc{gw}}(f)\stackrel{{\scriptstyle f\rightarrow 0}}{{\simeq}}{\Omega}_{\text{PBB}}\,(f_{\sigma})^{2}(f_{d})^{-2}(f_{s})^{-\left|3-2\beta\right|}(f_{1})^{\left|3-2\beta\right|-3}f^{3},$$ (A.9)

while for the smooth interpolation (A.8)

$$\Omega_{\textsc{gw}}^{\text{smooth}}(f)\stackrel{{\scriptstyle f\rightarrow 0}}{{\simeq}}A\,(f_{\sigma})^{2}(f_{d})^{-2}(f_{s})^{-\left|3-2\beta\right|}(f_{1})^{\left|3-2\beta\right|-3}f^{3},$$ (A.10)

which yields $A={\Omega}_{\text{PBB}}$.

Appendix B Parameter space of the PBB model

In this appendix, we show that the parameter space of the PBB model is $\{\beta,z_{s},z_{d},z_{\sigma}\}$ (section B.1) and we discuss the theoretical priors on these parameters (section B.2). We also determine the region in parameter space for which the promordial GWB amplitude is maximized (section B.3). In view of a numerical implementation of these conditions and to emphasize orders of magnitude, it may be useful to work with base-10 logarithmic expressions.

B.1 Reducing the number of parameters

First, we show that the transition scale $H_{1}$ is not independent and can be fixed by the other parameters,

$$H_{1}=H_{1}(\beta,z_{s},z_{d},z_{\sigma})\,.$$ (B.1)

By using previous results [58, 98], obtained under the natural assumption that the pivot scale belongs to the low-frequency band of the scalar spectrum (i.e., $k_{*}<k_{s}$), the condition on $H_{1}$ following from the normalization of the scalar spectrum can be written as

$$\left(\frac{H_{1}}{M_{\rm Pl}}\right)^{\frac{5-n_{\rm s}}{2}} = \frac{2\pi^{2}}{\mathcal{T}^{2}(\sigma_{i})}\mathcal{P}_{\rm s}(k_{*})\,z_{s}^{1-n_{\rm s}-2\beta}\left[\left(\frac{H_{*}}{M_{\rm Pl}}\right)^{-\frac{1}{2}}\left(\frac{mM_{\rm Pl}}{\sigma_{i}^{2}}\right)^{\frac{1}{3}}\right]^{n_{\rm s}-1}$$ (B.2)

$$= \frac{2\pi^{2}}{\mathcal{T}^{2}(\sigma_{i})}\mathcal{P}_{\rm s}(k_{*})\,z_{s}^{1-n_{\rm s}-2\beta}\left(\frac{H_{*}}{M_{\rm Pl}}\frac{z_{d}}{z_{\sigma}}\right)^{-\frac{n_{\rm s}-1}{2}},$$

where we have recast $m$ and $\sigma_{i}$ according to (4.10), $H_{*}$ is the curvature scale at the epoch in which the pivot mode $k_{*}$ re-enters the horizon and $\mathcal{T}(\sigma_{i})$ is the transfer function connecting the amplitude of the primordial axion fluctuations to the final amplitude of the scalar curvature modes of metric perturbations. A numerical integration of the scalar perturbation equations gives the simple result [97]

$$\mathcal{T}(\sigma_{i})\simeq 0.13\,\frac{\sigma_{i}}{M_{\rm Pl}}+0.25\,\frac{M_{\rm Pl}}{\sigma_{i}}-0.01\,,$$ (B.3)

where $\sigma_{i}$ can be expressed in terms of the $z$ parameters as in (4.10). To obtain $H_{*}$, we can conveniently refer to the equilibrium scale by noting that $k_{*}\simeq 5k_{\rm eq}$. This implies $H_{*}^{1/2}\simeq 5H_{\rm eq}^{1/2}$. On the other hand, it is known that the Hubble parameter at radiation-matter equality is given by $H_{\rm eq}\simeq 1.6\times 10^{5}H_{0}\approx 9.5\times 10^{-56}M_{\rm Pl}$. We thus obtain

$$\left(\frac{H_{*}}{M_{\rm Pl}}\right)^{1/2}\approx 1.5\times 10^{-27}\,.$$ (B.4)

Now we can then express the normalization (B.2) in terms of the four parameters $\{\beta,z_{s},z_{d},z_{\sigma}\}$ (and of known experimental numbers) as follows:

$$\log_{10}\left(\frac{H_{1}}{M_{\rm Pl}}\right) = \frac{2}{5-n_{\rm s}}\left\{\log_{10}\left[\frac{4.2\pi^{2}}{\mathcal{T}^{2}(\sigma_{i})}\right]-9+(1-n_{\rm s})(\log_{10}1.5-27)\right.$$ (B.5)

$$+\left.(1-n_{\rm s}-2\beta)\log_{10}z_{s}+\frac{n_{\rm s}-1}{2}\left(\log_{10}z_{\sigma}-\log_{10}z_{d}\right)\right\},$$

where we have used $\mathcal{P}_{\rm s}(k_{*})=2.1\times 10^{-9}$. It should be noted that $\mathcal{T}^{2}(\sigma_{i})$ also contains $H_{1}$ through (4.10) but the solution for $H_{1}$ can always be numerically obtained, in general, for any given set of values of the four independent parameters.

Finally, the other important quantity appearing in the GWB (4.17) is today’s value of the highest amplified frequency mode $f_{1}$, which is given by

$$f_{1} = \frac{\omega_{1}(\tau_{0})}{2\pi}=\frac{H_{1}a_{1}}{2\pi a_{0}}=\frac{H_{1}}{2\pi}\frac{a_{1}}{a_{\sigma}}\frac{a_{\sigma}}{a_{d}}\frac{a_{d}}{a_{\rm eq}}\frac{a_{\rm eq}}{a_{0}}\simeq\frac{H_{1}}{2\pi}\left(\frac{H_{\sigma}}{H_{1}}\right)^{\frac{1}{2}}\left(\frac{H_{d}}{H_{\sigma}}\right)^{\frac{2}{3}}\left(\frac{H_{\rm eq}}{H_{d}}\right)^{\frac{1}{2}}\left(\frac{H_{0}}{H_{\rm eq}}\right)^{\frac{2}{3}}=$$ (B.6)

$$= \frac{H_{1}^{\frac{1}{2}}}{2\pi}\left(\frac{z_{\sigma}}{z_{d}}\right)^{\frac{1}{2}}\frac{H_{0}^{\frac{2}{3}}}{H_{\rm eq}^{\frac{1}{6}}}\simeq\frac{3.9\times 10^{11}}{2\pi}\left(\frac{H_{1}}{M_{\rm Pl}}\right)^{\frac{1}{2}}\left(\frac{z_{\sigma}}{z_{d}}\right)^{\frac{1}{2}}\,{\rm Hz}\,.$$

B.2 Theoretical priors

Having thus determined that the parameter space is $\{\beta,z_{s},z_{d},z_{\sigma}\}$, let us now turn to the priors we can impose on it theoretically.

A first condition concerns the parameter $\beta$ controlling the power-law behaviour of the primordial GW spectrum at high frequencies, which is constrained to be in the range

$$0\leqslant\beta<3\,.$$ (B.7)

The lower limit is due to the assumption of growing string coupling (needed to implement a smooth bouncing transition [103, 60, 61]), while the upper limit has to be imposed to avoid background instabilities [124].

We have then a number of constraints following from the (already mentioned) hierarchy of the transition frequency scales, which must satisfy the conditions $f_{1}\gtrsim f_{\sigma}>f_{d}>f_{s}$. They imply

$$1\lesssim z_{\sigma}<z_{d}<z_{s}\,.$$ (B.8)

In addition, for an efficient implementation of the curvaton mechanism based on the oscillations of the Kalb–Ramond axion, the axion background field must be oscillating when it becomes dominant (at the curvature scale $H_{\sigma}$). From the axion dynamical equations, one finds [96, 97] that the oscillating regime starts at the scale $H_{m}\simeq m$. This leads to the condition $H_{m}\geqslant H_{\sigma}$, which implies $\sigma_{i}\leqslant M_{\rm Pl}$ and which, by using Eq. (4.10), can be written in logarithmic form as

$$\log_{10}\left(\frac{H_{1}}{M_{\rm Pl}}\right)+\frac{3}{2}\log_{10}z_{d}-\frac{7}{2}\log_{10}z_{\sigma}\leqslant 0\,.$$ (B.9)

Also, to be consistent with the established results of the post-inflationary scenario, we may expect that the reheating produced by the axion decay at the scale $H_{d}$, and marking the beginning of the standard cosmological evolution, occurs before the BBN scale, $H_{\textsc{bbn}}\simeq(1\,{\rm MeV})^{2}/M_{\rm Pl}$. This implies $H_{d}>H_{\textsc{bbn}}$ from which, using eqs. (4.5) and (4.10), we have the constraint

$$\log_{10}\left(\frac{H_{1}}{M_{\rm Pl}}\right)-3\log_{10}z_{d}+\log_{10}z_{\sigma}>-42-\log_{10}4\,,$$ (B.10)

where we have used $M_{\rm Pl}\approx 2\times 10^{18}$ GeV.

Finally, the conditions concerning the scalar perturbation spectrum must be imposed not only at the pivot frequency scale $k_{*}$ but also, in principle, to all frequency scales included into the multipole expansion of the CMB anisotropy, and constrained by observational data. This means, in other words, that also the highest frequency modes $k_{\textsc{lss}}$ presently constrained by large scale structure (LSS) observations must be below the lowest frequency branch of the axion perturbation spectrum [58, 98], and this implies $k_{\textsc{lss}}<k_{s}$, where $k_{\textsc{lss}}\sim 3\,{\rm Mpc}^{-1}\approx 60\,k_{*}$, namely $H_{\textsc{lss}}^{1/2}\simeq 60\,H_{*}^{1/2}$. The condition $k_{\textsc{lss}}/k_{s}=(k_{\textsc{lss}}/k_{1})z_{s}<1$ then leads to yet another constraint that can be written as follows. Since

$$\frac{k_{1}}{k_{\textsc{lss}}} = \frac{H_{1}a_{1}}{H_{\textsc{lss}}a_{\textsc{lss}}}=\frac{H_{1}}{H_{\textsc{lss}}}\frac{a_{1}}{a_{\sigma}}\frac{a_{\sigma}}{a_{d}}\frac{a_{d}}{a_{\textsc{lss}}}\simeq\frac{H_{1}}{H_{\textsc{lss}}}\left(\frac{H_{\sigma}}{H_{1}}\right)^{\frac{1}{2}}\left(\frac{H_{d}}{H_{\sigma}}\right)^{\frac{2}{3}}\left(\frac{H_{\textsc{lss}}}{H_{d}}\right)^{\frac{1}{2}}$$

$$= \left(\frac{H_{1}}{M_{\rm Pl}}\right)^{\frac{1}{2}}\left(\frac{H_{*}}{M_{\rm Pl}}\right)^{-\frac{1}{2}}\left(\frac{H_{*}}{H_{\textsc{lss}}}\right)^{\frac{1}{2}}\left(\frac{H_{\sigma}}{H_{d}}\right)^{-\frac{1}{6}},$$

from eqs. (4.5) and (4.10) we get $H_{\sigma}/H_{d}=(z_{d}/z_{\sigma})^{3}$ and from (B.4) we obtain

$$\log_{10}z_{s}<26-\log_{10}9+\frac{1}{2}\log_{10}\left(\frac{H_{1}}{M_{\rm Pl}}\right)+\frac{1}{2}\left(\log_{10}z_{\sigma}-\log_{10}z_{d}\right)\,.$$ (B.11)

B.3 Maximizing the signal

Given the condition (B.5) on $H_{1}/M_{\rm Pl}$, the amplitude and the frequency distribution of the GWB (4.17) or (4.18) are fully determined by $\beta,z_{s},z_{d},z_{\sigma}$. These four parameters are not completely free, as they must satisfy a non-trivial set of self-consistency conditions (Appendix B). Taking into account these constraints on the parameters, we can determine the maximal allowed region for the PBB signal in the spectral plane $(\Omega_{\textsc{gw}},f)$. Let us first notice that, thanks to the condition (B.7), the GWB (4.18) may be possibly decreasing only in the frequency branch $f_{d}\leqslant f\leqslant f_{\sigma}$. The peak of the spectrum may thus be located either at $f_{1}$ or at $f_{d}$, with corresponding amplitudes

$$\Omega_{\textsc{gw}}(f_{1})=\Omega_{{\rm r}0}\left(\frac{H_{1}}{M_{\rm Pl}}\right)^{2}\left(\frac{z_{\sigma}}{z_{d}}\right)^{2},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \Omega_{\textsc{gw}}(f_{d})=\Omega_{{\rm r}0}\left(\frac{H_{1}}{M_{\rm Pl}}\right)^{2}z_{d}^{|3-2\beta|-3}.$$ (B.12)

In the first case, given the constraints (B.8), the maximal amplitude can be reached for the limiting values $z_{d}\simeq z_{\sigma}$, which imply however $H_{d}\simeq H_{\sigma}$: hence, in that case, the axion starts decaying as soon as it becomes dominant, and there is not enough time for an efficient curvaton mechanism. Also, in that case, the maximal amplitude would correspond to a frequency range $f\sim f_{1}$, in general too high for the sensitivity of present detectors.

In the second case with the peak at $f_{d}$, given again the constraints (B.7) and (B.8), the maximal amplitude can be obtained either in the limit $z_{d}\simeq z_{\sigma}\rightarrow 1$ or in the limit $\beta\rightarrow 0$. Discarding the first possibility (for the same reasons as before), in order to find the allowed region for the GW signal of maximal intensity we will thus concentrate on the limiting case $\beta=0$ which, as we will see, automatically leads to a peak located in frequency ranges possibly accessible to third-generation detectors.

It should be noted, in addition, that the limiting amplitude reached at $f_{d}$ for $\beta=0$ is only controlled by the ratio $H_{1}/M_{\rm Pl}$, whose maximal allowed value is bounded by the constraints (B.9): hence, the amplitude of the GWB approaches its allowed maximum in the limit in which the condition (B.9) is saturated by the equality

$$\log_{10}z_{\sigma}=\frac{2}{7}\log_{10}\left(\frac{H_{1}}{M_{\rm Pl}}\right)+\frac{3}{7}\log_{10}z_{d}\,.$$ (B.13)

This result has two important consequences.

First of all, by using eq. (4.10), we can check that the above condition is equivalent to the condition $\sigma_{i}=M_{\rm Pl}$, and this uniquely fixes the transfer function (B.3) leading to the constant numerical value $\mathcal{T}^{2}\approx 0.137$. Second, by inserting into the above condition the general expression (B.5) for $H_{1}$, and solving for the variable $\log_{10}z_{\sigma}$, we can eliminate $z_{\sigma}$ everywhere and confine our discussion of the maximum allowed spectrum to a two-dimensional parameter space spanned by the variables $\{\log_{10}z_{s},\log_{10}z_{d}\}$, with $\beta=0$, $z_{\sigma}$ given by eq. (B.13) and $H_{1}$ given by

$$\log_{10}\left(\frac{H_{1}}{M_{\rm Pl}}\right) \simeq \frac{14}{37-9n_{\rm s}}\left[\log_{10}\left(\frac{4.2\pi^{2}}{0.137}\right)-9+(1-n_{\rm s})(\log_{10}1.5-27)\right.$$ (B.14)

$$+\left.(1-n_{\rm s})\log_{10}z_{s}+\frac{2}{7}(1-n_{\rm s})\log_{10}z_{d}\right].$$

We can now easily impose all constraints (B.8)–(B.11) and evaluate, in such a context, both the allowed region of parameter space and the maximal allowed value of the peak amplitude. It turns out (see Fig. 3), that the maximal value of eq. (B.14) compatible with the given constraints corresponds to $\log_{10}(H_{1}/M_{\rm Pl})\approx-3.29$, so that the expected maximum intensity of the primordial GWB is given by

$$\Omega_{\textsc{gw}}^{\rm max}=\Omega_{{\rm r}0}\left(\frac{H_{1}}{M_{\rm Pl}}\right)^{2}\approx 10^{-10.6}\,.$$ (B.15)

The allowed values of the parameters compatible with this maximal intensity (and with the imposed constraints) are in the range $18.7\leavevmode\nobreak\ \raise 1.72218pt\hbox{$<$}\kern-8.00003pt\lower 2.6694pt\hbox{$\sim$}\leavevmode\nobreak\ \log_{10}z_{s}\leavevmode\nobreak\ \raise 1.72218pt\hbox{$<$}\kern-8.00003pt\lower 2.6694pt\hbox{$\sim$}\leavevmode\nobreak\ 22.3$ and $2.19\leavevmode\nobreak\ \raise 1.72218pt\hbox{$<$}\kern-8.00003pt\lower 2.6694pt\hbox{$\sim$}\leavevmode\nobreak\ \log_{10}z_{d}\leavevmode\nobreak\ \raise 1.72218pt\hbox{$<$}\kern-8.00003pt\lower 2.6694pt\hbox{$\sim$}\leavevmode\nobreak\ 14.9$, as shown in Fig. 3. The corresponding value of $z_{\sigma}$ is given by eq. (B.13) and lies in the range $0\leavevmode\nobreak\ \raise 1.72218pt\hbox{$<$}\kern-8.00003pt\lower 2.6694pt\hbox{$\sim$}\leavevmode\nobreak\ \log_{10}z_{\sigma}\leavevmode\nobreak\ \raise 1.72218pt\hbox{$<$}\kern-8.00003pt\lower 2.6694pt\hbox{$\sim$}\leavevmode\nobreak\ 5.5$.