If dark matter DM consists of primordial black holes PBHs and particles simultaneously PBHs are generically embedded within particle DM halos. Such dressed PBHs dPBHs are subject to modified constraints compared to PBHs and can contribute to significant DM abundance in the mass range inlineformulammlmath displayinlinemmlmrowmmlmsupmmlmrowmmlmn10mmlmnmmlmrowmmlmrowmmlmommlmommlmn1mmlmnmmlmrowmmlmsupmmlmimmlmimmlmsupmmlmrowmmlmn10mmlmnmmlmrowmmlmrowmmlmn2mmlmnmmlmrowmmlmsupmmlmsubmmlmrowmmlmiMmmlmimmlmrowmmlmrowmmlmo stretchyfalsemmlmommlmrowmmlmsubmmlmrowmmlmathinlineformula. We show that diffractive lensing of chirping gravitational waves from binary mergers can not only discover but can also identify dPBH lenses and discriminate them from bare PBHs on the eventbyevent basis with potential to definitively establish the coexistence of subdominant PBHs and particle DM.
['Cosmology', 'Astrophysics', 'and Gravitation', 'Astrophysics - Cosmology and Nongalactic Astrophysics', 'Astrophysics - High Energy Astrophysical Phenomena', 'General Relativity and Quantum Cosmology', 'High Energy Physics - Phenomenology']
Coexistence Test of Primordial Black Holes and Particle Dark Matter from Diffractive Lensing
{'Coexistence Test of Primordial Black Holes and Particle Dark Matter from Diffractive Lensing': 'Han Gil Choi, Sunghoon Jung, Philip Lu, Volodymyr Takhistov  \nWe provide additional details of GW lensing, distinguishability of dressed and bare PBHs as well as sensitivity analysis.', '[astro-ph.HE].': "- [65] G. Bertone, A. R. A. C. Wierda, D. Gaggero, B. J. Kavanagh, M. Volonteri, and N. Yoshida, Dark Matter Mounds: towards a realistic description of dark matter overdensities around black holes, (2024), arXiv:2404.08731 [astro-ph.CO].\n- [66] J. Adamek, C. T. Byrnes, M. Gosenca, and S. Hotchkiss, WIMPs and stellar-mass primordial black holes are incompatible, Phys. Rev. D 100 , 023506 (2019), arXiv:1901.08528 [astro-ph.CO].\n- [67] B. C. Lacki and J. F. Beacom, Primordial Black Holes as Dark Matter: Almost All or Almost Nothing, Astrophys. J. Lett. 720 , L67 (2010), arXiv:1003.3466 [astroph.CO].\n- [68] M. P. Hertzberg, S. Nurmi, E. D. Schiappacasse, and T. T. Yanagida, Shining Primordial Black Holes, Phys. Rev. D 103 , 063025 (2021), arXiv:2011.05922 [hep-ph].\n- [69] S. Nurmi, E. D. Schiappacasse, and T. T. Yanagida, Radio signatures from encounters between neutron stars and QCD-axion minihalos around primordial black holes, JCAP 09 , 004, arXiv:2102.05680 [hep-ph].\n- [70] M. Oguri, V. Takhistov, and K. Kohri, Revealing dark matter dress of primordial black holes by cosmological lensing, Phys. Lett. B 847 , 138276 (2023), arXiv:2208.05957 [astro-ph.CO].\n- [71] R.-G. Cai, T. Chen, S.-J. Wang, and X.-Y. Yang, Gravitational microlensing by dressed primordial black holes, JCAP 03 , 043, arXiv:2210.02078 [astro-ph.CO].\n- [72] K. Eda, Y. Itoh, S. Kuroyanagi, and J. Silk, Gravitational waves as a probe of dark matter minispikes, Phys. Rev. D 91 , 044045 (2015), arXiv:1408.3534 [gr-qc].\n- [73] B. J. Kavanagh, D. A. Nichols, G. Bertone, and D. Gaggero, Detecting dark matter around black holes with gravitational waves: Effects of dark-matter dynamics on the gravitational waveform, Phys. Rev. D 102 , 083006 (2020), arXiv:2002.12811 [gr-qc].\n- [74] A. Coogan, G. Bertone, D. Gaggero, B. J. Kavanagh, and D. A. Nichols, Measuring the dark matter environments of black hole binaries with gravitational waves, Phys. Rev. D 105 , 043009 (2022), arXiv:2108.04154 [grqc].\n- [75] P. S. Cole, A. Coogan, B. J. Kavanagh, and G. Bertone, Measuring dark matter spikes around primordial black holes with Einstein Telescope and Cosmic Explorer, Phys. Rev. D 107 , 083006 (2023), arXiv:2207.07576 [astro-ph.CO].\n- [76] P. S. Cole, G. Bertone, A. Coogan, D. Gaggero, T. Karydas, B. J. Kavanagh, T. F. M. Spieksma, and G. M. Tomaselli, Distinguishing environmental effects on binary black hole gravitational waveforms, Nature Astron. 7 , 943 (2023), arXiv:2211.01362 [gr-qc].\n- [77] T. T. Nakamura, Gravitational lensing of gravitational waves from inspiraling binaries by a point mass lens, Phys. Rev. Lett. 80 , 1138 (1998).\n- [78] R. Takahashi and T. Nakamura, Wave effects in the gravitational lensing of gravitational waves from chirping binaries, The Astrophysical Journal 595 , 1039 (2003).\n- [79] S. Jung and C. S. Shin, Gravitational-Wave Fringes at LIGO: Detecting Compact Dark Matter by Gravitational Lensing, Phys. Rev. Lett. 122 , 041103 (2019), arXiv:1712.01396 [astro-ph.CO].\n- [80] K.-H. Lai, O. A. Hannuksela, A. Herrera-Mart'ın, J. M. Diego, T. Broadhurst, and T. G. Li, Discover-\n- g intermediate-mass black hole lenses through gravitational wave lensing, Physical Review D 98 , 083005 (2018).\n- [81] J. Urrutia and V. Vaskonen, Lensing of gravitational waves as a probe of compact dark matter, Mon. Not. Roy. Astron. Soc. 509 , 1358 (2021), arXiv:2109.03213 [astro-ph.CO].\n- [82] S. Basak, A. Ganguly, K. Haris, S. Kapadia, A. K. Mehta, and P. Ajith, Constraints on Compact Dark Matter from Gravitational Wave Microlensing, Astrophys. J. 926 , L28 (2022), arXiv:2109.06456 [gr-qc].\n- [83] R. Abbott et al. (LIGO Scientific, VIRGO), Search for Lensing Signatures in the Gravitational-Wave Observations from the First Half of LIGO-Virgo's Third Observing Run, Astrophys. J. 923 , 14 (2021), arXiv:2105.06384 [gr-qc].\n- [84] R. Abbott et al. (LIGO Scientific, VIRGO, KAGRA), Search for gravitational-lensing signatures in the full third observing run of the LIGO-Virgo network, (2023), arXiv:2304.08393 [gr-qc].\n- [85] H. Zhou, Z. Li, K. Liao, and Z. Huang, Constraints on compact dark matter from lensing of gravitational waves for the third-generation gravitational wave detector, Mon. Not. Roy. Astron. Soc. 518 , 149 (2022), arXiv:2206.13128 [astro-ph.CO].\n- [86] L. Dai, S.-S. Li, B. Zackay, S. Mao, and Y. Lu, Detecting lensing-induced diffraction in astrophysical gravitational waves, Physical Review D 98 , 104029 (2018).\n- [87] M. Oguri and R. Takahashi, Probing Dark Low-mass Halos and Primordial Black Holes with Frequencydependent Gravitational Lensing Dispersions of Gravitational Waves, Astrophys. J. 901 , 58 (2020), arXiv:2007.01936 [astro-ph.CO].\n- [88] H. G. Choi, C. Park, and S. Jung, Small-scale shear: Peeling off diffuse subhalos with gravitational waves, Phys. Rev. D 104 , 063001 (2021), arXiv:2103.08618 [astro-ph.CO].\n- [89] G. Tambalo, M. Zumalac'arregui, L. Dai, and M. H.Y. Cheung, Gravitational wave lensing as a probe of halo properties and dark matter, arXiv preprint arXiv:2212.11960 (2022).\n- [90] V. S. Berezinsky, V. I. Dokuchaev, and Y. N. Eroshenko, Formation and internal structure of superdense dark matter clumps and ultracompact minihaloes, JCAP 11 , 059, arXiv:1308.6742 [astro-ph.CO].\n- [91] M. Boudaud, T. Lacroix, M. Stref, J. Lavalle, and P. Salati, In-depth analysis of the clustering of dark matter particles around primordial black holes. Part I. Density profiles, JCAP 08 , 053, arXiv:2106.07480 [astroph.CO].\n- [92] S. Deguchi and W. D. Watson, Wave effects in gravitational lensing of electromagnetic radiation, Physical Review D 34 , 1708 (1986).\n- [93] T. T. Nakamura and S. Deguchi, Wave optics in gravitational lensing, Progress of Theoretical Physics Supplement 133 , 137 (1999).\n- [94] J.-P. Macquart, Scattering of gravitational radiation: Second order moments of the wave amplitude, Astron. Astrophys. 422 , 761 (2004), arXiv:astro-ph/0402661.\n- [95] S. Jung and S. Kim, Solar diffraction of LIGO-band gravitational waves, JCAP 07 , 042, arXiv:2210.02649 [astro-ph.CO].\n- [96] R. Takahashi, Quasi-geometrical optics approximation in gravitational lensing, Astronomy & Astrophysics  \n423 , 787 (2004).  \n- [97] S. Savastano, G. Tambalo, H. Villarrubia-Rojo, and M. Zumalacarregui, Weakly lensed gravitational waves: Probing cosmic structures with wave-optics features, arXiv preprint arXiv:2306.05282 (2023).\n- [98] C. Cutler and E. E. Flanagan, Gravitational waves from merging compact binaries: How accurately can one extract the binary's parameters from the inspiral waveform?, Physical Review D 49 , 2658 (1994).\n- [99] B. P. Abbott et al. (LIGO Scientific, Virgo), A guide to LIGO-Virgo detector noise and extraction of transient gravitational-wave signals, Class. Quant. Grav. 37 , 055002 (2020), arXiv:1908.11170 [gr-qc].\n- [100] J. D. Romano and N. J. Cornish, Detection methods for stochastic gravitational-wave backgrounds: a unified treatment, Living reviews in relativity 20 , 1 (2017).\n- [101] C. Talbot and E. Thrane, Measuring the binary black hole mass spectrum with an astrophysically motivated parameterization, The Astrophysical Journal 856 , 173 (2018).\n- [102] B. Abbott, R. Abbott, T. Abbott, S. Abraham, F. Acernese, K. Ackley, C. Adams, R. X. Adhikari, V. Adya, C. Affeldt, et al. , Binary black hole population properties inferred from the first and second observing runs of advanced ligo and advanced virgo, The Astrophysical Journal Letters 882 , L24 (2019).\n- [103] R. Abbott, T. Abbott, S. Abraham, F. Acernese, K. Ackley, A. Adams, C. Adams, R. Adhikari, V. Adya, C. Affeldt, et al. , Population properties of compact objects from the second ligo-virgo gravitational-wave transient catalog, The Astrophysical journal letters 913 , L7 (2021).\n- [104] R. Abbott et al. (KAGRA, VIRGO, LIGO Scientific), Population of Merging Compact Binaries Inferred Using Gravitational Waves through GWTC-3, Phys. Rev. X 13 , 011048 (2023), arXiv:2111.03634 [astro-ph.HE].\n- [105] B. P. Abbott, R. Abbott, T. Abbott, S. Abraham, F. Acernese, K. Ackley, C. Adams, V. Adya, C. Affeldt, M. Agathos, et al. , Prospects for observing and localizing gravitational-wave transients with advanced ligo, advanced virgo and kagra, Living reviews in relativity 23 , 1 (2020).\n- [106] S. Hild, M. Abernathy, F. e. Acernese, P. AmaroSeoane, N. Andersson, K. Arun, F. Barone, B. Barr, M. Barsuglia, M. Beker, et al. , Sensitivity studies for third-generation gravitational wave observatories, Classical and Quantum gravity 28 , 094013 (2011).\n- [107] V. Srivastava, D. Davis, K. Kuns, P. Landry, S. Ballmer, M. Evans, E. D. Hall, J. Read, and B. Sathyaprakash, Science-driven tunable design of cosmic explorer detectors, The Astrophysical Journal 931 , 22 (2022).\n- [108] J. Urrutia, V. Vaskonen, and H. Veermae, Gravitational wave microlensing by dressed primordial black holes, Phys. Rev. D 108 , 023507 (2023), arXiv:2303.17601 [astro-ph.CO].", 'A. Fitting geometric optics regime': 'In the geometric optics limit ( f →∞ ) of GW lensing, the solution of the Fresnel-Kirchhoff integral can be approximated by the usual two images [78, 93]  \nF ( f ) /similarequal | µ 1 | 1 / 2 + | µ 2 | 1 / 2 e 2 πif ∆ t -i π 2 , (S1)  \nwhere µ 1 and µ 2 are the magnifications of the images, and ∆ t is the arrival time difference between the lensing images. In GW observation of binary BH sources, the constant amplification due to lensing and the luminosity distance to the binary BH source are degenerate unless the source redshift is obtained from other observations. Therefore, the lensing amplification is reduced to  \nF ( f ) ∝ 1 + µ 1 / 2 r e 2 πif ∆ t -i π 2 , (S2)  \nand the only observables are the image magnification ratio µ r = | µ 1 /µ 2 | and ∆ t .  \nA point mass (PBH) lens always produces two lensing images. From its image properties, one can find that  \nµ r = 1 -( √ y 2 s / 4 + 1 + y s / 2 ) -1 1 -( √ y 2 s / 4 + 1 -y s / 2 ) -1 , (S3)  \nand  \n∆ t = 4 M l [ y s √ y 2 s / 4 + 1 + ln ( 1 + y 2 s / 2 + y s √ y 2 s / 4 + 1 ) ] , (S4)  \nwhere y s = x s /x E , and M l = M PBH (1 + z l ). Note that relationship between the pairs of numbers ( µ r , ∆ t ) and ( y s , M l ) is invertible, and one can find the inverse relation  \ny s = √ µ 1 / 2 r + µ -1 / 2 r -2 M l = ∆ t 2 ( √ µ r + µ -1 r -2 -ln µ r ) . (S5)  \nThis simple relation implies that any lensing signal induced by the interference of two lensing images can be interpreted as point mass lensing. Hence, in this situation, there is no distinguishability and different lens profiles are degenerate, unless additional information is given.', 'B. Strong diffraction versus weak diffraction': 'In our work, the weak diffraction (WD) and beyond Geometric Optics (bGO) signatures enable us to discriminate dPBHs from bare PBHs and others. Non-locally determined by lens profile, they induce nontrivial frequencydependent behaviors.  \nStrong diffraction is also a non-local phenomenon, which occurs for the case with x s < x E . But this effect does not contribute to distinguishability since it is related to the blurring of the Einstein ring, where any symmetric strong  \nFIG. S1. Same as Fig. 2, but for the case with x s < x E , which induces strong diffraction for the intermediate regime √ x s x E ≤ x F ≤ x E . This regime does not facilitate discrimination, once the GO regime is well fitted, but it still helps lensing detection due to the nontrivial frequency dependence over the unlensed F/A 0 = 1.  \n<!-- image -->  \nlensing system has similar properties [88]. This is illustrated in Fig. S1, which is produced with the same parameters as in Fig. 2 of the main text but in the case of x s < x E . Instead of bGO, strong diffraction occurs in the frequency range satisfying √ x s x E < x F < x E . In this regime, | F | ∝ f 1 / 2 for any symmetric lens profiles [88]. As a result, there is no distinction between dPBH lensing and PBH lensing (dashed) in this regime once data are matched at the geometric optics limit.', 'C. Observed events with discrimination capability': 'As a proxy of discrimination efficiency, we consider the fraction of observed GW lensing events that contain enough weak diffraction and bGO regimes in measured bands. Our analysis serves as a pilot study for a more dedicated analyses in the future. More specifically, we estimate the discrimination fraction of detectable lensing events by the following two conditions  \nx max F > x E , (S6a)  \nρ WD+bGO ≡ √ 4 ∫ f e f min df | h 0 ( f ) | 2 S n ( f ) > 10 , (S6b)  \nwhere x max F is the Fresnel length at the lowest sensitive frequency of the GW detector, and f e is the GW frequency at the x F = x E instance. x E is the Einstein radius of dPBH.  \nThe former condition Eq. (S6a) approximately corresponds to requiring that WD and bGO are included in the measurement band. Since the bGO effect outside the frequency range defined by Eq. (S6a) is small but non-zero, it could provide some discrimination power; Eq. (S6a) is a conservative requirement that can be further improved in dedicated analyses. The latter condition Eq. (S6b) approximately requires that the GW amplitude resolution in the weak diffraction plus bGO regime be better than ∼ 10% level. This is approximately the scale of differences observed in Fig. 2, which we also expect to be improved with a dedicated study.  \nWithin the region of the detectable lensing cross-section, we obtain the area of the region satisfying Eq. (S6a) and (S6b), which defines the discrimination cross-section, similarly to the lensing cross-section. Subsequently, we compute associated discrimination optical depths and discrimination rates. The discrimination rate divided by the lensing detection rate is the discrimination fraction of detectable lensing events.  \nFig. 4 of the main text shows the the discrimination fraction of detectable lensing events for Cosmic Explorer (CE, solid), Einstein Telescope(ET, dashed), aLIGO A+(dotted), and aLIGO(dot-dashed). We find that the discrimination fraction sharply drops near M PBH ∼ O (10 -10 2 ) M /circledot , where the Einstein radius of dPBH becomes comparable to the  \nFIG. S2. (Left panel) chirping GW spectrums assuming dPBH lensing (blue solid), bare PBH lensing (orange dashed), and no lensing (black solid). The characteristic strain h c ( f ) ≡ fh ( f ) is used. The GW source is 10 M /circledot -10 M /circledot binary BH at z s = 0 . 4. The lensing amplifications are given by M PBH = 20 M /circledot and z l = 0 . 17; The lensing parameters are the same as those in Fig. 2 of the main text. (Right panel) GW strain difference ∆ h c ≡ | h 1 c -h 2 c | between the curves given in the left panel. The difference between the dPBH and PBH (green solid) is a signal of discrimination, while the difference between the dPBH and no lensing (red dashed) is a signal of lensing detection. We also present the noise amplitude √ fS n ( f ) of CE (black dot-dashed) and ET (black dotted) for comparison.  \n<!-- image -->  \nFresnel lengths of GW spectrum. Above this mass range, lensing signals containing only geometric optics signatures are more frequent due to the large Einstein radius.  \nFigure S2 shows a typical example of the lensed GW strain that satisfies our discrimination criteria. We assume a chirping GW originated from a 10 M /circledot -10 M /circledot BH binary at z s = 0 . 4. The noticeable discrepancy between dPBH and PBH in the low-frequency region is due to diffractive lensing. Its magnitude is far above the noise amplitudes of CE and ET, hence it is measurable in those detectors. In the high-frequency region where lensing is described by GO, the discrepancy is buried under the detector noise and does not contribute to discrimination. Although the GO regime does not have discrimination power, it is a useful feature for lensing detection. Fig. S2 shows that the difference between dPBH lensing and no lensing is comparable at both high and low frequencies, equally contributing to lensing detection.', 'A. Likelihood ratio': 'Our definition of lensing detection significance is based on the likelihood ratio  \nΛ ≡ p ( d |H L ) p ( d |H 0 ) (S7)  \nbetween the lensing detection hypothesis H L and the no lensing hypothesis H 0 given the strain data d . The probability p ( d |H 0 ) and p ( d |H 1 ) are the marginal likelihoods [100]. The quantity Eq. (S7) is also known as the Bayes factor. In principle, the marginal likelihoods should be computed by the integral of posterior distribution times prior distribution, p ( d |H ) = ∫ d θ p ( θ | d, H ) p ( θ , H ), where θ is the model parameters given a hypothesis H . However, the computation of the integral for many lensing data realizations requires prohibitively large computing resources. Instead, we use an approximate version of Eq. (S7) [86, 100],  \nln Λ /similarequal 1 2 ( ρ 2 mL -ρ 2 m 0 ) , (S8)  \nwhere  \nρ 2 m ≡ max θ [( d | d ) -( d -h ( θ ) | d -h ( θ ))] (S9)  \nFIG. S3. The false-alarm probability for lensing detection criteria 〈 ln Λ 〉 > ln Λ c . In the cases of aLIGO (solid curves), 20 M /circledot binary BH with z s = 0 . 1 is assumed, and we assume the dressed PBH lensing signals are falsely identified with M PBH = 50 M /circledot (blue) and M PBH = 100 M /circledot (orange). z l is set to give the maximum d eff and x s is set to give 〈 ln Λ 〉 = ln Λ c . In the case of CE (dashed curves), only source redshift is changed to z s = 1.  \n<!-- image -->  \nis the square of the matched filter Signal-to-Noise Ratio(SNR) with respect to waveform h . Here, the inner product ( ·|· ) is defined as ( a | b ) ≡ 4Re ∫ ∞ 0 df a ∗ ( f ) b ( f ) /S n ( f ) for given strain data a and b under stationary Gaussian detector noise with the noise spectral density S n ( f ). We denote ρ m by ρ mL and ρ m 0 for a lensed h L ( θ L ) and unlensed h 0 ( θ 0 ) waveform, respectively. In these expressions, θ L and θ 0 are model parameters of H L and H 0 , respectively.  \nRelying on the detector noise, the matched filter SNR and the likelihood ratio are stochastic quantities. Assuming the lensed waveform signal h L and the detector noise is small d ∼ h L , the expected value of the log-likelihood ratio is  \n〈 ln Λ 〉 /similarequal 1 2 ( ρ 2 L -ρ 2 uL ) = 1 2 min θ 0 ( h L -h 0 ( θ 0 ) | h L -h 0 ( θ 0 )) , (S10)  \nwhere ρ 2 L ≡ 〈 ρ 2 mL 〉 /similarequal ( h L | h L ), and ρ 2 uL ≡ 〈 ρ 2 m 0 〉 /similarequal max θ 0 ( h 0 ( θ 0 ) | 2 h L -h 0 ( θ 0 )). We use Eq. (S10) to define the lensing cross-section.  \nIn our work, we assume that a lensing signal h L which gives 〈 ln Λ 〉 > ln Λ c is detectable. We should choose a large enough threshold ln Λ c to avoid the situations when ln Λ > ln Λ c is satisfied even if the strain data does not contain any lensing signal. We use the False-alarm probability P (ln Λ > ln Λ c | nolensing) to justify our choice of ln Λ c = 3. To compute the probability, we generate 10 5 realizations of d = n + h 0 with a fixed h 0 and the stationary Gaussian detector noise n and count the number of ln Λ > ln Λ c samples. In the computation of the ρ 2 mL , we choose h L ( θ L ) that leads to ( ρ 2 L -ρ 2 uL ) / 2 = ln Λ c . To minimize the computation resources for the maximizations, we use the analytic method described in the next subsection. The False-alarm probability results for some combinations of lensing parameters are shown in Fig. II A. We find that the False-alarm probabilities are not sensitive to lensing parameters or GW detectors and drop below 0 . 01 for ln Λ c = 3 which is small enough for the discussions in our work.', 'B. Optimization of the likelihood ratio': 'We consider only three model parameters: constant amplitude factor A , constant phase shift θ c , and time shift t c . We compute the matched filter SNR Eq. (S9) by maximizing the following function:  \n( d | d ) -( d -Ae i ( θ c + wt c ) h | d -Ae i ( θ c + wt c ) h ) (S11)  \nfor a given waveform h . Here we used w = 2 πf . The maximization of Eq. (S11) with respect to A and θ c is given by  \n1 ( h | h ) [( d | e iwt c h ) 2 +( d | ie iwt c h ) 2 ] . (S12)  \nMaximization over t c requires numerical computation in general. However, when the expression is maximized at a small t c , and we always inject the signal with t c = 0, one can obtain a leading order analytic expression for the t c  \nmaximization. The resulting analytic expression of ρ 2 m is  \nρ 2 m /similarequal 1 ( h | h ) [ ( d | h ) 2 +( d | ih ) 2 + { ( d | h )( d | iwh ) -( d | ih )( d | wh ) } 2 ( d | h )( d | w 2 h ) -( d | wh ) 2 +( d | ih )( d | iw 2 h ) -( d | iwh ) 2 ] . (S13)  \nFor a more detailed analysis, one can additionally consider binary BH intrinsic parameters such as total mass, mass ratio, black hole spins, etc, and lens intrinsic parameters such as lens mass, impact parameter, lens redshift, etc. For simplicity, we assume the parameters are not highly biased by the detector noise and are not correlated with each other. We leave a more detailed dedicated analysis for future works.', 'C. Unlensed GW waveform model': 'We use unlensed GW waveform in the Newtonian order to simplify our analysis. The waveform is given by [98]  \nh 0 ( f ) = -√ 5 24 π -2 / 3 A p ( M z ) 5 / 6 d L f -7 / 6 e i Ψ( f ) Ψ( f ) = θ c +2 πft c -π 4 + 3 128 ( π M z f ) -5 / 3 (S14)  \nwhere we used G = c = 1 units. Here, d L is the luminosity distance, and M z = ( m 1 m 2 ) 3 / 5 / ( m 1 + m 2 ) 1 / 5 (1 + z s ) is the detector-frame binary BH chirp mass defined by the component black hole masses m 1 and m 2 . The polarization angle and detector orientation dependence are included in the A p factor, and we set A p = 1 for simplicity.  \nThe upper cut-off of the GW frequency is set to the GW frequency at the inner-most-stable-circular-orbit frequency f isco = (6 3 / 2 π ( m 1 + m 2 )(1 + z s )) -1 . The lower cut-off of the frequency follows the lower frequency bound of the GW detector sensitivity range. We adopt 5 Hz for Cosmic Explorer, 2 Hz for Einstein telescope, and 10 Hz for LIGO.', 'III. LIMIT COMPUTATION': "We assume the occurrence of lensing events follows a Poisson distribution with the expected number of lensing detections N L = ∫ t max 0 dt ˙ N L , taking t max = 5 years. Under the lensing detection hypothesis with f PBH parameter, the probability of k lensing detection is  \nP ( k | f PBH ) = ( f PBH ν L ) k k ! e -f PBH ν L , (S15)  \nwhere ν L = N L ( f PBH = 1). Using Bayes' theorem, the posterior distribution of f PBH is given by  \np ( f PBH | k ) = π ( f PBH ) P ( k | f PBH ) P ( k ) , (S16)  \nwhere π ( f PBH ) is the prior distribution of f PBH , and P ( k ) ≡ ∫ df PBH π ( f PBH ) P ( k | f PBH ) is the marginalized likelihood. The lensing detection rate ˙ N L is given by  \n˙ N L = ∫ z c 0 dV c ( z s ) R 0 1 + z s 〈 τ ( z s ) 〉 , (S17)  \nwhere dV c ( z s ) is the comoving volume element at z s , and  \n〈 τ ( z s ) 〉 = ∫ ρ 0 ( M,η,z s ) > 8 dMdηp BBH ( M,η ) τ ( z s ; M,η ) (S18)  \nis the expected lensing optical depth. Here ρ 0 ( M,η,z s ) ≡ √ ( h 0 | h 0 ) is the optimal SNR [99] of unlensed GW signal h 0 given by binary total mass M = m 1 + m 2 , mass ratio η = ( m 1 m 2 ) /M 2 , and redshift z s .  \nThe optical depth τ is approximately the detectable lensing probability for the given source, 1 -e -τ /similarequal τ , and it is obtained by  \nτ /similarequal ∫ z s 0 dχ ( z l ) n PBH σ ( z l , z s ; M,η ) , (S19)  \nwith the comoving lensing cross-section,  \nσ = (1 + z l ) 2 ∫ 〈 ln Λ 〉 > ln Λ c d 2 x s , (S20)  \nmeasuring the comoving area of lens locations that can lead to detectable lensing. We used ln Λ c = 3 for the detection threshold. The comoving number density of dPBH lenses, n PBH = f PBH Ω DM ρ c /M PBH , is assumed to be constant. ρ c is the critical density of the Friedmann universe.  \nOur choice of p BBH ( M,η ) follows the Power law and Peak model [101], and we used R 0 = 28 . 3 Gpc -3 yr -1 . We take into account binary BH merger sources up to redshift z c = 10 and assume their populations do not evolve with redshift. We assume the flat ΛCDM model with the Hubble parameter H 0 = 67 . 74 km s -1 Mpc -1 , the matter density parameter Ω m = 0 . 3075, and the dark matter density parameter Ω DM = 0 . 2575.  \nWe consider the null detection k = 0 of lensing events to estimate the projected upper limit sensitivity of f PBH . We set a uniform prior on f PBH with the range from 0 to 1. In this case, P (0) = (1 -e -ν L ) /ν L , and  \np ( f PBH | 0) = ν L e -f PBH ν L (1 -e -ν L ) . (S21)  \nUsing the posterior, the 90% upper limit of f PBH is computed as  \nf 90% PBH = -1 ν L ln [ 1 -0 . 9(1 -e -ν L ) ] . (S22)"}