Pith. sign in

REVIEW 4 major objections 4 minor 28 references

CHIC: Caley-Hamilton, Invariants and Constants for Neutrino Oscillation Probabilities and Gradients

T0 review · 4 major / 4 minor · reviewed 2026-08-03 · deepseek-v4-flash

Pith's one-line read This paper claims that the Cayley–Hamilton theorem lets three-flavor neutrino oscillation probabilities and their derivatives be computed analytically from the Hamiltonian's matrix invariants, without ever diagonalizing it, and that the res

desk verdict The analytic gradient results are a real addition, but the core equations as printed contain multiple sign errors that break the formulas exactly where the paper claims generality, and there is no numerical validation to fall back on. read the letter →

arxiv 2512.16427 v2 pith:NEDH66QJ submitted 2025-12-18 hep-ph hep-ex

classification hep-phhep-ex PACS 14.60.Pq
keywords neutrinooscillationsCayley-Hamiltontheoremmatrixinvariantsmattereffectsoscillationprobabilitygradientsclosed-formsolutionoscillogradsconstant-densitypropagation
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

Neutrino oscillation probabilities are usually expensive to compute because they require exponentiating a 3x3 Hamiltonian at every energy and baseline. The paper shows that by shifting out the trace and applying the Cayley–Hamilton theorem, the amplitude matrix becomes a quadratic polynomial in the reduced Hamiltonian, with coefficients determined by just two invariants and the baseline. This separates the energy dependence (which lives only in the eigenvalues) from the baseline dependence (which lives only in exponentials), so energy-only quantities can be cached when scanning distances. The same polynomial structure yields closed-form derivatives with respect to every mixing parameter, and the paper implements both in a C++ library that performs on par with the fastest known code, with a 30% speedup in baseline-scanning use cases. These gradients enable fast interpolation and a new visualization—oscillograds—of local parameter sensitivity.

What carries the argument

The Cayley–Hamilton theorem: a 3x3 matrix satisfies its own characteristic polynomial, so exp(-ixHtilde) modulo the cubic p(lambda)=lambda^3 + (Tr(Htilde^2)/2) lambda - det(Htilde) equals A I + B Htilde + C Htilde^2. The coefficients A, B, C come from matching at the three eigenvalues, which are found by Viète's trigonometric formula for the reduced cubic, using sqrt(54 det(Htilde)^2 / Tr(Htilde^2)^3). This identity carries the whole argument: it converts a matrix exponential into arithmetic with the two invariants and the eigenvalues, cleanly separating energy from baseline, and its derivative analog reduces every parameter gradient to products of a precomputed matrix S with the derivative

What would settle it

Pick a set of physical parameters (e.g., standard NuFit values with an electron density such that det(Htilde)<0), compute P with the paper's closed-form formulas and with exact numerical diagonalization of the same Hamiltonian, and compare all nine channels over a grid of energies and baselines; any discrepancy beyond floating-point error in P or dP/dzeta would show the unstated det>=0 restriction is not merely formal.

Watch

Extended reading notes

Core claim

The central claim is that for a constant-density three-flavor Hamiltonian, the amplitude exp(-ixH) can be written exactly as A I + B Htilde + C Htilde^2, where Htilde is the traceless part of H and A, B, C are functions of the invariants Tr(Htilde^2), det(Htilde), and the baseline x. Using Viète's trigonometric solution of the reduced cubic, the eigenvalues of Htilde are expressed in closed form, and the coefficients become scalar products of the eigenvalue-dependent vector eta_j = exp(-i lambda_j x)/(2 lambda_j^2 - pi_j) with the invariant vectors pi, lambda, and 1. This factorization shows that all energy dependence sits in the eigenvalues and all baseline dependence in the simple exponent

Load-bearing premise

The whole computation hinges on the formula lambda_0 = sqrt(2 Tr(Htilde^2)/3) cos(theta_Htilde/3) with theta_Htilde involving sqrt(54 det(Htilde)^2 / Tr(Htilde^2)^3); this is a valid eigenvalue only if det(Htilde) >= 0, and the paper never states or tests this restriction, so for negative determinants the polynomial remainder in Eq. (11) may not equal the matrix exponential.

Editorial extensions

If this is right

  • When many baselines are evaluated at one fixed energy, caching the energy-dependent eigenvalues and coefficients lowers computation time by about 30% compared with the current fast implementation.
  • Closed-form gradients with respect to any parameter (mixing angles, CP phase, density, mass splittings, energy) are available without numerical differentiation, and a six-point Hermite interpolation using them reproduces the appearance probability to better than 10^-4 relative error.
  • Gradient maps ('oscillograds') expose regions where the probability is nearly insensitive to a parameter, revealing structure like the residual CP/ordering degeneracy around delta_CP ~ pi in long-baseline channels.
  • The same polynomial decomposition applies to any Hamiltonian that can be shifted to traceless form and extended to beyond-Standard-Model additions, so the approach is not tied to the standard three-flavor case.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • The method's dependence on the sign of det(Htilde) is never stated; in negative-determinant regions the trigonometric root is not an eigenvalue and the polynomial remainder formula would fail. A straightforward numerical comparison with direct diagonalization over the physical parameter space would show whether this is a real limitation or a benign edge case.
  • The same separation of invariants might be imported into two-flavor or four-flavor (sterile neutrino) scenarios; for four flavors the quartic characteristic polynomial is harder but the invariant structure is analogous, and analytic gradients would probably still be affordable.
  • Because the gradients are analytic, a fitting framework could use them directly for gradient-based minimization or Hessian approximations, potentially removing the need for finite-difference scans over nuisance parameters in future analyses.
  • Oscillograds are introduced as a tool; a natural extension is to plot second derivatives (curvature) to map where a parameter is best constrained, which could be used in experiment design to choose beam energies and baselines.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

4 major / 4 minor

Summary. The paper proposes an analytic method, implemented in the CHIC software, for computing three-flavor neutrino oscillation probabilities and their derivatives in constant-density matter. The method uses the Cayley-Hamilton theorem to express the evolution operator as a second-order polynomial in the traceless Hamiltonian, with coefficients constructed from matrix invariants and eigenvalues obtained via Vieta's trigonometric solution of the reduced cubic. The authors claim that this avoids explicit diagonalization, separates energy and baseline dependences, and gives a fast C++/Python implementation competitive with NuFast. They also introduce 'oscillograds' as a visualization of probability gradients. The central claim is that the formulas deliver correct probabilities and analytic gradients over the full physical parameter space.

Significance. If correct, the proposed factorization into energy-only invariants and baseline-dependent exponentials is a useful contribution for long-baseline and reactor analyses, and the analytic gradients would be a genuine practical advantage. The manuscript is self-contained and does not rely on fitted parameters; the algebraic strategy is elegant and the timing comparison suggests a potentially competitive implementation. However, the correctness of the central claim is currently undermined by several sign errors in the printed equations, and the absence of accuracy comparisons against exact diagonalization leaves these errors undetected. The significance of the contribution can only be assessed after the algebraic issues are repaired and validated.

major comments (4)
  1. [III, Eq. (14)] The arccos argument is written as sqrt(54 det(H̃)^2 / Tr(H̃^2)^3), i.e. |3√6 det(H̃)|/Tr(H̃^2)^{3/2}. The Vieta trigonometric solution of the reduced cubic requires the signed quantity 3√6 det(H̃)/Tr(H̃^2)^{3/2}. For det(H̃)<0, which occurs for inverted ordering in vacuum and in parts of antineutrino/matter parameter space, the printed definition gives θ_H = arccos(|a|) instead of arccos(a) = π − arccos(|a|). The resulting λ0 from Eq. (13) is then not an eigenvalue of H̃; λ1 and λ2 in Eq. (17) inherit the error, and the coefficients in Eqs. (11), (19), (21) and the gradients in Eq. (25) are all wrong. This is a load-bearing error in the central claim, not a stylistic issue.
  2. [III, Eq. (11)] The Cayley-Hamilton reduction as printed uses the polynomial H̃^3 + (Tr(H̃^2)/2)H̃ − det(H̃)I. For a traceless 3×3 Hermitian matrix the characteristic polynomial is λ^3 − (Tr(H̃^2)/2)λ − det(H̃), with a minus sign on the linear term. The subsequent eigenvalue formula (13) corresponds to the negative-sign polynomial, so Eq. (11) is internally inconsistent. As written, Eq. (11) does not annihilate H̃; e.g. for H̃ = diag(1,0,−1) the printed polynomial gives H̃^3 + H̃ ≠ 0. The sign must be corrected and propagated through the coefficient functions.
  3. [III, Eqs. (19)-(20); IV, Eq. (23)] The denominator 2λ_j^2 − π_j in the definition of J_j (and in η_j, I_jj, and I_jj') is incorrect. Lagrange interpolation for the distinct eigenvalues requires the denominator 2λ_j^2 + π_j = ∏_{k≠j}(λ_j − λ_k). With the printed minus sign, the J_j do not form a partition of unity. For the traceless example H̃ = diag(1,0,−1), Eq. (19) gives Σ_j J_j = diag(e^{−ix}, −1, e^{ix}) instead of diag(e^{−ix}, 1, e^{ix}). This error affects the positive-determinant sector as well, and it propagates through Eq. (23) into the gradient formula (25).
  4. [V, Table I; VI, Figures 1-3] The only quantitative validation provided is the timing comparison with NuFast in Table I. There is no comparison of CHIC probabilities or derivatives against exact diagonalization or an independent code for either mass ordering, for neutrinos and antineutrinos, or across the relevant matter-density range. Such a validation is essential here: the sign errors above are exactly the kind of defect that a numerical cross-check would have exposed. Please add accuracy tests reporting maximum absolute differences in P and ∂ζP against exact results over the physical parameter space considered.
minor comments (4)
  1. [II, Eq. (1)] The third matrix in the PMNS product appears to contain a typo: the (2,1) entry is printed as −s13 c12 but should be −s12 c12.
  2. [Appendix A, Eq. (A6)] The diag argument '0, δ_CP m2_21, Δm^2_31' contains an apparent typo; it should be diag(0, Δm^2_21, Δm^2_31).
  3. [Title and throughout] The name 'Caley-Hamilton' is misspelled; it should be 'Cayley-Hamilton'.
  4. [III, Eqs. (17)-(19)] The formulas assume non-degenerate eigenvalues. The denominators 2λ_j^2 + π_j vanish when two eigenvalues coincide (e.g. at MSW resonances). The paper should state this assumption explicitly and, if possible, provide a limiting prescription for degenerate cases.

Circularity Check

0 steps flagged · score 1.0 of 10

No meaningful circularity: the core derivation is self-contained and based on standard Cayley–Hamilton algebra; the one self-citation is an ancillary application claim, not a load-bearing input.

full rationale

The paper's central derivation uses the Cayley–Hamilton theorem to express the matrix exponential as a quadratic polynomial in the traceless Hamiltonian, with coefficients determined by the eigenvalues. Those eigenvalues are obtained from the characteristic polynomial via Viète's trigonometric formula and from the invariants Tr(H̃²) and det(H̃), which are themselves computed from the Hamiltonian by Eqs. (15)–(16). This is a closed mathematical construction: no parameter is fitted to the probabilities or gradients that are later computed, and no 'prediction' is defined in terms of the quantity it claims to predict. The cited uses of Cayley–Hamilton in Refs. [16–18] are prior applications of the same theorem, not a uniqueness claim that forces the present result, and the paper does not rely on them to establish its formulas; it re-derives the relevant expressions explicitly. The only self-citation, Ref. [23], is used to motivate the usefulness of analytic gradients in profiling, which is an ancillary application argument and not part of the derivation of the oscillation probabilities or their derivatives. Even if one doubted that motivation, the mathematical results stand independently. The potential correctness issue with Eq. (14), namely that the square root of det(H̃)² discards the sign of det(H̃) and may give a wrong eigenvalue for physical parameter regions, is a real validation concern but is not circularity: it is a bug in the implementation of a standard theorem, not a reduction of the output to the input. The paper also provides no numerical comparison to exact diagonalization, but absence of validation is not circularity. Overall, the derivation chain is self-contained and the result is not forced by definition, by fitted inputs, or by self-citation.

Assumptions & free parameters 0 free parameters · 4 assumptions · 0 invented entities

The paper introduces no new physical entities or free parameters. The central derivation uses standard mathematical theorems and the standard neutrino oscillation framework. The main assumptions are the constant-density approximation and the use of the traceless reduced Hamiltonian.

assumptions (4)
  • standard math Cayley-Hamilton theorem: every matrix satisfies its characteristic polynomial.
    Used in Eq. (11) to express the matrix exponential as a second-order polynomial in the Hamiltonian.
  • standard math Vieta's trigonometric formula for the roots of a cubic.
    Used in Eq. (13) to obtain the eigenvalues of the traceless Hamiltonian; validity depends on the condition for three real roots.
  • domain assumption Three-flavor neutrino oscillation formalism with PMNS matrix and constant matter density.
    The entire framework in Section II assumes standard three-flavor mixing and constant electron density; varying-density cases are approximated from this solution.
  • domain assumption The Hamiltonian is traceless after shifting by Tr(H)/3; global phase cancels in probabilities.
    The reduction to a traceless matrix is required for the reduced cubic; the global phase in Eq. (10) is stated to cancel in probabilities.

how reviews work

0 comments
Cite this review

Pith. "Pith review of CHIC: Caley-Hamilton, Invariants and Constants for Neutrino Oscillation Probabilities and Gradients." pith.science (2026). https://pith.science/paper/NEDH66QJ

@misc{pith2026251216427,
  author       = {Pith},
  title        = {Pith review of: CHIC: Caley-Hamilton, Invariants and Constants for Neutrino Oscillation Probabilities and Gradients},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/NEDH66QJ}},
  note         = {Machine review of arXiv:2512.16427}
}
read the original abstract

We use the Caley-Hamilton theorem to derive analytical solutions for the three-flavor neutrino propagation amplitude in a constant-density medium and their derivatives with respect to the mixing parameters. This approach avoids the diagonalization of the Hamiltonian and exploits precomputed matrix invariants to separate the dependence of oscillation probabilities on neutrino energy and propagation baseline. The results are implemented in the CHIC software, which provides simple, fast and efficient computation of oscillation probabilities and their derivatives. Finally, we demonstrate the value of probability gradients for neutrino data analyses and introduce a complementary visualization, the oscillograds, to probe underlying features of neutrino mixing.

Figures

Figures reproduced from arXiv: 2512.16427 by the authors.

Figure 2
Figure 2. FIG. 2. Two-dimensional plots showing the transition prob [PITH_FULL_IMAGE:figures/full_fig_p005_2.png] view at source ↗
Figure 3
Figure 3. FIG. 3. Two-dimensional plots showing the transition prob [PITH_FULL_IMAGE:figures/full_fig_p006_3.png] view at source ↗

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

28 extracted references · 2 canonical work pages

  1. [1]

    Pontecorvo, Mesonium and anti-mesonium, Sov

    B. Pontecorvo, Mesonium and anti-mesonium, Sov. Phys. JETP6, 429 (1957)

  2. [2]

    Z. Maki, M. Nakagawa, and S. Sakata, Re- marks on the unified model of elementary par- 7 ticles, Progress of Theoretical Physics28, 870 (1962), https://academic.oup.com/ptp/article- pdf/28/5/870/5258750/28-5-870.pdf

  3. [3]

    Fukudaet al.(Super-Kamiokande), Evidence for os- cillation of atmospheric neutrinos, Phys

    Y. Fukudaet al.(Super-Kamiokande), Evidence for os- cillation of atmospheric neutrinos, Phys. Rev. Lett.81, 1562 (1998), arXiv:hep-ex/9807003

  4. [4]

    Q. R. Ahmadet al.(SNO Collaboration), Measurement of the Rate ofν e +d→p+p+e − Interactions Pro- duced by 8BSolar Neutrinos at the Sudbury Neutrino Observatory, PRL87, 071301 (2001)

  5. [5]

    Huber, J

    P. Huber, J. Kopp, M. Lindner, M. Rolinec, and W. Win- ter, New features in the simulation of neutrino oscillation experiments with GLoBES 3.0: General Long Baseline Experiment Simulator, Comput. Phys. Commun.177, 432 (2007), arXiv:hep-ph/0701187

  6. [6]

    R. A. Wendell, Prob3++ (2018), https://github.com/rogerwendell/Prob3plusplus

  7. [7]

    C. A. Arg¨ uelles, J. Salvado, and C. N. Weaver, nusquids: A toolbox for neutrino propagation, Computer Physics Communications277, 108346 (2022)

  8. [8]

    A. A. Abudet al.(DUNE Collaboration), Low exposure long-baseline neutrino oscillation sensitivity of the dune experiment, Phys. Rev. D105, 072006 (2022)

Show all 28 references
  1. [9]

    Collaboration, Sensitivity of the hyper-kamiokande experiment to neutrino oscillation parameters using ac- celeration neutrinos (2025), arXiv:2505.15019 [hep-ex]

    K. Collaboration, Sensitivity of the hyper-kamiokande experiment to neutrino oscillation parameters using ac- celeration neutrinos (2025), arXiv:2505.15019 [hep-ex]

  2. [10]

    T. J. Collaboration, Potential to identify the neutrino mass ordering with reactor antineutrinos in juno (2024), arXiv:2405.18008 [hep-ex]

  3. [11]

    Bustamante, Nuoscprobexact: a general-purpose code to compute exact two-flavor and three-flavor neutrino os- cillation probabilities (2019), arXiv:1904.12391 [hep-ph]

    M. Bustamante, Nuoscprobexact: a general-purpose code to compute exact two-flavor and three-flavor neutrino os- cillation probabilities (2019), arXiv:1904.12391 [hep-ph]

  4. [12]

    P. B. Denton and S. J. Parke, Fast and accurate algo- rithm for calculating long-baseline neutrino oscillation probabilities with matter effects, Physical Review D110, 10.1103/physrevd.110.073005 (2024)

  5. [13]

    Cayley, A Memoir on the Theory of Matrices, Philo- sophical Transactions of the Royal Society of London Se- ries I148, 17 (1858)

    A. Cayley, A Memoir on the Theory of Matrices, Philo- sophical Transactions of the Royal Society of London Se- ries I148, 17 (1858)

  6. [14]

    W. Hamilton, On the existence of a symbolic and bi- quadratic equation which is satisfied by the symbol of linear or distributive operation on a quaternion, The Lon- don, Edinburgh, and Dublin Philosophical Magazine and Journal of Science (1862)

  7. [15]

    Fern´ andez, Chic,https://github.com/pabloferm/ CHIC(2025)

    P. Fern´ andez, Chic,https://github.com/pabloferm/ CHIC(2025)

  8. [16]

    Ohlsson and H

    T. Ohlsson and H. Snellman, Three flavor neutrino oscil- lations in matter, Journal of Mathematical Physics41, 2768 (2000)

  9. [17]

    Gr¨ onroos, T

    J. Gr¨ onroos, T. Ohlsson, and S. Vihonen, Revisiting se- ries expansions of neutrino oscillation and decay proba- bilities in matter, Phys. Rev. D111, 035003 (2025)

  10. [18]

    D. S. Chattopadhyay, K. Chakraborty, A. Dighe, and S. Goswami, Oscillation and decay of neutrinos in mat- ter: An analytic treatment, Physical Sciences Forum8, 10.3390/psf2023008066 (2023)

  11. [19]

    R. P. Feynman, An operator calculus having applications in quantum electrodynamics, Phys. Rev.84, 108 (1951)

  12. [20]

    Guennebaud, B

    G. Guennebaud, B. Jacob,et al., Eigen, https://libeigen.gitlab.io (2010)

  13. [21]

    Jakob, J

    W. Jakob, J. Rhinelander, and D. Moldovan, pybind11 — seamless operability between c++11 and python (2016), https://github.com/pybind/pybind11

  14. [22]

    Esteban, M

    I. Esteban, M. C. Gonzalez-Garcia, M. Maltoni, I. Martinez-Soler, J. P. Pinheiro, and T. Schwetz, Nufit-6.0: updated global analysis of three-flavor neu- trino oscillations, Journal of High Energy Physics2024, 10.1007/jhep12(2024)216 (2024)

  15. [23]

    C. A. Arg¨ uelles, P. Fern´ andez, I. Mart ´ ınez-Soler, and M. Jin, Measuring oscillations with a million atmospheric neutrinos, Phys. Rev. X13, 041055 (2023)

  16. [24]

    M. Maltoni, From ray to spray: augmenting ampli- tudes and taming fast oscillations in fully numerical neutrino codes, Journal of High Energy Physics2023, 10.1007/jhep11(2023)033 (2023). Supplemental Materials – S1 Appendix A: Derivatives of the Hamiltonian

  17. [25]

    Electron Density ,ne ∂ne ˜H= √ 2G F ·diag (1,0,0)− √ 2G F 3 I(A1) 2.∆m 2 21 ∂∆m2 21 ˜H= 1 2E U·   0 0 0 0 1 0 0 0 0   ·U † − 1 6E I(A2) 3.∆m 2 31 ∂∆m2 31 ˜H= 1 2E U·   0 0 0 0 0 0 0 0 1   ·U † − 1 6E I(A3) 4.θ 23 ∂θ23 U=   U20 U21 U22 −U10 −U11 −U12 −U00 −U01 −U...

  18. [26]

    Neutrino energy ,E ∂E ˜H=− 1 2E2 U·diag(0,∆m 2 21,∆m 2

  19. [28]

    d dx Ψ =−i ˜HΨ,(A9)

    Baseline,x The case of the baseline is different and simpler than the rest and we can simply use Schr¨ ondiger’s equation in Equation 2 without the need of applying de results in Section IV. d dx Ψ =−i ˜HΨ,(A9)

  20. [31]

    + ∆m2 21 + ∆m2 31 6E2 (A8)

Pith tools

Reviewed August 3, 2026 · model on record in the stance chip above.