REVIEW 4 major objections 7 minor
On the orbital eccentricities of primordial black hole binaries inside and outside of dark matter halos
T0 review · 4 major / 7 minor · reviewed 2026-07-30 · grok-4.5
Pith's one-line read LISA and DECIGO together can probe about a hundred stellar-mass primordial black hole binaries with eccentricity above 0.01; excluding them would tighten PBH abundance limits by an order of magnitude.
desk verdict Solid multi-band eccentricity maps for PBH binaries; the O(10^2) e>0.01 forecast and 10x abundance claim ride on inherited rates and a monotonic K-term that may overstate the high-e tail. read the letter →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
What carries the argument
Full orbital evolution of semi-major axis and eccentricity under the coupled equations that include both Peters–Mathews gravitational-wave emission and environmental binary-single terms (density and velocity dispersion inside growing dark-matter halo shells), followed by peak-harmonic characteristic-strain and SNR calculations that map each track onto detector bands and yield rescaled eccentricity histograms.
What would settle it
A multi-year LISA plus DECIGO search that finds zero stellar-mass black-hole binaries with measured eccentricity above 0.01, or that finds a substantially different number than the O(100) predicted at the current abundance ceiling, would directly test the claim.
Extended reading notes
Core claim
Given present limits on stellar-mass PBH abundance, LISA and DECIGO together are expected to detect O(10^2) PBH binaries that still have orbital eccentricity e > 0.01 when their waves enter the detector band; if those observatories can exclude such eccentric systems, the upper limit on the PBH abundance can be strengthened by an order of magnitude. Isolated early binaries circularize completely except for residual O(10^{-2}) eccentricities in LISA, while binaries that experience binary-single interactions inside dense dark-matter halos retain higher eccentricity even at late inspiral.
Load-bearing premise
The predicted detection counts rest on rescaling the simulated eccentricity histograms with the authors’ own earlier merger-rate calculations, normalized to a PBH abundance that saturates current LIGO-Virgo-KAGRA limits; if those rates or that abundance are too high, the O(100) forecast shrinks in proportion.
Editorial extensions
If this is right
- Ground-based detectors (aLIGO, ET, CE) will see essentially circular PBH binaries, so eccentricity cannot separate them from ordinary astrophysical black-hole binaries in those bands.
- LISA will retain residual eccentricities of order 0.01 for isolated PBH binaries, giving a low-frequency eccentricity window even without halo interactions.
- Dense inner halo shells can leave a few systems with e > 0.1 still visible to space-based detectors.
- A clean null result on e > 0.01 binaries from LISA and DECIGO would tighten the stellar-mass PBH dark-matter fraction by roughly a factor of ten.
- Direct-capture binaries contribute negligibly to the eccentric sample once current abundance limits are imposed.
Reading between the lines
- Eccentricity at millihertz-to-decihertz frequencies becomes a practical abundance diagnostic rather than only a formation-channel tag, linking space-based GW catalogs directly to dark-matter fraction limits.
- The same halo-environment machinery implies that any future tightening of LVK merger-rate bounds will linearly rescale the expected eccentric counts for LISA and DECIGO.
- Multi-band detections that catch the same binary first in LISA/DECIGO (still mildly eccentric) and later in ET/CE (circular) would furnish an independent check of the circularization timescale assumed here.
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. The manuscript computes the eccentricity distributions of stellar-mass primordial black hole (PBH) binaries at the frequencies of current and future GW detectors (LISA, DECIGO, ET, CE, aLIGO), tracking three channels: early-formed binaries evolving in isolation, early-formed binaries that enter dark matter halos and undergo binary-single interactions, and direct GW captures in halos. The authors integrate coupled semimajor-axis/eccentricity evolution equations (GW Peters–Mathews terms plus environmental hardening/eccentricity-pump terms) for samples of ~5×10^6 binaries per halo shell, compute characteristic strains via the peak harmonic, and rescale simulated histograms with channel merger rates from Ref. [36] at f_PBH = 3.4×10^-3. They find that unperturbed binaries circularize below detectability in all bands except LISA (residual e ~ 10^-2), that halo binaries retain a high-eccentricity tail, and that LISA+DECIGO together should probe O(10^2) binaries with e > 0.01, so that an exclusion of such eccentric binaries could improve PBH abundance limits by an order of magnitude.
Significance. If the results hold, this is a useful and timely forecast: it is, to my knowledge, the first study to follow PBH binary eccentricity self-consistently across all formation channels and all present/planned GW bands, with large Monte Carlo samples (5×10^6 binaries per shell/halo), explicit shell-resolved halo environments, and per-detector SNR/horizon calculations. The output is a falsifiable prediction — a concrete count of eccentric (e > 0.01) events for LISA and DECIGO under a stated abundance normalization — and a clear statement of what an exclusion would imply for f_PBH. The unperturbed-channel result (full circularization before ground-based bands, residual e ~ 10^-2 in LISA) is a robust, parameter-light check that agrees with prior analytic expectations and lends credibility to the pipeline. The main weaknesses are not in the GW-emission machinery (Eqs. 3–10 are standard and correctly applied) but in the imported environmental-interaction physics and in the normalization of the headline numbers.
major comments (4)
- [§IIA, Eqs. (1)–(2)] The headline O(10^2) count of e>0.01 detections is supplied almost entirely by the high-e tail of the binary-single channel (Fig. 6: DECIGO Ndet = 11840 with a tail extending to e ~ 1, versus the unperturbed channel which is below 0.01 in the DECIGO band). That tail is generated by the environmental term de/dt = + G H K(r,t) ρ_env a / v_disp, with H and K taken from Quinlan (1996) and Sesana et al. (2006) — N-body calibrations for massive black hole binaries hardening in dense stellar cusps, i.e., an extreme mass-ratio regime with a large reservoir of light perturbers. Three compounded issues: (i) as written, the K term is strictly positive, so while the environmental term dominates, eccentricity can only be pumped upward, in some shells for ~Gyr; real equal-mass binary-single encounters change e stochastically in both directions (Heggie-style scattering), and a monotonic secular pump is
- [§IIC, Eq. (12)] The SNR integral is taken over the overlap of the binary's full detector-frame frequency range [f_bin,min, f_bin,max] with the detector band, and t_obs enters only as a linear multiplier in Eq. (15). For a 30+30 M_sun binary, the inspiral time from 10^-3 Hz to merger is ~10^4 yr, and from 10^-2 Hz is ~20 yr — both far longer than the assumed t_obs = 5 yr (LISA) and 3 yr (DECIGO). If f_bin,max is the ISCO frequency rather than the frequency reached after t_obs of observation, the SNR, and hence z_max and N_det, are overestimated for the space detectors. The text does not state what f_bin,max is. Please clarify, and if the finite observation time is not accounted for, recompute z_max and N_det with f_start set by t_obs before merger (as in standard multiband treatments, e.g., Ref. [75]). This is load-bearing because the LISA numbers (N_det = 8 and 68 in the two channels) sit close to thres
- [§IIB, Eqs. (9)–(11)] The characteristic strain is approximated by the single peak harmonic, dropping the quadrature sum of Eq. (11). At moderate eccentricities the power is spread over several harmonics around n_peak, so this approximation misestimates h_c by a factor that varies with e; since z_max scales with h_c and N_det roughly with z_max^3 at low z, even a ~20–30% strain error is non-negligible for the event counts. The error could be quantified cheaply on a subsample by comparing Eq. (11) with the peak-harmonic approximation. Relatedly, the orientation-averaging factor F in Eq. (12) is introduced but its numerical value is never given; please state it.
- [§IIIB and §IV, Eq. (15) and Fig. 7] The two headline numbers — O(10^2) binaries with e>0.01 and the 'order of magnitude' improvement in f_PBH limits — are linear rescalings of the channel merger rates R_channel(z) imported from Ref. [36] (a companion preprint by the same authors) at the single normalization point f_PBH = 3.4×10^-3, f_PBH binaries = 0.5, m = 30 M_sun. This is not circular in a logical sense, but the manuscript should be more explicit that both numbers scale linearly with the assumed rate normalization, and should tabulate the central result directly: N_det(e>0.01) and N_det(e>0.1) per detector per channel. At present the O(10^2) figure can only be extracted by eye from the Fig. 7 histogram, which is unsatisfactory for the paper's main quantitative claim. In addition, DECIGO's N_det is dominated by the volume out to the (capped) z_max = 1000, where the halo mass function and the Ref. [36] rates are extrapola
minor comments (7)
- [Fig. 4 captions/legends] The legends quote 'DECIGO (10^-2 Hz, z ≤ 12)' while Table I gives z_max = 1000 for DECIGO. Presumably the z ≤ 12 reflects that halo-channel mergers only occur after halo formation; please make the legends consistent with Table I or explain the cut explicitly.
- [§IIIB, paragraph after Fig. 4] The sentence 'The fact that N1/N8 is much larger than one for DECIGO, ET, and CE, but less than one, is due to...' is incomplete — 'but less than one [for LISA and aLIGO]' appears to be intended.
- [Figs. 8 and 9] Each panel carries a duplicated halo-mass label (e.g., 'M(z=0) = 1.68×10^4 M⊙ M(z=0) = 1.68×10^4 M⊙').
- [§IIIB] Typo: 'withihn' → 'within' (paragraph above Eq. 16 reference to Fig. 5 bottom panel).
- [§IIA] The switch criterion |ȧ|env/|ȧ|GW ≤ 0.5 (and the same for ė) means the environmental term is still a ~50% correction at the handoff time t_GW; a brief note on the sensitivity of the final (a_GW, e_GW) to the choice 0.5 would be useful.
- [§IIC] The noise curves are taken from GWplotter [74]; please specify which curve configuration is used for each instrument (e.g., LISA mission duration/arm length assumption, DECIGO vs B-DECIGO, ET-D, CE 40 km), since z_max values in Table I depend on these choices.
- [§IIIA] The statement that alternative (f_PBH, f_PBH binaries) combinations consistent with LVK give 'the same number of detected events as in Table I' would benefit from one sentence of explanation, since the unperturbed rate generally scales as a non-trivial power of f_PBH.
Simulated Author's Rebuttal
We thank the referee for a careful and constructive report. We agree that the main weaknesses lie in the imported environmental-interaction physics and in the presentation/normalization of the headline numbers, rather than in the GW machinery. We will make four substantive revisions: (i) clarify the role and limitations of the Quinlan/Sesana hardening and eccentricity-pump coefficients and quantify the sensitivity of the high-e tail to the sign stochasticity of K; (ii) define f_bin,max explicitly and restrict the SNR integral to the portion of the inspiral observable within t_obs for LISA and DECIGO, recomputing z_max and N_det; (iii) quantify the error of the peak-harmonic approximation on a subsample and state the numerical value of F; and (iv) tabulate N_det(e>0.01) and N_det(e>0.1) per detector per channel, with explicit statements of the linear scaling with the rate normalization and the extrapolation involved at high z. We address each comment below.
read point-by-point responses
-
Referee: [§IIA] The high-e tail driving the O(10^2) count comes from the environmental term de/dt = + G H K(r,t) ρ_env a / v_disp with H, K from Quinlan (1996)/Sesana et al. (2006), calibrated on massive BH binaries in stellar cusps; K is strictly positive, so eccentricity is pumped monotonically for ~Gyr, whereas real equal-mass binary-single encounters change e stochastically in both directions.
Authors: We agree this is the weakest part of the modeling, and we thank the referee for pressing on it. Two defenses and one concession. First, the use of continuum hardening/eccentricity coefficients as effective orbit-averaged rates is a standard shortcut (e.g. Refs. [62–64, 73]), and our calibration in Ref. [36] — where the same formalism was benchmarked against direct N-body-style encounter rates for PBH halos — reproduces the binary-single merger-rate enhancement reasonably well. Second, the mass regime is less extreme than for MBH binaries in the sense that all perturbers have the same mass as the binary members, so the cusp-reservoir issue does not arise; the H and K coefficients enter only as an effective encounter-rate efficiency. However, we concede the referee's central physical point: our K is a mean positive drift, while Heggie-style equal-mass scattering randomizes e with a thermally-biased distribution, and a strictly monotonic pump over Gyr likely overproduces the extreme tail (e → 1). In the revision we will (a) state this limitation explicitly in §IIA, (b) perform a sensitivity test in which the environmental ė term is treated as a stochastic kick with zero mean drift but the same variance (or capped at the thermal-distribution expectation), and report how N_det(e>0.01) for DECIGO changes, and (c) if the O(10^2) headline count shifts materially, we will revise it and the abstract accordingly. We cannot, within this work, replace the effective equations with full few-body scattering experiments; we will say so plainly. revision: partial
-
Referee: [§IIC, Eq. (12)] f_bin,max is not defined. If it is the ISCO frequency rather than the frequency reached after t_obs, then for space detectors (inspiral times of 10^4 yr from 1e-3 Hz) the SNR, z_max, and N_det are overestimated. Please clarify and recompute with f_start set by t_obs before merger, as in Ref. [75]. This is load-bearing for the LISA counts.
Authors: The referee is correct that the text is ambiguous, and we apologize for the omission. In the submitted version f_bin,max is set by the frequency at ISCO; the only account of t_obs is the linear multiplier in Eq. (15). This indeed overestimates the accumulated SNR for LISA and DECIGO, whose sources dwell in band far longer than the mission duration. In the revision we will (i) define f_bin,min and f_bin,max explicitly, and (ii) restrict the integration in Eq. (12) for the space detectors to the window [f(t_ISCO − t_obs), f_ISCO] (equivalently the harmonics radiated during the observed segment), following the standard multiband treatment of Ref. [75], and recompute z_max and N_det for LISA and DECIGO. For ground-based detectors the inspiral through the band is much shorter than t_obs, so those numbers are unaffected. We note that the affected LISA counts (N_det = 8 and 68) will decrease, and the DECIGO counts will as well; we will update Tables I–II, Figs. 3, 6, 7, the abstract, and the O(10^2) claim accordingly. Since the high-e tail of the binary-single channel resides at the lower-frequency (earlier) portion of the band, the eccentric fraction within t_obs should be less suppressed than the total, but we will report the recomputed numbers rather than assert this. revision: yes
-
Referee: [§IIB, Eqs. (9)–(11)] The characteristic strain uses only the peak harmonic, dropping the quadrature sum of Eq. (11); at moderate e the power is spread over several harmonics, and since z_max ∝ h_c and N_det ∝ z_max^3, even 20–30% strain errors matter. Quantify the error on a subsample. Also, the value of the orientation-averaging factor F in Eq. (12) is never given.
Authors: We agree with both points. On the first: for the eccentricities relevant to the detectable populations (e ≲ 0.1 in the LISA/DECIGO bands for nearly all sources, as the high-e tail merges quickly after entering band), g(n,e) is strongly concentrated at n_peak, so we expect the quadrature correction to be modest; but 'we expect' is not a quantification. We will compute both Eq. (11) and the peak-harmonic approximation for a representative subsample spanning the simulated (a, e) grid and report the fractional difference in h_c, z_max, and N_det as a function of e; if the error exceeds ~10% anywhere in the population that dominates N_det, we will switch to the full sum, which is computationally cheap at our sample sizes. On the second: the value of F was omitted in error. We use F = 1/5 (sky- and orientation-averaged response for a single interferometer, following Refs. [68, 69]); we will state this explicitly in §IIC and note the rescaling of SNR had a different convention been adopted. revision: yes
-
Referee: [§IIIB, §IV] Both headline numbers are linear rescalings of channel merger rates from Ref. [36] at one normalization point (f_PBH = 3.4e-3, f_PBH binaries = 0.5, 30+30 M_sun). Be explicit about the linear scaling; tabulate N_det(e>0.01) and N_det(e>0.1) per detector per channel; and note that DECIGO's counts are dominated by z up to the capped z_max = 1000, where the halo mass function and rates are extrapolated.
Authors: We agree on all three requests. (i) We will add an explicit statement that N_det and the resulting f_PBH-improvement factor scale linearly with the assumed channel rate normalization, and give the simple rescaling formula so readers can translate to other (f_PBH, f_PBH binaries, mass) choices. (ii) We will add a table giving N_det(e>0.01) and N_det(e>0.1) per detector, per channel (unperturbed, binary-single, direct capture), replacing the current need to read Fig. 7 by eye — this is clearly the right way to present the paper's main quantitative claim. (iii) We will state that DECIGO's z_max = 1000 is a cap, that the Press–Schechter halo mass function and the Ref. [36] rates are extrapolations beyond their calibrated range at high z, and we will report a bounded variant — e.g. N_det computed with the integral truncated at z = 20 or 30 where the halo mass function is better constrained — so that the sensitivity of the DECIGO total to the high-z extrapolation is transparent. We note the eccentric (e>0.01) subset is less dominated by the highest redshifts than the total, since high-e sources are preferentially nearby, but the table will show this explicitly. revision: yes
- The referee's concern about the Quinlan/Sesana K-coefficient monotonic eccentricity pump cannot be fully resolved within the effective-equation framework; a definitive fix would require direct few-body scattering experiments for PBH-mass binaries in halo conditions, which is beyond the scope of this revision. We will quantify the sensitivity and state the limitation, but the high-e tail retains a model dependence we cannot eliminate.
Circularity Check
Eccentricity evolution is newly simulated and independent; absolute O(10^2) counts are linear rescalings of the authors' own prior merger rates and f_PBH normalization, a moderate but non-fatal self-citation burden.
-
self citation load bearing
[Sec. IIC, Eq. (15); Sec. IIIA (Table I); Sec. IIIB (N_det rescaling)]
"To compute N_det, we use the proper unperturbed comoving merger rate R_unperturbed from Ref. [36], assuming f_PBH = 3.4 × 10^{-3}, f_PBH binaries = 0.5, and m1 = m2 = 30 M_⊙. This combination of PBH parameters is in agreement with the current limits from the LVK Collaboration observations [16]. ... we need the comoving merger rate contributed by the halos that fall within its corresponding mass bin [M_k, M_{k+1}), which we take from Ref. [36]."
The headline forecast N_det = O(10^2) binaries with e>0.01 (and the claimed order-of-magnitude abundance improvement if none are seen) is obtained by rescaling newly simulated eccentricity histograms with R_channel(z) and f_PBH taken from the authors' own overlapping Refs. [36] and [16]. The absolute count is therefore linearly inherited from that prior normalization; only the eccentricity shape is independent in this work. This is load-bearing self-citation for the quantitative claim, not for the qualitative circularization result.
-
self citation load bearing
[Sec. IIA (binary-single channel setup); Sec. IIIB]
"Following the methodology of Ref. [35, 36], we classify PBH binaries into three evolutionary pathways... Based on earlier work in Ref. [36], we know for binary-single interactions inside halos with present masses between 10^9 M_⊙ 中 M 中 10^{15} M_⊙, that while they may have an effect on the total PBH merger rate, they are not common enough to increase very significantly the eccentricity of the PBH binaries."
Which halo-mass range is treated as 'effectively unperturbed' versus fully simulated, and the underlying environmental (rho_env, v_disp) and encounter framework that feeds Eqs. (1)–(2), are imported from the authors' prior papers rather than re-derived. The eccentricity integration itself is new, but the decision that only 10^4–10^9 M_⊙ halos matter for residual e, and the rate weights that turn shell histograms into N_det, rest on that self-cited scaffold.
full rationale
The paper's core dynamical content—integrating Peters–Mathews GW evolution plus environmental terms (Eqs. 1–2), building e(f_det) histograms for unperturbed, binary-single, and capture channels, and showing that only the halo binary-single channel retains a high-e tail into the LISA/DECIGO bands—is a genuine forward simulation, not a tautology or a fit renamed as a prediction. H and K are taken from external N-body calibrations (Quinlan 1996; Sesana et al. 2006), not from the target observable. What is load-bearing on overlapping authorship is only the absolute normalization: N_det (Eq. 15) multiplies those new eccentricity shapes by the channel-by-channel comoving merger rates R_channel(z) taken from the authors' Ref. [36] and by an f_PBH = 3.4e-3 choice tied to their Ref. [16] LVK limit. The O(10^2) e>0.01 forecast and the 'order-of-magnitude improvement' claim therefore scale directly with that inherited rate; the shape of the e distribution does not. Under the stated rules this is ordinary sequential self-citation with independent dynamical content remaining, scoring a 3—not a construction-by-definition result (6+) and not a clean 0.
Assumptions & free parameters
free parameters (6)
- f_PBH =
3.4e-3
- f_PBH binaries =
0.5
- m1=m2=m_single =
30 M_sun
- H, K interaction coefficients
- env-to-GW switch threshold =
0.5
- SNR detection threshold and tobs =
SNR≥8; tobs as in Table I
assumptions (6)
- domain assumption Orbit-averaged Peters–Mathews equations correctly describe GW-driven da/dt and de/dt outside dense environments.
- domain assumption Binary-single encounter effects on hard PBH binaries inside halos are captured by the Hρ/v and Kρ/v terms with coefficients from stellar-dynamics literature.
- ad hoc to paper Binary-single interactions are negligible before z=12 and inside present-day halos more massive than 10^9 M_sun for the purpose of eccentricity at detector frequencies.
- ad hoc to paper Characteristic strain may be approximated by the single peak harmonic n_peak(e) using the Hamers 2021 fit.
- domain assumption Comoving PBH merger rates R_channel(z) from Aljaf & Cholis (2025) Ref. [36] are accurate enough to rescale histogram counts into N_det.
- domain assumption Standard flat ΛCDM distances and the Press–Schechter-based halo mass function of Ref. [34] describe halo abundance and growth.
Cite this review
Pith. "Pith review of On the orbital eccentricities of primordial black hole binaries inside and outside of dark matter halos." pith.science (2026). https://pith.science/paper/S7IVVCOV
@misc{pith2026260723892,
author = {Pith},
title = {Pith review of: On the orbital eccentricities of primordial black hole binaries inside and outside of dark matter halos},
year = {2026},
howpublished = {\url{https://pith.science/paper/S7IVVCOV}},
note = {Machine review of arXiv:2607.23892}
}
abstract
Primordial black hole (PBH) binaries in the stellar mass range may still contribute a fraction of the detectable compact object binaries by LIGO and future GW observatories. PBH binaries at formation typically have very high eccentricities. In this paper, we study the eccentricity of stellar mass range PBH binaries from all formation channels and account for all evolutionary pathways. We simulate large samples of PBH binaries, tracking their full orbital evolution up to their merger or to the present day. For those that merge, we compute their GW strain, detectability, and eccentricity distributions for LISA, DECIGO, ET, CE, and aLIGO. We find that PBH binaries that evolve in isolation completely circularize by the time their GWs enter any GW band except for LISA's, where residual eccentricities of order $O(10^{-2})$ can exist. Binaries that become part of dark matter halos can have multiple binary-single interactions with other PBHs, especially if they reside in the more dense environments among them and can have higher eccentricities even at their late inspiral phase, probed by the GW observatories. Considering the current limits on the abundance of stellar mass range PBHs, we predict that LISA and DECIGO together would be able to probe $O(10^2)$ such binaries with $e>0.01$. If these future GW observatories in space can exclude such eccentric binaries, then limits on the PBH abundance can be improved by an order of magnitude.
Reviewed July 30, 2026 · model on record in the stance chip above.
Discussion (0). Continue with ORCID to comment.