{"id":"770344b7-169a-4eee-8589-3f2fc20bea1f","arxiv_id":"2411.17162","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":0,"one_line_summary":"A recursive tetrahedron refinement scheme computes Brillouin-zone integrals as weighted sums on the original grid, improving convergence for response functions with step, delta, and pole singularities.","lead":"This paper presents a recursive extension of the hybrid tetrahedron method for computing integrals over the Brillouin zone in crystalline solids. It expresses the integral as a weighted sum on the original momentum grid, enabling efficient and accurate treatment of singular integrands such as those in response functions.","discovery_kind":"extension","skeptic_critique":{"model":"deepseek-v4-flash","headline":"The central claim's practical reach depends on smooth, continuously labeled bands; standard sorted eigenvalues (and especially eigenvector-dependent matrix elements) are non-smooth at crossings, an issue the paper itself leaves unsettled in Section V.","rationale":"The reader's CONDITIONAL verdict is appropriate. The algorithmic core, Eq. (6) plus recursive weight collection, is internally coherent: if all input functions are smooth, quadratic interpolation and linear tetrahedron rules converge as refinement increases, and the code is available. The Lindhard test provides a solid validation for a single smooth band. The unresolved issue is the assumption of smooth, continuously labeled bands in multi-band applications. This is not a mere practical nuisance: the eigenvector matrix elements in the response functions are not even continuous at crossings, so the regular function F in Eq. (6) can be discontinuous along crossing manifolds. The authors knowingly flag the eigenvalue part in Section V, which is the reason the paper cannot claim general-purpose status without further work. Since the concern is about the scope of applicability rather than the internal validity of the algorithm, the verdict should remain CONDITIONAL: accept the method as a contribution, but require a demonstrated band-sorting strategy and an independent multi-band benchmark before it is advertised as a general tool. The reader identified the same root cause (band labeling), though I extend it to eigenvector-dependent numerators, so agreement is partial.","tokens_in":18367,"tokens_out":10989,"duration_ms":114688,"concrete_test":"Use the supplied BZIntegral.jl code on a two-band model with a genuine band crossing, e.g., H(k) = k_x sigma_x + k_y sigma_y + m sigma_z with m = 0 (Dirac point) or a two-band avoided-crossing model, and compute Im chi_+-(q, omega) from Eq. (15) using exact-diagonalization outputs sorted by energy. Compare the recursive method's results for nk = 8, 16 and nr = 0, 1, 2, 3 against a high-resolution reference obtained by direct summation on a 400 x 400 x 400 grid (or an analytic result where available). If the error does not decrease with nr at the expected rate, or if repeating with a smooth band-tracking labeling changes the result substantially, the band-sorting/eigenvector-continuity concern is confirmed.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The main result Eq. (6) approximates a BZ integral by a weighted sum over the initial grid, with the regular factor F(k) quadratically interpolated via the linear map Q and the singular weight W(k) built from interpolated eps(k) and D(k). This is valid only if the functions being interpolated are smooth enough for quadratic interpolation. In realistic multi-band electronic structure, band indices output by diagonalization are sorted by energy. The sorted eigenvalue functions are continuous but non-differentiable at band crossings, and the associated eigenvectors, which enter the numerator matrix elements in Eq. (15) and Eq. (18), are generically discontinuous. Quadratic interpolation over a tetrahedron containing such a crossing can overshoot or oscillate, corrupting both the approximation of F and the location of the step/delta/pole surfaces in W. Section V explicitly states that improper band labeling can cause uncertain discontinuities in band eigenvalue eps_nk and that higher-order interpolation usually causes overfitting in the case of discontinuous functions, concluding that an effective band sorting algorithm is needed for the method to work at its full potential. The fcc Co demonstration does not close this gap: the reference is the method's own high-resolution result, so convergence to that reference does not validate accuracy against an independent benchmark. Thus the paper's advertised application to realistic response-function calculations is conditional on a band-labeling solution that is neither provided nor tested.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"This paper proposes a recursive extension of the hybrid tetrahedron method for Brillouin-zone integration. The key result, Eq. (6), expresses a generic integral ∫ W(k)F(k) dk as a weighted sum of F on the initial k-grid, where the weights are obtained by propagating linear tetrahedron weights through the composite quadratic-interpolation map Q^(N,0). The weight function W(k) is assumed to depend on k through smooth functions (band eigenvalues ε(k), denominator D(k)) that can be quadratically interpolated; the method then handles step, delta, and 1/D singularities, including combinations such as Θ(ε_F−ε(k))/D(k). The Appendix provides closed-form linear tetrahedron weights for W=1/D with all limiting cases, as well as the standard step-function weights. Numerical demonstrations include the Lindhard function (compared to the exact result), the RPA transverse susceptibility of the honeycomb-lattice Hubbard model, and the Kohn-Sham susceptibility of fcc Co within TDDFT. A Julia implementation is released.","tokens_in":18576,"tokens_out":7586,"duration_ms":71446,"significance":"The central mathematical rearrangement in Eq. (6) is clean and correct, and the recursive weight-collection idea is a genuine practical improvement: it avoids storing exponentially many refined-grid values and makes the integral weights reusable for multiple integrands. The Lindhard-function test is a strong, falsifiable benchmark: the method converges systematically toward the exact analytic result with refinement, and the reported mean absolute error decreases rapidly with the number of refinements. The Appendix's closed-form weights for 1/D with all limiting cases is a useful reference for implementers, and the release of BZIntegral.jl is a concrete reproducibility asset. However, the paper's own Section V identifies an unresolved band-crossing problem that limits the method's practical reach to systems with smooth, continuously labeled bands; this must be addressed or clearly scoped in the manuscript before the broad claims in the abstract can be considered established.","major_comments":[{"comment":"Section V explicitly concedes that \"improper band labeling can cause uncertain discontinuities in band eigenvalue εnk\" and that \"an effective band sorting algorithm is needed for the method to work at its full potential.\" This is load-bearing for the advertised response-function applications, because Eqs. (15) and (18) mix band indices n,n′ and involve eigenvector-dependent matrix elements that are generically discontinuous at band crossings, while Table I and Eq. (6) rely on quadratic interpolation of smooth functions. The fcc Co benchmark in Fig. 5(c) is converged against the method's own 30×30×30, nr=1 result, so it does not independently validate accuracy in the presence of crossings. I request that the authors either supply a band-sorting algorithm, benchmark a crossing system against an independent exact or very-dense-grid reference, or explicitly restrict the claim of practical applicability to cases with continuous band labels.","section":"Section V and Eqs. (15), (18)"},{"comment":"The treatment of weight functions containing products of step functions, Eq. (8) and Eq. (10), applies identity (9) and then invokes the linear tetrahedron rule on Θ(−x1x2), with the justification that this is \"acceptable\" when the final tetrahedra are sufficiently tiny. The manuscript provides no error estimate and no numerical test for this case, despite the abstract's claim of \"simultaneously handling multiple singularities.\" I ask for a test of at least one product-of-step weight (e.g., a joint density of states or Eq. (10b)) or a clear statement that this part of the advertised functionality is heuristic and not benchmarked.","section":"Section III B, Eq. (9)"},{"comment":"The final paragraph states that the 1/D weights \"still hold after changing the arguments of logarithms from absolute values of D(k) to complex numbers D(k) itself,\" but the derivation preceding it treats real D with possible sign changes and absolute values. Since the response-function applications, Eqs. (13), (15), and (18), all use complex denominators with +iη, the analytic continuation of the closed-form weights to complex D should be justified, including the handling of branch cuts when D is near zero. Without this, a central ingredient of the flagship applications rests on an unproven assertion.","section":"Appendix VII C, final paragraph"}],"minor_comments":[{"comment":"Typographical errors should be corrected: \"respnse\" (Section IV, first paragraph), \"coloser\" (Section IV B), \"mehtod\" (Section IV B), and \"studys\" (Section IV B).","section":"Section IV"},{"comment":"The caption's \"log-log linear\" phrase should be clarified: the mean absolute error appears to scale as a power of Δk^3 for fixed nr, and the exponential decrease with nr is a separate claim; please specify the axes and the functional form being reported.","section":"Fig. 3(c,d) caption"},{"comment":"Eq. (12) introduces Q^(N,0)_ij with subscript order that is consistent with Eq. (6), but the index convention is not explicitly stated when the composite map is first defined in Eq. (4); please define the direction of the map (from refined to initial) at its introduction.","section":"Eq. (12)"},{"comment":"The limiting-case notation \"a&b → c\" and \"a&b&c → d\" is not self-explanatory; please define these as simultaneous limits (e.g., a→c and b→c) so that the formulas are unambiguous.","section":"Appendix VII C"}],"recommendation":"major_revision","confidential_remarks":"The manuscript is honest about its main limitation, which is commendable, but the abstract and introduction overstate the method's readiness for generic realistic response-function calculations. The fcc Co demonstration is weakened by the self-referential convergence check. The central derivation and the Lindhard benchmark are strong; I believe a major revision can bring the claims in line with the evidence."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Dear colleague,\n\nThe core of this paper is Eq. (6): a Brillouin-zone integral is rewritten as a weighted sum on the initial k-grid, with weights collected recursively through the composite quadratic-interpolation map Q^(N,0). That is a real algorithmic extension of MacDonald's hybrid tetrahedron method, and it removes the need to store an exponentially growing refined grid. The derivation is self-contained, the Appendix gives closed-form linear-tetrahedron weights for 1/D(k) with all limiting cases, and the tests against the exact Lindhard function show clear convergence improvement with each refinement. Credit where due: the authors also ship a Julia implementation and demonstrate the method on a honeycomb Hubbard model and fcc Co.\n\nThe soft spots are mostly ones the authors themselves flag. Section V states plainly that improper band labeling causes discontinuities in ε_nk and that higher-order interpolation overfits such discontinuities, so an effective band-sorting algorithm is needed for multi-band work. That is the main condition on the method's practical reach. The fcc Co test does not close this gap: the reference is the method's own high-resolution result, not an independent benchmark, so it shows self-consistency rather than external accuracy. The approximate treatment of composite step products via Eq. (9) is also justified only for sufficiently fine tetrahedra, which the paper acknowledges. Minor reproducibility gripe: the GitHub link lacks a commit hash and figure-reproduction scripts.\n\nNone of this undermines the central algorithmic claim. The method is a solid contribution for BZ integrals with singular weight functions, and the band-sorting problem is a known challenge in the field, not a hidden flaw. I'd send this to a competent referee. If it's accepted, I'd want the final version to either provide a working band-sorting strategy or state more carefully the class of systems for which the current implementation is reliable.\n\nFor you: worth reading if you work on tetrahedron methods or response-function codes. I would cite it.","headline":"A genuinely useful recursive weight-collection extension of the hybrid tetrahedron method, with the main caveat being the honestly flagged band-sorting problem.","tokens_in":19149,"tokens_out":2161,"would_cite":true,"duration_ms":20221,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":["71.15.-m"],"model":"deepseek-v4-flash","headline":"A recursive tetrahedron refinement scheme turns Brillouin-zone integrals into weighted sums on the original k-grid.","keywords":["Brillouin-zone integration","tetrahedron method","quadratic interpolation","weighted-sum integration","linear response functions","Lindhard function","Kohn-Sham susceptibility","singular integrands"],"falsifier":"Take a two-band model with a band crossing, compute a susceptibility with this method on a coarse grid using the code's sorted eigenvalues, and compare against a dense-grid reference computed with continuously labeled bands; if increasing the number of refinements does not shrink the error toward the dense-grid reference, the assumption of interpolable, continuous bands has failed.","tokens_in":18093,"feed_emoji":"📐","tokens_out":6821,"duration_ms":57042,"temperature":0.7,"pith_summary":"This paper proposes a recursive extension of the hybrid tetrahedron method that computes Brillouin-zone integrals as weighted sums over the original k-grid. The weights are generated by repeatedly dividing each quadratic tetrahedron into eight smaller tetrahedra using quadratic interpolation, then collecting linear-tetrahedron integration weights from the finest level back to the initial grid. The payoff is that refinement depth can be increased without storing the refined grid, and integrands carrying step, Dirac-delta, and pole singularities—the kind that appear in response and spectral functions—can be integrated in one framework. The authors demonstrate the convergence on the Lindhard function, the transverse spin susceptibility of the honeycomb-lattice Hubbard model, and the Kohn-Sham susceptibility of fcc cobalt.","feed_headline":"Recursive tetrahedron rule tames Brillouin-zone singular integrals","feed_subtitle":"Refinement adds accuracy while weights stay on the original k-grid, so response functions need no artificial broadening.","key_machinery":"The load-bearing object is the quadratic tetrahedron: a tetrahedron with values at its four vertices and six edge midpoints, which determines a unique quadratic interpolant. Dividing it into eight equal-volume subordinate quadratic tetrahedra and adding the new edge-midpoint values as weighted linear combinations of the parent values makes the refinement repeatable. Repeating this division and then applying analytic linear-tetrahedron integration on the finest grid, with all weights collected back through the linear maps $Q^{(p,p-1)}$, yields the main weighted-sum identity Eq. (6).","core_discovery":"On its own terms, the paper claims that the main identity, Eq. (6), holds: after any number of tetrahedron refinements, a Brillouin-zone integral can be expressed as $\\sum_m w_m F(k_m^{(0)})$, with the weights $w_m$ assembled by quadratic-interpolation maps from the finest grid back to the initial grid. The same recursive weight collection works when the singular weight function $W(k)$ contains a step function, a Dirac delta, a pole $1/D(k)$, or several of these at once, because the singular factors are evaluated on quadratically interpolated band eigenvalues and denominators and then treated by analytic linear-tetrahedron rules on the finest level. As a consequence, response and spectral functions, whose integrands contain denominators that can vanish, can be computed in the zero-broadening, zero-temperature limit rather than with a smearing parameter.","pith_inferences":["Because Eq. (6) only requires a linear relation between refined and initial grid values, the recursive weight-collection idea could in principle be applied with interpolation schemes other than the quadratic tetrahedron, such as adaptive or higher-order local interpolants.","The authors' flagged band-sorting problem suggests a natural paired extension: if eigenvalues are relabeled to be continuous across band crossings, the method should extend cleanly to degenerate multi-band response functions.","The appendix's two-dimensional triangle formulas make the same recursive treatment available for surface and 2D-material calculations, although the paper only demonstrates the 3D case.","Reusing the stored weights could make self-consistent linear-response loops much cheaper at equal accuracy, since the expensive interpolation and analytic-integration steps would not be repeated on every iteration."],"forward_implications":["With enough refinements, the method's accuracy approaches that of direct quadratic tetrahedron integration while using only the initial coarse k-grid values.","The same weighted-sum machinery handles step, delta, and pole singularities, so response functions can be computed without artificial broadening or smearing.","The integral weights depend only on the band structure and can be computed once and reused, for example inside a self-consistent density-functional perturbation theory loop.","In the demonstrated cases, refinements remove artificial peaks and jaggedness produced by broadening methods and reduce the overestimated magnon energies and Goldstone gap on coarse grids."],"supporting_citations":[{"why":"introduces the hybrid tetrahedron method that this paper extends to a recursive procedure","marker":"[9]"},{"why":"establishes that tetrahedron integrals can be written as weighted k-grid sums and supplies the closed-form step-function weights","marker":"[6]"},{"why":"provides the analytical linear-tetrahedron integration for a 1/D(k) denominator, adapted here to a symmetric form","marker":"[7]"},{"why":"defines the linear tetrahedron method for Fermi-surface step-function integrals that the recursive method uses at the finest level","marker":"[2, 3]"},{"why":"develops direct quadratic tetrahedron integration, the accuracy target that iterative refinement is claimed to approach","marker":"[4, 5]"},{"why":"is the prior response-function tetrahedron method with leveled linear approximants that the present pole-handling approach replaces","marker":"[8]"}],"fun_headline_variants":["Recursive tetrahedra enable zero-broadening singular Brillouin-zone integrals","Iterative tetrahedron refinement eliminates artificial broadening in k-space integrals","Recursive tetrahedron method computes singular integrals with zero smearing","Zero-broadening Brillouin-zone integrals from recursive tetrahedra"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The method assumes that band eigenvalues and denominator functions are smooth enough on the Brillouin zone that quadratic interpolation is faithful; the authors note that in practice eigenvalues from electronic-structure codes are sorted by size and become discontinuous at band crossings, and an effective band sorting algorithm is still needed.","fun_headline_variants_meta":{"raw":{"variants":["Recursive tetrahedra enable zero-broadening singular Brillouin-zone integrals","Iterative tetrahedron refinement eliminates artificial broadening in k-space integrals","Recursive tetrahedron method computes singular integrals with zero smearing","Zero-broadening Brillouin-zone integrals from recursive tetrahedra"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.00098,"raw_usage":{"total_tokens":4103,"prompt_tokens":831,"completion_tokens":3272,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":447,"completion_tokens_details":{"reasoning_tokens":3207}},"tokens_in":447,"tokens_out":3272,"duration_ms":20998,"temperature":1.0,"reasoning_tokens":3207,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-12T12:26:11.395174+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Take a two-band model with a band crossing, compute a susceptibility with this method on a coarse grid using the code's sorted eigenvalues, and compare against a dense-grid reference computed with continuously labeled bands; if increasing the number of refinements does not shrink the error toward the dense-grid reference, the assumption of interpolable, continuous bands has failed.","supporting_citations":[{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"introduces the hybrid tetrahedron method that this paper extends to a recursive procedure"},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"provides the analytical linear-tetrahedron integration for a 1/D(k) denominator, adapted here to a symmetric form"},{"cited_title":"[9] for readers’ convenience","cited_arxiv_id":null,"evidence_quote":"is the prior response-function tetrahedron method with leveled linear approximants that the present pole-handling approach replaces"}],"review_version":1}