{"id":"2f29cfff-d3d1-4eee-b8d8-863083fba1fc","arxiv_id":"1908.03222","paper_version":2,"verdict":"ACCEPT","confidence":"MODERATE","novelty_score":7.0,"correctness_risk":"low","formal_verification":"none","parameter_count":3,"one_line_summary":"Transit light curves and their parameter derivatives are computed in closed form for polynomial limb-darkening laws of any integer order, with improved numerical stability and speed.","lead":"This paper derives exact mathematical formulas for the brightness dip when a planet crosses a star with any polynomial edge-darkening law, including the slope of the dip with respect to each model parameter. The formulas are fast, numerically stable, and shipped as open-source code in four programming languages.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"No significant objection identified: the derivation is internally consistent and independently validated; remaining gaps are disclosed limitations, not threats to the central claim.","rationale":"The reader accepted the paper with moderate confidence, identifying the spherical/sharp-edge physical idealization as the weakest assumption. I agree that idealization limits real-data applicability, but it is explicitly disclosed and is not a correctness risk for the mathematical central claim. My independent check of the algebra found no internal inconsistency in the change of basis, the Green's theorem reduction, the recursion for M_n, or the derivative chain rule. The strongest independent evidence is the paper's cross-validation against multiple existing codes and its open-source implementations, which make the formulas practically reproducible. The only additional concern I found is that Sec. 6 does not spell out special cases for b = 0 and complete occultation for the higher-order terms; the k2-based branches in Eqs. (71)-(78) appear to assume partial overlap or an occultor fully inside the source. This is a genuine completeness gap in the text, but it is confined to measure-zero or geometrically distinct regimes, it does not undermine the derivation where the branches apply, and it is likely handled in the code. A single targeted test at (r, b) = (2.0, 0.5) can confirm whether the implementations cover the missing branch. Because the central claim is mathematically sound and the practical limitations are disclosed, I would keep the reader's ACCEPT verdict unchanged.","tokens_in":33718,"tokens_out":19049,"duration_ms":188656,"concrete_test":"Run the shipped Julia/C++ implementation at N = 10 with (r, b) = (2.0, 0.5), corresponding to complete occultation, comparing double-precision output against BigFloat; if the flux and all derivatives are zero to about 1e-14, the missing full-occultation branch is handled by the code and the Sec. 6 gap is only a documentation issue.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The paper's central claim is the exact construction F = I0 s^T A u with s built from the M_n recursion and analytic derivatives (Eqs. 93-98). I checked the change-of-basis algebra (Eqs. 10-18), the Green's-theorem reduction leading to Eq. (73), and the normalization I0 = 1/[pi(g0 + 2g1/3)]; these are internally consistent, and the paper's comparisons against Mandel & Agol (2002), Pal (2008), batman, and PyTransit provide independent numerical support. The physical idealizations (spherical bodies, azimuthal symmetry, sharp edges) are explicitly listed in Sec. 12 and are standard modelling assumptions, not hidden circularity. The self-reported precision loss for N > 25 (Sec. 10.1) is also disclosed. The one unstated fragility I find is in Sec. 6: s_n for n >= 2 is defined only through the P integral and the k2 branches, with no explicit special-case formulas for b = 0 or for complete occultation (b < r - 1). In those regimes the text's equations are not obviously valid, although the shipped code may handle them. This is a completeness gap in the exposition, not an error in the central derivation, and it does not change the verdict.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"This paper derives closed-form analytic expressions for the flux of a spherical star occulted by a circular, sharp-edged body, for limb-darkening profiles that are polynomials in the stellar elevation of arbitrary integer order. The flux is expressed as F = I0 s^T A u, where A is a fixed change-of-basis matrix, u contains the limb-darkening coefficients, and the solution vector s is computed by Green's theorem as a line integral plus recursion relations for the integrals M_n. The authors also derive analytic derivatives with respect to the radius ratio, impact parameter, and limb-darkening coefficients, and they compute exposure-time-averaged light-curve derivatives using adaptive Simpson quadrature. The paper validates the expressions against 256-bit BigFloat arithmetic and against Mandel & Agol (2002), Pál (2008), batman, and PyTransit, and it reports substantial improvements in numerical stability and speed.","tokens_in":33975,"tokens_out":13401,"duration_ms":153496,"significance":"If the derivation is correct, this is a substantial methodological contribution. The central construction is parameter-free in the sense that no fitted constants enter; the formulas follow from geometry and calculus, and the claims are supported by machine-checked derivations, reproducible code in multiple languages, and independent numerical benchmarks. The stability improvements near contact points and the availability of analytic derivatives are directly useful for gradient-based inference, Fisher-matrix calculations, and high-precision JWST-era transit modeling. The limitations imposed by the spherical, azimuthally symmetric, sharp-edged assumptions are disclosed in Section 12 and are standard for the field. The main gap I find is not in the core derivation but in the treatment of special geometric limits for the higher-order terms.","major_comments":[{"comment":"","section":"§6, Eqs. (71)-(78)"},{"comment":"","section":"§5, Eqs. (53)-(56)"}],"minor_comments":[{"comment":"","section":"§3.1"},{"comment":"","section":"Table 1"},{"comment":"","section":"§12"},{"comment":"","section":"§9-10"}],"recommendation":"major_revision","confidential_remarks":"The core derivation is sound and the paper is likely to be an important reference for transit modeling. The reason I recommend major revision rather than acceptance is the incomplete treatment of special geometric limits for higher-order limb darkening: the manuscript should either supply the missing formulas for b=0 and full occultation or clearly state that these cases are handled by the code with validated limiting procedures. I do not see circularity or hidden fitted parameters; the self-citations to starry are topical and appropriate."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"You should know this paper for what it is: the first closed-form treatment of transit light curves for polynomial limb darkening of arbitrary integer order, complete with analytic derivatives for radius ratio, impact parameter, and limb-darkening coefficients. The derivation is explicit and self-contained, building on Green's theorem and a well-chosen Green's basis. The change-of-basis algebra checks out, the recursion relations for the M_n integrals are internally consistent, and the formulas reduce properly to the known uniform, linear, and quadratic cases. The authors also deliver a genuinely improved overlap-area formula for the uniform case, stabilized for the numerically nasty limits near contact, and they provide code in Julia, Python, IDL, and C++ with tests and reproduction scripts.\n\nThe validation is a real strength. They compare against 256-bit BigFloat arithmetic near the contact points and against Mandel & Agol (2002), Pál (2008), batman, and PyTransit. The speed benchmarks are credible: roughly 35-45% faster than EXOFAST, 70% faster than Pál's Fortran code, and 5-30x faster than the Giménez/PyTransit recursion for high orders. The self-reported precision degradation around N = 25-30 is disclosed, and the paper correctly cautions users to check high-order cases.\n\nThe soft spots are minor. The stress-test note is right: in Sec. 6 the higher-order s_n are defined only through the P integral and the k^2 branches, with no explicit special-case formulas for b = 0 or complete occultation (b < r - 1). Those regimes are physically important and the equations as written are not obviously valid there; the shipped code likely handles them, but the paper should say so. Also, the exposure-time-averaged derivatives are computed with adaptive Simpson quadrature rather than analytically; that is a practical choice, but it is a numerical step in an otherwise analytic paper. Neither issue affects the central claim.\n\nOn the physical side, the assumptions of spherical bodies, azimuthal symmetry, and sharp edges are explicitly stated in Sec. 12 and are standard across the field. They limit real-data applicability at JWST precision, but the paper acknowledges this and explains why a numerically stable ideal model still matters for fitting and gradient-based inference.\n\nThis is a useful and reproducible methods paper. The target audience is anyone fitting transits or eclipses with high-precision photometry, especially JWST-era users who want fast gradients for Hamiltonian MCMC or Fisher information. I would cite it and bring it to a reading group. It deserves a serious referee: the derivation is original, the validation is thorough, and the code is a genuine community resource.","headline":"A carefully derived, well-validated analytic transit model with gradients for arbitrary integer-order polynomial limb darkening; the central result is solid and the minor exposition gaps do not undermine it.","tokens_in":34512,"tokens_out":2621,"would_cite":true,"duration_ms":29000,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":[],"model":"deepseek-v4-flash","headline":"For polynomial limb darkening of any integer order, the transit light curve is closed-form, stable, and differentiable.","keywords":["transit light curves","limb darkening","occultation","closed-form solution","elliptic integrals","numerical stability","analytic derivatives","exoplanet photometry"],"falsifier":"Compute the model at r = b + $10^{-12}$ in double precision and compare it with the same formulas evaluated in 256-bit arithmetic, for polynomial order N = 30 with coefficients chosen to alternate in sign; the paper's own error analysis predicts fractional errors near $10^{-6}$ at that order, so if double-precision evaluation instead shows larger or uncontrolled errors, the claim of uniform stability would fail at high order. For low order, N = 3 or less, the same test should stay within about $10^{-14}$ of the high-precision value.","tokens_in":33519,"feed_emoji":"🔭","tokens_out":8929,"duration_ms":88832,"temperature":0.7,"pith_summary":"This paper claims that the transit or occultation light curve of a star with polynomial limb darkening is not limited to the usual linear and quadratic cases: for any integer order N, the flux is a closed-form expression built from a solution vector that depends only on the radius ratio and impact parameter. The same construction yields analytic partial derivatives of the flux with respect to the radius ratio, impact parameter, and all limb-darkening coefficients. A sympathetic reader would care because the formulas are designed to stay numerically stable precisely where earlier analytic models lose precision, near ingress and egress contacts and when the planet's edge crosses the star's center, and because the derivatives make gradient-based fitting and Hamiltonian Monte Carlo practical. The paper also shows that finite-exposure, time-averaged fluxes and their derivatives can be computed accurately, and that the model applies to eclipsing binaries and multiple-star systems under spherical symmetry.","feed_headline":"Transit light curves become exact for any limb-darkening order","feed_subtitle":"Closed-form flux and analytic derivatives promise faster, numerically stable exoplanet fits for space-based photometry.","key_machinery":"The load-bearing construction is the Green's basis for the limb-darkening profile: the profile is re-expressed as a sum of functions with the constant term 1, the linear term z, and, for n >= 2, the combinations (n+2)z^n - n z^(n-2), so that each term becomes the exterior derivative of a known vector field and the disk integral turns into a boundary line integral around the occultor. That line integral is evaluated as a combination of integrals Mn(r,b) obeying upward or downward recursion relations, seeded by the general complete elliptic integral cel(kc,p,a,b) and, for the uniform and quadratic pieces, by a numerically stable circle-overlap area formula using a Kahan-style triangle construction and atan2 instead of arccos. The same Mn values are reused in the analytic derivatives, so gradients add little computational cost.","core_discovery":"The central claim, stated on the paper's own terms, is that when the specific intensity is a polynomial in the elevation coordinate z, the observed flux factorizes as F = I0 s^T A u, where u lists the limb-darkening coefficients, A is a fixed change-of-basis matrix from those coefficients to a Green's basis, and s is a solution vector depending only on the radius ratio r and impact parameter b. The uniform, linear, and quadratic terms s0, s1, and s2 are given in stabilized form, and every higher-order term sn is obtained by reducing the two-dimensional surface integral, via Green's theorem, to a line integral around the occulting body whose primitive integrals are linear combinations of recursively defined integrals Mn(r,b). The paper derives analytic partial derivatives of F with respect to r, b, and the limb-darkening coefficients, including small-impact-parameter expressions for the b-derivative that avoid division by b. The authors claim that these formulas evaluate at near-double precision across the parameter plane, with the documented caveat that precision degrades for very high polynomial order, roughly N = 25-30, because of cancellations among large alternating coefficients.","pith_inferences":["Non-integer limb-darkening laws, such as the power-2 law, are not directly solvable by this machinery, but the paper's framework suggests a practical pipeline: fit a high-order polynomial to the non-linear profile once, then fit the few non-linear parameters against the precomputed polynomial coefficients.","The documented precision ceiling near N = 25-30 implies that 'arbitrary order' has a real double-precision limit; users wanting very flexible profiles should regularize the coefficients or stay in the Green's basis, where high-order terms contribute zero net unocculted flux and are easier to keep small.","The same recursion-and-line-integral scheme could be reused for other two-circle boundary problems, such as ring shadows or spot crossings, whenever the integrand is polynomial in z.","Because the region near r = b is now accurate, perturbative transit-timing-variation searches that expand timing shifts as a linear response to the derivative of the light curve can drop their interpolation patches and use the derivative directly."],"forward_implications":["Transit fits can treat limb-darkening coefficients as free parameters at any integer polynomial order without resorting to numerical integration, which is what high-precision space photometry is expected to require.","The analytic derivatives make gradient-based optimization, Fisher-matrix forecasts, and Hamiltonian Markov Chain Monte Carlo sampling of transit parameters practical, since the gradient of the likelihood can be computed without finite differences.","Finite-exposure light curves are handled by adaptive quadrature that splits exposures at contact points, so time-averaged fluxes and their derivatives are available without re-binning or losing contact information.","The same closed-form machinery applies to eclipsing binaries and multiple-star occultations when both stars are treated as spherical and limb-darkened.","The stabilized overlap-area formula removes the need for special series patches near first-to-fourth contact, so the same code paths work for the whole parameter range rather than switching to limiting expressions."],"supporting_citations":[{"why":"Provides the standard uniform, linear, and quadratic light-curve formulas and the Lambda(r,b) elliptic-integral framework that this paper extends, corrects, and stabilizes.","marker":"Mandel & Agol (2002)"},{"why":"Previous infinite-series method for arbitrary polynomial limb-darkening order, used as the accuracy and speed baseline that the closed-form recursion supersedes.","marker":"Giménez (2006)"},{"why":"Earlier analytic derivatives for the quadratic limb-darkened model, which the paper extends to higher order and compares against.","marker":"Pál (2008)"},{"why":"Source of the Green's-basis transformation, notation, and recursion approach that the paper adapts to azimuthally symmetric limb darkening.","marker":"Luger et al. (2019a)"},{"why":"Defines the general complete elliptic integral cel used to remove logarithmic singularities and round-off cancellations in the linear and higher-order terms.","marker":"Bulirsch (1969)"},{"why":"Introduces the non-linear limb-darkening law used as the comparison model that high-order polynomial fits must approximate.","marker":"Claret (2000)"},{"why":"A widely used computational package that serves as the speed and accuracy baseline for quadratic and non-linear comparisons.","marker":"Kreidberg (2015)"}],"fun_headline_variants":["Exact transit light curves for any polynomial limb darkening","Closed-form transit flux with analytic derivatives up to any order","Exact transit light curves and derivatives for any limb darkening","Green's theorem yields exact transit light curves and derivatives","Arbitrary-order limb darkening: closed-form light curves and derivatives"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"If the two bodies are not effectively spherical with azimuthally symmetric, sharp-edged limb darkening, the closed-form model describes the idealized star rather than the actual observation.","fun_headline_variants_meta":{"raw":{"variants":["Exact transit light curves for any polynomial limb darkening","Closed-form transit flux with analytic derivatives up to any order","Exact transit light curves and derivatives for any limb darkening","Green's theorem yields exact transit light curves and derivatives","Arbitrary-order limb darkening: closed-form light curves and derivatives"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.002354,"raw_usage":{"total_tokens":9084,"prompt_tokens":975,"completion_tokens":8109,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":591,"completion_tokens_details":{"reasoning_tokens":8025}},"tokens_in":591,"tokens_out":8109,"duration_ms":49489,"temperature":1.0,"reasoning_tokens":8025,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-14T14:20:23.195850+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Compute the model at r = b + $10^{-12}$ in double precision and compare it with the same formulas evaluated in 256-bit arithmetic, for polynomial order N = 30 with coefficients chosen to alternate in sign; the paper's own error analysis predicts fractional errors near $10^{-6}$ at that order, so if double-precision evaluation instead shows larger or uncontrolled errors, the claim of uniform stability would fail at high order. For low order, N = 3 or less, the same test should stay within about $10^{-14}$ of the high-precision value.","supporting_citations":[],"review_version":1}