Implications of NANOGrav results and UV freeze-in in a fast-expanding Universe
IImplications of NANOGrav results and UV freeze-in in a fast-expanding Universe
Basabendu Barman, ∗ Amit Dutta Banik, † and Avik Paul ‡ Department of Physics, IIT Guwahati, Guwahati-781039, India Key Laboratory of Quark and Lepton Physics (MoE) and Institute of Particle Physics,Central China Normal University, Wuhan 430079, China Astroparticle Physics and Cosmology Division, Saha Institute of Nuclear Physics,HBNI, 1/AF Bidhannagar, Kolkata 700064, India
Recent pulsar timing data reported by the NANOGrav collaboration indicates the existence ofa stochastic gravitational wave (GW) background at a frequency f ∼ − Hz. We show that adark sector consisting of a Standard Model (SM) gauge singlet fermion χ and a singlet scalar φ ,both charged under a Z symmetry, is capable of generating such a low frequency GW via strongfirst order phase transition (SFOPT) through the modification of the standard cosmological history,where we assume faster-than-usual expansion at pre-BBN times driven by a new cosmological species ϕ whose energy density red-shifts with the scale factor as ρ ϕ ∝ a − (4+ n ) . Depending on the choice ofthe fast expansion parameters, reheat temperature and effective scale of the theory, it is also possibleto address correct dark matter (DM) relic abundance via freeze-in. We show that a successfulfirst order phase transition explaining NANOGrav results together with PLANCK observed DMabundance put bound on the fast expansion parameters requiring n (cid:46) I. INTRODUCTION
Recently, the North American Nanohertz Observa-tory for Gravitational Waves (NANOGrav) has reporteda possible detection of a Stochastic gravitational wavebackground (SGWB) with their 12.5 year data set ofthe pulsar timing array [1] that may be interpreted asa gravitational wave (GW) signal with amplitude of ∼ O (cid:0) − (cid:1) at a reference frequency of f yr = 1yr − .The NANOGrav collaboration has found no evidence ofmonopolar or dipolar correlations which may arise fromreference clock or solar-system ephemeris systematics,they also claim that there is no statistically significantevidence that this process has quadrupolar spatial cor-relations which can be a smoking gun for GW back-ground. However, several attempts were proposed, soas to interpret the NANOGrav observations which in-clude supermassive black hole merger events [1, 2], cosmicstrings [3–5], magneto-hydrodynamic (MHD) turbulenceat a first-order cosmological QCD phase transition [6] orphase transition from a hidden sector [7]. As shown in [8]the NANOGrav data can be well explained as a GW sig-nal from first order phase transition (FOPT) around thewarm dark matter (DM) physics scale that lies in thekeV-MeV range, well below the electroweak scale. In [9]it was shown that NANOGrav observation can be ex-plained by the first order GWs in the nonstandard ther-mal history with an early matter dominated (EMD) era.The expansion rate of the universe is parametrized bythe Hubble parameter which is directly related to thetotal energy density of the universe through standardFriedman equations [10]. In the standard case, it is as-sumed that the universe was radiation dominated before ∗ E-mail: [email protected] † E-mail: [email protected] ‡ E-mail: [email protected]
Big Bang Nucleosynthesis (BBN). However, there are nofundamental reasons to assume that the universe wasradiation-dominated prior to BBN at t ∼ freeze-in it has beenshown that in a fast expanding universe the DM produc-tion via decay and/or annihilation of the bath particlesis dramatically suppressed that demands a larger cou-pling between the dark and the visible sector to satisfythe observe relic abundance. This sheds some optimisticprospects on detecting DM produced via freeze-in whichotherwise is extremely difficult to detect.Motivated from all these in this work we propose asimple framework where we assume the presence of an-other species ϕ whose energy density red-shifts with thescale factor as ρ ϕ ∝ a − (4+ n ) . For n > ϕ en-ergy density dominates over radiation at early enoughtimes. However, the equality between ρ ϕ and ρ rad hap-pens at a temperature T R (cid:38) T BBN . This, as one canunderstand, alters the Hubble parameter leading to afaster-than-usual-expansion of the universe in the pre-BBN era. In such a modified cosmological backgroundwe consider a simple dark sector made up of a Stan-dard Model (SM) gauge singlet fermion χ and a singletscalar φ both charged under a Z symmetry. We assume χ is lighter than φ and hence it serves as a potentialDM candidate. The lack of any renormalizable interac-tion between the dark sector and the SM motivates theDM production via Ultra-violate (UV) freeze-in [13–15]. a r X i v : . [ a s t r o - ph . C O ] D ec The requirement of satisfying PLANCK observed [16, 17]relic abundance constraints the reheat temperature, ef-fective scale of the theory and also the fast expansionparameters n and T R . The singlet scalar, on the otherhand, triggers a strong first order electroweak phase tran-sition (SFOPT) that generates detectable GW signal.We show, in a fast expanding universe such GW sig-nals generated from electroweak phase transition canexplain the NANOGrav 12.5 yrs data within 2 σ confi-dence. The requirement of producing right DM abun-dance together with GW signal satisfying NANOGravresults put bound on the fast expansion parameters, typ-ically requiring n (cid:46) T R ∼ O (100) MeV.The paper is organized as follows: in Sec. II we describethe modification in the Hubble parameter due to fastexpansion; the effect of modified Hubble parameter onDM yield in the context of a toy model is detailed inSec. III; in Sec. IV we discuss the generation of detectableGW signal from strong first order phase transition that iscapable of explaining the NANOGrav 12.5 yrs data; wethen briefly touch upon a possible UV completion for thistoy model in Sec. V and finally we conclude in Sec. VI. II. A FAST EXPANDING UNIVERSE
As mentioned in Sec. I, we assume the universe beforeBBN has two different species: radiation and some otherspecies ϕ with energy densities ρ rad and ρ ϕ respectively.In presence of a new species ( ϕ ) along with the radiationfield, the total energy budget of the universe is ρ = ρ rad + ρ ϕ . One can always express the energy density of theradiation component as function of temperature T [10] ρ rad ( T ) = π g ∗ ( T ) T (1)with g ∗ ( T ) being the effective number of relativistic de-grees of freedom at temperature T . In the absence of en-tropy production per comoving volume i.e. sa = const.,radiation energy density redshifts as ρ rad ( t ) ∝ a ( t ) − . Incase of a rapid expansion of the universe the energy den-sity of ϕ field is expected to be redshifted faster than theradiation. Accordingly, one can assume ρ ϕ ∝ a ( t ) − (4+ n ) .Thus n > ϕ energy density dominates overradiation at early enough times.Now, the entropy density (per comoving volume) of theuniverse is expressed as s ( T ) = 2 π g ∗ s ( T ) T where g ∗ s ( T ) is the effective relativistic degrees of free-dom. A general form of ρ ϕ can then be constructed using GW signal in a modified cosmological scenario have also been stud-ied in details in [18–20] the entropy conservation g (cid:63) ( T ) / aT = constant in a co-moving frame as: ρ ϕ ( T ) = ρ ϕ ( T R ) (cid:18) g ∗ s ( T ) g ∗ s ( T R ) (cid:19) (4+ n ) / (cid:18) TT R (cid:19) (4+ n ) . (2)The temperature T R is an unknown variable and canbe considered as the point of equality where the two fluidshave equal energy densities: ρ ϕ ( T R ) = ρ rad ( T R ). Usingthis criteria, the total energy density at any temperature T reads [11, 12] ρ ( T ) = ρ rad ( T ) + ρ ϕ ( T ) (3)= ρ rad ( T ) (cid:34) g ∗ ( T R ) g ∗ ( T ) (cid:18) g ∗ s ( T ) g ∗ s ( T R ) (cid:19) (4+ n ) / (cid:18) TT R (cid:19) n (cid:35) (4)Now, from standard Friedman equation one can writethe Hubble parameter in terms of the energy density as H = √ ρ √ M Pl , (5)with M Pl = (8 πG ) − / = 2 . × GeV being the re-duced Planck mass. At temperature higher than T R withthe condition g ∗ ( T ) = ¯ g ∗ which can be considered to besome constant, the Hubble rate can approximately berecasted as [11] H ( T ) = H R ( T ) (cid:34) (cid:32) TT R (cid:33) n (cid:35) / ≈ π ¯ g / ∗ √ T M Pl (cid:18) TT R (cid:19) n/ , ( T (cid:29) T R ) , (6)where H R ( T ) is the Hubble parameter in the standardradiation dominated universe. In case of SM ¯ g ∗ ≡ g ∗ (SM) = 106 .
75. It is important to note from Eq. (6)that the expansion rate is larger than what it is sup-posed to be in the standard cosmological background for
T > T R and n >
0. It is noteworthy that T R can not betoo small such that it changes the standard BBN history.For a certain value of n , BBN constraints provide a lowerlimit on T R [11]: T R (cid:38) (15 . /n MeV . (7)As one can see, for larger n this constraint becomesmore and more lose. Before moving on to the next sectionwe would like to comment that in order to get valueslarger than n = 2, a negative scalar potential is to beconsidered. A specific structure of n > III. DARK MATTER YIELD VIA FREEZE-IN
In this section we will discuss the modification of stan-dard Boltzmann equation [10] governing the DM numberdensity due to the effect of fast expansion encoded inthe Hubble parameter H (Eq. (6)). The key for freeze-in [22, 23] DM production is to assume that DM was notpresent in the early universe. The DM is then producedvia annihilation and/or decay of the particles in the ther-mal bath. Due to the extremely feeble coupling of theDM with the visible sector particles the DM never reallyenters into thermal equilibrium. Now, the Boltzmannequation (BEQ) for DM production via annihilation for2 → a, b → , n + 3 H n = (cid:90) d Π d Π d Π a d Π b (2 π ) δ ( p + p − p a − p b ) |M| a,b → , f a f b , (8)where d Π j = d p j E j (2 π ) are Lorentz invariant phasespace elements, and f i is the phase space density of par-ticle i with corresponding number density being: n i = g i (2 π ) (cid:90) d pf i , (9)with g i is the internal degrees of freedom (DOF). Herewe have made two crucial assumptions: the initial abun-dance for particle 1 is negligible such that f ≈ ± f i ≈
1. The BEQ in Eq. (8) canbe rewritten as an integral with respect to the CM energyas [13, 22]:˙ n + 3 H n ≈ T π (cid:90) ∞ s =4 m χ ds d Ω P ab P |M| ab → √ s K (cid:18) √ sT (cid:19) , (10)where P ij = √ s (cid:112) s − ( m i + m j ) (cid:112) s − ( m i − m j ) → √ s in the limit m i,j →
0. The BEQ in terms of the yield Y χ = n /s can be written in the differential form as: − s ( T ) . H ( T ) .T. dY ann χ dT = T π (cid:90) ∞ s =4 m χ dsd Ω P ab P |M| ab → √ s K (cid:18) √ sT (cid:19) . (11)The total yield due then turns out to be: Y χ ( T ) = 1512 π (cid:90) T max T dT (cid:48) s ( T (cid:48) ) . H ( T (cid:48) ) (cid:90) ∞ s =4 m χ dsd Ω (cid:18) √ s (cid:19) |M| ab → √ s K (cid:18) √ sT (cid:48) (cid:19) (12) where we have put m i,j = 0 just for convenience. Theupper limit of the T integration is considered to be T max = T RH corresponding to the reheat temperatureof the universe. Eq. (12) should be divided into two eras:one before EWSB with T >
160 GeV and the other post-EWSB for T (cid:46)
160 GeV. For T RH (cid:29) T EW the DM pro-duction is dominant before EWSB. The lower limit onthe s integral can be approximated to zero if the DMmass is negligible compared to reheat temperature of theuniverse. For a lower reheat temperature, however, onecan not ignore the masses of the particle spectrum inthe theory and for T RH ∼ m i the IR nature of the yieldbecomes prominent [25]. Since in a fast expanding uni-verse the Hubble parameter gets modified according toEq. (6), the resulting DM yield also gets affected follow-ing Eq. (12). In the next section we shall quantify thisthrough analytical and numerical computation. Finally,the relic abundance of the DM at present temperaturecan be obtained via:Ω DM h = (cid:0) . × (cid:1) m DM GeV Y totalDM ( T ) , (13)where we find Y totalDM ( T ) by solving Eq. (12) numeri-cally. The important point here to note is that due tomodified Hubble rate H the DM yield and consequentlythe DM abundance is now a function of the fast expansionparameters n, T R on top of DM mass, reheat temperatureand the effective scale. A. A toy model
We consider a minimal dark sector comprising of aDirac fermion singlet χ and a scalar singlet φ . Thereexists a Z symmetry under which only the dark sectorparticles are charged as tabulated in Tab. I. We assumethe fermion χ to be the lightest particle charged underthe new symmetry and hence a potential DM candidate.With this symmetry assigned one can then write the ef-fective Lagrangian for the dark sector as: Symmetry χ φZ i -1 Table I. Transformation of the new particles under Z sym-metry. L ⊃ χ (cid:0) i /∂ − m χ (cid:1) χ − y χ χ c χφ + 1Λ χχ | H | − V ( H, φ ) (14)with the renormalizable scalar potential given by V ( H, φ ) ⊃ − µ H | H | + µ φ φ + λ H | H | + λ φ φ + λ Hφ | H | φ . (15)The Z symmetry prevents us from writing any oddterms in the scalar potential. Note that the portal cou-pling λ Hφ ∼ O (1) keeps the singlet scalar φ in equilib-rium with the thermal bath. As the new singlet scalar φ does not acquire any non-zero VEV, χ has only apure Dirac mass term m χ . Because of the Z symme-try it is also not possible to write Majorana mass termsfor χ that otherwise would have resulted in a pseudo-Dirac splitting [26]. The second term in Eq. (14), due toits renormalizable nature, contributes to the DM abun-dance via Infra-Red (IR) freeze-in due to φ → χχ decaywhich is kinematically accessible as long as m χ < m φ / HH → χχ contact interaction be-fore EWSB while after EWSB one can have both UVand IR freeze-in contribution due to non-zero VEV ofthe Higgs that results in renormalizable DM-SM inter-action. For UV freeze-in the abundance is sensitive tothe maximum temperature of the thermal bath, whichwe assume to be the reheat temperature T RH2 . This isin sharp contrast to the IR freeze-in scenario where thetwo sectors communicate via renormalizable operators,and the DM abundance is set by the IR physics i.e., the yield becomes maximum at low temperature, typi-cally at T ∼ m DM [22, 23]. Now, the reheat tempera-ture T RH of the universe is somewhat loosely bounded.Typically, the lower bound on T RH comes from the mea-surement of light element abundance during BBN, whichrequires T RH (cid:38) . T RH (cid:46) GeV to prohibit thermal grav-itino over production and (ii) simple inflationary scenar-ios that require at most T RH ∼ GeV [30, 31] for asuccessful inflation. In view of this, the maximum reheattemperature of the universe can be regarded as a freeparameter.
B. Relic abundance of the dark matter
A model-independent study of IR freeze-in in contextof fast expansion has been done in [12]. Here we wouldlike to see the effect of such fast expansion in the con-text of pure UV freeze-in before electroweak symmetrybreaking, together with IR freeze-in, once the electroweaksymmetry is broken. The requirement of obtaining rightDM relic abundance via UV freeze-in shall also put somebound on the reheat temperature, along with the cut-offscale Λ of the theory in the fast expanding scenario. Upto this end we note that the contribution to DM abun-dance via IR freeze-in through φ → χχ is negligible once In principle, the maximum temperature during reheating can belarger than T RH [27]. we choose the Yukawa coupling y χ (cid:46) −
12 3 . Such achoice of coupling might look like somewhat “fine-tuned”but the smallness of coupling for IR freeze-in is in itself arequirement to keep the DM out of thermal equilibrium.The small Yukawa coupling thus prevents the dark sectorto reach thermal equilibrium within itself without ham-pering the rest of our analysis. This makes UV freeze-inas the dominant process that can produce the total ob-served DM relic via decay and /or annihilation which wewill discuss below.In order to calculate the DM relic abundance by solvingEq. (12) we first point out the possible annihilation anddecay channels leading to DM pair production above andbelow the electroweak symmetry breaking temperature T (cid:39)
160 GeV. Before EWSB i.e.,
T >
160 GeV, for the SU (2) L scalar doublet, the propagating DOFs are theGoldstone bosons (GB): H = (cid:18) ϕ + ϕ (cid:19) , (16)which are then converted to the longitudinal DOF for themassive gauge bosons after EWSB in unitary gauge: H = (cid:18) h + v h √ (cid:19) , where h is the CP-even SM-like Higgs field. After elec-troweak symmetry breaking (EWSB) i.e., T (cid:46)
160 GeVthere can be several 2 → → m h > m χ (see Fig. 1). Another SM Higgs medi-ated 2 → φ in the initial state. Theamplitude for this process goes roughly as λ Hφ v h / Λ .For λ Hφ (cid:46) O (1) and Λ (cid:38) O (cid:0) (cid:1) GeV, however, thisprocess has negligible contribution to the DM yield viaIR freeze-in. Using Eq. (12) it is possible to derive anapproximate analytical expression for DM yield via UVfreeze-in (before EWSB). As all the SM particles are alsomassless before EWSB, the squared amplitude for the4-point processes simply become |M| (cid:39) s Λ . (17)Considering the DM and φ mass to be negligible com-pared to { Λ , T RH } , for T (cid:29) T R , we can analytically ob-tain the expression for DM yield before EWSB For example, for DM mass of 20 GeV and m φ = 70 GeV choosing y χ = 10 − gives rise to DM relic abundance of ∼ − dominantlyvia φ → χχ decay with a lifetime of ∼ . Figure 1. DM yield via freeze-in before EWSB (top left)andafter EWSB (top right and bottom). Here SM stands for allthe SM gauge bosons, quarks and leptons. Here H stands forthe SM-like Higgs field before EWSB and h is the CP-evenHiggs field after EWSB. Y χ = ξ (cid:16) T R T EW (cid:17) n/ T EW Λ ( n − (cid:34) − T RH T EW (cid:32) T EW T RH (cid:33) n/ (cid:35) , ( n (cid:54) = 2) ξT R log (cid:104) T RH T EW (cid:105) ( n = 2) (18)where ξ = M pl √ π is a constant of mass dimensionunity. After EWSB there will be IR dominated processesas well which are proportional to ( v/ Λ) . Now, since allthe SM particles become massive after EWSB (togetherwith the DM), hence one has to find the DM yield nu-merically in this regime. The only possible decay chan-nel contributing to DM abundance is h → χχ which isavailable if m χ < m h /
2. We can obtain an approximateanalytical expression for the yield due to decay [12] as Y decay χ = (cid:90) T max T min dT m h Γ h → χχ π K ( m h /T ) s ( T ) H ( T ) ≈ Γ h → χχ T n/ R π m n +2 h M pl A (cid:90) ∞ dxx n/ K ( x )= 2 M pl Γ h → χχ Aπ m h (cid:32) T R m h (cid:33) n/ (cid:34) n/ Γ (cid:16) n (cid:17) Γ (cid:16)
10 + n (cid:17)(cid:35) , (19)where A contains all the dimensionless constants and weassumed T max → ∞ (which, in reality, should be T EW )such that x max → T min = 0. The decay widthgoes roughly as Γ h → χχ ∝ ( v/ Λ) . In the case of 2 → Y → χ ∼ m χ M pl B Λ (cid:32) T R m χ (cid:33) n/ (cid:90) x max x min dxx n/ − , (20) with x max = m χ /T and x min = m χ /T EW , where B contains all the dimensionless constants. Taking all themasses of the particles into account, however, such a sim-ple analytical expression would be difficult to obtain andone has to resort to numerical computation. The asymp-totic nature of the yield after EWSB can be determinedfrom Eq. (19) and Eq. (20) which tells that the yield issuppressed for larger n as x → ∞ . Under the frame-work of faster expansion one can not vary T RH , T R and n arbitrarily as pointed out in [12]. The dominance of ρ ϕ at early universe sets an upper bound on the reheattemperature since beyond ρ ϕ ∼ M this framework fails.This, together with Eq. (6), sets an upper bound on thereheat temperature [12] H ∼ T M pl (cid:16) T RH /T R (cid:17) n/ (cid:46) M pl = ⇒ T RH (cid:46) M pl (cid:16) M − T R (cid:17) n/n +4 . (21)Note that the above bound gets stronger for larger n .Additionally, for the effective theory to remain validΛ > m χ , T RH is mandatory. For processes appearingbefore and after EWSB we compute the DM yield byperforming numerical integrations using Mathematica .In Fig. 2 we have illustrated how the DM yield varieswith x = m χ /T for a fixed T R = 40 MeV and reheattemperature T RH = 10 GeV for different choices of n -values for different DM masses m χ = { , } GeV. Forall these plots we have kept T RH = 10 GeV which im-plies n (cid:46) x → ∞ since all of them sat-isfy the DM relic abundance for a fixed DM mass of 50GeV (on left panel) and 500 GeV (on right panel). Alsonote that the DM yield builds up very early at a smaller x ( T ∼ T RH ) which is the typical feature of UV freeze-inwhen T RH (cid:29) m [13, 25]. For larger n , as one can no-tice from Eq. (12), the DM yield gets suppressed since H becomes larger according to Eq. (6) for a fixed T RH .This shows that the freeze-in production of DM in a fastexpanding universe is dramatically suppressed which hasalso been reported earlier [12]. Thus, a larger n requires asmaller Λ to achieve right DM abundance for a fixed T RH as Y χ ∝ . The orange curves ( n = 4) in Fig. 2 corre-spond to Λ ∼ GeV which is about five orders less thanthat for the red curves (standard case) corresponding towhich Λ ∼ GeV. This exactly what is reflected inFig. 3 where we show DM abundance satisfying contoursin Λ- T RH plane for DM mass of m χ = { , } GeV inthe left and in the right panel respectively consideringall processes appearing before and after EWSB. For a 50GeV DM h → χχ decay contributes significantly afterEWSB via IR freeze-in [32]. As a result we see the hori-zontal part in the red and green curves for T RH (cid:46)
10 TeV.This arises due to the IR nature of the freeze-in where theDM abundance depends no more on the reheat tempera-ture, rather on the mass of the decaying mother particle(which is Higgs in this case). For larger n although the (a) (b) Figure 2. (a): Variation of DM yield Y χ with x = m χ /T for different values of n = { , , , } shown in green, blue, black andorange respectively for m χ = 50 GeV, T RH = 10 GeV and T R = 40 MeV. Each curve satisfies observed relic abundance fordifferent values of the cut-off scale Λ ≈ { , , , } GeV respectively. (b): Same as (a) for m χ = 500 GeV. Thedashed black straight line in both plots corresponds to x EW = m χ /T EW . (a) (b) Figure 3. (a): Contours satisfying PLANCK observed relic abundance in the bi-dimensional plane of Λ- T RH for different choicesof n = { , , } shown in green, blue and orange respectively where the DM mass is fixed at 50 GeV and T R = 40 MeV. (b):Same as (a) for m χ = 500 GeV. curves look absolutely flat ( e.g., the blue one) but as onecan see from the inset figure their variation over Λ is al-most negligible compared to the standard case or n = 1case. For a heavier DM mass, typically m χ (cid:29) T EW , mostof the contribution to the DM yield comes from processesbefore EWSB. As a consequence we observe typical UVnature in the contours on the right panel of Fig. 3 wherethe cut-off scale rises linearly with the reheat tempera-ture. For small enough T RH (cid:28) T RH (cid:38) ∝ Λ − withrespect to the standard cosmological background in or-der to get the right DM abundance in a fast expandingbackground. This implies that for non-renormalizable in- teractions with dimension larger than five, it is possibleto bring down Λ ∼ TeV for n > n the DM yield gets moreand more suppressed, hence we stick only up to n = 4.For such a choice of n , as we have illustrated, it is possi-ble to obtain observed DM relic abundance for differentchoices of T R and T RH obeying Eq. (7) and Eq. (21) bysuitably modifying the cut-off scale Λ. As we shall see inthe next section, for such a choice of n, T R it is possibleto generate a strong first order phase transition that iscapable of explaining the NANOGrav 12.5 yrs data. IV. STOCHASTIC GRAVITATIONAL WAVESIGNAL IN A FAST EXPANDING UNIVERSE
In general, an extended scalar sector may give rise toa strong first order phase transition (SFOPT)resulting instochastic gravitational wave (GW) with detectable fre-quencies (see, for example [33]). In our case the presenceof an extra scalar singlet φ does that job. In the contextof the present model we would like to explore how thefrequency and amplitude of such a stochastic GW signalcan get modified due to the fast expansion of the uni-verse. Our particular interest is to explore if a simpledark sector can explain the NANOGrav data when thecosmological background is modified. As we shall see,the stochastic GW signal that can explain the reportedNANOGrav data [1] can only be generated in a fast ex-panding universe for some specific choices of T R , n .
1. Finite Temperature Effective Potential
To construct the finite temperature effective potential,we add two temperature corrected terms with the tree-level potential. The effective potential at a finite temper-ature T then can be written as [34] V eff = V tree-level + V T =01 − loop + V T (cid:54) =01 − loop , (22)where V tree-level , V T =01 − loop and V T (cid:54) =01 − loop are the zero tem-perature tree level potential, zero temperature Coleman-Weinberg (CW) one-loop potential and the finite tem-perature one-loop potential respectively. As mentionedearlier, the scalar sector under consideration contains areal scalar singlet on top of the SM-like Higgs doubletwhich leads to Eq. (15). The SM Higgs H can be repre-sented as H = G + √ (cid:0) v + h + iG (cid:1) , (23)where G + and G are the charged and the neutralGoldstone bosons after the spontaneous symmetry break-ing. To explore the electroweak phase transition we re-place the scalar fields H and φ with their VEVs v , v φ respectively. Now Eq. (15) can be expressed as V T =0tree − level = − µ H v + 12 µ φ v φ + 14 λ H v + 14 λ φ v φ + 14 λ Hφ v v φ . (24)Note that, although the scalar φ do not acquire anyvacuum expectation value (VEV) at T = 0 (to keepthe Z symmetry intact), but its VEV can be generatedat a finite temperature. Here the classical VEVs v, v φ change with temperature i.e., behave like a field and at the zero temperature it tends to the classical fixed values.The zero-temperature CW one-loop effective potential isgiven by [34, 35] V T =01 − loop = ± π (cid:88) i n i m i (cid:20) log m i Q − C i (cid:21) , (25)where the ‘+’ and ‘ − ’ symbols appear due to bosonsand fermions respectively. The quantity Q denotes therenormalization scale which we take Q = v = 246 . i is over all theparticle species associated in the scalar potential with i ∈ ( h, φ, G ± , G , W, Z, t ). In Eq. (25), the quantities n i , m i and C i represents the number of degrees of freedom,the field-dependent masses and renormalization-scheme-dependent numerical constant of the i th particle speciesrespectively. The associated degrees of freedom of theparticles are ( n W ± ) L = 4, ( n W ± ) T = 2, ( n Z ) L = 2,( n Z ) T = 1, n t = 12, n G ± = 2 and n h,φ,G = 1. Thevalues of the renormalization-scheme-dependent constant C i are ( C W,Z ) T = 1 / W, Z boson, ( C W,Z ) L = 3 / W, Z boson and for the other particles, C h,φ,G + ,G − ,G ,t = 3 /
2. We perform our calculationsby considering the Landau gauge where the Goldstonebosons are massless at T = 0 and also the theory isfree from ghost contributions [36]. With this, the finitetemperature one-loop effective potential can be expressedas [34] V T (cid:54) =01 − loop = T π (cid:88) i n i J ± (cid:20) m i T (cid:21) , (26)where J ± (cid:18) m i T (cid:19) = ± (cid:90) ∞ dy y log (cid:32) ∓ e − (cid:118)(cid:117)(cid:117)(cid:116) y + m i T (cid:33) . (27)We include the thermal correction to the boson massesby applying daisy resummation method [37] as µ (cid:48) H ( T ) = µ H + c T and µ (cid:48) φ ( T ) = µ φ + c T where c = 6 λ H + 2 λ Hφ
12 + 3 g + g (cid:48)
16 + y t c = 6 λ φ + 2 λ Hφ , (29)with the couplings g , g (cid:48) and y t as the SU (2) L gaugecoupling, U (1) Y gauge coupling and top Yukawa couplingof the SM respectively. We also include the thermal cor-rections to the W and Z boson masses with the modifiedthermal masses [36] m W ( T ) = 14 g v + 2 g T , (30)and m Z ( T ) = 18 (cid:16) g + g (cid:48) (cid:17) v + (cid:16) g + g (cid:48) (cid:17) T + 18 (cid:114)(cid:104)(cid:0) g − g (cid:48) (cid:1) (64 T + 16 T v ) + (cid:0) g + g (cid:48) (cid:1) v (cid:105) . (31)With this set-up we then proceed to analyze the GWproduction due to strong first order phase transition(SFOPT).
2. Gravitational wave production from SFOPT
Here we discuss the possible production mechanism ofGW from the strong first-order electroweak phase tran-sition originating mainly via bubble collisions [38–44],sound waves induced by the bubbles running throughthe cosmic plasma [45–48] and turbulence induced by thebubble expansions in the cosmic plasma [49–53]. The to-tal GW intensity Ω GW h for a particular frequency f can be estimated by considering the contributions of theabove three processes and can be expressed as [38–55]Ω GW h = Ω col h + Ω SW h + Ω turb h . (32)The expression for the GW intensity from the compo-nent of bubbles collision Ω col h is given byΩ col h = 1 . × − (cid:18) β H (cid:19) − . v w .
42 + v w (cid:18) κα α (cid:19) (cid:16) g ∗ (cid:17) − . (cid:18) ff col (cid:19) . . (cid:18) ff col (cid:19) . , (33)where the parameter β has the form β = (cid:20) H T ddT (cid:18) S T (cid:19)(cid:21) (cid:12)(cid:12)(cid:12)(cid:12) T n , (34)where S ( T ) is the Euclidean action of the critical bub-ble. In general, the nucleation of the bubbles occur at atemperature T n where the condition S ( T n ) /T n ≈ H denotes the Hubble pa-rameter at the nucleation temperature T n which is givenby Eq. (6) for the fast expanding universe. The bub-ble wall velocity v w is computed using the most generalexpression for the wall velocity [39, 43, 56, 57] v w = 1 / √ (cid:112) α + 2 α/
31 + α . (35)The quantity κ in Eq. (33) is given by κ = 1 − α ∞ α , (36)with [58, 59] α ∞ = 3024 π g ∗ (cid:18) v n T n (cid:19) (cid:34) (cid:16) m W v (cid:17) + 3 (cid:16) m Z v (cid:17) + 6 (cid:16) m t v (cid:17) (cid:35) . (37)where v n is the Higgs VEV at the nucleation temper-ature T n , m W , m Z and m t are the masses of the gaugebosons W , Z and the top quark t respectively. The phasetransition parameter α can be expressed as α = (cid:20) ρ vac ρ ∗ rad (cid:21) (cid:12)(cid:12)(cid:12)(cid:12) T n . (38)where ρ rad ∗ is the background energy density of theplasma and ρ vac is the energy density difference betweenfalse and true vacuum during the electroweak phase tran-sition. The quantities ρ vac and ρ rad ∗ have the followingform ρ vac = (cid:34)(cid:32) V higheff − T dV higheff dT (cid:33) − (cid:18) V loweff − T dV loweff dT (cid:19)(cid:35) , (39)and ρ ∗ rad = g ∗ π T n . (40)In Eq. (33) the peak frequency f col from the bubble col-lisions is given by f col = 16 . × − Hz (cid:18) . v w − . v w + 1 . (cid:19) (cid:18) β H (cid:19)(cid:18) T n
100 GeV (cid:19) (cid:16) g ∗ (cid:17) . (41)Next, the expression for GW intensity from the com-ponent of sound wave (SW) is given byΩ SW h = 2 . × − (cid:18) β H (cid:19) − v w (cid:18) κ v α α (cid:19) (cid:16) g ∗ (cid:17) − (cid:18) ff SW (cid:19)
74 + 3 (cid:18) ff SW (cid:19) , (42)where the parameter κ v κ v = α ∞ α (cid:20) α ∞ .
73 + 0 . √ α ∞ + α ∞ (cid:21) . (43)The peak frequency f SW from the sound wave contribu-tion reads f SW = 1 . × − Hz (cid:18) v w (cid:19) (cid:18) β H (cid:19) (cid:18) T n
100 GeV (cid:19) (cid:16) g ∗ (cid:17) . (44)For the contribution of SW to the total GW intensitydepends on the Hubble time scale. If it survives morethan a Hubble time then the expression in Eq. (42) willbe valid, otherwise we need to include a factor called suppression factor to the SW component of the GW in-tensity. Following [59–61] we compute the suppressionfactor H R ∗ ¯ U f (where ¯ U f is the root-mean-square (RMS)fluid velocity and R ∗ is the mean bubble separation) tocheck whether the SW components lasts more than aHubble time. Figure 4. Plot of calculated GW intensities as a functionof frequency for different choices of of { n, T R } . The sensitiv-ity curves of NANOGrav11, NANOGrav12.5 (purple colouredrectangular box), future GW detectors such as ALIA, BBO,DECIGO, LISA, aLIGO and aLIGO+ are also shown in thesame plane. Lastly, the expression for GW intensity from the tur-bulence in the plasma Ω turb h is given by Ω turb h = 3 . × − (cid:18) β H (cid:19) − v w (cid:18) (cid:15)κ v α α (cid:19) (cid:16) g ∗ (cid:17) − (cid:18) ff turb (cid:19) (cid:18) ff turb (cid:19) − (cid:18) πfh ∗ (cid:19) , (45)with h ∗ = 16 . × − Hz (cid:18) T n
100 GeV (cid:19) (cid:16) g ∗ (cid:17) . (46)Here (cid:15) = 0 . f turb from the tur-bulence mechanism has the form f turb = 2 . × − Hz (cid:18) v w (cid:19) (cid:18) β H (cid:19) (cid:18) T n
100 GeV (cid:19) (cid:16) g ∗ (cid:17) , (47) BP m h m φ λ H λ φ λ Hφ in GeV in GeV1 125.09 70 0.129 0.1 0.1 Table II. Chosen benchmark point to demonstrate strong first-order phase transition. BP v c T c v c T c v n T n α in GeV in GeV in GeV in GeV1 140.50 134 1.04 151.20 134.58 0.0028 Table III. Values of thermal parameters for the chosen BP.
To investigate the observational signatures of suchstochastic GW, our estimated model-dependent GWintensities are compared with the future space-basedand ground-based detectors such as ALIA, BBO, DE-CIGO, aLIGO, aLIGO+ and pulsar timing arrays suchas NANOGrav11 and NANOGrav12.5. We are typicallyinterested in the GW spectra from FOPT in a fast ex-panding universe which lie within the NANOGrav sensi-tivity. For calculating the thermal parameters relatedto the phase transition we use the publicly available
CosmoTransition package [34] and compute the totalGW intensities using Eqs. (32)-(46). We choose a bench-mark point (BP) with m φ = 70 GeV which is safe fromHiggs invisible decay constraint. Note that all the scalarcouplings giving rise to strong first order PT are λ i ∼ . { v c , T c , v n , T n , α } for thechosen BP are tabulated in Tab. III. The strength of theFOPT depends on the order parameter v c /T c and thetransition is said to be strong enough if it satisfies thecondition v c /T c (cid:62)
1, which in our case turns out to be0 v c /T c = 1 .
04. In Tab. IV we present the results for GWintensity and corresponding peak frequency for differentvalues of n and T R . From Tab. IV and Fig. 4 we ob-serve that the GW intensity increases and peak shiftsto the lower frequency region for smaller values of T R at a fixed n . The GW intensities as a function of fre-quency are plotted in Fig. 4 for different choices of n and T R . We compare our model-dependent GW intensitieswith the sensitivity plots of PTAs such as NANOGrav11,NANOGrav12.5 and future space-based, ground-basedgravitational wave detectors such as ALIA, BBO, DE-CIGO, LISA, aLIGO and aLIGO+ . From the Fig. 4one observes that for n = 3, T R = 40 MeV and n = 4, T R = 200 MeV the GW intensity lies within the sensitiv-ity curves of NANOGrav-12.5. There are other choices of { n, T R } that also give rise to detectable GW intensitieswhich fall within the reach of detectors like ALIA, BBO,DECIGO, LISA and NANOGrav11. Figure 5. The blue and black solid lines correspond to thepredictions of γ CP and A CP (with A CP the characteristic GWstrain amplitude at f = yr − and γ CP the spectral index of thepulsar timing-residual cross-power spectral density) obtainedfor different fixed values of { n, T R } which gives rise to Ω GW h required to explain the NANOGrav result. The dashed andsolid red curves correspond to to the 1 σ and 2 σ γ CP - A CP contours obtained by NANOGrav 12.5 yrs dataset. The results of PTA GW searches are generally reportedin terms of a pulsar timing-residual cross-power spectraldensity which has a frequency dependence of the form f − γ CP . Corresponding power spectrum of the character-istic GW strain is usually approximated as a power-lawof the form h c ( f ) = A CP (cid:32) ff yr (cid:33) α CP (48)with α CP = (3 − γ CP ) / In this work we use the power-law-integrated sensitivity curves fol-lowing [57, 62]. One can also use the other methods to representthe sensitivity curves [63–65] in circular orbits whose evolution is dominated by GWemission. Most importantly, The reference frequency isgiven by f = yr − . The power spectrum is related to theGW energy density via [66]Ω GW = 2 π H f h c ( f ) . (49) n T R β H f peak Ω GW h in MeV in Hz- - 5264.70 3.45 × − × − × − × − × − × − × − × − × − × − × − × − Table IV. The calculated values of peak frequencies and cor-responding GW intensity for the standard cosmology scenario(first row) and non-standard scenarios with different choicesof T R . The GW search results due to PTA are reported interms of joint A CP − γ CP posterior distributions or as theposterior distribution at a fiducial value of γ CP [1, 66].NANOGrav collaboration has obtained joint constraintson A CP and γ CP by fitting the power-law in Eq. (48) [1].Since Eq. (49) is a function of n, T R (through H ), hence itis possible to project limits on the fast expansion param-eters by exploiting Eq. (48). Now, as we have shown inFig. 4, corresponding to n = 3 , n − T R constraints from Fig. 4 onto the A CP − γ CP planewith corresponding 1 σ and 2 σ contours from NANOGravresult [1]. Here we observe that for a fixed n = 4 a larger T R shifts the straight line contour (orange) to the edgeof the 2 σ contour (red thick), while smaller n (black) iscloser to the 1 σ contour (red dashed). From Fig. 4 and 5we see that NANOGrav results not only support a fasterexpansion in the pre-BBN era, also favour n (cid:46) T R . Based on the analysis in Sec. IIIand Sec. IV we can infer that while right relic abundancefor the DM can be obtained for any n > T R butexplaining NANOGrav results from a SFOPT togetherwith the right DM abundance necessarily demands n (cid:46) T R (cid:46) O (100) MeV. V. POSSIBLE UV COMPLETION
Here we would like to elucidate a possible UV com-plete framework of the effective Lagrangian we discussedin Eq. (14). There might be several possibilities to ob-tain an effective interaction of the form in Eq. (14). Foran instance, the authors in [67] have discussed one suchalternative. However, in their case the Z is broken by1the non-zero VEV of the singlet scalar φ leading to aremnant Z symmetry under which the extra fermionsare odd (DM). On the other hand, in [68] a Z symme-try has been imposed under which all the fields in thedark sector are odd. Here we would like to furnish thesimplest possible model remembering the fact that in thepresent scenario the singlet scalar does not acquire anytree-level VEV. This is in contrast to the set-up discussedin [69] where the singlet fermionic WIMP-like DM is sta-bilized by a Z symmetry and the singlet scalar acquiresa VEV providing a strongly first order electroweak phasetransition. Fields SU (2) L U (1) Y Z χ iψ T : (cid:0) ψ , ψ − (cid:1) − iφ − Table V. Charges of different fields appearing in a possibleUV completion of the Lagrangian in Eq. (14).
We assume a dark sector that is made up of a vectorlike lepton (VLL) singlet χ and a SU (2) L VLL doublet ψ with the charges as mentioned in Tab. V. As mentionedearlier, we consider the new particles transform under a Z symmetry with their charges as in Tab. V. With thegiven SM × Z symmetry the renormalizable Lagrangiancan be written as L ⊃ m χ χχ + M ψ ψψ + Y ψ (cid:101) Hχ + y χ χ c χφ + H.c. (50)At a scale µ (cid:28) M ψ the heavy VLLs can be integratedout, resulting in operators of the form (cid:0) H † H (cid:1) ( χχ ) /M ψ .If the reheat temperature of the universe is smaller than M ψ then the additional VLLs are not present in the ther-mal bath and the resulting DM-SM interaction mimics adimension five operator with cut-off scale Λ ≡ M ψ . VI. CONCLUSION
In this paper we have studied a simple frameworkwhere it is possible to address both dark matter (DM)relic abundance via freeze-in and the recent NANOGravresult from a strong first order phase transition (SFOPT)in a fast expanding universe. In order to accomplishthem together, we consider a dark sector consisting ofa SM gauge singlet fermion χ and a scalar singlet φ bothcharged under a Z symmetry. We consider χ to be thelightest stable fermion in the dark sector that can serve asa viable DM candidate. Being singlet under SM gaugesymmetry the dark fermion does not have any interac-tions in the renormalizable level with the SM fields. Thescalar, on the other hand, remains in equilibrium withthe thermal bath via non-negligible portal coupling andis responsible in generating SFOPT.We focus on the DM genesis via freeze-in, typically UVfreeze-in, as the dark sector and the visible sector are con- joined through non-renormalizable interaction of dimen-sion five. The DM yield is computed by taking into ac-count all processes leading to DM production both beforeand after the electroweak symmetry breaking (EWSB).While the processes before EWSB (dominantly UV) areonly 2 → T R , n for fixed DM masses.This also put constraints on the other two relevant freeparameters of the model, namely the effective scale Λand the reheat temperature T RH which we consider to T RH (cid:38) n >
0. This, in turn, implies that to produce enoughDM via freeze-in to match observations, one has to bringdown the effective scale of interaction (equivalently, theDM-SM coupling) substantially with respect to the stan-dard cosmological scenario.The presence of the scalar singlet facilitates strong firstorder phase transition (SFOPT) which is otherwise im-possible to achieve within the SM particle content. In thepresence of fast expansion, we observe, such SFOPT givesrise to gravitational wave (GW) signal with intensitiesthat are detectable by several space-based and ground-based GW detectors. In a fast expanding universe suchstrong GW signals can be generated for λ Hφ (cid:38) . n (cid:46) σ posterior contours for ampli-tude and spectral slope obtained by the NANOGrav col-laboration. For all such n (cid:46) (cid:38) GeV that can well explain observedDM abundance for suitable choice of the temperature T R ∼ O (100) MeV. Thus, a simple framework of thekind discussed here connects reheat temperature of theuniverse with a SFOPT in a fast expanding cosmologicalbackground. A possible underlying UV complete theorywhich can give rise to such dimension-5 effective interac-tions have also been discussed briefly. It will however beinteresting, in the context of a UV complete theory, tolook into interplay of different parameters giving rise toDM relic abundance and a stochastic GW signal in a fastexpanding universe, together with more exotic physics.We shall address such issues in forthcoming works. ACKNOWLEDGEMENTS
BB would like to acknowledge discussions with Abhi-jit Kumar Saha, Debasish Borah and e-mail conversa-2tions Sunny Vagnozzi. AP would like to thank BiswajitBanerjee for the implementation of the
CosmoTransition package. Work of ADB is supported in part by the Na-tional Science Foundation of China (11775093, 11422545,11947235). [1] Z. Arzoumanian et al. (NANOGrav), (2020),arXiv:2009.04496 [astro-ph.HE].[2] A. Sesana, F. Haardt, P. Madau, and M. Volonteri, As-trophys. J. , 623 (2004), arXiv:astro-ph/0401543.[3] J. Ellis and M. Lewicki, (2020), arXiv:2009.06555 [astro-ph.CO].[4] S. Blasi, V. Brdar, and K. Schmitz, (2020),arXiv:2009.06607 [astro-ph.CO].[5] S. Chigusa, Y. Nakai, and J. Zheng, (2020),arXiv:2011.04090 [hep-ph].[6] A. Neronov, A. Roper Pol, C. Caprini, and D. Semikoz,(2020), arXiv:2009.14174 [astro-ph.CO].[7] Y. Nakai, M. Suzuki, F. Takahashi, and M. Yamada,(2020), arXiv:2009.09754 [astro-ph.CO].[8] A. Addazi, Y.-F. Cai, Q. Gan, A. Marciano, andK. Zeng, (2020), arXiv:2009.10327 [hep-ph].[9] S. Bhattacharya, S. Mohanty, and P. Parashari, (2020),arXiv:2010.05071 [astro-ph.CO].[10] E. W. Kolb and M. S. Turner, Front. Phys. , 1 (1990).[11] F. D’Eramo, N. Fernandez, and S. Profumo, JCAP ,012 (2017), arXiv:1703.04793 [hep-ph].[12] F. D’Eramo, N. Fernandez, and S. Profumo, JCAP ,046 (2018), arXiv:1712.07453 [hep-ph].[13] F. Elahi, C. Kolda, and J. Unwin, JHEP , 048 (2015),arXiv:1410.6157 [hep-ph].[14] S.-L. Chen and Z. Kang, JCAP , 036 (2018),arXiv:1711.02556 [hep-ph].[15] N. Bernal, F. Elahi, C. Maldonado, and J. Unwin, JCAP , 026 (2019), arXiv:1909.07992 [hep-ph].[16] P. A. R. Ade et al. (Planck), Astron. Astrophys. ,A16 (2014), arXiv:1303.5076 [astro-ph.CO].[17] N. Aghanim et al. (Planck), (2018), arXiv:1807.06209[astro-ph.CO].[18] N. Bernal and F. Hajkarim, Phys. Rev. D , 063502(2019), arXiv:1905.10410 [astro-ph.CO].[19] S. Bhattacharya, S. Mohanty, and P. Parashari, Phys.Rev. D , 043522 (2020), arXiv:1912.01653 [astro-ph.CO].[20] N. Bernal, A. Ghoshal, F. Hajkarim, and G. Lambiase,(2020), arXiv:2008.04959 [gr-qc].[21] B. Ratra and P. J. E. Peebles, Phys. Rev. D , 3406(1988).[22] L. J. Hall, K. Jedamzik, J. March-Russell, and S. M.West, JHEP , 080 (2010), arXiv:0911.1120 [hep-ph].[23] N. Bernal, M. Heikinheimo, T. Tenkanen, K. Tuomi-nen, and V. Vaskonen, Int. J. Mod. Phys. A32 , 1730023(2017), arXiv:1706.07442 [hep-ph].[24] J. Edsjo and P. Gondolo, Phys. Rev. D , 1879 (1997),arXiv:hep-ph/9704361.[25] B. Barman, S. Bhattacharya, and B. Grzadkowski,(2020), arXiv:2009.07438 [hep-ph].[26] A. De Simone, V. Sanz, and H. P. Sato, Phys. Rev. Lett. , 121802 (2010), arXiv:1004.1567 [hep-ph].[27] D. J. Chung, E. W. Kolb, and A. Riotto, Phys. Rev. D , 063504 (1999), arXiv:hep-ph/9809453. [28] T. Moroi, H. Murayama, and M. Yamaguchi, Phys. Lett.B , 289 (1993).[29] M. Kawasaki and T. Moroi, Prog. Theor. Phys. , 879(1995), arXiv:hep-ph/9403364.[30] L. Kofman, A. D. Linde, and A. A. Starobinsky, Phys.Rev. D , 3258 (1997), arXiv:hep-ph/9704452.[31] A. D. Linde, Particle physics and inflationary cosmology,Vol. 5 (1990) arXiv:hep-th/0503203.[32] B. Barman, D. Borah, and R. Roshan, JCAP , 021(2020), arXiv:2007.08768 [hep-ph].[33] M. Kakizaki, S. Kanemura, and T. Matsui, Phys. Rev.D , 115007 (2015), arXiv:1509.08394 [hep-ph].[34] C. L. Wainwright, Comput. Phys. Commun. , 2006(2012), arXiv:1109.4189 [hep-ph].[35] S. R. Coleman and E. J. Weinberg, Phys. Rev. D7 , 1888(1973).[36] P. Basler, M. Krause, M. Muhlleitner, J. Wittbrodt,and A. Wlotzka, JHEP , 121 (2017), arXiv:1612.04086[hep-ph].[37] P. B. Arnold and O. Espinosa, Phys. Rev. D47 ,3546 (1993), [Erratum: Phys. Rev.D50,6662(1994)],arXiv:hep-ph/9212235 [hep-ph].[38] A. Kosowsky, M. S. Turner, and R. Watkins, Phys. Rev.
D45 , 4514 (1992).[39] A. Paul, B. Banerjee, and D. Majumdar, JCAP ,062 (2019), arXiv:1908.00829 [hep-ph].[40] A. Kosowsky and M. S. Turner, Phys. Rev.
D47 , 4372(1993), arXiv:astro-ph/9211004 [astro-ph].[41] S. J. Huber and T. Konstandin, JCAP , 022 (2008),arXiv:0806.1828 [hep-ph].[42] A. Kosowsky, M. S. Turner, and R. Watkins, Phys. Rev.Lett. , 2026 (1992).[43] M. Kamionkowski, A. Kosowsky, and M. S. Turner,Phys. Rev. D49 , 2837 (1994), arXiv:astro-ph/9310044[astro-ph].[44] C. Caprini, R. Durrer, and G. Servant, Phys. Rev.
D77 ,124015 (2008), arXiv:0711.2593 [astro-ph].[45] M. Hindmarsh, S. J. Huber, K. Rummukainen, andD. J. Weir, Phys. Rev. Lett. , 041301 (2014),arXiv:1304.2433 [hep-ph].[46] J. T. Giblin, Jr. and J. B. Mertens, JHEP , 042 (2013),arXiv:1310.2948 [hep-th].[47] J. T. Giblin and J. B. Mertens, Phys. Rev. D90 , 023532(2014), arXiv:1405.4005 [astro-ph.CO].[48] M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J.Weir, Phys. Rev.
D92 , 123009 (2015), arXiv:1504.03291[astro-ph.CO].[49] C. Caprini and R. Durrer, Phys. Rev.
D74 , 063521(2006), arXiv:astro-ph/0603476 [astro-ph].[50] T. Kahniashvili, A. Kosowsky, G. Gogoberidze,and Y. Maravin, Phys. Rev.
D78 , 043003 (2008),arXiv:0806.0293 [astro-ph].[51] T. Kahniashvili, L. Campanelli, G. Gogoberidze, Y. Mar-avin, and B. Ratra, Phys. Rev.
D78 , 123006 (2008), [Er-ratum: Phys. Rev.D79,109901(2009)], arXiv:0809.1899[astro-ph]. [52] T. Kahniashvili, L. Kisslinger, and T. Stevens,Phys. Rev. D81 , 023004 (2010), arXiv:0905.0643 [astro-ph.CO].[53] C. Caprini, R. Durrer, and G. Servant, JCAP , 024(2009), arXiv:0909.0622 [astro-ph.CO].[54] B. Barman, A. Dutta Banik, and A. Paul, Phys. Rev. D , 055028 (2020), arXiv:1912.12899 [hep-ph].[55] A. Paul, U. Mukhopadhyay, and D. Majumdar, (2020),arXiv:2010.03439 [hep-ph].[56] P. J. Steinhardt, Phys. Rev.
D25 , 2074 (1982).[57] P. S. B. Dev, F. Ferrer, Y. Zhang, and Y. Zhang, JCAP , 006 (2019), arXiv:1905.00891 [hep-ph].[58] V. R. Shajiee and A. Tofighi, Eur. Phys. J.
C79 , 360(2019), arXiv:1811.09807 [hep-ph].[59] C. Caprini et al., JCAP , 001 (2016),arXiv:1512.06239 [astro-ph.CO]. [60] J. Ellis, M. Lewicki, and J. M. No, (2018), 10.1088/1475-7516/2019/04/003, [JCAP1904,003(2019)],arXiv:1809.08242 [hep-ph].[61] J. Ellis, M. Lewicki, J. M. No, and V. Vaskonen, JCAP , 024 (2019), arXiv:1903.09642 [hep-ph].[62] E. Thrane and J. D. Romano, Phys. Rev.
D88 , 124032(2013), arXiv:1310.5300 [astro-ph.IM].[63] K. Schmitz, (2020), arXiv:2002.04615 [hep-ph].[64] C. Moore, R. Cole, and C. Berry, Class. Quant. Grav. , 015014 (2015), arXiv:1408.0740 [gr-qc].[65] T. Alanne, T. Hugle, M. Platscher, and K. Schmitz,JHEP , 004 (2020), arXiv:1909.11356 [hep-ph].[66] S. Vagnozzi, (2020), arXiv:2009.13432 [astro-ph.CO].[67] A. Dutta Banik, A. K. Saha, and A. Sil, Phys. Rev. D , 075013 (2018), arXiv:1806.08080 [hep-ph].[68] P. Konar, A. Mukherjee, A. K. Saha, and S. Show,(2020), arXiv:2007.15608 [hep-ph].[69] M. Fairbairn and R. Hogan, JHEP09