Signature of Non-Minimal Scalar-Gravity Coupling with an Early Matter Domination on the Power Spectrum of Gravitational Waves

The signal strength of primordial gravitational waves experiencing an epoch of early scalar domination is reduced with respect to radiation domination. In this paper, we demonstrate that the specific pattern of this reduction is sensitive to the coupling between the dominant field and gravity. When this coupling is zero, the impact of early matter domination on gravitational waves is solely attributed to the alteration of the Hubble parameter and the scale factor. In the presence of non-zero couplings, on the other hand, the evolution of primordial gravitational waves is directly affected as well, resulting in a distinct step-like feature in the power spectrum of the gravitational wave as a function of frequency. This feature serves as a smoking gun signature of this model. In this paper, we provide an analytical expression of the power spectrum that illustrates the dependence of power spectrum on model parameters and initial conditions. Furthermore, we provide analytical relations that specify the frequency interval in which the step occurs. We compare the analytical estimates with numerical analysis and show they match well.


I. INTRODUCTION
Since the first detection of gravitational waves (GWs) at the Laser Interferometer Gravitational-Wave Observatory (LIGO) and Virgo, a new window to unravel the mysteries of the cosmos has opened [1][2][3][4].Even though the sources of the detected GWs have been astrophysical thus far [5], we hope to also detect the cosmological ones with the advance of the detectors [6].Among the possible sources of GWs, the stochastic gravitational wave background (SGWB) originating from inflation is of particular interest, because detecting it may shed light on the history of the early universe [7][8][9][10][11][12][13][14][15][16][17][18][19][20].The standard model of cosmology, under the assumption of radiation domination (RD) from inflation to matter-radiation equality, predicts an almost scale invariant power spectrum for the SGWB across all frequencies [21][22][23][24].Yet, many well-motivated proposals predict a transient phase, dominated by a different energy component, intervening between inflation and the Big Bang Nucleosynthesis (BBN) [25,26].Depending on the temperature range and other specific features of this transient period, the SGWB profile is expected to vary.Inspired by numerous extensions to the standard model of particle physics, we focus our attention on early matter domination [27][28][29][30][31][32][33][34][35][36][37][38][39].Specifically, we assume a scalar field, ϕ, which behaves like a pressureless fluid scaling as a −3 with a being the scale factor, causes a transient matter domination era.Once ϕ decays, the universe reverts to the RD era again.Assuming a minimal coupling between ϕ and gravity, SGWB is only affected indirectly and because of the alteration of the scale factor and the Hubble rate [25,26].We coin the term Indirect effect to denote this influence.In the context of early matter domination, the SGWB power spectrum is suppressed at high frequencies compared to standard cosmology, and the Equation of State (EoS) of the Universe governs the slope of the power spectrum [26].
In this paper, we study early matter domination where the dominating field has a non-minimal coupling with gravity; i.e., f (R, ϕ) = 1 8πG − ξϕ2 R. In the term ξRϕ 2 , ξ signifies the gravitational coupling constant and R denotes the Ricci scalar.This coupling stands as the sole feasible local, scalar interaction of its kind with the appropriate dimensions [40].A coupling between a scalar field and gravity is proposed in numerous models aiming to resolve some of the inherent problems with inflation, reheating, and baryogenesis [41][42][43][44][45][46][47][48].For instance, in the context of non-oscillatory inflationary models [49][50][51][52][53], the term ξRϕ 2 has been recently used in an efficient reheating scenario denoted as Ricci reheating [54][55][56][57].In the same context, a novel quintessential Affleck-Dine (AD) baryogenesis scenario has also been proposed [58] based on the term ξRϕ 2 , avoiding troublesome iso-curvature modes found in conventional AD scenarios.Moreover, the widely used self-interaction term for scalar fields, λϕ 4 , necessitates the inclusion of ξRϕ 2 in the Lagrangian for proper renormalization in curved space-time [40,59].Within the scalar sector of the Standard Model (SM) Lagrangian (Higgs) [54], ξRϕ 2 emerges as the missing term that upholds all the symmetries of both gravity and the SM.
As ξ is subject to running, it cannot be universally set to zero across all energy scales [60].Two values of ξ = 0 for minimal coupling and ξ = 1/6 for conformal coupling are of particular interest.The latter, for the case m ϕ = 0, leads to conformal invariance of the action and hence the field equation of motion [40,61].Nonetheless, no compelling reason supports the presence of minimal or conformal coupling in the real world, as no symmetry is enhanced [46].Furthermore, since ξ is dimensionless, there is no reason for it to be small.It could be non-minimal, i.e. of the order of unity or more.1 [46].Current experiments place a weak constraint on ξ (ξ < 10 15 ) due to the feeble gravitational interaction with the SM fields [60,62] However, GWs might offer insights into the existence and strength of such a coupling.
This paper delves into SGWB within a cosmological framework characterized by early matter domination, where the dominant field ϕ directly couples with gravity (ξRϕ 2 ).D'Eramo et al. previously examined the case of ξ = 0 [26].To highlight distinctions from [26], we specifically explore the non-minimal regime (ξ ≥ 1).Studies indicate that such gravity couplings manifest as additional terms in the evolution equation of GWs [41,[63][64][65][66].That is, alongside modifications to the scale factor and Hubble, i.e., indirect effect [67], an extra factor directly affects GW evolution.We term this extra impact the direct effect.Employing Wentzel-Kramers-Brillouin (WKB) analysis, we elucidate the power spectrum resulting from GW propagation in our cosmological setting.Our numerical findings demonstrate that the presence of the ξRϕ 2 term, with non-minimal coupling, deepens the kink shape in the power spectrum, compared with the case studied in [26].Furthermore, owing to the direct effect, an additional step-like feature emerges, resulting in an enhanced reduction in the power spectrum.We highlight several benchmarks that could be probed by upcoming gravitational wave experiments.
To deepen our understanding of the ξRϕ 2 term, it is crucial to provide an analytical interpretation of the resulting power spectrum.Hence, we use reasonable approximations to get an analytical expression for the power spectrum compared with the case of the standard cosmology and explain how each part of the spectrum in a specific frequency interval is shaped due to a specific physical effect that was dominant in the corresponding temperature interval.To this end, the high-frequency modes that became sub at high temperatures and thus got affected by the changes in the cosmological evolution are of particular interest.In particular, we determine the fraction of the spectrum of high frequencies to that of low-frequency modes that did not experience any changes in the cosmological history.In this context, we utilize a physical quantity called the dilution factor, which was introduced in [26].While this analytical interpretation was performed for the case of ξ = 0 in [26], exhibiting notable agreement with numerical results, our study expands this analysis to encompass the case of ξ ̸ = 0.The analytical formula we derive, expressed in terms of dilution and damping factors, closely aligns with numerical results, maintaining an acceptable level of accuracy.
Finally, we comment on the potential testability of this scenario, highlighting a selection of gravitational wave experiments that could probe these changes [68].These include the Laser Interferometer Space Antenna (LISA), which operates within the frequency range of 10 −5 − 1 Hz [69][70][71], the DECi-hertz Interferometer Gravitational wave Observatory (DECIGO) spanning 10 −3 − 10 Hz [72][73][74], the Einstein Telescope (ET) operating within 1 − 10 4 Hz [75][76][77], and the forthcoming Big Bang Observatory (BBO) encompassing 10 −3 −10 Hz [78][79][80].Other experiments will explore intermediate frequencies like square kilometer array or SKA, probing 10 −9 − 10 −6 Hz [81,82] and NANOGrav collaboration [83,84] that works based on pulsar timing arrays (PTA) measurements [64,68]. 2  This paper is organized as follows: In Sec.II, we introduce the model and the cosmological framework in detail.In Sec.III, we study the evolution of GWs, highlighting the effect of ξRϕ 2 .We present the WKB solution for the modes that become sub, in Sec.III A. We follow in Sec.III B, by introducing the power spectrum as the observable of GWs, and in Sec.III C, we highlight two of the pivotal frequencies in the power spectrum.Sec.IV is dedicated to numerical results for the power spectrum and a comparison with the analytical estimates.Finally, the concluding remarks are presented in Sec.V.

II. NON-MINIMALLY COUPLED SCALAR FIELD
Our theory is defined by a real scalar field ϕ that has feeble interactions with other fields, but has a non-minimal coupling with gravity: f (R, ϕ) = 1 8πG − ξϕ 2 R. The action of our theory is, thus, [56,60]: where g is the determinant of the metric g µν , G is the gravitational coupling constant, m ϕ is the mass of ϕ, and L M contains the kinetic component and the interactions of any other field in the cosmos, including its interaction with the scalar field, ϕ.Varying Eq. ( 1) with respect to g µν , one can obtain the gravitational equation as [60,63]: In the definition of Eq. ( 2), the energy-momentum tensor of the scalar filed ϕ is: where we have □ϕ ). Varying Eq. ( 1) with respect to ϕ yields the scalar field equation of motion [63]: Using Eqs.(2,4) and the Bianchi identity, = 0, the continuity equation of the matter part can be obtained, . In these equations, the field ϕ is only a function of time in order to respect the homogeneity and isotropy of the Universe, and L M describes the radiation.Therefore, the evolution equations for the scalar field ϕ and the energy density of radiation ρ R are: ρR + 4 where to obtain these equations, we have employed the flat Friedmann-Lemaitre-Robertson-Walker (FLRW) metric, ds 2 = −dt 2 +a(t) 2 δ ij dx i dx j ,4 with i, j = 1−3 specifying the spatial coordinates, explains the metric of a homogeneous and isotropic Universe, a is the cosmic scale factor, and accordingly the Ricci scalar is obtained as R = 6( ȧ2 + aä)/a 2 [63].Furthermore, Γ represents the total decay rate of the scalar.In this paper, we are oblivious to the exact nature of the fields into which ϕ decays, and rather we assume it is part of radiation.The total number of relativistic degrees of freedom contributing to the energy (entropy) density of radiation is represented by g ρ ⋆ (g s ⋆ ).The evolution of the Hubble rate, defined as H ≡ ȧ a , is described by the first Friedmann equation, obtained from Eq. ( 2) as: where M P l is the reduced Planck mass, and the energy density of ϕ can be obtained from the T (ϕ) 00 component of Eq. (3) [56,60]: Using Eqs. ( 5), (6), and (7), the system of coupled differential equations for the background would be complete for three unknowns, ϕ, ρ R , and H.The free parameters in these equations are In this project, we fix T in = 10 11 GeV, motivated by well-studied reheating scenarios which prefer T Reh ∈ (10 12 − 10 15 ) GeV [54,55,[89][90][91][92].The initial value of ϕ in has an upper bound of M P l to avoid undergoing the Universe to another inflationary period [26].In addition, since the most noticeable deviations from standard cosmology occur when ϕ in ∼ M P l , we choose ϕ in = 10 18 .The value of ξ, also, has an upper bound to make sure the Hubble rate stays real, i.e., ξ < (M 2 pl /ϕ 2 in ), which in our case gives rise to ξ max ≃ 5.95. 5 Since our setup is not sensitive to φin and any initial condition satisfying leads to the total equation of state ω ≃ 1/3; and hence, a negligible φ (see Appendix A), we fix φin = 0 for simplicity.Concentrating on one of the benchmarks of Ref. [26], i.e. m ϕ = 10GeV and Γ = 10 −8 GeV, we numerically solve the set of coupled equations and compare the results for the two cases of ξ = 0 (minimal coupling) [26] and ξ ≥ 1 (non-minimal coupling).6

A. General Evolution of the Background
Given our initial conditions at high temperatures, the friction term in Eq. 5 ( more specifically the Hubble term, H φ) forces ϕ to be stuck at its initial value, similar to the case of ξ = 0 [26].As the Hubble rate decreases, eventually the friction term becomes less efficient and ϕ starts oscillating when H ≃ m ef f ≡ ξR + m 2 ϕ .As demonstrated in Fig. 1, the effect of gravitational coupling causes ϕ to remain close to its initial value for a longer duration, because the Hubble rate is higher and m ef f , because of negative R, is lower for larger ξ: ) −1 (see Appendix A for this simplified expression of Hubble rate at high temperatures).Around the time when ϕ starts oscillating, the universe deviates from state of RD, or equivalently ω ≃ 1/3, or equivalently the state of Radiation Domination (RD).The temperatures at which this occurs is when H and m ef f lines meet (H ≃ m ef f ).The numerically determined factor i osc ≃ 6 has been implemented to obtain a more precise value for T osc , since ϕ oscillations start slightly after H ≃ m ef f (e.g., for ξ = 0, oscillations start at T ∼ 10 9 and for ξ = 5.95, they start at T ∼ 10 8 ).After some oscillations,7 ϕ behaves as sin(m ϕ t) m ϕ t .Hence, the energy density of ϕ redshifts like cold matter, leading to an era of matter domination (MD era).This process continues until the decay of ϕ becomes more efficient than the Hubble rate, at which point ϕ depletes.We refer to this era, the decay era (DE).When ρ ϕ becomes negligible with respect to ρ R , we return to radiation domination, which we will refer to as the Late RD era.The temperature at which we return to late RD era can be found by setting ρ ϕ ∼ ρ R [26]: when φ ≃ m ϕ ϕ and ρ ϕ ≃ φ2 , we replace Eqs. ( 5) and ( 6) with the following approximate Boltzmann Equations: Since g s ⋆ and g ρ ⋆ vary with temperature, especially at T EW ≃ 100GeV ≤ T ≤ T BBN , we evaluate them as a function of ρ R using Ref. [93].where α RD ≈ 0.64. 8Assuming a constant ω throughout the MD era, c RD can be numerically determined as c RD ≃ 1.07.At this point, the evolution of Universe is followed by the standard model of cosmology.
To investigate the differences between minimal (ξ = 0) and non-minimal (ξ ≥ 1) ϕ−gravity coupling, GWs are the best candidates.The imprint of ξ = 0 case on GWs has been investigated in Ref. [26].Due to the non-minimal coupling between ϕ and gravity, the evolution of ϕ significantly impacts in the evolution of GWs, which will be discussed in details in the following section.

III. THE EVOLUTION OF GWS AND THE POWER SPECTRUM
Gravitational waves, h ij , are perturbations in the space part of the FLRW metric [13,15,16,94]: with |h ij | ≪ 1. Taking h ij to be transverse and traceless h T T ii = k i h T T ij = 0 , there are two degrees of freedom for h ij indicating two polarizations, ζ = +, ×.By linearizing the gravitational equation, Eq. ( 2), the evolution equations of GWs is obtained: where It can be seen that an additional term, A, appears in the well-known evolution equations of the GWs [13][14][15][16], which comes from the ξRϕ 2 term.It is through this term that the evolution of ϕ directly impacts the evolution of the SGWB.
It is conventional to expand ij is the polarization tensor [13,16,94]. 9Substituting the polarization expansion and h 14), the evolution equation for the 8 Eq. 12 is derived assuming a constant EoS, or equivalently assuming a ∝ t α , with α being a constant. 9In our convention, e i e (1) i e (2) ) and e i e (2) i e (1) j ) [16].

FIG. 2:
This plot illustrates the rates of 2aH and aA, which represent the damping factors in the evolution of GWs, for the case of ξ = 5.95.At high temperatures, A ∝ φ is zero.However, as ϕ gains momentum, A rapidly grows and surpasses H within certain temperature intervals.In addition, the ripples of 2aH and aA become comparable at the onset of the oscillation era, and thus cancel each other out.
amplitude of GWs, h k is obtained as following [63]: where h k (t) represents the Fourier transform of h(x, t).For a better understanding of the behavior of GWs, Eq. ( 15) is expressed in terms of conformal time,η, given by η = dt a , as in Refs.[64][65][66]: with the prime denoting differentiation with respect to conformal time.In what follows we do not mention the subscript k, since the term k 2 can highlight that the equation is written for a specific momentum.The above equation is a harmonic oscillator with the damping term a 2 (2H + A).Using the numerical results of the background as obtained in Sec.II, 11 we distinguish different regimes in the behavior of the damping term, by comparing the evolution of the terms aA and 2aH as depicted in Fig. 2: • A ≳ 2H : Even though A is initially zero (A ∝ φ and φin = 0), it rapidly grows as ϕ gains momentum.
Depending on the value of ξ and our initial conditions, A may eventually exceed the Hubble rate.We refer to this era as A Domination era (AD era) for short. 12During AD era, the evolution of GWs is not only affected by the altered background but also directly by A. To highlight the effect of A, we refer to its contribution as the direct effect.Note that the direct effect goes to zero as ξ → 0. Conversely, for large values of ξ A remains comparable to H (A ∼ H) during the first few oscillations, and strongly damps the oscillations.
• 2H ≫ A ̸ = 0 : The effect of A rapidly diminishes, and the Hubble rate dominates the the damping term.In this case, the evolution of GWs is mainly influenced by the modified background.Since it is the effect of ξRϕ 2 coupling that indirectly affects the evolution of GWs through a and H, we denote this effect as indirect effect.
Eq. ( 15) is the same for both polarizations.Therefore, the index ζ is omitted in what follows.
Due to the fact that GWs are just perturbations on the background and do not have any back-reaction on it, it is justified to solve the evolution equations of background separately and then the resulting a, H, and A are inserted in the equation of GWs.In addition, this approach has two technical advantages: the background is solved only once, which significantly reduces the running time of the code, and fixes the precision of the numerical solution of the background for all modes.
Even though this Period effects the evolution of GWs, it does not change the EoS of the background.At the background level, the EoS parameter is still around 1  3 which denotes that we are still in RD era.
• A → 0 : This is the case in the late RD era, when ϕ is almost gone.Therefore, not only A diminishes but also the evolution of background becomes that of the standard cosmology.
Different regimes leave different signatures on the GWs power spectrum, as we will show in Sec.III B. For this purpose, let us first describe the behaviour of GWs at different frequency regimes.

A. Evolution of GWs
Recalling the well-known theory of damping harmonic oscillators, the general behavior of a specific mode in Eq. ( 16) can be understood by comparing its momentum, k, with the damping term, a 2 (2H + A).In this regard, we categorize different modes as: • Sup Modes: For a mode with k ≪ a 2 (2H + A), there exists a constant solution to Eq.( 16), h Prim k , where the superscript "Prim" stands for primordial.This solution describes the frozen behavior of the mode from the time it leaves the horizon during the inflationary epoch until it reaches the time η × when k ∼ a× 2 (2H × + A × ).The solution h Prim k serves as an initial value for the subsequent evolution of the mode after k ≫ a 2 (2H + A). • Sub Modes: For a mode with k ≫ a 2 (2H + A), it is conventional to parametrize the solution of Eq. ( 16) as , where χ k (η) represents the transfer function that describes the subsequent evolution of the mode from η × to the present time.By substituting this parametrization into Eq.(16), the evolution equation for the transfer function becomes: which has a general solution presented in Appendix B. By imposing the initial conditions χ = 1 and χ ′ = 0 for k ≪ a 2 (2H + A), the specific oscillatory solution for Eq. ( 17) can be expressed as: where with d k (η) denoting the damping factor for a specific mode k computed up to η.To calculate the observable, power spectrum of GWs, we set η = η 0 to obtain the WKB solution at present, where a 0 = 1.The modes that become sub prior to the AD era, experience the whole AD era, and thus their amplitudes behave as |χ k | = e −d a × .Conversely, the k modes that become sub after the AD era when A is negligible, do not experience any suppression from the d factor and they go as |χ k | = a × .If a k mode becomes sub during the AD era, the suppression factor on its amplitude will depend on its k value: Finally, let us comment on the behavior of a × as a function of k to gain an intuition about |χ k | in different k regimes.Since for the cosmological eras with constant ω, a ∝ η 2 1+3w , we can obtain a × given the functional form of 1+3w .However, this argument is only valid when the effect of A is ignored (before and after AD era when H ≫ A).For the modes that enter the horizon during AD, this simple analysis does not suffice and numerical analysis is required.

B. The Power Spectrum
To introduce the power spectrum of gravitational waves as an observable, we need their energy density given by [14,93]

T osc T osc T RD
FIG. 3: This plot demonstrates the evolution of a 2 (2H + A) with respect to temperature.The red line is for the case of ξ = 0, the blue line is for ξ = 5.95 while ignoring the A term, and the green line presents the case of ξ = 5.95 keeping the A term.Notice that since the A term is zero (ignored) in the red (blue) line, we have a 2 (2H + A) → aH.In this plot, T osc s are shown by stars, and T RD is represented by a triangle.As the definition of T osc does not depend on A, the temperature of oscillation for the blue and green lines coincide.Furthermore, T RD is independent of the value of ξ, and thus all of the lines in the plot share the same T RD .
where ⟨...⟩ av denotes spatial averaging.The power spectrum, Ω GW (t, k), is defined as the energy density of GWs per logarithmic frequency interval, divided by the critical energy density, ρ crit = 3H 2 (8πG) −1 , of the Universe [14,26,93]: The quantity P T (k) is the primordial tensor power spectrum which is a nearly scale-invariant power-law function of k around a pivot scale k * as P T (k) = A T (k/k * ) n T .Here, A T is the amplitude and n T is the tilt of the spectrum [16,26,94].Using the WKB solution, Eq. ( 18) for the sub modes k ≫ a 2 (2H + A), we get |χ ′ (η, k)| ≃ k|χ(η, k)| (see Appendix B).By substituting this derivative in Eq. ( 21) and setting η = η 0 , the present power spectrum for GWs is obtained as where f = k/(2πa 0 ), χ k , a 0 , and H 0 = 100h km/s/M pc are the present time quantities.In this work, we set A T ≃ 1.5 × 10 −10 , n T = 0.4 and k ⋆ = k CMB = 0.05 1 M pc in P T (k) as is done in [26] to access the maximal reach of future detectors.These values are consistent with the upper bound on the tensor-to-scalar ratio r = A T A R ≤ 0.07, where A R ≃ 2.1 × 10 −9 is the amplitude of scalar perturbations obtained from observational data [95]. 13Eq.( 22) elegantly shows how the parametrization h k (η) = h Prim k χ k (η) shows up in the power spectrum.Therefore, Ω 0 GW is generally made by concerning two factors: 1) the initial condition h P rim k and 2) the value of |χ k | 2 (influenced by the indirect and direct effects).Note that GWs start from after becoming sub, coded in the primordial tensor power spectrum, P T (k).The form of power spectrum Ω 0 GW for the high-frequency modes (becoming sub before AD) and low-frequency The power spectrum of GWs for standard cosmology (dashed black line), ξ = 0 (red line), ξ = 5.95 ignoring the A term (blue line), and ξ = 5.95 preserving the A term (green line) with respect to frequency is presented.The frequencies corresponding to the beginning of ϕ oscillation (Eq.( 26)), late RD era (Eq.( 25)), and the step (Eq.( 29)) are shown by star, triangle, and square, respectively.The ripples of the blue line are due to the intense oscillation of ϕ which is felt by the modes that become sub during the first RD era.In the case of the green line, the presence of A damps the ripples.Moreover, the A term leads to a noticeable drop in the power spectrum of the green line at high frequencies, manifesting a step-like feature.

C. Frequencies fosc and fRD
Before diving into the numerical results on power spectrum, let us discuss the relation between a k mode, or more precisely its corresponding frequency f = k/(2πa 0 ), and the temperature at which it becomes sub.Since this happens when k ≃ a 2 (2H + A), we can study Fig. 3 to find the transition frequencies.In this figure, the red line represents the ξ = 0 case, the green line is for ξ = 5.95, and we have included the blue line where the case of ξ = 5.95 is redrawn neglecting the effect of A. It is important to emphasize the blue line is not a physical case; however, it can illustrate the role of A more lucidly.Other than this numerical method using Fig. 3, we can use the following theoretical expression to find the frequency at T = T RD [26]: where T 0 ≃ 2.72548K ≃ 2.3 × 10 −13 GeV is the CMB temperature [96], the frequency f osc corresponding to the mode that becomes sub at T osc can be obtained using Eqs.( 10) and ( 12): Having these frequencies we can determine our scenario how narrates on power spectrum Ω 0 GW (f ).

IV. RESULTS
To obtain the power spectrum Ω 0 GW , we need to evolve GWs till today.However, given that after ϕ decay, we return to the standard cosmology, we can use simple scaling arguments to determine Ω 0 GW from the power spectrum after ϕ's effect has been diminished.To do this, we solve Eq. ( 17) from T in = 10 11 GeV to T f = 100 GeV, then scale |χ k (T f )| to its present value |χ k | using the conservation of entropy, S(T ) ∝ g s * (T )a 3 T 3 [93].In this section, our main goal is to investigate how the non-minimal (ξ ≥ 1) gravitational coupling of the scalar field, ξRϕ 2 , affects the power spectrum.In particular, we want to distinguish the direct and indirect effects on Ω 0 GW .This splitting can be done since they show up independently in the WKB solution, Eq. (18).To proceed, we study the evolution of χ k with and without the A term for ξ = 5.95, and contrast them with the case of ξ = 0. 14Fig.4 demonstrates the values of Ω 0 GW for the case of standard cosmology (dashed black line), the case of ξ = 0 (red line), and the non-minimal gravity coupling case ξ = 5.95, with (green line) and without (blue line) the A term.The frequencies that become sub during the RD era exhibit an almost scale-invariant spectrum.Introducing a transient MD era, causes high-frequency modes to be lower relative to the case of standard cosmology.That is because during the MD era, the damping term in Eq. ( 17) becomes a ′ a = aH ∝ 2 η , which is greater than that of the RD era, a ′ a = aH ∝ 1 η .Modes that become sub during the MD era, experience a weaker damping.
The aforementioned effects on Ω 0 GW due to the MD era are generally the same for both ξ = 0 and ξ > 0. However, the presence of the ξRϕ 2 term leaves a distinct imprint on the shape of the power spectrum.To better comprehend these effects, let us compare the blue line (ξ = 5.95, while ignoring the A term) in Fig. 4 with the red line (ξ = 0).The main difference between these two lines is the change in the Hubble rate due to the ξRϕ 2 term.The main points to highlight from the comparison of the blue and red lines are: 1) the oscillations of H appear in the power spectrum; 2) the MD era becomes longer15 making a deeper kink in Ω 0 GW ; 3) the temperature at which the Late RD era begins, coincides with that of the case ξ = 0.
To grasp the effect of the A term on Ω 0 GW , we compare the green line (ξ = 5.95, while preserving the A term) with the blue line (ξ = 5.95, ignoring the A term) in Fig. 4. As can be seen from these two lines, the A term has a notable effect on the shape of the power spectrum.One effect is that the oscillations on Ω 0 GW vanish, since the oscillations of aA cancel those of 2aH.Furthermore, an additional step-like feature appears in Ω 0 GW , which can be understood by the extra damping term in Eq. ( 17) during the AD era.We can find the temperature at which this step occurs, T step by setting A ≃ 2H/β, where the parameter β depends on ξ.Using this along with the form of A in terms of ϕ and Eq.(A3), the exact expression for φ2 at T step is obtained as: where γ = . Substituting A obtained from Eq. ( 27) and , where α step ≃ 0.57.This frequency is shown by a square in Fig. 4. To complete this part of our study, we obtain Ω 0

GW
for different values of ξ, and present the results in Fig. 5.As can be seen, the kink in Ω 0 GW deepens by increasing ξ, as expected.
Another way to interpret our results is by using dilution factor : D ≡ [S(T ≫ T osc )/S(T ≪ T RD )] , along with T osc and T RD given by Eqs.(10) and (12), the dilution factor becomes with α D ≈ 0.68 obtained from numerical matching.To interpret the results of Fig. 5, we employ the WKB solution (see Appendix B).For the sub mode with the wave number k, we have |χ k | ≃ a k .The parameter a k is the scale factor at which the k mode becomes sub.Replacing |χ k | in Eq. ( 22) with a k = k/H(a k ), for the modes with frequencies f osc and f RD , we obtain Ω 0 GW (f osc )/Ω 0 GW (f RD ) (ignoring A term) as: which demonstrates the indirect effect of the term ξRϕ 2 .Second, if we set the wave number as k = k step , we have where d = 1 tstep Adt.Eq. ( 32) exhibits the direct effect of the term ξRϕ 2 .Finally, implementing the correction parameter cr as the ratio of h 2 Ω 0 GW (f osc ) in the power spectrum with A to that of the power spectrum without A, we can write where and GW (f RD ) are given by Eqs. ( 32) and (31), respectively.To test our estimations, we compare some analytical results with numerical ones.Table I compares the analytical estimates and data driven pivotal frequencies f RD,osc,step for different values of ξ.Furthermore, the ratio of Ω 0 GW at these key frequencies according to the analytical expectations (Eqs.(31), (32), and ( 33)), and the data are provided in Table I as well.As it is illustrated, the predicted values are consistent with the numerical data.Albeit the analytical expressions depend on data driven parameters α RD,D,step and β, they prove to be valuable as they show the dependence of our observable Ω 0 GW (f ) on model parameters and initial conditions.For the sake of completeness, let us study the effect of m ϕ and Γ on the shape of Ω 0 GW , in the case of ξ = 0.Here we choose m ϕ = 10 GeV, 40 GeV and Γ = 10 −8 GeV, 4 × 10 −8 GeV.Since m ϕ determines the beginning of the MD era and Γ controls when it ends, only the position of the kinks change and the slopes in Ω 0 GW stay constant, as shown in Fig. 6.Finally, the effects of non-minimal gravitational coupling, ξRϕ 2 , of the scalar field on Ω 0 GW can be probed by future GWs experiments.Due to the shape of Ω 0 GW , we expect GWs to be detectable by some of the observatories focusing on low-frequencies but missed by high-frequency experiments.The shape of the GWs may also be probed by some experiments that can detect relatively lower-intensity GWs.In Fig. 7, we have picked some benchmarks that may be probed by LISA [97], BBO [98], and DECIGO [99].

V. CONCLUSION
In this paper, we investigated the evolution of the SGWB originating from inflation within a cosmological framework.We considered a scenario where a scalar field ϕ dominates the energy density of the Universe at high temperatures and has a non-minimal coupling with gravity, represented as ξRϕ 2 .Previous studies have demonstrated that early matter domination leads to a reduced signal strength of SGWB in the present.Various experiments such as LISA, BBO, and DECIGO can probe SGWB, while others like HLVK and ET may miss it.
By introducing a non-minimal coupling between ϕ and gravity, not only does this decrease amplify, but non-trivial features also emerge at certain frequencies.When the coupling parameter ξ is large, the damping term in the evolution of SGWB is enhanced due to a larger Hubble rate and the presence of a ξ-dependent term denoted as A in this study.The distinctive step-like feature observed in Fig. 7 is caused by the dominance of A over the Hubble rate.Analytically, this feature can be understood by examining the ratio Ω 0 GW (fstep) Ω 0 GW (f RD ) , which includes not only the dilution factor but also an additional damping term e − Adt compared to the case where ξ = 0. Additionally, when ξ ̸ = 0, the Hubble rate is relatively higher at the beginning of the MD era, resulting in a lower value of Ω GW h 2 , even if the A term is neglected, as the dilution factor discloses.We provided an analytical expression for Ω 0 GW at key frequencies and demonstrated that they align with the numerical results with acceptable precision.
In the scenario of early matter domination with ξ = 0, the power spectrum pattern exhibits deviations from a pure power law at frequencies determined by the mass and decay width of ϕ.Our study reveals that introducing the term ξRϕ 2 does not significantly alter the location of these deviations but introduces a distinct step feature that serves as a signature of this term in observational data.We identified a set of benchmarks that exhibit this feature within the parameter space where DECIGO and BBO experiments can probe.In general, the shape of the power spectrum of Primordial Gravitational Waves (PGWs) as a function of frequency provides valuable insights into the evolution of the early universe.FIG.7: The sensitivity curves of existing and future GW experiments and the power spectrum of GWs in ξ = 0 and ξ ̸ = 0 scenarios.There is hope that the effect of direct and indirect effects of modified gravity ξRϕ 2 can be captured for specific regions in the parameter space of the model by future detectors.The curves are resulted using different benchmarks of (ξ, m ϕ , Γ).Since we only can analyze the modes that become sub after T in = 10 11 GeV, the colored solid lines that are the results of the paper do not extend to arbitrarily large frequencies.This appendix shows that the benchmarks satisfying Eq. ( 9) have an attractor behavior.That is, even if we start with φin ̸ = 0, we will quickly merge to the path of ϕ = ϕ in and φin = 0 in the phase space.Our other objective in this appendix is to show that even in the case of non-minimal coupling where ξ ≥ 1, the start of evolution is effectively RD defined as ω = 1  3 , even though the energy density of the scalar field dominates over that of radiation ρ ϕ > ρ R .For better understanding, we rewrite the EoM of ϕ, Eq.( 5), with the assumption H ≫ Γ at the beginning as This equation has the form of a harmonic oscillator equation, with m ef f as the angular frequency and the friction term of 3H 2 .In the limit where 3H 2 ≫ m ef f , and with the initial conditions ϕ(T in ) = ϕ in and φ(T in ) = 0 the solution for ϕ becomes ϕ = ϕ in = const, which means that the field is stuck at its initial value.If the initial velocity of the harmonic oscillator is non-zero, the solution to the equation in this regime becomes which shows that in this case, the general solution does not have constant behavior and in fact is time-dependent.For any value of φ satisfying Eq. ( 9), the second term in Eq. (A2) becomes negligible compared with ϕ in and hence we can write the solution as ϕ(t) = ϕ in − φ 3H e −3Ht .In addition, the third term rapidly approaches zero with time, and and invoking the useful relation h(t) dt = h(f (x))f ′ (x) − h(g(x))g ′ (x), (B10) we can simply conclude that By substituting the above equation in Eq. ( 21), Eq. ( 22) is obtained.

FIG. 1 :
FIG.1:(a)The behavior of the scalar field, ϕ, in terms of temperature, T .(b) The evolution of the Hubble parameter and the effective mass as a function of temperature.In both plots, the red line is for ξ = 0, and the blue line is for ξ = 5.95.As shown in these plots, the Hubble rate increases and m ef f decreases with the value of ξ, and therefore the temperature at which H ≃ m ef f is lower for larger values of ξ.Thereby, the oscillation of ϕ is delayed with non-minimal coupling of ϕ with gravity.

10 − 5
FIG. 5: The power spectrum for different values of scalar-gravity coupling ξ is shown.Increasing ξ results in a decrease in the power spectrum at high frequencies.Since both the Hubble rate and A are proportional to 1 8πG − ξϕ 2 −1 , increasing ξ such that ξϕ 2 ≃ 1 8πG significantly increases the Hubble rate and A. Therefore, the drop at high frequencies for the case of ξ = 5.95 (green line) is considerably bigger than ξ = 5 (magenta line) or ξ = 4 (blue line).
ported by the Cluster of Excellence Precision Physics, Fundamental Interactions and Structure of Matter (PRISMA + EXC 2118/1) within the German Excellence Strategy (project ID 39083149).Also, F.E. is supported by the grant 05H18UMCA1 of the German Federal Ministry for Education and Research (BMBF), and H.M. is supported by the National Natural Science Foundation of China under Grants No. 12247107 and 12075007.Appendix A: Beginning of the evolution: Attractor Behavior, Hubble parameter, Ricci scalar, and EoS parameter 1/3, where S = a 3 s with s = 2π 2 45 g s * T 3 being the comoving entropy density.Using a ∝ t α ≃ ( α H ) α , and H osc ≃

TABLE I :
(32) table presents the values of (1) α RD,D,step and β obtained from numerical matching, (2) f RD,osc,step according to theoretical expectation and according to numerical analysis, (3) the dilution factor (4) the damping factor e −2d , and (5) the ratio of Ω 0 GW for each of the mentioned frequencies according to both theory and data analysis, for different values of ξ.Albeit one needs some numerical input for the theoretical prediction, this table shows that Eqs.(31),(32), and (33) match the numerical results with satisfactory accuracy.