A Unified Origin of Primordial Black Hole Dark Matter and Nanohertz Gravitational Waves

Abstract

Recent high-cadence observations by Subaru-HSC have identified a population of ultrashort-timescale microlensing events, providing a compelling window for planet-mass primordial black holes (PBHs) to constitute the entirety of dark matter. In this Letter, we demonstrate that this PBH population and the nanohertz stochastic gravitational-wave (GW) background reported by pulsar timing arrays (PTAs) can be naturally unified by a single primordial origin: a broad, nearly-flat enhancement of the curvature power spectrum with an amplitude of $O(10^{-2})$. The resulting PBH mass function spans the planet-to-solar mass range, while remaining consistent with all current observational constraints. This unified PBH--induced-GW framework makes concrete multi-messenger predictions, which can be decisively scrutinized by forthcoming microlensing surveys, next-generation PTAs, space-borne interferometers, precision astrometry, and laser ranging experiments.

Figures (5)

• Figure 1: The blue and red shaded contours show the 68% (1$\sigma$) and 95% (2$\sigma$) C.L. allowed regions of PBH abundance when all 12 microlensing candidates ("S12") or the 4 secure candidates ("S4") of Subaru-HSC are attributed to PBHs, respectively Sugiyama:2026kpv. The 95% C.L. allowed region from OGLE 6 events ("O6") is also shown as the gray shaded region Niikura:2019kqi. PBH mass functions are drawn for $f_\mathrm{tot}=1$ with $\langle M\rangle$ set at the 2$\sigma$ blue contour boundaries, while $k_\mathrm{min}$ is fitted to NANO-Grav data at 2$\sigma$. Mass functions without an IR cutoff are also shown for comparison (thin dotted). Other observational constraints include Subaru-HSC 2019 (light red dashed) Niikura:2017zjd, EROS (brown dashed) EROS-2:2006ryy, OGLE relaxed (purple shaded) and strict (purple dot-dashed) Mroz:2024mse, LIGO asymmetric merger (orange wiggling shaded) LIGOScientific:2025vwc, LIGO merger (orange shaded) Andres-Carcasona:2024wqk, LIGO stochastic background (orange dashed) Boybeyi:2024mhp, LIGO direct search (orange dotted) Kacanja:2026byy, and CMB disk accretion (steel blue shaded) Serpico:2020ehh.
• Figure 2: Posterior distributions of $k_{\min}$ (left) and $k^{\rm LC}_{\max}$ (right). The shaded regions indicate the 2-$\sigma$ bounds, and the dashed lines denote the best-fit values. In the right panel, the black curve corresponds to the Subaru-HSC sample with 4 secure candidates, while the gray curve represents the full sample of 12 candidates.
• Figure 3: GW spectrum $\Omega_\mathrm{GW}$ induced by a top-hat curvature power spectrum. The dark-blue line and shaded region represent the spectrum for the best-fit parameter and the 2-$\sigma$ C.L., respectively, inferred from the NANOGrav 15-yr data NANOGrav:2023ctt and the Subaru-HSC sample of all candidates Sugiyama:2026kpv, while the light-blue shaded region assumes 4 secure candidates. The first 14 bins of the NANO-Grav signal NANOGrav:2023ctt is shown in green violins together with some future observational projects, including the Roman Telescope (Seagreen) Wang:2022sxnPardo:2023cag, binary resonances (BRs) (Lightblue) Foster:2025nzfBlas:2026xws shown in LLR (dashed), eLO (dot-dashed), and eSLR (dotted) configurations, LISA (Red) LISA:2017pwj, Tai-ji (Orange) Luo:2019zal, and Tian-Qin (Pink) Luo:2025ewp. The gray-shaded region is the bound from relativistic DOF from CMB/BBN Cyburt:2004ycBinetruy:2012zePlanck:2018vygArbey:2021ysgGrohs:2023voo.
• Figure E1: The threshold of density contrast $\delta_\mathrm{c}(M)$ for different shape factor $q$, as a function of horizon mass $M_H$. See the text for details.
• Figure E2: PBH mass function (upper panel) and induced GW spectrum (lower panel) with $f_\mathrm{tot}=1$ and $\langle M\rangle=4\times10^{-7}~M_\odot$, for different threshold $\delta_\mathrm{c}$. The parameters are shown in the table. In PTA data, there is a pronounced correlation between the power spectrum amplitude $\mathcal{A}_{\mathcal{R}}$ and the IR cutoff $k_{\min}$, which are adjusted accordingly for different $\delta_\mathrm{c}$ (see text for details). The thin dotted curves are the mass functions and induced GW spectra without an IR cutoff, which exhibit strong tension with LVK constraints and PTA observations, respectively.