{"id":"03c05e40-9bc5-495e-9548-581b48a8a48e","arxiv_id":"2601.15198","paper_version":3,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":5.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":6,"one_line_summary":"Hierarchical inference on simulated LISA EMRI catalogues can tightly constrain sharp population features, separate mixtures with ~20 detections, and retain sensitivity under model misspecification.","lead":"Simulated extreme mass-ratio inspiral catalogues for LISA show that hierarchical Bayesian inference can tightly constrain black-hole population parameters and separate mixed formation channels with roughly 20 detections. The paper is a forecast that tells LISA analysts which population features will be measurable and which will remain uncertain.","discovery_kind":"extension","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Equation (5) defines Pdet as the complement of the detection probability; if the code follows the text, the selection function is inverted and every hierarchical result is biased.","rationale":"The reader's conditional verdict is reasonable, and their concern about Fisher-matrix likelihoods and SNR-only selection is a genuine limitation. However, the single most load-bearing issue I find is internal and precedes those concerns: Eq. (5) as written computes the complement of the detection probability. If that expression is what was used, the selection function is inverted and the entire hierarchical likelihood is wrong; no amount of added realism would fix it. If it is a typo, a one-line correction is needed, and the reader's external-realism concern remains the main caveat. I chose to disagree with the reader's weakest_assumption because they identified a modelling simplification while I identify a potentially fatal internal inconsistency in a central equation. The proposed test, checking the code and rerunning one case with the two definitions, settles the matter. If the code is correct and Eq. (5) is merely misprinted, the current CONDITIONAL verdict stands; if not, the verdict should be moved to REJECT or UNVERDICTED. I therefore leave the verdict unchanged pending this check.","tokens_in":25010,"tokens_out":14061,"duration_ms":152966,"concrete_test":"Clone the released analysis code (GitHub/Zenodo) and inspect the selection-function routine. Run one representative case, e.g. Model A with 10^3 injections, using the literal Eq. (5) expression P_det = 1 − P(ρ_n>ρ_t|ρ_opt), and compare the recovered x_c posterior with the paper's Figure 1. Then replace P_det with P(ρ_n>ρ_t|ρ_opt) and rerun. If the posteriors differ materially (e.g., x_c shifts by more than ~1σ), the pipeline is using the inverted selection function and the central results are invalid. If the posteriors are unchanged, Eq. (5) is a typo and should be corrected in the manuscript.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The selection function is the load-bearing part of the central claim: Eq. (3) requires α(Λ)=∫dθ P_det(θ)p_pop(θ|Λ), and Eq. (4) defines P_det=H(ρ_n−ρ_t). The expected detection probability for a source is therefore P(ρ_n>ρ_t|ρ_opt). However, Eq. (5) as written evaluates P_det = 1 − P(ρ_n>ρ_t|ρ_opt), which is the probability of non-detection, not detection. If the code implements this expression, then every Monte Carlo α(Λ) is the undetected fraction. Inserting such α into Eq. (2) inverts the selection correction and would systematically bias all recovered hyperparameters, so the reported unbiased recovery of injected values is not expected. This is the first thing to check because it is internal to the machinery that produces the headline claims about sharp features, ~20 detections, and misspecification resilience. Fixing event-level realism, as the reader suggests, would not matter if the selection function itself is inverted.","agreement_with_reader":"disagree"},"referee_report":{"model":"deepseek-v4-flash","summary":"This paper presents a hierarchical Bayesian inference framework for EMRI population inference with LISA, using an ML-based selection-function emulator (poplar), Fisher-matrix Gaussian event likelihoods, and nested sampling. The framework is applied to simulated catalogues drawn from two single-component population models (A and B), their homogeneous and heterogeneous mixtures, and deliberately misspecified models. The headline results are tight recovery of parameters controlling sharp features, disentangling of mixed populations with as few as ~20 detections, and the claim that key population features remain detectable even under model misspecification. The paper is structured as a simulation study with injection–recovery tests and posterior predictive checks.","tokens_in":25270,"tokens_out":8304,"duration_ms":89172,"significance":"If the results hold, the paper provides a useful proof-of-principle for EMRI population inference with LISA. Its strengths include end-to-end injection–recovery tests, posterior predictive checks, explicit treatment of selection effects, and reproducible code/data releases. The conceptual framework is standard hierarchical Bayesian inference; the novelty is the EMRI application with a fast emulator for the selection function. However, the headline claims about ~20 detections and percent-level branching-fraction recovery are stronger than the presented 90% intervals support, and the selection-function definition in Eq. (5) is internally inconsistent as written. The central derivation is standard and defensible once the selection-function typo is resolved, but the quantitative claims need to be either better supported or softened.","major_comments":[{"comment":"There is a sign error/inconsistency in the definition of Pdet used in the selection-function Monte Carlo. Eq. (4) defines Pdet = H(ρ_n − ρ_t), whose expectation is P(ρ_n > ρ_t | ρ_opt). However, the text after Eq. (5) states Pdet = 1 − P(ρ_n > ρ_t | ρ_opt). If the code implements the latter, α(Λ) in Eq. (2) is the undetected fraction, which would invert the selection correction and bias all recovered hyperparameters; the reported unbiased recovery of injected values would then be unexpected. This is load-bearing for all results. Please correct the equation/text and confirm, e.g. by reporting the computed α and the detected fraction for the simulated catalogues, that the implemented quantity is the detection probability. If this is only a typesetting error, it must still be fixed because Eq. (5) as written is what a reader would implement.","section":"§2, Eqs. (4)–(5)"},{"comment":"The claim that mixed populations can be disentangled with ~20 detections is stronger than the evidence. At 10^2 injections (19 detections) for Model A+B, the branching fraction is recovered as w = 0.70+0.06−0.07 (90%), i.e. roughly ±9% relative, not percent-level. More importantly, the subdominant B-component mass slopes are essentially unconstrained: λ_M^B = −2.1+1.0−0.8 and λ_μ^B = −2.5+0.9−1.3. At 10^3 injections (198 detections) these mass slopes remain broad. The data support detecting the presence of a second component and coarsely estimating its weight, but they do not support the statement that the two populations are 'disentangled' at ~20 detections if disentangling means recovering the subdominant component's mass distribution. Please either soften the claim, e.g. to 'detecting the presence of a second component,' or provide a quantitative criterion for disentangling that is ac","section":"§4.2, Fig. 3; Abstract; Conclusion"},{"comment":"The conclusion states that for the A+B mixture, 'the branching fraction is recovered with percent-level accuracy for 10^2 injections (19 detections).' This is inconsistent with Fig. 3, where the 90% interval is w = 0.70+0.06−0.07; the fractional uncertainty is ~9%, not percent-level. Percent-level recovery of w is only achieved at 10^3–10^4 injections (w = 0.73+0.02−0.02 and 0.696+0.007−0.007). Please correct this quantitative summary.","section":"Conclusion"},{"comment":"The forecast is built on an SNR-only selection function with a fixed threshold (ρ_t = 20) and on Fisher-matrix Gaussian event likelihoods. The authors acknowledge the Fisher approximation is valid in the high-SNR regime, but the quantitative uncertainties quoted throughout (e.g. 'within 1.5%' at 1886 detections) should be presented as idealized. Real LISA selection and parameter estimation may be non-Gaussian and parameter-dependent, which could change the forecasted constraints and the '~20 detections' conclusion. This is not a reason to reject, but the limitations should be stated more prominently in the abstract and conclusion, and the claims should be framed as conditional on the adopted idealised detection model.","section":"§2, selection function and event likelihood"}],"minor_comments":[{"comment":"The summation index runs from k=0 to N_t; presumably it should be k=1 to N_t.","section":"§2, Eq. (5)"},{"comment":"The text says 'for all analyses of mixed-population models, we fix the branching fraction at w=0.7,' but Table 3 and Fig. 3 treat w as a free hyperparameter. Please clarify that w=0.7 is the fixed value used to generate the injected catalogues, while w is later inferred.","section":"§3, mixed models"},{"comment":"The sentence 'For we adopt a modified Schechter-like distribution' is missing a model reference; likely 'For Model A we adopt...'.","section":"§3.1"},{"comment":"The final paragraph mentions 'strong constraints on spin–eccentricity correlations,' but the population models in Table 1 treat spin and eccentricity as independent distributions with no correlation parameter. If no correlation was modelled, this statement should be removed or rephrased as constraints on the marginal spin and eccentricity distributions.","section":"Conclusion"},{"comment":"The selection-function emulator is central to the analysis, but its validation is only cited to the authors' prior papers. Please add a brief summary of the validation (e.g. emulator error versus direct SNR evaluation) or an explicit statement that the end-to-end injection–recovery tests in this paper constitute the validation.","section":"§2, poplar emulator"}],"recommendation":"major_revision","confidential_remarks":"The manuscript relies heavily on the authors' own emulator and prior papers; this is not circular in itself, but given the Eq. (5) inconsistency, the editor should ask the authors to verify that the implemented selection function is the detection probability and to release or point to the code used for the injection–recovery runs. The overclaims about ~20 detections and percent-level branching-fraction recovery are fixable by rewording, but they are prominent enough that they should be corrected before publication."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Read this with the stress-test note in hand. The paper is a solid simulation study—the natural next step in a program the authors already built. What's genuinely new is systematic hierarchical inference tests for mixed EMRI populations (heterogeneous and homogeneous) and for deliberately misspecified models, using the poplar ML emulator and Fisher-matrix event likelihoods. Recovery on simulated catalogues is clean, and the posterior predictive checks are a real plus. The finding that sharply featured distributions—Schechter mass peak, narrow spin peak—are much better constrained than broad power laws is useful and well demonstrated. The misspecification results are also sensible: simple models flatten structure, flexible models mimic dominant components.\n\nThe soft spots, in order:\n\n1. Equation (5) as written is wrong or at least mis-stated. Eq. (4) defines Pdet as detection probability; Eq. (5) says Pdet = 1 − P(ρn > ρt | ρopt), the complement. If the code actually implements that, the selection function is inverted and every reported recovery would be biased. Since the recoveries are clean, I suspect a typo in the text, but this is the first thing to ask for: correct the equation, confirm the code uses H(ρn − ρt), or show a direct selection-function validation. You cannot accept the paper with that contradiction standing.\n\n2. The abstract's \"disentangled with ~20 detections\" is stronger than Fig. 3 supports. At 19 detections the branching fraction is recovered within about 10% and the subdominant component's mass parameters are barely constrained. The tight 2% numbers appear at 198–1980 detections. The claim should be reworded to say the dominant component and spin peak are recoverable at ~20, while full disentangling needs hundreds.\n\n3. The event-level realism is a legitimate caveat: Fisher-matrix Gaussian likelihoods and SNR-only selection with ρt=20 are standard shortcuts, but they're not cross-checked against full EMRI PE or a realistic detection pipeline. That's fine for a forecast as long as the numbers are presented as indicative, which they mostly are.\n\nThe paper is worth a serious referee. It will be cited by anyone building LISA EMRI population analyses. Fix the Eq. 5 issue and soften the abstract, and it's publishable.","headline":"Solid EMRI population forecast built on the authors' own pipeline; the '~20 detections' claim is overstated and Eq. 5 has a sign inconsistency that must be checked before acceptance.","tokens_in":25788,"tokens_out":2579,"would_cite":true,"duration_ms":28231,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":[],"model":"deepseek-v4-flash","headline":"A LISA-scale catalogue of extreme mass-ratio inspiral events can reveal the underlying massive black hole population, including sharp features and mixed formation channels, even when the assumed population model is imperfect.","keywords":["gravitational wave astronomy","extreme mass-ratio inspirals","LISA","hierarchical Bayesian inference","massive black hole populations","selection effects","population mixtures","model misspecification"],"falsifier":"Run the same hierarchical inference on simulated catalogues produced with a full EMRI detection pipeline and complete posterior sampling, instead of an SNR cut and Fisher-matrix likelihoods, and check whether the recovered mass-spectrum peak and mixture fraction fall within the forecast credible intervals; a systematic offset beyond those intervals would show the SNR-only and Fisher assumptions are the load-bearing simplification.","tokens_in":24884,"feed_emoji":"🛰️","tokens_out":3876,"duration_ms":46882,"temperature":0.7,"pith_summary":"Extreme mass-ratio inspirals—compact objects spiraling into massive black holes—will be observed in large numbers by LISA, and each event carries sub-percent measurements of the black hole's mass and spin. This paper asks whether those individual measurements can be combined into a trustworthy census of the massive black hole population, including which formation channels produced the events. The answer it argues for is yes: hierarchical Bayesian inference on simulated catalogues recovers the population's parameters, resolves mixtures of distinct formation channels with as few as roughly twenty detections, and remains sensitive to dominant population features even when the model used is deliberately misspecified. A sympathetic reader would care because this is the statistical machinery that will turn a stream of LISA detections into astrophysical conclusions about how massive black holes grow and how their nuclear environments shape inspirals.","feed_headline":"Twenty detections can separate black-hole formation channels","feed_subtitle":"Hierarchical inference on simulated LISA data recovers black-hole population features even when the model is wrong.","key_machinery":"The central object is the hierarchical hyperlikelihood, which combines per-event parameter-estimation likelihoods with a population model and normalises by a selection function that accounts for undetected sources. Detectability is modelled as a step function in signal-to-noise ratio above a threshold, approximated by Monte Carlo integration and made computationally tractable through machine-learning emulators. Event-level measurement uncertainty is handled with Fisher-matrix Gaussian likelihoods, appropriate at the high signal-to-noise ratios considered, and the population-level posterior is explored with nested sampling. This machinery lets the paper transform a catalogue of simulated dete","core_discovery":"The paper demonstrates that hierarchical population inference on simulated LISA EMRI catalogues can recover the hyperparameters of parametrised population models, and that the constraining power depends on the structure of the distribution being measured. Sharply defined features—such as the Schechter peak in the massive black hole mass function or a narrow spin distribution—are recovered with high precision, while smooth power-law slopes are harder to pin down. Mixed populations can be disentangled with roughly twenty detections, with the branching fraction of a 70/30 mixture recovered to a few percent. Under model misspecification, the recovered distributions track the dominant component,","pith_inferences":["Beyond the paper: the 'sharp features are easier to constrain than smooth slopes' pattern likely generalises beyond EMRIs to any gravitational-wave population study, including ground-based binary catalogues, where peaked mass or spin features should be the first targets for precision inference.","Beyond the paper: the roughly-twenty-detection threshold is contingent on an SNR-only selection function and Fisher-matrix measurement errors; a full detection pipeline and full posterior sampling could shift the threshold, so the number should be read as an order-of-magnitude guide rather than a fixed promise.","Beyond the paper: a directly testable extension is to replace the fixed branching fraction with a redshift-dependent formation-rate model, which would show whether the inference can separate channels that evolve differently across cosmic time.","Beyond the paper: the misspecification results suggest a practical analysis strategy—run a small suite of flexible phenomenological models and use posterior predictive checks to identify functional mismatch before interpreting recovered population parameters as physical."],"forward_implications":["With a few hundred EMRI detections, the peak of the massive black hole mass function can be constrained to roughly 1.5 percent, and narrow spin distributions to sub-percent precision.","A mixture of two formation channels can be separated with as few as about twenty detections, with the mixture fraction recovered to a few percent.","Even when two subpopulations share the same functional form but differ in their parameter values, the inference identifies bimodality rather than collapsing to a single averaged population.","When the assumed model is wrong, the resulting biases are systematic and directional: dominant population features are tracked, subdominant features are smoothed over, and posterior predictive checks flag the tension.","Key population trends remain detectable even under misspecification, so simplified phenomenological models can still yield meaningful astrophysical constraints when model selection is performed jointly with hierarchical inference."],"fun_headline_variants":["Twenty EMRI catches can separate black-hole birth channels","LISA's first 20 inspirals map black-hole populations","Sharp features, not slopes, betray black-hole origins in EMRI data","Even a wrong model still catches black-hole population features"],"cache_read_input_tokens":2304,"weakest_assumption_plain":"The forecasts assume that whether an EMRI is detected depends only on its signal-to-noise ratio exceeding a fixed threshold, and that per-event parameter uncertainties are Gaussian as described by the Fisher matrix; if real LISA detection and parameter estimation deviate from these assumptions, the predicted uncertainties and the roughly-twenty-detection threshold could change.","fun_headline_variants_meta":{"raw":{"variants":["Twenty EMRI catches can separate black-hole birth channels","LISA's first 20 inspirals map black-hole populations","Sharp features, not slopes, betray black-hole origins in EMRI data","Even a wrong model still catches black-hole population features"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000579,"raw_usage":{"total_tokens":2545,"prompt_tokens":701,"completion_tokens":1844,"prompt_tokens_details":{"cached_tokens":256},"prompt_cache_hit_tokens":256,"prompt_cache_miss_tokens":445,"completion_tokens_details":{"reasoning_tokens":1784}},"tokens_in":445,"tokens_out":1844,"duration_ms":13403,"temperature":1.0,"reasoning_tokens":1784,"cache_read_input_tokens":256,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-03T08:55:03.849022+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run the same hierarchical inference on simulated catalogues produced with a full EMRI detection pipeline and complete posterior sampling, instead of an SNR cut and Fisher-matrix likelihoods, and check whether the recovered mass-spectrum peak and mixture fraction fall within the forecast credible intervals; a systematic offset beyond those intervals would show the SNR-only and Fisher assumptions are the load-bearing simplification.","supporting_citations":[],"review_version":1}