{'fraction_non_alphanumeric': 0.04256410256410256, 'fraction_numerical': 0.014027149321266969, 'mean_word_length': 4.486759351208209, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Inference finds consistency between a neutrino flavor evolution model and Earth-based solar neutrino measurements

We continue examining statistical data assimilation (SDA), an inference methodology, to infer solutions to neutrino flavor evolution, for the first time using real - rather than simulated - data. The model represents neutrinos streaming from the Sun's center and undergoing a Mikheyev-Smirnov-Wolfenstein (MSW) resonance in flavor space, due to the radially-varying electron number density. The model neutrino energies are chosen to correspond to experimental bins in the Sudbury Neutrino Observatory (SNO) and Borexino experiments, which measure electron-flavor survival probability at Earth. The procedure successfully finds consistency between the observed fluxes and the model, if the MSW resonance - that is, flavor evolution due to solar electrons - is included in the dynamical equations representing the model.


I. Introduction
Neutrinos are ubiquitous in astrophysical environments. In stars like the Sun, they are a byproduct of the nuclear reactions that provide the energy flux counteracting gravity. Proto-neutron stars formed following the core collapse in a supernova cool by emitting copious neutrino pairs, and similar neutrino emission takes place following neutron star mergers. Observations of those neutrinos would provide not only valuable information about the interior of the star and collapse mechanism, but also the nature and properties of neutrinos. In particular, neutrino "flavor," a property that dictates neutrino interaction with matter, significantly affects the physics of these events [1][2][3][4][5][6][7][8][9][10][11][12][13][14].
A complete description of neutrino flavor evolution and transport in astrophysical environments is technically very involved. Powerful numerical integration codes exist for obtaining solutions to the flavor evolution problem in compact object environments [15][16][17][18][19][20][21][22]. Many of these codes, however, require adopting rather rigid physical assumptions regarding the symmetries of the problem, and it has been shown in recent years that relaxing these assumptions reveals physics that had been artificially hidden (e.g., see Ref. [23,24] and references therein). Further, because of the large density of neutrinos in these regimes, neutrino-neutrino interactions introduce nonlinear elements in the transport, rendering the problem computationally taxing for the existing large-scale codes. Additionally, these nonlinear "collective oscillations" could exhibit signatures of many-body quantum correlations among neutrinos (e.g., [25][26][27][28][29][30][31][32][33]). Hence it is helpful to explore the suitability of different computational tools to treat neutrino flavor evolution and transport.
In this paper we seek, for the first time, to test the SDA procedure using real neutrino data from experiments. Data on neutrinos from CCSN are limited to one event: SN1987A [55][56][57], wherein about 20 neutrinos were detected -not a sufficiently large sample to robustly investigate flavor evolution signatures. Hence we turn to an environment from which neutrino data are copious: the Sun. Our purpose here is to explore the efficacy of inference to describe the matter-enhanced neutrino oscillations in the Sun. Furthermore, it is worthwhile to investigate possible improvements that inference might offer for the analysis of solar neutrino data. In particular, inference might offer an independent check on new methods to remove cosmogenic-induced spallation in Super-Kamiokande, a recent effort to improve the precision of solar neutrino data [58].
To that end, we challenge the SDA procedure to find a solution -in terms of neutrino flavor evolution with radius, in a steady-state model -that is consistent with three elements: the standard solar model, a model of neutrino flavor evolution, and observed neutrino fluxes at Earth. In addition, we explore the significantly more complex problem of state prediction simultaneous with the estimation of unknown model parameters.

A. Neutrino flavor evolution
We consider two-flavor mixing of solar neutrinos, which is a good approximation because of the small value of the mixing angle θ 13 [59]. Hence the angle θ 12 describes the effective mixing between the electron neutrinos, ψ e , and a linear combination of the muon and tau neutrinos, ψ x . Hereafter we set θ = θ 12 , dropping the subscripts. We assume that flavor evolution is driven entirely by coherent forward-scattering. This arises from neutrino-matter interactions, and leads to an in-medium effective neutrino mass level crossing, referred to as the "MSW resonance." [60][61][62]. The MSW resonance is associated with an enhanced e ↔ x flavor conversion probability. The neutrinos are produced at the core of the Sun as electron neutrinos, then they undergo flavor oscillations while also forward-scattering off the background electrons as they travel through the Sun. For each neutrino, this propagation is described either by the equation: or by the equation [63,64]: where ω = δm 2 /(2E) is the vacuum oscillation frequency of a neutrino with energy E, δm 2 being the mass-squared difference in vacuum. The unit vector B = sin(2θ)x − cos(2θ)ẑ represents flavor mixing in vacuum, with mixing angle θ. The function V(r) arises from neutrino interactions with the background electrons. Assuming that flavor evolution is driven entirely by coherent forward-scattering, V(r) is given by the Wolfenstein correction to the neutrino mass: V(r) = √ 2G F N e (r) [62], where N e is the electron density, taken from the Standard Solar Model [65]. Finally, the relation between the dynamical variables in Eq. (1) and Eq. (2) is: In Eq. (2), the flavor state of the neutrino is characterized by the "polarization vector" P. Eq. (3) shows that the components of the polarization vector are real numbers. In the inference method adopted in this work, the dynamical variables and input parameters need to real. For that reason, in this paper we adopt Eq. (2) to describe neutrino evolution. The density matrix ρ of the system can also be written in terms of the polarization vector as: Here, theẑ component of the neutrino polarization vector denotes the net flavor content of electron flavor minus the superposition of muon and tau flavors. Hence in the Sun, Eq. (1) is solved with the initial conditions ψ e (0) = 1, ψ x (0) = 0, or equivalently, Eq. (2) is solved with the initial conditions P z = +1, P x = P y = 0. In this configuration, the density matrix as defined in Eq. (4) represents a pure quantum state. In the forwardscattering regime, this purity (i.e. ρ 2 = ρ) is preserved through the course of flavor evolution, and one therefore has | P| 2 = 1 throughout the evolution.
We solve Eq. (2) inside the Sun, using the inference procedure to be described in later sections. In particular, the solution includes a prediction of P at the surface of the Sun. As our model does not extend beyond the Sun, we need to use the time evolution operator to relate P at the endpoint of our model (i.e. the solar surface) to the measured electron flavor survival probability, P survival , at Earth's surface 1 .
Neutrino propagation from the solar surface (denoted by ' s') to Earth (denoted by '⊕') can be analytically calculated to obtain the density matrix on Earth: where the operator U is: In Eq. (6) the entries are calculated by solving Eq. (1) in vacuum: ψ e (L) = cos ωL − i cos 2θ sin ωL, ψ x (L) = i sin 2θ sin ωL. (8) where L is the Earth-Sun distance. If one uses only the 'day' data (so that the measured neutrinos haven't passed through the earth), as we do, then P ⊕ represents the flavor state of the neutrino at the detector. Hence for the z component of P ⊕ on or near the surface of Earth, we have: P ⊕,z = 1 − 2 sin 2 2θ sin 2 ωL P s,z − 2 cos 2θ sin 2θ sin 2 ωL P s,x − 1 2 (sin 2ωL sin 2θ) P s,y .
Averaging over ωL (that is, over multiple oscillation cycles in vacuum), we obtain: P ⊕,z can be related to the measured P survival using: and hence we have the following constraint on the 'final' (i.e., solar surface) values of our state variables: 1 − sin 2 2θ P s,z − cos 2θ sin 2θP s,x = 2P survival − 1.
(12) Eq. (12) yields a linear relation between P z and P x at the surface of the Sun.

B. Neutrino data
To test the inference procedure, we use 8 B day-time neutrino flux observed by the SNO [66] and Borexino [67] experiments.
For the Borexino data, we use only the observed ppchain neutrinos, and not the CNO neutrinos. This is a reasonable choice for the Sun, as its core temperature is relatively low, meaning that few CNO neutrinos are produced. In addition, for simplicity we use day data only. The Borexino survival probabilities are listed for three discrete energies [67].
The SNO collaboration used an analytic fitting formula for the survival probability P survival [66] as a function of neutrino energy E ν . For daytime data only, the formula is: where the best-fit values of the coefficients c 0 , . . . , c 2 , along with the uncertainties, are given in Table I. For the neutrino oscillation parameters δm 2 and θ, we use the following values: θ = 0.5838 radians (15)

A. General formulation
Statistical data assimilation is an inference procedure wherein measured quantities are assumed to arise from a dynamical physical model. It is designed for cases wherein only a subset of the model state variables can be experimentally accessed. We write this model F in D ordinary differential equations: where r is the parameterization -for example, distance or time. Components x a of vector x are the model state variables. The p are any unknown parameters to be estimated, and they may vary with r. Measured quantities comprise a subset L of the D state variables. We seek to estimate the evolution of all state variables that is consistent with measurements, and to predict model evolution at parameterized locations where measurements do not exist.

B. Optimization formulation
We use a path-integral formulation of SDA, which can be summarized in three equations. The path integral is an integral representation of the master equation for the stochastic process represented by Eq. (16). We seek the probability of obtaining a path X in the model's state space given observations Y : This expression states: the path X for which the probability -given Y -is greatest is the path that minimizes the quantity A 0 , which we call our action. A formulation for A 0 will permit us to obtain the expectation value of any function G(X) on a path X: Expectation values are the quantities of interest when the problem is statistical in nature. For many estimation problems, the quantity of interest is the path itself: G(X) = X. The action is written in two terms: The first term describes Markov-chain transition probabilities governing the model dynamics. The second term is the conditional mutual information (CMI) [68], which asks: "How much information, in bits, is learned about event x(n) upon observing event y(n), conditioned on having previously observed event(s) Y (n -1)?" 2 . Simplifications are then made to write a computationallyfunctional form of A 0 , and model-specific equality constraints may be added. See Ref. [69] for a derivation of Eq. (19).
The SDA problem is then cast as an optimization, where the action is a cost function -a succinct and powerful equivalency. The cost function of the optimizer is equivalent to the action on paths in the state space that is searched. Generally, the action surface is ((D + p) × (N + 1))-dimensional, where N + 1 is the number of discretized model locations, taken to be independent dimensions. One seeks the path X 0 = {x(0), . . . , x(N), p(0), . . . , p(N)} in state space on which A 0 attains a minimum value. Minima are found by requiring that small variations to the action vanish under small perturbations [70], thereby enforcing the Euler-Lagrange equations of motion upon any path. We extremize the cost function via the variational method.
After many simplifications (see Appendix A of Ref [50]), the Markov-chain term (first term of Eq. (19)) reduces to a "model error", which describes the divergence of the prediction from model dynamics. The mutual information (second term of Eq. (19)) reduces to a "measurement error", describing the divergence of the prediction from measurements. (For a pedagogical treatment, see Ref. [69].) The action A 0 used in this paper is written as: The model error, A model , imposes adherence to the model evolution of all D state variables x a . The outer sum on n runs through all odd-numbered discretized locations. The sum on a runs through all D state variables. The terms within the first and second sets of curly brackets represent the errors in the first and second derivatives, respectively, of the state variables.
The measurement error, A meas , imposes adherence to measurements. The variables y l , for l = 1, . . . , d, represent the d quantities measured at locations r m ∈ {meas}, where N meas is the total number of locations. These are to be compared against the quantities h l,m (x), where h l,m are transfer functions that relate the state variables in the model to the quantities being measured, at each location. In our optimization design, the measured quantities are the values of P z for each neutrino energy, at two locations: (i) the center of the Sun, and (ii) the surface of Earth. At the Sun's center, the "measurement" of P z can be compared directly against the model P z at the same location, rendering the transfer functions trivial: h 0 ( P) = P z for each neutrino. On the other hand, since our model grid does not extend beyond the solar surface to the Earth, the P z measurement at Earth is compared against an extrapolated P z value derived from the Polarization vector at the Sun's surface, as shown in Eq. (10). In other words, measuring P z at Earth is equivalent to measuring a linear combination of P z and P x at the surface of the Sun. Therefore, the transfer function at the Sun's surface becomes: for each neutrino energy 3 . The measurement term can then be written as where k = {1, . . . , N ν } are the neutrino energy bins, P k is the polarization vector for the k th energy bin in the model equations of motion (Eq. (2)), with components P x,k , P y,k , P z,k , and P meas z,k is the measurement of P z at the specified location, directly associated with the survival probability at that location (Eq. (11)). R and L are the solar radius and the earth-sun distance, respectively.

C. Identifying a lowest minimum of the action
The action surface of a nonlinear model will be nonconvex. To identify a lowest minimum, we perform an iterative annealing in terms of the ratio of model and measurement error, R f and R m , respectively 4 [71]. It works as follows.
We define R m to be a constant (in this paper it is 1.0), and R f as: R f = R f ,0 α β , where R f ,0 = 10 −3 , α = 2.0, and β -the annealing parameter -is initialized at zero. Relatively free from model constraints, the action surface is convex. Then we increase β recursively, each time 3 Note that in our previous papers [50][51][52][53][54], the transfer function was always trivial, since we simply employed simulated data at the far boundary of the model domain. 4 More generally, R m and R f are inverse covariance matrices for the measurement and model errors, respectively. In this paper the measurements are taken to be mutually independent, rendering these matrices diagonal.
recalculating the action, toward the deterministic limit of R f R m . The aim is to remain sufficiently near to the lowest minimum so as not to become trapped in a local minimum as the model dynamics resolve the surface 5 .

IV. The task for optimization
Our model beams were assigned energies corresponding to the experimentally detected neutrino energies. We placed two constraints on each beam. One was an assumed "measurement" of P z at the center of the Sun: 1.0 for each beam, or: pure ν e flavor. The other measurement was P survival , the electron flavor survival probability of neutrinos measured at Earth's surface, related to the P z at earth through Eq. (11).
The task for the procedure, depicted in Fig. 1, was to take those two constraints, together with the potential V(r) from the standard solar model [65] and the dynamics of flavor evolution (Eq. (2)), to find a solution -in terms of the radial evolution of the polarization vectors -that is consistent with model and data. We expected this inference task to be a challenge, specifically for estimating the polarization vector at the solar surface. From  FIG. 1. Schematic of the inference task. The constraints provided to the procedure are an assumed initial condition at the solar center (blue circle) and a measured survival probability at Earth (yellow circle). Meanwhile, the model consists of the potential V(r) and flavor transformation dynamics (Eq. (2)) through the Sun. The inference task is to consider both the constraints and the model, to predict the flavor evolution within the Sun (red region), while also accounting for the transformation (grey wave) of the Earth-based measurement (yellow circle) to the solar surface (red circle).
Eq. (12), one can see that a single measured value of survival probability at Earth corresponds to a (P x ,P z ) pair at the solar surface. In principle, many possible pairs could satisfy that relation. The SDA procedure is tasked with finding a pair that also obeys model dynamics (Eq. (2)).

A. Details
The procedure is given full knowledge of model parameter values (except for the parameter estimation of Section V B), all listed in Table I. The three energies E B , i correspond to three discrete bins in the Borexino experiment [67]; SNO employed an energy range, which we discretized into eight values [66]. To assign a unique value of the matter potential V(r) at each radial location, we used a piecewise linear interpolation: the V(r) term adopted from the solar model in Ref. [65] contains values at 1,219 discrete radial locations, while our model contains 121,901 6 .
The optimization was performed by the open-source Interior-point Optimizer (Ipopt) [72]. Ipopt employs a Hermite-Simpson method of discretization and a constant step size. We employed 121,901 steps and a step size of δr of 2.85237. The discretization of state space, calculations of the model Jacobean and Hessian matrices, and the annealing procedure are performed via a Python 6 Our choice to increase the number of discrete locations by a factor of 100 was motivated by our aim to resolve the oscillations; namely, the step size had to be sufficiently small.
interface [73] that generates C code to be read by Ipopt.
Simulations were run on a computing cluster equipped with 201 GB of RAM and 24 GenuineIntel CPUs (64 bits), each with 12 cores.
To compare solutions to the model dynamics, we generated a simulation via forward integration, with all beams initialized at [P x ,P y ,P z ]=[0,0,1] at the solar core, and using the same discretized grid as the optimization. This integration was performed by Python's odeINT package, which uses an adaptive step. Our complete procedure, including forward-integration codes, codes to interface with Ipopt, instructions for designing and running experiments on our supercomputing cluster, and examples for new users, can be found in a publicly available repository [74].
For all experiments to be described in Section V, four independent paths were initialized randomly. That is, each initialization consisted of as many random choices as there are dimensions in the model 7 . The user-defined search range for state variables (P x ,P y , and P z ) was: [-1.0:1.0], the full dynamical range for each. For the search range used in parameter estimation, see Section V B.

V. Results
Key findings are as follows: • When the MSW transition within the Sun is included in the model dynamics, the procedure finds a solution consistent with both model dynamics and measured survival probabilities P survival at Earth, within the published experimental errors on P survival , for both SNO and Borexino data. (Alternatively, when the MSW transition is ignoredor, V(r) is set to zero -a solution compatible with both model and measurements is not found.) • A preliminary parameter estimation shows that the measured survival probabilities contain information about model mixing angle θ (Eq. (2)).

A. State prediction without parameter estimation
The left panel of Fig. 2 shows the logarithm of the action as a function of annealing parameter β, over the 7 As noted, generally the dimensionality of the action surface is ((D + p) × (N + 1)), where D, N + 1, and p are the number of state variables, discretized model locations, and parameters, respectively. All optimization procedures described in this paper, however, contained at most one unknown parameter, which was assumed to have the same value at all (N+1) discretized locations. Thus, the dimensionality of these procedures was: D × (N + 1) + p.  4 10.0 E SNO, 5 10.5 E SNO, 6 11.0 E SNO, 7 11.5 E SNO, 8 12 course of annealing, for the case in which the matter potential V(r) is taken from the standard solar model (Section III), a scenario wherein the electron number density effects a significant MSW transition. Beginning at β = 0, the action increases as the weight of model error is increased. Gradually, however, it levels off. This "plateau," around β ∼ 16 to 20, indicates that a solution consistent with both model and measurements has been found 8 (for a detailed study of these action(β) plots, see Ref. [51]). All four randomly-initialized paths converged to an identical solution.
Indeed, the predictions of state variable evolution well match the true model evolution. The left panel of Fig. 3 shows the true (blue) versus predicted (red) 8 The further increase in the action beyond β ∼ 21 is most likely due to discretization error, as the optimization technique uses a different method of discretization compared to the forward integration.
state variable evolution of P x (top), P y (middle), and P z (bottom), for the first of the three beams corresponding to the Borexino data; the result is representative of all beams in both SNO and Borexino models. Fig. 4 shows detail that cannot be discerned by eye on Fig. 3. At top are three segments of the evolution of P z : the first thousand steps, beginning at the solar center (left), middle thousand (middle), and final thousand ending at the solar surface (right); each section is roughly one hundredth of the full 121,901-step series. The blue solid dot at far left denotes the location of the assumed initial conditions on P z , and the red dot at far right indicates the first location of prediction, given the survival probability measured at Earth (Earth is not depicted in the figure).
At bottom are the Fourier transforms of P z corresponding to each segment. The predicted oscillation frequency throughout the Sun is perfect, to within the resolution permitted by the density of sampled locations. These results are insensitive to the addition of the published maximum experimental errors on survival probability, for both SNO and Borexino.
As a test that the procedure recognizes consistency between measurement and model, we set the matter potential V(r) to zero and repeated the optimization. Setting V(r) to zero essentially tells the procedure that the electron number density inside the Sun is not sufficiently high to effect an appreciable MSW transition.
The right panel of Fig. 2 shows the resulting plot of action-versus-β. In contrast to our original result, here the action never attains the signature "plateau" indicative of a successful optimization. Rather, it increases exponentially. The corresponding predictions appear in the right panel of Fig. 3: agreement with true (blue) model is poor. As with the original experiments, all four randomly-initialized paths converged to this solution. Together, Figures 2 and 3 and convince us that the MSW effect within the Sun cannot be neglected if one seeks to account for the measured Earth-based survival probabilities.
We conducted a further test that the SDA procedure recognizes consistency between measurement and model. As noted, one vital component of a successful prediction was a (P x , P z ) pair at the solar surface that was consistent with both the survival probability measured at earth and the model dynamics (Eq. (2)), together with the Earth-Sun transformation function of Eq. (21).
With this in mind, we replaced the true published measured survival probabilities with "test" survival probabilities that are not physically possible; specifically: values above 1.0 and below 0.0. For these cases, the action(β) plot behaved similarly to the right panel of Fig. 2: an exponential increase (not shown). That is, the model was unable to find a solution that satisfied both model and measurements -as is the expected outcome.

B. State prediction with simultaneous parameter estimation
Adding parameter estimation to the inference task renders it significantly more challenging [75,76]. The root of the difficulty is that parameters -unlike state variables -do not obey a known dynamical law, so there is no straightforward way to correlate state variable evolution with parameter estimate error 9 . With that in mind, we took the preliminary step of setting the mixing angle θ of Eq. (2) as an unknown parameter to be estimated. In these tests, the true value of θ was 0.58, and the permitted search range was 0.001 to 1.571. Fig. 3. Top: from left, three 1,000-point segments of the P z evolution, beginning at the solar center (far left) -where the blue dot denotes the assumed initial condition, and ending at the solar surface (far right)where the red dot denotes the location where P x and P z must be consistent with the measurement at Earth (Earth is not depicted). Bottom: Fast Fourier decomposition for the respective regions at top, showing that the oscillation frequency is well predicted throughout the Sun.

FIG. 4. Details of the predictions shown in
Our first parameter attempt failed (not shown), and to ascertain whether computational expense might be the cause, we temporarily removed the Earth-to-Sun transformation of Eq. (21) from the model dynamics. To be clear: our aim here was to offset the increased computational expense incurred with parameter estimation by decreasing the complexity of the equations of motion, to identify whether computational expense was the culprit in the failed first attempt. Across four randomly-initialized paths, the resulting estimates of θ were, respectively: [1.57, 0.98, 0.01, 0.59] radians. That is, one of the four paths identified a solution near the true value of 0.58. Fig. 5 shows the corresponding state prediction for that estimate. Encouraged by this mild success, we restored the Earth-Sun transformation (Eq. (21)) -and in addition, we amended the requirement regarding a match to the Earth-based measurement of survival probability. Specifically, instead of using only the final values of P x and P z at the solar surface, we averaged P x and P z over the last 1000 radial locations in the model and compared those averages to the measured survival probability, using the Sun-Earth transformation described earlier. The motivation for this change was that, without a fixed mixing angle, the rapid oscillations in the predicted path could, in principle, be out of phase with the model path. With this change, the results improved: the estimates of θ across four independent paths were: [0.79, 0.78, 0.56, 0.89] radians, respectively. Fig. 6 shows the cor-responding state predictions for the closest estimate of 0.56.
There is much honing to be done to yield reliable results from parameter estimation. We need to ascertain why the problem invites multiple solutions across independent paths searched, and aim to increase the fraction of successful paths. We note, however, that based on this preliminary study, inference can indeed extract information about the mixing angle in Eq. (2) from Earth-based measurements.

VI. Conclusion
For the first time we have applied inference to astrophysical neutrino data. The ability of the SDA procedure to identify state predictions consistent with both model and Earth-based measurements is encouraging, particularly given the multiple unphysical solutions at the solar surface that those measurements permit. Further, while the parameter estimation needs honing, the results presented in Section V B indicate that inference-based procedures do possess the capability to recognize the dependence of Earth-based flavor survival probability measurements on the neutrino mixing angle θ. Finally, for the interested reader, the Appendix addresses a topic noted briefly earlier: the phenomenon of decoherence through neutrino propagation.   6. State predictions akin to Fig. 3, again with mixing angle θ recast as an unknown parameter to be estimated, and including the Earth-Sun transformation. Importantly, in this version of the experiment, the measurements were taken to correspond not only to the (P x ,P z ) pair at the solar surface R , but rather to the average of those two values over the final 1,000 discretized model locations in the Sun. See text for explanation. One out of four paths converged to this solution, which corresponds to a value of θ = 0.56, taken at a value of β of 20.

Appendix: Coherent versus incoherent detection
In Section II A, we wrote a transformation matrix to carry the polarization vector between the Earth and Sun. In doing so, we omitted an accounting of possible decoherence. It might be valuable in the future to investigate whether including decoherence in the model improves state and parameter estimation. Here we explain the phenomenon to the interested reader.
When a neutrino arrives at earth from a sufficiently distant source, it is in general a coherent superposition of more than one mass (i.e., propagation) eigenstate. However, since the different propagation eigenstates have different masses and therefore propagate at different velocities, their wavepackets eventually become spatially separated-to the extent that any neutrino interaction in a detector would only involve the participation of one of these mass eigenstates at a time.
Since the individual mass eigenstates themselves do not oscillate in time (only their coherent superpositions do), the probability of detecting a particular mass eigenstate in a certain flavor (e.g., P 1e = ν e |ν 1 ) is independent of time or position, once the neutrino is propagating in vacuum (i.e., if there are no matter effects). For the purposes of detection, since only one mass eigenstate participates at a time (assuming sufficient wavepacket separation), one may treat the neutrino flux arriving at the earth as an incoherent mixture of mass eigenstates rather than as a coherent superposition. Mathematically, this can be described as follows.
Consider a single neutrino existing as a coherent superposition of mass eigenstates. Its density matrix in the mass basis may be written as follows: where n ν i = a † ν i a ν i are expectation values of the number operators for the mass eigenstates i ∈ {1, 2}, and ρ 12 depends on the relative phase between the mass eigenstates. For a neutrino propagating in vacuum, the n ν i remain invariant in time, and only the off-diagonal entries are time-dependent.
In contrast, for a neutrino that can be considered to have essentially devolved into an incoherent mixture of mass eigenstates, the density matrix is simply In each case, one may ask the question of how much can be learned from a measurement (i.e., a detection). First, let us define the quantities n ν e and n ν x , which are the expectation values of the flavor-basis number operators. These can be calculated as follows: n ν e = a † ν e a ν e = Tr ρ ν a † ν e a ν e = n ν 1 cos 2 θ + n ν 2 sin 2 θ + 2 Re{ρ 12 } sin θ cos θ, n ν x = a † ν x a ν x = Tr ρ ν a † ν x a ν x = n ν 1 sin 2 θ + n ν 2 cos 2 θ − 2 Re{ρ 12 } sin θ cos θ.
Here, we have used the definitions from earlier in the section, and the familiar unitary transformation between the flavor and mass eigenbasis: a ν e a ν x = cos θ sin θ − sin θ cos θ a ν 1 a ν 2 .
Suppose, for instance, one detects enough neutrinos at a certain energy so as to statistically obtain a sufficiently accurate determination of n ν e and n ν x at a given location. Then, in case of neutrinos that are still coherent superpositions by the time of arrival at the detector, one cannot uniquely determine n ν 1 and n ν 2 from n ν e and n ν x , unless one also measures ρ 12 , which is physically impossible. Instead, a possible workaround may be to measure n ν e and n ν x at multiple locations, as was explored in Ref. [53] in the context of supernova neutrino detection.
On the other hand, if the neutrinos that arrive at the detector are an incoherent mixture, then ρ 12 = 0, and a measurement of n ν e and n ν x at a single location is sufficient to uniquely determine n ν 1 and n ν 2 .