Cosmic superstrings, metastable strings and ultralight primordial black holes: from NANOGrav to LIGO and beyond
Satyabrata Datta Rome Samanta

=1


 
        











 

 
 =10000
=10000





[],n,






    addtoresetequationsectionsection
        
        
        amisatyabrata703@gmail.com
         samanta@na.infn.it§ ABSTRACT 
           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 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.
           
        Cosmic superstrings, metastable strings and ultralight primordial black holes: from NANOGrav to LIGO and beyond
            Rome Samanta
            
        ===============================================================================================================
        
        
        
        § INTRODUCTION
        
        Ultralight primordial black holes PBHs <cit.> with initial mass M_ BH≲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 <cit.>. 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 <cit.> and produce density fluctuation to induce GW background testable in the planned GW detectors <cit.>. 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 <cit.> 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 <cit.> and metastable strings <cit.>[Effect of PBHs on GWs from stable cosmic strings have been recently explored in Refs.<cit.>]. 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., <cit.>.
        
        The 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) <cit.>, 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 the 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)  <cit.>. Unlike stable cosmic strings <cit.> 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  <cit.>, cosmic superstrings and metastable cosmic strings fit the data better <cit.> 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 <cit.>. 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 <cit.>.  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. 
        
        Focusing 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∈[10^6g,10^8 g] 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) <cit.> could probe the PBH-imprinted 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 <cit.> as large as ∼ 10^5. 
        
        In the companion theory paper by NANOGrav collaboration <cit.>, 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∈[10^6g,10^8 g], 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.<cit.> (subsequent references are in Sec.<ref>). For the M_ BH∈[10^6g,10^8 g], 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 <cit.> accommodating neutrino mass, BAU via leptogenesis, 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. 
        
        
        
        
        
        
        
        
        § 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.
         Stable 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)<cit.>. Long strings are described by a correlation length L=√(μ/ρ_∞), with ρ_∞ is the long string energy density and μ is the string tension defined as <cit.>μ=π v_Φ^2 h(λ,g^'),h(λ,g^')≃ 1, unless the coupling λ and g^' are not strongly hierarchical <cit.>. 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 Γ≃ 50 <cit.>, α≃ 0.1 <cit.>, 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=2k/l_k=a(t_0)/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 <cit.>Ω_ GW= ∑_k=1^k_ max2k ℱ_α Gμ^2 Γ_k/fρ_c∫_t_i^t_0[a(t)/a(t_0)]^5 n_ω(t,l_k) dt  ,
        where ρ_c is the critical energy density of the universe, ℱ_α≃ 0.1 is an efficiency factor <cit.> and n_ω(t,l_k) is the loop number density which can be computed from the velocity-dependent-one-scale model as <cit.>
        n_ω(t,l_k)=A_β/α(α+Γ G μ)^3(1-β)/[l_k(t)+Γ G μ t]^4-3βt^3β . 
        In Eq.(<ref>), β=2/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 <cit.>. The quantity Γ_k=Γ k^-δ/ζ(δ) quantifies the emitted power in k-th mode, with δ=4/3 (δ=5/3) for loops containing small-scale structures such as cusps (kinks) <cit.> 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., PTA 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.
        
        The integral in Eq. (<ref>) is subjected to the Heaviside functions Θ(t_i-t_ fric)Θ(t_i-l_c/α) 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 <cit.>. For l>l_c, the GW emission dominates over particle production as shown by recent numerical simulations <cit.>. The critical length can be approximately computed as l_c≃δ_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. 
        
          Without any intermediate matter-dominated epoch, Eq.(<ref>) 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.<ref>) which for k=1 mode is given by <cit.>Ω_ GW^ plt = 128πℱ_α Gμ/9ζ(δ)A_R/ϵ_RΩ_R[(1+ϵ_R)^3/2-1]  . 
         In Eq. (<ref>), 
        ϵ_R=α/Γ Gμ≫ 1, A_R≡ A_β≃ 5.4 and Ω_R∼ 9× 10^-5. Note that the plateau amplitude Ω_ GW^ plt∝√(μ)∝ 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.<ref>). This spectral break frequency f_ brk can be estimated as  <cit.>
        f_ brk≃ 0.45  Hz(10^-12/Gμ)^1/2(T_ brk/ GeV) ,
        where 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., Ω_ GW^ plt,∞≃Ω_ GW^ pltζ(δ) and the spectrum falls as f^1-δ instead of f^-1<cit.>.
        
        Given 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/T_ brk>1, where we derive T_ cusp  as 
        
        T_ cusp≃ 2.3× 10^4  GeV(G μ/10^-12)^5/4 .
        
        Similar to Eq.(<ref>), an expression for f_ cusp can therefore be obtained as 
        
        f_ cusp≃ 1.01× 10^4  Hz(G μ/10^-12)^3/4 
        which we find to be consistent with Eq.(50) of Ref.<cit.>. While analyzing the GW spectrum with PBHs, one has to make sure  that f_ cusp>f_ brk is satisfied (see Fig.<ref>, top panel).
         Cosmic 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<cit.> (see Ref.<cit.> for a different scaling relation). This method has also been adopted by the NANOGrav collaboration <cit.>. Therefore, we shall use all the methodologies considered for stable strings but rescaling the loop number density by P^-1, with P∈[10^-3,1]<cit.>. Note that, this is not a rescaling of Gμ that one obtains for hierarchical coupling (h(λ,g^')≠ 1) <cit.>, thereby making the spectral shape, e.g., in the PTA band, different from the stable cosmic strings.  
         Metastable cosmic strings: Metastable cosmic strings decay into string segments connecting monopole-antimonopole 
        pairs by quantum tunneling <cit.>. In the semiclassical approximation, the network is characterized  by the decay rate per string unit length as 
        Γ_d = μ/2πexp(-πκ),  √(κ) = m_M/μ^1/2≃Λ_ GUT/v_Φ,
        where m_M  is the monopole mass which is of the order of Λ_ GUT≃ 10^16 GeV. From Eq.(<ref>), one can define a time t_s=1/Γ_d ^1/2, 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.(<ref>) receives a correction factor <cit.>:
        
        D(l,t)= e^-Γ_d[l(t-t_s)+1/2Γ Gμ(t-t_s)^2]Θ(α t_s-l(t_s)),
        
        where the Θ function with l(t_s)≃ 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.(<ref>) 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 <cit.>. 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.<ref>). As such, for sufficiently small values of κ, the spectrum reaches the plateau at a very high frequency f≳ 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 <cit.>. ii) The infra-red tail of the GW spectrum behaves as Ω_ GW∝  f^2<cit.> 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., <cit.>. 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. 
        
        
        We conclude the discussion by showing the GW spectra in Fig.<ref> 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.<ref> (left) we show the GW spectrum for cosmic super (stable) strings with the green solid (dashed) curve. For the superstrings, we choose the benchmark values P=6× 10^-3 and Gμ=2× 10^-12 while for the stable cosmic strings, we rescale Gμ as Gμ→ Gμ/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 blue-tilted 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) <cit.> ). 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[With the blue curve we show the GW spectra for √(κ)=6 evading all the bounds with the dominant contribution at very high-frequencies. Although this spectrum is not so relevant for the present study, we find it interesting to mention it as it is generally not present in the case of stable (and super) strings.]. In both plots, the black dotted lines represent the analytical expression for the plateau amplitude and the spectral break frequency given in Eq.(<ref>) and (<ref>), 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.
        
        An 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.<cit.> 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. 
        
        
        
        
        § 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 <cit.>:
        dρ_R/dz+4/zρ_R   =   0,
        dρ_ BH/dz+3/zH/H̃ρ_ BH   -   Ṁ_ BH/M_ BH1/zH̃ρ_ BH=0,
        where H is the Hubble parameter (H_0∼ 1.44× 10^-42 GeV),  z=T_ Bf/T, with T_ Bf being the PBH formation temperature and the time derivative Ṁ_ BH is given by
        Ṁ_ BH=𝒢g_*B(T_ BH)/30720πM_Pl^4/M_ BH^2,
        where M_Pl=1.22× 10^19 GeV, 𝒢≃ 3.8 is a grey body factor and g_*B(T_ BH)≃ 100 being the number of relativistic degrees of freedom below the PBH temperature T_ BH with Standard Model (SM) containing left-handed neutrinos <cit.>. The quantity H̃ and the scale factor a evolve as
        H̃=(H+𝒦), da/dz=(1-𝒦/H̃)a/z,
        
        where 𝒦=Ṁ_ BH/M_ BHρ_ BH/4ρ_R. To derive Eq.(<ref>)-Eq.(<ref>), 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 
        
        T_ ev=(5 M_Pl^2/π^3g_*(T_ ev)τ^2)^1/4,
        
        where τ is the PBH lifetime which can be computed by integrating Eq.(<ref>). This allows also to express T_ ev in terms of M_ BH as
        
        T_ ev=0.1(M_ BH/4.1× 10^7 g)^-3/2  GeV.
        Eq.(<ref>)  shows that with T_ ev≡ T_ brk, all the plots in Fig.<ref> can be reproduced for PBHs with initial mass M_ BH≃ 4.1× 10^7  g. 
        
        
        
        It is crucial to note the following: i) An eternal falling of the GW spectrum in the range f_ brk≲ f≲ f_ LIGO, (see Fig.<ref>), 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 <cit.>, 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.<ref>. ii) The PBHs should evaporate before the BBN at a temperature T≳ 10 MeV. Therefore, from Eq.(<ref>), M_ BH is generically constrained as M_ BH≲ 2× 10^8 g, which we use as the absolute upper bound on M_ BH. 
        
        We 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.(<ref>) and Eq.(<ref>), we obtain the spectral break frequency as
        
        f_ brk^ PBH≃ 0.045  Hz(10^-12/Gμ)^1/2(M_ BH/4.1× 10^7 g)^-3/2 .
        When a large number of modes are summed, GW energy density beyond f_ brk^ PBH can be written as 
        Ω_ GW(f)≃Ω_ GW^ plt,∞(f/f_ brk^ PBH)^-1/3Θ[ f-f_ brk^ PBH],
        where 
        Ω_ GW^ plt,∞≃ 3.1× 10^-8 Υ^-1(Gμ/10^-7)^1/2,  Υ=P   for    superstrings
        Υ=1   for    metastable  strings.
        Using Eq.(<ref>), Eq.(<ref>) can be re-expressed as 
        Ω_ GW(f)≃Ω_ GW^ plt,∞Ω_ GW^ mod,PBH,
        
        where Ω_ GW^ mod,PBH, given by
        Ω_ GW^ mod,PBH≃ 0.1 (10^-12/Gμ)^1/6(M_ BH/4.1× 10^7 g)^-1/2(f/25 Hz)^-1/3
         captures the modulation in the GW spectra because of PBH domination. The absence of PBHs corresponds to Ω_ GW^ mod, PBH=1. In which case, Eq.(<ref>) and the non-observation of SGWB by LIGO-O3 run; i.e., Ω_ GW^ any model≲ 1.7× 10^-8 at f_ LIGO=25 Hz, implies the following: 
         without  the  PBH,  LIGO-O3 excludes  
            P≲ 1.8 (Gμ/10^-7)^1/2   for    superstrings,
        
            Gμ≳ 3×10^-8   for    metastable strings.
        Though Eq.(<ref>) 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):
        
        
        Category A: Ω_ GW^ plt,∞≳Ω_ GW^ LIGO-O3 and Ω_ GW^ LIGO-O3≳Ω_ GW^ plt,∞Ω_ GW^ mod,PBH(f=25  Hz)≳Ω_ GW^ LVK-D.
        
        
        Category B: Ω_ GW^ LIGO-O3≳Ω_ GW^ plt,∞≳Ω_ GW^ LVK-D and Ω_ GW^ LIGO-O3≳Ω_ GW^ plt,∞Ω_ GW^ mod,PBH(f=25  Hz) ≳Ω_ GW^ LVK-D.
        
        Category C: Ω_ GW^ LVK-D≳Ω_ GW^ plt,∞Ω_ GW^ mod,PBH(f=25  Hz) for all Ω_ GW^ plt,∞.
        
        
        Here Ω_ GW^ LVK-D=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.(<ref>), 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_ brk^ PBH≲25  Hz and Ω_ GW^ LIGO-O3≳Ω_ GW^ plt,∞Ω_ GW^ mod,PBH(f=25  Hz) ≳Ω_ GW^ LVK-D. 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 at ET. In Fig.<ref>, 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 <ref> (Superstrings) and Table <ref> (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. 
        
        
         §.§ Numerical results: PBH mass ranges for PTA GW signals evading LIGO-O3 bound 
        
        On the top panel in Fig.<ref>, 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<cit.> (one may also look at the Appendix <ref> of this article). The benchmark points tabulated in Table <ref> corresponding to the GW spectra in Fig.<ref> (left) are denoted with the symbols   BP1 (Cat-A),   BP2 (Cat-B), and 2  BP3 (Cat-C). A similar description holds for metastable strings on the Gμ and √(κ) plane (Fig.<ref>, top) with the benchmarks tabulated in Table <ref> and corresponding GW spectra shown in Fig.<ref> (right). In the bottom panel of both figures, the GW spectra presented in Fig.<ref> 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∼[10^6 g, 10^8 g] 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 Ω_ GW^ LVK-D with f_ brk  at lower frequencies. On the other hand, for Category B, only lighter PBHs are allowed for the opposite reason, i.e., Ω_ GW^ LVK-D≲Ω_ GW(f=25  Hz) requiring f_ brk≲ 25  Hz 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 <cit.>. As such, from Fig.<ref>, 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 <cit.>. Therefore, one can produce Category-B signals for metastable strings with  M_ BH∼ 10^5.5 g. 
        
        
         §.§ 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) <cit.>.  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 as <cit.> SNR = √(t_ obs∫_f_ min^f_ maxdf (Ω_ GW(f)/Ω_ noise(f))^2) ,
        where Ω_ noise(f) is the noise spectrum of the detector and f_ min(f_ max) are minimum (maximum) accessible frequency. In Fig.<ref>, 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 owing 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. 
        
        
        Given 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. 
        
        
        
        §  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.<ref> and Fig.<ref>, 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 <cit.>. 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 <cit.> based on the seesaw mechanism <cit.> of neutrino masses stemming also the possibility of baryogenesis via leptogenesis <cit.> and superheavy dark matter in the presence of PBHs<cit.>. In the scenario, it was reported for the first time, that if there are ultralight PBHs in the early universe, seesaw Lagrangian can accommodate the novel right-handed neutrino (RHN) mass spectrum Λ_ GUT∼ M_R3>M_R2>M_R1, where M_Ri being the mass of the ith RHN, with heaviest one (N_R3) being super heavy dark matter and the other two generating baryon asymmetry. The analysis in Ref.<cit.>  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)_B-L^ PBH seesaw model. 
        
        
         §.§ 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., <cit.>), ultralight PBHs could possibly merge, releasing  GWs <cit.>, PBHs emit gravitons, which constitute GWs at ultra-high frequencies <cit.>, and the inhomogeneous distribution of PBHs, which results in density fluctuations, triggers the production of GWs <cit.>. 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.<cit.>. 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
        
            Ω_ GW^ PBH,I(t_0,f)≃Ω^ peak,I_ GW( f/f_ peak)^11/3Θ (f_ peak,I-f),
        
        where the peak amplitude is given by
        
            Ω^ peak,I_ GW≃ 2× 10^-6(β/10^-8)^16/3(M_ BH/10^7 g)^34/9,
        
        and the peak frequency is given by
        
            f_ peak,I≃ 5.4  Hz( M_ BH/10^7 g)^-5/6.
        
        The ultra-violet cutoff frequency f_ peak corresponds to the comoving scale representing the mean separation of PBHs at the formation time.
        
        On the other hand, the present spectral energy density of GWs from gravitons is given by <cit.>
            Ω_ GW(f) h^2= (2π)^4βγ^2f^4/4H_0^2M_ BH∫_t_ Bf^t_ ev(1-t-t_ Bf/τ_ ev)^2/3( a(t_ Bf)/a(t_0))^-3/exp( 16π^2 f a(t_0)M_ BH/a(t)M_ Pl^2)-1
        
        By maximizing the Eq.(<ref>) w.r.t the present-day frequency f, one obtains the peak frequency and consequently the peak GW amplitude as 
        
        
        f_ peak,g≃ 2.8 a(t_ ev)T_ BH≃  2 × 10^17 Hz( M_ BH/10^7 g)^1/2,  Ω_ GW^ peak,g∼ 10^-7.
        In Fig.<ref> and Fig.<ref> (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.(<ref>)) 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≲ 5.6× 10^-6Δ N_ eff, with Δ N_ eff≲ 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.<ref> 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.<ref>. 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.<ref>, 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.<ref> and Fig.<ref> extracted for Ω_ GW∝ f^-1/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.<ref> and Fig.<ref>, the duration of PBH domination would be such that from f_ brk until the f_ LIGO, the spectrum would fall as Ω_ GW∝ f^-1/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. 
        
        Besides 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.<ref> and <ref>, where the GW strain h_c is related to Ω_ GW as
        
        h_c(f)=√(3 H_0^2/2π^2Ω_ GW)f^-1.
        
        Because Eq.(<ref>) characterizes the GW background, for the allowed mass ranges of PBH presented in Sec.<ref>, 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 <cit.> are not presently competitive with the Δ N_ eff bound from Planck <cit.> (according to Ref.<cit.>, the electromagnetic resonant cavities can potentially reach below
        the dark radiation projections around f∼ 10^8 Hz; the black oscillating curve in Fig.<ref> and Fig.<ref>), 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. 
        
        
         §.§ Implication on a realistic BSM model of leptogenesis and super-heavy dark matter
        
        In the seesaw mechanism <cit.>, the light neutrino mass m_i∼ 0.01 eV inferred by neutrino oscillation data and cosmological measurements <cit.>, can be explained as a ratio of two scales: m_i≃Λ_ EW^2/Λ_ GUT, with Λ_ EW and Λ_ GUT being the Electroweak and GUT scale, respectively. At the Lagrangian level, these two scales are introduced as ℒ_ seesaw∼Λ_ EWL N_R+Λ_ GUTN^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., <cit.>. Such attempts are theoretically well-built and predictive in particle physics experiments, however, at the expense of fine-tuning the L N_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 <cit.> consider N_R to be extremely light 𝒪( GeV,  keV), therefore tiny L N_R coupling. In Ref.<cit.>, 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≃Λ_ EW^2/Λ_ GUT, i.e., considering all the RHNs to be 𝒪(Λ_ 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 <cit.>. 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 
        
        M_ DM≃ 4.7×10^12 (M_ BH/10^7 g)^-5/2  GeV.
        The 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≃ 6.3× 10^-10<cit.> is taken into account, the parameter β and scale of thermal leptogenesis get correlated as <cit.>β≃ 1.2× 10^-9(10^7g/M_ BH) (T_ lepto/10^14 GeV),
        where we have assumed leptogenesis is dominated by the decay of N_1 so that M_1∼ T_ lepto. Eq.(<ref>) can be shown to be accurate, matching the numerical results obtained from Eq.(<ref>),(<ref>) and the modified leptogenesis equations in the presence of PBHs <cit.>. The PBHs can also populate the decaying RHNs to produce lepton asymmetry, but if Eq.(<ref>) 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.<ref>. 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.<ref>) in a way that Eq.(<ref>) prevails. A detailed discussion related to this can be found in Ref.<cit.>. Note now the features of the model. First, The DM and the PBH masses are anti-correlated (cf. Eq.(<ref>)). 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.(<ref>)), β increases, consequently the duration of PBH domination. Therefore, as the RHN mass scale (including the DM) approaches Λ_ GUT,  remarkably, besides maintaining the  m_i≃Λ_ EW^2/Λ_ GUT relation to a great extent, the GWs produced from cosmic strings automatically evades LIGO bound owing to large β or the long duration of PBH-domination. Assuming the hierarchy in the RHN masses[This assumption is not necessary, but more convenient from the model building perspective. As such, one can choose to impose a simple condition  M_ DM>T_ RH>M_Ri(=1,2) on the parameter space, to avoid DM production from scattering while maintaining the condition for thermal leptogenesis T_ RH>M_Ri(=1,2). ] as M_ DM>M_Ri(=1,2) we can re-write Eq.(<ref>) as
        μ=π M_DM^2 f_ DM^-2h(λ,g^'),
        where 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 𝒪(1) couplings in Eq.(<ref>), and using Eq.(<ref>) we obtain an expression for Gμ as 
        
        Gμ≃ 4.7× 10^-8(M_ BH/10^6 g)^-5. 
        With this relation, we now redo the PTA fit for metastable strings performed in the previous section. The result is shown in Fig.<ref>. Note that because of this extra condition in Eq.(<ref>), this U(1)_B-L^ PBH model fits the PTA data for an extremely constrained range of PBH mass and produces  GW signals 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. 
        
        
        
        
        § 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 GW signals 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 non-observation 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∈[10^6g,10^8 g] 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. 
        
        
        § 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).
         
        
        § SUPERSTRINGS AND METASTABLE STRINGS FIT TO THE NANOGRAV DATA
        
         In Fig.<ref> we present the superstring and metastable fit to the PTA data using <cit.>. The 2D posteriors have been used in the main text.
