ournal of C osmology and A stroparticle hysics J The rise of the primordial tensor spectrum from an early scalar-tensor epoch

. Primordial gravitational waves (PGW) produced during inﬂation span a large range of frequencies, carrying information on the dynamics of the primordial universe. During an early scalar-tensor dominated epoch, the amplitude of the PGW spectrum can be enhanced over a wide range of frequencies. To study this phenomenon, we focus on a class of scalar-tensor theories, well motivated by high energy theories of dark energy and dark matter, where the scalar is conformally and disformally coupled to matter during the early cosmological evolution. For a conformally dominated epoch, the PGW spectrum has a ﬂat step-like shape. More interestingly, a disformally dominated epoch is characterised by a peaked spectrum with a broken power-law proﬁle, with slopes depending on the scalar-tensor theory considered. We introduce a graphical tool, called broken power-law sensitivity curve, as a convenient visual indicator for understanding whether a given broken power-law proﬁle can be detected by GW experiments. We then analyse the GW spectra for a variety of representative conformal and disformal models, discussing their detectability prospects with the Einstein Telescope (ET), Laser Interferometer Space Antenna (LISA), DECi-hertz Interferometer Gravitational wave Observatory Big Bang Observer


Introduction
Scalar-tensor theories are ubiquitous in extensions of the Standard Model (SM) of particle physics and cosmology. For example, string theory approaches to particle physics and inflationary model building generically predict the presence of several new ingredients, in particular, new particles such as scalar fields with clear geometrical interpretations. These scalars couple conformally and disformally to matter living on branes, extended objects where matter is localised. In string D-brane constructions, longitudinal string fluctuations are identified with the matter fields such as the SM and/or dark matter (DM) particles, while transverse fluctuations correspond to scalar fields. Constraints on scalar fields and its couplings to matter today are extremely tight, for example from solar system tests [1][2][3][4], and recently from strict bounds on the speed of gravitational waves [5].

JCAP08(2022)010
On the other hand, the thermal history of the early Universe remains uncertain. Its evolution involves a sequence of epochs, each characterised by a certain expansion rate H, which is key to understanding the physics of processes occurring during these eras. According to current observational evidence, the universe was radiation dominated at the time of Big Bang Nucleosynthesis (BBN) and, most likely, there was an early era of vacuum domination known as inflation. However, there might have been a non-standard cosmological history between the end of inflation and the onset of BBN.
An epoch of a scalar-tensor modification to general relativity (GR) domination before the onset of BBN, can change the cosmological expansion rate, without violating present constraints on scalar fields and its couplings to matter. For example, a modified expansion rate felt by matter can change the standard predictions for the dark matter relic abundance as studied in [6][7][8][9][10][11][12][13][14][15][16][17][18][19] for the conformal case, and in [20] for the disformal case.
In this work, we investigate the imprints of a scalar-tensor dominated epoch with conformal and disformal couplings, on the primordial gravitational wave spectrum produced during inflation. This spectrum spans a large range of frequencies with an almost scale-invariant profile, whose amplitude is too small to be detected by future gravitational wave experiments. We show how this signal is enhanced during an era of scalar-tensor domination with distinct signatures in the conformal and disformal cases.
The conformal case, which characterises the Brans-Dicke class of scalar-tensor theories, was recently discussed in [21] for a particular choice of conformal coupling. In this case, the PGW flat spectrum is enhanced to a flat spectrum with a larger amplitude during the scalar-tensor epoch. We demonstrate that this step-like enhancement depends on the conformal factor (and initial conditions), offering the possibility of probing its existence using correlations between different experiments.
Remarkably, during a disformally dominated epoch, the PGW spectrum has a characteristic peak with a distinctive frequency profile, which offers a smoking-gun signature of the early scalar-tensor epoch. In order to easily understand in a visual manner whether any given broken power-law GW spectrum can be detected by a GW experiment, we introduce the notion of broken power-law sensitivity curve (BPLS).
The paper is organised as follows. In section 2, we review the calculation of the primordial gravitational wave spectrum in standard cosmology, following closely [21][22][23][24]. In section 3, we introduce the general modified cosmological setup due to scalar-tensor theories motivated by D-branes, as well as a more phenomenological set-up, popular in the literature. We focus on the calculation of the modified expansion rate for different cases of interest: a purely conformal case, a purely disformal case, and a conformal-disformal case in D-brane-like scenarios. In section 4, we discuss the rise of the primordial gravitational wave spectrum due to a scalar-tensor epoch in the three cases discussed in section 3. Specifically, subsection 4.1 focuses on the conformal rise of the PGW spectrum. Subsection 4.2 discusses the interesting disformal case, which gives rise to a characteristic peaked spectrum. In passing, we explain how our results differ with a similar peaked spectrum arising from a short kination era due to a spinning axion. In subsection 4.3, we introduce and discuss the concept of the broken power-law sensitivity curves. After our conclusions in section 5, appendix A contains more details on the calculations in section 3.

Primordial gravitational waves in standard cosmology
In this section, we briefly review the PGW spectrum evolution in standard cosmologies [21][22][23][24]. Given a model of inflation, a stochastic background of GWs arises inevitably from the primordial tensor fluctuations, whose equation of motion is given bÿ where H is the Hubble expansion rate in GR and we consider that the evolution of the primordial tensor fluctuations is source-free for the frequencies we are interested in, f 10 −10 Hz, corresponding to temperatures T 4 MeV. To solve the tensor perturbation equation, one can write it in Fourier space as with λ = +, × corresponding to the two independent polarisations, and λ being the spin-2 polarisation tensor satisfying the normalisation condition λ ij λ * ij = 2δ λλ . Using (2.2), the solution to (2.1) can be written as where T is a transfer function and h λ inf ( k) is the amplitude of the tensor perturbation. The energy density of the relic GW today is given by (2.4) The relic density of the PGWs from the tensor perturbation is calculated from dρ GW (t, k) d ln k , (2.5) where ρ c is the critical density of the Universe. This can be further rewritten using (2.3) as [21][22][23][24] Ω GW (τ, k) = P t (k) (T τ (τ, k)) 2 12a 2 (τ )H 2 (τ ) 1 24 where the subscript 'τ ' denotes derivatives with respect to conformal time, i.e. dτ = dt/a, and the subindex 'hc' denotes horizon crossing [21], i.e. k = a hc H hc . The primordial tensor spectrum P t (k) is determined by (2.7) The fractional energy density in primordial gravitational waves, as observed today, can be expressed as Ω 0 GW (k) 1 24

JCAP08(2022)010
where H 0 denotes the Hubble parameter today. Using entropy conservation, we can write a hc a 0 = g * s,0 g * s,hc where g * s represents the entropy degrees of freedom (see subsection 3.1.1 for more details).
Combining the above two equations, we obtain (2.10) The amplitude of the primordial tensor power spectrum is given by P t = r A S , where A S = 2.1 × 10 −9 is the amplitude of the primordial scalar power spectrum [25], and we consider the latest upper bound for the tensor-to-scalar ratio reported by BICEP/Keck, r = 0.036 [26]. Now, the frequency, f 0 , of GWs observed today, at which the corresponding mode k = 2πf 0 a 0 (where a 0 is the scale factor today) reenters the horizon, can be related to the temperature T hc via the following relation [22,23,27,28]: where ρ hc is the total energy density of the universe at horizon crossing. Using this equation, we can obtain the evolution of Ω GW as a function of f 0 in section 4.

Scalar-tensor theories
In this section we introduce the general scalar-tensor setup describing a D-brane like scalartensor theory with conformal and disformal couplings motivated by D-brane world scenarios, as well as a more phenomenological setup, where these couplings are in principle not related to each other. We follow closely [19,20], extending and generalising the results found there. The starting action is given by: where where κ 2 = M −2 Pl = 8 π G, b = 0, 1, depending on whether we are in a phenomenological (b = 1) or D-brane (b = 0) set-up. We will see below, how each different set-up affects the results. M is a mass scale, which can be related to the tension of the D-brane [19,20], and we take M = 0 when we consider the phenomenological model.

JCAP08(2022)010
The disformally coupled metric,g µν , is given by 1 where C(φ) and D(φ) are the conformal and disformal couplings of the scalar to the metric, respectively. 2 These functions are in principle arbitrary up to a causality constraint, which requires C > 0 and C + 2DX > 0 [29]. In a D-brane motivated set-up, these functions are related via M 4 CD = 1 (see appendix C of [19] for details) and automatically satisfy the causality constraint.
Cosmological equations. Consider an homogeneous and isotropic FLRW metric given by with a(t) being the scale factor. In this background, the equations of motion become

5)
where H =ȧ/a, the dots indicate derivatives with respect to t, and the primes indicate derivatives with respect to the field φ. The Lorentz factor is given by and we have used the equation of state, P = ωρ, with P and ρ being the pressure and energy density, respectively. For the scalar field, the energy density ρ φ and the pressure P φ are given by [19,20] for details). The energy-momentum conservation equation gives ∇ µ T µν tot = ∇ µ T µν φ + T µν = 0. That is, the scalar field and matter are not separately conserved in the Einstein frame. The time component of this constraint yields the equationṡ

JCAP08(2022)010
Using the last equation (3.12), Q 0 can be rewritten as Plugging this into the (non)conservation equation for matter (3.12), we obtaiṅ We are interested in the disformal, or Jordan, frame, which is defined as the frame where matter and entropy are conserved. That is,∇ µT µν = 0, where∇ µ is computed with respect to the disformal metric, and the energy-momentum tensor is defined bỹ (3.15) where the tilde denotes the Jordan frame. This can be related to the energy-momentum tensor in the Einstein frame through the following relation: Then the energy densities, pressures, equations of state, and the scale factors in the Einstein and Jordan frames are related through the following expressions: The Hubble parameter in the Jordan frame is given bỹ with dt = C 1/2 γ −1 dt, so that it is computed from H and φ. As indicated above, in the Jordan frame, the continuity equation for matter takes the standard form, that is [19]: Note that if we consider a purely disformal case, with C = 1, then the scale factor in both frames coincide. Moreover, at the onset of BBN, whrn C = γ = 1, the two frames coincide.

Modified expansion rate
We are interested in computing the modified expansion rate in the Jordan frame to compare it with the standard GR evolution. For this purpose, it is convenient to introduce the dimensionless scalar ϕ = κ φ, and swap time derivatives with derivatives with respect to the -6 -

JCAP08(2022)010
number of e-folds, i.e. N = ln(a/a 0 ), so that dN = Hdt. With these changes, we can rewrite equations (3.5)-(3.7) as follows where the subscript N denotes a derivative with respect to the number of e-folds, we have used (3.17) to replace the energy density and equation of state in the Jordan frame, and we have defined and λ = V/ρ (=Ṽ/ρ). Also, in terms of ϕ and N -derivatives, the Lorentz factor is given by Now, the expansion rate in General Relativity (GR), which we express is given by (κ GR will be specified below) using (3.17) we can write H GR entirely as a function of H, ϕ, ϕ N as follows (see [6,19,20]): Therefore, once we find a solution for H and ϕ, we can compare the expansion ratesH with H GR , by introducing the following parameter ξ, which measures the departure from the standard expansion: with κ 2 GR being the gravitational constant and κ 2 being the value of the gravitational constant measured by local experiments for conformally coupled theories, and ϕ 0 is the value of the scalar field at the present time. Notice that ξ can be larger or smaller than one, indicating an enhancement or reduction ofH with respect to H GR . This means thatH can grow during the cosmological evolution. 3 Below we solve the system of equations (3.20), (3.21), and (3.22) numerically. For this, we first change the number of e-folds to the Jordan frame,Ñ , which is given by [19] N ≡ ln a a 0 = ln Thus, in terms ofÑ , the derivatives of ϕ can be expressed as

Degrees of freedom and the kick function
When the conformal factor is turned on, α(ϕ) = 1 and thus the equation of motion for the scalar field (3.22), has a term, α(ϕ)(1 − 3ω), which depends non-trivially on the equation of stateω, and which acts as an effective potential, V eff = ln C 1/2 , whenω = 1/3. Thus, deep in the radiation dominated era, whenw ∼ 1/3, this term drops out. As the universe cools down, when the temperature of the universe drops below the rest mass of a particle species, the particle becomes non-relativistic, and there arise small departures in the value ofw from w = 1/3. This leads to non-zero contributions due to this term, (1 − 3w) in eq. (3.22), which generates a kick in the scalar field as we will see. This factor is referred to as the 'kick function'. We now compute the kick function referred to above, which we shall use to solve eqs. (3.21) and (3.22) numerically. For a particle species 'A' that is in thermal equilibrium with the radiation bath during the radiation dominated era, its energy density and pressure are given by [22,23,30,31]  where x = (E/T ) 2 − y 2 A , y A = m A /T , and g A is the number of degrees of freedom for each species 'A'. The total energy density and pressure is equal to the sum of energy densities and pressures of all these particles:ρ The total number of entropy degrees of freedom is given by (see figure 1) The evolution of this quantity is essential to relate the scale factor and the temperature of the universe through entropy conservation (2.9). We evaluate the energy density, pressure, and the total number of entropy degrees of freedom numerically for all the particles in table 1. We then add each of these quantities to the corresponding quantities of the relativistic particles. Before the QCD phase transition, the relativistic species of interest are the gluons (g A = 16), photons (g A = 2), light quarks (g A = 36), and neutrinos (g A = 6). After the QCD phase transition, i.e. below 170 MeV, the relativistic species that remain are the photons and neutrinos. Also, after neutrino decoupling, i.e. below 1 MeV, we take into account the fact that the neutrino temperature evolves differently compared to the photon temperature. This can be estimated using conservation of entropy.
The kick function,Σ, during the radiation, matter dominated and Λ dominated epochs is given byΣ  where we have consideredp M = 0,p Λ = −ρ Λ , and the equation of state parameter is given byw = (1 −Σ)/3. So the Hubble parameter in GR that we use below, is given by the total energy density as andρ R is given by eq. (3.34a), while ρ 0 , Ω M 0 , and Ω Λ 0 are the energy density, the matter density parameter, and the dark energy density parameter evaluated today, respectively. Note that there is no contribution from the scalar field energy density in the Jordan frame. In the Einstein frame on the other hand, ρ total = ρ R + ρ M + ρ Λ + ρ ϕ .
We now focus on the two cases of interest separately, namely (b = 1, M = 0) and (b = 0, M = 0), to study the expansion rate modification during an early scalar-tensor theory domination in the universe's evolution.

Phenomenological case
In this case, we take M = 0 and b = 1 in (3.2b). That is, the scalar field has a standard kinetic term and there is in principle no particular relation between C and D, except the causality constraint [29]. Below we consider the cosmological evolution for the following cases: a purely conformal case, that is D = 0 (γ = 1), and a purely disformal case, that is C = 1, D = 0.

Purely conformal enhancement
In the purely conformal case, we have D(ϕ) = 0, hence γ = 1. We further take λ ∼ 0, such that the dark energy today is fully dominated by a cosmological constant. In this case, the As we discussed above, deep in the radiation dominated era,w ∼ 1/3, hence the last term in (3.39) drops out and the equation can be solved analytically, to give ϕ N ∝ e −N i.e. the field velocity decreases rapidly in the radiation dominated era. Further integration shows ) . Thus the field settles to a constant value after a few e-folds [6,31,32]. This behaviour holds, even for larger values of ϕ i N , however for values close √ 6, the approximation ∆ϕ ϕ i N stops being valid, and the constant value to which the field settles differs largely from ϕ i N . We shall encounter this behaviour explicitly in the full numerical solution of (3.39) (see figure 6).
As the universe cools down, when the temperature of the universe drops below the rest mass of a particle species, the particle becomes non-relativistic, and there arise small departures in the value ofw fromw = 1/3. This leads to non-zero contributions due to (1 − 3w) in the last term in eq. (3.39), which generate a kick in the scalar field with α(ϕ) acting as an effective potential: V eff = ln C 1/2 .
We now consider a set of conformal functions motivated by the choice in [6], which satisfy the requirement that the standard cosmological evolution is recovered at the onset of BBN. Moreover, this choice will allow us to demonstrate the dependence on this choice (and the initial conditions) on the enhancement of the expansion rateH, compared to the standard GR case, H GR . In section 4, we will discuss how this impacts the PGW signal and the potential for its detection.
Conformal expansion rate enhancement. We consider a suitable modification of the conformal factor used in [6,19], given by with b = 0.1 and β = 8. In [6,19], n = 1 was chosen such that the evolution of the Hubble parameter matches that of the standard GR evolution after BBN. Here we choose n = 1, 2, 4 in order to demonstrate the dependence of the signal on the conformal factor. The effective potential is a runaway of the form V eff = n ln(1 + b e −βϕ ). As discussed above, deep in the radiation era, any initial velocity ϕ N , goes rapidly to zero and ϕ → constant value. As soon asω differs slightly from 1/3, the effective potential kicks in, and the field rolls along it untilω ∼ 1/3 again [6,19,31,32]. During the following matter and dark energy dominated eras,ω = 1/3 and the field keeps rolling down its effective potential, weighted by 2(1 − 3ω) (see (3.39)).
The enhancement of the expansion rate with respect to GR depends on the initial conditions and the conformal factor, and thus can in principle be probed depending on its effects, e.g. on the dark matter relic abundance [6,19], or on the stochastic gravitational wave background, as we shall explore in the next section. In general, both the initial position and velocity of the scalar field can take any value, positive or negative. We choose initial conditions for the scalar field and its velocity to be less than the Planck scale, that is, in the range (ϕ i , ϕ i N ) ∈ (±1, ±1). In the runaway effective potential dictated by ln C 1/2 for the conformal factor (3.40), there are the following possibilities: (a) The scalar field starts somewhere up in the runaway effective potential with zero or positive initial velocity, thus staying at or reaching a constant value until the first particles become non-relativistic, giving a kick to the scalar field, so that it rolls down its potential untilω ∼ 1/3 again, and thereafter stays constant until the kick function releases it again during the radiation era. It then evolves rapidly until C reaches 1, well before BBN. It is thus clear that the largest enhancement will occur when the field starts as high as possible in the effective potential, where the conformal factor will be the largest.

(b)
The second possibility arises when the initial velocity is negative. In this case, deep in the radiation era, the field quickly reaches a constant negative value until the first particles become non-relativistic, turning on the effective potential V eff . At this point, the field starts rolling down its potential with kicks dictated byΣ (3.36). The field subsequently rolls down its effective potential during the matter and dark energy dominated eras. The enhancement is thus dictated by the smallest constant negative value reached by the scalar deep in the radiation era. The subsequent evolution proceeds as in the previous case.
We numerically calculate the evolution and enhancement using both types of initial conditions, which are summarised in tables 2 and 3 below. In the first case, we choose the initial conditions with ϕ iÑ = 0 at the initial temperature of 10 15 GeV, with ϕ i as in table 2. We start the evolution of the scalar field deep in the radiation era. The evolution of ϕ and the corresponding behaviour of the conformal factor (3.40) are shown in figure 3. As discussed before, deep in the radiation era, the scalar stays at a constant value, after which, when the first particles become non-relativistic, the effective potential turns on asω = 1/3, and the field starts to roll down the runaway. The conformal factor starts at a large value set by the initial value of ϕ, eventually dropping back to unity well before the onset of BBN. We see that larger powers of n, give a larger value of C and thus larger enhancement. This is also reflected in the Hubble parameter, figure 4. Table 2. Type (a) initial conditions with ϕ iÑ = 0, for the purely conformal case with conformal function (3.40). The initial temperature in all the cases isT i = 10 15 GeV.  Using the solutions for ϕ and C, we compute the modified expansion rate in the Jordan frame and compare it to H GR (see (3.29)) to get: 41) where ϕ 0 is the value of the field today, and B is given by (3.23) with M = 0. The behaviour of the Hubble parameters is shown in figure 4. The notch observed in this plot is due to the fact that, in the Jordan frame, the Hubble parameter can become smaller than H GR for a brief period of time when (1 + α(ϕ) ϕ N )/ √ B < 1, as observed in [19] (see figure 5). As we can see from figure 4, the relative enhancement due to the different powers of n in (3.40) is not very prominent.
Next, we consider the initial conditions with non-zero initial velocity ϕ iÑ = 0, and initial temperature of 10 15 GeV as shown in table 3. The resultant evolution of ϕ and the corresponding behaviour of the conformal factor (3.40) are shown in figure 6. After an initial negative velocity, the scalar field settles down to a constant value, 4 after which, when the first particles become non-relativistic, the effective potential turns on asω = 1/3, and the field starts to roll down the effective runaway potential. The conformal factor reaches a larger maximum value than in the case with zero initial velocity, during the radiation era, eventually dropping back to unity well before the onset of BBN. 4 Note that the approximation ∆ϕ ∼ ϕ i N is not valid as the initial velocity is of order √ 6.    less than one for a brief period of time, thus making the Hubble parameter in the Jordan frame smaller than H GR [19]. From eq. (3.41), it is evident that the amplitude of the Hubble parameter at any temperature is directly proportional to the conformal factor. This effect directly impacts the enhancement of the amplitude of the gravitational waves as we discuss this in the next section. From figure 7, one can easily see that the enhancement is largely amplified. Further, we can see that the larger the conformal factor, the larger the enhancement of the expansion rate.

Purely disformal case
In the purely disformal case, we have C(ϕ) = 1. We study two different scenarios -one wherein the disformal factor is a constant, and one wherein the disformal factor is a function of the scalar field ϕ. Since in the phenomenological case, the only relation between the conformal and disformal factors is given by the causality constraint, we have to ensure that we choose the relevant parameters such that the following condition is always satisfied: In terms of φ = κϕ andÑ derivatives, this simplifies to which will need to be satisfied by our choice of D(ϕ). The equations of motion (see eqs. (3.21), (3.22)) in this case simplify to (note that in the pure disformal case, N =Ñ ): Constant disformal factor: D = D 0 . We start by considering the simplest possibility of a constant disformal function, D = D 0 . In this case, δ = 0 and the last term in (3.44b) vanishes. The initial condition for the Hubble parameter is found by finding a real positive solution to the cubic equation for H in (3.28) for which γ i ∼ 1 (see [19] and appendix A for details). This imposes the following condition on D 0 : which is further complemented by the causality condition (3.43) on D 0 . Using these two constraints, we choose the parameters and initial conditions as given in table 4. The evolution of the scalar field ϕ and the Lorentz factor γ are shown in figure 8 (dashed lines). As we see there, the scalar field stays constant for a few e-folds, but quickly evolves towards larger values, causing γ to increase, before going back to one well before the end of BBN. This causes an enhancement of the Hubble parameter, according to ξ = γ 3/2 /B 1/2 , as shown in figure 9. It is important to point out that the choice of the initial temperature impacts the moment at which the enhancement occurs. Thus, for a larger initial temperature, the enhancement will occur earlier [20]. Interestingly, as we shall discuss in the next section, we can probe this with gravitational waves. Field dependent disformal factor: D = D 0 ϕ 2 . In this case, δ(ϕ) = 0, and we cannot neglect the last term in (3.44b). We again proceed to set the initial condition for the Hubble parameter by finding the real positive solution of (3.28) for which γ i ∼ 1 (see [19] and appendix A for details). This imposes the following condition Further, the causality condition (3.43) implies that we must have Taking into account these constraints, we choose the parameters and initial conditions as given in table 4. As in the previous case, by starting the evolution at a higher temperature, the enhancement of the Hubble parameter will occur earlier, which will impact when the GW background will be enhanced. We shall see this explicitly in the next section. The evolution of the field ϕ and the Lorentz factor γ are shown in figure 8 (dotted line). The Hubble parameters for the constant and field dependent cases are shown in figure 9, (dashed and dotted lines respectively). Compared to the constant case, the field dependent example results in a slight increase in the maximum of γ, which thus is reflected in a larger enhancement in the Hubble parameter with a slightly different profile, which will be reflected also in the PGW effect.

D-brane scalar tensor theories
This case corresponds to b = 0 in (3.2b), and it arises in D-brane cosmology scenarios.
The scalar field has a non-standard kinetic term dictated by the Dirac-Born-Infeld (DBI) action, (3.2b) and the modifications to the expansion rate and its effect on the dark matter relic abundance was discussed in [20]. This type of scenario can arise from a post-string inflationary scenario. At this scale, the universe is already four-dimensional and moduli associated to the compactification have been stabilised. 5 However, in what follows, our JCAP08(2022)010  study is purely phenomenological and can be used as a first step to understand the effects for gravitational waves in D-brane scalar tensor theories in the early universe. As shown in appendix C of [19], the canonical normalisation of φ, obtained by expanding the DBI action, implies a relation between the conformal and disformal factors through M 4 CD = 1. Thus, in this section, we study the solutions for the D-brane conformally and disformally coupled matter with the choice above, which implies δ(φ) = −α(φ) (see eqs. (3.24) and (3.25)).

Evolution equations.
In this case, the evolution equations (3.21), (3.22) simplify to where γ is given by (3.26), and here B is given by In what follows we consider the expansion rate modification in the purely disformal case, wherein C = 1 and therefore D = 1/M 4 , and in the conformal plus disformal case wherein C = 1 and D = 1/CM 4 with C given by (3.40). Remember that, in this case, the conformal and disformal couplings are related and thus we cannot have a purely conformal case as in the phenomenological example above. Moreover, the causality constraint is always satisfied.

Purely disformal case
In the purely disformal case, C = 1 and thus D = 1/M 4 (= D 0 ). As in the phenomenological case, the initial condition for the Hubble rate H can be obtained by looking for the positive real solution of the cubic equation for H (see appendix A and [19,20] for details), given the initial values for ϕ i , ϕ iÑ , M . These initial conditions can be used together to obtain the initial value of the Lorentz factor, γ i , which needs to be of order of one, in order to satisfy the constraints on the present value of the Hubble parameter.
With this in mind, we have chosen M = 1.200 × 10 5 GeV, as well as the other parameter values and initial conditions as given in table 4. The evolution of ϕ and γ are shown in figure 8 (solid line), while the Hubble parameter is shown in figure 9 (solid line). We see that, compared to the phenomenological examples above, the enhancement in the pure Dbrane disformal scenario is larger, for the same initial conditions. Moreover, the γ profile also differs slightly. This behaviour has interesting implications for the PGW spectrum discussed in the next section.

Conformal-disformal case
We now turn on a non-trivial conformal factor, C(ϕ). We choose the same conformal function as in the phenomenological case, namely, (3.40). In the present case, this immediately fixes D through the condition M 4 C D = 1. That is, in this case, both conformal and disformal factors are turned on. Therefore, we shall have to solve the full set of coupled equations (3.48). Initial conditions and evolution. As in the phenomenological case, we can either choose a suitable initial field value and a zero initial velocity (see discussion in section 3.2.1) or a nonzero initial velocity. In the present case, if the initial velocity is zero, the initial condition for the Hubble rate With this in mind, we again consider two sets of initial conditions: (a) Zero initial velocity: analogous as to the phenomenological case, we can start with a zero velocity and a suitable value for the scalar field. We choose the same initial values as in the phenomenological case (see table 2 .21)). Note that, in order to solve the equations (3.48) in terms ofÑ , we have to convert the initial conditions obtained by the above procedure into initial conditions on the set of values of ϕ i , ϕ iÑ , M . With this in mind, we consider the initial conditions as mentioned in table 5. The initial conditions on ϕ and ϕ N are different from what we had chosen for the same conformal factor in the phenomenological case. This is because the equations (3.48) are different from the phenomenological case, and the same choice of initial conditions is not numerically viable in both cases. Since starting with a large negative initial velocity leads to a greater enhancement in the conformal factor, we choose the initial velocity accordingly and then choose the initial condition on the field such that the equations (3.48) can be numerically solved. The evolution of the field and the conformal factor is shown in figure 13, while the Hubble expansion rates are shown in figure 14. As we can see, compared to the previous choice of starting with zero initial velocity, in this case, C ∼ 1 at the start of the evolution. It subsequently increases to its maximum value and then drops back quickly to one well before the onset of BBN. The behaviour is thus very similar to the phenomenological case, and it is dominated by the conformal evolution. In fact, during the whole evolution, γ ∼ 1.
Comparing figures 10 and 13 with figures 3 and 6, we find that the evolution of the field, the conformal factors, and the expansion rates are similar overall in the purely conformal phenomenological scenario and the conformal-disformal D-brane scenario. This implies that, even in the presence of a disformal factor, the field evolution is dominated by the conformal factor. In order to better illustrate the different behaviour in the two scenarios, in figure 15, we show the evolution of the scalar field and the conformal factor for the purely conformal phenomenological case (solid line), and the conformal-disformal D-brane case (dashed line), for n = 1 and initial conditions as in tables 3 and 5 respectively. As we see, the phenomenological case seems to be more effective in enhancing the Hubble rate and thus the PGW as we will see.

The rise of the primordial tensor spectrum
In section 2, we briefly described the estimation of the fractional energy density of gravitational waves in a standard cosmological scenario, given by eq. (2.8). In the previous section 3, we saw that during an epoch of scalar-tensor domination, the expansion rate is modified, thus inducing a non-trivial modification in the PGW energy density as follows: whereã andH are the modified scale factor and Hubble parameter in the Jordan frame, which is the frame in which energy density and entropy are conserved. As we have seen in the previous section, the Jordan frame is indistinguishable from H GR at the onset of BBN, when the conformal and/or disformal factors become one. Using (3.29) and (3.30), (4.1) becomes where γ is given by (3.26) and B by (3.23). We now discuss the rising of the PGW spectrum due to the modified expansion rate during a period of scalar-tensor domination epoch, in the post-inflationary evolution. In subsection 4.1, we consider the characteristic step-like enhancement of the PGW due to a conformal coupling in the phenomenological and D-brane scalar-tensor theories discussed in subsections. 3.2.1 and 3.3.2. Then, in subsection 4.2.2, we discuss the non-trivial rise of the PGW spectrum due to a disformal coupling in the phenomenological and D-brane scalar-tensor theories discussed in subsections 3.2.2 and 3.3.1.

Conformal enhancement of the PGW spectrum
In this section, we consider the effect of the conformal coupling on the expansion rate and thus in the PGWs. As discussed previously, in the phenomenological model, the disformal and conformal couplings are not related, except via the causality condition, (3.42) and this case corresponds to setting D = 0. On the other hand, for the D-brane case, these couplings are related through M 4 CD = 1. We have also seen in section 3.3.2 that, when both couplings are turned on, the disformal effect is negligible and the behaviour of the expansion rate is dominated by the conformal term.
As also discussed in the previous section, in the presence of a conformal coupling, an effective potential is turned on when the temperature of the universe drops below the rest mass of a particle species so that it becomes non-relativistic, thereby causing a small departure fromw = 1/3, which generates a 'kick' in the scalar field. The temperature when the kick function starts to trigger a change in the evolution of the field is approximately given Hence, for frequencies above this, there is a rise and eventual enhanced flat spectrum for Ω 0 GW (k) h 2 for both sets of initial conditions discussed in sections 3.2 and 3.3.2, namely, the scalar field starting with zero and non-zero initial velocities. As we have already illustrated in these sections, the enhancement in the Hubble parameter can be increased (or decreased) depending on the choice of the conformal factor.

Phenomenological case
In the phenomenological conformal scenario, D = 0 and thus γ = 1. Therefore, eq. (4.2) simplifies toΩ where B is given by (3.23) with b = 1, M = 0. From this equation, we see that an enhancement is determined mostly by the amplitude of the conformal factor C(ϕ). The quantity (1 + α(ϕ) ϕ N ) 2 /B reaches a large value ∼ 29 at very early temperatures but then drops rapidly to O(1). As shown in figure 5, this quantity exhibits a dip in its value when the equation of state starts changing, which also results in a dip in the PGW spectrum before it gets amplified due to the conformal factor. In figure 16 (left panel), we show the resulting step-like rise of the gravitational wave amplitude for the three choices of n and the corresponding initial conditions in table 3.

D-brane conformal-disformal enhancement
When the conformal coupling is turned on in the D-brane scalar-tensor case, the disformal coupling is also turned on according to the relation M 4 CD = 1. Thus, in this case, γ = 1, and (4.2) can be rewritten as where B is given by (3.23) with b = 0. However, as we have seen in section 3.3.2, the modification of the expansion rate is driven mostly by the conformal coupling. As we saw, throughout the entire evolution, γ ∼ 1, and thus the rise in the PGW spectrum is driven by the conformal factor. The dip in the value of (1 + α(ϕ)ϕ N ) 2 /B (cf. figure 11) is reflected in a dip also in the corresponding PGW spectrum. In figure 16 (right panel), we show the rise of the PGW spectrum for the three choices of n in (3.40) and the corresponding initial conditions in table 5. As in the phenomenological case, the enhancement has a step-like flat spectrum. Thus the profile for both, phenomenological and D-brane like is the same. Note that in [21] the profile for the PGW had also a step-like flat behaviour, for a the different conformal function, C = e βφ 2 . We thus expect a step-like flat behaviour for other choices of the conformal function.
From figure 16, we can see that, in the phenomenological pure conformal example, where C and D are unrelated (except via the causality condition), the amplification of the PGW spectrum is larger compared to the D-brane model, where C and D are related through M 4 CD = 1. Recall that the choice of the conformal function is exactly the same in both cases. This effect may be caused by the disformal contribution.
-24 -JCAP08(2022)010 Figure 16. The conformally enhanced gravitational wave spectra have been plotted as functions of frequency. On the left panel, we show the PGW spectra for the phenomenological pure conformal case for the set of initial conditions in table 3. On the right panel, we show the PGW spectra for the D-brane conformal-disformal case for the set of initial conditions in table 5. For both cases, the solid, dashed, and dotted lines correspond to n = 1, n = 2, and n = 4 respectively. The power-law sensitivity curves have been plotted using the method described in [34] (more on this in section 4.3).

Disformal rise of the PGW spectrum
When only the disformal coupling is turned on, that is, C = 1, eq. (4.2) simplifies to: where γ is given by (3.26) and B by (3.23). Thus, the amplitude of the gravitational waves depends only on the Lorentz factor γ and B as does the modified Hubble parameter. As we discussed in sections 3.2.2 and 3.3.1, the position of the maximum enhancement of the Hubble parameter -and thus the gravitational wave amplitude -depends on the initial temperature. The scalar field quickly evolves, producing a large peak in γ, eventually settling to a constant value well before the onset of BBN. Thus, the earlier the initial temperature is chosen, the earlier the enhancement of H occurs, which is reflected in the characteristic peaked enhancement at a higher frequency in the gravitational wave spectrum. As we shall discuss now, the disformal coupling in both these cases gives rise to a characteristic peak in the PGW spectrum, with a peculiar frequency dependence. Moreover, the phenomenological and disformal cases give rise to different frequency profiles, thus making the two cases potentially distinguishable by future experiments.

Phenomenological disformal enhancement
In section 3.2.2, we discussed two choices of the disformal function: a constant D = D 0 and a field dependent one, D = D 0 ϕ 2 (see table 4). As we have seen, the field dependent coupling gives rise to a larger maximum value of the Lorentz factor γ compared to the constant coupling, thus producing a bigger effect on the modified expansion rate. This is reflected in the amplitude of the PGWs as well. In figure 17, we compare the PGW spectra for D = D 0 (dashed line) and D = D 0 ϕ (dotted line) for the initial conditions specified in table 4.
As we can see from figure 17, the phenomenological disformal couplings give rise to very characteristic peaks in the PGW spectrum. In particular, the spectrum rises with a slope pro- portional tof 2 0 at the beginning, subsequently changing slope tof 25 0 for the case with constant disformal factor, andf 20 0 for the case with field dependent disformal factor. The spectrum then drops asf −3 0 for both cases. Note that this behaviour, namely the difference in the slopes for the two cases, has important implications as it makes these two couplings in principle distinguishable by future GW experiments. Interestingly, the slopes in the D-brane case are very different (see below), thus again being distinguishable from a phenomenological case.
The shifting effect of the initial temperature condition is illustrated in figure 18. In this plot, we use the same initial conditions as mentioned in table 4, but with an initial temperature given byT i = 10 11 GeV, instead ofT i = 10 7 GeV. The slopes of the curves are preserved, but the peak occurs at larger frequencies, which are relevant for the Einstein telescope [35]. Note that, for earlier initial temperatures, the peak can access ultra high frequencies, accessible by future experiments [36].

D-brane disformal enhancement
In the purely disformal D-brane-like case, C = 1 and thus D = 1/M 4 . As we have seen, the enhancement of the Hubble parameter is larger than in the phenomenological case (see figure 8).
This is thus reflected in the rise of the PGW spectrum, as can be seen in figure 17, where the solid line corresponds to the D-brane disformal case. Indeed in this case, the rise of the PGW spectrum is much larger than that in the constant and field dependent phenomenological cases. It could be possible that a power law disformal coupling, D = D 0 ϕ r with r > 2 for the phenomenological case might reach to and above the D-brane case. However, interestingly, the frequency dependence of the PGW enhancement is very different  in both cases. The D-brane disformal coupling produces a very characteristic enhancement with a slope proportional tof 2 0 at the beginning, subsequently changing tof 5 0 , in contrast to thef 25 0 orf 20 0 behaviour in the phenomenological cases. Moreover, the spectrum then drops as f −3 0 , i.e. with the same slope as in the phenomenological case. This behaviour can be directly understood from the behaviour of the Lorentz factor γ for the three cases. We can see from figure 8 that the rise in γ differs between the phenomenological cases (dashed and dotted lines), and the D-brane case (solid line), while it drops down with the same slope in all the cases. In figure 18, we show the rise in the PGWs when the initial condition for temperature is set at a higher value. Again, the D-brane case gives the largest enhancement, thus crossing the sensitivity curves of Einstein Telescope. Setting even higher initial temperatures will be relevant for the ultra high frequency experiments. Thus, a detection of this characteristic peaked spectrum by either ET, LISA or ultra high frequency experiments will tell us the epoch of scalar-tensor domination in the early universe.

Einstein frame rise and spinning axions
We finish this section with a discussion on the peaked enhancement of the PGW spectrum, when computed using the Einstein frame Hubble parameter (3.28). In this case, the spectrum grows and falls with a frequency slope that goes as Ω 0 and ω the Einstein frame equation of state of the total energy density ρ total . The peak arises since the energy density of the universe in the Einstein frame (see (3.17) with C = 1), behaves as matter with ω = 0 (ρ ∝ a −3 ) at early times, subsequently changing to a period of 'kination' domination due to the scalar field energy density with ω = 1 (ρ ∝ a −6 ), to finally start behaving as radiation when γ drops to one. This non-standard evolution of the energy density in the Einstein frame causes the peaked enhancement in the PGW spectrum. In figure 19 (left panel), we show the energy density of the scalar field, ρ ϕ , and the background, ρ bg , in the Einstein frame (see (3.17) with C = 1), while in the right panel we show the peaked PGW spectrum. A similar behaviour, although in a very different set-up, has been recently considered in [37,38], where a short kination era is implemented due to a spinning axion in field space. Let us stress however that the systems are very different, and in the present case, it is the Jordan frame where entropy and energy density are conserved (see (3.19)), and thus in this frameρ ∝ a −4 before the onset of BBN, as it should.

The concept of broken power-law sensitivity curve
We introduce the concept of broken power-law sensitivity curve (BPLS) as a graphical tool to show the sensitivity of GW experiments to scenarios producing stochastic GW background (SGWB) spectra with broken power-law profiles, as the ones we met in section 4.2. The BPLS generalizes the methods used to obtain the nominal and the power-law (PLS) sensitivity curves.
Nominal sensitivity curves -see e.g. [39] for a general discussion -provide a visual understanding on whether a GW event, characterised by a given frequency and amplitude, can or can not be detected with sufficiently high signal-to-noise ratio (SNR) by a GW experiment.
Power-law sensitivity curves (PLS) [34] make visually manifest the fact that, by integrating over frequencies, we can exploit the broadband nature of a SGWB signal, by accumulating SNR and making a detection of the SGWB more feasible even when its amplitude does not fit within the nominal noise curves. Assuming that the SGWB is described by a power-law profile in frequency, scaling as f β with β being a constant slope, the PLS is made by the envelope of experimental limits which can be placed for each slope β, varying β over a given interval, and integrating the signal over frequencies. As clearly explained in [34], this allows one to gain orders of magnitude in sensitivity, at the price of making hypothesis on the slope of the SGWB frequency profile. Moreover, additional sensitivity to a SGWB is accumulated integrating the signal over time.
However, the frequency profiles of various examples of SGWBs encountered in the literature are often not well described by a single power-law -they can instead be described by a broken power-law (see e.g. the analysis in [40,41]). A broken-power law profile changes slope at a value of GW frequency within the band of a given GW experiment. In this case, by using the standard PLS of [34] one tends to overestimate the sensitivity of an experiment to the SGWB signal. In fact, when assuming that the frequency profile is a single powerlaw over the entire instrumental band, one misses the possibility that the frequency profile changes slope and drops within the experiment frequency range. In this way, the signal 'leaves' the sensitivity region of the experiment sooner than what is expected for a single power-law. For this reason, we propose to generalise the same method of [34] to the case of broken power-laws, and obtain more realistic sensitivity curves adapted for this class of -29 -JCAP08(2022)010 models. We follow the discussion of [34] step by step, adapting it to our context. We start from the noise spectral densities of specific GW experiments. We assume that the signal has a broken power-law functional form, denoted by Ω β 1 ,β 2 (f /f ). The broken power-law is characterised by two slopes parameterised by the indices β 1 and β 2 , as well as a reference frequency f at which the slope of the spectrum changes. For any given choice of the set S = (β 1 , β 2 , f ) we determine the overall amplitudeΩ S of the spectrum which is necessary for achieving an SNR=10 by integrating over the entire instrumental frequency band. Then, the broken power-law sensitivity curve is the envelope of the profiles Ω β 1 ,β 2 (f /f ), with am-plitudeΩ S , drawn as a function of the frequency f . This method for building the BPLS requires an extremisation process over the three-dimensional space of the set S = (β 1 , β 2 , f ) (while the PLS requires extremisation only over a single parameter). We wrote a parallelised Fortran code for performing the procedure quickly and efficiently, which is available upon request. We plot the resulting BPLS curve in figure 20 for the case of four GW experiments: Einstein Telescope (we use ET-B version presented in [42,43] as nominal curve), LISA (we use the noise curves of [44]), BBO, and DECIGO (we use [45]). Figure 20 shows that the BPLS curves lose sensitivity with respect to the PLS ones, by a factor of order one. However, the loss in sensitivity is quite dependent on the interval over which we allow the parameters β 1,2 to vary. The wider the interval of variation for β 1,2 , the larger is the loss in sensitivity with respect to the PLS curves, because the broken power-laws can be more peaked, and we are left with a smaller frequency band wherein to integrate and accumulate SNR. We have chosen β 1,2 to vary over the interval (−25, 25) because it is the range of slopes we met in the scenarios of section 4.2.
Notice that PLS and BPLS in general tend to converge at the extrema of the frequency band of each experiment. This is understandable, since if the change in slope of BPL occurs at those extrema, then the signal effectively acts as a single power law over the rest of the frequency band, with the slight change in slope close to the extrema contributing very little to the SNR.
As with the PLS, the BPLS offers a practical visual aid for understanding whether the predictions of a given model can be detected by a GW experiment. If the model profile enters within the BPLS, then further sophisticated tools are needed for understanding whether the fine details of the model profile can be reconstructed -see e.g. the proposals and analyses in [44,46,47].

Conclusions
Cosmological inflation predicts the existence of a primordial gravitational wave (PGW) background spanning a large range of frequencies with an almost scale-invariant profile, whose amplitude is too small to be detected by future gravitational wave experiments. The fractional energy density in primordial gravitational waves, as observed today, depends on the primordial tensor spectrum generated during inflation, and crucially, on the expansion rate of the Universe from the end of inflation until today.
In this work we showed that if a period of well motivated scalar-tensor theories dominated epoch occurs before the onset of BBN, a large enhancement of the PGW signal can occur at frequencies probed by future GW interferometers such as LISA, ET, BBO and DECIGO.
We considered the most general physically consistent relation between two metrics and a single scalar field, which includes a conformal as well as a disformal (or derivative) cou--30 -JCAP08(2022)010 pling [29],g µν = C(ϕ) g µν +D(ϕ) ∂ µ ϕ ∂ ν ϕ. These couplings arise naturally from (D-)branes in (D-)brane models of cosmology [48]. In this case, the functions C and D are not independent (M 4 CD = 1) and thus, they cannot simply be turned on and off as in a phenomenological treatment of these theories.
We first considered the effect of a conformal coupling (conformal-disformal in the Dbrane case) on the PGW spectrum. In the phenomenological case, a conformal coupling was first considered in [21]. In the present work we explored different conformal functions (3.40) motivated by previous work in [6,19,20], and different initial conditions, to illustrate the conformal effect more clearly. Generically, the enhancement of the PGW has a step-like behaviour, similar to that found in [21] for the conformal function, C = e β φ 2 , and depends on the choice of conformal factor and initial conditions (see figures 3, 6, 10, 13). In the phenomenological case, we saw that a suitable conformal factor can enhance the flat spectrum, reaching the sensitivity curves of LISA and ET, and thus can be tested by performing correlations between these experiments. Interestingly, the conformal-disformal effect in D-brane-like scalar tensor theories is less effective at enhancing the PGW spectrum (see figure 15). This may be due to the disformal effect, which although subdominant could affect the enhancing effect.
We next considered the more interesting purely disformal effect, again both in a phenomenological set-up and in the non-trivial D-brane-like scalar tensor set-up. In the phenomenological approach, one can simply switch off the conformal factor, C = 0, while keeping D = 0. The pure disformal factor gives rise to a distinct peaked spectrum. We showed that the rising slope of this peaked spectrum behaves as (f 2 0 →f 25 0 ) for the case with constant disformal factor, while it behaves as (f 2 0 →f 20 0 ) for the case with field dependent disformal factor. However, for both cases, the spectrum falls asf −3 0 . The amplitude can be enhanced by a suitable choice of disformal function (see figure 17).
On the other hand, the D-brane-like pure disformal coupling, obtained by setting C = 1 and thus D = M −4 , gives rise to a substantial enhancement, larger than the phenomenological case, and more remarkably with a characteristic peaked spectrum with very different rising slopes (f 2 0 →f 5 0 ), but same falling slopes (f −3 0 )! In this case, the pure disformal rise reaches well within the LISA and ET sensitivity curves, depending on the initial conditions (see figures 17, 18). Indeed, as discussed in [20], the enhancement of the expansion rate, and thus the PGW, can be shifted to earlier/larger times/frequencies, by changing the initial conditions. This peaked-like spectrum with distinct slopes in the frequency is a feature of the disformal effect, which does not arise in other modified cosmologies and can thus be a smoking gun signal of a period of D-brane disformal scalar-tensor dominated epoch. In order to study the sensitivity of GW experiments to spectra with these features more accurately, we used the tool of broken power-law sensitivity curves, introduced in section 4.3.
Throughout this paper we have assumed λ ∼ 0, namely, the contribution of the scalar potential to the energy density being negligible during the whole evolution. In particular, we assume that the present day acceleration is due to a pure cosmological constant Λ. This does not need to be the case and the scalar field could contribute to the dark energy density, as quintessence. We leave the exploration of this possibility for future work. Another interesting aspect for future exploration is the possible implications of the scalar field displacement in the disformal case. As we saw, the disformal coupling drives the field quickly to large values, which may be imply the presence of additional fields, according to recent quantum gravity constraints (see [49][50][51] for recent reviews). with The two other solutions to eq. (A.13) are given by 2 ∆ Since we set the initial conditions on the field during the radiation dominated era, we know thatρ(T ) = π 2 /30 g * (T )T 4 . Therefore, for a particular set of initial conditions ϕ i , ϕ iÑ , the value of M is bounded by the following interval: (A.20) Using the above condition, for a particular set of values of ϕ i , ϕ iÑ , M , we can obtain the real positive solutions for eq. (A.13). This gives us the initial condition H i for the Hubble parameter. These initial conditions can be used together to obtain the initial value of the Lorentz factor, γ i . We find that, in order to satisfy the constraints on the present value of the Hubble parameter, we need to choose suitable initial conditions such that γ i ∼ 1.

A.4 D-brane conformal-disformal
The cubic equation for H 2 that needs to be solved in this case is give by [19,20]:

JCAP08(2022)010
One of the solutions to eq. (A.21) is given by The two other solutions to eq. (A.21) are given by 2 ∆ Since we set the initial conditions on the field during the radiation dominated era, we know thatρ(T ) = π 2 /30 g * (T )T 4 . Therefore, for a particular set of initial conditions ϕ i , ϕ i N , the value of M is bounded by the following interval: Using the above condition, for a particular set of values of ϕ i , ϕ i N , M , and by using the normalization condition M 4 CD = 1, we can obtain the real positive solutions for eq. (A.21). Note that, in order to solve the equations (3.48) in terms ofÑ , we have to convert the initial conditions obtained by the above procedure into initial conditions on the set of values of ϕ i , ϕ iÑ , M .
For the scenario wherein the field evolution is started with zero initial velocity, we find that A 1 = 0. In this case, the cubic equation for H 2 reduces to the following quadratic equation: This equation also illustrates that the choice of initial condition for H is independent of M and depends only on the initial temperature throughρ i and the initial field value through C i .