From ray to spray: augmenting amplitudes and taming fast oscillations in fully numerical neutrino codes

In this note we describe how to complement the neutrino evolution matrix calculated at a given energy and trajectory with additional information which allows to reliably extrapolate it to nearby energies or trajectories without repeating the full computation. Our method works for arbitrary matter density profiles, can be applied to any propagation model described by an Hamiltonian, and exactly guarantees the unitarity of the evolution matrix. As a straightforward application, we show how to enhance the calculation of the theoretical predictions for experimentally measured quantities, so that they remain accurate even in the presence of fast neutrino oscillations. Furthermore, the ability to “move around” a given energy and trajectory opens the door to precise interpolation of the oscillation amplitudes within a grid of tabulated values, with potential benefits for the computation speed of Monte-Carlo codes. We also provide a set of examples to illustrate the most prominent features of our approach.

A Exploiting the symmetries of the system 24 1 Motivation The discovery of neutrino oscillations have finally provided robust observational evidence that the Standard Model of particle physics is not the ultimate theory of nature.Lepton flavor conversion requires neutrinos to be massive, something which was not accounted for in the original formulation of the Standard Model.Adding neutrino masses through the usual Higgs mechanism is of course possible, but involves the introduction of right-handed neutrino states which, being gauge singlets, are not prevented by any known symmetry to acquire a Majorana mass.Hence, one way or another, New Physics seems to be at work in the neutrino sector, either through the appearance of something fundamentally new such as Majorana particles, or through some unknown mechanisms which prevents them.It is therefore understandable that during the last decades an intense experimental neutrino program has been carried out, and that even more powerful experiments are being developed for the coming years.
Neutrino experiments, by their very nature, aim at reconstructing neutrino properties by observing the effects of flavor conversion during propagation from source to detector.Such conversion depends of course on the assumed theoretical model (standard three neutrinos, extra sterile states, non-standard interactions, etc.) and on the specific value of its parameters.But it also depends on a set of "dynamical variables" characterizing the neutrino state (such as its energy E) or its trajectory (such as the path length L, or more generically the matter profile encountered along the path).In what follows we will sometimes refer to a specific instance of these variables (i.e., a concrete choice of energy, trajectory, etc.) as a "ray".Although in principle the experimental setup aims at determining the dynamical variables as accurately as possible, so to minimize their impact on the oscillation pattern and therefore extract the maximum information on the neutrino properties, in practice some amount of uncertainty is unavoidable.For example, the energy spectrum of the neutrinos emitted by the source is usually non-monochromatic, and the energy resolution of the detector is finite, so neutrino energy is never perfectly known.For atmospheric neutrinos [1] the imperfect reconstruction of the arrival direction implies that the traveled length and crossed matter profile are uncertain, furthermore the altitude of the production point is totally unobserved.Similarly, for solar neutrinos [2] it is impossible to determine at which point of the solar core the neutrino was produced, and therefore the exact profile of the traversed matter.All this implies that the comparison of experimental results with the theoretical expectations (be it the χ 2 of the "number of events" in a given data bin defined in terms of reconstructed quantities, or just the likelihood function for each individual event) always implies integrals (or averages) over dynamical variables such as the neutrino energy E or the path length L. Now, these integrals can be performed in many ways.Let us focus on the neutrino energy for definiteness, as it is ubiquitous to all experiments.One can resort to Riemann integration and divide the relevant range into a number of small bins, choose a representative value in each of them, calculate the corresponding conversion probabilities (for the given theoretical model and parameter values), and sum.Or one can use Monte-Carlo techniques, generating a random set of energy values according to some appropriate distribution.Either way, the conversion probabilities for each sampled point are calculated assuming a specific energy value «E», but are then used for the whole interval «dE» which such energy represent.This procedure implicitly assumes that the conversion probabilities will "stay the same" over the interval «dE», or at least that their variation can be reliably inferred from nearby sampled points without the need of further information.This is certainly true if the integration bins are "small enough", but sometimes this is prohibitively difficult to achieve: for example, when oscillations are "fast", which requires to choose very small intervals and therefore a very large number of integration points.
In brief, any point explicitly sampled from the integration domain (i.e., any "ray") actually represent a small but extended region (a "spray") around it.In view of this, when calculating neutrino propagation for a given ray one should make sure to recompile enough information to describe the conversion probabilities in the whole neighborhood it represents.In this note we present a method to address this issue in full generality, i.e., without assuming a specific oscillation model or matter density profile.We focus on scenarios where the evolution of the neutrino state is unitary and can be described in terms of an hermitian Hamiltonian by means of a Schrödinger equation.We do not consider here dissipative processes such as those commonly accounted by a Lindblad equation, nor neutrino selfinteractions which become relevant in dense media such as supernovas.A peculiar feature of our approach is that it is entirely formulated in terms of generic matrices which play a specific role in the description of neutrino propagation (such as the Hamiltonian H, the evolution matrix S ≡ T exp[−i H dx], etc.) but whose concrete expression as a function of the parameters of the model is never taken into account.This happens because we do not aim at providing analytic formulas valid for specific oscillation scenarios (which is already a widely studied topic in the literature, see refs.for an incomplete list), but rather at developing a model-independent framework which could serve as a guideline to enrich existing algorithms and extend their range of applicability "from the ray to the spray".This work is organized as follows.In chapter 2 we describe how to complement the conversion amplitudes calculated at a given energy and trajectory with extra information which allows to extrapolate them accurately to a neighborhood of such ray.In chapter 3 we apply this formalism to the case of "fast" neutrino oscillations, showing how the corresponding averaging effects can be implemented in fully numerical calculations.In chapter 4 we provide a set of examples to illustrate the advantages and limitations of the proposed approach, and in chapter 5 we summarize our conclusions.Finally, in appendix A we briefly discuss how specific symmetries of the neutrino system can be efficiently exploited within our formalism.

Formalism
The simplest and best known scenario accounting for leptonic flavor conversion consists in mass-induced neutrino oscillations in vacuum.From the phenomenological point of view, the fundamental properties of such model are: a) the evolution Hamiltonian H 0 is inversely proportional to the neutrino energy E.
This in particular implies that [H 0 (E 1 ), H 0 (E 2 )] = 0 even for E 1 ̸ = E 2 , so that it exists a basis in flavor space (the mass basis) where H 0 is diagonal for all energies; b) due to translational invariance of vacuum, H 0 is independent of the neutrino position in space, so that the evolution matrix S 0 ≡ exp(−iH 0 L) depends on the neutrino trajectory only through its total length L.
In this case the oscillation probabilities take a very simple form, essentially reducing to the sum of terms proportional to cos(γ i L/E) or sin(γ i L/E) where γ i generically denotes appropriate functions of the oscillation parameters.This simplicity allows for an analytic treatment of various oscillation effects which would otherwise be hard to handle in a fully numerical framework.For example, the GLoBES software [89,90] implements a feature called "low-pass filter" which averages probabilities and suppresses aliasing in the presence of very fast neutrino oscillations, for neutrino trajectories exhibiting translational invariance (i.e., vacuum and constant density).A similar functionality is also provided by the nuSQuIDS toolbox [91], again assuming a spatially uniform environment.In general, the two properties a) and b) above are responsible for the particularly simple dependence of S 0 on the neutrino energy and position, respectively.In the rest of this chapter we will show how, by means of suitable first-order Taylor expansions of the generic Hamiltonian H and the evolution matrix S, it is possible to attain such simplicity also in the general case of fully numerical neutrino propagation in an arbitrary matter profile, so that some analytic techniques commonly used for vacuum oscillations become applicable.

The constant-density case
Let's start by considering neutrinos propagating over a distance L in a constant matter potential, so that the evolution Hamiltonian H(E) does not depend on the position.It is reasonable to assume that H is a smooth function of E, so we can expand it at first order in a neighborhood of a reference energy Ē: With this, the matrix S(E) describing the neutrino propagation is: Now, it would be convenient to factorize S( Ē + ξ E ) into the product of two terms, one describing the evolution at the reference energy, S ≡ S( Ē), and the other accounting for the perturbation induced by the energy shift ξ E .In other words, we are seeking an expression of the kind: for a suitable matrix K E , which can formally be defined as: . (2.4) If the matrices H and H ′ E commute, as is the case in vacuum, it is immediate to see that K E is proportional to the derivative of the Hamiltonian, K E = H ′ E L. In the general case the expression for K E reads [92]: where we have introduced the unitary matrix Ū relating the flavor basis to the effective mass basis at reference energy Ē: In eq.(2.5) the operator ⊙ denotes the Hadamard product of two matrices, which consists in an element-wise multiplication of the corresponding elements: The matrix C is hermitian and its diagonal entries are equal to 1, while the non-diagonal entries are no larger than 1 in modulus: The matrix K E is also hermitian, which ensures that S( Ē + ξ E ) is exactly unitary for any value of ξ E .If [ H, H ′ E ] = 0 then it is possible to choose Ū so that H and H ′ E are simultaneously diagonalized, in which case the element-wise multiplication by C has no effect and K E = H ′ E L as previously stated.From the computational point of view, the most time-consuming step in the determination of K E is the diagonalization of H [93,94].However, such a step renders the computation of S (which requires to perform a matrix exponential) essentially trivial, so that the time spent in diagonalizing H is recovered from its exponentiation.We find therefore that computing also K E does not result in significant slow-down with respect to computing only S.

Multiple layers and arbitrary matter profile
In the previous section we have seen that, in the case of constant matter density, neutrino oscillations around a central energy Ē can be described in terms of two matrices: a unitary one S, describing the evolution for the specific value E = Ē, and an hermitian one K E , allowing to extrapolate S to nearby energies E = Ē + ξ E .Let us now consider the case in which the neutrino crosses two consecutive layers, each one with its own matter density.We will denote by ( S1 , K E 1 ) the evolution and perturbation matrices of the first layer, and by ( S2 , K E 2 ) those of the second layer.The combined evolution reads: The expression for S ≡ S( Ē) is readily obtained by setting ξ E = 0 in the previous formula, which yields S = S2 • S1 as expected.As for the combined K E , it can be found by means of eq.(2.4), which gives: (2.8) These two expressions can be summarized in a single "multiplication rule" among the pairs of ( S, K E ) matrices characterizing each layer: Interestingly, this product is associative, it has (1, 0) as identity element, and every pair ( S, K E ) has ( S † , − SK E S † ) as inverse, so that the set of pairs form a group.This is not a surprise, since it is clear from eq. (2.3) that the introduction of K E does not alter the unitarity (and therefore the algebraic structure) of the evolution matrix S(E).From eq. (2.9) it is evident that the extension to trajectories with varying matter potential of the formalism developed in the previous section follows the same line as the usual "fixed-energy" calculation of the evolution matrix S, except that this one is now replaced by a ( S, K E ) pair.Concretely, we proceed as follows: • we divide the trajectory into N smaller layers, in such a way that the variation of the matter potential within each of them is small compared to its average value; • we calculate the evolution pair ( Sn , K E n ) for all the n = 1, . . ., N layers, under the hypothesis of constant matter density; • we merge together all the individual pairs using the product defined in eq.(2.9), so that the overall evolution pair is ( S, As for the constant-density case, no significant amount of extra time is required to compute also K E as compared to computing only S, provided that the time for matrix multiplication is negligible with respect to that of matrix exponentiation.

Perturbations of the neutrino trajectory
So far we have only considered perturbations of the evolution matrix S around a reference energy Ē.However, sometimes the calculation of the event rates may require additional integration over other dynamical variables: for example, for atmospheric neutrinos the arrival direction (parametrized by the zenith angle Θ) plays a key role.In these cases, the same first-order expansion which we have just presented for the neutrino energy can be repeated for the other integration variables X, by considering the corresponding K X matrices.For instance, in the case of atmospheric neutrinos, for each reference energy Ē and zenith angle Θ we will write: Notice that, although strictly speaking the final result depends on the ordering in which we introduce the perturbation factors e −iK E ξ E and e −iK Θ ξ Θ (as in the general case the matrices K E and K Θ may not commute), the effect of interchanging them is of order ξ E ξ Θ , and can therefore be neglected in our first-order expansion.
For what concerns neutrinos crossing the Earth, the construction of K Θ is particularly simple.Given the spherical symmetry of the Earth, it is possible to approximate its density profile with a large number of constant-density shells.A given trajectory will cross a specific sequence n = 1, . . ., N of such shells, each one with a length L n (Θ) determined by the geometry.A large variation of Θ will cause shells to drop in or out of the reference path Θ (i.e., a change in N ), and this is a non-analytic effect which no Taylor expansion can reproduce.But for trajectories close enough to the central one, the sequence of shells (and therefore of the Hamiltonians Hn used within each of them) will not change, only the traveled length in each shell will be affected.Therefore: (2.11) where we have taken advantage of the fact that the "perturbation" commutes with the Hamiltonian Hn , so that no Hadamard product with a matrix C is required in this case.The ( Sn , K E n , K Θ n ) matrices for the various layers can then be merged together using the composition rule in eq.(2.9), trivially extended to accommodate also the Θ derivative: It is clear that further integration variables X which may affect neutrino propagation can be handled in the same way by simply appending their own perturbation matrices K X to the tuples of eq.(2.12).A particularly simple situation occur when X accounts for a longitudinal extension of the neutrino trajectory at one of its extremes, as implied by averaging over an extended production region or a non-negligible detector volume.In this case the matrix K X for the entire trajectory is directly related to the concrete value assumed by the evolution Hamiltonian at the relevant extreme, namely K X = H src at the source or K X = S † H det S at the detector.

Improved evolution within a definite layer
Till now we have discussed how to account for small deviations of dynamical variables (such as neutrino energy E or zenith angle Θ) from a central value used to perform the actual calculations.In this section we will instead concentrate on the reference ray itself (defined as E = Ē and Θ = Θ) and in particular on the construction of its evolution matrix (denoted as S in previous sections), showing how first-order Taylor expansion can be used to improve its computation as well.
As seen in section 2.2, a generic way to handle neutrino propagation in an arbitrary matter profile it to divide the trajectory into a number of layers, small enough so that the variation of the matter density within each of them can be considered small.We will focus here on one of these layers, assumed to have length L, and parametrize by x ∈ [0, L] the instantaneous neutrino position inside it.Let us denote by S(x) the unitary matrix describing the transition and survival amplitudes of the neutrino state from the beginning of the layer to position x, so that S(0) = 1 whereas S(L) corresponds to the evolution over the entire layer.The matrix S(x) satisfies the same equation as the state vector: By construction, the matter density varies little within the layer.A zero-order approximation is therefore to assume that it is perfectly constant, as we did in section 2.1.In this case H(x) = H and eq.( 2.13) can be solved immediately: where ω = diag{ω i } = Ū † H Ū as defined in eq.(2.6).For convenience we have denoted by S(x) the solution in the constant density approximation.The evolution matrix S(L) of the entire layer is: in full agreement with the formalism of section 2.1, and in particular with the definition of S appearing in eqs.(2.3) and (2.6).
The purpose of this section is to go beyond the constant-density approximation.To this aim, let us now introduce a small perturbation ∆(x) and decompose the Hamiltonian H(x) within our layer as: where H can be chosen as the value of H(x) at some specific location (such as x = 0 or x = L/2) or be defined by the condition ⟨∆⟩ ≡ L 0 ∆(x) dx L = 0. Following the approach of ref. [36], we seek the solution of eq.(2.13) in the form (2.17) with K(x) satisfying |K ab (x)| ≪ 1. Inserting eq.(2.17) into eq.(2.13) and keeping only the first order terms in ∆(x) and K(x), we find: At this point it is convenient to switch to the effective mass basis, so that H and S(x) become diagonal.Defining: we see that the differential equation for K(x) separates into individual components, and can therefore be solved by ordinary integration: This expression is just a special case of the general formalism presented in ref. [42], and corresponds to the truncation of the Magnus series to its first term.Specific derivations for concrete matter density profiles can be found in the literature, for example the Earth structure predicted by the PREM model [95] involves density shells which are well described by eq.(2.16).Accounting for the perturbation ∆(x) on top of the constant part H allows to compute the neutrino evolution inside each Earth shell in a single shot, without the need to further subdivide it into smaller layers.Examples of this approach can be found, e.g., in refs.[4,96] for solar neutrinos and in ref. [39] for atmospheric neutrinos.
In the present work, however, we do not want to stick to concrete matter density profiles, but we are interested instead in formulas which can be applied to generic situations.If the constant-density limit can be regarded as a zero-order approximation, then the natural first-order generalization is to assume that ∆(x) is a linear function of x within the given layer [40]: In such case the integral in eq. ( 2.20) can be computed analytically and we get: where sinc ′ (x) denotes the first derivative of the unnormalized sinc(x) function: Switching back to the flavor basis and using matrix notation, eq.(2.21) becomes: In principle the factor e i( φi − φj ) in eq.(2.21) could have been included into the definition of Ĉij , in which case the phase matrices e ±i φ in eq.(2.23) would not have appeared.In either case the matrix Ĉ is hermitian and with zero diagonal entries; but with the present choice Ĉ is also purely imaginary, which helps to speed up calculations when both H (and therefore Ū ) and H ′ are real matrices.In our convention the expression for the evolution matrix S(L) of the entire layer, including the first-order correction K, reads:

24)
Notice that the matrix K defined in eq.(2.23) plays a different role than the perturbation matrices K E and K Θ introduced in the previous sections.As detailed in eq.(2.10), (K E , K Θ ) describe how to alter the evolution matrix S when the dynamical variables (E, Θ) deviate from their reference values ( Ē, Θ) by finite amounts (ξ E , ξ Θ ).In turn, K encodes a correction to the constant-density approximation which is not controlled by any tunable quantity, and therefore there is no reason to keep it separated from S. For this reason, the correct way to implement K into the formalism developed in the previous sections is simply to replace S ≡ S(L) from eq. (2.15) with S ≡ S(L) from eq. (2.24) in the construction of the tuple ( S, K E , K Θ ) characterizing neutrino propagation.As for trajectories with multiple layers, one simply repeats this procedure for each layer and then combines them together using the multiplication rule of eq.(2.12).
3 Averaging In chapter 2 we have presented a formalism which allows to calculate the neutrino transition amplitudes in an extended neighborhood of a reference energy and trajectory.Here we will make use of these results to derive expressions for the flavor conversion probabilities, which are the key ingredient in the calculation of the theoretical predictions for experimentally measured quantities.In particular, we will show how our approach ensures that the integrals over the dynamical variables (such as the neutrino energy or trajectory) remain accurate even in the presence of fast neutrino oscillations, avoiding aliasing without the need to increase the density of integration points.The number of events observed by a neutrino experiment can be usually written as the sum over many oscillation channels (corresponding to initial flavor, final flavor, neutrino chirality, and so on) of expressions of the form: where P (E) is the neutrino conversion probability for the given oscillation channel, and N (E) denotes the "unoscillated number of events" which takes into account the neutrino flux at the source, the cross-section of the process, the efficiency and finite resolution of the detector, the number of targets and running time, and in general every factor or function which is required to properly describe the experimental setup.In principle, the integral in eq.(3.1) should extend over all the dynamical variables which affect neutrino propagation, such as the arrival direction or the production point for extended sources, but for definiteness we focus here only on the neutrino energy E.
In order to evaluate the integral numerically, it is useful to divide the integration domain into small intervals [E i , E i+1 ], so that N ch = i N i ch with: where ⟨ ⟩ denotes the average over the given bin, i.e., the integral itself divided by the bin's width ∆ E .In what follows we will assume that the function N (E) is relatively "smooth", in the sense that within each energy interval it is well approximated by a straight line: In turn, while we assume that the probability P (E) is continuous and differentiable, we do not require that it exhibits such slow variation, at least not so for every point in the parameter space.With this, eq.(3.2) becomes: ) and can therefore consistently be neglected with respect to the first one.It should be noted, however, that for probability functions which are slowly varying on the bin's energy range (so that their first-order expansion, P ( Ē + ξ E ) ≈ P + P ′ E ξ E , is a good approximation as we assumed for N ), then the second term in eq.(3.4) is of order O(∆ 3 E ) (because ⟨ξ E /∆ E ⟩ = 0, so that the leading P contribution vanishes).This suggests that keeping only the first term in the expansion, E ) approximation at least in some case.We will return on this later on.

Average over energy
As we have just seen, calculating the integral in eq.(3.1) requires estimating the average value of both functions N (E) and P (E) over the range of each energy interval.The former is pretty easy, as under the assumption that the first-order expansion is accurate enough within the bin, we can simply use the value of N in the central point of the bin: N = N ( Ē).Alternatively, if N does not depend on the parameters of the model (as it is the case when the physics model under consideration only affect neutrino propagation), we can afford estimating N = ⟨N ⟩ numerically by subdividing the bin into smaller parts, as in any case this is a one-time-only calculation.
In order to calculate the average probability ⟨P ⟩, we take advantage of eq.(2.3).Let α and β denote the initial and final neutrino flavor state, so that P (E) ≡ |S βα (E)| 2 .Then: where we have introduced the matrix V E diagonalizing K E : Expanding in components: Averaging P (E) over the bin's energy range reduces to calculating e i(λ It is interesting to notice that for ∆ E = 0 this expression immediately reduces to ⟨P ⟩ = P ( Ē), so a numerical code which implements averaging as described here can also trivially provide unaveraged results.Such situation also arises when the eigenvalues of K E are "small", (λ E j − λ E i )∆ E ≪ 1, in which case the oscillation probabilities vary slowly over the bin's energy range.In turn, in the limit of very fast oscillations, (λ , so interference among different (i ̸ = j) effective mass eigenstates is suppressed leading to full decoherence.
For completeness, we also provide the expression of the higher-order term ⟨P • ξ E /∆ E ⟩ which appears in eq.(3.4): Notice that sinc ′ (x) ∼ −x/3 for x → 0, as can be seen in eq.(2.22), so in the limit is suppressed by one power of ∆ E , making the second term of eq.(3.4) of order O(∆ 3 E ) -as already inferred in the introduction of this chapter.However, if (λ E j − λ E i )∆ E ∼ 1 such suppression does not take place, and the corresponding correction -although still subleading withe respect to the ⟨P ⟩ contributionis simply of order O(∆ 2 E ).For simplicity, in the rest of this note we will neglect this term.

Average over trajectory
The generalization of these results to integrals over multiple dynamical variables follows the same line.Let us consider the case of atmospheric neutrinos described in section 2.3.Now in addition to the neutrino energy E we should integrate also over the zenith angle Θ, and the expression of S(E, Θ) is given by eq.(2.10).In order to calculate the average probability ⟨P ⟩, we first diagonalize the perturbation matrices K E and K Θ : With this, denoting by α and β the initial and final neutrino flavor state, we can write: which, after averaging over their respective bin intervals ∆ E and ∆ Θ , yield: This expression, albeit correct, is not very illuminating.Things become clearer if we make use instead of the following algorithm, which reproduce eq.(3.12) by applying a chain of transformations (ρ 0 → ρ 1 → ρ 2 → ρ) to the density matrix describing the neutrino state: a) we begin by setting the density matrix to the projector onto the initial neutrino state: b) we rotate it to the basis where K Θ is diagonal, multiply it element-wise by a matrix G Θ , and rotate it back to the flavor basis: c) we rotate it to the basis where K E is diagonal, multiply it element-wise by a matrix G E , and rotate it back to the flavor basis: d ) we apply the evolution operator S, thus obtaining the density matrix ρ at the detector: The average probability is then given by ⟨P ⟩ = Tr ρ Π (β) = ρ ββ .In any case, it should be noticed that this approach involves the construction of the entire neutrino density matrix, so that it is readily at hand in situations when the bare probabilities do not suffice (for example, in the presence of flavor-changing neutrino interactions in the detector, as is the case for NSI with electrons [97]).The sequential approach described above makes it manifest the way averaging acts.Each matrix K X associated to a dynamical variable X gets decomposed into two parts: its eigenvalues, which induce decoherence by suppressing the off-diagonal elements of the density matrix, and its diagonalizing matrix V X , which determines in which basis the aforesaid suppression takes place.Averaging over different variables (E and Θ in our example) results in subsequent decoherence applied in different basis.As already noted at the beginning of section 2.3, the order in which we average over E and Θ affects the final result, but only at subleading order ∆ E ∆ Θ .
From the mathematical point of view, decoherence is introduced through element-wise multiplication of the density matrix ρ by a suitable matrix G X .This process does not spoil the hermiticity of ρ since G X is itself hermitian.Furthermore, the condition Tr(ρ) = 1 is unaltered as the diagonal entries of G X are identically 1 by construction.Finally, the property Tr(ρ 2 ) ≤ 1 is preserved since |G X ij | ≤ 1 for any i ̸ = j pair.

Integral over production point
Sometimes the neutrino source has a sizable spatial extension, so that the integration over the production region cannot be neglected.Such integral can be formally decomposed into two components: the longitudinal one, corresponding to the direction of neutrino propagation, and the transversal one, which is orthogonal to it.Integration over the transversal variables X is just a special realization of the "average over trajectory" discussed in the previous section, and can therefore be described in terms of suitable perturbation matrices K X and smearing matrices G X .The same is certainly true also for the longitudinal integral, with the extra benefit that the corresponding changes to the trajectory only affect a little portion at its beginning and are therefore straightforward to implement.On the other hand, the integral over the longitudinal direction becomes even simpler when neutrino propagation exhibits translational invariance inside the production region, as stated in property b) at the start of chapter 2. A typical example is provided by atmospheric neutrinos, for which the oscillation probabilities depend not only on the energy and zenith angle, but also on the altitude of the production point in the atmosphere.Usually the air matter density can safely be neglected, so that propagation proceeds as in vacuum and is described by an Hamiltonian H 0 (be it the usual vacuum term, or a different one if New Physics is present) independent of the position.In the rest of this section we will focus on this particular case, using atmospheric neutrinos as a guideline.
Let us begin by fixing the neutrino energy and zenith angle to reference values Ē and Θ, and neglecting their variation at first.In this case, denoting by ℓ the slant height of the production point (i.e., the distance to the ground level as measured along the neutrino trajectory, which is not necessarily vertical), we have: where ξ ℓ is the distance from a reference position l, which may or may not coincide with ground level.The first thing to notice is the formal similarity of this expression with eq. ( 2.3), the main conceptual difference being that eq.(3.17) is exact for any ξ ℓ and not the outcome of a first-order expansion.Letting U 0 be the matrix which diagonalizes H 0 , so that U † 0 H 0 U 0 = ω 0 = diag{ω 0 i }, we can write: The next step is to average over the altitude of the production point.Denoting by π ℓ (ξ ℓ ) the probability density of creating a neutrino at slant height ( l + ξ ℓ ), we get: Hence, thanks to the assumption of translational invariance of H 0 , the integral over the neutrino production point can be performed in a single shot, without the need of splitting the integration domain into smaller steps.Furthermore, eq.(3.19) is impressively similar to the formalism presented in the previous sections, which suggests that its numerical implementation can be easily merged with the average over the neutrino energy and direction.Indeed, this is accomplished by modifying the algorithm in section 2.3 as follows: a-c) we proceed as before until the construction of the density matrix ρ 2 ; d ) we rotate it to the basis where H 0 is diagonal, multiply it element-wise by a matrix G ℓ , and rotate it back to the flavor basis: e) we apply the evolution operator S, thus obtaining the density matrix ρ at the detector: As can be seen, this approach treats averaging over neutrino energy, arrival direction and production altitude on the same footing.It should be noted, however, that our formalism is based on a first-order expansion, and therefore relies on the assumption that the three dynamical variables ξ E , ξ Θ and ξ ℓ are sufficiently small.While ξ E and ξ Θ can be kept under control by suitably choosing their corresponding intervals ∆ E and ∆ Θ , the range of ξ ℓ is determined by the properties of the Earth's atmosphere (or more in general of the neutrino source), and cannot be changed.Yet this does not spoil the accuracy of the calculations when ∆ E = ∆ Θ = 0 as long as translational invariance ensures that eq.(3.17) is exact.
In other words, our procedure may fail to account for terms of order ξ E ξ ℓ or ξ Θ ξ ℓ when deviating from the reference ray, but while the smallness of ξ ℓ is not a priori guaranteed for arbitrary sources, such terms are still subleading due to the smallness of ξ E and ξ Θ .This issue can be further mitigated by tuning the reference altitude l, for example ensuring that the mean of π ℓ (ξ ℓ ) is zero.Anyway, the size of the Earth's atmosphere is indeed small compared to the overall radius of the Earth, hence in this case the validity of the computation is secured by the physical system.And of course, for non-uniform sources (such as the core of the Sun) or when the interplay between the source's overall extension and the size of the energy and zenith integration bins cannot be neglected, one always have the option of splitting the production range into steps small enough to circumvent these issues, and handle the longitudinal integral numerically.
Looking at the definition of G ℓ in eq.(3.20), we see that its elements are strictly related to the Fourier transform πℓ of the altitude distribution function π ℓ : G ℓ ij ≡ πℓ (ω 0 i − ω 0 j ).This is also the case for G E or G Θ , since they were constructed assuming uniform priors within their respective ranges: π E (ξ E ) ≡ rect(ξ E /∆ E )/∆ E and similarly for π Θ (ξ Θ ), whose Fourier transform is indeed the sinc(x) function.This suggests that the flat averaging over the bin's range which we have performed so far can be generalized by assuming alternative distributions for ξ E and ξ Θ .For example, a Gaussian prior for π E (ξ E ) would yield: This is precisely the idea behind the low-pass filter in ref. [90], and can be useful to describe, e.g., the smearing of the oscillation probabilities induced by a finite energy resolution ∆ E of the detector -provided that ∆ E is small enough for our first-order expansion to hold.Alternatively, in Monte-Carlo calculations where the dynamical variables E And Θ are chosen randomly by an integrator routine and no energy or angular grid are defined, it may be convenient to introduce exponential smearing on scales ∆ E and ∆ Θ well below the resolution of the detector, so to properly handle fast oscillations without spoiling the reliability of the simulation.

Tabulation and interpolation
To conclude this chapter, let's briefly comment on a trivial extension of the techniques described so far.In eq. ( 2.10) we have illustrated how the perturbation matrices (K E , K Θ ) can be used to "shift" the evolution matrix S from its central value S calculated at ( Ē, Θ) to a nearby position ( Ē + ξ E , Θ + ξ Θ ).In the previous sections we have used this formula to derive accurate averages over energy and zenith angle, assuming some distribution π E (ξ E ) and π Θ (ξ Θ ) (either plain rectangular functions with widths ∆ E and ∆ Θ , or more general ones such as Gaussian priors) around the central value ( Ē, Θ).However, our formalism trivially allows to perform averages also around shifted values, ( Ē + δ E , Θ + δ Θ ).This is accomplished by means of shifted priors, π E (ξ E − δ E ) and π Θ (ξ Θ − δ Θ ), which leads to a rephasing of the G E and G Θ matrices: This simple observation opens the door to efficient tabulation and interpolation of oscillation amplitudes.Consider the case where a Monte-Carlo generator needs to compute the neutrino conversion probabilities for a very large number of (E, Θ) rays.A well-known technique to speed up computations is to first tabulate the probabilities on a representative grid of ( Ēi , Θj ) values, and then extract the actual (E, Θ) ray by interpolation.The problem in doing so, however, is that in the presence of fast oscillations a fixed ( Ēi , Θj ) grid may fail to reproduce the oscillation pattern accurately enough.The solution is to tabulate instead the ( S, K E , K Θ ) matrices (which for further convenience can also be factorized at this stage into unitary ( Ū , V E , V Θ ) and diagonal ( ω, λ E , λ Θ ) components) for each ( Ēi , Θj ) node, and later use this information to reconstruct the probabilities once the required (E, Θ) value is known.The most straightforward way to perform this last step is to find the closest ( Ēi , Θj ) node and use it for extrapolation.A more refined approach is to locate the [ Ēi , Ēi+1 ] × [ Θj , Θj+1 ] cell containing (E, Θ), derive an estimate of the conversion probabilities from each of its vertices, and then produce a weighted average of such estimates as in ordinary interpolation.In this second case we can also evaluate the reliability of the result by comparing the probabilities obtained from the various vertices, as for accurate calculations they should all be similar among them.

Examples
In this chapter we will present a number of examples to illustrate the main features of the formalism just introduced.Concretely, we will focus on three aspects: Taylor expansion in energy and trajectory (described in sections 2.1, 2.2 and 2.3), improving the accuracy of S within a definite layer (described in section 2.4), and averaging in the presence of fast oscillations (described in chapter 3).

Taylor expansion in energy and trajectory
In figures 1 and 2 we plot the oscillation probabilities in various channels for atmospheric neutrinos (solid lines) and antineutrinos (dashed lines) crossing the Earth matter.We assume standard three-neutrino oscillations and set the corresponding parameters to the NuFIT-5.2best-fit value [98,99].We fix Ē = 0.3 GeV (figure 1) or Ē = 3 GeV (figure 2) as reference value for the neutrino energy, as well as cos Θ = −0.9 for the zenith angle of the arrival direction, and compute the matrices ( S, K E , K Θ ) defined in chapter 2. We then plot the dependence of the probabilities on the neutrino energy E = Ē + ξ E (left panels) and zenith angle Θ = Θ + ξ Θ (right panels), and compare the exact calculations (thick colored lines) with the extrapolation based on eq.(2.10) (thin black lines).
As can be seen, all black lines in all panels match the value of their colored counterpart at zero shift.This is by construction, as that corresponds precisely to the reference value used for the calculation of S ≡ S( Ē, Θ).However, in addition to point-like coincidence the black lines are also tangent to the colored ones, and this is a consequence of taking into account also the first-order terms.To make it clear, if we had neglected K E and K Θ in our calculations -thus sticking just to the usual zero-order approximation -all the black lines would have been perfectly horizontal.
Another feature of our correction terms is that they capture the relevant oscillation "frequencies" (in E and Θ) of the system, so that the black lines can often "track" the exact calculations for a sizable interval around zero.This is particularly the case at low energy (see figure 1), when oscillations are dominated by the vacuum term for which our formalism becomes exact.Notice, however, that even at Ē = 0.3 GeV matter effects still play an important role, as the clear differences between same-channel neutrino and antineutrino probabilities demonstrate, yet this does not spoil the accuracy of the extrapolation.
of the non-commutativity of the system, so they cannot help in estimating the "curvature" of the lines (which is a second-order effect) beyond the simple periodic oscillation pattern.This is clearly visible in the red, orange and cyan neutrino lines, where our extrapolation sizably deviates from the exact result already for energy shifts at the few-percent level.This underlines that our method is not intended for large-scale extrapolations, and in particular the size of the integration bin should be kept small enough to ensure that any non-oscillatory effect is properly accounted for numerically.In other words, the formalism described here takes care of potentially fast oscillations stemming from large derivatives of the evolution Hamiltonian, but the features of the slow oscillation pattern still require a dense grid in the (E, Θ) plane.

Improved in-layer calculation
Neutrino propagation in arbitrary matter profiles can be handled by dividing the path into a number of sufficiently small layers.As described in section 2.4, within each layer we can either assume plain constant density, or add a correction proportional to the first derivative of the matter potential.To illustrate the advantage of the second choice, in figure 3 we plot the P ee survival probability for a neutrino produced in the center of the Sun and detected at infinite distance.For definiteness we assume two-neutrino oscillations with sin 2 θ = 0.3 and ∆m 2 = 7.4 × 10 −5 eV 2 , and the solar matter distribution and chemical composition given in ref. [2].For such a model the MSW effect takes place [100,101] and neutrino probabilities can be calculated analytically using the adiabatic approximation.This result is an excellent benchmark to check the accuracy of our formalism, therefore we have plotted it in both panels as a thick dashed orange line.
The solid lines in figure 3 have been computed in a fully numerical way.Concretely, we have divided the trajectory inside the Sun into as many layers as indicated in the legend, and we have obtained the overall evolution matrix S for the full path by multiplying together Figure 3. Asymptotic survival probability for a neutrino produced in the center of the Sun.We assume a 2ν oscillation model with sin 2 θ = 0.3 and ∆m 2 = 7.4 × 10 −5 eV 2 .In both panels the orange dashed line is obtained with the adiabatic formula.The colored solid lines correspond to fully numerical calculations, with the trajectory inside the Sun divided into as many layers as indicated in the legend.In the left panel we assume that the matter density is constant within each layer, while in the right panel we account for the first-order correction described in section 2.4.
the contribution of the various layers.In the left panel we have assumed constant density within each layer, as described in eq.(2.15).As can be seen, in order for the numerical calculation to reproduce the analytic result accurately enough over the relevant energy range, one need at least O(10 4 ) layers.Qualitatively, this can be understood as follows.
A requirement for MSW conversion is that the vacuum and matter term of the evolution Hamiltonian become comparable at some point along the trajectory.These terms can be conveniently quantified through the oscillation length they induce, l osc = 2π/(ω 2 − ω1 ) where ωi are the eigenvalues of the Hamiltonian (see eq. (2.6)).For ∆m 2 = 7.4 × 10 −5 eV 2 we get l vac osc = 33 km • E/MeV in vacuum, while for E → ∞ we have l mat osc ≥ 160 km in solar matter.The condition l vac osc ∼ l mat osc requires E ≳ few MeV (and indeed the transition between the vacuum-dominated and matter-dominated regime occur in this range, as clearly visible in figure 3) and implies oscillation lengths l osc ≳ O(100 km).The numerical computation is accurate when the layer size does not exceed the oscillation length, and this is ensured in all the MSW region only if the layer length is smaller than O(100 km): hence, the number of layers should be at least O(R ⊙ /100 km) ≃ O(10 4 ), with R ⊙ = 7 × 10 5 km being the solar radius.Empirically, one may say that the break of adiabaticity induced by the "jumps" in the potential among consecutive constant-density layers should occur at scales well below the oscillation length.In the left panels we fix sin 2 θ 23 = 0.572 and plot ∆χ 2 as a function of ∆m 2  31 for various accelerator neutrino experiments.In the right panel we show the allowed regions (at 1σ, 90%, 2σ, 99%, 3σ CL for 2 d.o.f.) in the (θ 23 , ∆m 2 3ℓ ) plane from the global analysis of Super-Kamiokande atmospheric data.The colored lines or regions are based on the averaging procedure described in chapter 3, while the black lines do not.See text for details.
On the other hand, in the right panel of figure 3 we have taken into account the linear variation of the matter potential inside each layer, as encoded in eq.(2.24).This effectively removes the artificial "jumps" introduced by the ladder-like schematization of the potential in the constant-density limit, and leads to impressively accurate results with as little as a few tens of layers.It should be remembered, however, that our formula for the evolution in a linearly-varying potential is not exact (unlike the constant-density case) but rather obtained through a perturbative expansion, so the layer length should always be kept sufficiently small for the approximation to hold.

Averaging fast oscillations
In chapter 3 we have shown how the perturbation matrix K E (and K Θ for atmospheric neutrinos) can be used to improve the calculation of the energy integral (and zenith-angle one) commonly required to estimate the theoretical prediction of a given measurement.In particular, our approach naturally handles fast neutrino oscillations, yielding properly av-eraged results without the need of ad-hoc solutions.To illustrate this feature, in figure 4 we consider various experiments and compare the fits obtained with and without the inclusion of K E and K Θ .For definiteness we assume standard 3ν oscillations and fix the undisplayed parameters to the NuFIT-5.2best-fit value [98,99].
In the left panels we focus on accelerator experiments MINOS [102], NOvA [103] and T2K [104] and plot the overall ∆χ 2 (defined with respect to the local minimum) as a function of ∆m 2  31 .The energy integral for each experiment is converted into a sum by subdividing the relevant range into uniform bins in logarithmic scale, with a density of 100 bins per decade.The central point of each energy bin is chosen as representative value for the entire bin and used to calculate neutrino propagation, encoded in the evolution matrix S.This simple integration method is straightforward to implement, and produce accurate results until well beyond the boundaries of the experimentally allowed region.However, for the sake of illustration we are interested here in the domain ∆m 2  31 ≳ 10 −1 eV 2 , even though it is completely ruled out by the data.In this limit oscillations become so fast that the conversion probabilities can no longer be regarded as "constant" within an energy bin.If this fact is ignored and the probabilities are still extracted solely from the bin's representative S matrix, as is the case for the black lines in figure 4, then the calculation becomes unreliable due to aliasing effects.Conversely, if the perturbation matrix K E is also taken into account as described in eqs.(3.6) and (3.8), then fast oscillations are automatically averaged and accuracy is recovered.As a further example, in the right panel of figure 4 we plot the allowed region in the (θ 23 , ∆m 2 3ℓ ) plane from our own global analysis of Super-Kamiokande atmospheric data [105].In this case the energy integral is estimated with a density of 50 points per decade in logarithmic scale, while the neutrino arrival direction is discretized into 100 points uniformly distributed in cos Θ ∈ [−1, +1].Despite the very large number of sampled points (more than 30 000 rays in the full (E, Θ) plane) the results obtained with calculations based exclusively on S are inaccurate, as illustrated by the black lines.This is driven by sub-GeV data, for which E ≲ 1 GeV so that the oscillation probabilities of neutrinos coming from below the horizon are "fast" (i.e., they vary a lot even for small energy and zenith changes).Once again, taking into account the information encoded in the perturbation matrices K E and K Θ fixes the issue, as can be deduced from the colored regions.
Of course, one of the reasons behind the failure of calculations based solely on S is that we have chosen a regular grid of sampling points in both energy and zenith angle, which favors aliasing effects: randomizing our grid would have mitigated the problem.Also, for atmospheric neutrinos we have verified that doubling the density of points (100 per decade in energy, and 200 overall in zenith) significantly improves the quality of the fit (at least in the standard 3ν case), at the cost of a factor ∼ 4 in computer time.In general, various methods exists to handle fast oscillations, but they all have some kind of drawback.For example, one can use an adaptive integration routine which "detects" poor accuracy and adds extra points to compensate for it, but this usually implies substantial slow-down in difficult regions.Alternatively, one can introduce a "low-pass filter" as described in ref. [90], but this requires to choose an "averaging length" according to the details of the experiment under consideration, furthermore its implementation is only feasible in limited situations (e.g., in ref. [90] this option is provided just for constant density).Integration over a spatially uniform production region can be performed analytically, as implemented in ref. [91] and discussed in detail in section 3.3, but it is only effective for averaging purposes when the oscillation length is smaller than the source's size, which is not very often unless the source occupies a considerable fraction of the overall baseline.Finally, in specific scenarios where fast oscillations are known to occur (such as in solar neutrinos due to large ∆m 2  31 , in atmospheric neutrinos at sub-GeV energy, or in long-baseline experiments when extra eV-scale sterile states are considered), it may be possible to factor them out analytically while leaving the treatment of non-fast oscillations numerical (see, e.g., appendices C and D of ref. [106]) or semi-analytical [78], but this approach requires the derivation of appropriate formulas for each propagation model, and it relies on oscillation frequencies being "infinitely large" so that it cannot handle smooth transitions between "fast" and "slow" oscillations.In contrast, our method generically works for any oscillation model which can be described in terms of an evolution Hamiltonian, requires a fixed amount of computation time irrespective of the specific point in parameter space being simulated (i.e., it is not affected by whether fast oscillations arise or not), and does not require the choice of an "averaging length" because the finite extension of the integration bins take care of that (and the result is independent of it, as long as the bin is small enough for the first-order approximation to hold).

Summary
In this note we have presented a general formalism which allows to considerably enhance the accuracy and performance of numerical neutrino codes needed to calculate the theoretical predictions for experimentally measured quantities.In particular: • our approach does not make any assumption on the underlying theory determining neutrino propagation, hence it can be applied to a vast set of models such as standard three-neutrino oscillations, extra sterile neutrinos, non-standard neutrino-matter interaction, violation of fundamental symmetries, and so on.Furthermore, it works for arbitrary matter density profiles; • our method relies on a first-order Taylor expansion of the neutrino evolution matrix S(E, Θ) around a reference energy Ē and trajectory Θ.As described in eq.(2.10), S( Ē + ξ E , Θ + ξ Θ ) is related to its central value S ≡ S( Ē, Θ) through suitable perturbation matrices (K E , K Θ ), and its unitarity is guaranteed for any value of (ξ E , ξ Θ ).The set of ( S, K E , K Θ ) tuples with the multiplication rule in eq.(2.12) forms a group, and provides the building blocks to compute neutrino propagation on trajectories comprising multiple density layers; • our formalism ensures that the integrals over neutrino energy and trajectory embedded in the theoretical predictions of experimental measurements remain accurate even in the presence of fast neutrino oscillations, avoiding aliasing without the need to increase the density of integration bins.This is achieved through element-wise multiplication of the neutrino density matrix with smearing matrices (G E , G Θ ) in a way entirely controlled by the perturbation matrices (K E , K Θ ), as seen in eqs.(3.14) and (3.15).In Riemann integration the bin's width naturally acts as a low-pass filter, yet its specific value does not affect the final result as long as it is small enough for the first-order approximation to hold.In Monte-Carlo simulations a suitable cutoff can be introduced by appropriate priors π E (ξ E ) and π Θ (ξ Θ ) such as Gaussian functions; • our method also allows for efficient tabulation and interpolation of oscillation amplitudes.The ( S, K E , K Θ ) matrices can be pre-computed on a representative grid of ( Ēi , Θj ) values, which are then used to reconstruct the conversion probabilities once the required (E, Θ) ray is known.This procedure avoids the loss of accuracy and aliasing effects which usually appear in the presence of fast oscillations when the probabilities themselves are tabulated and interpolated; • finally, for atmospheric neutrinos (and in general extended sources) our approach unifies averaging over neutrino energy and direction with the integral over the neutrino production point, also described by a suitable smearing matrix G ℓ as shown in eq.(3.20).Furthermore, it naturally leads to the construction of the density matrix at the detector, which is convenient when considering scenarios where the plain oscillation probabilities do not suffice.
On the technical side, a pre-existing object-oriented code accounting for neutrino propagation solely in terms of S can be adapted to incorporate our formalism by replacing the matrix S and its product with the tuple ( S, K E , K Θ ) and the multiplication rule in eq.(2.12).The addition of the first-order terms does not result in significant slow-down, as long as the computation time required for matrix exponentiation is comparable to that of matrix diagonalization and overwhelms that of matrix multiplication.Concretely, one extra diagonalization for each K E or K Θ matrix is required to perform averaging, which results in doubling the computation time for constant-density paths (such as accelerator neutrinos) but negligible impact on trajectories with a large number of different layers (such as atmospheric neutrinos).In brief, our approach provides a lossless enhancement with respect to computations based on S alone, which can thus be regarded as its zero-order limit.

A Exploiting the symmetries of the system
Let us consider a system described by an Hamiltonian H(E) ≡ O H(E) O † where O is a unitary matrix.This situation occur, for example, in standard three-neutrino oscillations, where the matrix O accounts for the θ 23 and δ CP parameters while the "reduced" Hamiltonian H(E) depends solely on θ 12 , θ 13 , ∆m 2 21 , ∆m 2 31 .In this case, it is immediate to see that S(E) ≡ O S(E) O † for all energies, so that one can perform the bulk of calculations in the so-called "propagation basis" using the reduced matrix H(E) (which is simpler and depends on less parameters) and then reintroduce O at the end.The formalism developed in this note is completely transparent with respect to this factorization.In particular: The same happens when multiple derivatives (such as K E and K Θ ) or the altitude of the production point for atmospheric neutrinos are considered.
As for the actual averaging, the algorithm described in section 3.3 is trivially modified in its first and last step to incorporate O: a) the initial matrix ρ 0 must be rotated to the propagation basis, so that: In summary, models where a group of parameters can be factorized out and reintroduced at the end are perfectly compatible with our formalism.Finally, we want to comment on a rather common feature of many oscillation models which is sometimes exploited to speed-up computations.It is not infrequent that the oscillation probabilities depend on the neutrino energy only through a particular combination of it with the parameters of the model (here collectively denoted as ⃗ Ω), so that P is invariant under a simultaneous rescaling of the energy E → αE and suitable transformation of the parameters ⃗ Ω → ⃗ Ω α .This is the case, for example, in standard three-neutrino oscillations, where the probabilities depend on the energy and the mass-squared differences through the combined ratios ∆m 2 ij /E.Such situation allows to reuse the probability spectrum P (E) tabulated for a given point ⃗ Ω in parameter space, for all the other points related to it by the transformation ⃗ Ω → ⃗ Ω α , as these would require a simple "shift" E → αE of the tabulated energy values.When energy averaging is introduced into the game, one should be careful not to spoil the invariance of the system.Concretely, denoting by ∆[ Ē] the range of the bin with central energy Ē, we should make sure that ∆[α Ē] = α ∆[ Ē].This is trivially achieved if ∆[ Ē] Ē = constant, i.e., a uniform spacing in log(E) is used to define the energy grid.

Figure 2 .
Figure 2. Same as figure 1 but for Ē = 3 GeV.In the left (right) panel we show the effects of modifying the energy (direction).See text for details.
α) ij = δ αi δ αj ; (A.1) b-d ) remain the same as before, expressed in terms of the matrices "without the O"; e) the matrix O is reintroduced in the construction of the final density matrix ρ: ρ 3 → ρ ≡ O S ρ 3 S † O † .(A.2) is no larger than 1/2 in absolute value, so the second term in eq.(3.4) is at most of order O(∆ 2 E • the full Hamiltonian H(E) can be decomposed asH( Ē + ξ E ) ≈ H + H′ E ξ E with H ≡ O H O † and H′ E ≡ O H ′ E O † .Similarly, S ≡ O S O † ; • the Hamiltonian H is diagonalized by a matrix Ů ≡ O Ū , and its eigenvalues are the same as H. Hence, Ů † H Ů = Ū † H Ū = ω.Consequently, the matrix C used to construct K E is unaffected by O, and we get KE = O K E O † as expected; • in brief, the pair ( S, KE ) accounting for the full Hamiltonian H(E) is related to the reduced one ( S, K E ) by an overall rotation of each individual matrix: ( S, KE ) ≡ O ( S, K E