While topologically stable cosmic strings are disfavored by the recent observation of nHz stochastic gravitational waves GW by Pulsar Timing Arrays PTA e.g. NANOGrav cosmic metastable strings and superstrings are not. However because the gravitational waves from all classes of strings generally span a wide range of frequencies they contradict LIGOs nonobservation of stochastic gravitational waves at the fsim  25 Hz band for a substantial stringparameter space favoured by the PTA data. Suppose ultralight primordial black holes MBHlt109 g existed in the early universe. In this case they reduce the amplitude of the GWs at higher frequencies by providing an early matterdominated phase alleviating the tension between LIGO observation and PTA data. We show that the recent PTA data complemented by future LIGOVirgoKAGRA LVK runs plus detectors such as LISA and ET would be able to dapple the properties and further search strategies of such ultralight primordial black holes which are otherwise fairly elusive as they evaporate in the early universe by Hawking radiation.
['General Relativity and Quantum Cosmology', 'Astrophysics - Cosmology and Nongalactic Astrophysics', 'High Energy Physics - Phenomenology']
Cosmic superstrings metastable strings and ultralight primordial black holes from NANOGrav to LIGO and beyond
{'Cosmic superstrings, metastable strings and ultralight primordial black holes: from NANOGrav to LIGO and beyond': "Satyabrata Datta 1, a and Rome Samanta 2, 3, b  \n1 Department of Physics and Institute of Theoretical Physics,  \nNanjing Normal University, Nanjing, 210023, China  \n2 Scuola Superiore Meridionale, Universit'a degli studi di Napoli  \n'Federico II', Largo San Marcellino 10, 80138 Napoli, Italy  \n3 INFN - Sezione di Napoli, Complesso Universitario Monte S. Angelo, I-80126 Napoli, Italy  \nWhile topologically stable cosmic strings are disfavored by the recent observation of nHz stochastic gravitational waves (GW) by Pulsar Timing Arrays (PTA), e.g., NANOGrav, cosmic metastable strings and superstrings are not. However, because the gravitational waves from all classes of strings generally span a wide range of frequencies, they contradict LIGO's non-observation of stochastic gravitational waves at the f ∼ 25 Hz band for a substantial string-parameter space favoured by the PTA data. Suppose ultralight primordial black holes ( M BH < 10 9 g) existed in the early universe. In this case, they reduce the amplitude of the GWs at higher frequencies by providing an early matter-dominated phase, alleviating the tension between LIGO observation and PTA data. We show that the recent PTA data complemented by future LIGO-Virgo-KAGRA (LVK) runs plus detectors such as LISA and ET would be able to dapple the properties and further search strategies of such ultralight primordial black holes which are otherwise fairly elusive as they evaporate in the early universe by Hawking radiation.  \nb  \nsamanta@na.infn.it", 'I. INTRODUCTION': "Ultralight primordial black holes PBHs [1-3] with initial mass M BH /uni2272 10 8 g are elusive objects as they evaporate in the early universe, therefore, they can not be constrained by astrophysical means like heavier black holes [4, 5]. Nonetheless, if they existed abundantly and dominated the energy density in the early universe, plenty of direct or indirect observable signatures are achievable in the form of gravitational radiation. For instance, they emit gravitons to constitute a high-frequency gravitational wave background [6-9] and produce density fluctuation to induce GW background testable in the planned GW detectors [10-14]. In addition, a PBH domination imprints the spectrum of stochastic gravitational waves background (SGWB) originating from other independent sources. Such imprints on the GWs could be probed with the spectral shape reconstruction techniques [15] planned for the GW detectors. Therefore, a combined study of the GW spectrum associated with the ultralight PBHs-the GW spectrum obtained directly from PBHs plus inquiring into the imprints of a PBH-dominated phase on an independent GW spectrum could be a riveting way to understand the properties of ultralight PBHs. The latter case is the central theme of this article. We study imprints of PBHs on GW background produced by cosmic superstrings [16-24] and metastable strings [25-33] 1 . We classify the other GW backgrounds obtained directly from PBH evaporation or their dynamics in the early universe as complementary probes that could be envisaged to distinguish the effect of PBH domination on GWs from cosmic strings from any other early matter domination (EMD) scenarios, see, e.g., [24, 36].  \nThe motivation to consider the mentioned classes of strings is twofold. First, because the ultralight PBHs evaporate at early times, before the Big Bang Nucleosynthesis (BBN) [37-39], they generically affect any propagating GW spectrum during the PBH-domination phase at higher frequencies. Cosmic strings radiate GWs spanning a wide range of frequencies; generally with a strong scale-invariant amplitude at higher frequencies. Therefore, ultralight PBHs leave possibly  \nthe cleanest detectable imprints on such GW spectrum. In addition, a wide-spanning GW signal has the following advantage: if a measurement at a certain frequency band fixes the parameter space describing the GW source, signals at other frequency bands get automatically constrained. This redirects us to the second motivation. Cosmic superstrings and metastable strings are speculated to be among the possible sources of nHz SGWB recently discovered by the Pulsar Timing Arrays (PTA) [40-45]. Unlike stable cosmic strings [46-49] described only by the parameter: Gµ , where G is the Newton's constant and µ ∼ v 2 Φ is the string tension with v Φ being the vacuum expectation value of the scalar field describing the strings [50-52], cosmic superstrings and metastable cosmic strings fit the data better [24, 44, 53-60] because the last two being a two-parameter string model. In addition to Gµ , the inter-commutation probability P for superstrings and the parameter κ describing the separation between the Grand Unified Theory (GUT) scale and the string scale for metastable strings provide additional degrees of freedom to match the required GW amplitude and spectral slope at nHz frequencies. On the other hand, the high-frequency scale-invariant behaviour of the GW spectrum remains similar for all string classes. Fitting PTA data with GWs from superstrings and metastable strings predicts extremely strong scale-invariant amplitude at higher frequencies if the universe undergoes standard cosmological expansion. In fact, the high-frequency amplitudes are so strong that a significant string parameter space consistent with the PTA data gets ruled out by the LIGO-O3 run which did not observe any SGWB [61]. This puts superstring and metastable string interpretation to the PTA data perhaps in jeopardy unless an additional GW contribution from supermassive black holes (SMBH) is considered to relax the tension somewhat [44]. A PBH-dominated phase in the early universe, however, can cause a characteristic fall-off of the spectrum at higher frequencies, alleviating the tension with the PTA and LIGO completely. Besides, because the PTA data constrains the string parameters and the PBH mass is bounded from above so that they evaporate before BBN, even though reduced, the predicted overall signal strength could be pretty strong.  \nFocusing exclusively on the parameter range of superstrings and metastable strings favoured by the PTA data, and considering non-rotating PBHs with monochromatic initial mass, we find that PBHs with an initial mass within the range M BH ∈ /bracketleft.alt1 10 6 g , 10 8 g /bracketright.alt would imprint the GWs at higher frequencies so that they evade the LIGO bound. We broadly classify such viable signals into two classes: whether or not the next LIGO run (LVK-Design) [62] could probe the PBHimprinted GW signal obtained from superstrings and metastable strings. We find that for most of the mentioned mass range, PBHs can produce signals at LVK-D while being consistent with the PTA data. Heavier PBHs with an initial mass around 10 8 g reduce the GW amplitude such that the signals fall short of the sensitivity reach of LVK-D, however strong enough to produce a Signal-to-Noise-Ratio at ET [63] as large as ∼ 10 5 .  \nIn the companion theory paper by NANOGrav collaboration [44], the effect of an EMD on GWs from superstrings and metastable strings has been mentioned to evade the LIGO-O3 bound. Therefore, root-idea-wise, this article does not offer anything new. We rather report if the PBHs provide such EMD, the scenario becomes pretty distinct compared to any other EMD, e.g., provided by a scalar field. This is because, first, the PBHs provide additional GW signatures. For the obtained mass range M BH ∈ /bracketleft.alt1 10 6 g , 10 8 g /bracketright.alt , although the GW in the form of gravitons appears at extremely high frequency, GWs, e.g., induced by PBH density fluctuation appear in the detectable range. It is not trivial to reconstruct such a combined GW signal by any other means. Second, lately, there has been an extensive effort to understand the effects of ultralight PBHs on Beyond the Standard Model (BSM) physics such as particle dark matter (DM) production and baryon asymmetry of the universe (BAU), see, e.g., Refs.[64-67] (subsequent references are in Sec.IV B). For the M BH ∈ /bracketleft.alt1 10 6 g , 10 8 g /bracketright.alt , such effects could be quite significant, specifically, BSM physics occurring around the Electroweak (EW) scale, even down to the BBN (MeV) scale. In this context, we discuss a recently proposed BSM model [68, 69] accommodating neutrino mass, BAU via leptogen-  \n, super-heavy DM, and metastable strings (the last feature is new). We show how a PBH-string landscape could be capable of probing untrod parameter regions of BSM models.", 'II. STABLE STRINGS, SUPERSTRINGS, METASTABLE STRINGS, AND NON-STANDARD COSMOLOGY': "In this section, we briefly discuss/review the foremost aspects of cosmic superstrings and metastable strings required to compute GW spectra in the presence of PBHs. The gravitational wave spectral features obtained from the superstring and metastable string network differ from the standard stable strings mostly at the low frequencies, e.g., in the PTA band. On the other hand, imprints of early matter domination (by PBHs in our case) on the GW spectra at higher frequencies are similar to all string classes. It is therefore useful to highlight first the properties and the technicalities required to compute the GW spectrum from stable cosmic strings. The methodology used can be generalized to the case of superstrings and metastable strings. In addition, it would be opportune to understand why the GWs from stable cosmic strings do not fit the recent PTA data, unlike the other two string classes.  \nStable cosmic strings : Gravitational waves are radiated from cosmic string loops chopped off from the long strings resulting from the spontaneous breaking of the gauged U ( 1 ) [70, 71]. Long strings are described by a correlation length L = /radical.alt1 µ /slash.left ρ ∞ , with ρ ∞ is the long string energy density and µ is the string tension defined as [72]  \nµ = πv 2 Φ h /parenleft.alt1 λ, g ' /parenright.alt1 , (II.1)  \nh ( λ, g ' ) /uni2243 1, unless the coupling λ and g ' are not strongly hierarchical [73, 74]. The time evolution of a radiating loop of initial size l i = αt i is given by l ( t ) = l i -Γ Gµ ( t -t i ) , where Γ /uni2243 50 [70, 71], α /uni2243 0 . 1 [75, 76], G is the Newton's constant, and t i being the initial time of loop production. The total energy loss from a loop can be decomposed into a set of normal-mode oscillations with instantaneous frequencies f k = 2 k /slash.left l k = a ( t 0 )/slash.left a ( t ) f , where k = 1 , 2 , 3 ...k max , f is the present-day frequency at t 0 , and a is the scale factor. The total GW energy density can be computed by summing all the k modes giving [75, 76]  \nΩ GW = k max /summation.disp k = 1 2 k F α Gµ 2 Γ k fρ c /integral.disp t 0 t i /bracketleft.alt4 a ( t ) a ( t 0 ) /bracketright.alt4 5 n ω ( t, l k ) d t , (II.2)  \nwhere ρ c is the critical energy density of the universe, F α /uni2243 0 . 1 is an efficiency factor [75] and n ω ( t, l k ) is the loop number density which can be computed from the velocity-dependent-one-scale model as [77-80]  \nn ω ( t, l k ) = A β α ( α + Γ Gµ ) 3 ( 1 -β ) [ l k ( t ) + Γ Gµt ] 4 -3 β t 3 β . (II.3)  \nIn Eq.(II.3), β = 2 /slash.left 3 ( 1 + ω ) with ω is the equation of state parameter of the universe, and A β = 5 . 4 ( A β = 0 . 39) for radiation-dominated (matter-dominated) universe [80]. The quantity Γ k = Γ k -δ /slash.left ζ ( δ ) quantifies the emitted power in k -th mode, with δ = 4 /slash.left 3 ( δ = 5 /slash.left 3) for loops containing small-scale structures such as cusps (kinks) [81] and makes it evident that k = 1 (fundamental) mode is the dominant one contributing to the GWs. In principle, the GW spectrum corresponding to the fundamental mode broadly captures most of the underlying features of a given BSM model producing cosmic strings. However, when it comes to comparing the spectrum with real data, e.g.,  \nPTA data, and based on that, one intended to forecast future predictions, it is desirable to sum a large number of modes to obtain accurate results.  \nThe integral in Eq. (II.2) is subjected to the Heaviside functions Θ ( t i -t fric ) Θ ( t i -l c /slash.left α ) that set cut-offs on the GW spectrum at high frequencies f ∗ above which the spectrum falls as f -1 for the fundamental mode. The quantity t fric represents the time until which the string network is damped by friction [82]. For l > l c , the GW emission dominates over particle production as shown by recent numerical simulations [83, 84]. The critical length can be approximately computed as l c /uni2243 δ w ( Γ Gµ ) -γ where δ w = ( √ λv Φ ) -1 is the width of the string and γ = 2 ( γ = 1) for loops containing cusps (kinks). In most cases, the cut-offs owing to friction and kinks are weaker (occur at higher frequencies). We therefore focus on the string loops containing cusps only.  \nWithout any intermediate matter-dominated epoch, Eq.(II.2) can be solved to obtain the GW spectra taking into account two main contributions: i) from the loops created in the radiation domination but emitting in the standard matter domination, ii) from the loops created as well as emitting in the radiation domination. The former corresponds to a low-frequency peak while the latter produces a scale-invariant plateau at high frequencies (see Fig.1) which for k = 1 mode is given by [75, 76, 85]  \nΩ plt GW = 128 π F α Gµ 9 ζ ( δ ) A R ϵ R Ω R /bracketleft.alt2( 1 + ϵ R ) 3 /slash.left 2 -1 /bracketright.alt2 . (II.4)  \nIn Eq. (II.4), ϵ R = α /slash.left Γ Gµ /uni226B 1, A R ≡ A β /uni2243 5 . 4 and Ω R ∼ 9 × 10 -5 . Note that the plateau amplitude Ω plt GW ∝ √ µ ∝ v Φ , implying strong amplitude GWs are produced for larger symmetry breaking scale. In scenarios that feature an early matter-dominated epoch (in this work by PBHs) instead, the plateau breaks at a high frequency f brk , beyond which the spectrum falls as Ω GW ( f > f brk ) ∼ f -1 for the fundamental mode (green and red curves on the top panel in Fig.1). This spectral break frequency f brk can be estimated as [86]  \nf brk /uni2243 0 . 45Hz /parenleft.alt4 10 -12 Gµ /parenright.alt4 1 /slash.left 2 /parenleft.alt3 T brk GeV /parenright.alt3 , (II.5)  \nwhere T brk ≡ T end is the temperature corresponding to the end of early matter domination. When summed over a large number of modes ( k →∞ ), two important changes occur in the GW spectrum. First, the plateau amplitude increases approximately by a factor of ζ ( δ ) , i.e., Ω plt , ∞ GW /uni2243 Ω plt GW ζ ( δ ) and the spectrum falls as f 1 -δ instead of f -1 [49, 87].  \nGiven an observed spectrum with a plateau followed by a fall-off, one has to compare two parallel scenarios to pinpoint the origin of such a spectrum-whether the EMD or the particle production cut-off causes the spectral shape. Considering the former is equivalent to imposing a constraint on the parameter space as T cusp /slash.left T brk > 1, where we derive T cusp as  \nT cusp /uni2243 2 . 3 × 10 4 GeV /parenleft.alt3 Gµ 10 -12 /parenright.alt3 5 /slash.left 4 . (II.6)  \nSimilar to Eq.(II.5), an expression for f cusp can therefore be obtained as  \nf cusp /uni2243 1 . 01 × 10 4 Hz /parenleft.alt3 Gµ 10 -12 /parenright.alt3 3 /slash.left 4 (II.7)  \nwhich we find to be consistent with Eq.(50) of Ref.[84]. While analyzing the GW spectrum with PBHs, one has to make sure that f cusp > f brk is satisfied (see Fig.1, top panel).  \nCosmic superstrings : Cosmic superstrings are characterized by a weaker inter-commutation probability P because of their quantum-mechanical nature. The lack of availability of numerical simulations with larger simulation times makes the actual GW spectrum produced by superstrings pretty uncertain. Nonetheless, the standard paradigm to compute the spectrum has been to multiply the loop number density by a factor of P -1 [21, 88] (see Ref.[23] for a different scaling relation). This method has also been adopted by the NANOGrav collaboration [44]. Therefore, we shall use all the methodologies considered for stable strings but rescaling the loop number density by P -1 , with P ∈ /bracketleft.alt1 10 -3 , 1 /bracketright.alt [44]. Note that, this is not a rescaling of Gµ that one obtains for hierarchical coupling ( h ( λ, g ' ) ≠ 1) [73, 74], thereby making the spectral shape, e.g., in the PTA band, different from the stable cosmic strings.  \nMetastable cosmic strings : Metastable cosmic strings decay into string segments connecting monopole-antimonopole pairs by quantum tunneling [25, 26]. In the semiclassical approximation, the network is characterized by the decay rate per string unit length as  \nΓ d = µ 2 π exp (-πκ ) , √ κ = m M µ 1 /slash.left 2 /uni2243 Λ GUT v Φ , (II.8)  \nwhere m M is the monopole mass which is of the order of Λ GUT /uni2243 10 16 GeV. From Eq.(II.8), one can define a time t s = 1 Γ 1 /slash.left 2 d , before which the cosmic string network behaves similarly to a stable cosmic string network, around which the monopole nucleation starts, and after which, the network behaves as a decaying network with new loops no longer produced. For t > t s the loop number density in Eq.(II.3) receives a correction factor [30]:  \nD ( l, t ) = e -Γ d /bracketleft.alt1 l ( t -t s )+ 1 2 Γ Gµ ( t -t s ) 2 /bracketright.alt Θ /parenleft.alt1 αt s -¯ l ( t s )/parenright.alt1 , (II.9)  \nwhere the Θ function with ¯ l ( t s ) /uni2243 l + Γ Gµt , represents the fact that only the loops produced before t s contribute to the number density. In principle, unless the κ is quite large (which is the case for stable strings), it is sufficient to take the upper limit of the integration in Eq.(II.2) as t ∼ t eq , because during t s < t < t eq the loop number density encounters rapid exponential suppression. In this work, we consider GW emission only from loops. Nonetheless, when monopoles do not carry unconfined flux, cosmic string segments formed owing to monopole nucleation can contribute to GWs [30]. Compared to the stable cosmic strings, two distinct features of the GW spectrum from metastable strings are i) for smaller values of κ , the frequency at which the spectrum reaches the plateau gets exponentially right-shifted on the Ω GW vs. f plane (the blue curve in Fig.1). As such, for sufficiently small values of κ , the spectrum reaches the plateau at a very high frequency f /uni2273 10 Hz, thus avoiding all the existing bound on Gµ , barring the one expected from BBN bound on the effective number of degrees of freedom [89]. ii) The infra-red tail of the GW spectrum behaves as Ω GW ∝ f 2 [30, 31] favored by the PTA data when fitted to a power-law GW spectrum. Markedly, only a few other models exhibit Ω GW ∝ f 2 power-law in the infra-red tail, see, e.g., [90, 91]. Note that, compared to the stable cosmic strings, in both cases, we have an extra parameter in addition to Gµ : the intercommutation probability P for superstrings, and κ for metastable strings. This additional degree of freedom (i.e., considering a two-parameter string model) makes it easier to fit the PTA data which are characterized also by two parameters; the amplitude and the spectral index.  \nWe conclude the discussion by showing the GW spectra in Fig.1 for superstrings and metastable strings comparing them to the ones from stable strings. We choose a benchmark value of T brk = 0 . 1 GeV and consider only the k = 1 mode, leaving a more realistic analysis by summing a large number of modes for the next section. In Fig.1 (left) we show the GW spectrum for cosmic super (stable) strings with the green solid (dashed) curve. For the superstrings, we choose the  \nFIG. 1. All GW spectra are produced for k = 1 mode. Top-left: GW spectra for cosmic super (stable) strings shown with green solid (dashed) curve. Two benchmark string parameters chosen for superstrings are P = 6 × 10 -3 and Gµ = 2 × 10 -12 . For the stable strings, we rescale Gµ as Gµ → Gµ /slash.left P 2 so that comparable plateau amplitude is achieved in both cases. The horizontal and vertical dotted black lines represent the analytical expression for the plateau amplitude (multiplied by P -1 for superstrings), and the spectral break frequency f brk corresponding to a long-lasting matter domination ending at T = 0 . 1 GeV. The green vertical lines correspond to the spectral break frequencies where the spectrum would have fallen due to particle production if there were no matter domination. Top-right: GW spectra for cosmic metastable (stable) strings are shown with the red solid (dashed) curve. In this case, two benchmark string parameters are chosen as √ κ = 8 and Gµ = 2 × 10 -7 . The dotted horizontal and vertical lines represent the same as in the left panel figure. The red vertical line represents the frequency at which the spectrum would have fallen if there was no matter domination. The GW spectrum represented by the blue curve is for √ κ = 6 assuming no matter domination. The density reach of SKA [92], LISA [93], DECIGO [94] and LVK-D [62] are shown with the filled coloured curves. Bottom: The benchmark GW spectra in the PTA band with the recent posteriors reported by NANOGrav recently.  \n<!-- image -->  \nbenchmark values P = 6 × 10 -3 and Gµ = 2 × 10 -12 while for the stable cosmic strings, we rescale Gµ as Gµ → Gµ /slash.left P 2 such that comparable plateau amplitude can be obtained in both cases. Note that, for superstrings, one has the required blue-tilted GW spectrum in the PTA band (A bluetilted distribution of fourteen NANOGrav posteriors is shown in black; in this work, we use only NANOGrav data which are consistent with all other PTAs (IPTA) [95] ). On the other hand, even though the required GW amplitude can be obtained for the stable strings, the spectrum is nearly scale-invariant in the PTA band and thereby does not fit the data. A similar situation occurs in the case of metastable strings as shown in the top right panel with the red curves 2 . In both plots, the black dotted lines represent the analytical expression for the plateau amplitude and the spectral break frequency given in Eq.(II.4) and (II.5), respectively. The green and the red vertical lines represent the spectral break frequency f cusp at which the spectrum would have fallen if there was no matter domination. In the bottom panel, we zoomed in on the GW spectra within the PTA band to show the consistency of the blue-tiled GW spectrum obtained from cosmic superstrings and metastable strings with the data.  \nAn important aspect of this study is considering a falling GW spectrum at higher frequencies to avoid the LIGO bound. Although we zero in on early matter domination for this purpose, as discussed earlier, the spectrum may also fall because of particle production. With a straightforward generalization of Eq.(11) of Ref.[83] with δ w = ( √ λv Φ ) -1 , it can be shown that for making the spectrum fall before the LIGO frequency band while being consistent with the PTA data, an extremely small value of λ is required. For the standard case with λ ∼ 1 considered in the numerical simulations, the spectrum always falls at very high frequency owing to particle production. This makes the case of an early matter-dominated scenario pragmatic while interpreting the PTA data as cosmic metastable (super) strings radiated GWs evading the LIGO-O3 bound.", 'III. IMPRINTS OF ULTRALIGHT PBHS ON GWS FROM COSMIC (SUPER/METASTABLE) STRINGS': 'We focus on non-rotating PBHs with monochromatic initial mass M BH . The energy density of the PBHs ( ρ BH ) and radiation ( ρ R ) evolve according to the following Friedmann equations [49, 96, 97]:  \ndρ R dz + 4 z ρ R = 0 , (III.1)  \ndρ BH dz + 3 z H ˜ H ρ BH -˙ M BH M BH 1 z ˜ H ρ BH = 0 , (III.2)  \nwhere H is the Hubble parameter ( H 0 ∼ 1 . 44 × 10 -42 GeV), z = T Bf /slash.left T , with T Bf being the PBH formation temperature and the time derivative ˙ M BH is given by  \n˙ M BH = G g ∗ B ( T BH ) 30720 π M 4 Pl M 2 BH , (III.3)  \nwhere M Pl = 1 . 22 × 10 19 GeV, G /uni2243 3 . 8 is a grey body factor and g ∗ B ( T BH ) /uni2243 100 being the number of relativistic degrees of freedom below the PBH temperature T BH with Standard Model (SM) containing left-handed neutrinos [98]. The quantity ˜ H and the scale factor a evolve as  \n˜ H = ( H +K) , da dz = /parenleft.alt3 1 -K ˜ H /parenright.alt3 a z , (III.4)  \nwhere K = ˙ M BH M BH ρ BH 4 ρ R . To derive Eq.(III.1)-Eq.(III.4), we assume the entropy ( g ∗ s ) and the energy ( g ∗ ρ ) degrees of freedom are equal and constant. For a given value of β ≡ ρ BH ( T Bf ) ρ R ( T Bf ) and M BH , the above equations can be solved to determine the duration of PBH domination and the evaporation temperature. In the case where the PBH dominates the energy density, an analytical expression for the evaporation temperature is obtained as  \nT ev = /parenleft.alt4 5 M 2 Pl π 3 g ∗ ( T ev ) τ 2 /parenright.alt4 1 /slash.left 4 , (III.5)  \nwhere τ is the PBH lifetime which can be computed by integrating Eq.(III.3). This allows also to express T ev in terms of M BH as  \nT ev = 0 . 1 /parenleft.alt3 M BH 4 . 1 × 10 7 g /parenright.alt3 -3 /slash.left 2 GeV . (III.6)  \nEq.(III.6) shows that with T ev ≡ T brk , all the plots in Fig.1 can be reproduced for PBHs with initial mass M BH /uni2243 4 . 1 × 10 7 g .  \nIt is crucial to note the following: i) An eternal falling of the GW spectrum in the range f brk /uni2272 f /uni2272 f LIGO , (see Fig.1), requires longer duration of PBH domination and therefore large β -we consider this to constrain PBH mass with PTA data and the LIGO bound plus projection. However, in scenarios where β can be constrained (can be bounded from above, resulting in a short duration of PBH domination) by GWs from PBH density fluctuations [11], the falling of the GW spectrum might stop before the LIGO frequency band and another subsequent plateau with suppressed amplitude (originating out of the contribution from the loops in the first radiation epoch before the PBH start to dominate) might show up. Therefore one needs to modify our results accordingly, see, e.g., sec.IV. ii) The PBHs should evaporate before the BBN at a temperature T /uni2273 10 MeV. Therefore, from Eq.(III.6), M BH is generically constrained as M BH /uni2272 2 × 10 8 g , which we use as the absolute upper bound on M BH .  \nWe now proceed to the detailed discussion of the impact of PBH-induced matter domination on GWs from cosmic superstrings and metastable strings. First, from Eq.(II.5) and Eq.(III.6), we obtain the spectral break frequency as  \nf PBH brk /uni2243 0 . 045Hz /parenleft.alt4 10 -12 Gµ /parenright.alt4 1 /slash.left 2 /parenleft.alt3 M BH 4 . 1 × 10 7 g /parenright.alt3 -3 /slash.left 2 . (III.7)  \nWhen a large number of modes are summed, GW energy density beyond f PBH brk can be written as  \nΩ GW ( f ) /uni2243 Ω plt , ∞ GW /parenleft.alt4 f f PBH brk /parenright.alt4 -1 /slash.left 3 Θ /bracketleft.alt1 f -f PBH brk /bracketright.alt , (III.8)  \nwhere  \nΩ plt , ∞ GW /uni2243 3 . 1 × 10 -8 Υ -1 /parenleft.alt3 Gµ 10 -7 /parenright.alt3 1 /slash.left 2 , /uni23A7 /uni23AA /uni23AA /uni23A8 /uni23AA /uni23AA /uni23A9 Υ = P for superstrings Υ = 1 for metastable strings . (III.9)  \nUsing Eq.(III.7), Eq.(III.8) can be re-expressed as  \nΩ GW ( f ) /uni2243 Ω plt , ∞ GW Ω mod , PBH GW , (III.10)  \nwhere Ω mod , PBH GW , given by  \nΩ mod , PBH GW /uni2243 0 . 1 /parenleft.alt4 10 -12 Gµ /parenright.alt4 1 /slash.left 6 /parenleft.alt3 M BH 4 . 1 × 10 7 g /parenright.alt3 -1 /slash.left 2 /parenleft.alt3 f 25Hz /parenright.alt3 -1 /slash.left 3 (III.11)  \nFIG. 2. When PBHs dominate in the early universe: Left (right): Gravitational wave spectra for three benchmark points [see Table I, II] from each category discussed in the text are shown for superstrings (metastable strings). We consider GWs radiated from loops only for metastable strings, discarding contribution from string segments. Each benchmark GW spectrum is consistent with PTA data and evades LIGO-O3 bound on stochastic gravitational waves.  \n<!-- image -->  \ncaptures the modulation in the GW spectra because of PBH domination. The absence of PBHs corresponds to Ω mod , PBH GW = 1. In which case, Eq.(III.9) and the non-observation of SGWB by LIGO-O3 run; i.e., Ω anymodel GW /uni2272 1 . 7 × 10 -8 at f LIGO = 25 Hz, implies the following:  \nwithoutthePBH , LIGO -O3excludes /uni23A7 /uni23AA /uni23AA /uni23AA /uni23A8 /uni23AA /uni23AA /uni23AA /uni23A9 P /uni2272 1 . 8 /parenleft.alt2 Gµ 10 -7 /parenright.alt2 1 /slash.left 2 for superstrings , Gµ /uni2273 3 × 10 -8 for metastable strings . (III.12)  \nThough Eq.(III.12) is the most interesting region to explore with non-standard cosmological evolution, e.g., with PBHs, for a comprehensive discussion, we present our numerical results by dividing parameter space into three different categories (all cases are consistent with PTA data):  \nCategory A: Ω plt , ∞ GW /uni2273 Ω LIGO -O3 GW and Ω LIGO -O3 GW /uni2273 Ω plt , ∞ GW Ω mod , PBH GW ( f = 25Hz ) /uni2273 Ω LVK -D GW .  \nCategory B: Ω LIGO -O3 GW /uni2273 Ω plt , ∞ GW /uni2273 Ω LVK -D GW and Ω LIGO -O3 GW /uni2273 Ω plt , ∞ GW Ω mod , PBH GW ( f = 25Hz ) /uni2273 Ω LVK -D GW .  \nCategory C: Ω LVK -D GW /uni2273 Ω plt , ∞ GW Ω mod , PBH GW ( f = 25Hz ) for all Ω plt , ∞ GW .  \nHere Ω LVK -D GW = 2 × 10 -9 at f = 25 Hz is the maximum projected sensitivity (minimum testable value of Ω GW ) reach of the LVK-D. Note that Category A corresponds to recovering the excluded parameter space in Eq.(III.12), for which one must require a non-standard cosmological evolution, e.g., by PBH. For Category B, PBHs are not essential because even the scale-invariant amplitude does not contradict the LIGO-O3 bound. A PBH domination however can be tested for f PBH brk /uni2272 25Hz and Ω LIGO -O3 GW /uni2273 Ω plt , ∞ GW Ω mod , PBH GW ( f = 25Hz ) /uni2273 Ω LVK -D GW . For Category C, though Ω GW is consistent with the PTA data, the LVK-D cannot test the spectrum because the amplitude is smaller than its sensitivity reach. However, detectors like the Einstein Telescope (ET) can potentially test such signals. We quantify this class of signals by computing Signal-to-Noise-Ratio  \nTABLE I. Benchmarks: SuperstringsTABLE II. Benchmarks: Metastable strings  \nat ET. In Fig.2, we show the benchmark GW spectrum (summing over larger modes) from each category belonging to both string classes. The benchmark parameters are tabulated in Table I (Superstrings) and Table II (Metastable strings). In what follows, we present detailed numerical results including a more rigorous parameter space scan to find an approximate mass range of PBHs that can produce such signals.', 'A. Numerical results: PBH mass ranges for PTA GW signals evading LIGO-O3 bound': 'On the top panel in Fig.3, we show the allowed regions for fixed PBH mass that would produce GW signals in each category belonging to cosmic superstrings. The gray contours in each plot are the NANOGrav 1 σ and 2 σ posteriors of Gµ and P [44] (one may also look at the Appendix A of this article). The benchmark points tabulated in Table I corresponding to the GW spectra in Fig.2 (left) are denoted with the symbols ♣ BP1 (Cat-A), ♠ BP2 (Cat-B), and /uni22C6 BP3 (Cat-C). A similar description holds for metastable strings on the Gµ and √ κ plane (Fig.4, top) with the benchmarks tabulated in Table II and corresponding GW spectra shown in Fig.2 (right). In the bottom panel of both figures, the GW spectra presented in Fig.2 are zoomed in to the NANOGrav and the LIGO frequency band. Note that for all categories in both figures, only a constrained range of PBH mass M BH ∼ /bracketleft.alt1 10 6 g, 10 8 g /bracketright.alt would produce a viable GW signal consistent with the PTA data. For category C, slightly heavier PBHs are allowed; to make the spectrum fall below Ω LVK -D GW with f brk at lower frequencies. On the other hand, for Category B, only lighter PBHs are allowed for the opposite reason, i.e., Ω LVK -D GW /uni2272 Ω GW ( f = 25Hz ) requiring f brk /uni2272 25Hz at much higher frequencies. Our parameter space scanning does not account for a correction factor due to the possible variation of effective degrees of freedom, which would make plateau amplitude slightly less. Therefore, in this analysis, the LIGO constraints on the string parameter space are slightly stronger than they should be [44]. As such, from Fig.4, it seems that without PBH the NANOGrav 2 σ region is excluded by LIGO-O3. But when the correction factor is considered, a small portion of 2 σ region around Gµ ∼ 10 -7 does not contradict the LIGO-O3 bound [44]. Therefore, one can produce Category-B signals for metastable strings with M BH ∼ 10 5 . 5 g .', 'B. Numerical results: Expected Signal-to-Noise-Ratio at the Einstein Telescope': 'The standard way to assess the detectability of predicted GW signal in the future planned detectors is to compute the Signal-to-Noise-Ratio (SNR) [99, 100]. Because the PTA GWs are of strong amplitudes, a generic expectation from any such signal spanning a wide frequency range would be to obtain a large SNR at the high-frequency detectors. In the present case, because PBH domination reduces the signal amplitude at higher frequencies, it is useful to quantify the detectability of such signals with SNR. Note that, in our case, the higher the frequency smaller the amplitude of GWs. One gets the most reduced signal at higher frequencies for Category C. We therefore assess those signals by computing SNR at the Einstein Telescope. The SNR is defined  \nFIG. 3. All plots are for superstrings. Top panel, left: The coloured regions marked with different PBH mass M BH on the Gµ -P plane, correspond to the GW spectra with plateau amplitude Ω plt , ∞ GW /uni2273 Ω LIGO -O3 GW ≡ 1 . 7 × 10 -8 and testable with LVK-D at 25 Hz. In the absence of PBHs, the regions would be excluded by LIGO-O3 because in that case, the GW spectra do not fall at f /uni226A 25 Hz, and the scale-invariant amplitudes contradict the LIGO bound. The regions inside the gray contours are consistent with the NANOGrav 15yrs data. Middle: The colored regions correspond to the GW spectra Ω LIGO -O3 GW /uni2273 Ω plt , ∞ GW /uni2273 Ω LVK -D GW ≡ 2 × 10 -9 . Even in the absence of PBHs, these regions do not contradict the LIGO-O3 bound. But in that case, LVK-D would probe a scale-invariant spectrum. Right: These regions correspond to GWs with smaller amplitude at the LIGO frequency band so that LVK-D can not test them, but detectors such as ET can. For the chosen benchmark point from each category (Table. I) the GW spectrum is shown in Fig.2 (left). Bottom panel: Zoomed in GW spectra in the PTA band against the NANOGrav posteriors (the violins) and in the LIGO frequency band.  \n<!-- image -->', 'as [99, 100]': 'SNR = /roottop /rootmod /rootmod /rootbot t obs /integral.disp f max f min df /parenleft.alt4 Ω GW ( f ) Ω noise ( f ) /parenright.alt4 2 , (III.13)  \nwhere Ω noise ( f ) is the noise spectrum of the detector and f min ( f max ) are minimum (maximum) accessible frequency. In Fig.5, we present the SNR@ET for PTA GWs radiated from superstrings (metastable strings) in the upper (lower) panel for the benchmark values M BH = 10 7 , 8 g. Note that even though at the ET band the signal strength reduces, the predicted SNR is pretty high. This is  \nFIG. 4. All plots are for metastable strings. Top panel, left: The coloured regions marked with different PBH mass M BH on the Gµ -κ plane, correspond to the GW spectra with plateau amplitude Ω plt , ∞ GW /uni2273 Ω LIGO -O3 GW ≡ 1 . 7 × 10 -8 and testable with LVK-D at 25 Hz. In the absence of PBHs, the regions would be excluded by LIGO-O3 because in that case, the GW spectra do not fall at f /uni226A 25 Hz, and the scale-invariant amplitudes contradict the LIGO bound. The regions inside the gray contours are consistent with the NANOGrav 15yrs data. Middle: The colored regions correspond to the GW spectra Ω LIGO -O3 GW /uni2273 Ω plt , ∞ GW /uni2273 Ω LVK -D GW ≡ 2 × 10 -9 . Even in the absence of PBHs, these regions do not contradict the LIGO-O3 bound. Right: These regions correspond to GWs with smaller amplitude at the LIGO frequency band so that LVK-D can not test them, but detectors such as ET can. For the chosen benchmark point from each category (Table II) the GW spectrum is shown in Fig.2 (right). Bottom panel: Zoomed in GW spectra in the PTA band against the NANOGrav posteriors (the violins) and in the LIGO frequency band.  \n<!-- image -->  \nowing to the required strong amplitude GWs at the PTA band plus the maximally allowed value of PBH mass which is around M BH ∼ 10 8 g -the BBN constraint.  \nGiven the above discussion and the parameter scan rendering an approximate (fairly accurate) range of PBH masses that make string-radiated GWs consistent with the LIGO bound, we now discuss the phenomenological advantage of having a PBH domination over any other early matter domination.  \nFIG. 5. Top [bottom]: Expected Signal-to-Noise-Ratio (SNR) at the Einstein Telescope for the PTA GWs radiated from cosmic superstrings [metastable strings] for M BH ∼ 10 7 g (left), and M BH ∼ 10 8 g (right). The SNR in each plot has been computed considering a f -1 /slash.left 3 fall of the GWs beyond f brk .  \n<!-- image -->', 'IV. PBHS OVER ANY OTHER EMD: ADDITIONAL PHYSICS': "As mentioned in the introduction, the advantage is broadly twofold. First, PBHs themselves are additional sources of GWs, thereby making the overall stochastic GW spectrum peculiar. As such, PBHs inevitably produce high-frequency gravitons, plus, as shown recently, the inhomogeneity in the distribution of PBHs over the spaces could induce tensor modes if the PBHs dominate the universe's energy density. As we show below, remarkably, for the obtained PBH mass range in Fig.3 and Fig.4, the peak of such induced GWs from PBH density fluctuations appears in the LISA-DECIGO-ET band. Second, lately, there have been considerable efforts to discuss the effect of ultralight PBHs on particle dark matter and baryogenesis models [64-67, 101-125]. The PBHs with the predicted mass range in this article would have profound effects (perhaps the best testable effects) on BSM particle physics models. To this end, we discuss a recently proposed BSM scenario [69] based on the seesaw mechanism [126] of neutrino masses stemming also the possibility of baryogenesis via leptogenesis [127] and superheavy dark matter in the presence of PBHs[68, 69]. In the scenario, it was reported for the first time, that if there are ultralight  \nPBHs in the early universe, seesaw Lagrangian can accommodate the novel right-handed neutrino (RHN) mass spectrum Λ GUT ∼ M R 3 > M R 2 > M R 1 , where M Ri being the mass of the i th RHN, with heaviest one ( N R 3 ) being super heavy dark matter and the other two generating baryon asymmetry. The analysis in Ref.[69] was made with stable cosmic strings produced owing to the spontaneous breaking of an U ( 1 ) B -L (embeddable in GUT) responsible also for generating the RHN masses. Here we present an analysis of the same model considering a metastable string network incorporating the new PTA data. In what follows, we present a concise discussion with numerical results for the additional GWs from PBHs and the U ( 1 ) PBH B -L seesaw model.", 'A. Complementary GW signatures: density fluctuations and graviton emission': "The ultralight PBHs are many-faced objects concerning the stochastic GWs. For instance, the initial curvature perturbations leading to the formation of PBHs give rise to GWs (see, e.g., [128130]), ultralight PBHs could possibly merge, releasing GWs [131], PBHs emit gravitons, which constitute GWs at ultra-high frequencies [6-9], and the inhomogeneous distribution of PBHs, which results in density fluctuations, triggers the production of GWs [10-12]. Here we do not discuss any specific production mechanisms of PBHs, nor the possibility of merger formation. We rather focus on the last two options, being slightly biased to the fact that graviton emission is inevitable and independent of the production mechanism, the GWs from PBH-density fluctuation seem to provide a realistic stochastic background with spectral shape dependent on their spatial distribution which has been assumed to be like poison like in Refs.[10, 11]. Even though the PBH gas, on average, behaves like pressure-less dust, the inhomogeneous spatial distribution results in density fluctuations, which are isocurvature in nature. When PBHs dominate the universe's energy density, the isocurvature component gets converted to curvature perturbation, which subsequently induces secondary GWs. Due to the significant density fluctuations at smaller scales (equivalent to the average separation of PBHs at T Bf ), substantial GWs are generated. These GWs are further amplified due to the nearly instantaneous evaporation of PBHs. The amplitude of such induced GWs in the present day is given by  \nΩ PBH , I GW ( t 0 , f ) /uni2243 Ω peak , I GW /parenleft.alt4 f f peak /parenright.alt4 11 /slash.left 3 Θ ( f peak , I -f ) , (IV.1)  \nwhere the peak amplitude is given by  \nΩ peak , I GW /uni2243 2 × 10 -6 /parenleft.alt3 β 10 -8 /parenright.alt3 16 /slash.left 3 /parenleft.alt3 M BH 10 7 g /parenright.alt3 34 /slash.left 9 , (IV.2)  \nand the peak frequency is given by  \nf peak , I /uni2243 5 . 4 Hz /parenleft.alt3 M BH 10 7 g /parenright.alt3 -5 /slash.left 6 . (IV.3)  \nThe ultra-violet cutoff frequency f peak corresponds to the comoving scale representing the mean separation of PBHs at the formation time.  \nOn the other hand, the present spectral energy density of GWs from gravitons is given by [6-8]  \nΩ GW ( f ) h 2 = ( 2 π ) 4 βγ 2 f 4 4 H 2 0 M BH /integral.disp t ev t Bf ( 1 -t -t Bf τ ev ) 2 /slash.left 3 /parenleft.alt2 a ( t Bf ) a ( t 0 ) /parenright.alt2 -3 exp /parenleft.alt3 16 π 2 fa ( t 0 ) M BH a ( t ) M 2 Pl /parenright.alt3 -1 (IV.4)  \nFIG. 6. Top-left: Combined GW spectrum (Metastable strings+density fluctuation) for the benchmark parameters belonging to Category C, see Table II. The spectrum has been produced for β = 5 × 10 -9 and summing over 10 3 modes of cosmic string loop oscillations. Top-right: GWs only from metastable cosmic strings for the same benchmark parameters but for different values of β , showing the impact of the finite duration of PBH domination. Bottom: GW strain h c from density fluctuation (dashed), from Metastable strings (solid), and graviton radiation from PBH (dotted) with the projected sensitivities (including the black curly curve) of various high-frequency GW detectors. All the benchmark parameters are the same as in the top-left plot. The region above the red dashed line is excluded from ∆ N eff bound by Planck.  \n<!-- image -->  \nFIG. 7. Top-left: Combined GW spectrum (Metastable strings+density fluctuation) for the benchmark parameters belonging to Category C, see Table I. The spectrum has been produced for β = 3 × 10 -8 and summing over 10 3 modes of cosmic string loop oscillations. Top-right: GWs only from metastable cosmic strings for the same benchmark parameters but for different values of β , showing the impact of the finite duration of PBH domination. Bottom: GW strain h c from density fluctuation (dashed), from Metastable strings (solid), and graviton radiation from PBH (dotted) with the projected sensitivities (including the black curly curve) of various high-frequency GW detectors. All the benchmark parameters are the same as in the top-left plot. The region above the red dashed line is excluded from ∆ N eff bound by Planck.  \n<!-- image -->  \nBy maximizing the Eq.(IV.4) w.r.t the present-day frequency f , one obtains the peak frequency and consequently the peak GW amplitude as  \nf peak , g /uni2243 2 . 8 a ( t ev ) T BH /uni2243 2 × 10 17 Hz /parenleft.alt3 M BH 10 7 g /parenright.alt3 1 /slash.left 2 , Ω peak , g GW ∼ 10 -7 . (IV.5)  \nIn Fig.6 and Fig.7 (top-left) we show the total GW spectrum: density fluctuation spectrum combined with the one from metastable strings and superstrings for the benchmark parameters belonging to Category C in each string class. Instead of large β , here we chose a benchmark value such that the duration of PBH domination is relatively short and the density fluctuation GW spectrum is visible (the sharp peak at higher frequencies). Note that because of the strong dependence of the peak GW amplitude on β (see Eq.(IV.2)) in this scenario, the duration of PBH domination can not be arbitrarily large, otherwise, the GWs would saturate the BBN bound on ∆ N eff : ∫ d f f -1 Ω GW ( f ) h 2 /uni2272 5 . 6 × 10 -6 ∆ N eff , with ∆ N eff /uni2272 0 . 2. Consequently, the falling of the GW spectrum from cosmic strings would stop at higher frequencies, potentially showing up a subsequent GW plateau originating from the loops in the radiation epoch before the PBH domination. This has been shown in the top-right panel of Fig.6 for metastable strings for three benchmark values of β . Nonetheless, the appearance of the second plateau is not always the case if the particle production cut-off is strong enough. We show it for superstrings on the top-right panel of Fig.7. Compared to the metastable string case, here owing to the strong particle production cut-off, the GW spectra closely follow the large β case. As discussed in Sec.II, for the GWs with amplitude at the level of PTAs, the particle production cut-off appears mostly at frequencies higher than the LIGO frequency band. Therefore, for shorter-duration PBH-domination scenarios, the regions in Fig.3 and Fig.4 extracted for Ω GW ∝ f -1 /slash.left 3 falling, requires modification now to account for a flattening of the spectrum. Although we do not present such a scan, to observe GW peak from density fluctuation on top of the strong amplitude cosmic string radiated GWs required by PTAs, reasonably large β (respecting, of course, the BBN bound) is required. For this reason, for most of the parameter space presented in Fig.3 and Fig.4, the duration of PBH domination would be such that from f brk until the f LIGO , the spectrum would fall as Ω GW ∝ f -1 /slash.left 3 , keeping the results of the previous section intact. We expect modification to the parameter space for heavier PBH masses around 10 8 g leading to smaller values of f brk , and thereby potential flatting of the spectrum at the LIGO frequency band due to the shorter duration of PBH domination.  \nBesides these low-frequency peculiarities, a PBH+superstrings (metastable strings) system provides a characteristic high-frequency GW spectrum constituting gravitons as shown in the bottom panels of Fig.6 and 7, where the GW strain h c is related to Ω GW as  \nh c ( f ) = /radical.alt4 3 H 2 0 2 π 2 Ω GW f -1 . (IV.6)  \nBecause Eq.(IV.5) characterizes the GW background, for the allowed mass ranges of PBH presented in Sec.III, a considerable GW background constituting gravitons can be obtained within a constrained frequency range spanning approximately, only a decade. Although the proposed high-frequency GW detectors [132-142] are not presently competitive with the ∆ N eff bound from Planck [143] (according to Ref.[144], the electromagnetic resonant cavities can potentially reach below the dark radiation projections around f ∼ 10 8 Hz; the black oscillating curve in Fig.6 and Fig.7), optimistically, in the future graviton background may be detected. In any case, these marked spectral features in the total GW background are described with three characteristic frequencies f brk , f peak , I , and f peak , g make the PBH + string scenario more appealing than an ordinary EMD+string case.  \nFIG. 8. A tentative timeline for the events occurring in the U ( 1 ) PBH B -L model. Without discussing any GUT scheme, this article assumes a metastable string network formation after the U ( 1 ) B -L breaking which also makes the right-handed neutrinos (one of them being a super-heavy dark matter) massive.  \n<!-- image -->", 'B. Implication on a realistic BSM model of leptogenesis and super-heavy dark matter': "In the seesaw mechanism [126], the light neutrino mass m i ∼ 0 . 01 eV inferred by neutrino oscillation data and cosmological measurements [143, 145], can be explained as a ratio of two scales: m i /uni2243 Λ 2 EW Λ GUT , with Λ EW and Λ GUT being the Electroweak and GUT scale, respectively. At the Lagrangian level, these two scales are introduced as L seesaw ∼ Λ EW LN R + Λ GUT N c R N R , where L is the Standard Model (SM) lepton doublet and N R is a RHN field. Besides generating light neutrino mass, the seesaw Lagrangian can explain the entire DM relic and the observed BAU, if there are three RHN fields, with one being the DM and the other two generating BAU through leptogenesis. A plethora of studies is dedicated to exploring this three-dimensional aspect of the seesaw-the neutrino mass+BAU+DM, see, e.g., [146-152]. Such attempts are theoretically wellbuilt and predictive in particle physics experiments, however, at the expense of fine-tuning the LN R Yukawa coupling and treating the RHN mass as a free parameter, thereby giving up the supposition that active neutrino masses arise from the GUT and EW scale physics. As such, perhaps the most studied scenarios based on [146, 147] consider N R to be extremely light O( GeV , keV ) , therefore tiny LN R coupling. In Ref.[69], it was proposed that if there are PBHs in the early universe, all three aspects of the seesaw can be addressed without significant deviation from the relation m i /uni2243 Λ 2 EW Λ GUT , i.e., considering all the RHNs to be O( Λ GUT ) , with one of them being the GUT-scale (super heavy) DM. In this case, the scenario becomes testable in GW experiments. The idea is based on dark matter production from ultralight PBH evaporation [64-67]. If such DM makes up 100% relic density without contradicting structure formation, the DM mass must be in the super-heavy range following the relation  \nM DM /uni2243 4 . 7 × 10 12 /parenleft.alt3 M BH 10 7 g /parenright.alt3 -5 /slash.left 2 GeV . (IV.7)  \nThe thermal leptogenesis mechanism (thermal scatterings populate the decaying RHNs producing the lepton asymmetry) generates the BAU. But in this case, the final asymmetry gets modified to account for an entropy dilution owing to the PBH domination. Once the condition for the observed baryon-to-photon ratio η B /uni2243 6 . 3 × 10 -10 [143] is taken into account, the parameter β and scale of  \nFIG. 9. Fit to the PTA data with U ( 1 ) PBH B -L model for neutrino mass, super-heavy DM and leptogenesis. The produced GW signal from metastable strings belongs to Category A discussed in the previous section, i.e., the GW spectrum consistent with the PTA data must produce a signal at the LVK-D.  \n<!-- image -->  \nthermal leptogenesis get correlated as [69]  \nβ /uni2243 1 . 2 × 10 -9 /parenleft.alt4 10 7 g M BH /parenright.alt4 /parenleft.alt3 T lepto 10 14 GeV /parenright.alt3 , (IV.8)  \nwhere we have assumed leptogenesis is dominated by the decay of N 1 so that M 1 ∼ T lepto . Eq.(IV.8) can be shown to be accurate, matching the numerical results obtained from Eq.(III.1),(III.2) and the modified leptogenesis equations in the presence of PBHs [49]. The PBHs can also populate the decaying RHNs to produce lepton asymmetry, but if Eq.(IV.7) is to be satisfied, such PBHs evaporate after Spaleron freeze-out temperature, and the lepton asymmetry does not get processed to BAU. Therefore we discard this option. Because the seesaw Lagrangian can be embedded in GUT groups with the RHN mass term generated after the breaking of an intermediate (While going to GUT → SM) symmetry U ( 1 ) B -L , cosmic strings are formed. Not exploring any explicit GUT scheme, we consider the formation of a metastable string network after the U ( 1 ) B -L breaking. The tentative timeline of this U ( 1 ) B -L model in the presence of PBHs looks like Fig.8. It is crucial to note that like the BAU, dilution of any pre-existing DM relic (e.g., from scattering by B -L gauge bosons) should be taken care of (the red line in Fig.8) in a way that Eq.(IV.7) prevails. A detailed discussion related to this can be found in Ref.[69]. Note now the features of the model. First, The DM and the PBH masses are anti-correlated (cf. Eq.(IV.7)). As such, ∼ 10 15 GeV DM gets produced from ∼ 10 6 g PBHs. Lighter PBHs are required to produce heavier DM. Second, as the leptogenesis scale increases (cf. Eq.(IV.8)), β increases, consequently the duration of PBH domination. Therefore, as the RHN mass scale (including the DM) approaches Λ GUT , remarkably, besides maintaining the m i /uni2243 Λ 2 EW Λ GUT relation to a great extent, the GWs produced from cosmic strings automatically evades LIGO bound owing to large β or the long duration of PBHdomination. Assuming the hierarchy in the RHN masses 3 as M DM > M Ri (= 1 , 2 ) we can re-write  \nEq.(II.1) as  \nµ = πM 2 DM f -2 DM h /parenleft.alt1 λ, g ' /parenright.alt1 , (IV.9)  \nwhere f DM is the DM Yukawa coupling which determines the interaction strength of N R field with B -L scalar; N c R N R Φ B -L . Working with O( 1 ) couplings in Eq.(IV.9), and using Eq.(IV.7) we obtain an expression for Gµ as  \nGµ /uni2243 4 . 7 × 10 -8 /parenleft.alt3 M BH 10 6 g /parenright.alt3 -5 . (IV.10)  \nWith this relation, we now redo the PTA fit for metastable strings performed in the previous section. The result is shown in Fig.9. Note that because of this extra condition in Eq.(IV.10), this U ( 1 ) PBH B -L model fits the PTA data for an extremely constrained range of PBH mass and produces GWsignals belonging only to Category A. Therefore any non-observation of SGWB in the LVK-D would rule out the model/the best parameter choice of the model. Besides, because the PBH mass window is narrow, the other GW signals, e.g., from density fluctuation plus in the form of gravitons are peaked in a stringent range in frequencies. In an upcoming publication, we shall provide a detailed discussion of the model for more general parameter space with all possible GW signatures.", 'V. CONCLUSION': 'An abundant presence of ultralight PBHs leading to their domination on energy densities in the early universe could produce plenty of observable signatures specifically in the form of gravitational radiation. In addition to producing direct GW background, e.g., in the form of gravitons, a PBH domination gets imprinted on the spectral shape of GWs originating from other independent sources. Such a spectral shape reconstruction with the GW detectors and studying complementary GWsignals originating directly from the PBHs could be an interesting way to probe the properties of ultralight PBHs. Here we study the impact of PBH domination on the GWs emitted from cosmic superstrings and metastable strings which are speculated to be among the sources of nHz GWs observed by the PTA. Such an interpretation might become perilous, as a GW background originating from cosmic strings spans a wide range in frequencies thereby contradicting the nonobservation of stochastic GW background by LIGO dealing with comparable signal strength. In fact, the LIGO constraint rules out a significant string-parameter space favoured by the recent PTA data. A PBH-dominated early universe, however, affects the propagation of such GWs to reduce the amplitude at the high frequencies thereby evading the LIGO constraints. Assuming monochromatic and non-rotating PBHs we work out the mass range of such ultralight PBHs which makes the string-radiated GWs consistent with LIGO bound and the PTA data. We show that PBHs with mass in the range M BH ∈ /bracketleft.alt1 10 6 g , 10 8 g /bracketright.alt would produce such viable signals that can be tested in the mid-band detectors such as LISA, high-frequency detectors such as LVK-D plus ET with a high signal-to-noise ratio. We discuss why a PBH-domination is distinguishable from any other early matter-domination scenarios that reduce the amplitude of string-radiated GWs at higher frequencies. In this context, we discuss a BSM model based on U ( 1 ) B -L -seesaw providing neutrino mass, leptogenesis, super-heavy dark matter, and metastable cosmic strings. We show how a PBH-cosmic string scenario could be extremely predictive in probing unexplored parameter space of such BSM scenarios.  \nFIG. 10. Triangular plot on the left (right) from the MCMC runs for superstrings (metastable strings) with the NANOGrav 15-year dataset using PTArcade [153].  \n<!-- image -->', 'ACKNOWLEDGEMENTS': 'The work of SD is supported by the National Natural Science Foundation of China (NNSFC) under grant No. 12150610460. The work of RS is supported by the research project TAsP (Theoretical Astroparticle Physics) funded by the Istituto Nazionale di Fisica Nucleare (INFN).', 'Appendix A: Superstrings and Metastable strings fit to the NANOGrav data': "In Fig.10 we present the superstring and metastable fit to the PTA data using PTArcade [153]. The 2D posteriors have been used in the main text.  \n- [1] S. Hawking, Mon. Not. Roy. Astron. Soc. 152 , 75 (1971).\n- [2] B. J. Carr and S. W. Hawking, Mon. Not. Roy. Astron. Soc. 168 , 399 (1974).\n- [3] A. Escriv'a, F. Kuhnel, and Y. Tada, (2022), 10.1016/B978-0-32-395636-9.00012-8, arXiv:2211.05767 [astro-ph.CO].\n- [4] B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, Rept. Prog. Phys. 84 , 116902 (2021), arXiv:2002.12778 [astro-ph.CO].\n- [5] M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama, Class. Quant. Grav. 35 , 063001 (2018), arXiv:1801.05235 [astro-ph.CO].\n- [6] R. Anantua, R. Easther, and J. T. Giblin, Phys. Rev. Lett. 103 , 111303 (2009), arXiv:0812.0825 [astro-ph].\n- [7] A. D. Dolgov and D. Ejlli, Phys. Rev. D 84 , 024028 (2011), arXiv:1105.2303 [astro-ph.CO].\n- [8] R. Dong, W. H. Kinney, and D. Stojkovic, JCAP 10 , 034 (2016), arXiv:1511.05642 [astro-ph.CO].\n- [9] A. Ireland, S. Profumo, and J. Scharnhorst, Phys. Rev. D 107 , 104021 (2023), arXiv:2302.10188 [gr-qc].\n- [10] T. Papanikolaou, V. Vennin, and D. Langlois, JCAP 03 , 053 (2021), arXiv:2010.11573 [astro-ph.CO].  \n- [11] G. Dom'enech, C. Lin, and M. Sasaki, JCAP 04 , 062 (2021), [Erratum: JCAP 11, E01 (2021)], arXiv:2012.08151 [gr-qc].\n- [12] G. Dom'enech, V. Takhistov, and M. Sasaki, Phys. Lett. B 823 , 136722 (2021), arXiv:2105.06816 [astro-ph.CO].\n- [13] T. Papanikolaou, JCAP 10 , 089 (2022), arXiv:2207.11041 [astro-ph.CO].\n- [14] G. Dom'enech, Universe 7 , 398 (2021), arXiv:2109.01398 [gr-qc].\n- [15] C. Caprini, D. G. Figueroa, R. Flauger, G. Nardini, M. Peloso, M. Pieroni, A. Ricciardone, and G. Tasinato, JCAP 11 , 017 (2019), arXiv:1906.09244 [astro-ph.CO].\n- [16] G. Dvali and A. Vilenkin, JCAP 03 , 010 (2004), arXiv:hep-th/0312007.\n- [17] E. J. Copeland, R. C. Myers, and J. Polchinski, JHEP 06 , 013 (2004), arXiv:hep-th/0312067.\n- [18] M. Yamada and K. Yonekura, Phys. Lett. B 838 , 137724 (2023), arXiv:2204.13125 [hep-th].\n- [19] M. Yamada and K. Yonekura, Phys. Rev. D 106 , 123515 (2022), arXiv:2204.13123 [hep-th].\n- [20] M. G. Jackson, N. T. Jones, and J. Polchinski, JHEP 10 , 013 (2005), arXiv:hep-th/0405229.\n- [21] M. Sakellariadou, JCAP 04 , 003 (2005), arXiv:hep-th/0410234.\n- [22] A. Avgoustidis and E. P. S. Shellard, Phys. Rev. D 73 , 041301 (2006), arXiv:astro-ph/0512582.\n- [23] P. P. Avelino and L. Sousa, Phys. Rev. D 85 , 083525 (2012), arXiv:1202.6298 [astro-ph.CO].\n- [24] J. Ellis, M. Lewicki, C. Lin, and V. Vaskonen, Phys. Rev. D 108 , 103511 (2023), arXiv:2306.17147 [astro-ph.CO].\n- [25] A. Vilenkin, Nucl. Phys. B 196 , 240 (1982).\n- [26] J. Preskill and A. Vilenkin, Phys. Rev. D 47 , 2324 (1993), arXiv:hep-ph/9209210.\n- [27] L. Leblond, B. Shlaer, and X. Siemens, Phys. Rev. D 79 , 123519 (2009), arXiv:0903.4686 [astroph.CO].\n- [28] A. Monin and M. B. Voloshin, Phys. Rev. D 78 , 065048 (2008), arXiv:0808.1693 [hep-th].\n- [29] W. Buchmuller, V. Domcke, and K. Schmitz, Phys. Lett. B 811 , 135914 (2020), arXiv:2009.10649 [astro-ph.CO].\n- [30] W. Buchmuller, V. Domcke, and K. Schmitz, JCAP 12 , 006 (2021), arXiv:2107.04578 [hep-ph].\n- [31] W. Buchmuller, V. Domcke, and K. Schmitz, JCAP 11 , 020 (2023), arXiv:2307.04691 [hep-ph].\n- [32] M. A. Masoud, M. U. Rehman, and Q. Shafi, JCAP 11 , 022 (2021), arXiv:2107.09689 [hep-ph].\n- [33] G. Servant and P. Simakachorn, Phys. Rev. D 109 , 103538 (2024), arXiv:2312.09281 [hep-ph].\n- [34] A. Ghoshal, Y. Gouttenoire, L. Heurtier, and P. Simakachorn, JHEP 08 , 196 (2023), arXiv:2304.04793 [hep-ph].\n- [35] D. Borah, S. Jyoti Das, R. Roshan, and R. Samanta, Phys. Rev. D 108 , 023531 (2023), arXiv:2304.11844 [hep-ph].\n- [36] S. Antusch, K. Hinze, S. Saad, and J. Steiner, Phys. Lett. B 856 , 138924 (2024), arXiv:2405.03746 [hep-ph].\n- [37] R. H. Cyburt, B. D. Fields, K. A. Olive, and T.-H. Yeh, Rev. Mod. Phys. 88 , 015004 (2016), arXiv:1505.01076 [astro-ph.CO].\n- [38] M. Kawasaki, K. Kohri, and N. Sugiyama, Phys. Rev. Lett. 82 , 4168 (1999), arXiv:astro-ph/9811437.\n- [39] T. Hasegawa, N. Hiroshima, K. Kohri, R. S. L. Hansen, T. Tram, and S. Hannestad, JCAP 12 , 012 (2019), arXiv:1908.10189 [hep-ph].\n- [40] G. Agazie et al. (NANOGrav), Astrophys. J. Lett. 951 , L8 (2023), arXiv:2306.16213 [astro-ph.HE].\n- [41] J. Antoniadis et al. (EPTA, InPTA:), Astron. Astrophys. 678 , A50 (2023), arXiv:2306.16214 [astroph.HE].\n- [42] D. J. Reardon et al. , Astrophys. J. Lett. 951 , L6 (2023), arXiv:2306.16215 [astro-ph.HE].\n- [43] H. Xu et al. , Res. Astron. Astrophys. 23 , 075024 (2023), arXiv:2306.16216 [astro-ph.HE].\n- [44] A. Afzal et al. (NANOGrav), Astrophys. J. Lett. 951 , L11 (2023), arXiv:2306.16219 [astro-ph.HE].\n- [45] J. Antoniadis et al. (EPTA), (2023), arXiv:2306.16227 [astro-ph.CO].\n- [46] S. Blasi, V. Brdar, and K. Schmitz, Phys. Rev. Lett. 126 , 041305 (2021), arXiv:2009.06607 [astroph.CO].\n- [47] J. Ellis and M. Lewicki, Phys. Rev. Lett. 126 , 041304 (2021), arXiv:2009.06555 [astro-ph.CO].\n- [48] R. Samanta and S. Datta, JHEP 05 , 211 (2021), arXiv:2009.13452 [hep-ph].\n- [49] S. Datta, A. Ghosal, and R. Samanta, JCAP 08 , 021 (2021), arXiv:2012.14981 [hep-ph].\n- [50] T. W. B. Kibble, J. Phys. A 9 , 1387 (1976).\n- [51] M. B. Hindmarsh and T. W. B. Kibble, Rept. Prog. Phys. 58 , 477 (1995), arXiv:hep-ph/9411342.  \n- [52] R. Jeannerot, J. Rocher, and M. Sakellariadou, Phys. Rev. D 68 , 103514 (2003), arXiv:hepph/0308134.\n- [53] D. G. Figueroa, M. Pieroni, A. Ricciardone, and P. Simakachorn, Phys. Rev. Lett. 132 , 171002 (2024), arXiv:2307.02399 [astro-ph.CO].\n- [54] J. Ellis, M. Fairbairn, G. Franciolini, G. Hutsi, A. Iovino, M. Lewicki, M. Raidal, J. Urrutia, V. Vaskonen, and H. Veermae, Phys. Rev. D 109 , 023522 (2024), arXiv:2308.08546 [astro-ph.CO].\n- [55] D. G. Figueroa, M. Pieroni, A. Ricciardone, and P. Simakachorn, Phys. Rev. Lett. 132 , 171002 (2024), arXiv:2307.02399 [astro-ph.CO].\n- [56] S. Antusch, K. Hinze, S. Saad, and J. Steiner, Phys. Rev. D 108 , 095053 (2023), arXiv:2307.04595 [hep-ph].\n- [57] S. Antusch, K. Hinze, and S. Saad, (2024), arXiv:2406.17014 [hep-ph].\n- [58] G. Lazarides, R. Maji, and Q. Shafi, Phys. Rev. D 108 , 095041 (2023), arXiv:2306.17788 [hep-ph].\n- [59] R. Maji and W.-I. Park, JCAP 01 , 015 (2024), arXiv:2308.11439 [hep-ph].\n- [60] D. Marfatia and Y.-L. Zhou, JHEP 07 , 204 (2024), arXiv:2312.10455 [hep-ph].\n- [61] R. Abbott et al. (KAGRA, Virgo, LIGO Scientific), Phys. Rev. D 104 , 022004 (2021), arXiv:2101.12130 [gr-qc].\n- [62] J. Aasi et al. (LIGO Scientific), Class. Quant. Grav. 32 , 074001 (2015), arXiv:1411.4547 [gr-qc].\n- [63] B. Sathyaprakash et al. , Class. Quant. Grav. 29 , 124013 (2012), [Erratum: Class.Quant.Grav. 30, 079501 (2013)], arXiv:1206.0331 [gr-qc].\n- [64] T. Fujita, M. Kawasaki, K. Harigaya, and R. Matsuda, Phys. Rev. D 89 , 103501 (2014), arXiv:1401.1909 [astro-ph.CO].\n- [65] O. Lennon, J. March-Russell, R. Petrossian-Byrne, and H. Tillim, JCAP 04 , 009 (2018), arXiv:1712.07664 [hep-ph].\n- [66] D. Hooper, G. Krnjaic, and S. D. McDermott, JHEP 08 , 001 (2019), arXiv:1905.01301 [hep-ph].\n- [67] L. Morrison, S. Profumo, and Y. Yu, JCAP 05 , 005 (2019), arXiv:1812.10606 [astro-ph.CO].\n- [68] R. Samanta and F. R. Urban, JCAP 06 , 017 (2022), arXiv:2112.04836 [hep-ph].\n- [69] D. Borah, S. Jyoti Das, R. Samanta, and F. R. Urban, JHEP 03 , 127 (2023), arXiv:2211.15726 [hep-ph].\n- [70] A. Vilenkin, Phys. Lett. B 107 , 47 (1981).\n- [71] T. Vachaspati and A. Vilenkin, Phys. Rev. D 31 , 3052 (1985).\n- [72] C. T. Hill, H. M. Hodges, and M. S. Turner, Phys. Rev. D 37 , 263 (1988).\n- [73] W. T. Emond, S. Ramazanov, and R. Samanta, JCAP 01 , 057 (2022), arXiv:2108.05377 [hep-ph].\n- [74] M. Chianese, S. Datta, G. Miele, R. Samanta, and N. Saviano, (2024), arXiv:2406.01231 [hep-ph].\n- [75] J. J. Blanco-Pillado, K. D. Olum, and B. Shlaer, Phys. Rev. D 89 , 023512 (2014), arXiv:1309.6637 [astro-ph.CO].\n- [76] J. J. Blanco-Pillado and K. D. Olum, Phys. Rev. D 96 , 104046 (2017), arXiv:1709.02693 [astro-ph.CO].\n- [77] C. J. A. P. Martins and E. P. S. Shellard, Phys. Rev. D 54 , 2535 (1996), arXiv:hep-ph/9602271.\n- [78] C. J. A. P. Martins and E. P. S. Shellard, Phys. Rev. D 65 , 043514 (2002), arXiv:hep-ph/0003298.\n- [79] L. Sousa and P. P. Avelino, Phys. Rev. D 88 , 023516 (2013), arXiv:1304.2445 [astro-ph.CO].\n- [80] P. Auclair et al. , JCAP 04 , 034 (2020), arXiv:1909.00819 [astro-ph.CO].\n- [81] T. Damour and A. Vilenkin, Phys. Rev. D 64 , 064008 (2001), arXiv:gr-qc/0104026.\n- [82] A. Vilenkin, Phys. Rev. D 43 , 1060 (1991).\n- [83] D. Matsunami, L. Pogosian, A. Saurabh, and T. Vachaspati, Phys. Rev. Lett. 122 , 201301 (2019), arXiv:1903.05102 [hep-ph].\n- [84] P. Auclair, D. A. Steer, and T. Vachaspati, Phys. Rev. D 101 , 083511 (2020), arXiv:1911.12066 [hep-ph].\n- [85] L. Sousa, P. P. Avelino, and G. S. F. Guedes, Phys. Rev. D 101 , 103508 (2020), arXiv:2002.01079 [astro-ph.CO].\n- [86] Y. Cui, M. Lewicki, D. E. Morrissey, and J. D. Wells, JHEP 01 , 081 (2019), arXiv:1808.08968 [hep-ph].\n- [87] S. Blasi, V. Brdar, and K. Schmitz, Phys. Rev. Res. 2 , 043321 (2020), arXiv:2004.02889 [hep-ph].\n- [88] J. J. Blanco-Pillado, K. D. Olum, and X. Siemens, Phys. Lett. B 778 , 392 (2018), arXiv:1709.02434 [astro-ph.CO].\n- [89] A. Peimbert, M. Peimbert, and V. Luridiana, Rev. Mex. Astron. Astrofis. 52 , 419 (2016), arXiv:1608.02062 [astro-ph.CO].  \n- [90] E. Babichev, D. Gorbunov, S. Ramazanov, R. Samanta, and A. Vikman, Phys. Rev. D 108 , 123529 (2023), arXiv:2307.04582 [hep-ph].\n- [91] S. Balaji, G. Dom'enech, and G. Franciolini, JCAP 10 , 041 (2023), arXiv:2307.08552 [gr-qc].\n- [92] A. Weltman et al. , Publ. Astron. Soc. Austral. 37 , e002 (2020), arXiv:1810.02680 [astro-ph.CO].\n- [93] P. Amaro-Seoane et al. (LISA), (2017), arXiv:1702.00786 [astro-ph.IM].\n- [94] S. Kawamura et al. , Class. Quant. Grav. 23 , S125 (2006).\n- [95] G. Agazie et al. (International Pulsar Timing Array), Astrophys. J. 966 , 105 (2024), arXiv:2309.00693 [astro-ph.HE].\n- [96] G. F. Giudice, E. W. Kolb, and A. Riotto, Phys. Rev. D 64 , 023508 (2001), arXiv:hep-ph/0005123.\n- [97] I. Masina, Eur. Phys. J. Plus 135 , 552 (2020), arXiv:2004.04740 [hep-ph].\n- [98] J. H. MacGibbon, Phys. Rev. D 44 , 376 (1991).\n- [99] M. Maggiore, Phys. Rept. 331 , 283 (2000), arXiv:gr-qc/9909001.\n- [100] K. Schmitz, JHEP 01 , 097 (2021), arXiv:2002.04615 [hep-ph].\n- [101] Y. F. Perez-Gonzalez and J. Turner, Phys. Rev. D 104 , 103021 (2021), arXiv:2010.03565 [hep-ph].\n- [102] N. Bernal, C. S. Fong, Y. F. Perez-Gonzalez, and J. Turner, Phys. Rev. D 106 , 035019 (2022), arXiv:2203.08823 [hep-ph].\n- [103] A. Cheek, L. Heurtier, Y. F. Perez-Gonzalez, and J. Turner, Phys. Rev. D 105 , 015023 (2022), arXiv:2107.00016 [hep-ph].\n- [104] P. Gondolo, P. Sandick, and B. Shams Es Haghi, Phys. Rev. D 102 , 095018 (2020), arXiv:2009.02424 [hep-ph].\n- [105] G. Franciolini and P. Pani, Phys. Rev. D 108 , 083527 (2023), arXiv:2304.13576 [astro-ph.CO].\n- [106] B. Barman, D. Borah, S. Das Jyoti, and R. Roshan, JCAP 08 , 068 (2022), arXiv:2204.10339 [hep-ph].\n- [107] B. Barman, D. Borah, S. Jyoti Das, and R. Roshan, Phys. Rev. D 107 , 095002 (2023), arXiv:2212.00052 [hep-ph].\n- [108] D. Borah, S. Jyoti Das, and R. Roshan, Nucl. Phys. B 1002 , 116528 (2024), arXiv:2208.04965 [hepph].\n- [109] B. Barman, S. Jyoti Das, M. R. Haque, and Y. Mambrini, Phys. Rev. D 110 , 043528 (2024), arXiv:2403.05626 [hep-ph].\n- [110] D. Borah, N. Das, S. Jyoti Das, and R. Samanta, (2024), arXiv:2403.02401 [hep-ph].\n- [111] R. Calabrese, M. Chianese, J. Gunn, G. Miele, S. Morisi, and N. Saviano, Phys. Rev. D 109 , 103001 (2024), arXiv:2311.13276 [hep-ph].\n- [112] K. Schmitz and X.-J. Xu, Phys. Lett. B 849 , 138473 (2024), arXiv:2311.01089 [hep-ph].\n- [113] K.-Y. Choi, J. Kim, and E. Lkhagvadorj, JCAP 02 , 020 (2024), arXiv:2309.16122 [hep-ph].\n- [114] R. Calabrese, M. Chianese, J. Gunn, G. Miele, S. Morisi, and N. Saviano, Phys. Rev. D 107 , 123537 (2023), arXiv:2305.13369 [hep-ph].\n- [115] T. C. Gehrman, B. Shams Es Haghi, K. Sinha, and T. Xu, JCAP 10 , 001 (2023), arXiv:2304.09194 [hep-ph].\n- [116] N. Bhaumik, A. Ghoshal, and M. Lewicki, JHEP 07 , 130 (2022), arXiv:2205.06260 [astro-ph.CO].\n- [117] N. Bernal, Y. F. Perez-Gonzalez, Y. Xu, and O. Zapata, Phys. Rev. D 104 , 123536 (2021), arXiv:2110.04312 [hep-ph].\n- [118] Y.-P. Wu, E. Pinetti, and J. Silk, Phys. Rev. Lett. 128 , 031102 (2022), arXiv:2109.09875 [astroph.CO].\n- [119] P. Sandick, B. S. Es Haghi, and K. Sinha, Phys. Rev. D 104 , 083523 (2021), arXiv:2108.08329 [astro-ph.CO].\n- [120] N. Bernal, F. Hajkarim, and Y. Xu, Phys. Rev. D 104 , 075007 (2021), arXiv:2107.13575 [hep-ph].\n- [121] A. Cheek, L. Heurtier, Y. F. Perez-Gonzalez, and J. Turner, Phys. Rev. D 105 , 015022 (2022), arXiv:2107.00013 [hep-ph].\n- [122] I. Masina, Grav. Cosmol. 27 , 315 (2021), arXiv:2103.13825 [gr-qc].\n- [123] A. Ambrosone, R. Calabrese, D. F. G. Fiorillo, G. Miele, and S. Morisi, Phys. Rev. D 105 , 045001 (2022), arXiv:2106.11980 [hep-ph].\n- [124] T. Kitabayashi, Int. J. Mod. Phys. A 36 , 2150139 (2021), arXiv:2101.01921 [hep-ph].\n- [125] M. Y. Khlopov, A. Barrau, and J. Grain, Class. Quant. Grav. 23 , 1875 (2006), arXiv:astroph/0406621.\n- [126] M. Gell-Mann, P. Ramond, and R. Slansky, Conf. Proc. C 790927 , 315 (1979), arXiv:1306.4669 [hep-th].  \n- [127] M. Fukugita and T. Yanagida, Phys. Lett. B 174 , 45 (1986).\n- [128] R. Saito and J. Yokoyama, Phys. Rev. Lett. 102 , 161101 (2009), [Erratum: Phys.Rev.Lett. 107, 069901 (2011)], arXiv:0812.4339 [astro-ph].\n- [129] K. Inomata, M. Kawasaki, K. Mukaida, T. Terada, and T. T. Yanagida, Phys. Rev. D 101 , 123533 (2020), arXiv:2003.10455 [astro-ph.CO].\n- [130] R.-g. Cai, S. Pi, and M. Sasaki, Phys. Rev. Lett. 122 , 201101 (2019), arXiv:1810.11000 [astro-ph.CO].\n- [131] D. Hooper, G. Krnjaic, J. March-Russell, S. D. McDermott, and R. Petrossian-Byrne, (2020), arXiv:2004.00618 [astro-ph.CO].\n- [132] N. Seto, S. Kawamura, and T. Nakamura, Phys. Rev. Lett. 87 , 221103 (2001), arXiv:astroph/0108011.\n- [133] N. Aggarwal, G. P. Winstone, M. Teo, M. Baryakhtar, S. L. Larson, V. Kalogera, and A. A. Geraci, Phys. Rev. Lett. 128 , 111101 (2022), arXiv:2010.13157 [gr-qc].\n- [134] N. Aggarwal et al. , Living Rev. Rel. 24 , 4 (2021), arXiv:2011.12414 [gr-qc].\n- [135] A. S. Chou et al. (Holometer), Phys. Rev. D 95 , 063002 (2017), arXiv:1611.05560 [astro-ph.IM].\n- [136] V. Domcke and C. Garcia-Cely, Phys. Rev. Lett. 126 , 021104 (2021), arXiv:2006.01161 [astro-ph.CO].\n- [137] T. Akutsu et al. , Phys. Rev. Lett. 101 , 101101 (2008), arXiv:0803.4094 [gr-qc].\n- [138] A. Ito, T. Ikeda, K. Miuchi, and J. Soda, Eur. Phys. J. C 80 , 179 (2020), arXiv:1903.04843 [gr-qc].\n- [139] A. Ito and J. Soda, Eur. Phys. J. C 80 , 545 (2020), arXiv:2004.04646 [gr-qc].\n- [140] V. Domcke, C. Garcia-Cely, and N. L. Rodd, Phys. Rev. Lett. 129 , 041101 (2022), arXiv:2202.00695 [hep-ph].\n- [141] A. Ejlli, D. Ejlli, A. M. Cruise, G. Pisano, and H. Grote, Eur. Phys. J. C 79 , 1032 (2019), arXiv:1908.00232 [gr-qc].\n- [142] A. Berlin, D. Blas, R. Tito D'Agnolo, S. A. R. Ellis, R. Harnik, Y. Kahn, and J. Schutte-Engel, Phys. Rev. D 105 , 116011 (2022), arXiv:2112.11465 [hep-ph].\n- [143] Y. Akrami et al. (Planck), Astron. Astrophys. 641 , A10 (2020), arXiv:1807.06211 [astro-ph.CO].\n- [144] N. Herman, L. Lehoucq, and A. F'uzfa, Phys. Rev. D 108 , 124009 (2023), arXiv:2203.15668 [gr-qc].\n- [145] I. Esteban, M. C. Gonzalez-Garcia, M. Maltoni, T. Schwetz, and A. Zhou, JHEP 09 , 178 (2020), arXiv:2007.14792 [hep-ph].\n- [146] T. Asaka, S. Blanchet, and M. Shaposhnikov, Phys. Lett. B 631 , 151 (2005), arXiv:hep-ph/0503065.\n- [147] T. Asaka, S. Blanchet, and M. Shaposhnikov, Phys. Lett. B 631 , 151 (2005), arXiv:hep-ph/0503065.\n- [148] S. Dodelson and L. M. Widrow, Phys. Rev. Lett. 72 , 17 (1994), arXiv:hep-ph/9303287.\n- [149] X.-D. Shi and G. M. Fuller, Phys. Rev. Lett. 82 , 2832 (1999), arXiv:astro-ph/9810076.\n- [150] P. Di Bari, P. O. Ludl, and S. Palomares-Ruiz, JCAP 11 , 044 (2016), arXiv:1606.06238 [hep-ph].\n- [151] P. Di Bari, K. Farrag, R. Samanta, and Y. L. Zhou, JCAP 10 , 029 (2020), arXiv:1908.00521 [hep-ph].\n- [152] A. Datta, R. Roshan, and A. Sil, Phys. Rev. Lett. 127 , 231801 (2021), arXiv:2104.02030 [hep-ph].\n- [153] A. Mitridate, D. Wright, R. von Eckardstein, T. Schroder, J. Nay, K. Olum, K. Schmitz, and T. Trickle, (2023), arXiv:2306.16377 [hep-ph]."}