Testing Stochastic Gravitational Wave Signals from Primordial Black Holes with Optical Telescopes

Primordial black holes (PBHs) can constitute the predominant fraction of dark matter (DM) if PBHs reside in the currently unconstrained"sublunar'' mass range. PBHs originating from scalar perturbations generated during inflation can naturally appear with a broad spectrum in a class of models. The resulting stochastic gravitational wave (GW) background generated from such PBH production can account for the recently reported North American Nanohertz Observatory for Gravitational Waves (NANOGrav) pulsar timing array data signal, and will be testable in future GW observations by interferometer-type experiments such as Laser Interferometer Space Antenna (LISA). We show that the broad mass function of such PBH DM is already being probed by Subaru Hyper Suprime-Cam (HSC) microlensing data and is consistent with a detected candidate event. Upcoming observations of HSC will be able to provide an independent definitive test of the stochastic GW signals originating from such PBH DM production scenarios.

Primordial black holes (PBHs), formed in the early Universe prior to any galaxies and stars, are a viable dark matter (DM) candidate (e.g. ). PBHs could also play a central role in a variety of astrophysical phenomena, such as progenitors [36][37][38][39][40][41][42][43] for the Laser Interferometer Gravitational-Wave Observatory (LIGO) gravitational wave (GW) events [44][45][46], seeds for formation of super-massive black holes [37,47,48], as well as the source of new signals [43,49,50] from compact star disruptions due to PBH capture, among others. Depending on the formation time, the resulting PBHs can span many decades of orders of magnitude in mass. PBHs formed with mass above the Hawking evaporation limit of ∼ 10 15 g survive until present and contribute to the DM abundance. PBHs residing in the smaller "sublunar" mass range of ∼ 10 −16 − 10 −10 M can constitute the entirety of DM [51][52][53].
Recently, the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) collaboration has reported a putative signal from a 12.5 years analysis of pulsar timing array data [54]. The signal is consistent with a stochastic GW background (SGWB) with a GW strain amplitude of A ∼ 10 −15 and a frequency of f ∼ 10 −9 Hz, resulting in a radiation-dominated postinflationary Universe SGWB of Ω GW h 2 ∼ 10 −10 with h being the dimensionless Hubble parameter. While the reported NANOGrav signal is in tension with their own previous 11 years data analysis [55], as well as with the Parkes Pulsar Timing Array (PPTA) [56] analysis, it has been argued that previous pulsar timing array results have been over-interpreted in astrophysical terms, leading to suppressed potential signals. In the nHz frequency range, a significant source of GWs is expected to originate from supermassive black hole binaries (e.g. [57,58]). However, predictions of this SGWB contribution suffer from significant uncertainties. The reported signal has been also interpreted in terms of cosmic strings [59][60][61][62], domain walls [63], phase transitions [64][65][66][67], and tensor perturbations from inflation [68].
During PBH production in the early Universe from enhanced curvature perturbations, a SGWB will also be generated at second-order. This allows to naturally associate the reported NANOGrav observations with PBH formation [69][70][71]. Particularly intriguing is the possibility of a broad PBH mass function [72] arising in a class of motivated models [73][74][75][76] that can account for the entirety of DM in the open lower-mass parameter space window as well as SGWB signal simultaneously [70].
In this work, we demonstrate that a broad PBH mass function, responsible for DM and that simultaneously can account for the reported SGWB signal, is already being probed by Subaru Hyper Suprime-Cam (HSC), and will be definitively tested with upcoming observations of HSC and other optical telescopes. In contrast to other measurements designed to probe the SGWB, such as those of the future Laser Interferometer Space Antenna (LISA), this constitutes an independent test of the reported signal. The potential of HSC to probe PBH models with a broad mass function has been recently stressed in the context of other general PBH formation scenarios [34].
We now briefly describe PBH and SGWB production from a broad mass function that could arise in a class of models such as [73][74][75][76], following Ref. [70]. We consider the broad and flat power spectrum of curvature perturbations of the form [72] where k s k l , Θ is the Heaviside step function and A ζ is the amplitude. This also gives the spectrum of overdensity δ perturbations P δ (k).
PBHs form when a sufficient overdense region, corresponding to the density fluctuations with a sufficiently large amplitude at a certain scale, enters the Hubble horizon. Hence, given the perturbation spectrum, the energy density in PBHs at formation can be found from Press-Schechter formalism as [26] where δ c is the critical density contrast threshold for collapse, and near the collapse threshold one has to take into account the scaling of the PBH mass with respect to the horizon and γ c 0.36 [91][92][93][94][95]. Here, P (δ, σ PBH ) is the probability distribution of density fluctuations entering the horizon, assumed to be Gaussian with a variance where W (k, R H ) is a window function for smoothing over the horizon scale R H ∼ 1/aH, and T (k, R H ) is a transfer function smoothing over sub-horizon modes [72]. After matter-radiation equality the DM fraction consisting of PBHs can be expressed as a mass function (e.g. [96]) (4) where M eq 3 × 10 17 M is the horizon mass at equality, and Ω DM is cold dark matter density. The total contribution of PBHs to DM f PBH,tot is found by integrating f PBH over ln M PBH .
Enhanced curvature perturbations, responsible for PBH generation, also induce GWs at second order (see e.g. Ref. [97]). The contribution today of this SGWB is given by [70,98,99] where k = 2πf , Ω r,0 is the radiation abundance today, I(x, y) is the kernel function, and the integration region S covers x > 0 and |1 − x| ≤ y ≤ 1 + x. The factor c g accounts for the change in number of degrees of freedom of thermal radiation during time evolution. In Ref. [70] it was shown, by fitting Eq. (5), that the reported NANOGrav signal can result from the perturbation spectrum of Eq. (1), peaking at M PBH ∼ 10 −14 M and responsible for PBH DM, with A ζ 5.8 × 10 −3 , k s = 10 9 k l 1.6 Hz, and f PBH,tot = 1. We display this result in Fig. 1. The peak of the mass function denotes the horizon mass when the shortest scale ∼ 1/k s re-enters. The sub-dominant peak around solar-mass results from a change in the equation of state at the QCD phase transition 2 [96]. The PBH mass function between peaks follows ∼ M −1/2 functional dependence 3 . Since the SGWB spectrum from PBHs extends over many orders of magnitude, this scenario can be tested with future GW observations such as those of LISA.
Here we suggest that optical telescopes can provide an independent test of this explanation for the NANOGrav signal and will provide a definitive probe in the near future. Even though the mass function of PBHs constituting DM peaks at much smaller masses, where microlensing effects are negligible, the broad mass function has a tail overlapping with the HSC sensitivity range. This is actually consistent with HSC observations of Andromeda galaxy (M31) [77], which reported a candidate event consistent with PBH at f PBH (M ∼ 10 −9 M ) ∼ 10 −2 . While these observations lend credibility to the hypothesis of PBHs being the source of the SGWB detected by NANOGrav, additional HSC observations will be able to test this scenario.
Following Ref. [34], we estimate the reach of HSC, exploiting results from HSC Monte Carlo simulations as well as their analysis tools [77]. The expected number of observed microlensing events reads PBH consistent with the HSC candidate event reported after 7 hours of observations [77]. The thick purple line denotes the best fit normalization. The HSC constraint on the total PBH DM fraction f PBH,tot , assuming a monochromatic mass function, is also shown in blue. We account for the finite source effects [52,78]. The green line shows the mass function which is consistent with the reported NANOGrav SGWB signal and normalized to be f PBH,tot = Ω PBH /Ω DM = 1 [70].
[Right] Forecast of the exclusion region for the normalization of a PBH mass function ∝ M −1/2 PBH spectrum (shown in red), combining existing 7 hours of observation and additional 4 hours of observation (assuming 1 and 0 event for each, respectively). The blue region shows the constraint on the total PBH DM fraction f PBH,tot , assuming a monochromatic mass function, obtained with the same observation time (4 + 7 = 11 hours in total). Constraints from extragalactic γ-rays from BH evaporation [79] (additional constraints in this region have been also recently suggested [80][81][82][83][84]), microlensing Kepler data [85], MACHO/EROS/OGLE microlensing [86], the accretion effects on the CMB observables [87][88][89] as well as dwarf galaxy heating [90] are also displayed.
Here, (Ω PBH /Ω DM ) = f PBH,tot is the mass fraction of DM in the form of PBHs, dN event /d log(t FWHM ) is the expected differential number of PBH microlensing events per logarithmic interval of the fullwidth-at-halfmaximum (FWHM) microlensing timescale t FWHM for a single star in M31, dN s /dm r is the luminosity function of source stars in the photometric r-band magnitude range [m r , r + dm r ], and (t FWHM , m r ) is the detection efficiency (which quantifies the probability that a microlensing event for a star with magnitude m r and light curve timescale t FWHM is detected by HSC event selection procedures). We normalize the PBH mass function f PBH (M ) by enforcing the condition that PBHs constitute the entirety of DM,  Fig. 1. Notice that as the PBH spectrum is not monochromatic, the allowed range does not reach the differential HSC exclusion region for a monochromatic mass function.
Based on Poisson statistics, detecting a single event after the 7 hours of observations already carried out by HSC is compatible with the PBH DM hypothesis and a spectrum scaling as M −1/2 at ∼ 16% CL. New detections are expected with only 2.4 hours of observations. The scenario of f PBH,tot = 1 and spectrum f PBH (M PBH ) given by Eq. (4) can be excluded at a 2-σ level (95% CL) with additional 4 hours of observation, when combined with the existing 7 hours of observation and assuming null event detection in future observations. The red shaded region in right panel of Fig. 1 is the exclusion region after 11 hours of observation in total.
In conclusion, we have highlighted that SGWB signals from PBHs with an extended mass function and comprising DM can be well probed by optical telescopes, independently of other GW-specific experiments. In particular, SGWB from PBH DM formation due to a broad perturbation spectrum and consistent with the signal detected by the NANOGrav collaboration is already being tested by HSC. Near future observations by HSC and other optical telescopes will provide a definitive indepen-dent test of this possibility.
The work of A.K., V.T. and E.V. was supported by the U.S. Department of Energy (DOE) Grant No. DE-SC0009937. A.K. was also supported by the World Premier International Research Center Initiative (WPI), MEXT, Japan. S.S. is supported by International Graduate Program for Excellence in Earth-Space Science (IGPEES), World-leading Innovative Graduate Study (WINGS) Program, the University of Tokyo. This research was also supported in part by the National Science Foundation under Grant No. NSF PHY-1748958. This work was supported in part by JSPS KAKENHI Grant Numbers 15H03654, 15H05887, 15H05893, 15H05896, 15K21733, 19H00677, 19H01895 and 20H04727.