Implications of Pulsar Timing Array Data for Scalar-Induced Gravitational Waves and Primordial Black Holes: Primordial Non-Gaussianity $f_{\mathrm{NL}}$ Considered

Multiple pulsar-timing-array collaborations have reported strong evidence for the existence of a gravitational-wave background. We study physical implications of this signal for cosmology, assuming that it is attributed to scalar-induced gravitational waves. By incorporating primordial non-Gaussianity $f_{\mathrm{NL}}$, we specifically examine the nature of primordial curvature perturbations and primordial black holes. We find that the signal allows for a primordial non-Gaussianity $f_{\mathrm{NL}}$ in the range of $-4.1\lesssim f_{\mathrm{NL}} \lesssim 4.1$ (68\% confidence intervals) and a mass range for primordial black holes $m_{\mathrm{pbh}}$ spanning from $\sim10^{-5}M_{\odot}$ to $\sim10^{-2}M_{\odot}$. Furthermore, we find that the signal favors a negative non-Gaussianity, which can suppress the abundance of primordial black holes. We also demonstrate that the anisotropies of scalar-induced gravitational waves serve as a powerful tool to probe the non-Gaussianity $f_{\mathrm{NL}}$. We conduct a comprehensive analysis of the angular power spectrum within the nano-Hertz band. Looking ahead, we anticipate that future projects, such as the Square Kilometre Array, will have the potential to measure these anisotropies and provide further insights into the primordial universe.


I. INTRODUCTION
Multiple collaborations in the field of pulsar timing array (PTA) observations have presented strong evidence for a signal exhibiting correlations consistent with a stochastic gravitational-wave background (GWB) [1][2][3][4].The strain has been measured to be on the order of 10 −15 at a pivot frequency of 1 yr −1 .Though this GWB aligns with expectations from astrophysical sources, specifically inspiraling supper-massive black hole (SMBH) binaries [5], it is important to note that the current datasets do not rule out the possibility of cosmological origins or other exotic astrophysical sources, which have been explored in collaborative accompanying papers [6,7].Notably, several cosmological models have demonstrated superior fits to the signal compared to the SMBH-binary interpretation.If these alternative models are confirmed in the future, they may provide compelling evidence for new physics.
In this study, our focus lies on the cosmological interpretation of the signal, specifically the existence of scalarinduced gravitational waves (SIGWs) [8][9][10][11][12][13].This possibility had been used for interpreting the NANOGrav 12.5year dataset [14] in Refs.[15][16][17][18][19][20][21][22][23][24].It was recently revisited by the PTA collaborations [6,7], but the statistics of primordial curvature perturbations was assumed to be Gaussian.However, it was demonstrated that primordial non-Gaussianity f NL significantly contributes to the energy density of SIGWs [25][26][27][28][29][30][31][32][33].This indicates noteworthy modifications to the energy-density spectrum, which is crucial for the data analysis of PTA datasets.On the other hand, it has been shown that primordial non-Gaussianity f NL could generate initial inhomogeneities in SIGWs, leading to anisotropies characterized by the angular power spectrum [33].Related studies can be found in Refs.[34][35][36][37][38][39][40][41][42].Our analysis will encompass a comprehensive examination of the angular power spectrum within the PTA band.Moreover, this spectrum is capable of breaking the degeneracies among model parameters, particularly leading to possible determination of f NL , and playing a crucial role in distinguishing between different sources of GWB.Therefore, by interpreting the signal as originating from SIGWs, we aim to study physical implications of PTA datasets for the nature of primordial curvature perturbations, including their power spectrum and angular power spectrum.
The remaining context of this paper is arranged as follows.In Section II, we will provide a brief summary of the homogeneous and isotropic component of SIGWs.In Section III, we will show implications of the PTA data for the power spectrum of primordial curvature perturbations and then for the mass function of PBHs.In Section IV, we will study the inhomogeneous and anisotropic component of SIGWs and show the corresponding angular power spectrum in PTA band.In Section V, we make concluding remarks and discussions.

II. SIGW ENERGY-DENSITY FRACTION SPECTRUM
In this section, we show a brief but self-consistent summary of the main results of the energy-density fraction spectrum in a framework of SIGW theory.
To quantify contributions of f NL to the energy density, we express the primordial curvature perturbations ζ in terms of their Gaussian components ζ g , i.e., [69] Here, f NL represents the non-linear parameter that characterizes the local-type primordial non-Gaussianity.To simplify the subsequent analytic formulae, we introduce a new quantity as follows It is worth noting that perturbation theory requires the condition A S F 2 NL < 1, where A S will be defined later.We define the dimensionless power spectrum of ζ g as ⟨ζ g (q)ζ g (q ′ )⟩ = δ (3) In this work, we assume that ∆ 2 g (q) follows a log-normal FIG. 1. Unscaled (or equivalently, AS = 1 and F NL = 1) contributions to the energy-density fraction spectrum of SIGWs.
distribution with respect to ln q [21, 31, 70-72] where A S represents the spectral amplitude at the spectral peak wavenumber q * , and σ denotes the standard deviation that characterizes the width of the spectrum.The wavenumber q is straightforwardly converted into the frequency ν, namely, q = 2πν.Through a detailed derivation process based on Wick's theorem, we can decompose Ωgw ∼ ⟨ζ 4 ⟩ into three components depending on the power of f NL .However, the complete derivations have been simplified by employing a Feynman-like diagrammatic approach [25,28,[30][31][32][33].
Here, we present the final results

Ωgw
where we provide the analytic expressions for Ω(n) gw , which are proportional to A 2 S (A S F 2 NL ) n with n = 0, 1, 2, in Appendix A. They were evaluated using the vegas package [73], while their numerical results for σ = 1/2, 1, 2 are reproduced in Fig. 1.Specifically, Ω(0) gw corresponds to the energy-density fraction spectrum in the case of Gaussianity, while Ω(1) gw and Ω(2) gw fully describe the contributions of local-type primordial non-Gaussianity f NL .
The energy-density fraction spectrum of SIGWs at the present conformal time η 0 can be expressed as Ωgw (η, q) .
(6) In the above equation, Ω rad,0 h 2 = 4.2 × 10 −5 represents the physical energy-density fraction of radiations in the present universe [74].T and T eq correspond to the cosmic temperatures at the emission time and the epoch of matter-radiation equality, respectively.ν can be related to T , g * ,ρ (T ), and g * ,s (T ) as follows [21] Here, g * ,ρ and g * ,s represent the effective relativistic degrees of freedom in the universe, which are tabulated functions of T as provided in Ref. [75].To illustrate the interpretation of current PTA data in the framework of SIGWs, we depict Ωgw,0 (ν) with respect to ν in Fig. 2, using three specific sets of model parameters.

III. IMPLICATIONS OF PTA DATA FOR NEW PHYSICS
In this section, we investigate the potential constraints on the parameter space of the primordial power spectrum and PBHs using the NANOGrav 15-year (NG15) data.While it is possible to obtain constraints from other PTA datasets using the same methodology, we do not consider them in this study, as they would not significantly alter the leading results of our current work.
By performing a comprehensive Bayesian analysis [7], we could gain valuable insights for the posteriors of four independent parameters, i.e., F NL , A S , σ, and ν * , for which the priors are set to be F NL ∈ [−30, 30], log 10 A S ∈ [−3, 1], σ ∈ [0, 5], and log 10 (ν * /Hz) ∈ [−9, −5].Here, we also adopt the aforementioned condition of perturbativity, namely, A S F 2 NL < 1.The inference results within 68% confidence intervals are given as We can also recast Eq. ( 8) into constraints on f NL , i.e., Fig. 3 shows two-dimensional contours in log 10 A S − F NL plane at 68% (dark blue regions) and 95% (light blue regions) confidence levels.There is a full degeneracy in the sign of primordial non-Gaussianity f NL , as the energydensity fraction spectrum is dependent of only the absolute value of F NL , as demonstrated in Fig. 1.The above results indicate that the PTA observations have already emerged as a powerful tool for probing physics of the early universe.
We can further recast the constraints on the primordial curvature power spectrum into constraints on the nature of PBHs, which is characterized by their mass function.Due to significant uncertainties in the formation scenarios of PBHs (as discussed in reviews such as Ref. [42]), we adopt a simplified scenario [61] to illustrate the importance of primordial non-Gaussianity f NL .The initial mass function of PBHs is described by ) where P (ζ) represents the probability distribution function (PDF) of primordial curvature perturbations, σ g is the standard variance of the Gaussian component ζ g in the PDF, and ζ c stands for the critical fluctuation.We further find σ 2 g = ⟨ζ 2 g ⟩ = ∆ 2 g (q)d ln q = A S by considering the power spectrum defined in Eq. ( 4).Additionally, it is known that ζ c is of order O(1), with specific values of 0.7 and 1.2, as suggested by Ref. [76].
To evaluate Eq. ( 13), we devide F NL into two regimes, i.e., F NL > 0 and F NL < 0. In the case of F NL > 0, we solve the equation ζ(ζ g ) = ζ c , yielding a relation By substituting it into Eq.( 13), we gain where erfc(x) is the complementary error function.Similarly, in the case of −(4ζ c ) −1 < F NL < 0, we gain (16) In contrast, in the case of F NL < −(4ζ c ) −1 , no PBHs were formed in the early universe, since the curvature perturbations are expected to never exceed the critical fluctuation.As a viable candidate for cold dark matter, the abundance of PBHs is determined as [77] f pbh ≃ 2.5 × 10 8 β g * ,ρ (T f ) 10.75 where m pbh represents the mass of PBHs, and T f denotes cosmic temperature at the formation occasion.Roughly speaking, m pbh can be related to the horizon mass m H and then the peak frequency ν * , namely, [17] Based on Eq. ( 11), we could infer that the mass range of PBHs is the order of O(10 −5 − 10 −2 )M ⊙ .However, the inferred abundance of PBHs exceeds unity in the case of a sizable positive F NL , indicating an overproduction of PBHs.This is because the inferred value of A S is typically one order of magnitude larger than the value of A S that leads to f pbh = 1.To illustrate this result more clearly, we include into Fig. 3 two solid curves corresponding to m pbh = 10 −2 M ⊙ and f pbh = 1 in the cases of ζ c = 0.7 (purple curve) and ζ c = 1.2 (rose curve), respectively.For comparison, we mark the critical value F NL = −(4ζ c ) −1 with dotted lines.Therefore, we find that a negative F NL is capable of alleviating the overproduction of PBHs, especially when considering a sizable negative F NL , namely, F NL < −(4ζ c ) −1 , which prevents the formation of any PBHs.However, due to large uncertainties in model buildings, it remains challenging to exclude the PBH scenario through analyzing the present PTA data.
In summary, it is crucial to measure the primordial non-Gaussianity or at least determine the sign of F NL in order to assess the viability of the PBH scenario.However, it is impossible to determine the sign of F NL through measurements of the energy-density fraction spectrum of SIGWs, due to the sign degeneracy.In the next section, we will propose that the inhomogeneous and anisotropic component of SIGWs has the potential to break the sign degeneracy, as well as other degeneracies in model parameters, opening up new possibilities for making judgments about the PBH scenario in the future.

IV. SIGW ANGULAR POWER SPECTRUM
In this section, we investigate the inhomogeneities and anisotropies in SIGWs via deriving the angular power spectrum in the PTA band, following the research approach established in our previous paper [33].
The inhomogeneities in SIGWs arise from the longwavelength modulations of the energy density generated by short-wavelength modes.As discussed in Section II, SIGWs originate from extremely high redshifts, corresponding to very small horizons.However, due to limitations in the angular resolution of detectors, the signal along a line-of-sight represents an ensemble average of the energy densities over a sizable number of such horizons.Consequently, any two signals would appear identical.Nevertheless, the energy density of SIGWs produced by short-wavelength modes can be spatially redistributed by long-wavelength modes if there are couplings between the two.The local-type primordial non-Gaussianity f NL could contribute to such couplings.
Similar to the temperature fluctuations of relic photons [78], the initial inhomogeneities in SIGWs at a spatial location x can be characterized by the density contrast, which is denoted as δ gw (η, x, q), given by where the energy-density full spectrum ω gw (η, x, q) is defined in terms of the energy density, namely, ρ gw (η, x) = ρ crit (η) d 3 q, ω gw (η, x, q)/q 3 .We specifically get ω gw ∼ ⟨ζ 4 ⟩ x , where the subscript x denotes an ensemble average within the horizon enclosing x [33,34] [74].Using Feynman-like rules and diagrams, we get an expression for δ gw (η, x, q), i.e., [33] where we introduce a quantity of the form The present density contrast, denoted as δ gw,0 (q), can be estimated analytically using the line-of-sight approach [80][81][82].It is contributed by both the initial inhomogeneities and propagation effects, given by [34] δ gw,0 (q) = δ gw (η, x, q) + [4 − n gw,0 (ν)] Φ(η, x) .(22) Here, n gw,0 (q) denotes the index of the present energydensity fraction spectrum in Eq. ( 6), given by For the propagation effects, we consider only the Sachs-Wolfe (SW) effect [83], which is characterized by the Bardeen's potential on large scales Φ(η, x) = 3 5 We assume the statistical homogeneity and isotropy for the density contrasts on large scales, similar to the study of cosmic microwave background (CMB) [84].
The anisotropies today can be mapped from the aforementioned inhomogeneities.The reduced angular power spectrum is useful to characterize the statistics of these anisotropies.It is defined as the two-point correlator of the present density contrast, namely, where δ gw,0 (q) has been expanded in terms of spherical harmonics, i.e., Roughly speaking, we get C ℓ ∼ δ 2 gw,0 ∝ ⟨ζ gL ζ gL ⟩ ∼ ∆ 2 L .Detailed analysis using Feynman-like rules and diagrams was conducted in our previous paper [33].We summarize the main results as follows which can be recast into the angular power spectrum Analogous to CMB, for which the root-mean-square (rms) temperature fluctuations is determined by [ℓ(ℓ + 1)C CMB ℓ /(2π)] 1/2 , the rms density contrast for SIGWs is determined by [ℓ(ℓ + 1)C ℓ (ν)/(2π)] 1/2 , which represents the variance of the energy-density fluctuations.It is vital to note that the rms density contrast is constant with respect to multipoles ℓ, but depends on frequency bands.
In Figure 4, we present the rms density contrast as a function of gravitational-wave frequency.We also include the energy-density fraction spectrum for comparison.Roughly speaking, we find that Cℓ is the order of 10 −4 , depending on specific sets of model parameters.It is worth noting that the angular power spectrum can break degeneracies among these parameters.For instance, based on Fig. 4, we observe a coincidence in the energy-density fraction spectra for three different parameter sets.However, the angular power spectrum breaks this coincidence, particularly in the case of the sign degeneracy of f NL .This result suggests that measurements of the anisotropies in SIGWs have the potential to determine the primordial non-Gaussianity [33].Recently, 10 9  10 8 10 7 [Hz] an upper limit of Cℓ < 20% was inferred from the NG15 data [86].However, this limit is not precise enough to test the theoretical predictions of our present work.In contrast, based on Fig. 4, we anticipate that the Square Kilometre Array (SKA) program [85] will offer sufficient precision to measure the non-Gaussianity f NL .

V. CONCLUSIONS
In this study, we examined the implications of recent PTA datasets for understanding the nature of primordial curvature perturbations and primordial black holes (PBHs).Specifically, we investigated the influence of primordial non-Gaussianity f NL on the inference of model parameters, and vice versa, by analyzing the recent NG15 data.In particular, at 68% confidence level, we inferred |f NL | < 4.1, which is competitive with the constraints from measurements of CMB.Even when considering the non-Gaussianity f NL , we found that the PBH scenario is in tension with the NG15 data, except when a sizable negative f NL is considered, which can significantly suppress the abundance of PBHs.Our results indicated that the PTA observations have already emerged as a powerful tool for probing physics of the early universe and dark matter.Moreover, we proposed that the anisotropies of SIGWs serve as a powerful probe of the non-Gaussianity f NL in the PTA band.For the first time, we conducted the complete analysis of the angular power spectrum in this frequency band and found that it can effectively break potential degeneracies among the model parameters, particularly the sign degeneracy of f NL .Additionally, we explored the detectability of the anisotropies in SIGWs in the era of the SKA project.
Notes added.-During the preparation of this paper, a related study [87] appears, which examines the posteriors of NG15 data.The authors suggest that the Gaussian scenarios for SIGWs are in tension with the current PTA data at a 2σ confidence level, but non-Gaussian scenarios that suppress the abundance of PBHs can alleviate this tension.Given the significant uncertainties in the formation scenarios of PBHs (as discussed in reviews such as Ref. [42]), the main focus of our research is to simultaneously examine the energy-density fraction spectrum and the angular power spectrum of SIGWs, by incorporating the complete contributions arising from primordial non-Gaussianity f NL .We also address the importance of pri-mordial non-Gaussianity to SIGWs through a Bayesian analysis over the NG15 data.dφ 12 cos 2φ 12 J(u 1 , v 1 , x → ∞)J(u 2 , v 2 , x → ∞) (A2) The calculation for the average of the squared oscillation J(u, v, x → ∞) has been provided in Ref. [33], as well as in earlier studies referenced in Refs.[12,13,30,31], i.e., J(u i , v i , x → ∞)J(u j , v j , x → ∞) = 9 1 − s 2 i 1 − s 2 j t i (t i + 2) t j (t j + 2) s 2 i + t 2 i + 2t i − 5 s 2 j + t 2 j + 2t j − 5 8 (−s i + t i + 1) 3 (s i + t i + 1) 3 (−s j + t j + 1) 3 (s j + t j + 1) The equations presented in this appendix can be utilized to numerically calculate the energy density of SIGWs in a self-consistent manner.

2 ,
FIG.2.Energy-density fraction spectra of SIGWs for different sets of independent parameters.The NG15 data are also shown for comparison.