Constraints on Primordial Black Hole Dark Matter from the Stochastic Gravitational-Wave Background Gokdeniz Baydar∗ University of Minnesota, Twin Cities (Dated: April 20, 2025) This work investigates the contribution of primordial black hole (PBH) binaries to the stochastic gravitational-wave background (SGWB) using data from the LIGO–Virgo–KAGRA (LVK) collabo- ration. A Bayesian inference framework is employed to model the PBH population with a log-normal mass distribution, incorporating early binary formation and suppression effects from environmental interactions such as tidal disruptions and clustering. Observational limits from the LVK O1–O3 runs are used to constrain the PBH dark matter fraction. The results place stringent upper bounds on PBH abundances in the 0.05–10M⊙ mass range and demonstrate that suppression effects signif- icantly reduce the allowed parameter space. These refined constraints are essential for assessing the role of PBHs as dark matter candidates and emphasize the importance of future observations with instruments such as LIGO O4, LISA, and DECIGO. Keywords: Primordial black holes, stochastic gravitational-wave background, LIGO, suppression effects I. INTRODUCTION Gravitational waves (GWs), ripples in spacetime predicted by Einstein’s General Theory of Relativity, have created a new frontier in cosmology. The first direct detection of GWs by the LIGO-Virgo-Kagra Collaboration (LVK) in 2015 made direct investiga- tion of this phenomenon possible. Among various po- tential sources, primordial black holes (PBHs) have emerged as candidates contributing to the stochastic gravitational-wave background (SGWB). SGWB is a pervasive and continuous GW signal arising from nu- merous unresolved sources distributed throughout the universe. Primordial black holes, theoretically proposed in 1966 by Yakov Zeldovich and Igor Novikov, are thought to form from the collapse of large curvature perturbations generated during inflation in the early universe [3, 4]. Unlike stellar-mass black holes orig- inating from collapsed massive stars, PBHs span an extensive mass range and can form independently of stellar evolutionary processes. Their versatile mass spectrum, potentially ranging from subatomic-scale black holes to supermassive ones, makes PBHs also a viable candidate for dark matter. Thus, constraining their abundance and mass distribution is crucial, not only for GW astrophysics but also for understanding cosmological models and the nature of dark matter. The SGWB generated by PBHs predominantly arises from their formation mechanisms, early- universe interactions, and binary coalescence. PBH binaries can form dynamically through gravitational capture events in dense environments or through pri- mordial gravitational decoupling during the universe’s radiation-dominated era [2, 8]. Each formation path- way leads to distinct observable signatures in the SGWB, represented by its spectral shape and ampli- tude. As a result, careful analysis of SGWB data from current and future GW observatories offers constraints on PBH populations. ∗ Correspondence email address: bayda002@umn.edu The LIGO collaboration has placed upper bounds on the SGWB in the 20–100Hz frequency band using data from their first three observing runs (O1–O3). This has limited the contribution of PBHs to the to- tal dark matter abundance [1]. Interpreting these con- straints requires a framework that includes the PBH mass distribution, binary formation channels, redshift evolution, and suppression effects [2, 5]. This study adopts a minimal modeling approach based on a log-normal PBH mass function. The GW energy spectrum ΩPBH(f) is computed by integrating the merger rate over redshift and mass, including rel- evant suppression factors. The predictions from this model are then rigorously compared with the LIGO upper limits to constrain the dark matter fraction of PBHs, denoted as fPBH. The computational framework used here is based on the existing methodologies developed by Töre Boy- beyi, who also supervised this project. By using sev- eral statistical methodologies, including Bayesian in- ference and Markov Chain Monte Carlo (MCMC) sim- ulations, this work quantitatively assesses constraints on the PBH fraction of dark matter. Key results in- clude a posterior contour plot in the µ–fPBH plane, portraying the allowed parameter space, and a con- straint curve illustrating upper limits on fPBH across different PBH masses [1, 8]. This analysis aims to improve our understanding of PBHs and the parameters that constrain them, while also helping to guide observational strategies and data-analysis efforts for future gravitational-wave detectors. II. THEORETICAL FRAMEWORK The abundance and mass distribution of PBHs are key to cosmological models [3, 4]. The present-day fraction of dark matter composed of PBHs, fPBH(M), can be calculated from the initial mass fraction, β(M), which quantifies the fraction of the mass of the uni- verse collapsing into PBHs at formation. The relation between these two fractions is expressed as follows: mailto:bayda002@umn.edu 2 fPBH(M) = 2.7× 108 ( γ 0.2 )1/2 ( g∗,form 10.75 )−1/4 ( M M⊙ )−1/2 β(M) (1) Here, γ represents the collapse efficiency factor, typi- cally around 0.2 for standard formation scenarios, and g∗,form is the effective number of relativistic degrees of freedom at the epoch of PBH formation. The initial mass fraction, β(M), can be approximated by: β(M) ≈ γ√ 2π ν(M) exp [ −ν(M)2 2 ] , (2) where ν(M) ≡ δth/σM represents the ratio of the threshold density contrast δth required for collapse, to the variance of the density perturbations σM at a mass scale M [10]. This parameter directly influences the likelihood of collapse and thus the initial abundance of PBHs. To characterize PBH populations, a log-normal mass distribution is used and given by: ψ(m) = 1√ 2πσm exp ( − (lnm− lnµ)2 2σ2 ) , (3) where µ is the central mass around which PBH masses cluster, and σ defines the width of the distribu- tion. This form naturally emerges from enhanced cur- vature perturbations in numerous inflationary models [4, 6]. Its shape allows for both narrow and broad mass spectra, making it a flexible tool for exploring a wide range of formation scenarios. In contrast, a monochromatic mass function in which all PBHs have the same mass is a highly idealized case. Although this method is much more simplistic, it often leads to overly sharp constraints which do not reflect com- plex, extended mass distributions predicted by realis- tic PBH formation models [1, 8]. PBHs predominantly form binaries through two principal channels: early-universe gravitational de- coupling and late-time dynamical interactions within dark matter halos. The merger rate of PBH bina- ries, defined as the number of mergers per unit vol- ume and unit time, significantly impacts the stochas- tic gravitational-wave background (SGWB). It can be mathematically expressed as: dRE2 d lnm1d lnm2 ≈ 1.6× 106 Gpc−3yr−1 f 53/37 PBH ( M M⊙ )−32/37 × ( t t0 )−34/37 η−34/37S1S2 ψ(m1)ψ(m2) (4) In this formula, M = m1 +m2 represents the total binary mass, η = m1m2/(m1+m2) 2 is the symmetric mass ratio, and t is the cosmic time with t0 being the current age of the Universe. The suppression factors S1 and S2 account for the disruption of PBH bina- ries. The factor S1 models binary disruption caused by matter fluctuations, while S2 captures the disrup- tion due to early PBH clustering from initial Poisson fluctuations when PBHs are abundant [1, 2]. S1 ≈ 1.42 [ ⟨m2⟩/⟨m⟩2 N̄ + C + σ2 M f2PBH ]−21/74 e−N̄ , (5) S2 ≈ min ( 1, 9.6× 10−3f−0.65 PBH e(0.03 ln2 fPBH) ) . (6) Late-time binary formation through halo-based capture provides a secondary contribution and is char- acterized by [2]: dRL2 d lnm1d lnm2 ≈ 3.4× 10−6 Gpc−3yr−1 f2PBH ( σv km/s )−11/7 × η−5/7 δeff ψ(m1)ψ(m2) (7) Here, δeff describes the effective halo density con- trast, and σv is the characteristic velocity dispersion within dark matter halos. Finally, the GW energy density can initially be de- fined as: ΩGW(f) = 1 ρc dρGW d ln f . (8) This equation can now be extended by including the contributions from merging PBH binaries, leading to this full expression [1]: ΩGW(f) = f ρc,0 ∫ zmax 0 dz ∫ d lnm1d lnm2 p(m1)p(m2) (1 + z)H(z) × d2R d lnm1d lnm2 (z,m1,m2) dEGW dfr (9) where ρc,0 is the current critical density of the Uni- verse, H(z) the Hubble parameter, and dEGW/dfr the GW energy emitted by the spectrum per merger event in the source frame. The integration of this equation allows for comparison with empirical constraints from GW detectors such as LVK [8, 10]. III. COMPUTATIONAL METHODS In order to test this theoretical framework against GW observations, a comprehensive computational pipeline was used. This computational framework specifically computes the SGWB arising from PBH binary mergers by combining Monte Carlo sampling methods, semi-analytical integrations, and empirical gravitational wave constraint data from the LVK col- laboration. Initially, the PBH mass function is defined. Then, a log-normal mass distribution parameterized by a cen- tral mass µ and width σ is adopted. These two pa- rameters, along with the overall PBH fraction of dark matter, fPBH, and the contributions from the astro- physical background, Ωcbc, form the parameter space explored through Bayesian inference techniques. The analysis uses the Bayesian inference library bilby to utilize the Markov Chain Monte Carlo (MCMC) algo- rithm provided by the emcee sampler. 3 Log-uniform priors were assigned to all model pa- rameters to reflect the wide uncertainty in their values and to ensure unbiased exploration across several or- ders of magnitude. The prior distributions and their corresponding ranges are summarized in Table I. Parameter Prior Type Range µ [M⊙] LogUniform 10−2 to 120 σ LogUniform 10−2 to 100 fPBH LogUniform 10−5 to 1 Ωcbc LogUniform 10−10 to 10−7 Table I. Log-uniform prior distributions used in the Bayesian inference. A critical component of the computational frame- work is the construction of the likelihood func- tion. This is implemented via a custom PBHModel class, which computes the predicted SGWB spectrum ΩGW(f) and compares it to the 95% upper limits pub- lished by the LVK Collaboration from the O1–O3 ob- serving runs. This comparison is performed through- out the frequency range in the observational band, with a particular focus around ∼ 25Hz where the de- tectors are most sensitive. The likelihood penalizes models that exceed observational bounds, effectively constraining the parameter space. Formally, the log- likelihood is evaluated as logL = −1 2 ∑ i ( ΩGW(fi)− ΩUL(fi) σ(fi) )2 , where ΩUL(fi) denotes the observed upper limit and σ(fi) is an effective uncertainty scale. In the absence of a detected SGWB signal, this likelihood acts as a conservative filter, ruling out models that are incon- sistent with current bounds. To interpret and visualize the statistical results of the Bayesian analysis, log-scaled corner plots are used. These graphs provide a detailed representation of the correlations and constraints of the parameters µ, σ, and fPBH. Furthermore, two-dimensional density es- timates are produced in the µ–fPBH plane as a contour plot to visualize the credible parameter regions and di- rectly compare with the existing observational limits from the LVK data. In addition, a complementary ap- proach is used by fixing the distribution parameters (µ and σ) and systematically calculating the maximum allowed fPBH at each mass scale. The computational framework is implemented in Python and is organized into modules for mass func- tions, merger rates, suppression effects, and likelihood evaluations. A configuration system enables quick switching between modeling scenarios, while input datasets, including LVK constraints and PBH lim- its, are dynamically loaded during execution. This setup enables efficient exploration of PBH parame- ter space and refines dark matter constraints using gravitational-wave observations. IV. RESULTS AND DISCUSSION The analysis constrains PBH populations with a log-normal mass distribution, focusing on the param- eters (µ, σ), fPBH, and the astrophysical background Ωcbc. The log-scaled corner plots obtained from the LVK O1–O3 SGWB data are shown below. Figure 1. Corner plots showing log-scaled posterior distri- butions for the parameters, without suppression (top) and with suppression (bottom). The log-scaled corner plots (Figure 1) illustrate parameter correlations and posterior distributions for (µ, σ, fPBH,Ωcbc) with and without suppression. These plots highlight the constraints placed by cur- rent observational data, clearly showing how includ- ing suppression effects significantly narrows the viable parameter space. A concise summary of the estimated best-fit param- eter values obtained from our Bayesian inference anal- ysis is provided in Table II: 4 Parameter Without Suppression With Suppression µ [M⊙] 0.33164 0.46184 σ 0.10340 0.09556 fPBH 5.84× 10−4 5.50× 10−4 Ωcbc 5.938× 10−10 6.999× 10−10 Table II. Best-fit parameter estimates from Bayesian in- ference analysis for suppressed and unsuppressed binary formation scenarios. The values of the best-fit parameters for the sup- pressed and unsuppressed binary formation scenar- ios are summarized in Table II. The central mass µ is found to be in the sub-solar range having val- ues around 0.33M⊙ for the unsuppressed case and 0.46M⊙ for the suppressed case. These results suggest that, if PBHs contribute to the dark matter content, their masses are likely concentrated around these in- termediate scales. The width σ of the log-normal dis- tribution remains relatively narrow (∼ 0.1), indicat- ing a tight clustering around the central mass rather than a broader mass spread. The inferred dark mat- ter fractions fPBH are small, approximately 5 × 10−4 in both scenarios, reaffirming that only a subdomi- nant PBH contribution is allowed by current SGWB limits. Finally, the values of Ωcbc are of the order 10−10, consistent with expectations for the residual as- trophysical binary background. Overall, these results demonstrate the strong constraining power of SGWB observations on PBH dark matter scenarios. Figure 2 presents the joint posterior distributions of log10(µ) and log10(fPBH), constructed using a smooth kernel density estimate (KDE) method. Credible re- gions corresponding to confidence levels 30%, 68%, 95%, and 99.7% are indicated, with the cyan dashed curve representing the LIGO upper limit. In both suppression and non-suppression scenarios, the poste- riors favor low fPBH values across a range of µ between ∼ 0.05 and 10M⊙. The suppression model leads to a slightly more concentrated distribution, indicating that early binary disruption mechanisms tighten the allowed parameter space. In contrast, the unsup- pressed case exhibits broader credible regions, espe- cially at higher masses. The posterior density stays well below the LIGO limit in both cases, confirming the strength of the constraints from stochastic GW observations. Figure 3 shows the predicted fPBH values required to match the LVK SGWB constraint at 25 Hz, com- pared to the upper limits. The predicted curves for both suppression and non-suppression binary forma- tion scenarios lie several orders of magnitude below the observational bounds across the entire PBH mass range. This confirms that the parameter spaces ex- plored here are consistent with existing SGWB obser- vations. The close overlap between the suppression and non-suppression predictions also indicates that suppression effects have a relatively small impact on the allowed fPBH values in this frequency regime. The results broadly agree with previous stud- ies [2, 7], but demonstrate tighter constraints due to the inclusion of comprehensive suppression modeling. Figure 2. Kernel Density Estimate (KDE) plots illustrat- ing credible regions without suppression (top) and with suppression (bottom). Figure 3. Exclusion curves for PBH dark matter fraction, comparing suppressed and unsuppressed models. The inferred posterior distributions overlap with mass ranges predicted by PBH formation scenarios involv- ing phenomena such as the QCD phase transition [6], although the corresponding dark matter fractions re- main significantly below unity. In conclusion, the analysis strongly suggests that PBHs within the mass range 0.05–10M⊙ cannot con- stitute the entire dark matter under current GW back- ground constraints. However, a subdominant PBH component remains plausible, particularly when envi- ronmental suppression effects are considered. These results lay the foundation for future studies based on data from the LIGO O4 observation run. 5 V. CONCLUSION Bayesian inference results indicate that primordial black holes (PBHs) in the 0.05–10M⊙ mass range can account for only a small fraction of dark matter. Suppression effects significantly reduce the predicted merger rates and tighten the resulting constraints. Al- though broadly consistent with previous studies, the inclusion of suppression mechanisms provides a more refined exclusion of PBH scenarios. This framework supports future extensions to ex- plore alternative PBH mass distributions, spin dy- namics, and formation channels. Sensitivity improve- ments from LIGO O4 and future missions such as LISA and DECIGO are expected to further tighten constraints, solidifying SGWB measurements as a powerful probe of dark matter and early-universe physics. [1] Boybeyi, T., Clesse, S., Kuroyanagi, S., & Sakel- lariadou, M. Search for a gravitational wave background from primordial black hole binaries using data from the first three LIGO-Virgo- KAGRA observing runs. arXiv 2412.18318 (2024). https://arxiv.org/abs/2412.18318 [2] Raidal, M., Vaskonen, V., & Veermäe, H. For- mation of primordial black hole binaries and their merger rates. arXiv 2404.08416 (2024). https://arxiv.org/abs/2404.08416 [3] Escriva, A., Kuhnel, F., & Tada, Y. Pri- mordial Black Holes. arXiv 2211.05767 (2025). https://arxiv.org/abs/2211.05767 [4] Carr, B., & Kuhnel, F. Primordial Black Holes as Dark Matter Candidates. arXiv 2110.02821 (2022). https://arxiv.org/abs/2110.02821 [5] Escriva, A., Germani, C., & Sheth, R. K. Analyt- ical thresholds for black hole formation in general cosmological backgrounds. arXiv 2007.05564 (2020). https://arxiv.org/abs/2007.05564 [6] Inomata, K., Kawasaki, M., Mukaida, K., & Yanagida, T. K. Double inflation as a single ori- gin of primordial black holes for all dark matter and LIGO observations. arXiv 1711.06129 (2019). https://arxiv.org/abs/1711.06129 [7] Clesse, S., & García-Bellido, J. Detecting the gravitational wave background from primordial black hole dark matter. arXiv 1812.11011 (2020). https://arxiv.org/abs/1812.11011 [8] Iovino, A. J. Cosmic Whispers of the Early Uni- verse: Gravitational Waves and Dark Matter from Primordial Black Holes. Ph.D. Thesis, Sapienza University of Rome (2025). arXiv 2501.03065. https://arxiv.org/abs/2501.03065 [9] Ali-Haïmoud, Y., Kovetz, E. D., & Kamionkowski, M. Merger rate of primordial black-hole binaries. arXiv 1709.06576 (2017). https://arxiv.org/abs/1709.06576 [10] Sasaki, M., Suyama, T., Tanaka, T., & Yokoyama, S. Primordial black holes—perspectives in gravita- tional wave astronomy. arXiv 1801.05235 (2018). https://arxiv.org/abs/1801.05235 https://arxiv.org/abs/2412.18318 https://arxiv.org/abs/2404.08416 https://arxiv.org/abs/2211.05767 https://arxiv.org/abs/2110.02821 https://arxiv.org/abs/2007.05564 https://arxiv.org/abs/1711.06129 https://arxiv.org/abs/1812.11011 https://arxiv.org/abs/2501.03065 https://arxiv.org/abs/1709.06576 https://arxiv.org/abs/1801.05235 Constraints on Primordial Black Hole Dark Matter from the Stochastic Gravitational-Wave Background Abstract Introduction Theoretical Framework Computational Methods Results and Discussion Conclusion References