{"id":"0fe2d3d8-9cc5-4c86-be8c-fc1966779794","arxiv_id":"2607.18406","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":0,"one_line_summary":"Feynman integrals can be solved globally over phase space by rewriting them as second-order PDE equilibrium problems and discretizing with finite elements.","lead":"Physicists show that hard multi-loop Feynman integrals can be computed by rewriting them as partial differential equations and solving those equations over a whole phase-space region at once with finite-element meshes. The approach is a proof-of-concept on a two-loop four-point integral and could eventually speed up precision predictions for future electron-positron colliders.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"End-to-end claim is untested: reported 14 s and 1.4e-4 use exact analytic boundary and load data, excluding the cost of the inputs the method depends on.","rationale":"The reader's weakest assumption identifies the same load-bearing dependency on externally supplied boundary and lower-sector data. I agree, and I sharpen it by pointing out that the FEM load function f requires lower-sector integrals over the entire domain, not just on the boundary, so the input-data requirement is even larger than 'boundary conditions.' The PDE construction itself is internally consistent: taking derivatives of the first-order system yields a second-order equation whose coefficients are obtained from the same IBP reduction, and triangularity in sectors is preserved. The variational formulation and FEM discretization are standard. The convergence results are plausible and the algebraic rates are consistent with piecewise-linear FEM on domains with mild singularities. The concern is therefore about the interpretation and practical reach of the central claim, not about a mathematical contradiction. The paper itself concedes that the numerical experiments are a proof of concept on an analytically well-controlled example and that boundary conditions are input data, so the appropriate disposition remains conditional: the method is promising and worth testing end-to-end, but the headline efficiency claim is not yet demonstrated. No adjustment to the reader's verdict is needed.","tokens_in":10017,"tokens_out":9043,"duration_ms":84344,"concrete_test":"Run an end-to-end benchmark on the same 12-MI family: compute γ_n on the boundary nodes and all lower-sector MIs on the interior mesh using a path-based method (e.g., AMFlow or DiffExp) with controlled accuracy; then assemble f and solve the top-sector FEM system. Compare the resulting top-sector errors against AsyInt at the 75,008 interior nodes, and report total wall time including boundary/load generation, mesh construction, assembly, and LinearSolve. If the total time is comparable to or larger than direct AMFlow evaluation at the same interior points, the 'once and for all' practical claim is not supported.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central practical claim is that FEM solves Feynman integrals over a region 'once and for all,' but the benchmark measures only the discretization error of a manufactured problem with external analytic input. In Eq. (4), the Dirichlet data γ_n are stated to be 'input data.' In Eq. (11), the load vector F is assembled from the function f, which for each sector n depends on all lower-sector master integrals over the whole domain (Eq. (2)). In the numerical experiments, the paper explicitly says: 'The analytic solutions of these 12 MIs are taken from AsyInt and serve as the boundary conditions and as benchmarks.' Thus both γ_n and f are exact analytic inputs, not outputs of the proposed pipeline. The reported runtime is only LinearSolve time (Table I caption), excluding mesh generation, IBP reduction, assembly, and—crucially—the cost of computing boundary values and lower-sector integrals by the path-based methods the paper invokes. Without quantifying that input cost, the 14 s and 1.4e-4 average error do not demonstrate an end-to-end computational advantage over direct evaluation. This is load-bearing because the paper's motivation is to avoid path-based transport over the interior; if boundary/load generation is as expensive as direct evaluation, the method is essentially interpolation of someone else's computation. The summary statement 'requiring only a limited number of boundary conditions' understates the data requirement: every sector needs boundary values, and every sector except the first needs all lower-sector solutions over the whole domain to build f. This concern is acknowledged in the text but not addressed quantitatively, so the central claim's practical validity remains conditional.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper proposes to compute multi-loop Feynman integrals by reformulating the system of first-order differential equations as a second-order PDE, deriving a variational form, and solving it with the finite element method over a bounded two-dimensional phase-space region. The method is demonstrated on a massless two-loop non-planar four-point integral family with 12 master integrals. Using analytic results from AsyInt for boundary conditions and lower-sector master integrals, the authors report algebraic convergence of the O(ϵ^0) top-sector integral with average relative errors near 1e-4 and runtimes of up to 14.2 s for the sparse linear solve. The paper claims this is a new 'equilibrium' approach that solves the integral over a broad phase-space region 'once and for all', in contrast to one-dimensional path-based transport methods.","tokens_in":10308,"tokens_out":4059,"duration_ms":38263,"significance":"If the second-order PDE construction is correct and the required boundary/load data can be obtained cheaply, the method could complement existing path-based differential-equations approaches by providing a genuinely global solution over a region. The paper is a proof of concept with concrete numerical benchmarks, and the absence of fitted parameters is a strength: the Galerkin system is assembled directly from IBP-derived coefficients. However, the practical value hinges on two premises that the manuscript does not establish: the PDE system in Eq. (1) is asserted rather than derived, and the numerical experiments rely on exact analytic inputs for all boundary conditions and lower-sector integrals. The reported runtimes exclude the cost of generating those inputs, so the end-to-end computational advantage over direct evaluation is not demonstrated.","major_comments":[{"comment":"The construction of M(s,t,ϵ) is asserted without a derivation or an explicit example. The claim that 'the IBP relations required for constructing the second-order PDEs are identical to the ones for first-order DEs' needs justification: differentiating a first-order system introduces products of coefficient matrices and their derivatives, and the resulting M may be singular or lead to a non-elliptic operator. No ellipticity or well-posedness check is given. This is load-bearing because the Galerkin convergence observed in the benchmarks cannot be expected for indefinite or degenerate operators. Please provide the explicit derivation of Eq. (1) for the top-sector 2×2 system (or a small worked example) and discuss the properties of the operator -Δ - V over the chosen domain.","section":"Methodology, Eq. (1)"},{"comment":"The central claim of solving the integrals 'once and for all' is not tested end-to-end. The runtimes in Tables I and II are only LinearSolve times, while the paper states that 'the analytic solutions of these 12 MIs are taken from AsyInt and serve as the boundary conditions and as benchmarks.' Thus both γ_n in Eq. (4) and f in Eq. (2) are exact analytic inputs. The cost of obtaining these inputs via path-based transport or asymptotic methods is never quantified. Without that cost, the 14 s solve cannot be compared with direct evaluation, and the method may amount to interpolating someone else's computation. Please report the total wall-clock time including IBP reduction, mesh generation, assembly, and the numerical determination of boundary conditions and lower-sector integrals, or clearly state that the proof of concept only measures the discretization error of a manufactured problem.","section":"Numerical experiments, Tables I–II and captions"},{"comment":"Only the O(ϵ^0) term of the top-sector integral I_NPL is benchmarked (with the statement that I_NPL^(num) behaves similarly). However, the method solves 55 PDEs over five ϵ-orders, and the load functions for higher sectors depend on the FEM solutions of lower sectors. Since the lower-sector solutions are analytic inputs in this experiment, the error propagation through the sector hierarchy is not tested. Please report per-sector and per-ϵ-order errors for the full family, at least for the largest mesh, to support the claim that the whole system is solved reliably.","section":"Numerical experiments, error reporting"}],"minor_comments":[{"comment":"The decomposition In = u + u0 is confusingly written. 'u0|Ω/∂Ω = 0' is not a standard notation; please clarify that u0 is a finite-extension of the boundary data into the interior, e.g. by harmonic extension or by nodal interpolation on the boundary layer.","section":"Eq. (6)"},{"comment":"The phrase 'the absolute value is equivalent to the L2-norm for complex numbers' is imprecise; for a complex number the absolute value is the modulus, not an L2 norm.","section":"Page 3, error definition"},{"comment":"The statement 'requiring only a limited number of boundary conditions as input' understates the data requirement: every sector needs boundary values on ∂Ω and the load function f requires all lower-sector master integrals over the entire domain. Please qualify this sentence.","section":"Summary"},{"comment":"The assembly of the Galerkin matrices is described only in words. A short pseudo-code or a reference to a standard FEM assembly routine would improve reproducibility.","section":"Methodology, Eq. (9)–(11)"}],"recommendation":"major_revision","confidential_remarks":"The paper's idea is genuinely interesting and the numerical benchmark is clean as far as it goes, but the headline claim of a global method that solves integrals 'once and for all' is not yet supported because the input generation is delegated to existing tools without any cost or error analysis. I would like to see the following before reconsidering: a concrete derivation of M(s,t,ϵ) with a well-posedness discussion, and a comparison where boundary/load data are obtained numerically from a path-based method with total runtime reported. This can be done within the scope of a revised manuscript, so I recommend major revision rather than rejection."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"The two-line take: this paper introduces a genuinely new angle on multi-loop integrals — recast the first-order DE system as a second-order PDE equilibrium problem and solve it with finite elements over a bounded phase-space region. The two-loop non-planar benchmark is clean and convergent, with no fitted parameters. The soft spot is that the end-to-end claim is still untested: the boundary and load inputs are exact analytic data from AsyInt, and the reported 14 s is only LinearSolve time. So the practical advantage over path-based transport is not yet demonstrated.\n\nWhat is actually new: prior work solves first-order DEs along one-dimensional paths (Chebyshev spectral, AMFlow, DiffExp). This paper fomulates a global boundary-value problem on a bounded phase-space region and solves it with a standard Galerkin/FEM discretization in one pass. That framing is absent from the cited literature. The numerical experiments are honest — they deliberately use a known analytic integral to validate the method, report algebraic convergence rates, and show near-linear runtime scaling with mesh nodes. That is solid evidence the discretization works on this example.\n\nWhere the soft spots are: the derivation of M(s,t,ε) in Eq. (1) is asserted, not shown. A referee should ask for at least one explicit construction or a coercivity argument for the variational form; the numerical success helps but is not a substitute. The bigger issue is the benchmark's input cost. Both γ_n and the load f are built from exact analytic solutions from AsyInt, and only the linear solver is timed. The paper frames itself as avoiding interior transport, but if generating boundary and lower-sector data costs as much as direct evaluation, the method is essentially interpolating someone else's computation. That concern is real and not addressed quantitatively. It is not fatal — boundary data is lower-dimensional, so real savings may exist — but the paper should release code/meshes and run an end-to-end comparison that includes the cost of all inputs. Finally, max errors at ~1% at 75k nodes are fine for a proof of concept but far from FCC-ee precision; the suggested hp-FEM improvement is plausible but undemonstrated.\n\nWho is this for: people who do multi-loop numerics and care about phase-space interpolation methods. It deserves serious referee time — the idea is new and the proof-of-concept is real — but I would ask for a derivation check, an end-to-end cost benchmark, and artifact release before trusting it as a production method.","headline":"Genuinely new idea — FEM on a second-order PDE for Feynman integrals — but the benchmark uses exact analytic boundary/load data, so the practical claim is not yet tested.","tokens_in":10838,"tokens_out":4094,"would_cite":true,"duration_ms":36641,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["65N30","35J15","81Q30"],"pacs":[],"model":"deepseek-v4-flash","headline":"Multi-loop Feynman integrals can be computed as solutions of a second-order PDE in variational form, and the finite element method solves them over an entire phase-space region in one computation.","keywords":["Feynman integrals","second-order PDEs","finite element method","variational formulation","master integrals","integration-by-parts reduction","two-loop four-point integrals","phase-space evaluation"],"falsifier":"Run the identical 12-master-integral finite-element solve with boundary values produced by a path-based transport code instead of analytic formulas, then compare interior values against an independent high-precision evaluation. If the interior error is much larger than the reported 1.4×10⁻⁴ while the analytic-boundary solve stays accurate, the method's practical value collapses.","tokens_in":9878,"feed_emoji":"🧮","tokens_out":8433,"duration_ms":68238,"temperature":0.7,"pith_summary":"This paper aims to change how Feynman integrals are evaluated: instead of tracing each kinematic point along a one-dimensional path from known boundary values, it proposes to treat the whole region of phase space as a single equilibrium problem. Each master integral obeys a second-order PDE, and a finite-element discretization of the corresponding variational form turns the calculation into one sparse linear solve that returns the integrals everywhere in the region at once. The proof of concept is a two-loop non-planar four-point integral family; the top-sector integral is reproduced to average relative error 1.4×10⁻⁴ in 14 seconds on a single CPU core. If it holds up, this would give collider phenomenologists a fast, global alternative to point-by-point transport for precision predictions at future e⁺e⁻ machines. The open question the paper leaves is the cost of the boundary data, which currently come from analytic or path-based methods.","feed_headline":"Finite elements compute two-loop integrals in a single pass","feed_subtitle":"A variational second-order PDE turns phase-space evaluation into one sparse solve: 1.4e-4 average error on 75k nodes.","key_machinery":"The central object is the second-order PDE system ∇²I = M(s,t,ε)I for the vector of master integrals, rewritten in variational form as −∫∇Iₙ·∇v − ∫V Iₙ v = ∫f v, with prescribed boundary values γₙ on the boundary of the (s,t) domain. Here V is a potential and f is a load built from lower-sector master integrals; the same integration-by-parts identities used for first-order differential equations produce the system, so no new reduction cost is incurred. The workhorse is a piecewise-linear nodal-basis finite-element discretization: it assembles sparse matrices element by element and solves the resulting linear system once. The variational form only requires the solution and test functions to h","core_discovery":"The paper establishes that multi-loop Feynman integrals can be treated as the equilibrium solution of a second-order PDE over the whole (s,t) phase-space region, rather than as a one-dimensional transport problem. The evidence is a proof-of-concept: the massless two-loop non-planar four-point family is reduced to 12 master integrals, 55 PDEs after ε-expansion, and solved over a bounded region in one sparse linear solve. On a graded mesh with 75,008 interior nodes, the O(ε⁰) top-sector integral matches the known analytic result to average relative error 1.4×10⁻⁴ and maximal error 8.1×10⁻³, in about 14 seconds, with algebraic convergence rate around 0.8 in the average. The authors stress the e","pith_inferences":["The practical bottleneck the paper leaves implicit is boundary data: if accurate boundary values can only come from the same expensive analytic machinery the method is meant to bypass, the one-shot interior solve gains less than the headline numbers suggest. A quantitative cost/accuracy comparison of boundary input versus interior error would settle this.","An immediate testable extension is to feed the boundary data from the path-based transport methods cited in the paper rather than from analytic solutions, and compare interior errors on the same 75k-node mesh.","Since the benchmark uses a known analytic result, applying the same FEM machinery to a case where no analytic expression exists—for instance, a two-loop four-point integral with internal masses—would show whether the boundary-data requirement can be met in practice.","The same variational formulation might be used not just for master integrals but for the amplitudes themselves, if integration-by-parts reduction can be performed pointwise on the mesh; that would skip the separate assembly of master integrals."],"forward_implications":["If the convergence rate holds at higher loops, the method supplies the full (s,t)-dependence of master integrals from one linear solve, replacing per-point evaluations.","The second-order system is built from the same integration-by-parts reduction as first-order DEs, so no additional reduction work is needed to apply the method.","Higher-order finite elements and spectral bases are natural upgrades that the paper says could raise the algebraic convergence rate, potentially to exponential in smooth regions.","Massive thresholds, where integrals have kinks but remain continuous, stay in the square-integrable-derivative space and can be handled by partitioning the phase space or using complex masses.","Coupled master integrals are covered by off-diagonal potential blocks, as demonstrated for the two coupled top-sector integrals solved simultaneously."],"fun_headline_variants":["One sparse solve computes two-loop Feynman integrals","Galerkin solve yields two-loop phase-space integrals","PDE equilibrium method resolves two-loop four-point","Two-loop integrals: one PDE solve covers phase space"],"cache_read_input_tokens":2304,"weakest_assumption_plain":"The load-bearing premise is that accurate values of every master integral on the boundary of the phase-space region can be supplied at acceptable cost; the paper treats those boundary values as given input and does not measure their price or error.","fun_headline_variants_meta":{"raw":{"variants":["One sparse solve computes two-loop Feynman integrals","Galerkin solve yields two-loop phase-space integrals","PDE equilibrium method resolves two-loop four-point","Two-loop integrals: one PDE solve covers phase space"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.001294,"raw_usage":{"total_tokens":5067,"prompt_tokens":639,"completion_tokens":4428,"prompt_tokens_details":{"cached_tokens":256},"prompt_cache_hit_tokens":256,"prompt_cache_miss_tokens":383,"completion_tokens_details":{"reasoning_tokens":4375}},"tokens_in":383,"tokens_out":4428,"duration_ms":27231,"temperature":1.0,"reasoning_tokens":4375,"cache_read_input_tokens":256,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-01T15:29:33.809392+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run the identical 12-master-integral finite-element solve with boundary values produced by a path-based transport code instead of analytic formulas, then compare interior values against an independent high-precision evaluation. If the interior error is much larger than the reported 1.4×10⁻⁴ while the analytic-boundary solve stays accurate, the method's practical value collapses.","supporting_citations":[],"review_version":1}