{"id":"9b14b810-26e3-4c90-aa6b-67ad8e9614ad","arxiv_id":"2412.14852","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":0,"one_line_summary":"Optimal estimators and variance formulas are derived for the Legendre harmonic coefficients of the pulsar timing array Hellings-Downs correlation, with effective degrees of freedom computed for real arrays.","lead":"This paper derives optimal estimators for the harmonic coefficients of the Hellings-Downs pulsar timing correlation, expressing each coefficient's variance as the product of an effective angular and an effective frequency factor. It provides explicit formulas and a table of those factors for current pulsar timing arrays such as NANOGrav, EPTA, and IPTA.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Approach 1's matched-filter estimator (Eq. 4) is not derived and is not shown to be minimum-variance for the stated model; it is only unbiased for the fiducial expected coefficients, so the claim that both approaches yield optimal estimators is unsupported for Approach 1.","rationale":"The reader's weakest assumption focuses on the Gaussian spectral model being correct. That is a standard and shared assumption for all PTA optimal-statistic work; if it fails, both approaches fail, but it is not specific to this paper's contribution. The more load-bearing issue is that Approach 1 is presented as an 'optimal' estimator and the paper's variance formulas are advertised as the variances of optimal estimators, yet the optimality is never proven and the estimator is not the GLS solution for the stated model. The paper itself flags this ('we do not know how to derive it from first principles'), and the unproven inequalities about 2L_eff+1 reinforce the sense that the geometric interpretation is not fully controlled. Approach 2 is a standard GLS result and its variance derivation is sound, so the paper has real value; however, the abstract's claim to derive optimal estimators for both approaches is not supported. The concrete Monte Carlo or analytic check would settle whether Eq. (4) is merely an unproven but valid estimator or actually suboptimal and biased for the realized coefficients. Given that, the appropriate verdict remains conditional: accept only if the optimality claim for Approach 1 is either proven within the stated objective or removed and reframed as a matched-filter heuristic. This is consistent with the reader's CONDITIONAL verdict, but for a more specific and internal reason than the spectral-model assumption.","tokens_in":9207,"tokens_out":19495,"duration_ms":162757,"concrete_test":"Use a small synthetic PTA with N_pul=10 and N_freq=1. Compute C from AR19 and the full Fisher matrix F_{ll'} = (mu_l H)^T C^{-1} mu_{l'} H for l,l'=2..6, and verify that F is non-diagonal. Generate 10^4 Gaussian realizations of ZZ for a true coefficient vector c that equals <c> plus a perturbation in one high-l component. Compare three quantities: (i) the Approach 1 estimator (4), (ii) the GLS estimator (F^{-1} d)_l, and (iii) the variances predicted by (14) and (36). If (4) shows bias when other c_{l'} differ from their fiducial values, or if its empirical mean squared error relative to the true c_l exceeds that of the GLS estimator, then Eq. (4) is not the minimum-variance estimator of the realized c_l and the optimality claim must be revised.","verdict_should_be":"CONDITIONAL","load_bearing_attack":"The central claim is that two optimal estimators of c_l are derived and that their variances factor as <c_l>^2/[(2L_eff+1)N_freq]. The algebraic variance computation for each proposed estimator is straightforward, but Approach 1's optimality is not established. In the linear model y=ZZ with mean S c = sum_l c_l mu_l H and covariance C, the minimum-variance unbiased estimator of c_l (for fixed, unknown c) is the l-th component of (S^T C^{-1} S)^{-1} S^T C^{-1} y, i.e., Approach 2's (F^{-1} d)_l. Equation (4) instead uses only the single template mu_l H and normalizes by u_l, the inner product of that template with the fiducial mean mu H. Unless F is diagonal, the expectation of (4) equals c_l only when the true coefficients equal the fiducial <c>; otherwise it is biased for the realized harmonic coefficient. The variance reported in (14) is the scatter of this estimator about the ensemble mean, not the variance of an unbiased estimator of the actual c_l. The paper acknowledges it does not know how to derive (4) from first principles, and the assertion that it is optimal because it whitens is a detection matched-filter argument, not a minimum-variance estimation result. This does not invalidate the covariance algebra or Approach 2, but it undercuts the abstract's 'derive optimal estimators' and the claimed variance-minimizing property of Approach 1. A related admitted limitation is that the expected inequality 2L_eff+1 <= 2l+1 is unproven and Approach 1's values are not monotonic in pulsar count; that too suggests the geometric degrees-of-freedom interpretation is not fully controlled.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper develops harmonic-space estimators for the Hellings-Downs angular correlation in pulsar timing arrays. Assuming a Gaussian ensemble with known gravitational-wave-background and pulsar-noise spectra, it proposes two estimators of the Legendre coefficients c_l: Approach 1 is a per-harmonic matched-filter estimator (Eq. 4) whose variance is factored as <c_l>^2 / [(2L_eff+1) N_freq] (Eq. 14), and Approach 2 is a chi-squared-minimizing set of estimators ≈ F^{-1}d whose covariance is (F^{-1})_{ll} (Eqs. 30-36). Explicit formulas for the effective angular and frequency degrees of freedom are given in Eqs. (18)-(19) and (38)-(39), with the many-pulsar limit reducing to 2L_eff+1 → 2l+1. Numerical values of 2L_eff+1 are tabulated for current and fictional PTAs.","tokens_in":9575,"tokens_out":10868,"duration_ms":76429,"significance":"If properly qualified, the paper provides a useful and explicit variance decomposition for harmonic coefficients of the HD correlation, and the Approach 2 derivation is a clean application of standard Gaussian least-squares / Fisher-matrix theory. The many-pulsar limit correctly recovers the known 2l+1 degrees of freedom and connects to earlier work by Roebber and Holder. The paper is also transparent about its limitations: it states that it does not know how to derive Eq. (4) from first principles and that the expected inequalities for 2L_eff+1 are unproven. The main weakness is that Approach 1 is presented as an 'optimal estimator' without a supporting derivation, and the variance reported for it is not the mean-square error for estimating the realized coefficient c_l. This affects the central claim in the abstract and conclusion, though it is fixable by reframing or by supplying a derivation under an explicit optimality criterion.","major_comments":[{"comment":"The estimator in Eq. (4) is not derived from first principles. The paper itself says 'we do not know how to derive it from first principles,' and its optimality is asserted rather than proven. In the stated linear model y = sum_l c_l μ_l H with covariance C, the minimum-variance unbiased estimator of a fixed coefficient c_l is the l-th component of F^{-1}d, Eq. (33), not Eq. (4). For Eq. (4), E[ħ_c_l] = <c_l> (sum_{l'} c_{l'} F_{ll'}) / u_l, which equals c_l only in special cases, such as diagonal F or c = <c>. Hence Eq. (4) is not an unbiased estimator of the realized harmonic coefficient c_l; it is unbiased only for the ensemble mean <c_l>. Consequently the quantity called 'variance' in Eqs. (13)-(14) is the scatter of the estimator around <c_l>, not the mean-square error of an estimator of c_l. The abstract's claim that the paper derives optimal estimators for the c_l is therefore not supported for Approach 1 as written. I recommend either deriving Eq. (4) from an explicit optimality criterion (e.g., a random-effects prior with a specified prior covariance) or rephrasing Approach 1 as a matched-filter detection statistic whose ensemble variance is given by Eq. (14), and adjusting the abstract and conclusion accordingly.","section":"Approach 1: Matched Filter, Eqs. (4)-(14)"},{"comment":"The Approach 2 derivation is sound: minimizing χ^2 in Eq. (27) gives the normal equations (30), and ħ_c = F^{-1}d is unbiased with covariance (F^{-1})_{ll}. However, the notation σ^2_{ħ_c_l} is used for two different quantities. In Approach 1 it is the variance of the estimator about the ensemble average <c_l>, whereas in Approach 2 it is the variance of an unbiased estimator about the true coefficient c_l. The paper should make this distinction explicit and should not present both approaches as 'optimal estimators of c_l' with directly comparable variances. If the intended target in Approach 1 is <c_l>, that should be stated; if the target is c_l, the estimator is biased and the reported variance understates the estimation error.","section":"Approach 2: Best Fit, Eqs. (27)-(36)"},{"comment":"The split of the variance into (2L_eff+1) and N_freq is a convention, not a derived property of the estimators. The paper acknowledges this and fixes the split by two requirements, which is acceptable. However, the conclusion's phrasing that 'the variance is written as a ratio' with a denominator that is 'an effective number of degrees of freedom' should be clearly labeled as a convention in the abstract and introduction as well, so that readers do not interpret Eqs. (18)-(19) and (38)-(39) as unique physical factorizations. This is a presentational point, but it bears on how the central formulas are likely to be quoted.","section":"General framework, Eqs. (14), (36)"}],"minor_comments":[{"comment":"The symbol c_l is used both for the realized coefficient in Eq. (1) and, via ħ_c, for the estimated set in Approach 2; please introduce distinct notation to avoid ambiguity between actual coefficients, their ensemble means, and their estimators.","section":"Throughout"},{"comment":"The paper depends heavily on equations from the companion paper [7] that are not reproduced here, including the explicit form of C in (AR19), the definition of the crossover frequency between (AR34) and (AR35), and the inverse covariance in (AR35). Please include a brief appendix or explicitly list the needed definitions, since the paper cannot be fully evaluated without the companion paper at hand.","section":"Introduction and Eq. (9)"},{"comment":"The numerical entries in Table I are hard to read: for example, '2 .015 1.725' appears to combine two numbers with a space after the decimal point. Please reformat the table into separate columns with standard decimal alignment.","section":"Table I"},{"comment":"The sentence 'if not, then (18) is the effective number of degrees of freedom which could be observed with the given set, if enough SNR were available' is vague; please define what 'enough SNR' means or remove the phrase to avoid an untestable claim.","section":"After Eq. (25)"},{"comment":"The paper states that 2L_eff+1 ≤ 2l+1 is expected but unproven, and that adding pulsars does not always increase 2L_eff+1 for Approach 1. Since these statements affect the interpretation of Table I, please flag them at the first definition of 2L_eff+1 rather than only near the table.","section":"Near Table I and Conclusion"}],"recommendation":"major_revision","confidential_remarks":"The paper's Approach 2 is essentially the standard Fisher-matrix solution applied to harmonic-space PTA analysis; its value lies in the explicit variance factorization and the numerical tables. The main issue for the editor is the unsupported optimality claim for Approach 1, which appears in the abstract and conclusion. This is fixable by reframing or by adding a derivation, so I recommend major revision rather than rejection. A second editorial concern is that the paper is not self-contained and depends on the companion paper [7] for essential definitions; if [7] is not yet accepted, that dependency should be factored into the review process."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Here is my read. The paper's real contribution is a pair of estimators for the Legendre coefficients c_l of the Hellings-Downs correlation and a variance decomposition into an effective angular factor and an effective spectral factor. The second approach is a textbook generalized least squares solution of a chi-squared minimization, and the covariance algebra leading to (F^{-1})_{ll} is clean and correct. The first approach is a matched filter borrowed from the authors' earlier position-space work, and the paper honestly says it does not know how to derive it from first principles. That honesty is warranted: the estimator is not proven to be minimum-variance among unbiased estimators of the realized coefficient, and in general it is only unbiased for the ensemble mean under the fiducial model. The paper's later distinction—that Approach 1 minimizes variance between realizations while Approach 2 minimizes mismatch for our realization—partly fixes the framing, but the abstract and the \"unbiased estimator of c_l\" wording go further than the math supports. That is the main soft spot, and it is real.\n\nThe variance factorization itself is definitional rather than deep: two requirements fix the split into (2L_eff+1) and N_freq, so it is a convenient convention rather than a derived result. The genuinely new material is the effective angular degrees-of-freedom table for real PTA arrays and the demonstration that 2L_eff+1 approaches 2l+1 in the many-pulsar limit. The paper leans heavily on the authors' own previous papers for the covariance model, but the harmonic estimators and the variance formulas are not in that prior work, so the reliance is legitimate. The undemonstrated inequalities—2L_eff+1 <= 2l+1 and monotonicity in pulsar count—are explicitly flagged, and the observed non-monotonicity for Approach 1 suggests the geometric interpretation is not fully controlled.\n\nWho gets value from this? PTA practitioners who want a ready-made harmonic-space estimator with error bars. It is written in the authors' telegraphic style and genuinely requires the companion paper at hand. It is not a paradigm shift, but it is a solid methods contribution. I would send it to a serious referee. My recommendation: accept after a revision that either derives Approach 1 properly, restricts its claims to what is actually proven, or demotes it to a heuristic comparison case. The Approach 2 material alone justifies publication.","headline":"Useful harmonic-space estimators for the HD curve, with a clean GLS approach and a heuristic matched filter whose optimality claim needs tempering.","tokens_in":10063,"tokens_out":3583,"would_cite":true,"duration_ms":27377,"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":"This paper derives optimal estimators for the harmonic coefficients of the Hellings-Downs correlation and gives their variances as $\\langle c_l\\rangle^2$ divided by an effective number of angular and frequency degrees of freedom.","keywords":["pulsar timing arrays","Hellings-Downs correlation","harmonic analysis","Legendre polynomials","gravitational-wave background","optimal estimators","covariance matrix","cosmic variance"],"falsifier":"Generate simulated PTA data from a known Gaussian background and pulsar noise with a deliberately misspecified noise spectrum, then compare the empirically measured variance of the matched-filter estimator across many realizations with (14) using N_freq from (19); a systematic mismatch would show the spectral factorization is not valid outside the assumed ensemble.","tokens_in":8974,"feed_emoji":"📡","tokens_out":5425,"duration_ms":41696,"temperature":0.7,"pith_summary":"The paper is trying to show that the angular correlation pattern of a gravitational-wave background seen by pulsar timing arrays, the Hellings-Downs curve, can be reconstructed optimally in harmonic form as a sum of Legendre polynomials. It derives two families of estimators for the coefficients $c_l$: one that minimizes the variance of each coefficient in isolation and one that globally fits the whole curve, and it computes the variance of both. The variance takes the form $\\sigma^2 = \\langle c_l\\rangle^2 / [(2L_{\\rm eff}+1) N_{\\rm freq}]$, where $2L_{\\rm eff}+1$ is an effective number of angular degrees of freedom set by the pulsar sky geometry and $N_{\\rm freq}$ is an effective number of frequency bins set by the signal and noise spectra. In the limit of many pulsars uniformly spread over the sky, $2L_{\\rm eff}+1$ becomes $2l+1$, matching the usual multipole degeneracy. A sympathetic reader would care because the variance formula tells how accurately future PTA datasets can measure each harmonic and thus how well deviations from the HD prediction can be probed.","feed_headline":"Optimal variance found for each HD harmonic coefficient","feed_subtitle":"New formula splits the error into angular and frequency degrees of freedom; many pulsars recover the familiar 2l+1.","key_machinery":"The central object is the covariance matrix $F_{ll'} = (\\mu_l H)^t C^{-1} (\\mu_{l'} H)$, built from the inverse covariance $C^{-1}$ of the pulsar-pair cross-correlation data, the spectral matrix $H$, and the Legendre templates $\\mu_{l,ab} = P_l(\\cos\\gamma_{ab})$. Its diagonal entries set the variance of the matched-filter estimator, and its inverse $F^{-1}$ gives the variance of the global-fit estimator. The identity that carries the argument is the factorization of the variance into angular and frequency degrees of freedom: $2L_{\\rm eff}+1$ is defined via the geometry-only matrix $G_{ll'}$ in the crossover-frequency limit and then shown to hold generally, while $N_{\\rm freq}$ absorbs the spectral dependence. In the many-pulsar limit, $G_{ll'}$ becomes diagonal $\\frac{2l+1}{2\\langle c_l\\rangle^2}\\delta_{ll'}$ because Legendre polynomials are orthogonal over the sky, which yields $2L_{\\rm eff}+1 \\to 2l+1$.","core_discovery":"On the paper's own terms, the discovery is that the harmonic coefficients $c_l$ of the HD correlation can be estimated from PTA data by two optimal linear procedures, and that the estimation variance of each coefficient is exactly the square of its expected value divided by the product of an effective number of angular degrees of freedom (depending only on pulsar sky positions) and an effective number of frequency degrees of freedom (depending on the gravitational-wave and pulsar-noise spectra). For the matched-filter approach the variance is $\\sigma^2_{\\hat{c}_l} = \\langle c_l\\rangle^2 F_{ll}/u_l^2$, written as (14); for the global chi-square approach it is $(F^{-1})_{ll}$, written as (36). In the many-pulsar, uniform-sky limit the angular factor reduces to $2l+1$, so the variance becomes $\\langle c_l\\rangle^2 / [(2l+1) N_{\\rm freq}]$, generalizing the single-frequency-bin result of Ref. [10] to arbitrary frequency content.","pith_inferences":["Beyond the paper: because the angular factor $2L_{\\rm eff}+1$ depends only on geometry, the same table can be used to compare different arrays' sensitivity to high-$l$ harmonics without specifying a noise model.","Beyond the paper: the dirty-map/clean-map analogy suggests that a deconvolution step like $F^{-1}$ could be applied iteratively to PTA maps, and the variance $(F^{-1})_{ll}$ supplies the natural error bars for such cleaned maps.","Beyond the paper: a testable extension is to apply the two estimators to publicly released PTA timing residuals and check whether the measured $c_l$ for $l\\ge2$ follow the predicted $1/l^3$ falloff within the quoted variances.","Beyond the paper: if the conjectured inequality $2L_{\\rm eff}+1 \\le 2l+1$ holds, then the many-pulsar limit is an upper bound on the angular information any PTA can extract; testing it numerically for random arrays would be a quick check."],"forward_implications":["PTA collaborations can compute per-harmonic error bars for their reconstructed HD curves using only the pulsar sky positions and assumed spectra, via (18) and (19) or (38) and (39).","Because the harmonic coefficients with $l<2$ are predicted to be zero, measured values of $c_0$ and $c_1$ provide consistency checks of the Gaussian-background model.","Adding more pulsars or increasing observation time reduces the variance by increasing $N_{\\rm freq}$ and the effective angular degrees of freedom, which can suppress cosmic-variance-like effects.","In the many-pulsar limit, the variance scales as $1/(2l+1)$ per frequency bin, matching and generalizing the earlier single-frequency result of Ref. [10]."],"supporting_citations":[{"why":"Companion paper that supplies the position-space optimal estimators, the covariance matrix C, and the crossover-frequency limit used throughout.","marker":"[7]"},{"why":"Provides the geometry matrix G and its many-pulsar inverse, used to evaluate 2L_eff+1 in the uniform-sky limit.","marker":"[8]"},{"why":"Source of the expected harmonic coefficients ⟨c_l⟩ for the HD curve, which enter the variance formulas as the numerator.","marker":"[9]"},{"why":"Earlier single-frequency-bin harmonic-space variance result that the many-pulsar formula (26) generalizes.","marker":"[10]"},{"why":"Background harmonic analysis of PTA angular correlations that motivates expressing the HD curve as a Legendre sum.","marker":"[11]"},{"why":"Related variance calculation for the HD correlation in position space, providing context for the frequency-weighted estimators.","marker":"[13]"}],"fun_headline_variants":["Optimal variance formulas for each HD harmonic","PTA harmonic errors split into sky and frequency","Generalized 2l+1 law for arbitrary frequency bins","New estimators for pulsar timing array correlations","Dirty and clean maps for PTA angular harmonics"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The load-bearing premise is that the gravitational-wave signal and each pulsar's noise are random processes whose statistical properties (their spectra) are known, so the covariance matrix C and its inverse used to weight the estimators are correct; if those spectral models are wrong, the estimators are not optimal and the quoted variances are wrong.","fun_headline_variants_meta":{"raw":{"variants":["Optimal variance formulas for each HD harmonic","PTA harmonic errors split into sky and frequency","Generalized 2l+1 law for arbitrary frequency bins","New estimators for pulsar timing array correlations","Dirty and clean maps for PTA angular harmonics"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000209,"raw_usage":{"total_tokens":1422,"prompt_tokens":977,"completion_tokens":445,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":593,"completion_tokens_details":{"reasoning_tokens":372}},"tokens_in":593,"tokens_out":445,"duration_ms":3632,"temperature":1.0,"reasoning_tokens":372,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-11T11:50:51.786859+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Generate simulated PTA data from a known Gaussian background and pulsar noise with a deliberately misspecified noise spectrum, then compare the empirically measured variance of the matched-filter estimator across many realizations with (14) using N_freq from (19); a systematic mismatch would show the spectral factorization is not valid outside the assumed ensemble.","supporting_citations":[{"cited_title":"Allen and J","cited_arxiv_id":null,"evidence_quote":"Provides the geometry matrix G and its many-pulsar inverse, used to evaluate 2L_eff+1 in the uniform-sky limit."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Source of the expected harmonic coefficients ⟨c_l⟩ for the HD curve, which enter the variance formulas as the numerator."},{"cited_title":"Roebber and G","cited_arxiv_id":null,"evidence_quote":"Earlier single-frequency-bin harmonic-space variance result that the many-pulsar formula (26) generalizes."},{"cited_title":"Allen, Pulsar timing array harmonic analysis and source angular correlations, Phys","cited_arxiv_id":null,"evidence_quote":"Background harmonic analysis of PTA angular correlations that motivates expressing the HD curve as a Legendre sum."},{"cited_title":"Allen, Variance of the Hellings-Downs correlation, Phys","cited_arxiv_id":null,"evidence_quote":"Related variance calculation for the HD correlation in position space, providing context for the frequency-weighted estimators."}],"review_version":1}