First-order Phase Transition interpretation of PTA signal produces solar-mass Black Holes

We perform a Bayesian analysis of NANOGrav 15yr and IPTA DR2 pulsar timing residuals and show that the recently detected stochastic gravitational-wave background (SGWB) is compatible with a SGWB produced by bubble dynamics during a cosmological first-order phase transition. The timing data suggests that the phase transition would occur around QCD confinement temperature and would have a slow rate of completion. This scenario can naturally lead to the abundant production of primordial black holes (PBHs) with solar masses. These PBHs can potentially be detected by current and advanced gravitational wave detectors LIGO-Virgo-Kagra, Einstein Telescope, Cosmic Explorer, by astrometry with GAIA and by 21-cm survey.


I. INTRODUCTION
By measuring cross-correlations in the arrival times of pulses emitted by rotating neutron stars, Pulsar Timing Arrays (PTAs) have been established as a mean to detect nano-Hertz (nHz) frequency Gravitational Waves (GW).In 2020, a common low-frequency noise has been identified in the datasets of NANOGrav [1], EPTA [2], and PPTA [3], and confirmed in 2022 by IPTA [4] (IPTA2) which combines data from the former.To distinguish a GW origin from systematic effects requires timing delay correlations to have a quadrupolar dependence on the angular separation between pulsars [5].In June 2023, following the analysis of their most recent data, the collaborative efforts of NANOGrav, EPTA and PPTA (NG15, EPTA2 and PPTA3) have identified compelling statistical evidence for such interpulsar correlations [6][7][8], with Bayes factors of 600, 60, and 11, respectively.The primary expected source of GWs at low frequencies is believed to be from supermassive black holes binaries (SMBH) [9][10][11].The stochastic GW background (SGWB) inferred from PTA data corresponds to the upper limit of the astrophysical predicted interval, see Fig. 1.This could suggest that SMBH binaries are slightly more massive and more numerous than initially anticipated [12][13][14][15].Alternatively, the PTA SGWB might originate from new physics taking place in the early universe [15][16][17][18].The last hypothesis however comes with its own set of challenges.For instance, ascribing the SGWB to inflation necessitates unnaturally large values for the spectral tilt n t ≃ 1.8 and a low reheating temperature T reh ≲ 10 GeV [19].GW induced by a Gaussian spectrum of curvature perturbation would results in excessive PBHs production [20][21][22].A SGWB resulting from PBH mergers would not align with structure formation [23,24].A cosmic strings network, when arising from a global symmetry is excluded by Big-Bang Nucleosynthesis (BBN) [25][26][27][28], while when arising from a local symmetry is not favoured by the Bayesian analysis [16,29].To evade BBN bound, a first-order phase transition (1stOPT) sourcing PTA signal would necessitate the latent heat to be released dominantly to the Standard Model (SM), e.g.[15][16][17][18][30][31][32][33].Interestingly however, the 1stOPT interpretation of PTA SGWB requires a reheating temperature around the scale of QCD confinement 100 MeV, with a rather low completion rate β/H ≲ 12 and a large latent heat fraction α ≳ 0.5 [16].This overlaps with the region where 1stOPT have been recently found to produce PBHs in observable amount [34].The PBH prior has been omitted in all previous analysis of the 1stOPT interpretation of PTA data [15-18, 30-33, 35-46].Similarly, domain wall networks annihilating into SM degrees of freedom are credible early-universe interpretation of the PTA signal associated with the production of of multi-solar-mass PBHs [47,48].
In this letter, we perform a Bayesian search for SGWB from 1stOPT in NANOGrav 15-year (NG15) and IPTA DR2 (IPTA2) timing residuals, including both BBN-N eff -bound and PBH-overproduction constraints as priors in the analysis.To simplify the numerical strategy, we focus on the region α ≫ 1 of strong supercooling, where PBH production is the most efficient.In this region, the dependency of the GW signal on both the wall velocity (v w = 1) and the latent heat fraction α disappears. 1 We argue that the SGWB from 1stOPT is given by the bulk flow model independently of whether the latent heat is still stored in bubble walls at percolation or has been released to the plasma before.We find that PBH formation does not exclude the 1stOPT interpretation of PTA signal.Instead, a SGWB from supercooled PT is favoured with respect to the SMBH binary hypothesis by a Bayes factor of 15 in NG15 data set.We point for the first time, the existence of a multi-messenger window: the NG15 posterior contains a region producing [0.1 − 10] solar-mass PBHs, see Fig. 3.The merging of such PBHs could source GWs with kHz frequencies in the range of LIGO-Virgo [49][50][51][52][53], and ET/CE [54,55].Additionally, their presence could be detected from lensing in GAIA [56][57][58] or from heating in 21-cm survey [59][60][61].
We also consider the negative hypothesis in which the SGWB observed in PTA would not result from a supercooled PT and derive lower limits on the rate of completion β/H ≳ [10 − 20], implying that the universe could not have boiled longer than [5% − 10%] of a Hubble time during the QCD phase transition.

II. GRAVITATIONAL WAVES FROM FIRST-ORDER PT
PT parameters -The strength of a 1stOPT is characterized by the ratio of its latent heat ∆V , defined as the vacuum energy difference between the two minima of the potential driving the transition, to the radiation energy density ρ rad (T n ) at the nucleation temperature where we have neglected a ratio of number of relativistic degrees of freedom.A 1stOPT is said supercooled when α ≳ 1, in which case the universe enters a stage of vacuumdomination at temperature T eq which ends at T n when bubble growth converts the latent heat into radiation energy density.The rate at which nucleation takes place is controlled by the time derivative of the tunneling rate per unit of volume After the phase transition completes, the universe is reheated back to the temperature T eq up to changes in number of degrees of freedom which we again neglect.Energy budget -The dynamics of weak phase transition α < 1 is rather well understood [63,64].The non-relativistic motion of bubble walls, γ w ≃ 1, converts the latent heat into thermal and kinetic energy of the plasma, which propagate under the form of long-lasting sound waves [65], and ultimately turn into turbulence [66][67][68][69][70]. GWs sourced by sound waves have been intensively simulated on the lattice in the recent years [71][72][73][74][75], and analytical modelling have been proposed [76,77].The dynamics of supercooled phase transition α > 1 is more complex due to the large Lorentz factor γ w ≫ 1 of bubble walls [78,79].In the relativistic limit, the acceleration of bubble walls with tension σ is set by the pressure balance [34] The friction pressure P fric is dominantly induced by transition radiation [78], which resummed at leading-logs, reads [79] where g D is a gauge coupling and v ϕ is the vev of the scalar field driving the phase transition.As bubble walls accelerate, the retarding pressure P fric grows linearly with γ w .Scalar field gradient -It is necessary to distinguish two scenarios according to whether the retarding pressure stops the walls from accelerating before collision P fric = ∆V or not P fric ≪ ∆V [79,81].In the later case, bubble walls runaway γ w ↗, and the latent heat is dominantly kept in terms of bubble wall kinetic energy which is the main source of GWs.This occurs for very large supercooling GWs from scalar field gradient were first computed in the socalled "envelop" approximation where walls are infinitely thin and collided parts are neglected [82][83][84][85][86]. Later, collided parts were added to the computation in the so-called "bulk flow" model at the analytical [87] and numerical level [88][89][90][91].It was found that the long-lasting propagation of the infinitely thin shells produces an IR enhancement of the GW spectrum as Ω PT ∝ f 1 instead of Ω PT ∝ f 3 .For relativistic wall velocities, the bulk flow model predicts [88] with the spectral shape S PT (f ) peaked on f ϕ and the redshift factor between percolation " * " and today "0" a * /a 0 = 1.65×10 −2 mHz T eq 100 GeV We added the correction factor with c * = O(1) to impose an f 3 scaling for emitted frequencies smaller than the Hubble factor H * /(2π) as required by causality [92][93][94][95].We fix c * = 1 and leave the determination of c * for future studies.Plasma dynamics -If Eq. ( 5) is not satisfied, bubble walls reach a constant Lorentz factor γw = 0, and the latent heat of the phase transition is dominantly transferred to the plasma, which is the main source of GWs.Friction-dominated bubble wall motion is expected to generate extremely thin and relativistic fluid configurations, which become long-lasting shock waves after bubble collisions [96].The large hierarchy between the bubble radius and the thickness of the shock front is a major challenge to numerical treatment.However, from a gravitational viewpoint an extremely peaked momentum distribution carried by the plasma should be indistinguishable from an extremely peaked momentum distribution carried by the scalar field.Hence we expect the GW signal in both situation to be similar.A second difficulty in modelling plasma dynamics is the possibility for bubble walls to be followed ↑ f < 1 yr The violin diagrams depict the posterior probability distribution of the SGWB energy density in each frequency bins of NG15 and IPTA2 data sets.We overlay with solid lines the SWGB from 1stOPT, obtained using Eq. ( 6), using mean posterior value for the PT parameters.The dotted lines illustrate the SGWB originating from SMBH binaries, employing the mean posterior value for the amplitude and fixing the power-law index to β = 2/3.The gray band represents the 90% confidence interval for the projected SGWB based on a Monte Carlo simulation of a binary population of SMBHs [62].
β/H 68% and 95% confidence contours NG15 (14 bins) IPTA2 (13 bins) Left: Colored regions are posterior distributions in term of the reheating temperature T reh and rate of completion β/H of a strong 1stOPT (α ≫ 1).They are obtained after performing a Bayesian analysis of PTA dataset.We overlay the CMB, LIGO/Virgo and microlensing (EROS) constraints on PBHs produced during such 1stOPT.Right: Lower limit on the rate of completion β/H in the negative hypothesis in which the PTA SGWB would not arise from a strong PT (α ≫ 1).We cast 68% and 95% lower limit using Bayesian inference as explained in App.A (orange and blue), or using the Power Law Integrated Curve in [80] (gray) assuming a signal-to-noise ratio (SNR) threshold of 5 and 10.
by relativistic shells of free-streaming particles [79,[97][98][99], breaking down the fluid description.A recent study in the moderately relativistic regime γ w ≲ 10 [100] suggests that the GW spectrum again resembles the one predicted in bulk flow model.For the two aforementioned reasons, in the present work we assume the GW signal to be given by the bulk flow model in Eq. ( 6) in the whole strongly supercooled regime T n ≪ T eq , independently of whether Eq. ( 5) is satisfied or not. 2

Numerical strategy -
We searched for GW from 1stOPT in two open-access datasets, NG15 [6] and IPTA2 [4].The released data are presented in terms of the timing-residual cross-power spectral density S ab (f ) ≡ Γ ab h 2 c (f )/(12π 2 )f −3 , where h c (f ) ≃ 1.26 • 10 −18 (Hz/f ) h 2 Ω GW (f ) signifies the characteristic strain spectrum [101] and Γ ab denotes the Overlap Reduction Function (ORF) between pulsars 'a' and 'b' within a given PTA [102].We used the software packages enterprise [103] and enterprise extensions [104] to compute the likelihood of observing given timing residuals assuming the presence of the SGWB from 1stOPT given in Eq. ( 6).We used PTMCMC [105] to generate the posterior distribution.For IPTA2, we marginalized over white, red and dispersion measure noises as prescribed in [4,21,106].For NG15, we instead used the handy wrapper PTArcade [107] with "enterprise" mode in which marginalization over noise parameters is automatized.We used GetDist [108] tool to plot the results.To circumvent pulsar-intrinsic excess noise at high frequencies, the SGWB search was confined to the lowest 14 and 13 frequency bins of the NG15 and IPTA2 datasets, respectively.We included the BBN constraints assuming that the 1stOPT sector reheates dominantly into Standard Model degrees of freedom and, when specified, the one from PBH overproduction discussed in Sec.IV, to infer the prior distribution of 1stOPT parameters.Detailed information regarding data analysis and prior choices can be found in App.A. Table I.Mean parameter values with 68% confidence interval of the probability distribution for distinct prior information (rows) and observed data (columns).The GW spectrum from 1stOPT is taken from Bulk Flow model Eq. ( 6) in the supercooled limit α ≫ 1.
Supercooled PT -We conducted searches for GW from strong 1stOPT (α ≫ 1) in isolation, GW from SMBH binaries individually, as well as a combined analysis of 1stOPT and SMBH binaries.In Fig. 1, we show the GW spectra with parameters set to their mean posterior values given in Tab.I.The 68% and 95% confidence contours are depicted in Fig. 2left.The posterior for the combined analysis of 1stOPT and SMBH is reported to Figs. 5 and 6 in the appendix.We assumed a flat prior on the strain amplitude of the SGWB from SMBH binaries, as well as the spectral slope of 13/3 associated with GW-driven inspirals.To quantify the evidence provided by the observed PTA data, denoted as D, in favor of one model, say X, versus another, say Y , we employ the Bayesian factor BF Y,X ≡ P(D|Y ) / P(D/X), (10) which we compute using the product-space sampling method [102] implemented in enterprise extensions [104].
Here, P(D/X) is the likelihood probability of observing data D given the model X.The outcomes of the Bayesian model comparison presented in Tab.II, according to Jeffrey's scale [109,110], suggests that NG15 data 'substantially' favours the presence of a GW signal from 1stOPT aside to the one from SMBHB.Instead, IPTA2 data remains inconclusive.Exclusion bounds -Under the assumption that the PTA signal does not arise from 1stOPT, we have derived upper limits on the GW signal emanating from 1stOPT.As depicted in Fig. 2-right, these limits correspond to lower bounds on the rate of completion, going up to β/H ≲ 20.As discussed in App.A, these lower limits are conservative as the GW spectrum from SMBH was not included in the analysis.

IV. PRIMORDIAL BLACK HOLES
Supercooled late-blooming mechanism -In [34], it was demonstrated that PBHs could be produced in observable amount during supercooled PT through a process termed "late-blooming".During 1stOPT, the nucleation sites of bubbles are randomly dispersed across the entire volume of the false vacuum.As the universe gets close to the point of percolation, there remains a non-zero probability of identifying Hubble-sized regions where nucleation has not yet initiated.Throughout the supercooled PT, these delayed regions maintain a constant vacuum energy, while the energy density in their vicinity redshifts like radiation.Upon completion of percolation, these "late-bloomers" evolve into overdense regions.If these regions are Hubble-sized and exceed a certain density threshold δρ/ρ ≳ 0.45, they collapse into  PBHs.We direct the reader to [34] for the precise analytical formula to estimate the abundance and mass of those PBHs. 3he mass distribution of those PBHs, left for future studies in [34], is assumed to resemble a delta function in the present work.We included the PBH overproduction constraints as a prior in the Bayesiasn analysis.The Bayes factors shown in Tab.II is unaffected for IPTA2 and only decreases from 24 to 15 for NG15.We have plotted the contour lines representing the PBH fraction of dark matter f PBH in Fig. 2 and the PBH mass in Fig. 3.In addition, we overlay cosmological and astrophysical constraints on this population of PBHs.Excluded regions and detection prospects -With solid lines, we show current constraints.In yellow, we have the exclusion regions arising from distortion of the Cosmic Microwave Background (CMB) caused by X-rays from accretion which modify the ionization history between recombination and reionizaton [116][117][118].In purple, we show the constraints using the search for photometric magnification (strong lensing) of stars in the Magellanic clouds conducted on Eros data [119].The solid cyan-colored region represents constraints derived from the data collected by LIGO/Virgo in-terferometers [49][50][51][52][53].With dashed lines, we show future prospects.In green, we have the reach of 21 cm surveys due to heating and ionization of the intergalactic medium via X-rays produced during accretion [59][60][61].In red, we have the forecast from the search for transient astrometric deviation (weak lensing) of single or multiple stars in GAIA time-series data [56][57][58].Finally, in dashed cyan we show the prospect for detecting GW from PBH binaries with Einstein telescope and Cosmic Explorer [54,55].

V. CONCLUSION
We conducted a Bayesian analysis of the NANOGrav 15-yr (NG15) and IPTA DR2 (IPTA2) timing residuals.Our findings indicate that NG15 indicate a substantial preference for the presence of a strong first-order phase transitions (1stOPT) in isolation or combined with SGWB from SMBH binaries, while IPTA2 remains inconclusive on which scenario is preferred.The phase transition is characterized by a remarkably low completion rate, e.g.β/H ≃ 12.6 and 10.7 for NG15 with and without astrophysical signal from SMBH binaries.From a theoretical perspective, such a value is typical of supercooled phase transitions, characterized by a strong firstorder phase transition with a parameter α significantly larger than 1, e.g.[63,64,120], which motivates the choice of prior α ≫ 1 done in this work.These cosmological scenarios have been demonstrated to produce primordial black holes (PBHs) in considerable quantities when β/H ≲ 7 [34].The Bayes factor of the strong 1stOPT interpretation with respect to SMBH binary one is only reduced from 24 to 15 in NG15 after including the PBH prior, while it is not affected in IPTA2.
However, we showed that the 1stOPT interpretation of the PTA signal might be associated with the presence of solarmass PBHs in our universe today.We further assessed the potential for detecting these PBHs using different observational techniques, including 21 cm cosmological hydrogen line observations, astrometry with the GAIA mission and next-generation kilohertz frequency GW interferometers such as the Einstein Telescope (ET) and Cosmic Explorer (CE).We conclude that 1stOPTs can be ranged alongside domain wall networks [47,48] and scalar induced GW [20][21][22] in the category of the early-universe interpretations of PTA signal capable of producing multi-solar-mass PBHs in quantities that are potentially observable.
In the event that an astrophysical explanation becomes definitive, we established 68% and 95% exclusion constraints on the parameter space of 1stOPT, up until β/H ≳ 20.Under these conditions, it would effectively preclude any possibility of detecting PBHs from supercooled PTs within the mass range We must emphasize that our current comprehension of the GW spectrum resulting from supercooled phase transitions is still in its early stages.The assumptions are founded on the bulk flow model, in which GWs are sourced by the expansion of an infinitely thin distribution of the stress-energy momentum tensor.Future investigations are necessitated to probe potential modifications of the GW spectrum that could be induced by non-linear effects, such as those arising from relativistic shock waves, or deviations from a fluid description.
Acknowledgements.-The author is grateful to Iason Baldes, Ryusuke Jinno, Marius Kongsore, Fabrizio Rompineve, Miguel Vanvlasselaer and Tomer Volansky for fruitful discussions and to the Azrieli Foundation for the award of an Azrieli Fellowship.This work was conducted using the high performance computing cluster resources of Tel Aviv University.We compare the posterior resulting from the Bayesian analysis realized with our own wrapper of enterprise extensions [104] as described in App.A (orange) with the one using the public software PTArcade [107] (green) .Right: Numerous analysis in the literature were performed with NANOGrav 12.5-yr (NG12) with 5 frequency bins (purple), which is very different from the recent NG15 14-bin posterior (blue).We also compare the posterior obtained with the fast method of [126], called "ceffyl", which relies on a fit of the violin plot (red).All posteriors in the right panel were generated with PTArcade [107].

A. Data analysis
The purpose of this Appendix is to delineate the Bayesian search methodology employed in our study.We started rely on the NG15 dataset [127] and on Version B of the IPTA2 dataset [128].To ascertain noise parameters of IPTA2, we closely follow the approach adopted by IPTA collaboration [4], see also [21,106].We then checked that we obtained consistent result with the software PTArcade [107] in which noise marginalization has been automatised, see Fig. 4-left.Instead the Bayesian analysis of NG15 was done solely with the "enterprise" mode of PTArcade [107].We perform the search for SGWB in the first 13 and 14 frequency bins of IPTA2 and NG15, respectively.IPTA2 analysis.-We now describe the Bayesian analysis of IPTA2 which we performed ourselves without the use of PTArcade [107].We adapted the software packages enterprise [103] and enterprise extensions [104] to incorporate GW spectra from 1stOPT in terms of the power spectrum in timing residual, and used them to compute the likelihood function, symbolized as P(D|θ).This function encapsulates the probability of observing the data D given a specific set of model parameters θ.The posterior distribution, P(θ|D), which illustrates the probability distribution of model parameters θ given the observed data D, is linked to the likelihood function via Bayes's theorem Within this equation, P(θ) is the prior distribution, representing preliminary knowledge of the parameters prior to data observation, while P(D) is the marginal likelihood or evidence, functioning as a normalization constant to ensure that the posterior distribution integrates to 1.The parallel-tempering Markov Chain Monte-Carlo sampler PTMCMC [105] was employed to reconstruct the posterior distribution P(θ|D) using an enhanced version of the Metropolis-Hastings algorithm [102].The GetDist tool [108] was subsequently used to plot the posterior distributions and upper limits.The pulsar noise parameters employed in the likelihood function can be classified into three distinct categories: white noise, red noise, and dispersion measures (DM).The white noise parameters are grouped into three sets for each backend/receiver associated with a given pulsar: ).The values of the white noise parameters are fixed to the mean posterior values obtained by performing single pulsar analysis devoid of GW signals.We only kept pulsars with more than 3 years of observation time which corresponds to 53 pulsars.Instead the Bayesian analysis of NG15 data performed via PTArcade contains 68 pulsars with more than 3 years of observation.We employ the Jet Propulsion Laboratory Development Ephemeris DE438 and the Terrestrial Time reference timescale of the International Bureau of Weights and Measures BIPM18.Next, for the multi-pulsar analysis incorporating the GW signals, we account for two power-law red noise parameters per pulsar, specifically the amplitude at the reference frequency of yr −1 denoted as A red , and the spectral index denoted as γ red .Additionally, we incorporate power-law errors associated with dispersion measures (DM).We note that the treatment of DM noise as a Gaussian process is specific of IPTADR2 dataset.Instead, in the analysis of NG15 data performed via PTArcade, but also in the analysis of NANOGrav 12.5-year (NG12) done in [1], pulse dispersion is modelled by a set of "per-epoch" parameters describing the DM offset from a nominal fixed value [129,130].These can add dozens of additional parameters per pulsar [131].In the individual pulsar analysis of PSR J1713+0747 (in IPTA2 but also in NG15), we extend our consideration to encompass a DM exponential dip parameter, following the methodology described in [4].The priors for the noise parameters are reported in Tab.III, along with the priors for the parameters for the GW spectra from 1stOPT and SMBH binaries.To economize on computational time, we adopt the methodology of previous studies [4,132] and in our search for a GW background we utilize only auto-correlation terms I = J in the Overlap Reduction Function (ORF) Γ IJ , rather than the complete Hellings-Downs ORF with I ̸ = J.We acquire 10 6 samples per analysis presented in this study and discard 25% of each chain as burn-in.We could replicate the posteriors of [132] and [4] for a power-law model with excellent concurrence.The violin features shown in Figs. 1 and 5 are obtained with the free-spectrum approach described in [133].We do not repeat this analysis and instead take the data directly from NG15 and IPTA2.
Our study encompasses two types of analyses.The first, a detection analysis, identifies the region of parameter space in which GWs from 1stOPT can account for the common-spectrum process in the datasets.Here, we use a uniform prior on the logarithm of each parameter and adopt a prior on β/H due to the BBN bound and -when mentioned -PBH overproduction.The second, an lower-limit analysis, seeks to constrain the rate of completion of the phase transition β/H.There, we use a uniform prior on H/β instead of log 10 (β/H) as described in [102,134].We made the conservative choice to not include the GW spectrum from SMBH in the lower-limit analysis, see the related discussion in [21].All prior choices are given in Tab.III.
BBN prior.-As a sub-component of the total energy density of the universe, the latent heat ∆V can impact the expansion rate of the universe which is strongly constrained by BBN and CMB.Its effect can be encoded in the effective number of extra neutrino relics where ρ γ is the photon number density.The total number of effective degrees is constrained by CMB measurements [135] to N eff = 2.99 +0.34 −0.33 and by BBN predictions [136,137] to N eff = 2.90 +0.22 −0.22 whereas the SM prediction [138,139] is N eff ≃ 3.045.The latent heat parameter of a generic 1stOPT reads where T is the photon temperature and g * (T ) contains eventual dark degrees of freedom.The maximal contribution to N eff occurs at reheating after percolation The BBN bound ∆N eff ≲ 0.3 [140,141] applies after neutrino decouples below the temperature T dec where g * (T < T dec ) ≡ 2 + (7/8) • 6 • (4/11) 4/3 ≃ 3.36.We obtain Two scenarios must be distinguished.The first one is when reheating after percolation occurs in a dark sector, in which case Eq.A5 is the BBN constraints.The second one is when reheating after the 1stOPT occurs into the Standard Model, in which case Eq.A5 applies only if the reheating temperature is below the neutrino decoupling temperature T reh ≲ 1 MeV.The last case is the scenario we consider in this work.Note that stronger BBN constraints have been considered in the literature, e.g.T reh ≲ 3 MeV in [32], or T reh ≲ 2 MeV and 4 MeV for cases of electromagnetic and hadronic decays respectively [142,143].PBH prior.-The condition of not producing PBH with an energy density larger than the one of observed dark matter, f PBH < 1, implies a lower bound on the rate of completion of a 1stOPT [34] β/H ≳ 5. 54  where we have introduced an analytical function fitted on numerical results of [34].When specified, we include the constraint in Eq. (A6) as prior information on β/H and T eq .Due to the exponential dependence of the PBH abundance on β/H, the precise PBH constraints due to astrophysical and cosmological constraints, as shown in e.g.Fig. 3, make little difference with respect to simple criterion f PBH < 1.

B. Combined GW from 1stOPT and SMBH binaries
The characteristic strain spectrum of a population of circular GW-driven SMBH binaries is a red-tilted power-law [144]    where A SMBH is the strain amplitude at 1 yr −1 ≃ 3.2×10 −8 .In terms of the fractional energy density, it corresponds to the blue-tilted power-law

Figure 3 .
Figure 3.The ellipses are the posterior distributions obtained after a Bayesian search of SGWB sourced by a supercooled 1stOPT in NG15 and IPTA2 data sets.We overlay the region producing PBHs detectable by different observatories, see Sec.IV for details.

Figure 4 .
Figure 4. Left: We compare the posterior resulting from the Bayesian analysis realized with our own wrapper of enterprise extensions[104] as described in App.A (orange) with the one using the public software PTArcade [107] (green) .Right: Numerous analysis in the literature were performed with NANOGrav 12.5-yr (NG12) with 5 frequency bins (purple), which is very different from the recent NG15 14-bin posterior (blue).We also compare the posterior obtained with the fast method of[126], called "ceffyl", which relies on a fit of the violin plot (red).All posteriors in the right panel were generated with PTArcade[107].

1 − 2 / 3 , 1 ↑Figure 5 .
Figure 5.We show the combined GW signal from 1stOPT and SMBH binaries with mean posterior values for the rate of completion β/H and reheating temperature T reh (blue and orange).The gray band shows the SGWB from the incoherent superposition of a population of Monte-Carlo-simulated SMBH binaries[62].The band brackets 90% of the simulated population.

Figure 6 .
Figure 6.We performed a Bayesian analysis of a superposition of SGWB from 1stOPT and SMBH binary mergers.

Table II
. Bayesian factors BFY,X with values significantly exceeding 1 indicate support for interpretation Y with respect to X. Conversely, values approaching 1 suggests no discernible preference between X and Y .We can see that the 1stOPT interpretation is favored with respect to SMBH binaries in NG15 data and that the PBH prior only slightly worsens the fit.

Table III .
Prior assumptions on the parameters used in the Bayesian analysis of this work.