Model-independent test of T violation in neutrino oscillations

We propose a method to establish time reversal symmetry violation at future neutrino oscillation experiments in a largely model-independent way. We introduce a general parametrization of flavour transition probabilities which holds under weak assumptions and covers a large class of new-physics scenarios. This can be used to search for the presence of T-odd components in the transition probabilities by comparing data at different baselines but at the same neutrino energies. We show that this test can be performed already with experiments at three different baselines and might be feasible with experiments under preparation/consideration.

Introduction. The violation of time reversal (T) and charge-parity (CP) symmetries are central topics in particle physics. CP violation (CPV) is one of the necessary conditions to generate a matter-antimatter asymmetry in the early Universe [1], and under the well founded assumption of CPT conservation, CPV is equivalent to T violation (TV). A particularly active field is the search for CPV in neutrino oscillations [2][3][4]. Unfortunately, the experimental signature is rather indirect, and it is not possible to construct model-independent CP-asymmetric observables in neutrino oscillation experiments. This is related to the fundamental obstacle that experiments and detectors are made out of matter (and not antimatter). Moreover, the passage of the neutrino beam through Earth matter introduces environmental CPV due to matter effects [5].
The standard approach to this problem is to perform a model-dependent fit to data. This involves the assumptions that neutrino production, detection and propagation is fully understood in terms of Standard Model (SM) interactions, that neutrino mixing is unitary, and only the three SM neutrino flavours exist. In this case oscillation physics can be parametrized in terms of a unitary 3 × 3 lepton-mixing matrix [6,7] and two neutrino masssquared differences. CPV is then described by a complex phase δ in the mixing matrix [2,8] which can be fitted against data. "Observation of CPV" is considered equivalent to establishing that δ is different from 0 and π at a certain confidence level. Within this restricted framework, current data start to provide first indications of preferred regions for the parameter δ [9][10][11][12][13].
Large activity is devoted to study the impact of nonstandard scenarios on the search for CPV in neutrino oscillations. Examples are non-unitary mixing [14,15], non-standard neutrino interactions [16][17][18], or the presence of sterile neutrinos [19][20][21]. In such new-physics scenarios, additional complex phases appear, which can act as new sources for CP and T violation. Typically one adopts a specific parameterization of new-physics and again performs a parametric fit in the extended model. Our aim in this letter is to go a step beyond such approaches and develop a largely model-independent test, covering a wide class of non-standard scenarios. Our approach is based on fundamental principles about TV noted in the seminal paper by N. Cabibbo [2].
Model-independent description of flavour evolution. Here we specify our approach to describe the neutrino survival and transition probabilities P αβ , with α, β = e, µ, τ . P αβ is the probability for a neutrino ν α produced at the neutrino source to arrive as ν β at the detector. We adopt the following assumptions: (i) Propagation of the three SM neutrino states is described by a hermitian Hamiltonian H(E, x), which depends on neutrino energy E and in general on the matter density at the position x along the neutrino path.
(ii) We assume that for the experiments of interest, medium effects can be described to sufficient accuracy by a constant matter density which is approximately the same for all considered experiments. This is a good approximation for experiments with baselines less than several 1000 km [22,23]. In an accompanying paper [PRD] we show that for the baselines relevant for our test, the effect of non-constant density is negligible.
Assumption (ii) implies that the matter effect does not introduce environmental TV by itself [24,25] [PRD]. Therefore, any observation of TV can be related to fundamental TV of the theory. (General discussions about TV in neutrino oscillations can be found e.g. in [2,[26][27][28][29][30][31][32][33][34].) Furthermore, the Hamiltonian becomes position independent and we can diagonalize it as H(E) = W λW † , with W being a unitary matrix and λ = (λ i ) is a diagonal matrix of the real eigenvalues λ 1 , λ 2 , λ 3 of H. Both W and λ depend on the neutrino energy. Note that we allow for arbitrary non-standard matter effects. In general, W and λ will be different for neutrinos and antineutrinos.
(iii) We allow for arbitrary (non-unitary) mixing of the energy eigenstates ν i with the flavour states ν α relevant for detection and production, We make no specific assumption on the complex coefficients N αi . In particular, we do not relate them to the unitary matrix W , we allow them to be arbitrary (suf-ficiently smooth) functions of energy, and they can be different for neutrino production and detection. But we do assume that they are the same for different experiments (at the same energy). Below we focus on the experimentally relevant ν µ → ν µ disappearance and ν µ → ν e appearance channels. Under the assumptions (i), (ii), (iii) the corresponding probabilities are obtained as [PRD] with and ω ij ≡ λ j − λ i . Similar expressions are obtained in the context of non-unitary mixing, e.g., [35][36][37]. As usual we have traded the time dependence in the evolution equation with the space coordinate, leading to the appearance of the baseline L in the above probabilities. Therefore, T reversal is formally equivalent to L → −L [2,25]. The first line of Eq. (3) is invariant under T, whereas the second line is T-odd. Fundamental TV can be established by proving the presence of the L-odd term in the probability.
For practical reasons we will introduce one more assumption. We list it here for completeness and provide further discussion and motivation below: (iv) We impose that the oscillation frequencies ω ij deviate only weakly from the ones corresponding to the standard three-flavour oscillation case.
Note that our assumptions (i), (ii), (iii) are rather general and cover a large class of new physics scenarios, including non-standard interactions in production, detection, and propagation [38], generic non-unitarity [35][36][37], or the presence of sterile neutrinos, as long as their associated oscillation frequencies are large compared to ω ij [39].
The TV test. The strategy we propose to probe the L-odd terms is to measure the oscillation probability as a function of L at a fixed energy and check whether Leven terms are enough to describe the data or if TV is required. Under these conditions, the effective frequencies and mixings in the Hamiltonian are the same, and so the data at different baselines (but at the same energy) can be consistently combined. Notice that antineutrino data cannot be analyzed together with neutrino data, as their effective frequencies and mixings are in general different from the neutrino's; nevertheless, separate tests could be made for neutrino and antineutrino data.
In the absence of TV, all c α i are real and the data points could be described by the L-even part of the oscillation probability. We define (c α i real) For the two relevant channels, these probabilities depend on 8 parameters, which we collectively denote by θ: 6 real coefficients c µ i , c e i (i = 1, 2, 3) and two independent ω ij , e.g., ω 21 and ω 31 . We assume now that the probabilities P µµ and P µe are measured at a fixed energy at several baselines L b . We denote the corresponding measured values by p dis b and p app b with the uncertainties σ dis b and σ app b , respectively. Below we are going to assume that p dis b and p app b correspond to the values predicted by standard three-flavour neutrino (3ν) oscillations in matter.
We now ask the question if we can exclude the hypothesis of T conservation parametrised by Eq. (4), if the data correspond to 3ν oscillations with TV, i.e., for a CP phase δ different from 0 or π. To this aim we construct the χ 2 function The best-fit T-conserving model is obtained by considering χ 2 min (E) = min θ χ 2 even (E; θ) . We will take the value of χ 2 min (E) as a rough indication of how strongly T conservation can be excluded by data, and leave a more detailed statistical analysis for future work. Considering that each baseline provides 2 data points (appearance and disappearance) and that the T-even model has 8 parameters, it is clear that we need more than 4 experiments at different baselines. Let us note, however, that our parameterization includes so-called zero-distance effects, due to the non-unitary mixing in Eq. (1). Therefore, the near-detector(s) of long-baseline experiments provide already two data points at L ≈ 0 and effectively only more than 3 experiments are needed.
This requirement can even be further relaxed if we impose one additional assumption, which can be motivated by the fact that we have overwhelming evidence that the standard three-flavour scenario is approximately correct and any new physics effect can only be subleading. Therefore, we introduce assumption (iv) mentioned above: we assume that the oscillation frequency ω 21 deviates only weakly from the one corresponding to the standard 3ν case. Technically we impose this requirement by calculating the effective mass-squared difference in matter, ∆m 2 21 (E), assuming the standard matter effect and add the following prior to Eq. (5): Estimated number of appearance signal events at future accelerator experiments, assuming normal mass ordering and true δ = 90 • . Data from Refs. [47,48] (DUNE), [43] (T2HK), [44] (T2HKK), and [46] (ESSνSB).
The frequencies ω ij are determined by the effective evolution Hamiltonian, assumption (i), and are independent of the mixing in Eq. (1). A general parameterization of the Hamiltonian is provided by the non-standard neutrino interaction scenario, see e.g., Ref. [40]. Using the results of Ref. [40] we estimate the possible deviation to σ 21 = 0.1∆m 2 21 , see the appendix for more details. It turns out that the other independent frequency, ω 31 , is effectively constrained by the long-baseline data used in our fit, and therefore it is not necessary to impose an analogous prior for it. The prior in Eq. (6) acts as an additional data point for each energy bin (note that also the prior is energy dependent). Therefore, under this additional assumption, we come to the remarkable result that our model-independent test can be performed already with 3 experiments at different baselines plus near detectors.
The crucial requirement, however, is sufficient overlap in neutrino energy. If experiments have overlapping energy ranges, we can combine information from different energies. However, to be completely modelindependent, the minimization has to be done individually for each energy, since we do not want to make any assumptions about the energy dependence of the unknown new physics. This is an important difference to usual model-dependent analyses.
Realistic baselines and energies. Let us now consider planned long-baseline accelerator experiments in order to see if such a test realistically can be carried out in the future. We consider the following experiments: the DUNE project in USA (L = 1300 km) [41,42], T2HK in Japan (L = 295 km) [43], with the option of a second detector in Korea, T2HKK (L = 1100 km, 1.5 • off axis) [44], and a long-baseline experiment at the European Spalation Source in Sweden, ESSνSB (L = 540 km) [45,46].
Expected event numbers are obtained from Design Re- ports or detailed studies of the physics potential and are shown for the appearance channel in the case of 3ν oscillations and δ = 90 • in Fig. 1. In practice, we will see that only the two energy bins between 0.7 and 0.9 GeV provide relevant sensitivity, as data points with sufficient statistics are needed at 1st and 2nd oscillation maxima. We note that the energy spectrum from the NOνA experiment [10] has no overlap with the T2K beam and therefore it cannot be used for this analysis. We use the information from Fig. 1 (and the corresponding data for the disappearance channel) to estimate the statistical uncertainties in Eq. (5) as σ br /P even (L b , E r ) = √ S br + B br /S br at baseline b and energy bin r. We take the background events B br directly from the experimental studies and estimate the number of signal events from the N br in the Figure assuming S br = N br × P even (L b , E r ; θ)/P 3ν (L b , E r ). For the near detector data points, we assume the standard P αβ (L → 0) = δ αβ with σ = 0.01.
In Fig. 2 we show the data points for the appearance and disappearance probabilities as a function of the baseline for the 0.7-0.8 GeV energy bin. We can see that the disappearance data points essentially fix the oscillation frequency, whereas the appearance data are crucial for the TV test. The "true" oscillation probability assumed to generate the data points correspond to standard 3ν oscillations with maximal TV (δ = 90 • ) and normal mass ordering. We find that no satisfactory L-even fit is possible for the 4L and 3L (HKK) combinations at this energy. The essential information is obtained from the relative heights of the first and second appearance oscillation peaks, see the appendix for further discussion. Note that disappearance probabilities can reach values larger than one in our fit, since we do not impose unitarity in our effective parameterization of the T-even transitions.
In order to connect our test with experiments, one should take into account the fact that finite energy resolution effectively changes the L-dependence in their measurements, which will in turn affect the sensitivity of the TV test. We assume a given energy resolution ∆E around the central bin energy E 0 , and smear the transition probability by convoluting it with a Gaussian with mean E 0 and width ∆E. To illustrate the effect we assume here ∆E = 0.1E 0 . In order to perform the convolution one must assume a certain energy dependence of the transition probability. Our assumption is that the energy dependence of the amplitudes c α i is slow enough, such that it can be neglected within an interval of few ∆E. The only significant energy dependence would thus be in the oscillation phases ω ij . According to assumption (iv) introduced above, we assume that ω 31 ∝ 1/E, as in the standard 3ν oscillation case. We have checked that our results are independent of energy smearing of ω 21 terms. The impact of the finite energy resolution is illustrated in the right panels of Figure 2.
Our results for maximal TV are summarized in Tab. I, which shows the χ 2 min values for the various energy bins for different experiment combinations, with and without including the energy smearing. We observe that 0.75 GeV is the most relevant energy bin, whereas the one at 0.85 GeV still provides some sensitivity. The strong impact of the energy resolution is apparent. We also find that the detector in Korea is essential, whereas both DUNE and ESS provide little sensitivity but at least one of them is needed to fix the ω ij from disappearance data.
In Fig. 3 we show the summed χ 2 min contributions from the 0.75 and 0.85 GeV bins as a function of the value of the 3ν CP phase δ assumed to calculate the "data" to which the T-even model is fitted. In addition to the features mentioned above, we see from Fig. 3 that the test is sensitive only to δ 90 • , whereas no sensitivity appears around 270 • . This behaviour stems from the enhancement of the second oscillation maximum in the latter case (contrary to its suppression around 90 • ): only when the second oscillation maximum is smaller than the first one does the P even µe (L) fail to fit the data. Bins with E > 1 GeV are not useful in the test because of the absence of measurements at both maxima. See the appendix for further discussion. For illustration purpose we show in Fig. 3 the effect doubling the event numbers in DUNE. This shows that there is significant potential to increase the sensitivity of the test by suitable optimizations. The increased sensitivity emerges from the 0.85 GeV bin, since at this energy the DUNE baseline is close to the 2nd oscillation maximum.
The results for inverted mass ordering (IO) are qualitatively similar to the one from normal ordering (for IO we show only the relevant range of δ in Fig. 3). Further details on IO are given in the appendix. If antineutrino data are assumed (instead of neutrino data) the result is roughly obtained for δ → 2π − δ in Fig. 3, with highest sensitivity around δ 270 • . This is to be expected, since antineutrino oscillation probabilities are obtained from the neutrino ones by replacing δ → −δ (in addition to the sign-flip of the matter potential).

Summary.
We propose a largely model-independent test to search for T violation in neutrino oscillations by comparing transition probabilities at the same energy and different baselines. The test can be done under rather general assumptions covering a wide range of new physics scenarios. Within some modest assumptions, the test can be performed already with experiments at three different baselines plus near detectors. The crucial requirements are sufficient event numbers in the neutrino energy overlap region between the experiments and good neutrino energy reconstruction [48,49]. Our estimates show that with the planned long-baseline experiments DUNE, T2HK, and T2HKK, this test can be potentially carried out. In order to cover all T-violating values of δ, data for neutrinos and antineutrinos are necessary. We stress that a detector at the Tokai-Korea baseline is required in addition to DUNE and T2HK. Some optimization studies, especially in the low-energy region of the DUNE and high-energy region of the T2HKK beams, may be required. The results presented here warrant more detailed sensitivity studies based on realistic experiment simulations and statistical analyses, which we leave for future work.

ACKNOWLEDGMENTS
This project has received support from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 860881-HIDDeN, and from the Alexander von Humboldt Foundation.

SUPPLEMENATRY MATERIAL
A. Estimation of the ∆m 2 21 prior If the neutrino data come from a theory with T violation, the L-even fit converges to oscillation frequencies different from the real ones trying to compensate the wrong functional form. These different oscillation frequencies thus show either the presence of T violation or an underlying BSM mechanism that shifts the effective mass-squared differences ∆m 2 ij from their (energydependent) standard values in matter.
We characterize the possible size of this second effect using the constraints on non-standard neutrino interactions from current global oscillation data [36]. The NSI effects can be parametrized in the matter part of the neutrino Hamiltonian as where E αβ (x) = e αβ + p αβ + Y n (x) n αβ , Y n is the neutronto-electron number ratio, and q αβ are the non-standard couplings to electron, protons, neutrons for q = e, p, n, respectively.
Considering a typical value q αβ 0.1 from the global fit results in Ref. [36], we generate random values for the non-standard couplings and diagonalize the Hamiltonian to study their effect on the neutrino mass squared differences. We find that ∆m 2 21 (E) is shifted around its standard value in matter with σ = 0.09. Therefore, we conservatively choose σ 21 (E) = 10%∆m 2 21 (E) to describe the allowed deviations of the smallest oscillation frequency from the SM value via the prior term Eq. (6).
In the case of the large oscillation frequency, we find that such effect is negligible: the T-even disappearance data always fix its best-fit point to a value within the range allowed by NSIs.
Let us emphasize that this approach remains largely model independent. We consider the NSI as an effective parametrization to study how much the eigenvalues of the Hamiltonian are allowed to deviate from their standard values. This argument does not depend on the specific new physics scenario and is independent of the mixing coefficients in Eq. (1).
B. Discussion of particular behaviours of the L-even fit The consistency of the test requires that χ 2 min = 0 in the absence of T violation, and we checked that we indeed recover this result whenever data are generated for 3ν oscillations with sin δ = 0. However, depending on the specific values of the data points to be fitted, it may happen that one obtains a small χ 2 min even in the presence of T violation. In this Section, we present some cases that illustrate the typical behaviour of the L-even fit.
We show in Figure 4 the data and best-fit curves for δ = 90 • at the energies E = 0.65, 0.75, 0.85 GeV. For the clarity of the following discussion, we show the results without the smearing effect due to a finite energy resolution; the qualitative features apply also to the case including energy smearing. The middle panels of Figure 4 are identical to the left panels of Figure 2. In general, the disappearance data fix the largest oscillation frequency, whereas appearance data provide the sensitivity to TV. The two higher-energy bins are the most sensitive ones, where the fit is unable to describe the data from the T-violating oscillation probability. On the other hand, we see that for E = 0.65 GeV (left panels) the L-even model can provide a reasonable fit to the 5 appearance data points, even if the fitted function looks nothing like the true oscillation probability. It remains an interesting question, whether such a drastic deviation from the standard behaviour could be tested model-independently by other observations, for instance atmospheric neutrinos. Notice also the relatively large zero-distance effect for the fit with all baselines, which show the importance of the near detector measurement in constraining nonstandard P αβ (L → 0) = δ αβ . Increasing the error bar of this data point by one order of magnitude to 0.1 decreases χ 2 min (0.75 GeV) from 7.40 to 5.05.
The analogue plots for δ = 270 • are shown in Fig. 5, where one sees that the shape of the oscillation probability is much more oscillatory-like than for δ = 90 • and the L-even fit can quite accurately reproduce the true oscillation probability. This feature explains the lack of sensitivity of the TV test in the region around the maximally T-violating value sin δ = −1. From this Figure it becomes also clear that this behaviour does not depend on the specific baselines used here; even many more relatively accurate data points at different baselines would not be sensitive to this case.
To understand why this happens, one should look into the functional form f (L) of the fitting function and the true probability. We provide here a qualitative discussion. For simplicity we set θ 23 = π/4 and consider oscillations in vacuum; the arguments remain valid also when the matter effect is included. The 3ν appearance probability is then given by P 3ν ≈ 2s 2 13 sin 2 ∆ + 2s 13α ∆ sin ∆ cos(∆ + δ) +  for neutrinos (antineutrinos).
In performing the TV test, we are attempting to fit the probability (9) with an L-even function. According to our requirement (iv) it will have a similar ∆-dependence as Eq. (8) with δ = 0 or π: where a, b, c are three fit amplitudes with a, c ≥ 0. Disappearance data require |∆ | ≈ |∆|. The expression in Eq. (10) is equivalent to our L-even fit assuming that zero-distance effects are small, as it happens with the data we used in the analysis.
Comparing this with Eq. (9), it becomes immediately clear that it is not possible to fit the T-violating probabilities when the second oscillation maximum is suppressed with respect to the first one, i.e., for δ = +(−)90 • for neutrinos (antineutrinos). For the opposite case, i.e., δ = −(+)90 • for neutrinos (antineutrinos), the second oscillation maximum is enhanced and consistent with the condition (11). Indeed, one can show that in this case, choosing three reference values for ∆ close to (but not exactly at) the first and second oscillation maximum and the first oscillation minimum, all three coefficients, a, b, c, can be fitted.
Hence, this simple considerations confirm the behaviour visible in figs. 4 and 5 (which are based on full oscillation probabilities including matter effects): the suppression (enhancement) of the second oscillation maximum due to the sign of the term proportional to sin δ is responsible for the (lack of) sensitivity around δ = 90 • (270 • ) for neutrinos, and the opposite for antineutrinos. Higher energy bins have no sensitivity to TV due to the fact that there is no longer any data point at the second oscillation maximum.
Furthermore, it becomes also clear why this behaviour is independent of the mass ordering. This follows from the observation that the sign of the termα∆ in Eq. (9) is   [7.59] independent of the mass ordering. We give χ 2 min results for inverted ordering in Tab. II, and figs. 6 and 7 show some of the corresponding best-fit curves, which are qualitatively similar to their normal-hierarchy counterparts. Matter effects induce deviations from the above analytical arguments; however, at the energies under consideration, the matter potential is small enough that these deviations are subleading.