{"id":"7b29b488-bf69-4bef-92a4-6358daab2752","arxiv_id":"2608.08419","paper_version":1,"verdict":"ACCEPT","confidence":"MODERATE","novelty_score":5.0,"correctness_risk":"low","formal_verification":"none","parameter_count":0,"one_line_summary":"FORTVSH is a validated Fortran toolbox that computes vector spherical harmonics accurately and efficiently, with batch routines stable to ℓ=2000 and built-in poloidal-toroidal magnetic field representation.","lead":"This paper presents FORTVSH, a new open-source Fortran library that computes vector spherical harmonics, including the poloidal-toroidal forms used in magnetohydrodynamics simulations. The package ships with validation tests showing machine-precision accuracy and batch routines that run about an order of magnitude faster than single-mode loops.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"High-ℓ stability claim rests on in-house tests; external benchmark at ℓmax=2000 would settle it.","rationale":"The reader's weakest_assumption identifies exactly the same load-bearing concern: the high-ℓ stability of the Holmes-Featherstone recurrence underlies the whole batch VSH stack, and it is validated only in-house. I agree this is the most important risk to the central claim. However, I do not think it changes the verdict. The paper provides thorough internal validation, machine-precision results at moderate ℓ, explicit stability studies, and a reproducible code archive. The stability concern is a verification gap rather than a demonstrated failure; the appropriate disposition remains ACCEPT with moderate confidence, possibly noting that an external high-ℓ benchmark would strengthen the claim. My proposed concrete test is exactly the kind of independent check that would either confirm the ℓmax=2000 range or reveal its boundary. I also note the paper's self-reported limitation (no inverse transform) and the reproducibility of the build across two toolchains as positive evidence. No ad hominem, no manufactured objection: the central argument is sound as far as it goes, but the most load-bearing premise deserves an external check.","tokens_in":13751,"tokens_out":1653,"duration_ms":14381,"concrete_test":"Run an independent external benchmark: compute normalized associated Legendre functions P̄ℓm(cosθ) with FORTVSH's ASSOC_LEGENDRE_NORM_ALL and compare against SHTOOLS' PlmBar (or an independent arbitrary-precision mpmath reference) for ℓmax=2000, with a dense sweep of colatitudes including θ=0.001, π−0.001, and a scan over all m for several fixed ℓ (e.g., ℓ=500, 1000, 1500, 2000). Report the maximum absolute error and the worst (ℓ,m,θ) combination. If any error exceeds ~1e-12, the claimed stability range is not supported.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central claim that every VSH batch routine is numerically stable to ℓmax=2000 depends entirely on the Holmes and Featherstone (2002) modified forward-column recurrence as implemented in ASSOC_LEGENDRE_NORM_ALL (Appendix A). The paper's validation compares this routine against mpmath reference values, but only for the diagonal m=ℓ case in Figure 1 (unnormalized) and for aggregate maximum absolute differences in Figure 2. The normalized associated Legendre functions are known to be numerically delicate at high degree and order, especially near the poles (X near ±1), because the recurrence coefficients grow and cancellation can occur. The paper does not report a fixed-θ, worst-case error scan across all m for ℓmax=2000, nor does it compare against an independent high-ℓ external benchmark such as SHTOOLS or a direct ratio-of-factorials evaluation at selected (ℓ,m,θ) points. If the recurrence loses accuracy for some m near ℓmax or for colatitudes very close to 0 or π, the advertised stable range of the batch VSH routines would be overstated, and the headline claim of machine precision over ℓmax=2000 would weaken. This is not an internal inconsistency—the tests as reported are self-consistent—but the extrapolation to 'stability inherited by every VSH routine' is only as strong as the underlying Legendre recurrence, which is asserted rather than independently demonstrated.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"FORTVSH is an open-source Fortran module that computes vector spherical harmonics in the L2-eigenfunction and polar (lambda) representations, along with scalar spherical harmonics, associated Legendre functions, Clebsch-Gordan and Wigner coefficients, and Geppert-Wiebicke coupling coefficients. The paper derives the relevant analytic identities (curl formulas, inner products, rotation relations), describes single-mode and batch routines, and presents validation against mpmath reference values and closed-form Geppert-Wiebicke test cases, together with application examples for dipole synthesis, boundary matching, and spectral decomposition. The headline claims are machine-precision evaluation of the implemented VSH forms, numerical stability of normalized Legendre batch routines to lmax=2000, order-of-magnitude speedups over single-mode loops, and suitability as a dependency-free foundation for poloidal-toroidal MHD spectral codes.","tokens_in":13963,"tokens_out":7917,"duration_ms":82127,"significance":"Should the claims hold, FORTVSH fills a concrete gap: SHTOOLS and SPHEREPACK do not natively provide the radial-poloidal-toroidal VSH triad used in MHD poloidal-toroidal decomposition. The manuscript's strengths are its open-source release with an archived version, validation against external mpmath references and independent closed-form analytic results, explicit statement of package limitations (forward-transform only, no inverse transforms yet), and the demonstration of batch speedups. The derivation of the analytic identities is consistent with cited references and no circularity was found. The principal weakness is that the high-ell stability claim, which is central to the advertised advantage over single-mode routines, rests on in-house comparisons over a sampled grid rather than an independent high-degree benchmark.","major_comments":[{"comment":"The advertised stability to lmax=2000 for the normalized associated Legendre recurrence and the Discussion's statement that this stability is 'directly inherited by every SSH and VSH routine' are not sufficiently supported. Figure 2 reports aggregate maximum absolute differences over 15 uniformly spaced interior colatitudes theta in (0,pi), with no fixed-theta scan across all m near lmax and no comparison near the poles, where the Holmes and Featherstone (2002) forward-column recurrence is most sensitive. I request either an independent high-degree benchmark (for example SHTOOLS or direct ratio-of-factorials evaluation at selected (l,m,theta) points at lmax=2000) or a documented worst-case error scan over m and theta, including theta close to 0 and pi, for ASSOC_LEGENDRE_NORM_ALL and at least one representative batch VSH routine. Without this, the load-bearing claim of machine precision over the full advertised range is an extrapolation.","section":"Section 3.1 / Appendix A / Discussion"},{"comment":"There is an ambiguity about what Figure 2 plots. The text says each batch_ALL subroutine is compared pointwise against its single-mode counterpart, but the Figure 2 caption says the plotted quantity is the maximum absolute difference between batch_ALL subroutines and reference values computed in mpmath. These are different validations. If Figure 2 shows the mpmath comparison, then the batch-versus-single-mode consistency check for the vector VSH routines is not displayed; if it shows the consistency check, then the independent reference validation for VSH_TOR, VSH_POL_UP, and VSH_POL_DN is missing. Please clarify the figure content and, if both tests were run, show both.","section":"Section 3.1, 'Batch consistency' paragraph and Figure 2"},{"comment":"The stable ranges of the two evaluation tiers are stated only for the Legendre routines as a group: unnormalized routines to lmax=200 and normalized routines to lmax=2000. Because single-mode VSH routines (VSH_TOR, VSH_POL_UP, VSH_POL_DN, PVSH_*) are implemented through ASSOC_LEGENDRE_AND_DERIV, which is built on the unnormalized Bonnet recurrence, they are presumably limited to the shorter range, while only the batch routines inherit the Holmes-Featherstone stability. The manuscript should state these per-routine limitations explicitly; otherwise a user calling a single-mode VSH function at l=1000 might wrongly assume the advertised lmax=2000 stability applies.","section":"Appendix A / Section 3.1"}],"minor_comments":[{"comment":"Equation (45) is empty; there is no equation content between Eqs. (44) and (46). Please either provide the missing expression or renumber the equations.","section":"Section 2.1, Eq. (45)"},{"comment":"The caption says 'Dotted line at 2.22 x 10 16' but should read 2.22 x 10^-16, which is the machine epsilon for double precision.","section":"Figure 5 caption"},{"comment":"The text says the VSH and Geppert-Wiebicke routes agree 'to within a few units in the last place of double precision', but Figure 4 axes are labeled 'relative deviation (%)'. These two metrics are incompatible; please state the reported quantity consistently.","section":"Section 3.1, paragraph after Eq. (67)"},{"comment":"The sentence 'All tests were generated with suitable l max determined by stability studies' is not reproducible; please give the actual lmax values and grid sizes used in Figures 1 and 2, either in the captions or in a table.","section":"Section 3.1, first paragraph"},{"comment":"The scalar spherical harmonic orthonormality test is quadrature-limited: with N=10,000 midpoint-rule points, the reported residuals of about 10^-7 are consistent with O(N^-2) truncation error, as the text acknowledges. This test should be labeled as a convergence/consistency check rather than a numerical-precision validation, or supplemented by a convergence study demonstrating that the residual decreases to machine precision as N increases.","section":"Section 3.1, SSH orthonormality test"},{"comment":"The statement 'All original weights are perfectly recovered' lacks a tolerance; please report the maximum absolute deviation between the recovered and input coefficients.","section":"Section 3.2, Part A"}],"recommendation":"major_revision","confidential_remarks":"I agree with the reader's overall positive assessment, but the high-degree stability evidence should be strengthened before acceptance. The empty Eq. (45) and the Figure 2 caption/text ambiguity are fixable presentation issues, but the external benchmark or worst-case error scan at lmax=2000 is a substantive validation gap for the paper's central claim. The author's disclosure of AI assistance in the acknowledgements is transparent and does not affect my technical assessment."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Short version: this is the first dependency-free Fortran VSH library I know of that implements the full radial-poloidal-toroidal triad, and the paper does an honest job validating it. The math is standard; the contribution is the consolidated package plus its validation suite, and that is enough to deserve referee time.\n\nWhat is new: not the identities — those come from Varshalovich, Geppert-Wiebicke, Holmes and Featherstone, and others — but the implementation, the batch routines over all l,m up to lmax, and the explicit poloidal-toroidal basis for MHD work. The tests are real: single-mode vs batch consistency, SSH orthonormality, poloidal-toroidal orthogonality, rotation inversion between the polar and L2 bases, and two Geppert-Wiebicke inner products checked against closed forms. Agreement sits around 1e-15, which is at double-precision rounding. The code is public on GitHub and Zenodo, BSD-3 licensed, and the build was verified on two toolchains. That is reproducible evidence, and I weight it.\n\nSoft spots: the lmax=2000 stability claim rests entirely on the Holmes-Featherstone recurrence as implemented and on in-house error plots against mpmath. The paper does not show a worst-case scan over all m at fixed theta near the poles, nor does it compare against an independent high-degree library like SHTOOLS. That is a real gap, but not a fatal one, because the moderate-l machine-precision claim is backed by multiple external references. The empty Equation (45) is a typo, not a mathematical hole. No inverse transform is included, but the paper states that openly.\n\nThe math holds up. The derivations are standard and the cross-checks are meaningful. The criticisms of SHTOOLS and SPHEREPACK for not supporting poloidal-toroidal decomposition look accurate. Self-citations are not a problem here; the cited works are properly credited and relevant.\n\nWho this is for: people building spectral MHD codes in spherical geometry, neutron-star and dynamo modelers, and researchers who want VSH without pulling in SHTOOLS or SPHEREPACK. The package is lightweight and embeddable.\n\nMy recommendation: send it to review. A referee with access to the code should run the high-l tests independently and ask for a fixed-theta, all-m error scan at lmax=2000 plus a comparison against SHTOOLS. Those are modest requests, not grounds for rejection.","headline":"Dependency-free Fortran VSH library with genuine validation; the high-l stability claim needs one independent benchmark, but the package deserves serious review.","tokens_in":14472,"tokens_out":2027,"would_cite":true,"duration_ms":21433,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":["95.30.Qd"],"model":"deepseek-v4-flash","headline":"FORTVSH computes vector spherical harmonics to machine precision in all standard forms, with batch routines stable to degree 2000.","keywords":["vector spherical harmonics","poloidal-toroidal decomposition","associated Legendre functions","numerical stability","magnetohydrodynamics","spectral methods","Fortran numerical library","Clebsch-Gordan coefficients"],"falsifier":"Evaluate ASSOC_LEGENDRE_NORM_ALL and the VSH batch routines at $\\ell_{\\max}=2000$ on a dense grid of colatitudes that includes points extremely close to the poles, and compare every $(\\ell,m)$ component against mpmath reference values; any deviation above the roughly $10^{-15}$ double-precision floor would refute the stability claim. A second independent check is the orthogonality test $Y^{(+1)}_{JM} \\cdot Y^{(0)}_{JM} = 0$ at the same grid points.","tokens_in":13511,"feed_emoji":"🧲","tokens_out":7951,"duration_ms":74150,"temperature":0.7,"pith_summary":"The paper presents FORTVSH, a dependency-free Fortran module for evaluating vector spherical harmonics (VSH) in the polar form used for poloidal-toroidal decomposition and in the standard angular-momentum form. The central claim is that the package computes every implemented VSH mode to double-precision rounding, preserves the analytic identities that connect the two representations, and supplies batch routines that stay stable to angular degree $\\ell_{\\max} = 2000$ while running about an order of magnitude faster than loops over single-mode calls. The motivation is spectral MHD: a magnetic field written as poloidal and toroidal scalar potentials can be advanced mode-by-mode only if the VSH basis, its curl identities, and its inner products are all computed reliably. If the stability claim holds, the package gives simulation codes a small numerical core that existing spherical-harmonics libraries do not provide, since those handle only curl-free scalar-potential fields or purely tangential two-dimensional flows.","feed_headline":"Vector spherical harmonics reach machine precision at degree 2000","feed_subtitle":"Batch routines run ten times faster and preserve the poloidal-toroidal identities MHD codes need.","key_machinery":"The load-bearing object is the normalized associated Legendre batch evaluator ASSOC_LEGENDRE_NORM_ALL, which implements the modified forward-column recurrence of Holmes and Featherstone: the $4\\pi$ normalization is absorbed into the recurrence coefficients so all intermediate values stay of order $1/\\sqrt{4\\pi}$, extending stable degree from about 200 to about 2000. Every scalar and vector spherical harmonic batch routine calls this module once per evaluation point, so its stability is inherited. Around this core, FORTVSH assembles the three polar VSH families $Y^{(-1)}$, $Y^{(0)}$, and $Y^{(+1)}$ from the scalar harmonics via radial, gradient, and cross-product formulas, provides the Clebsch-Gordan and Wigner 6-j machinery that rotates between polar and $L^2$-eigenfunction forms, and packages the curl identities needed to propagate modes through the induction equation.","core_discovery":"The central discovery is that one normalized Legendre recurrence, applied in batch mode, can carry the entire VSH construction to machine precision and to $\\ell_{\\max} = 2000$. Specifically, FORTVSH's ASSOC_LEGENDRE_NORM_ALL folds the $4\\pi$ normalization into the modified forward-column recurrence coefficients, and every scalar harmonic and VSH batch routine is built on that single module. The validation shows maximum absolute errors at or below about $10^{-15}$ relative to mpmath reference values, exact reproduction of the bilinear orthogonality $Y^{(+1)}_{JM} \\cdot Y^{(0)}_{JM} = 0$, exact reconstruction of the Clebsch-Gordan rotation between polar and $L^2$-eigenfunction bases, and agreement of the Geppert-Wiebicke coupling coefficients with closed-form analytic expressions. The paper claims these results hold for $\\ell$ up to 2000 for normalized routines and up to 200 for unnormalized Legendre routines, with batch evaluation outperforming single-mode loops by an order of magnitude.","pith_inferences":["The stability claim rests on in-house comparisons at 15 colatitude points and one azimuth, so a user who needs all modes at arbitrary points should run an independent high-$\\ell$ benchmark before trusting $\\ell_{\\max}=2000$ in production.","If the recurrence really is stable to 2000, the same code likely tolerates somewhat higher degrees, but with error growth the paper does not quantify; a simple extension would measure the crossover degree where errors exceed machine precision.","The package currently covers forward transforms only; adding an inverse transform that synthesizes fields from coefficient arrays would turn the validation examples into a full simulation-capable spectral transform library.","The Geppert-Wiebicke coupling coefficients implemented here are tailored to Hall-MHD two-mode interactions, so the natural next test is a small spectral Hall-cascade simulation that checks whether coefficient-based mode coupling reproduces direct evaluation of the nonlinear terms."],"forward_implications":["A magnetic field in poloidal-toroidal form can be evaluated, rotated between representations, and integrated in the same package, so the full VSH triad needed by the induction equation is available in one call.","Batch evaluation at $\\ell_{\\max}=2000$ makes high-degree spectral transforms practical on a single CPU core; the order-of-magnitude speedup over single-mode loops carries directly into time-advance loops.","Machine-precision preservation of orthogonality means a spectral decomposition will not leak power into spurious modes, keeping MHD mode amplitudes clean over long runs.","Because the package has no FFT or BLAS/LAPACK dependency, it can be embedded directly into purpose-built spherical-shell simulation codes.","Existing spherical-harmonics tools cannot express toroidal fields or the radial-tangential coupling of poloidal fields; FORTVSH's explicit radial-poloidal-toroidal triad fills that structural gap."],"supporting_citations":[{"why":"Supplies the tensor-spherical-harmonic definitions, the Clebsch-Gordan series, and the Wigner 6-j identities on which the VSH formulas are built.","marker":"Varshalovich et al., 1988"},{"why":"Provides the modified forward-column recurrence that ASSOC_LEGENDRE_NORM_ALL implements and that extends stability to $\\ell_{\\max}=2000$.","marker":"Holmes and Featherstone (2002)"},{"why":"Defines the I and J coupling coefficients used for the Hall-MHD inner products and cross-validated against closed forms.","marker":"Geppert and Wiebicke (1991)"},{"why":"Introduces the $Y^{(\\lambda)}_{JM}$ alternate VSH notation used for poloidal-toroidal decomposition.","marker":"Akheizer and Berestetskii (1965)"},{"why":"Establishes the alternative VSH notation compared in Table 1 and used for magnetostatic applications.","marker":"Barrera et al. (1985)"},{"why":"Provides the SHTOOLS comparison and the batch-output indexing convention adopted for FORTVSH arrays.","marker":"Wieczorek and Meschede, 2018"},{"why":"Provides the SPHEREPACK comparison, showing why a purely tangential 2D vector transform cannot represent radial-poloidal coupling.","marker":"Adams and Swarztrauber, 1999"},{"why":"Establishes the poloidal-toroidal decomposition of magnetic fields that motivates the package's design.","marker":"Krause and Rädler, 1980"},{"why":"Supplies the standard angular-momentum definitions of the spin eigenfunctions used in the VSH construction.","marker":"Edmonds, 1960"}],"fun_headline_variants":["VSH hits machine precision at degree 2000","Batch VSH: 10x speed, 1e-15 accuracy","FORTVSH: normalized recurrence to degree 2000","Machine precision vector spherical harmonics","VSH routines reach 2000 with roundoff error"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The advertised range holds only if the Holmes-Featherstone normalized Legendre recurrence stays accurate to $\\ell_{\\max}=2000$ for every m and every colatitude; the paper supports this with in-house error plots against mpmath rather than an independent high-degree benchmark.","fun_headline_variants_meta":{"raw":{"variants":["VSH hits machine precision at degree 2000","Batch VSH: 10x speed, 1e-15 accuracy","FORTVSH: normalized recurrence to degree 2000","Machine precision vector spherical harmonics","VSH routines reach 2000 with roundoff error"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000317,"raw_usage":{"total_tokens":1748,"prompt_tokens":854,"completion_tokens":894,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":470,"completion_tokens_details":{"reasoning_tokens":815}},"tokens_in":470,"tokens_out":894,"duration_ms":9247,"temperature":1.0,"reasoning_tokens":815,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-14T04:37:34.957484+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Evaluate ASSOC_LEGENDRE_NORM_ALL and the VSH batch routines at $\\ell_{\\max}=2000$ on a dense grid of colatitudes that includes points extremely close to the poles, and compare every $(\\ell,m)$ component against mpmath reference values; any deviation above the roughly $10^{-15}$ double-precision floor would refute the stability claim. A second independent check is the orthogonality test $Y^{(+1)}_{JM} \\cdot Y^{(0)}_{JM} = 0$ at the same grid points.","supporting_citations":[],"review_version":1}