NANOGrav results and LIGO-Virgo primordial black holes in axion-like curvaton model

We discuss a possible connection between the recent NANOGrav results and the primordial black holes (PBHs) for the LIGO-Virgo events. In particular, we focus on the axion-like curvaton model, which provides a sizable amount of PBHs and GWs induced by scalar perturbations around the NANOGrav frequency range. The inevitable non-Gaussianity of this model suppresses the induced GWs associated with PBHs for the LIGO-Virgo events to be compatible with the NANOGrav results. We show that the axion-like curvaton model can account for PBHs for the LIGO-Virgo events and the NANOGrav results simultaneously.

Introduction.-TheNANOGrav collaboration has recently reported the results of the 12.5-year pulsar timing observation, which shows strong evidence of a stochastic process around O (nHz) with a common amplitude and a common spectral index across all pulsars [1].Although quadrupolar spatial correlations of the stochastic process, which should exist for gravitational wave (GW) signals, have not been found yet, it is worth investigating potential implications of the process in terms of stochastic GWs.
There are already some works discussing the NANOGrav results in the context of the induced GWs associated with PBHs.Ref. [32] showed that the NANOGrav signal is too small to be compatible with the induced GWs associated with PBHs for the LIGO-Virgo events in their setups, but on the other hand it can be consistent with the PBHs as the primordial seeds of supermassive BHs.Also, possible connections of the NANOGrav results with the PBHs for DM [33], O (1)M BHs [34], and planet-mass BHs [35] were discussed (see also Ref. [36]).
The main aim of this letter is to show that the axion-like curvaton model [19,37,38] can generate PBHs for the LIGO-Virgo events that are consistent with the NANOGrav results.The key difference from the setups in Ref. [32] is that we take into account the primordial non-Gaussianity of curvature perturbations, which is inevitably produced in the curvaton model.The primordial non-Gaussianity suppresses the induced GWs with the PBH abundance fixed, which enables our model to explain PBHs for the LIGO-Virgo events and the NANOGrav results simultaneously.
Axion-like curvaton model.-First,we summarize the basic properties of the axion-like curvaton model (see Refs. [19,38,39] for detail).This model was originally introduced in the framework of supersymmetry [40,41].The relevant field corresponds to the flat direction, denoted as Ψ = ϕe i θ / 2. After taking into account the supergravity effect, we can approximate the effective potential for ϕ as where H is the Hubble parameter during the inflation, c is the coefficient coming from the supergravity effect, and f is a symmetry breaking scale.We assume that c is of O (1) and the initial value of ϕ is much larger than f , which makes ϕ roll down the potential to the minimum f during inflation.Once ϕ reaches the minimum f , the curvaton is given by σ ≡ θ f .Similarly to the QCD axion, we assume that a certain non-perturbative dynamics generates a mass to the curvaton at some point after reheating.Its potential is given by where Λ is the scale of the non-perturbative dynamics and the curvaton mass is obtained from m σ ≡ Λ 2 / f .The second equality is valid when the curvaton is near its minimum.
Once the curvaton acquires its mass, it starts to oscillate around the minimum and behaves as the non-relativistic matter component.During a radiation-dominated (RD) era, its energy fraction in the total energy continues to grow.
In this work, we consider the case where the curvaton decays to radiation before it dominates the Universe.
In the axion-like curvaton model, the power spectrum of curvature perturbations on small scales, relevant to the PBH production, can be approximated by where η is the conformal time and k * is the inverse of the horizon scale at the arrival of ϕ to the minimum f .Note that, on large scales observed by CMB, the power spectrum is dominated by the inflaton fluctuations and hence different from Eq. (3).A ζ is given by where H | η=1/k * denotes its value at η = 1/k * and θ is the misalignment angle of the curvaton.r stands for the ratio between the energy densities of the curvaton and radiation until the curvaton decay, and it is frozen after its decay, i.e., where the subscript D indicates the value at the curvaton decay.The tilt n 1 is related to the coefficient c in Eq. ( 1): which originates from the evolution of the radial direction ϕ during the inflation.Since a typical size of fluctuations of the phase-direction field at the horizon exit is H /(2π) regardless of the value of ϕ at that time, the fluctuation of the misalignment angle is given by δθ = H /(2πϕ exit ), with ϕ exit being the value at the horizon exit.The curvaton perturbation, δσ = f δθ, gets more suppressed for a larger value of ϕ exit corresponding to perturbations at a larger scale, i.e., a smaller k.On the other hand, the other tilt n 2 stems from the time-dependence of H , which was not taken into account in our previous work [19].As we consider slow-roll inflation, n 2 − 1 should be smaller than unity.Note that n 2 is the tilt of the small-scale perturbations and therefore can be different from that on larger scales, determined by the Planck collaboration [42]. 1BH abundance.-Next,we introduce equations connecting the power spectrum and the PBH abundance.The scale of the perturbation for the PBH production is related to the PBH mass as [7], where γ is the fraction of the PBH mass in the horizon mass at the production and g is the effective degrees of freedom of radiation at that time.As fiducial values, we take γ = 0.2 [4] and g = 10.75.We do not include the effect of the critical collapse [43][44][45][46] for simplicity because Ref. [47] implies that it does not modify the PBH abundance so much.We adopt the Press-Schechter formalism throughout this work. 22 In this formalism, the production rate of PBHs at the time of formation, β, is given by where δ c is the threshold of the PBH formation.Here, we have assumed the perturbations follow the Gaussian statistics and we will explain how to take into account the non-Gaussianity later.The variance of the smoothed density contrast is denoted by σ 2 .(Hereafter, it does not indicate the curvaton field.)The variance is obtained from where M denotes the PBH mass, R = 1/k(M ) is the smoothing scale for the PBH production of mass M , W is a window function in Fourier space, and D(x) is the transfer function of the gravitational potential, Φ, with the normalization of Φ = − 3 2 ζ in the superhorizon limit.Remember that the wavenumber is related to the PBH mass through Eq. ( 7).We take the following transfer function regardless of the scales: Although the transfer function for k k D (≡ 1/η D ) could be modified [19], we expect such modification does not change our results so much because the main contribution comes from the scales of k * < k < k D .In this letter, we take the real-space top-hat window function as a fiducial example and take δ c = 0.51 for a pure RD era [51] (see also Ref. [52] for the discussion on the window function dependence).We also take into account a small modification of δ c during the QCD phase transition based on the results in Ref. [53].Note that, with the real-space top-hat window function, σ 2 (M ) P ζ (k) holds for a scale-invariant spectrum [19,52].The current fraction of PBHs in DM over a logarithmic interval is where Ω DM is the current energy density parameter of DM and Ω PBH (M ) is the density parameter of PBHs whose masses are smaller than M .
Induced gravitational waves.-Here,we briefly review the fundamental equations of GWs induced by scalar perturbations.We take the conformal Newtonian gauge where a is the scale factor, Φ and Ψ are scalar perturbations, and h i j is the tensor perturbation corresponding to GWs.Since we focus on the early Universe, we may assume the perfect fluid condition, Φ = Ψ.See Refs.[54][55][56][57][58][59][60][61][62][63][64][65] for the recent discussion on the gauge (in)dependence of the induced tensor perturbations.
The energy density parameter of the induced GWs is where H is the conformal Hubble parameter.P h is the time-averaged power spectrum of the induced GWs, given by [19] The integrand I (u, v, k, η) is obtained from where x = kη and x = k η.Since the curvaton energy is assumed to be subdominant in the Universe in our scenario, we may use the Green function in a RD era, namely where Θ(x) is the Heaviside step function.The function F in Eq. ( 16) can be expressed by means of the transfer function of the gravitational potential, T (η, k): where the prime means a derivative with respect to η.A concrete expression of T (η, k) is given in Sec.IV of Ref. [19].
For scales that enter the horizon after the curvaton decay, T (η, k) is almost the same as D(kη), given in Eq. ( 11), except for the normalization.On the other hand, for scales that enter the horizon before the curvaton decay, we need to pay attention to some issues overlooked in our previous work [19].
In the previous work, we implicitly assumed that the curvaton decay occurs instantaneously, and impose continuity on T and T before and after the curvaton decay.However, in a realistic situation, the curvaton decay occurs gradually on a time scale of its decay rate and such an approximation might not be valid, given the results in Ref. [66].Although the situation of this work is not exactly the same as that of Ref. [66], some suppression of gravitational potential during the curvaton decay might occur similarly.The study of this effect is beyond the scope of this letter and we use the expression of T (η, k) in Ref. [19] for simplicity.For this reason, we should keep in mind that the GWs spectrum of this work could possibly be more suppressed on scales k > k D .During a RD era, the density parameter in Eq. ( 14) finally asymptotes to a constant value after the scalar perturbations on the peak scale (k k * ) enter the horizon because the induced GWs behave as radiation.Taking into account the subsequent matter-dominated and dark-energydominated era, we obtain the following expression of the density parameter at present [19]: 75 where η const is the conformal time when Ω GW becomes constant at the RD era, g const is degrees of freedom at η const , and Ω r,0 is the current radiation density parameter.Non-Gaussianities.-In the axion-like curvaton model, the curvature perturbations that produce PBHs inevitably have the primordial non-Gaussianity.Besides, the nonlinear relation between the density perturbations and the curvature perturbations induces the intrinsic non-Gaussianity for the density perturbations [67][68][69].The non-Gaussianity is characterized by where δ g follows the Gaussian statistics and σ 2 is its variance.Using the equations in Ref. [67], we numerically evaluate the skewness µ 3 for the scale-invariant power spectrum as where the first term is the contribution from the primordial non-Gaussianity and the second term is from the intrinsic non-linear relation.In the curvaton model, f NL is related to r D as The probability distribution function (PDF) of the Gaussian part is given by Solving Eq. ( 20) with respect to δ g , we obtain The PDF of the non-Gaussian part is now written as Then, we can express β as δ g+ (δ c ) P G (δ g )dδ g for µ 3 < 0 . ( Since Ω GW is roughly proportional to P 2 ζ , we define the suppression factor of Ω GW as Here P ζ | µ 3 =0 is the power spectrum with µ 3 = 0 which gives the same β as that for P ζ with non-zero µ 3 in the numerator.We take into account the non-Gaussianity effect on the induced GWs multiplying this suppression factor to the Ω GW that is calculated with P ζ | µ 3 =0 .See also Refs.[22,25] for more detailed discussion. Results.-As a fiducial example, we take the following parameter values:  (28).The shaded regions are constrained by EROS/MACHO [70,71], caustic crossing events [72], CMB anisotropy with the assumption of the spherical and the disk accretion onto PBHs [73], and GWs from mergers [74].See also Ref. [12] for other constraints.
Figure 1 shows the PBH mass spectrum with this parameter set.We can see f PBH ∼ 10 −3 at M = 30M , which is consistent with the merger rate estimated by the LIGO-Virgo collaboration [10,11].The bump around M ∼ 1M is due to the suppression of the PBH threshold during the QCD with the tilt of Ω GW ∝ f −0.4 .The constraints are also shown from the previous results of EPTA [75] and PPTA [76], and the future sensitivity of SKA [77].The black dotted lines show the spectra on k > k D , which might be suppressed more because the curvaton decay is not instantaneous (see below Eq. ( 18)).
phase transition [53].Although the mass spectrum seems inconsistent with the constraint from the CMB anisotropy in Ref. [73], we should keep in mind that it still has uncertainties in accretion models, as the reference itself mentions.
Figure 2 shows the spectra of GWs induced by the scalar perturbations that realize the mass spectrum in Fig. 1.Note again that the non-Gaussianity, dependent on r D , changes the relation between the abundance of PBHs and the power spectrum of curvature perturbations.From this figure, we can see that if we take r D 0.6, the induced GWs are consistent with the 2σ region of the NANOGrav stochastic process.On the other hand, for a small r D 0.6, corresponding to f NL 2.1, the non-Gaussianity suppresses the induced GWs too much to explain the NANOGrav results.For comparison, we also provide the induced GWs neglecting the primordial non-Gaussianity from the curvaton.We can see that, if there is no primordial non-Gaussianity, the induced GWs are too large to be compatible with the NANOGrav results.This is qualitatively consistent with Ref. [32].
Conclusion.-In this letter, we discuss the implications of the recent NANOGrav results in terms of PBHs that explain the BHs detected by the LIGO-Virgo collaboration.Such PBHs for the LIGO-Virgo events predict large GWs, which are induced by scalar perturbation, around the NANOGrav frequency.We pursue the possibility to explain the NANOGrav results by these induced GWs.We study the axion-like curvaton model as a concrete example, which inevitably produces the primordial non-Gaussianity.We have demonstrated that the curvaton model can account for the LIGO-Virgo PBH scenario and the NANOGrav results at the same time, which is enabled by the non-Gaussianity sup-pressing the induced GWs with the PBH abundance fixed.Our result implies that the NANOGrav results could be a signal from the LIGO-Virgo PBHs in the case where the primordial non-Gaussianity exists.

FIG. 2 .
FIG.2.The spectra of the induced GWs with r D = 0.8 (black, top, thick) ∼ 0.4 (bottom, thin).The difference of r D between two adjacent lines is 0.1.We also show the spectrum without the non-Gaussianity from curvaton (µ 3 /σ = −9/4, brown dot-dashed).The green-shaded region shows the 2σ region of the NANOGrav results with the tilt of Ω GW ∝ f −0.4 .The constraints are also shown from the previous results of EPTA[75] and PPTA[76], and the future sensitivity of SKA[77].The black dotted lines show the spectra on k > k D , which might be suppressed more because the curvaton decay is not instantaneous (see below Eq. (18)).