Pulsar timing array PTA data releases show strong evidence for a stochastic gravitationalwave background in the nanohertz band. When the signal is interpreted by a scenario of scalarinduced gravitational waves SIGWs we encounter overproduction of primordial black holes PBHs. We wonder if varying the equation of state EoS of the early Universe can resolve this issue and thereby lead to a consistent interpretation of the PTA data. Analyzing a data combination of PTA bigbang nucleosynthesis and cosmic microwave background we find that an epoch with EoS inlineformula texmathwsimcalO102texmath inlineformula between the end of inflation and the onset of radiation domination can significantly suppress the production of PBHs leading to alleviation of the PBHoverproduction issue. With the inferred interval inlineformula texmathw0.440.400.52texmath inlineformula at 95 confidence level our scenario can interpret the PTA data just as well as the conventional scenario of SIGWs produced during the radiation domination. QingHua Zhu is Supported by the National Natural Science Foundation of China 12305073. ZhiChao Zhao is supported by the National Key Research and Development Program of China 2021YFC2203001 and the National Natural Science Foundation of China 12005016. Sai Wang is partially supported by the National Natural Science Foundation of China 12175243 the National Key RD Program of China 2023YFC2206403 the Science Research Grants from the China Manned Space ProjectCMSCSST2021B01 and the Key Research Program of the Chinese Academy of Sciences XDPB15. Xin Zhang is supported by the National SKA Program of China 2022SKA0110200 2022SKA0110203 and the National Natural Science Foundation of China 11975072 11835009 11805031. This work is supported by Highperformance Computing Platform of China Agricultural University
['physics of the early universe', 'gravitational waves', 'primordial black hole', 'pulsar timing arrays', 'Astrophysics - Cosmology and Nongalactic Astrophysics', 'Astrophysics - High Energy Astrophysical Phenomena', 'General Relativity and Quantum Cosmology']
Unraveling the early universes equation of state and primordial black hole production with PTA BBN and CMB observations
{"Unraveling the early universe's equation of state and primordial black hole production with PTA, BBN, and CMB observations": "Qing-Hua Zhu  \n1 Zhi-Chao Zhao 2 Wang Sai 3 Xin Zhang 4 , 5 , 6  \n- 1 Department of Physics and Chongqing Key Laboratory for Strongly Coupled Physics, Chongqing University, Chongqing 401331, People's Republic of China\n- 2 Department of Applied Physics, College of Science, China Agricultural University, Qinghua East Road, Beijing 100083, People's Republic of China\n- 3 Theoretical Physics Division, Institute of High Energy Physics, Chinese Academy of Sciences, Beijing 100049, People's Republic of China\n- 4 Key Laboratory of Cosmology and Astrophysics (Liaoning) & College of Sciences, Northeastern University, Shenyang 110819, People's Republic of China\n- 5 National Frontiers Science Center for Industrial Intelligence and Systems Optimization, Northeastern University, Shenyang 110819, People's Republic of China\n- 6 Key Laboratory of Data Analytics and Optimization for Smart Industry (Ministry of Education), Northeastern University, Shenyang 110819, People's Republic of China  \nAbstract: Pulsar timing array (PTA) data releases showed strong evidence for a stochastic gravitational-wave background in the nanohertz band. When the signal is interpreted by a scenario of scalar-induced gravitational waves (SIGWs), we encounter overproduction of primordial black holes (PBHs). We wonder if varying the equation of state (EoS) of the early Universe can resolve this issue and thereby lead to a consistent interpretation of the PTA data. Analyzing a data combination of PTA, big-bang nucleosynthesis, and cosmic microwave background, we find that an epoch with EoS w ∼O (10 -2 ) between the end of inflation and the onset of radiation domination can significantly suppress the production of PBHs, leading to alleviation of the PBH-overproduction issue. With the inferred interval w =0 . 44 +0 . 52 -0 . 40 at 95% confidence level, our scenario can interpret the PTA data just as well as the conventional scenario of SIGWs produced during the radiation domination.  \nKey words: Physics of the early universe, gravitational waves, primordial black hole, pulsar timing arrays", '1 Introduction': 'Recently, multiple data releases of pulsar timing arrays (PTAs) showed significant evidence for a stochastic gravitational-wave background in the nanohertz band [14]. The signal has been interpreted to arise from scalarinduced gravitational waves (SIGWs) by the authors of Refs. [5-22]. Though this scenario fits the PTA data better than other scenarios, it results in an issue of overproduction of primordial black holes (PBHs) [23, 24]. In these works, the authors have assumed that SIGWs were produced during an epoch of radiation domination with the equation of state (EoS) w =1 / 3. However, the EoS of the early Universe might deviate from the radiation domination [25-39]. We wonder if the overproduction issue can be resolved via considering an arbitrary EoS of the early Universe during an epoch between the end of inflation and the onset of radiation domination.  \nIt is well known that SIGWs could be nonlinearly produced by linear cosmological curvature perturbations,  \nwhen the latter reentered the Hubble horizon after the end of inflation [40-45]. The energy-density spectrum of SIGWs would be significantly affected by the EoS of the Universe at the production time of SIGWs [28, 39, 4552]. The corresponding imprints may be left on the PTA signal, since the Universe is transparent to the gravitational waves [53]. Therefore, we anticipate that SIGWs is an effective probe of the early-Universe physics, particularly for the EoS parameter.  \nIt is also known that PBHs were produced due to gravitational collapse of these perturbations if the latter excess a critical overdensity when they reentered into the Hubble horizon [54]. The mass function of PBHs, and thereby the abundance, would also be significantly impacted by the EoS of the Universe at the formation time of PBHs [51, 55-57]. Therefore, the issue of PBH overproduction might be resolved when an appropriate EoS of the early Universe is taken into account.  \nIn this work, we study a scenario of SIGWs produced during the epoch with an arbitrary EoS between the end  \nof inflation and the onset of radiation domination. Differing from studies assuming that all of SIGWs were produced during the w -domination before the onset of radiation domination [49, 52], we consider that SIGWs were produced during both epochs. Our scenario is inspired by an enhancement mechanism for SIGWs caused by a sudden transition from the early-matter domination to the radiation domination [46, 47]. We expect that a similar enhancement mechanism can also exist in our scenario. By fitting the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) 15year data release [3] to our scenario, we study whether varying EoS of the early Universe can improve the ability of SIGWs to interpret the PTA signal. Meanwhile, we propose that the enhancement of SIGW spectrum can be significantly constrained by further taking into account the big-bang nucleosynthesis (BBN) [58] and cosmic microwave background (CMB) [59] data. It is worth noting that our results would be different from those of Refs. [23, 24], where the authors considered SIGWs produced during the radiation domination only and analyzed the PTA data only. Moreover, we study the production process of PBHs during the w -dominated epoch. We will show that the PBH abundance is significantly changed, probably leading to a resolution of the PBHoverproduction issue.  \nThe remaining context of this paper is arranged as follows. In Section 2 , we develop the scenario of SIGWs by taking into account the arbitrary EoS of the early Universe. In Section 3 , we demonstrate the method of data analysis and show the inferred parameter region. In Section 4 , we study the formation of PBHs in our scenario and the issue of PBH overproduction. In Section 5 , we show conclusions and discussion.', '2 SIGW energy-density spectrum': "It has been known that a sudden transition from the early-matter domination to the radiation domination can result in a significant enhancement of the energy-density spectrum of SIGWs [46, 47]. In this work, we replace the early-matter domination with an epoch dominated by 'matter' with an arbitrary EoS parameter w ∈ [0 , 1], following the research method of Refs. [48, 49, 52]. However, we still consider the sudden transition, since in general, it is too complicated to study a finite duration of the transition, which is left to future works. We will study the imprints of such a scenario on the energy-density spectrum of SIGWs in the following. In fact, there are several mechanisms to generate the sudden transition from the early-matter domination to the radiation domination [29, 30, 39, 47, 57, 60-65].  \ni.e., [47]  \na ( η ) a R =    ( η η R ) 2 1+3 w for η < η R 2 (1+3 w ) ( η η R -1 ) +1 for η ≥ η R , (1)  \nwhere η is the conformal time. We introduce a R with η R at the occasion of transition, i.e.,  \na R = a ( η R ) = ( 1+3 w 2 ) H 0 η R √ Ω m, 0 a eq , (2)  \nη R = 8 × 10 9 s ( 2 1+3 w )[ g ∗ ( T R ) 10 . 75 ] -1 6 ( 1MeV T R ) , (3)  \nwhere g ∗ is an effective number of relativistic species as a function of the cosmic temperature T R corresponding to η R , H 0 is the Hubble constant, Ω m,0 is the present-day energy-density fraction of dark matter, and a eq denotes the scalar factor at the epoch of matter-radiation equality. The conformal Hubble parameter H = a ' /a is  \nH = { 2 (1+3 w ) η for η < η R 1 η -1 2 (1 -3 w ) η R for η ≥ η R , (4)  \nwhere the prime denotes a derivative with respect to η .  \nThe equation of motion for SIGWs is shown as [40-  \n45]  \nh '' ij +2 H h ' ij -∆ h ij = -4Λ ab ij S ab , (5)  \nwhere h ij is the strain of SIGWs, Λ ab ij is the transversetraceless operator, and S ab denotes the source term consisting of nonlinear couplings between the linear cosmological perturbations ψ , namely, S ab ∼ ψ 2 . Analytic formulas of S ab can be explicitly obtained via calculating the non-linear cosmological perturbations with xpand [66].  \nThe energy-density fraction spectrum of SIGWs for the sub-horizon modes ( kη ≫ 1) is given by [44, 45]  \nΩ GW ( k,η ) = 1 48 ( k H ) 2 [ P h ( k,η )+ ¯ P h ( k,η ) ] , (6)  \nwhere k is a wavenunber, and P h ( k,η ) and ¯ P h ( k,η ) stand for the two-point correlators of h ij and h ' ij , respectively, defined as  \n⟨ h ij, k h ij, ¯ k ⟩ = 2(2 π ) 3 δ ( k + ¯ k ) 2 π 2 k 3 P h ( k,η ) , (7)  \n⟨ h ' ij, k h ' ij, ¯ k ⟩ = 2(2 π ) 3 δ ( k + ¯ k ) 2 π 2 k 3 ¯ P h ( k,η ) . (8)  \nIn fact, Eq. (8) is redundant, since it can be derived from Eq. (7), but it lets us express Eq. (6) in a simple form. Semi-analytic formulas for h ij, k and h ' ij, k can be obtained explicitly. According to Eq. (5), the stochastic nature of h ij, k is attributed to that of Ψ k , i.e., the initial conditions of ψ . The two-point correlator of Ψ k is defined as  \nWe adopt a sudden transition of the scale factor a ,  \n⟨ Ψ k Ψ ¯ k ⟩ = (2 π ) 3 δ ( k + ¯ k ) [ 3(1+ w ) 5+3 w ] 2 2 π 2 k 3 P ζ ( k ) , (9)  \n10  \nwhere P ζ ( k ) is a power spectrum of primordial curvature perturbations ζ , and we have related ψ to ζ . For simplicity, we consider a monochromatic spectrum, i.e., P ζ = Ak ∗ δ ( k -k ∗ ), with k ∗ being a reference wavenumber. Such a spectrum may be related to the formation of PBHs [67-69], and has been extensively studied in the literature (see reviews in Refs. [70, 71] and references therein). It would be interesting to study a finite-width spectrum [72, 73], which may change parts of our results. However, such a study is left to future works.  \n<!-- image -->  \nFigure 1. Energy-density fraction spectrum of SIGWs produced during the w -domination (dashed) versus the total spectrum (solid). For comparison, we show the conventional scenario of SIGWs (dotted).  \n<!-- image -->  \nIn Fig. 1, we show the energy-density fraction spectra of SIGWs for several different sets of model parameters. We take k = 2 πf and k ∗ = 2 πf ∗ throughout this work. For the small w , i.e., w → 0, the spectrum experiences an enhancement due to the sudden transition. In other words, when the early Universe transits from the almost early-matter domination to the radiation domination, the SIGWs produced during the radiation domination dominates the spectrum. This result is consistent with the enhancement mechanism proposed in Refs. [46, 47], where the authors considered a sudden transition from the early-matter domination to the radiation domination.", '3 Data analysis and results': 'To compare the theoretical results with the observations, we should know the energy-density fraction spectrum of SIGWs in the present Universe, i.e.,  \nΩ SIGW ( f ) ≃ Ω r , 0 × Ω GW ( k,η ) , (10)  \nwhere Ω GW has been shown in Eq. (6), and the physical energy-density fraction of radiations in the present-day Universe is h 2 Ω r , 0 ≃ 4 . 2 × 10 -5 with the dimensionless Hubble constant being h = 0 . 6766, as measured by the Planck satellite [74].  \nBesides the above SIGW component, the PTA signal may be contributed by an astrophysical component arising from the binary supper-massive black holes (SMBHs). The corresponding energy-density spectrum  \nis parameterized as a power-law [24, 75]  \nΩ BHB ( f ) = 2 π 2 f 2 yr 3 H 2 0 A 2 BHB ( f f yr ) 5 -γ BHB , (11)  \nwhere A BHB is the spectral amplitude, γ BHB is the spectral index, and f yr is a pivot frequency corresponding to one year. Therefore, the total spectrum is expressed as a sum of Eq. (10) and Eq. (11).  \nUsing the same statistical tools used by the PTA collaborations, we analyze the NANOGrav 15-year data [3] by using the publicly-available PTArcade [76, 77], which implements the Bayes parameter inferences for PTA data by providing a wrapper of ENTERPRISE [78] and Ceffyl [79]. Here, the total spectrum would be fitted to the PTA data during our data analysis.  \nWe take the BBN [58] and CMB [59] constraints on the effective number of relativistic species into account, which can further pin down the parameter space of SIGWs. In terms of the integrated energy-density fraction, as defined by ∫ ∞ k min d ln k h 2 Ω SIGW ( k ), the upper limits at 95% confidence level are given as 1 . 3 × 10 -6 for BBN [58] and 2 . 9 × 10 -7 for CMB [59]. The lower bound of the integral is denoted as k min =2 πf min , where f min is 1 . 5 × 10 -11 Hz for BBN and 3 × 10 -17 for CMB [80]. The BBN and CMB data have also been taken into account in studies of SIGWs [16] and phase-transition gravitational waves [81]. Here, we take each of them into account by assuming the integrated energy-density spectrum to follow a Gaussian distribution, whose variance is given by one-half of the aforementioned upper limit. Therefore, we adopt a modified version of Ceffyl [79] via adding two likelihoods to account for the BBN and CMB constraints, respectively.  \nThe parameter space under investigation is spanned by log 10 A , log 10 ( f ∗ / Hz), w , log 10 ( f ∗ η R ), log 10 A BHB , and γ BHB . We adopt the uniform priors for the SIGW parameters, i.e., log 10 A ∈U [ -8 , 1], log 10 ( f ∗ / Hz) ∈U [ -9 , -4 . 5], w ∈ U [0 , 1], and log 10 ( f ∗ η R ) ∈ U [0 , 5]. The priors for the last two parameters indicate a prior for η R within [10 9 , 10 14 ] seconds. According to Eq. (3), this prior promises T R to be higher than ∼ 5 MeV, i.e., a temperature at which the w -domination should transition to the radiation domination, as suggested in Refs. [8284]. In addition, the priors for the SMBH parameters, i.e., log 10 A BHB and γ BHB , are the same as those of Ref. [23], which follow a bivariate Gaussian distribution, as wrapped in PTArcade .  \n)  \nη  \n,  \nf  \nπ  \n2  \n(  \nW  \nG  \nΩ  \n2  \n∗  \n∗  \nFigure 2. Posteriors of six independent parameters.  \n<!-- image -->  \nTable 1. Posterior and prior distributions for the parameters in the analysis. Here, U denotes the uniform distribution. We show the results at 95% confidence level.  \n(  \n)  \nFigure 3. Comparison between the best-fit spectrum and the NANOGrav signal.  \n<!-- image -->  \nWe depict the posteriors of six independent parameters in Fig. 2. Compared with Ref. [23], which analyzed the NANOGrav data only, we combine the NANOGrav data with the BBN and CMB data. Hence, the inferred posteriors in our work are changed significantly. In particular, we find that at 95% confidence level, A is bounded from the upper, but not from the lower, i.e., A ≤ 0 . 39, and f ∗ is within 1 . 2 × 10 -8 -1 . 9 × 10 -6 Hz. In contrast, we get looser constraints on w and f ∗ η R . We find  \nthat most of the considered EoS fit the joint dataset just as well as w =1 / 3 within 95% confidence interval. In addition, our constraints on A BHB and γ BHB are compatible with those of Ref. [23]. We summarize the information of priors and posteriors in Tab. 1. Furthermore, we compare the best-fit spectrum with the observed signal in Fig. 3.  \nFigure 4. Comparison of the probability distribution functions of the EoS parameter w between the SIGW-only model (red) and the SIGW+SMBH model (blue).  \n<!-- image -->  \nFollowing the same approach mentioned above, we perform an additional analysis of the SIGW-only model, meaning that we consider the SIGW parameters only and disregard the SMBH parameters. In Fig. 4, we compare the inferred probability distribution functions of w between the SIGW-only model (red curve) and the combined SIGW+SMBH model (blue curve). We find that they are compatible with each other.', '4 PBH production': 'The PBH formation during the radiation domination results in the overproduction of PBHs [23, 24]. In contrast, we expect that PBHs were produced during an epoch instead of the radiation domination. The PBH abundance can be changed to some extent, possibly leading to avoidance of the overproduction issue.  \nWhen the k modes reentered into the Hubble horizon, i.e., k ≃ H = 2 / [(1 + 3 w ) η ], and collapsed to form PBHs, the PBH fraction of dark matter is given by [71]  \nf pbh ( k ) = a 3 Ω CDM ( H H 0 ) 2 γβ ( k )  \n(12) = Ω m, 0 Ω CDM [ a eq a R ]( η η R ) -6 w 1+3 w γβ ( k ) ∣ ∣ ∣ ∣ η = 2 ,  \n(1+3 w ) k  \nwhere H is the Hubble parameter, and we have γ =0 . 356 [23], Ω m, 0 =0 . 3089, and Ω CDM =0 . 2589 [85]. The proba-  \n7  \nbility of PBH production is defined as [86]  \nβ ( k ) = ∫ ∞ δ c 1 √ 2 πσ ( k ) exp ( -δ 2 2 σ ( k ) ) d δ = 1 2 erfc ( δ c √ 2 σ ( k ) ) , (13)  \nwhere the critical overdensity for gravitational collapse is δ c = [3(1+ w ) / (5+3 w )] sin 2 [( π √ w ) / (1+3 w )] [55]. The coarse-grained overdensity is defined as  \nσ 2 ( k ) = ∫ ∞ 0 dln qW 2 ( q,R ) P δ ( q,η ) ∣ ∣ ∣ R = η = 2 (1+3 w ) k , (14)  \n̸  \nwhere W ( k,R ) = exp[ -( kR ) 2 / 2] is a window function, and the power spectrum of overdensity is related to that of primordial curvature perturbations ζ , i.e., P δ = [2(1+ w ) / (5+3 w )] 2 ( k/ H ) 4 T 2 ( kη ) P ζ ( k ) with T ( kη ) being the transfer function of ψ . One should note that the transfer function for the case of w =1 / 3 is different from that for the radiation domination.  \nTo define the PBH abundance, we should derive a relation between the PBH mass M PBH and k via evaluating k/k ∗ = H / H ∗ by making use of M PBH = γM H = γ (4 π/ 3) ρH -3 and Eq. (1). Such a relation is  \nk k ∗ = ( M PBH M PBH , ∗ ) -1+3 w 3(1+ w ) , (15)  \nwhere M PBH , ∗ corresponds to k ∗ . Finally, the PBH abundance is defined as [23, 87]  \nf PBH = ∫ M PBH,R 0 f pbh ( k ( M PBH ))dln M PBH = 3(1+ w ) 1+3 w ∫ ∞ 2 (1+3 w ) η R f pbh ( k )dln k , (16)  \nwhere M PBH,R corresponds to the k modes reentering into the Hubble horizon at η R .  \nFigure 5. The PBH overproduction versus the SIGW scenario with a sudden transition from the epoch of w =1 / 3 to the radiation domination. We denote f PBH = 1 as dashed curves, above which PBHs are over-produced. For comparison, posteriors are shown as shaded regions. In the left panel, w =1 / 3 is shown in the dotted line. In the right panel, f PBH =1 for w =1 / 3 is also shown in the dotted curve.  \n<!-- image -->  \n<!-- image -->  \n̸  \n∗  \nIn Fig. 5, we show to what extent we could resolve the issue of PBH overproduction via comparing f PBH =1, as predicted by our scenarios, with the inferred posteriors. As suggested in Ref. [23], where only the PTA data was analyzed, the issue of PBH overproduction is inevitable in the case of w = 1 / 3. Though being changed significantly, as shown in Sec. 3 , the posteriors associating with the combination of PTA, BBN, and CMB still can not avoid the PBH overproduction. In this sense, we pay attention to scenarios of the varying EoS in the early Universe. It was surprising to find that the scenarios of w ∼O (10 -2 ), which are allowed by our joint dataset, can alleviate the overproduction issue.', '5 Conclusions and discussion': 'In this work, we have studied the influence of the EoS of the early Universe on the production of SIGWs, and further interpreted the recent PTA data releases as well as the overproduction issue of PBHs. Via generalizing the scenario of Refs. [45-47] and following the method of Refs. [48, 49, 52], we computed the energy-density fraction spectrum of SIGWs for a sudden transition of the early Universe from an arbitrary EoS w ∈ [0 , 1] to the radiation domination w =1 / 3. We showed that the SIGW spectrum exhibits a notable enhancement as the parameter w approaches zero, as shown in Fig. 1. To get precise constraints on the model parameters, we analyzed the data combination incorporating the NANOGrav 15-year PTA data with the BBN and CMB data. To our knowledge, it is the first time to perform Bayes data analysis for the considered topic. It is also the first time to combine the PTA data with the BBN and CMB data for the Bayes data analysis. The posteriors of model parameters have been shown in Tab. 1. In addition, studying the PBH production in the w -dominated epoch, we found that the PBH-overproduction issue can be alleviated to some extent, as demonstrated by Fig. 5.  \nIn fact, the production of SIGWs in the early Universe with an arbitrary EoS was studied in Refs. [4850, 52]. In these works, the gravitational waves were produced during the w -domination, without further production during the subsequent radiation domination. However, it was shown that for SIGWs produced during both the early-matter domination and subsequent radiation domination, the SIGW spectrum could be significantly enhanced due to a sudden transition of EoS [46, 47]. We have found the similar enhancement mechanism for our scenario.  \nWe have proposed that the SIGW spectrum is tightly constrained by the PTA data and the total energydensity fraction by the BBN, and CMB data. In Tab. 1, we obtained the allowed parameter region, e.g., the  \nspectral amplitude A ≤ 0 . 39, spectral peak frequency f ∗ ∼ 10 -8 -10 -6 Hz, EoS w =0 . 44 +0 . 52 -0 . 40 , etc. Compared with those of Ref. [23], the posteriors have been significantly changed due to combining the PTA data with the BBN and CMB data and considering different scenarios of SIGWs. We revealed the importance of BBN and CMB for data analysis. Our results can be further tested with future observations. Other works relevant to the background EoS can be found in Ref. [88].  \nThrough exploring the production of PBHs during the epoch of w -domination, we have found that our scenario of w ∼ O (10 -2 ) can suppress the abundance of PBHs, leading to alleviation of the overproduction issue [23, 24]. One should note that we derived the systematic formulas for the PBH abundance, implying a consistent study of PBH production. In addition, the speed of transition from the w -domination to the radiation domination is another concern to investigate. Introducing a finite speed of transition into our model may change our results [57, 89]. However, it is extremely complicated to take the effect of transition speed into account, following the existing research approach. On the other hand, to our knowledge, there are not any works along this re-  \narch direction until now. We think that such a study warrants for another paper, which is left to future works.  \nWe acknowledge Zucheng Chen, Shi Pi and YanHeng Yu for helpful discussion. Q.H.Z is supported by the National Natural Science Foundation of China (Grant NO. 12305073). Z.C.Z. is supported by the National Key Research and Development Program of China Grant No. 2021YFC2203001 and the National Natural Science Foundation of China (Grant NO. 12005016). S.W. is partially supported by the National Natural Science Foundation of China (Grant No. 12175243), the National Key R&D Program of China No. 2023YFC2206403, the Science Research Grants from the China Manned Space Project with No. CMSCSST-2021-B01, and the Key Research Program of the Chinese Academy of Sciences (Grant No. XDPB15). X.Z. is supported by the National SKA Program of China (Grants Nos. 2022SKA0110200 and 2022SKA0110203) and the National Natural Science Foundation of China (Grants Nos. 11975072, 11835009, and 11805031). This work is supported by High-performance Computing Platform of China Agricultural University.', 'A Appendix: SIGWs due to a sudden transition of the early Universe': "In conformal Newtonian gauge, the perturbed spatially-flat Friedmann-Robertson-Walker (FRW) metric is  \nd s 2 = a 2 { -(1+2 ϕ )d t 2 + [ (1 -2 ψ ) δ ij + 1 2 h ij ] d x i d x j } , (17)  \nwhere η is a conformal time, a is a scale factor of the Universe, both ϕ and ψ are the linear cosmological scalar perturbations, and h ij denotes the second-order tensor perturbations, i.e., SIGWs. In the absence of anisotropic stress, we have ψ = ϕ , exactly. Here, the evolution of ψ is determined by the master equation,  \nψ '' +3(1+ w ) H ψ ' -w ∆ ψ = 0 , (18)  \nwhere ' denotes the derivative with respect to η , and we have adopted an adiabatic sound speed c 2 s = w with 0 ≤ w ≤ 1 for a phenomenological study. By solving Eq. (18), we obtain the Fourier modes of ψ in the w -domination and radiation domination, respectively, namely,  \nψ ( w ) k = 0 F 1 ( 7+9 w 2(1+3) w ; -1 4 wk 2 η 2 ) Ψ k , (19) ψ k = [ -√ 3 kη w cos ( kη w / √ 3 ) +3sin ( kη w / √ 3 ) ( kη w ) 3 ] ψ (A) k + [ -3cos ( kη w / √ 3 ) -√ 3 kη w sin ( kη w / √ 3 ) ( kη w ) 3 ] ψ (B) k , (20)  \nwhere the superscripts ( w ) denote the w -domination, we introduce η w = η -(1 -3 w ) η R / 2, Ψ k stands for a random variable related to the primordial curvature perturbations, and 0 F 1 is the confluent hypergeometric function. According  \nto the continuity and differentiability of ψ k , we can determine ψ (A) k and ψ (B) k as follows  \nψ (A) k = [ ( y 2 -9) √ 3 cos ( y √ 3 ) +3 y sin ( y √ 3 )] ψ ( w ) k ∣ ∣ ∣ η R + y [ √ 3cos ( y √ 3 ) + y sin ( y √ 3 )] k -1 ψ ( w ) k ' ∣ ∣ ∣ η R ,  \n(21) ψ (B) k = [ -√ 3 y cos ( y √ 3 ) +( y 2 -9)sin ( y √ 3 )] ψ ( w ) k ∣ ∣ ∣ η R -y [ y cos ( y √ 3 ) -√ 3sin ( y √ 3 )] k -1 ψ ( w ) k ' ∣ ∣ ∣ η R , (22)  \nwhere we have introduced y = kη w .  \nThe equation of motion for SIGWs has been presented in Eq. (5) in main text, where the source term S ab consisting of the linear perturbations is given by  \nS ab = 2(5+3 w ) 3(1+ w ) ∂ a ψ∂ b ψ + 4 3(1+ w ) H ( ∂ a ψ∂ b ψ ' + ∂ a ψ ' ∂ b ψ ' )+ 4 3(1+ w ) H 2 ∂ a ψ ' ∂ b ψ ' . (23)  \nTherefore, in the Fourier space, we obtain the subhorizon modes of h ij, k based on Eq. (5) in the main text, namely,  \nh ij, k = sin( kη ) kη w h (A) ij, k + cos( kη ) kη w h (B) ij, k + sin( kη ) kη w ∫ η η R d¯ η { cos( k ¯ η )Λ ab ij ( ˆ k ) S ab, k (¯ η ) ∣ ∣ w =1 / 3 [ ¯ η -1 2 (1 -3 w ) η R ]} -cos( kη ) kη w ∫ η η R d¯ η [ sin( k ¯ η )Λ ab ij ( ˆ k ) S ab, k (¯ η ) ∣ ∣ w =1 / 3 ( ¯ η -1 2 (1 -3 w ) η R )] , (24)  \nwhere S ab, k stands for the Fourier mode of S ab with the wavevector k . Due to the continuity and differentiability of h ij, k , the time-independent quantities h (A) ij, k and h (B) ij, k can be determined by SIGWs produced during the w -domination (denoted as h ( w ) ij, k ), namely,  \nh (A) ij, k = η w cos( kη R ) h ( w ) ij, k ' ∣ ∣ ∣ η = η R +[cos( kη R )+ kη w sin( kη R )] h ( w ) ij, k ∣ ∣ ∣ η = η R , (25)  \nh (B) ij, k = [ kη w cos( kη R ) -sin( kη R )] h ( w ) ij, k ∣ ∣ ∣ η = η R -η w sin( kη R ) h ( w ) ij, k ' ∣ ∣ ∣ η = η R . (26)  \nwhere we define  \nh ( w ) ij, k = ( kη ) α Y α ( kη ) ∫ η 0 d¯ η k [ π 2 ( k ¯ η ) 1 -α J α ( k ¯ η )Λ ab ij ( ˆ k ) S ab, k ] -( kη ) α J α ( kη ) ∫ η 0 d¯ η k [ ( k ¯ η ) 1 -α Y α ( k ¯ η )Λ ab ij ( ˆ k ) S ab, k ] , (27)  \nand we introduce α = -3(1 -w ) / [2(1 + 3 w )], Y α ( x ) is the Bessel function of second kind, and J α ( x ) is the Bessel function of first kind.", 'References': "- 1 H. Xu et al. , Res. Astron. Astrophys. 23 , 075024 (2023), arXiv:2306.16216 [astro-ph.HE] .\n- 2 J. Antoniadis et al. (EPTA, InPTA:), Astron. Astrophys. 678 , A50 (2023), arXiv:2306.16214 [astro-ph.HE] .\n- 3 G. Agazie et al. (NANOGrav), Astrophys. J. Lett. 951 (2023), 10.3847/2041-8213/acdac6, arXiv:2306.16213 [astro-ph.HE] .\n- 4 D. J. Reardon et al. , Astrophys. J. Lett. 951 (2023), 10.3847/2041-8213/acdd02, arXiv:2306.16215 [astro-ph.HE] .\n- 5 G. Franciolini, A. Iovino, Junior., V. Vaskonen, and H. Veermae, Phys. Rev. Lett. 131 , 201401 (2023), arXiv:2306.17149 [astro-ph.CO] .\n- 6 K. Inomata, K. Kohri, and T. Terada, Phys. Rev. D 109 , 063506 (2024), arXiv:2306.17834 [astro-ph.CO] .\n- 7 Y.-F. Cai, X.-C. He, X.-H. Ma, S.-F. Yan, and G.-W. Yuan, Sci. Bull. 68 , 2929 (2023), arXiv:2306.17822 [gr-qc] .\n- 8 S. Wang, Z.-C. Zhao, J.-P. Li, and Q.-H. Zhu, Phys. Rev. Res. 6 , L012060 (2024), arXiv:2307.00572 [astro-ph.CO] .\n- 9 L. Liu, Z.-C. Chen, and Q.-G. Huang, Phys. Rev. D 109 , L061301 (2024), arXiv:2307.01102 [astro-ph.CO] .\n- 10 K. T. Abe and Y. Tada, Phys. Rev. D 108 , L101304 (2023), arXiv:2307.01653 [astro-ph.CO] .\n- 11 R. Ebadi, S. Kumar, A. McCune, H. Tai, and L.-T. Wang, Phys. Rev. D 109 , 083519 (2024), arXiv:2307.01248 [astroph.CO] .\n- 12 D. G. Figueroa, M. Pieroni, A. Ricciardone, and P. Simakachorn, Phys. Rev. Lett. 132 , 171002 (2024), arXiv:2307.02399 [astro-ph.CO] .\n- 13 Z. Yi, Q. Gao, Y. Gong, Y. Wang, and F. Zhang, Sci. China Phys. Mech. Astron. 66 , 120404 (2023), arXiv:2307.02467 [grqc] .\n- 14 E. Madge, E. Morgante, C. Puchades-Ib'a˜nez, N. Ramberg,  \n- W. Ratzinger, S. Schenk, and P. Schwaller, JHEP 10 , 171 (2023), arXiv:2306.14856 [hep-ph] .\n- 15 H. Firouzjahi and A. Talebian, JCAP 10 , 032 (2023), arXiv:2307.03164 [gr-qc] .\n- 16 S. Wang, Z.-C. Zhao, and Q.-H. Zhu, Phys. Rev. Res. 6 , 013207 (2024), arXiv:2307.03095 [astro-ph.CO] .\n- 17 Z.-Q. You, Z. Yi, and Y. Wu, JCAP 11 , 065 (2023), arXiv:2307.04419 [gr-qc] .\n- 18 G. Ye and A. Silvestri, Astrophys. J. Lett. 963 , L15 (2024), arXiv:2307.05455 [astro-ph.CO] .\n- 19 S. A. Hosseini Mansoori, F. Felegray, A. Talebian, and M. Sami, JCAP 08 , 067 (2023), arXiv:2307.06757 [astroph.CO] .\n- 20 S. Balaji, G. Dom'enech, and G. Franciolini, JCAP 10 , 041 (2023), arXiv:2307.08552 [gr-qc] .\n- 21 J.-H. Jin, Z.-C. Chen, Z. Yi, Z.-Q. You, L. Liu, and Y. Wu, JCAP 09 , 016 (2023), arXiv:2307.08687 [astro-ph.CO] .\n- 22 B. Das, N. Jaman, and M. Sami, Phys. Rev. D 108 , 103510 (2023), arXiv:2307.12913 [gr-qc] .\n- 23 A. Afzal et al. (NANOGrav), Astrophys. J. Lett. 951 (2023), 10.3847/2041-8213/acdc91, arXiv:2306.16219 [astro-ph.HE] .\n- 24 J. Antoniadis et al. (EPTA, InPTA), Astron. Astrophys. 685 , A94 (2024), arXiv:2306.16227 [astro-ph.CO] .\n- 25 M. Drees, F. Hajkarim, and E. R. Schmitz, JCAP 06 , 025 (2015), arXiv:1503.03513 [hep-ph] .\n- 26 K. Saikawa and S. Shirai, JCAP 1805 , 035 (2018), arXiv:1803.01038 [hep-ph] .\n- 27 B. Carr, S. Clesse, J. Garc'ıa-Bellido, and F. Kuhnel, Phys. Dark Univ. 31 , 100755 (2021), arXiv:1906.08217 [astro-ph.CO]  \n.  \n- 28 I. Dalianis and C. Kouvaris, JCAP 07 , 046 (2021), arXiv:2012.09255 [astro-ph.CO] .\n- 29 G. Dom'enech, C. Lin, and M. Sasaki, JCAP 04 , 062 (2021), [Erratum: JCAP 11, E01 (2021)], arXiv:2012.08151 [gr-qc] .\n- 30 K. D. Lozanov and V. Takhistov, Phys. Rev. Lett. 130 , 181002 (2023), arXiv:2204.07152 [astro-ph.CO] .\n- 31 N. Bhaumik and R. K. Jain, Phys. Rev. D 104 arXiv:2009.10424 [astro-ph.CO] .\n- , 023531 (2021),\n- 32 M. R. Haque, D. Maity, T. Paul, and L. Sriramkumar, Phys. Rev. D 104 , 063513 (2021), arXiv:2105.09242 [astro-ph.CO] .\n- 33 T. Papanikolaou, V. Vennin, and D. Langlois, JCAP 03 , 053 (2021), arXiv:2010.11573 [astro-ph.CO] .\n- 34 H.-L. Huang, Y. Cai, J.-Q. Jiang, J. Zhang, and Y.-S. Piao, Res. Astron. Astrophys. 24 , 091001 (2024), arXiv:2306.17577 [gr-qc] .\n- 35 M. C. Johnson and M. Kamionkowski, Phys. Rev. D 78 , 063010 (2008), arXiv:0805.1748 [astro-ph] .\n- 36 M. S. Turner, Phys. Rev. D 28 , 1243 (1983).\n- 37 V. Poulin, T. L. Smith, D. Grin, T. Karwal, and M. Kamionkowski, Phys. Rev. D 98 , 083525 (2018), arXiv:1806.10608 [astro-ph.CO] .\n- 38 A. Vikman, Phys. Rev. D 71 , 023515 (2005), arXiv:astroph/0407107 .\n- 39 G. Dom'enech, V. Takhistov, and M. Sasaki, Phys. Lett. B 823 , 136722 (2021), arXiv:2105.06816 [astro-ph.CO] .\n- 40 K. N. Ananda, C. Clarkson, and D. Wands, Phys. Rev. D75 , 123518 (2007), arXiv:gr-qc/0612013 [gr-qc] .\n- 41 D. Baumann, P. J. Steinhardt, K. Takahashi, and K. Ichiki, Phys. Rev. D76 , 084019 (2007), arXiv:hep-th/0703290 [hep-th] .\n- 42 S. Mollerach, D. Harari, and S. Matarrese, Phys. Rev. D69 , 063002 (2004), arXiv:astro-ph/0310711 [astro-ph] .\n- 43 H. Assadullahi and D. Wands, Phys. Rev. D81 , 023527 (2010), arXiv:0907.4073 [astro-ph.CO] .\n- 44 J. R. Espinosa, D. Racco, and A. Riotto, JCAP 09 , 012 (2018), arXiv:1804.07732 [hep-ph] .\n- 45 K. Kohri and T. Terada, Phys. Rev. D97 , 123532 (2018), arXiv:1804.08577 [gr-qc] .\n- 46 K. Inomata, K. Kohri, T. Nakama, and T. Terada, J. Phys. Conf. Ser. 1468 , 012002 (2020).\n- 47 K. Inomata, K. Kohri, T. Nakama, and T. Terada, Phys. Rev. D 100 , 043532 (2019), arXiv:1904.12879 [astro-ph.CO] .\n- 48 F. Hajkarim and J. Schaffner-Bielich, Phys. Rev. D 101 , 043522 (2020), arXiv:1910.12357 [hep-ph] .\n- 49 G. Dom'enech, Int. J. Mod. Phys. D 29 , 2050028 (2020), arXiv:1912.05583 [gr-qc] .\n- 50 G. Dom'enech, S. Pi, and M. Sasaki, JCAP 08 , 017 (2020), arXiv:2005.12314 [gr-qc] .\n- 51 G. Dom'enech and S. Pi, Sci. China Phys. Mech. Astron. 65 , 230411 (2022), arXiv:2010.03976 [astro-ph.CO] .\n- 52 G. Dom'enech, Universe 7 , 398 (2021), arXiv:2109.01398 [gr-qc] .\n- 53 R. Flauger and S. Weinberg, Phys. Rev. D 99 , 123030 (2019), arXiv:1906.04853 [hep-th] .\n- 54 B. J. Carr and S. W. Hawking, Mon. Not. Roy. Astron. Soc. 168 , 399 (1974).\n- 55 T. Harada, C.-M. Yoo, K. Kohri, K.-i. Nakao, and S. Jhingan, Astrophys. J. 833 , 61 (2016), arXiv:1609.01588 [astro-ph.CO]  \n.  \n- 56 T. Nakama, Phys. Rev. D 101 , 063519 (2020).\n- 57 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- 58 R. Cooke, M. Pettini, R. A. Jorgenson, M. T. Murphy, and C. C. Steidel, Astrophys. J. 781 , 31 (2014), arXiv:1308.3240 [astro-ph.CO] .\n- 59 T. J. Clarke, E. J. Copeland, and A. Moss, JCAP 10 , 002 (2020), arXiv:2004.11396 [astro-ph.CO] .\n- 60 I. Dalianis and G. P. Kodaxis, Galaxies 10 , 31 (2022), arXiv:2112.15576 [astro-ph.CO] .\n- 61 G. White, L. Pearce, D. Vagie, and A. Kusenko, Phys. Rev. Lett. 127 , 181601 (2021), arXiv:2105.11655 [hep-ph] .\n- 62 S. Kasuya, M. Kawasaki, and K. Murai, JCAP 05 , 053 (2023), arXiv:2212.13370 [astro-ph.CO] .\n- 63 M. Kawasaki and K. Murai, JCAP 01 , 050 (2024), arXiv:2308.13134 [astro-ph.CO] .\n- 64 M. M. Flores, A. Kusenko, L. Pearce, Y. F. PerezGonzalez, and G. White, Phys. Rev. D 108 , 123002 (2023), arXiv:2308.15522 [hep-ph] .\n- 65 K. Harigaya, K. Inomata, and T. Terada, Phys. Rev. D 108 , L081303 (2023), arXiv:2305.14242 [hep-ph] .\n- 66 C. Pitrou, X. Roy, and O. Umeh, Class. Quant. Grav. 30 , 165002 (2013), arXiv:1302.6174 [astro-ph.CO] .\n- 67 S. Hawking, Mon. Not. Roy. Astron. Soc. 152 , 75 (1971).\n- 68 M. Khlopov and A. Polnarev, Physics Letters B 97 , 383 (1980).\n- 69 M. Khlopov, B. A. Malomed, and I. B. Zeldovich, Mon. Not. Roy. Astron. Soc. 215 , 575 (1985).\n- 70 B. Carr and F. Kuhnel, Ann. Rev. Nucl. Part. Sci. 70 , 355 (2020), arXiv:2006.02838 [astro-ph.CO] .\n- 71 M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama, Class. Quant. Grav. 35 , 063001 (2018), arXiv:1801.05235 [astroph.CO] .\n- 72 T. Papanikolaou, JCAP 10 , 089 (2022), arXiv:2207.11041 [astro-ph.CO] .\n- 73 S. Basilakos, D. V. Nanopoulos, T. Papanikolaou, E. N. Saridakis, and C. Tzerefos, Phys. Lett. B 850 , 138507 (2024), arXiv:2307.08601 [hep-th] .\n- 74 N. Aghanim et al. (Planck), Astron. Astrophys. 641 , A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .\n- 75 G. Agazie et al. (NANOGrav), Astrophys. J. Lett. 952 , L37 (2023), arXiv:2306.16220 [astro-ph.HE] .\n- 76 A. Mitridate, (2023), 10.5281/zenodo.7876430.\n- 77 A. Mitridate, D. Wright, R. von Eckardstein, T. Schroder, J. Nay, K. Olum, K. Schmitz, and T. Trickle, (2023), arXiv:2306.16377 [hep-ph] .  \n- 78 J. A. Ellis, M. Vallisneri, S. R. Taylor, and P. T. Baker, 'Enterprise: Enhanced numerical toolbox enabling a robust pulsar inference suite,' Zenodo (2020).\n- 79 W. G. Lamb, S. R. Taylor, and R. van Haasteren, Physical Review D 108 , 103019 (2023).\n- 80 M. Maggiore, Gravitational Waves. Vol. 2: Astrophysics and Cosmology (Oxford University Press, 2018).\n- 81 T. Bringmann, P. F. Depta, T. Konstandin, K. SchmidtHoberg, and C. Tasillo, JCAP 11 , 053 (2023), arXiv:2306.09411 [astro-ph.CO] .\n- 82 A. C. Sobotka, A. L. Erickcek, and T. L. Smith, Phys. Rev. D 107 , 023525 (2023), arXiv:2207.14308 [astro-ph.CO] .\n- 83 P. F. de Salas, M. Lattanzi, G. Mangano, G. Miele, S. Pastor, and O. Pisanti, Phys. Rev. D 92 , 123534 (2015),", 'arXiv:1511.00672 [astro-ph.CO] .': '- 84 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- 85 P. A. R. Ade et al. (Planck), Astron. Astrophys. 594 , A13 (2016), arXiv:1502.01589 [astro-ph.CO] .\n- 86 W. H. Press and P. Schechter, Astrophys. J. 187 , 425 (1974).\n- 87 K. Ando, K. Inomata, and M. Kawasaki, Phys. Rev. D97 , 103528 (2018), arXiv:1802.06393 [astro-ph.CO] .\n- 88 V. K. Oikonomou, Phys. Rev. D 108 , 043516 (2023), arXiv:2306.17351 [astro-ph.CO] .\n- 89 K. Inomata, K. Kohri, T. Nakama, and T. Terada, JCAP 10 , 071 (2019), [Erratum: JCAP 08, E01 (2023)], arXiv:1904.12878 [astro-ph.CO]'}