Probing high scale seesaw and PBH generated dark matter via gravitational waves with multiple tilts

We propose a scenario where a high scale seesaw origin of light neutrino mass and gravitational dark matter (DM) in MeV-TeV ballpark originating from primordial black hole (PBH) evaporation can be simultaneously probed by future observations of stochastic gravitational wave (GW) background with multiple tilts or spectral breaks. A high scale breaking of an Abelian gauge symmetry ensures the dynamical origin of the seesaw scale while also leading to the formation of cosmic strings responsible for generating stochastic GW background. The requirement of a correct DM relic in this ballpark necessitates the inclusion of a diluter as PBH typically leads to DM overproduction. This leads to a second early matter dominated epoch after PBH evaporation due to the long-lived diluter. These two early matter dominated epochs, crucially connected to the DM relic, lead to multiple spectral breaks in the otherwise scale-invariant GW spectrum formed by cosmic strings. We find interesting correlations between DM mass and turning point frequencies of GW spectrum which are within reach of several near future experiments like LISA, BBO, ET, CE, etc.


I. INTRODUCTION
The observation of neutrino oscillation has been the most compelling experimental evidence suggesting the presence of beyond standard model (BSM) physics.This is because the standard model (SM) alone can not explain the origin of neutrino mass and mixing, as verified at neutrino oscillation experiments [1,2].Different BSM frameworks have, therefore, been invoked to explain non-zero neutrino mass.Conventional neutrino mass models based on seesaw mechanism [3][4][5][6][7][8][9][10][11][12][13] typically involve the introduction of heavy fields like right handed neutrinos (RHN).Typically, such canonical seesaw models have a very high seesaw scale keeping it away from the reach of any direct experimental test.This has led to some recent attempts in finding ways to probe high scale seesaw via stochastic gravitational wave (GW) observations [14][15][16][17][18][19][20][21][22][23].In these scenarios, the presence of additional symmetries was considered whose spontaneous breaking not only leads to the dynamical generation of RHN masses but also leads to the generation of stochastic GW from topological defects like cosmic strings [14][15][16][17], domain walls [18] or from bubbles generated at first order phase transition [19][20][21][22][23].
In addition to the origin of neutrino mass, the SM also fails to explain the origin of dark matter (DM), a mysterious, non-luminous and non-baryonic form of matter comprising around 27% of the present universe.The present DM abundance is often reported in terms of density parameter Ω DM and reduced Hubble constant h = Hubble Parameter/(100 km s −1 Mpc −1 ) as [24] Ω DM h 2 = 0.120 ± 0.001 (1) at 68% CL.Among a variety of BSM proposals to explain the origin of DM, the weakly interacting massive particle (WIMP) has been the most well-studied one.In the WIMP framework, a particle having mass and interactions around the electroweak ballpark can lead to the observed DM abundance after undergoing thermal freeze-out in the early universe [25].While WIMP has promising direct detection signatures, contrary to high scale seesaw, the null results at direct detection experiments [26] have motivated us to consider other alternatives.Such alternative DM scenarios may have suppressed direct detection signatures but can leave signatures in future GW experiments.There have been several recent works on such GW signatures of DM scenarios [27][28][29][30][31][32][33][34][35][36].
Motivated by this, we consider a high scale seesaw scenario with gravitational DM gen-erated from evaporating primordial black holes (PBH) assuming the latter to dominate the energy density of the universe at some stage.While the seesaw scale is dynamically generated due to the spontaneous breaking of a U (1) gauge symmetry at a high scale, DM can be generated during a PBH-dominated epoch.Unlike sub-GeV or superheavy DM (with mass ≥ O(10 9 GeV)) considered in most of the PBH-generated DM scenarios like [37][38][39][40][41][42][43], we consider DM around MeV-TeV ballpark [44].Since DM in this ballpark gets overproduced during a PBH-dominated era, we also consider the presence of a long-lived diluter whose late decay brings the DM abundance within limits.The diluter is also dominantly produced from PBH as its coupling with the SM remains suppressed in order to keep it long-lived.
The high-scale breaking of U (1) gauge symmetry leads to the formation of cosmic strings (CS) [45,46], which generate stochastic GW background with a characteristic spectrum within the reach of near future GW detectors if the scale of symmetry breaking is sufficiently high [47,48].In the present setup, we show that the PBH domination followed by the subsequent diluter domination era lead to multiple spectral breaks in the GW spectrum formed by CS network.In a recent work, [32] the correlation between such spectral break due to non-standard cosmological era [49][50][51] of diluter domination and DM properties were studied within the framework of a realistic particle physics model with high scale seesaw.
Here we extend this to a PBH-generated gravitational DM scenario with the unique feature of multiple spectral breaks in scale-invariant GW background generated by CS, not considered in earlier works.Depending upon DM mass, the three spectral breaks or turning point frequencies can be within reach of different future GW experiments.The novel feature of our study is in finding the role of multiple early matter domination eras in the otherwise scale-invariant GW spectrum of CS network leading to multiple turning point frequencies and the connection of these eras as well as the turning point frequencies of GW spectrum with mass of gravitational DM produced from PBH evaporation.
This paper is organized as follows.In section II, we briefly discuss our setup followed by a discussion of dark matter relic from PBH in section III.In section IV, we discuss the details of the gravitational wave spectrum generated by cosmic strings and subsequent spectral breaks due to PBH and diluter domination for different DM masses.In section V, we discuss briefly the gravitational waves generated from PBH density fluctuations.Finally, we conclude in section VI.

II. THE FRAMEWORK
The present framework can be embedded in any Abelian extension of the SM that guarantees the formation of cosmic strings due to the spontaneous symmetry breaking with an era of primordial black hole domination.A singlet scalar with no coupling to the SM or to the new Abelian gauge sector is assumed to be the DM having purely gravitational interactions.A diluter in the form of a heavy singlet Majorana fermion or a heavy RHN is also incorporated while keeping its coupling to the SM leptons tiny.In order to ensure that the PBH domination and diluter domination arise as separate epochs for simplicity, we consider the diluter to be singlet under the new Abelian gauge symmetry as well.
Here, we outline the minimal interaction terms to realise our scenario with one possible UV completion being given in Appendix A. The two RHNs having sizeable couplings to leptons thereby generating light neutrino masses via type I seesaw mechanism1 [3-5, 7, 8, 12] are denoted by N 1,2 while the one playing the role of diluter and having tiny Yukawa coupling is denoted by N 3 .Dark matter is assumed to be a real singlet scalar with only gravitational interactions.The RHNs N 1,2 are charged under a U (1) gauge symmetry such that the spontaneous U (1) symmetry breaking by a singlet scalar S dynamically generates the seesaw scale.The relevant Lagrangian of the newly introduced fields is given by In order to ensure N 3 to be long-lived Y α3 ≪ 1. Due to tiny coupling with the bath particles, the diluter is dominantly produced only from PBH evaporation.Similarly, due to the absence of any DM coupling with the bath, it will also be solely produced from PBH.Additionally, the singlet scalar S can acquire a non-zero vacuum expectation value (VEV), leading to the generation of N 1,2 as well as U (1) gauge boson masses.The same symmetry breaking also leads to the formation of cosmic strings.We discuss the details of DM relic generation and cosmic string dynamics leading to gravitational waves in the upcoming sections.

III. DARK MATTER AND DILUTER FROM PBH
As mentioned above, the DM in the present framework only has gravitational interactions and hence can be dominantly produced from PBH evaporation.On the other hand, as a result of its feeble interaction with the SM bath, N 3 can also be assumed to be produced purely from the PBH evaporation.While the PBH itself can act as a DM if its mass is chosen appropriately, in this work we focus on the ultra-light mass regime of PBH where it is not cosmologically long-lived.As is well known from several earlier works, the DM in a mass range: a few keV to around 10 10 GeV, when produced from the PBH evaporation are always overabundant if PBHs dominate the total energy density of the Universe [31,52].
A recent study [44] has shown that the overproduction of such DM can be prevented if one considers a diluter that injects entropy at a late time into the thermal bath as a result of its slow decay.Following [44], in the present work we also consider such multiple early matter-dominated eras and show how their presence can affect the GW spectrum generated by cosmic strings in a unique way.
PBH is assumed to be formed in the era of radiation domination after inflation.In such a situation, the PBH abundance is characterized by a dimensionless parameter β defined as where ρ BH (T in ) and ρ R (T in ) represent the initial PBH energy density and radiation energy density respectively while T in denotes the temperature at the time of PBH formation.The expression of T in can be found in appendix B. Once PBHs are formed, the radiation-dominated universe eventually turns into a matter-dominated universe and remains matter-dominated till the epoch of PBH evaporation.The temperature corresponding to the PBH evaporation is denoted as T ev and can be obtained from Eq. (B8).The condition determining PBH's evaporation during the radiation domination is given as [53], where M pl indicates the Planck mass and with ω i = 2 s i + 1 for massive particles of spin s i , ω i = 2 for massless species with s i > 0 and ω i = 1 for s i = 0. β s = 2.66, 4.53, 6.04, 9.56 for s i = 0,1/2,1,2 respectively, whereas , where m i is the mass of the i th species [53].β c is the critical PBH abundance that leads to an early matter-dominated era.In our scenario, we consider β to be large enough such that PBH dominates the energy density of the universe at some epoch.In the above inequality, γ ≃ 0.2 is a numerical factor that contains the uncertainty of the PBH formation, and G ∼ 4 is the grey-body factor.Here, M pl denotes the Planck mass while m in and T BH denote the mass of PBH at the time of formation and instantaneous Hawking temperature of PBH (see appendix B for the details).Additionally, we also assume the PBHs to be of Schwarzschild type without any spin and charge and having a monochromatic mass spectrum implying all PBH have identical masses.It is important to mention that a PBH can dominantly evaporate only into particles lighter than its instantaneous Hawking temperature.
An upper bound on the PBH mass is obtained by demanding its evaporation temperature T ev > T BBN ≃ 4 MeV, as the evaporation of PBH also results in the production of the radiation that can disturb the successful prediction of the big bang nucleosynthesis (BBN).
While the upper bound on the PBH mass can be obtained by comparing its evaporation temperature with the BBN temperature, a lower bound on its mass can be obtained from the cosmic microwave background (CMB) bound on the scale of inflation, i.e.H(T in ) ≤ H I ≤ 2.5 × 10 −5 M pl , where H(T in ) = 1  2 t in with t(T in ) = m in M 2 pl γ (as obtained from Eq. (B1)).Using these BBN and CMB bounds together, we have a window for allowed initial mass for PBH that reads 0.1 g ≲ m in ≲ 3.4 × 10 8 g.The range of PBH masses between these bounds is unconstrained from BBN, CMB and other such bounds at low redshifts [54].However, Refs.[55,56] have recently put constraints on such low PBH mass as well, from the gravitational wave background produced from PBH density fluctuations, since these GW can lead to extra relativistic degrees of freedom during BBN.We discuss this GW spectrum in Section V.
Since in this work, our primary goal is to investigate the effect of a PBH-dominated Universe on the production of dark matter and its signature on gravitational waves, we remain agnostic about the formation of PBH.However, our choice of the initial energy density of PBH β, which we consider to be a free parameter satisfying β > β c ensures that we have a PBH domination.Now, the production of PBH with the mass and β ranges we have considered can take place, say from perturbations generated during inflation, when they re-enter the horizon, or from phase transitions.For inflationary production mechanisms, the mass of PBH formed can be related to the scale of the perturbation as [57][58][59][60][61] where M ⊙ denotes the solar mass and g * (T k ) indicates the relativistic degrees of freedom when the scale 'k' enters the horizon.Hence, for forming PBH with a mass function peaked around m in ∼ 10 3 −10 6 g, which, as we will see, is the range we have used, an enhancement to the power spectrum is needed at extremely small scales k ∼ 10 19 −10 21 Mpc −1 , corresponding to the last stages of inflation.We leave inflationary potentials giving rise to such features in the power spectrum for future studies.Now, the initial abundance of PBH β is determined by the amplitude of the power spectrum P R [60,62], which we can roughly estimate to be given by P R ∼ 0.2 ln(1/β) [60].The amplitude is thus only logarithmically dependent on β and we can have our required β > β c with P R ∼ 10 −2 .PBH can also be formed during a first order phase transition (FOPT) [63][64][65][66][67][68][69] a result of bubble collisions, Fermi-ball collapse among several others.For example, in the collapsing Fermi-ball scenario [68], a dark sector fermion χ with an asymmetry between χ and χ leads to the formation of Fermi-balls after getting trapped in the false vacuum.Due to the strong Yukawa force among the dark fermions, such Fermi-balls can collapse into PBH.The energy fraction for such PBH is given by [68] with T * , v b , β PT /H being the nucleation temperature, bubble wall velocity and time scale ratio of universe's expansion and the FOPT respectively.For PBH mass of our interest, say m in = 10 7 g, we can have The initial PBH mass of 10 7 g can be obtained for the same choice of parameters if dark sector asymmetry is ∼ O(10 −4 ).
There are other PBH formation mechanisms discussed in the literature namely, from the collapse of topological defects [70,71] which we do not elaborate here.Thus, it is possible to generate PBH with the desired initial mass and energy fraction if we incorporate additional physics as mentioned above.In the spirit of minimality, we stick to our minimal setup to highlight the key feature of our proposal namely, the effect of multiple early matter domination on cosmic string generated GW spectrum with a non-trivial connection to dark matter genesis from evaporating PBH.

A. Particle production from PBH
Once formed, the PBH can evaporate by Hawking radiation [72,73].Among the products of PBH evaporation, there might be a stable BSM particle that can contribute to the observed DM abundance [31, 37-44, 53, 74-89] or there might be a new particle produced which is responsible for generating the present matter-antimatter asymmetry of the universe [38,44,52,72,[90][91][92][93][94][95][96][97][98][99][100][101][102].The number of any particle X produced during the evaporation of a single PBH can be estimated by considering the differential mass decrease , where T BH is the Hawking temperature and dE is the corresponding energy emitted.Considering the mean energy of the emitted particles as 3T BH , the differential number of emitted particles is found to be Integrating the above expression gives us the number of particles X as Here, T in BH denotes the instantaneous PBH temperature at the time of the formation.While g ⋆,H accounts for all emitted particles (Eq.5), g X,H corresponds to the particle X.Note that for m X < T in BH , the emission of particles takes place from the beginning of PBH formation, while in the opposite case of m X > T in BH , the emission occurs only when T BH reaches m X .This leads to different expressions in these two limits after integration.The yield of a particle produced during the evaporation of a PBH is related to the PBH abundance (n BH ) and is expressed as, where s denotes the entropy density at evaporation.PBH abundance at the time of its evaporation can be obtained using the first Friedmann equation as where τ is the PBH lifetime (cf.Eq. (B7)).Substituting Eq. ( 11) in Eq. ( 10), one can easily obtain the asymptotic yield of the particles produced during the PBH evaporation.
Following this, one can easily calculate the abundance of the DM and the diluter (N 3 ) produced during the evaporation of the PBH.Once the number of DM particles produced from PBH evaporation is known, the DM relic abundance Ω DM h 2 at the present epoch can be calculated as with 10640 π .Here ζ parametrizes a possible entropy production after PBH evaporation until now, i.e., ζ (sa 3 ) ev = (sa 3 ) 0 .Looking at the above equation one finds that the heavier the DM mass (for m DM > T in BH ), the lesser will be the DM relic abundance.On contrary, in a scenario with m DM < T in BH , the DM relic abundance becomes proportional to the DM mass.The parameter space for DM produced during PBH evaporation is highly constrained.The correct relic abundance is satisfied only if m DM ≳ 10 10 GeV and m in ≳ 10 6 g [31,52].In the rest of the parameter space, the DM is mostly overabundant if one assumes that there exists no entropy injection into the thermal plasma after the evaporation of the PBH, i.e. ζ = 1.This tension can be relaxed considerably if one assumes ζ ∼ O (10) or O (100).While DM of mass around keV or lighter do not get overproduced, it suffers from strong Lyman-α bounds [40,103,104] severely restricting the parameter space in β − m in plane.We do not consider such light or superheavy DM in our setup and focus on the intermediate mass range which typically gets overproduced during a PBH dominated phase.
Analogous to the DM, the diluter N 3 is also produced during the PBH's evaporation whose abundance can be calculated using Eq.(10).Once produced, it decays very slowly as a result of its very feeble interaction with the first generation of the SM lepton doublet and the Higgs.Subsequent to the PBH evaporation, the universe enters into a radiationdominated era for a short duration following which the energy density of N 3 overtakes the radiation energy density and the universe re-enters the matter-dominated phase in a similar manner as was shown in [44].The universe experiences a change in temperature evolution and cools slower as the diluter starts to decay.This is because, in the process of its decay, the diluter also injects entropy into the thermal bath.Finally, after the decay of the N 3 is complete the universe again becomes radiation-dominated.A schematic of the scenario is shown in Fig. 1 and the evolution of the energy densities of the different components is shown in Fig. 2 for a particular choice of benchmark parameters.As a result of the entropy injection, the number density of DM is brought within observed limits thereby alleviating the tension of DM relic abundance produced from PBH evaporation in the intermediate mass range mentioned before.
It should be noted that in addition to PBH evaporation, DM and diluter can also be produced from 2 to 2 scatterings mediated by massless graviton.Adopting the approach discussed in earlier works [42,44,[105][106][107][108][109], we find that such productions are negligible for the choices of m in , M DM and N 3 mass M 3 considered in our numerical analysis to be discussed below.
Finally, we would like to mention here that one can also have leptogenesis from the decay of the RHNs produced from PBH evaporation.However, as studied in details in Refs.[44,96,97,99], the baryon asymmetry produced for the PBH mass range considered in our analysis, the asymmetry turns out to be smaller than the observed value.As we will see in the next section, decreasing the PBH mass is not favorable from the gravitational wave detection perspective as it would increase the turning point frequency to a higher value beyond the reach of near-future detectors and would also result in particle production from cosmic strings to be dominant compared to gravitational wave emission.Moreover, the asymmetry that might be produced thermally at a high temperature would be diluted because of two epochs of entropy dilution arising from PBH evaporation and N 3 decay.

IV. GRAVITATIONAL WAVES FROM COSMIC STRINGS
Cosmic strings [45,46], are one-dimensional objects which appear as topological defects after the spontaneous breaking of a symmetry group containing a vacuum manifold that is not simply connected [46].The simplest group that exhibits such a feature is U (1) which naturally appears in many BSM frameworks.Numerical simulations [110,111] based on Nambu-Goto action indicate that dominant energy loss from a string loop is in the form of GW radiation if the underlying symmetry is gauged.Thus, being one of the potential sources of primordial GW, they have gained a great deal of attention after the recent finding of a stochastic common spectrum process across many pulsars [17,[112][113][114][115].If the symmetry breaking scale (Λ CS ) is high enough (Λ CS ≳ 10 9 GeV), the resulting GW background is detectable.This makes CS an outstanding probe of super-high scale physics [14,16,32,[116][117][118][119][120][121][122][123][124][125].The properties of CS are described by their normalised tension Gµ ∼ GΛ 2  CS with G being the Newton's constant.Unless the motion of a long-string network gets damped by thermal friction [126], shortly after formation, the network oscillates (t osc ) and enters the scaling regime [111,127,128] which is an attractor solution of two competing dynamicsstretching of the long-string correlation length due to cosmic expansion and fragmentation of the long strings into loops which oscillate to produce particle radiation or GW [47,48,129].
A set of normal-mode oscillations with frequencies f k = 2k/l constitute the total energy loss from a loop, where the mode numbers k = 1, 2, 3....∞2 .Therefore, the GW energy density parameter is defined as Ω GW (t 0 , f ) = k Ω (k) GW (t 0 , f ), with t 0 being the present time and f ≡ f (t 0 ) = f k a(t 0 )/a(t).Present day GW energy density corresponding to the mode k is computed with the integral [130] where n (t, l k ) is a scaling loop number density which can be computed analytically using Velocity-dependent-One-Scale (VOS) [131][132][133] model3 , ρ c is the critical energy density of the universe and Γ k = Γk −δ ζ(δ) depends on the small scale structures in the loops such as cusps (δ = 4/3) and kinks (δ = 5/3).In this article, we consider only cusps to compute GW spectrum.A similar analysis can also be done considering kinks.Although for higher modes, contributions to the GWs from cusps are dominant, the number of cusps and kinks per loop oscillation cannot be straightforwardly determined with numerical simulations, which do not include gravitational wave back-reaction in general.A preference for cusps over kinks comes from the so-called smoothing mechanism, where the kinky loops are expected to be smoothed out due to the gravitational wave back-reaction [134].This mechanism, however, has been challenged in [135,136].Therefore, we would like to place the consideration of cusps over kinks in this article as a choice rather than a preference.
A typical feature of GWs from CS is a flat plateau due to loop formation and decay during radiation domination, with an amplitude given by where ϵ r = α/ΓGµ with α the initial (at t = t i ) loop size parameter, Ω r ≃ 9 × 10 −5 and A r = 5.4 [133].In our analysis, we have considered α ≃ 0.1 [130,137] (as suggested by simulations), Γ ≃ 50 [129,130,137], and CMB constraint Gµ ≲ 10 −7 [138] which lead to α ≫ ΓGµ.In this limit, Eq.( 14) implies Ω (k=1) GW (f ) ∼ Λ CS , a property that makes models with larger breaking scales more testable with GWs from CS.The observed frequency today of the GW spectrum is connected to the temperature at which the loops responsible for that particular frequency have been formed.The major contribution to the GW emission is at a time when the loop reaches its half-life t E = t H ≃ αt i 2ΓGµ , where t i indicates the time of loop formation.Now, since the frequency emitted at the time t E is related to the length of the loop as (f (t E ) = 2k/l(t E )), we get (considering k=1) where t eq (z eq ), t 0 indicate the standard matter-radiation equality time (redshift) and the present time respectively.
The dependence of the GW amplitude on the frequency can be qualitatively understood as follows [51].The frequency dependence receives contribution from two factors.Because of the Universe's expansion, the GW energy density dilutes as a −4 .Hence, the GW amplitude emitted by loops formed at higher temperatures and subsequently at higher frequencies undergoes more suppression until the present time.However, at earlier times the loop production rate, which varies as ∝ t −4 i [49,130] also increases, leading to enhancement of the GW amplitude.Thus, the total dependence on frequency is the effect of these two factors.In the case of radiation domination, these two contributions exactly cancel each other, leading to a flat spectrum.If there exist cosmological epochs other than radiation at higher temperatures, there would be a deviation from flatness, which would first start to appear because of the contribution of those loops that have been formed at a time say t ∆ (with corresponding temperature T ∆ ), when the standard radiation era began or the nonstandard era ended.Thus, the frequency at which the spectral break occurs or the turning point frequency f ∆ is given by Eq. ( 15) with t i = t ∆ [49][50][51], which in terms of temperature gives us Beyond f ∆ , the spectrum goes as Ω GW ∼ f −1 for k = 1 mode, which can be estimated by considering the corresponding change in the loop number density and the cosmic evolution because of the matter-dominated background (Please refer to [51] for a detailed derivation).Now, as shown in Fig. 1, in our scenario we have two phases of early matter domination, one because of N 3 domination just before BBN (MD2) and the other because of PBH (MD1) at much earlier epochs.This leads to some interesting patterns in the GW spectrum as we will show below.First, let us discuss the impact of the recent matter domination because of N 3 (MD2), which ends at temperature T 5 .Note that the lifetime of N 3 should be large enough such that it does not decay before leading to matter-domination.We attempt to find an analytical estimate for T 5 in terms of the other parameters in our scenario like M DM , m in , M 3 such that the turning point frequency f ∆ given by equation ( 16) (with T ∆ = T 5 ), can be related to these parameters.For this purpose, we need to calculate the entropy dilution arising because of N 3 decay, given by S = s(T 5 )a(T 5 ) 3 s(T 4 )a(T 4 ) 3 , which can be estimated to be [25,139,140] where Y 3 = n N 3 /s is the initial yield of N 3 arising from PBH and Γ 3 is its decay width.
This should be equal to = Ω DM h 2 /0.12 to produce the observed DM abundance.T N 3 ≃ T 5 represents the temperature at the end of N 3 decay and assuming instantaneous decay of N 3 , we write Note that Y N 3 , in turn depends on m in , M 3 (cf.Eq. ( 10)).Using Eq. ( 10) for Y N 3 , Eq. ( 12) for Ω DM h 2 , Eq. ( 18) for Γ 3 and substituting in Eq. ( 17), we find For both the cases above, T 5 decreases with a higher value of dark matter mass M DM .This is expected since the dark matter relic from PBH (M DM < T in BH ) increases with M DM and hence a longer duration of N 3 domination would be required to produce the observed relic, which decreases T 5 .Now, for the case M 3 > T in BH , increasing the initial PBH mass decreases the DM relic (see Eq. ( 12)), and also the N 3 yield, with the latter effect being more dominant.Hence, N 3 starts dominating much later which results in a smaller T 5 .For the other case M 3 < T in BH , the dependence on m in does not appear because of a similar dependence of the N 3 yield on m in as that of DM.Finally, increasing the diluter mass M 3 decreases its yield for M 3 > T in BH and leads to a decrease in T 5 .The reverse is true for the case M 3 < T in BH .In our analysis, we consider M 3 > T in BH , and fix its value at M 3 = 10 11 GeV4 .Thus, the turning point frequency (Eq.( 16)) varies as 12.The experimental sensitivities of SKA [141], GAIA [142], EPTA [143], THEIA [142], µARES [144], LISA [145], DECIGO [146], BBO [147], ET [148], CE [149] and aLIGO [150]  Next, we look into the effects of the matter-dominated era arising from PBH (MD1).
Here, the expression for the turning point frequency f ∆ given by Eq. ( 16) is changed because of the other intermediate epochs before BBN.This can be understood as follows.
A loop formed at the time t 3 would contribute to the spectrum much later when it reaches its half-life which is given by [51] tM The frequency emitted is related to the length of the loop at its half-life time f ( tM ) = 2k/l( tM ) and hence we get where f is the frequency observed today.Now, usually, tM is in the radiation-dominated era after BBN, which gives us Eq. ( 16).However, in our case, tM would be in the radiation or early matter-dominated era before BBN, which gives us 5 Substituting tM from Eq. ( 20), we find the turning point frequency, which we call f ′ ∆ to be given by , we get Thus, the TPF f ′ ∆ depends not only on the PBH evaporation temperature T 3 , but also on T 5 , i.e, the temperature when N 3 dominated era ends.Using Eq. ( 19) and B8, we find Similar turning point frequencies can also be obtained for loops below a critical length (ΓGµ) 2 [151], for which particle production becomes dominant over the GW emission.This translates into the bound Gµ > 2.4×10 −16 T 4/5 ∆ .In our scenario, we choose Gµ = 10 −11 which is high enough such that the GW emission always dominates.This gives an upper bound on the TPF as T ∆ ≲ 5 × 10 5 GeV, and considering T ∆ = T 3 = T ev , we get a lower bound on the initial PBH mass as m in ≳ 2 × 10 3 g. 5For most of our parameters used, the radiation-dominated era after PBH domination (RD2) lasts for a very small duration.Hence, tM is in the N 3 dominated era.However, for large PBH masses, the radiation era is prolonged and loops reach their half-life tM in the radiation era itself, which modifies Eq. ( 22) as . (right panel), along with the sensitivities of future experiments like LISA [145], BBO [147], ET [148], CE [149], µARES [144].Here, Gµ = 10 −11 and M 3 = 10 11 GeV.
We now compute the GW spectra by integrating Eq. ( 13) over the different cosmological epochs.All the temperatures (T 2 , T 3 , T 4 , T 5 ) corresponding to the different cosmological backgrounds (see Fig. 1) which would be required in this integration are found by numerically solving the following coupled Boltzmann equations [95,97,99] for ρ BH , ρ R , ρ DM and Here, ∆ = 1 + T 3g * s(T ) dT .ρ R and ρ BH denote the comoving energy densities of radiation and PBH respectively, H is the Hubble parameter, whereas n N 3 and n DM represent the comoving number densities of N 3 and DM respectively.Note that the terms on the R.H.S of ρ R evolution equation denote entropy dilution from PBH evaporation and N 3 decay, with the corresponding second and third terms in the temperature evolution equation.Since we are considering a scenario of DM with only gravitational interactions, it gets produced only from PBH.Similarly, production of N 3 from the thermal bath also remains suppressed because of its tiny Yukawa couplings required for a longer lifetime.Γ BH→N 3 , Γ BH→DM denote the production rate from PBH [95], whereas Γ 3 indicates the decay rate of N 3 .
In Fig. 3, we show the GW spectra for three DM masses (left panel) and three PBH masses (right panel), differing by orders of magnitude.The GW spectral behavior can be connected to the different cosmic regimes represented in Fig. 1.We show this for one of the BPs in the right panel of Fig. 3. Here, the label indicates the contribution to the GW spectrum from loops emitting in different cosmic eras.We refrain from showing the contribution coming from RD1 since it occurs at ultra high frequencies with sub-dominant amplitude out of reach of near-future detectors.The decay rate of N 3 is chosen such that Ω DM h 2 = 0.12.As anticipated above, the spectral breaks due to the two early matterdominated eras can be clearly seen.The first break (f ∆ ) at the lower frequency, from flat to power-law, depends on the temperature T 5 which marks the end of N 3 domination.It is evident here that f ∆ decreases with DM mass M DM and initial PBH mass m in .The break from flat to power-law at higher frequencies (f ′ ∆ ) is related to PBH evaporation at the temperature T 3 , and a higher evaporation temperature (lower m in ) implies break at a higher frequency (right panel).However, as discussed above, this break also depends on the other intermediate eras before BBN, specifically on the temperature T 5 (see Eq. ( 24)).This is apparent from the left panel where even for the same initial PBH mass m in (and hence the same T 3 ), the break occurs at different frequencies.Note that there appears another break in the spectrum in the middle, from power-law to flat, which can also be probed in the detectors.
In the left panel of Fig. 4, we show the contours in the M DM − f ∆ plane (grey-shaded) which satisfies the correct DM relic, along with the sensitivities of the GW experiments which can probe them.We find our numerical results to agree with our analytical estimates derived earlier up to a difference of around O(1).As expected, the contours move up when we increase the PBH mass m in .The upper bound on the DM mass appears since N 3 has to decay before BBN.Note that this bound changes with m in , since a heavier PBH evaporates at a lower temperature, closer to BBN.Decreasing DM mass leads to Ω DM h 2 < 0.12, which gives us a lower bound on DM mass.In the right panel of Fig. 4, for the same PBH masses we show the corresponding predictions at the higher frequencies, i.e. the M DM − f ′ ∆ plane.As we can see, for lighter PBH and heavier DM, f ′ ∆ falls outside the sensitivities.Finally, in Table I, we provide some benchmark values of DM mass M DM and PBH mass m in along with the GW experiments which can probe the corresponding spectral breaks.Finally, recall that we are using only the k = 1 mode for computing the GW spectra.
Although summing upto high modes is computationally heavy and beyond the scope of this work, here we would like to mention the expected change and the consequences.From Eq. ( 13) one can see that the spectra for the k th mode is related to the fundamental k = 1 mode through [15,51] For the part of the spectra behaving as Ω GW ∝ f 0 , the total spectrum just gets re-scaled as Ω GW (f ) ≃ ζ(δ)Ω GW ∝ f −1 , summing over all modes changes this slope to f −1/3 [51].This effect of the higher k modes is illustrated in Fig. 5. Thus, we can conclude that including the effects of the higher k modes keeps the turning point frequencies nearly unchanged, while leading to an increase in the overall amplitude of the GW spectra, and hence still being in favor of being detected by the various GW experiments.I. Some benchmark points (BP) of M DM and m in along with the GW experiments which can probe the 1 st TPF(f ∆ ) from flat to power-law, 2 nd TPF from power-law to flat and 3 rd TPF(f ′ ∆ ) from flat to power-law.

V. GRAVITATIONAL WAVES FROM PBH
Apart from the GW background generated from Cosmic Strings, PBH themselves might be involved in the production of GW in several ways.The evaporation of PBH can produce gravitons which might constitute an ultra-high frequency GW spectrum [152].PBH can also form mergers which can lead to GW emission [153].Next, the scalar perturbations leading to the formation of PBH can induce GWs at second-order [154], which can be enhanced during PBH evaporation [155].Finally, the inhomogeneity in the distribution of PBH may also induce GW at second order, as recently studied in Ref. [55,56,156].We concentrate on this last possibility since such GW generated are independent of the formation history of PBH and especially, as we will see, for the PBH mass range we are considering, it can lead to GW signals within the sensitivity of near future GW detectors.
The distribution of PBH after they are formed are random and follow Poissonian statistics [156].These inhomogeneities induce GW at second order, when PBH dominate the energy density of the Universe, and are enhanced during PBH evaporation [55].The dominant contribution to the GW amplitude observed today can be written as [55,157,158] where Here the factor S (c.f.Eq. ( 17)) accounts for any entropy injection into the SM bath after the PBH evanescence, which in our case is provided by the decay of the long-lived diluter field.The GW spectrum has an ultraviolet cutoff at frequencies corresponding to comoving scales representing the mean separation between PBH [55,156], which is given by Note that the peak of the GW amplitude (Eq.( 28)) increases with β, which has been a free parameter in our analysis so far.In Fig. 6, we show the GW spectra generated from PBH density fluctuations for the benchmark points given in Table I, considering β = 10 −8.5 .
We ensure that β > β c holds for all the BPs, which is required for our analysis to be valid.
Moreover, we also ensure that , where β max is the upper bound on β which comes from the contribution of this GW spectrum to extra relativistic degrees of freedom during BBN [55,56].It is pertinent to note that so far our analysis of dark matter and GW spectrum from cosmic strings has been independent of the parameter β, provided it satisfies β > β c .Hence, in the study of GW from PBH density fluctuations, β appears to be a free parameter, which we consider to be small such that the GW peak amplitude which varies as β 16/3 (cf.Eq. ( 28)) remains subdominant compared to that from cosmic strings.The motivation for particularly focusing on the GW spectrum from cosmic strings is that it has distinct signatures of multiple matter dominated eras because of the presence of multiple turning point frequencies, which is usually absent for GW from other sources.Also, as discussed briefly in section III, it is indeed possible to generate PBH with such initial energy fraction and mass within some realistic scenarios.The most general scalar potential involving the different scalars can be expressed as, The scalar S breaks the U (1) Lµ−Lτ symmetry once its CP even component develops a non-zero VEV v µτ .This breaking also results in an additional non-zero mixing between the RHNs as can be seen from Eq. (A2).Once the Electroweak Symmetry is broken, the Higgs doublet (H) also develops a non-zero vev v = 246 GeV.As a result of electroweak symmetry breaking, the different scalars will mix with each other.We will not go into the details of the scalar mixing as it remains unimportant for our analysis and refer the reader to [162] for the details.In fact, due to the restrictive nature of Dirac neutrino and RHN mass matrices, fitting light neutrino data requires another singlet scalar, as pointed out in [162,163].The details of these extensions do not affect our analysis and hence we skip them here.

FIG. 1 .
FIG. 1.A schematic of the dominant phases of evolution of the universe at different epochs from a reheating temperature T 1 to the present temperature T 0 (neglecting recent dark energy domination).T 3 is the black hole evaporation temperature (equation (B8)) and T 6 represents the standard matter-radiation equality temperature of ∼ 0.75 eV.

FIG. 3 .
FIG. 3. Gravitation wave spectra for the fundamental (k = 1) mode for different values of DM mass (left panel) and initial PBH mass (right panel), satisfying the condition Ω DM h 2 = 0.12.The are shown as shaded regions of different colours.The different cosmic regimes of Fig. 1 are indicated for one of the BPs in the right panel (see text).

FIG. 5 .
FIG.5.Effect of summing over the high k modes on the GW spectrum.

Finally, we consider 1 τ 2 ∼ 2 pl 16 π 2 1/ 4 .√ 3 T
scalar singlet dark matter ϕ to have only gravitational interactions and hence it has a bare mass term (denoted by m DM ) only in the Lagrangian.Another heavy RHN, similar to N e with vanishing U (1) Lµ−Lτ charge but tiny couplings to leptons is considered to be negligible.This diluter RHN, to be denoted as N 3 , also couples more feebly to other RHNs so that N 3 dominantly decays into L e , H at late epochs.While N e could, in principle, play the role of diluter as well, it will lead to difficulties in generating correct light neutrino data due to restricted mass matrix structures.Due to tiny couplings of N 3 with other particles in the thermal bath, it is produced only non-thermally from the evaporation of the PBH.The consequence of this slow decay is a N 3 dominated epoch after the evaporation of PBH.as the PBH lifetime.One can calculate the evaporation temperature by considering H(T ev ) ∼ ρ rad (T ev ), T ev ≡ 45 M 3 g ⋆ (T ev ) τ the PBH component dominates at some point the total energy density of the universe, the SM temperature just after the complete evaporation of PBHs is: T ev = 2/ ev[40].For PBH to dominate the energy density of the Universe before they evaporate, one should have the following condition on β[53] β > β c = γ −1/2 G g ⋆,H (T BH )