{"id":"f7ac846a-e4f5-44eb-bd6b-10fe4c955cc2","arxiv_id":"2504.21133","paper_version":2,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":7.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":1,"one_line_summary":"Closed-form analytic covariance templates for the 2, 3, and 4 point galaxy correlation functions, built from a 1/k power-law power spectrum, reproduce Boltzmann-code results at percent level and trace sparsity to closed-triangle configurations.","lead":"This paper derives closed-form formulas for the covariance of galaxy clustering measurements using a simplified power-law model of the matter power spectrum. The formulas match numerical calculations at roughly percent level and explain why high-order clustering covariances are sparse, which matters for surveys like DESI, Euclid, Roman, and SPHEREx.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Equation (5.43) omits the triangular-support Heavisides on its shot-noise term, so the analytic f-integral is nonzero in regions where the defining integral (5.1) vanishes.","rationale":"The reader's weakest assumption concerned the sketched regularization of divergent equal-argument contributions. My review found a more concrete, checkable defect in the same family of manipulations: when the Dirac-delta closure is used in Eq. (5.24), the sifting property is applied without enforcing that s lie within the u-integration interval, so the shot-noise term in Eq. (5.43) loses its triangle-inequality support. This is not a question of UV regularization; it is an algebraic omission that makes the headline analytic f-integral disagree with the integral that defines it. Because the f-integrals are the building blocks of the 3PCF and 4PCF covariance matrices, the analytic covariance expressions built on Eq. (5.43) inherit the error. The paper's numerical comparisons of model versus camb power spectra remain informative, but they do not test the closed-form Eq. (5.43) itself, since the figures appear to use direct numerical evaluation of Eq. (5.1) with the model P(k). The fix is straightforward: restore the Heaviside support factor to the shot-noise term (and, correspondingly, to Eq. 5.25). For that reason I do not think the reader's conditional verdict should be changed to rejection; the concern is substantive but addressable, and a single point check settles it. I mark agreement as partial because the reader identified a related but distinct regularization worry, not the missing triangular support in the shot-noise term.","tokens_in":67704,"tokens_out":15164,"duration_ms":167408,"concrete_test":"Evaluate Eq. (5.43) for {ℓ,ℓ′,ℓ″}={0,0,0} at ri=r′_j=1 h^-1 Mpc and s=10 h^-1 Mpc with nbar=3×10^-4 h^3 Mpc^-3 and A=277 h^-2 Mpc^2, and compare with the exact shot-noise term of Eq. (5.1): (1/2π² nbar) ∫_0^∞ dk/k sin(k)² sin(10k), which is analytically zero because no closed triangle can be formed. A nonzero value from Eq. (5.43) confirms that the missing Heaviside factor H(R_+−1)H(1−R_−) must be restored to the 2π/(A nbar s) term. A simple numerical quadrature of the same configuration over k=10^-3 to 10 h Mpc^-1 should reproduce the analytic zero; if instead it tracks Eq. (5.43), the printed formula is invalid.","verdict_should_be":"UNCHANGED","load_bearing_attack":"In deriving the shot-noise part of the analytic f-integral, the paper applies the closure/sifting step from Eq. (5.24) to Eq. (5.25) without carrying the integration bounds: the u-integral of δ(s−u) over [|ri−r′_j|, ri+r′_j] is zero unless s lies in that interval. The correct Eq. (5.25) should contain the factor H(s−|ri−r′_j|)H(ri+r′_j−s), i.e. H(R_+−1)H(1−R_−). Instead, Eq. (5.43) places the resulting term 2π/(A nbar s) after all three Heaviside terms, so it contributes for every s. This makes the closed-form f-integral disagree with its definition even at the most basic level. For example, take ℓ=ℓ′=ℓ″=0, ri=r′_j=1, s=10 (no closed triangle). The shot-noise contribution to Eq. (5.1) is proportional to ∫_0^∞ (dk/k) sin(k)² sin(10k), which vanishes by Dirichlet integrals, but Eq. (5.43) gives a nonzero value ∝ 1/(8π nbar s ri r′_j). Thus the central building block for the 3PCF and 4PCF covariances is wrong as written outside the triangular region, directly contradicting the paper's sparsity mechanism. This is also an instance of the paper's broader 'measure-zero' treatment of delta/equal-argument terms not being safe: after sifting, the delta's support constraint is not measure-zero in the final integration variable s.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"This paper develops an analytic Gaussian-random-field covariance model for galaxy 2-, 3-, and 4-point correlation functions. The physical power spectrum is replaced by the model P(k)=A/k+1/\\bar n; using standard spherical-Bessel integral identities, the authors obtain a closed form for the 2PCF covariance (Eq. 3.16), closed forms for the f-integrals that are the building blocks of the 3PCF and 4PCF covariances (Eqs. 5.6 and 5.43), and approximate sparsity of the resulting covariance matrices by the requirement that ri, r'j, and s form closed triangles. The model is validated against camb-based numerical integrals and against the full 3PCF and 4PCF covariance matrices, and the paper proposes a rank-one correction scheme for inversion together with a complexity discussion.","tokens_in":68116,"tokens_out":9215,"duration_ms":98539,"significance":"If the derivations are correct, this is a useful and potentially practical contribution: it provides an interpretable, fast, fully analytic template for covariance matrices that are otherwise too high-dimensional to estimate from mocks, and it gives a concrete structural explanation for their sparsity. The paper's strengths are the breadth of the derivations, the explicit closed forms, the extensive numerical comparisons against camb, and the clear discussion of triangular-region support. The weaknesses are that one of the two central closed forms, Eq. (5.43), is wrong as written outside the triangular region, and that the treatment of divergent equal-argument terms is asserted rather than rigorously derived; both issues affect the central claims of the paper, so they must be addressed before the results can be accepted.","major_comments":[{"comment":"The sifting step from Eq. (5.24) to Eq. (5.25) drops the integration bounds of the u-integral. Since the u-integral in Eq. (5.24) runs over u \\in [|ri - r'j|, ri + r'j], the integral \\int du \\delta(s-u) (...) is nonzero only when |ri - r'j| \\le s \\le ri + r'j. The correct Eq. (5.25), and hence the shot-noise term 2\\pi/(A \\bar n s) in Eq. (5.43), must carry the factor H(s-|ri-r'j|)H(ri+r'j-s), equivalently H(R_+ - 1)H(1 - R_-). As printed, that term sits outside all three Heaviside factors and contributes for every s. For the simple case \\ell=\\ell'=\\ell''=0, ri=r'j=1, s=10, Eq. (5.43) gives a nonzero shot-noise contribution proportional to 1/(8\\pi \\bar n s ri r'j), whereas the defining integral (5.1) vanishes because \\int_0^\\infty dk \\, \\sin^2 k \\, \\sin(10k)/k = 0 by the Dirichlet integral. Thus the closed-form f-integral disagrees with its definition outside the triangular region, and the sparsity mechanism derived from products of Eq. (5.43) is not supported as written.","section":"Section 3.1 and Section 5.1"},{"comment":"The treatment of divergent equal-argument spherical-Bessel contributions is only sketched. Footnote 5 states that binning eliminates the ultraviolet divergence of I^{[2,lin]}_\\ell(r,r') at r=r', and Section 5.1 drops the ri=s contribution to Eq. (5.6) because \"the integral of a function at one point vanishes.\" This is not sufficient once products of f-integrals are integrated over s: a term supported at a single point can become a finite contribution after multiplication by a Dirac-delta shot-noise term (as in Eq. (5.6)) or after binning. The missing support constraint identified in Eq. (5.43) is an explicit example of a pointwise/support issue that is not measure zero in the final integration variable. The authors need to provide a rigorous regularization or an explicit binned derivation for these equal-argument terms before the f-integral results and the covariance comparisons built on them can be considered complete.","section":"Section 3.1 and Section 5.1"},{"comment":"The accuracy summary is presented as the mean, standard deviation, and maximum of the absolute difference D between model and camb f-integrals, yet the text and table captions conclude \"single-digit percent level accuracy\" and quote the largest difference as \"1.19%\" or \"2.22% away from a perfect match.\" As printed, the columns are dimensionful absolute differences, not relative percentages; for example, Table 2's Max(AbsVal)=1.19\\times10^{-2} does not by itself establish 1.19% agreement. Please report relative errors with an explicit normalization (for example, relative to the typical or maximum |f| over the plotted region), or provide percent-error maps, so that the central accuracy claim is directly verifiable.","section":"Tables 2 and 3, Section 6.1.3"}],"minor_comments":[{"comment":"The amplitude A=277 h^{-2} Mpc^2 is hand-tuned to one power spectrum and one number density, and Appendix A notes that it should be re-tuned for other number densities. Please state explicitly how the quoted validation accuracy depends on this choice and whether the reported f-integral accuracy is partly in-sample.","section":"Section 2.1"},{"comment":"It would be helpful to state explicitly whether the \"model\" f-integrals and covariance matrices in the figures are evaluated from the closed forms (5.6) and (5.43) or by direct numerical integration of Eq. (5.1) with P(k)=A/k+1/\\bar n. This matters for assessing the implications of the support error in Eq. (5.43).","section":"Section 6.1, Figures 7-13"},{"comment":"The phrase \"the largest difference is 1.19% away from a perfect match\" is inconsistent with the numerical column label Max(AbsVal). If these are absolute differences, the percentage wording should be removed or replaced with a properly normalized relative error.","section":"Page 44, Table 2 caption"},{"comment":"After correcting the missing Heaviside factors, the statement that Eq. (5.25) matches equation 3.21 of [205] should be rechecked, since the cited result presumably includes the appropriate support constraint in its domain of validity.","section":"Section 5.2.1, Eq. (5.25)"}],"recommendation":"major_revision","confidential_remarks":"The manuscript is within scope for JCAP and the topic is timely for upcoming spectroscopic surveys. The main technical issue, the missing triangular support in the shot-noise term of Eq. (5.43), is substantial but appears fixable by carrying the integration bounds through the sifting step and recomputing the affected comparisons. I recommend major revision rather than rejection, provided the authors also address the equal-argument regularization and the normalization of the accuracy metrics."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"The 2PCF covariance closed form and the general f-integral machinery are genuinely new, and the derivations are mostly careful. The structural insight that products of f-integrals suppress non-overlapping triangle regions is nice and likely correct in spirit. But the stress-test note is right, and it lands on the main equation.\n\nThe problem is in the sifting step from Eq. (5.24) to (5.25). The u-integral of δ(s−u) over [|ri−r′j|, ri+r′j] is only nonzero when s lies in that interval. That constraint has to appear as Heaviside factors in Eq. (5.25), and it must survive into Eq. (5.43). Instead, the 2π/(A nbar s) shot-noise term sits outside all the H(1−R+) and H(R±−1) factors, so it contributes for every s. For ℓ=ℓ′=ℓ″=0, ri=r′j=1, s=10, the defining integral vanishes by the Dirichlet integral, but Eq. (5.43) gives a nonzero value. So the claimed closed-form f-integral disagrees with its own definition outside the triangle region, and the sparsity story is contradicted by the shot-noise piece.\n\nThis is a fixable error, not a fatal one. But it means Eq. (5.43) cannot be used as published, and the validation tables do not currently test whether this term is correct away from the triangle region.\n\nOther soft spots are proportionate: the equal-argument ultraviolet regularization is only sketched (footnote 5 and Section 5.1), and may deserve more care; the amplitude A is fitted to the same camb spectrum used for validation, so the percent-level metrics are partly in-sample; and no code or data are provided, which limits independent checking.\n\nWhat the paper does well: the 2PCF covariance (3.16) is a clean closed form; the f-integral derivation is systematic and uses legitimate special-function identities; the comparison against camb covers many multipole channels and separations; and the triangle-overlap explanation of covariance sparsity is a useful contribution to the literature.\n\nThis is a paper for LSS methodologists working on analytic covariance templates. It deserves a serious referee, but the referee should push for a corrected shot-noise term, a re-validation that checks the support behavior explicitly, and ideally a public implementation.\n\nRecommendation: send to peer review, but expect major revision before acceptance.","headline":"The f-integral formulas are elegant but Eq. (5.43) drops the triangular-support Heavisides on the shot-noise term, so the paper's central analytic result is wrong as written.","tokens_in":68569,"tokens_out":5936,"would_cite":false,"duration_ms":64221,"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":"Using the power spectrum model $P(k)=A/k+1/\\bar n$, the paper derives closed-form expressions for the Gaussian-random-field covariance matrices of the 2-, 3-, and 4-point correlation functions and traces their sparsity to the triangle…","keywords":["covariance matrices","N-point correlation functions","Gaussian random field","spherical Bessel functions","f-integrals","triangle inequality sparsity","shot noise","power-law power spectrum"],"falsifier":"Numerically evaluate equation (5.1) directly for $P(k)=A/k+1/\\bar n$ on a very fine $k$ grid for fixed $\\{\\ell,\\ell',\\ell''\\}$, without dropping the $r=r'$, $s=u$, and $r_i=s$ points, and compare the binned result with equations (5.6) and (5.43); a mismatch at the few-percent level at small separations would indicate that the discarded pointwise terms survive binning.","tokens_in":67533,"feed_emoji":"📐","tokens_out":7290,"duration_ms":68034,"temperature":0.7,"pith_summary":"The paper aims to replace the numerically integrated power-spectrum integrals in analytic covariance templates for galaxy clustering with fully closed forms. Using the model $P(k)=A/k+1/\\bar n$, it derives an exact expression for the 2PCF covariance in terms of the geometric mean and ratio of pair separations, and closed forms for the f-integrals that are the building blocks of the 3PCF and 4PCF covariances. It reports single-digit-percent agreement with the Boltzmann-solver-based integrals across many multipole channels and separations. The structural payoff is an explanation of sparsity: each f-integral is only large when the two side lengths and the pair separation can form a closed triangle, so products of f-integrals inherit a smaller overlap region. If the model holds, future surveys can compute and invert covariance templates without expensive oscillatory integrals or many mock catalogs.","feed_headline":"One power-law spectrum turns galaxy covariances into closed forms","feed_subtitle":"Percent-level agreement with numerically integrated power spectra also explains why these matrices are sparse.","key_machinery":"The load-bearing object is the f-integral, a triple product of spherical Bessel functions weighted by the power spectrum, defined in equations (4.3) and (5.1). The paper evaluates it after substituting $P(k)=A/k+1/\\bar n$, splitting into integrals $I^{[3,\\mathrm{lin}]}$ and $I^{[3,\\mathrm{quad}]}$ with powers $k$ and $k^2$. The $k^2$ piece collapses to a Dirac delta via the spherical Bessel closure relation; the $k$ piece is reduced, using an orthogonality identity to insert a free angular order, to finite sums over Wigner 3-j and 6-j symbols, binomial coefficients, and hypergeometric and Meijer G-functions. The mechanism that produces sparsity is the set of Heaviside functions $H(1-R_{+,ij})$ and $H(R_{-,ij}-1)$ in equation (5.43), which confine the dominant contribution to the region where the three lengths form a closed triangle. Multiplying f-integrals then shrinks the nonzero region to the overlap of their triangular supports.","core_discovery":"The central discovery is that the $1/k$ power-law plus shot-noise model makes every integral needed for the leading-order (Gaussian-random-field) covariance of the 2-, 3-, and 4-point correlation functions analytically tractable. For the 2PCF, the covariance reduces to equation (3.16), expressed through $g=\\sqrt{rr'}$, $\\chi=r/r'$, and the scale $\\eta=1/(A\\bar n)$ where cosmic variance balances shot noise. For the higher-order functions, the f-integral $f_{\\ell,\\ell',\\ell''}(r_i,r'_j,s)=\\int_0^\\infty (k^2 dk/2\\pi^2)P(k)j_\\ell(kr_i)j_{\\ell'}(kr'_j)j_{\\ell''}(ks)$ is evaluated in closed form in two cases: one spherical Bessel function with zero order and argument (equation 5.6), and all arguments nonzero (equation 5.43), the latter expressed through finite sums over Wigner symbols, hypergeometric functions, and Meijer G-functions. The paper argues that these closed forms match integrals against the true power spectrum at the single-digit-percent level, and that the Heaviside functions in the closed form enforce the triangle inequalities $|r_i-r'_j|\\le s\\le r_i+r'_j$. This triangle constraint is the origin of the covariance sparsity: products of f-integrals are significant only where their triangular regions overlap.","pith_inferences":["If the triangle-support picture is generic, sparsity patterns in high-order covariances could be predicted from pure geometry (side-length triangle inequalities) rather than from evaluating the full integrals, suggesting targeted compression schemes that only store overlapping-triangle blocks.","The $A/k$ model has no baryon acoustic oscillation wiggle; the residual errors at low $\\{\\ell,\\ell'\\}$ and the diagonal features in the half-inverse tests may be correctable by adding a small $k$-dependent correction to $A$ while preserving closed forms, for example a sum of power laws $P(k)=\\sum_a A_a k^{-\\alpha_a}$, since the same Bessel integral machinery applies term-by-term.","A direct testable extension is to compare the analytic precision matrix from the inversion lemma, with the correction estimated from a small number of mocks, against the precision matrix from thousands of mocks; the paper's sparsity argument predicts the correction is dominated by a few low-order multipole blocks."],"forward_implications":["The 2PCF covariance in equation (3.16) can be evaluated element-by-element from $g$, $\\chi$, $\\bar n$, and $V$ without numerical integration over $k$.","The closed-form f-integrals make the GRF covariance of the 3PCF and 4PCF computable by evaluating special functions at precomputed ratios $R_{\\pm,ij}=|r_i\\pm r'_j|/s$, with cost dominated by the $s$-grid rather than by a $k$-grid.","Covariance sparsity is explained geometrically: an element is appreciable only when the relevant side lengths and separations satisfy the triangle inequalities, and products of f-integrals are appreciable only in the overlap of their triangular regions.","Because the analytic model is invertible, the true covariance can be written as the analytic template plus a low-rank correction, enabling use of a matrix inversion lemma to obtain the precision matrix without inverting a huge mock covariance.","The same f-integral formulas extend to N-point correlation functions beyond the 4PCF, since higher-order GRF covariances are built from products of these same building blocks."],"supporting_citations":[{"why":"Supplies the GRF 2PCF covariance integral that the paper's analytic evaluation starts from.","marker":"[28]"},{"why":"Supplies the isotropic-basis 3PCF and 4PCF covariance formulas and defines the f-integrals the paper evaluates.","marker":"[176]"},{"why":"Provides the isotropic N-point basis functions used to write the covariance expansions.","marker":"[177]"},{"why":"Provides the analytic triple-spherical-Bessel integral used for the $k^2$ piece of the f-integrals.","marker":"[205]"},{"why":"Integral tables used for the double-spherical-Bessel integrals with constant and linear weights.","marker":"[202]"},{"why":"Gives the orthogonality trick that rewrites $k^n j_\\ell(kr)$ and reduces triple Bessel integrals to double integrals.","marker":"[206]"},{"why":"Matrix inversion lemma used to invert the template-plus-correction covariance.","marker":"[207]"},{"why":"Supplies the claim that precision matrices are sparse, which this paper's triangle argument explains.","marker":"[126]"}],"fun_headline_variants":["Closed forms for 2-4pt galaxy covariances from a 1/k spectrum","Triangle inequalities explain galaxy covariance sparsity","Analytic galaxy covariances to percent accuracy via power-law","One power law makes 2-4pt covariance matrices tractable","Power-law spectrum yields closed-form covariance for galaxy surveys"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The derivation assumes that the divergent contributions that occur when two spherical Bessel arguments coincide (and orders match) can be discarded because the covariance is binned in separation, so these points occupy measure zero; if binning or subsequent integration over $s$ gives those points finite weight, the closed-form f-integrals and the covariance built from them miss real pieces.","fun_headline_variants_meta":{"raw":{"variants":["Closed forms for 2-4pt galaxy covariances from a 1/k spectrum","Triangle inequalities explain galaxy covariance sparsity","Analytic galaxy covariances to percent accuracy via power-law","One power law makes 2-4pt covariance matrices tractable","Power-law spectrum yields closed-form covariance for galaxy surveys"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000302,"raw_usage":{"total_tokens":1850,"prompt_tokens":1163,"completion_tokens":687,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":779,"completion_tokens_details":{"reasoning_tokens":601}},"tokens_in":779,"tokens_out":687,"duration_ms":6530,"temperature":1.0,"reasoning_tokens":601,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-16T05:13:09.089150+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Numerically evaluate equation (5.1) directly for $P(k)=A/k+1/\\bar n$ on a very fine $k$ grid for fixed $\\{\\ell,\\ell',\\ell''\\}$, without dropping the $r=r'$, $s=u$, and $r_i=s$ points, and compare the binned result with equations (5.6) and (5.43); a mismatch at the few-percent level at small separations would indicate that the discarded pointwise terms survive binning.","supporting_citations":[{"cited_title":"On a General Method for Resolving Integrals of Multiple Spherical Bessel Functions Against Power Laws into Distributions","cited_arxiv_id":"2112.07809","evidence_quote":"Gives the orthogonality trick that rewrites $k^n j_\\ell(kr)$ and reduces triple Bessel integrals to double integrals."}],"review_version":1}