{"id":"d42ae282-79d5-4c5f-a7f0-3a1e06925e9c","arxiv_id":"2608.18281","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":9,"one_line_summary":"CosmoPyro, a fully differentiable hierarchical Bayesian code, uses one- and two-dimensional Gaussian processes for the black hole mass distribution and infers h=0.66 and h=0.57 from GWTC-5, consistent with LVK results.","lead":"A new computer code, CosmoPyro, measures how fast the universe is expanding by using the shapes of black hole merger signals without needing a visible light counterpart. It reports values consistent with LIGO-Virgo-KAGRA within 1 sigma, and shows the method works on simulated future data.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"GP-2D's single shared correlation length in S and δ is an admitted but unquantified systematic that may bias the headline h=0.57; an anisotropic-GP rerun would settle it.","rationale":"The GP-2D shared correlation length is the most load-bearing concern because it is an admitted, unquantified systematic in the central h=0.57 measurement. The paper's own Sec. 4.2 and App. G flag it, but no analysis quantifies the bias. The simulated O5 validation is not designed to detect it, since the injected Multi Peak population has no anisotropic GP structure. A direct test with an anisotropic GP is feasible and would settle whether the reported h shift between GP-1D and GP-2D is astrophysical or a modeling artifact. I credit the paper for the careful validation against the LVK Multi Peak result (App. F), the marginalization over power-spectrum hyperparameters, and the multiple-seed convergence checks; these make the concern a condition for acceptance rather than a reason for rejection. The reader's weakest assumption identifies the same issue, so agreement is 'agree' and the recommended verdict is unchanged (CONDITIONAL).","tokens_in":30481,"tokens_out":8424,"duration_ms":89642,"concrete_test":"Implement a separable GP-2D power spectrum with independent k_cut,S and k_cut,δ (e.g., P2D(k_S,k_δ) = P1D(k_S; A_PS, k_cut,S) × P1D(k_δ; A_PS, k_cut,δ)), re-run the GWTC-5 analysis, and compare the h posterior to the isotropic result. If the anisotropic h posterior shifts by more than roughly 0.1 relative to h=0.57, the shared-length assumption is load-bearing for the headline measurement.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The load-bearing concern is the isotropic correlation structure of the GP-2D mass model, which uses a single power spectrum P2D(k) with one k_cut for both the log-total-mass coordinate S and the log-mass-ratio coordinate δ (Eq. C.3, App. C). Because the model cannot assign different smoothness to the two dimensions, the inferred k_cut is a compromise: if the true mass-ratio distribution is smooth while the primary-mass distribution has sharp features (as the GP-1D reconstruction and LVK parametric models suggest), the shared length can smear out features in S and bias the inferred mass distribution, hence h. The paper states this explicitly in Sec. 4.2: 'If the true correlation lengthscale is different for the two dimensions (S and δ), this could bias the final result.' This is load-bearing because the GP-2D headline value h=0.57 is 0.09 lower than the GP-1D value, and the simulated O5-like validation does not test the scenario: its injected truth is the Multi Peak model (Sec. 3.2), which has no anisotropic GP structure, so the isotropic GP-2D cannot be shown to be unbiased under the very misspecification the authors flag. App. G also cautions that the h posterior can be more affected for more extreme power-spectrum settings, particularly for GP-2D. The concern is not that the model is wrong, but that the central h measurement may inherit a systematic whose magnitude is currently unquantified.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"This paper introduces CosmoPyro, a JAX/NumPyro hierarchical Bayesian code for spectral-siren Hubble-constant inference. It models the BBH source-frame mass distribution either with a one-dimensional Gaussian process for the primary mass plus a running power-law mass-ratio distribution (GP-1D) or with a two-dimensional Gaussian process in the coordinates S = log(m1+m2) and δ = -log q (GP-2D), including the standard selection-effect integral evaluated via injected events. Applying the code to 231 BBH events from GWTC-5 gives h = 0.66^{+0.17}_{-0.20} (GP-1D) and h = 0.57^{+0.20}_{-0.15} (GP-2D), both consistent with the LVK spectral-siren results within 1σ. A simulated 1000-event O5-like catalog with Multi Peak injections yields h = 0.58 ± 0.09/0.11 for both GP models, recovering the injected h = 0.67 at the 1σ level. The code is validated against the LVK Multi Peak GWTC-5 analysis in App. F, and the paper includes robustness studies for power-spectrum hyperparameters (App. G) and mass-ratio misspecification (App. H). The central methodological contribution is a differentiable, GPU-accelerated nonparametric population model for spectral-siren cosmology.","tokens_in":30951,"tokens_out":5774,"duration_ms":61574,"significance":"If the central claims hold, CosmoPyro is a valuable addition to the spectral-siren toolkit. It demonstrates that NUTS-based sampling can handle the O(25,000)-dimensional parameter spaces arising from two-dimensional Gaussian-process population models, that the GP power-spectrum hyperparameters can be marginalized over, and that the resulting h constraints are statistically consistent with the LVK parametric results on the same catalog. The public code, the App. F validation against LVK pipelines, and the extensive robustness appendices are clear strengths. The most interesting scientific result is the GP-2D value h = 0.57, but that result depends on a model whose isotropic correlation structure in S and δ is acknowledged by the authors to be a possible source of bias. The significance of the paper is therefore conditional on quantifying this systematic. The framework itself is sound, and the issue appears addressable within the manuscript's scope.","major_comments":[{"comment":"The GP-2D power spectrum P2D(k) in Eq. (C.3) uses a single k_cut for both the S and δ dimensions. The paper explicitly states in Sec. 4.2 that if the true correlation lengths differ in these two dimensions, 'the inferred correlation lengthscale would be a weighted average of the two,' which can bias the final result. This is load-bearing because the GP-2D headline value h = 0.57 is 0.09 below the GP-1D value, and because the O5-like validation of Sec. 4.3 injects a Multi Peak population (Tab. 2) that has no anisotropic GP structure, so the validation cannot demonstrate unbiasedness under the misspecification the authors themselves flag. App. H tests mass-ratio misspecification only for GP-1D. Please add either an anisotropic GP-2D variant or an injection/recovery study with different correlation lengths in S and δ, and report how h shifts. Without this, the magnitude of the potential systematic on the headline h = 0.57 result is unquantified.","section":"Sec. 4.2, App. C, App. G"},{"comment":"The simulated O5-like validation uses a simplified single-event likelihood described in App. E (Gamma posteriors in detector-frame chirp mass, mass ratio, and luminosity distance, with the SNR factor sampled exactly), rather than the full gravitational-wave parameter estimation used for the real GWTC-5 events. In this setup both GP models recover h = 0.58 ± ~0.1 from an injected h = 0.67, so the injected value sits at the upper boundary of the 68% credible interval, not near the median. The claim that the pipeline 'recovers the injected Hubble constant at the 1σ level' is therefore weaker than it appears, and the single different-seed check mentioned in Sec. 4.3 is not shown quantitatively. Please provide an ensemble of independent simulated catalogs (or at least the distribution of recovered medians over several seeds), and ideally a validation on a few events with realistic PE or an end-to-end injection campaign, so the reader can assess how much of the offset is statistical noise versus a systematic of the simplified likelihood.","section":"Sec. 4.3, App. E"}],"minor_comments":[{"comment":"The description of the GP-2D window function w(S,δ) says 'For each component, it is the same as Eq. (2.6), where m_min, m_max, σf,low, and σf,high are all inferred.' It is not clear how the two component-wise windows are combined in (S,δ) coordinates, nor how the priors on m_min and m_max in Tab. 3 map onto the total-mass and mass-ratio coordinates. Please spell out the joint window definition explicitly.","section":"Sec. 2.2.1 (GP-2D)"},{"comment":"The symmetry GF(S,δ) = GF(S,-δ) is imposed by mirroring the white-noise field 'before applying the power spectrum', but the FFT-based construction with a finite grid and periodic boundaries should explain how the mirrored field is defined at δ = 0 and at the grid edge, and how the zero-slope condition at q = 1 is enforced exactly.","section":"Eq. (2.13)"},{"comment":"The simplified O5-like catalog is calibrated with a fixed amplitude A and concentration parameters k_Md ~ U(10,1000), k_q ~ U(3,7), k_dL ~ U(1,5). A few sentences comparing the resulting distance and mass-ratio error distributions to the expected O5 projections, beyond the mean relative uncertainties quoted in App. E, would help the reader judge the realism of the validation.","section":"App. E"},{"comment":"The sentence 'To adjust the LVK results with a prior of p(h) = U(0.1,2.0) to the results presented here, we reduce the LVK posteriors to this prior range' is redundant with footnote 4. Consider merging the two statements.","section":"Footnotes 4 and main text"},{"comment":"Reference [118] is cited as 'In preparation'. If the work is not publicly available, please mark it clearly as such or remove it from the main-text citation list.","section":"References"},{"comment":"The truncated-normal mass-ratio misspecification test is performed only for GP-1D. Since the main systematic concern is the GP-2D mass-ratio flexibility, it would be useful to state explicitly whether the conclusion is expected to hold for GP-2D or to run the analogous test.","section":"App. H"}],"recommendation":"major_revision","confidential_remarks":"The manuscript is technically competent and unusually honest about its limitations, but the headline GP-2D result should not be presented as a measurement until the anisotropic correlation-length systematic is quantified. The fix is feasible within the manuscript's scope - an anisotropic GP-2D test or an injection study with distinct S and δ correlation lengths - so I do not recommend rejection. The App. F validation against LVK and the public code are strong points that should be highlighted in the final version."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"CosmoPyro is a genuinely useful contribution: a differentiable JAX/NumPyro spectral-siren pipeline with a new two-dimensional Gaussian-process mass model in (log total mass, -log mass ratio), applied to GWTC-5 to get h = 0.66 (+0.17-0.20) for GP-1D and h = 0.57 (+0.20-0.15) for GP-2D. The validation against the LVK Multi Peak result in App. F shows good agreement, and the GPU implementation brings a 1000-event production analysis down to 15-25 hours. That is real, reproducible method work.\n\nWhat I think is actually new: the GP-2D model with the mirrored symmetry in delta, the joint S-delta coordinate choice, and pushing non-parametric spectral sirens to GWTC-5. The mass-ratio structure they see with GP-2D, a q distribution peaking at unity near 10 solar masses and roughly flat at 20, is a nice example of what the flexible model buys you over the running-power-law GP-1D.\n\nNow the soft spots, in order of size.\n\nWeakest point: the GP-2D model uses one isotropic power spectrum in S and delta, so a single k_cut controls smoothness in both directions. The paper says in Sec. 4.2 that if the true correlation lengths differ, the result could be biased. That is the right caveat, but it is left unquantified, and the simulated O5-like validation does not test it: the injected truth is the Multi Peak model, which has no anisotropic GP structure. So the case where GP-2D is misspecified is precisely the one not stress-tested. The headline h=0.57 sits 0.09 below the GP-1D value, and some of that gap could be this systematic. Not necessarily wrong, but the magnitude is unknown.\n\nSecond: the O5-like validation uses a simplified single-event likelihood (Gamma approximations, analytic SNR) rather than a real PE pipeline, and the recovered h = 0.58 is a 1-sigma downward shift from the injected 0.67. The paper says this shift is absent for a different random seed, but that run is not shown. That claim should be documented or shown.\n\nThird: the code link is broken. The text says the code is available at \"/githubCosmoPyro\", which is not a URL. For a methods paper that leans on a public code, that is a real reproducibility gap.\n\nThese are fixable. The central conclusions, that both GP models recover h consistent with LVK within 1 sigma and that the code works, likely survive. The anisotropic-GP concern is the one that could bite at higher precision, and I would want it quantified in the revision rather than just acknowledged. This deserves a serious referee, and I would accept it for peer review with the expectation of major revision on the validation and reproducibility points.","headline":"Solid, useful spectral-siren tool with new GP-2D mass model and new GWTC-5 h values, but the GP-2D headline number carries an admitted, unquantified systematic and the validation is weaker than it first looks.","tokens_in":31427,"tokens_out":1791,"would_cite":true,"duration_ms":20825,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"pith_extraction":{"msc":[],"pacs":[],"model":"deepseek-v4-flash","headline":"Flexible mass models recover the Hubble constant from gravitational waves.","keywords":["gravitational-wave cosmology","spectral sirens","Hubble constant","Gaussian processes","hierarchical Bayesian inference","binary black hole populations","non-parametric mass models","automatic differentiation"],"falsifier":"Run the GP-2D analysis on the same catalog with two independent correlation lengthscales, one for total mass and one for mass ratio; if the inferred $h$ shifts by more than the quoted $1\\sigma$ uncertainty, the shared-lengthscale assumption is the source.","tokens_in":30284,"feed_emoji":"🌌","tokens_out":11208,"duration_ms":95778,"temperature":0.7,"pith_summary":"This paper introduces CosmoPyro, a fully differentiable hierarchical Bayesian pipeline that infers the Hubble constant from gravitational-wave events while modeling the source-frame mass distribution as a Gaussian process rather than a fixed parametric shape. Applied to the 231 confident black-hole mergers in the fifth gravitational-wave transient catalog, the one-dimensional model returns $h = 0.66^{+0.17}_{-0.20}$ and the two-dimensional model returns $h = 0.57^{+0.20}_{-0.15}$. Both are consistent with the latest measurements from the same catalog within $1\\sigma$, even though the two models reconstruct noticeably different mass distributions. On a simulated catalog of 1000 events at anticipated next-run sensitivity, both models recover the injected $h = 0.67$ at the $1\\sigma$ level. A sympathetic reader would take this as evidence that spectral-siren cosmology can proceed without knowing the true mass distribution in advance.","feed_headline":"Flexible mass models recover the Hubble constant from gravitational waves","feed_subtitle":"With 231 black-hole mergers, two flexible mass models give h=0.66 and 0.57, matching earlier estimates.","key_machinery":"The load-bearing object is a Gaussian random field $\\mathrm{GF}$ defined on the logarithm of the mass distribution, $\\log p(\\mathrm{mass}) \\sim \\mathrm{GF}$, with smoothness controlled by a power spectrum through an amplitude $A_{\\mathrm{PS}}$ and a cutoff wavenumber $k_{\\mathrm{cut}}$. GP-1D applies the field to the primary mass $m_{1,s}$ and gives the mass ratio a running power-law distribution; GP-2D applies it jointly to $S = \\log(m_{1,s}+m_{2,s})$ and $\\delta = -\\log q$, so that correlations between total mass and mass ratio are learned rather than assumed. The field is generated on a grid in Fourier space and interpolated to event parameters, and automatic differentiation makes gradient-based sampling of the roughly 25,000 field parameters feasible. Selection effects enter through Monte Carlo estimates of the detection probability, with a variance penalty guarding against noisy likelihood estimates.","core_discovery":"The paper's claim is that a Gaussian random field placed on the logarithm of the source-frame mass distribution is a workable replacement for parametric mass models in spectral-siren cosmology, and that the Hubble constant can be recovered without committing to a particular functional form. With the one-dimensional Gaussian process (GP-1D), the inferred value is $h = 0.66^{+0.17}_{-0.20}$; with the two-dimensional process (GP-2D) that jointly models total mass and mass ratio, it is $h = 0.57^{+0.20}_{-0.15}$. The two mass distributions differ markedly, yet both $h$ posteriors remain consistent with the latest measurements within $1\\sigma$. In a simulated 1000-event O5-like catalog with a known injected value $h = 0.67$, both models return $h = 0.58^{+0.09}_{-0.09}$ and $h = 0.58^{+0.11}_{-0.10}$, recovering the truth at the $1\\sigma$ level. The paper reads this as evidence that the method does not need the true mass-distribution family to be known in advance.","pith_inferences":["Inference: the roughly $0.09$ gap between the GP-1D and GP-2D $h$ medians on the same data is a plausible measure of the current systematic budget from mass-ratio modeling; a Bayes-factor comparison, which the paper leaves for future work, would indicate which model the data favor.","Inference: the paper's shared-lengthscale caveat points to an immediate extension, a two-dimensional Gaussian process with independent correlation lengths in $S$ and $\\delta$; if that model shifts $h$, the shift quantifies the misspecification.","Inference: because the machinery scales to tens of thousands of inferred field parameters, it could absorb redshift-dependent mass evolution or mass-spin correlations, systematics the paper flags as increasingly important for larger catalogs.","Inference: the GP-2D reconstruction's flat mass-ratio distribution around $20\\,M_\\odot$ is a testable prediction for fixed-$h$ population analyses with existing catalogs."],"forward_implications":["Spectral-siren $H_0$ measurements no longer have to commit to a parametric mass model, giving a direct check on a known source of systematic bias.","If the richer correlation structure favored by the two-dimensional model is real, mass-ratio modeling choices become a leading systematic in current $h$ measurements rather than a negligible detail.","At next-run catalog sizes of about 1000 events, either flexible model constrains $h$ to roughly $\\pm 0.1$ and recovers the injected value, so the method is ready for O5-scale data.","Marginalizing over the Gaussian-process power-spectrum hyperparameters is the recommended default, because fixing them changes the reconstructed mass distribution more than it changes $h$.","The same differentiable pipeline extends to non-parametric tests of modified gravitational-wave propagation without re-parameterizing the mass model."],"supporting_citations":[{"why":"Establishes that gravitational-wave observations directly measure luminosity distance, the physical basis of spectral-siren cosmology.","marker":"[1]"},{"why":"Shows that source-population modeling choices matter for gravitational-wave cosmology, motivating non-parametric mass models.","marker":"[5]"},{"why":"Introduces the spectral-siren method using the full mass distribution, the approach this paper implements with Gaussian processes.","marker":"[6]"},{"why":"Demonstrates a one-dimensional binned Gaussian-process mass model for spectral sirens, the direct predecessor of GP-1D.","marker":"[33]"},{"why":"Models the logarithm of the heavier-mass distribution with a Gaussian process, the technique adapted to the new pipeline.","marker":"[35]"},{"why":"Provides the reference Hubble-constant measurement from the same catalog that the paper compares against.","marker":"[86]"}],"fun_headline_variants":["Gaussian process mass models deliver H0 from GWTC-5","No need to know mass shape: GPs get H0 from GW","Flexible mass models recover H0 without fixed shape","GP mass distributions yield h=0.66 and 0.57 from GW","CosmoPyro: GP mass model fits H0 without parametric shape"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The load-bearing assumption is that one shared correlation lengthscale can describe how smooth the true mass distribution is in both the total-mass and mass-ratio directions of the two-dimensional model, so if the true smoothness differs between those directions, the inferred Hubble constant can be biased.","fun_headline_variants_meta":{"raw":{"variants":["Gaussian process mass models deliver H0 from GWTC-5","No need to know mass shape: GPs get H0 from GW","Flexible mass models recover H0 without fixed shape","GP mass distributions yield h=0.66 and 0.57 from GW","CosmoPyro: GP mass model fits H0 without parametric shape"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.001387,"raw_usage":{"total_tokens":5701,"prompt_tokens":1119,"completion_tokens":4582,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":735,"completion_tokens_details":{"reasoning_tokens":4489}},"tokens_in":735,"tokens_out":4582,"duration_ms":32500,"temperature":1.0,"reasoning_tokens":4489,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-27T20:01:15.851688+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run the GP-2D analysis on the same catalog with two independent correlation lengthscales, one for total mass and one for mass ratio; if the inferred $h$ shifts by more than the quoted $1\\sigma$ uncertainty, the shared-lengthscale assumption is the source.","supporting_citations":[],"review_version":1}