Microscopic and Macroscopic Effects in the Decoherence of Neutrino Oscillations

We present a generic structure (the layer structure) for decoherence effects in neutrino oscillations, which includes decoherence from quantum mechanical and classical uncertainties. The calculation is done by combining the concept of open quantum system and quantum field theory, forming a structure composed of phase spaces from microscopic to macroscopic level. Having information loss at different levels, quantum mechanical uncertainties parameterize decoherence by an intrinsic mass eigenstate separation effect, while decoherence for classical uncertainties is typically dominated by a statistical averaging effect. With the help of the layer structure, we classify the former as state decoherence (SD) and the latter as phase decoherence (PD), then further conclude that both SD and PD result from phase wash-out effects of different phase structures on different layers. Such effects admit for simple numerical calculations of decoherence for a given width and shape of uncertainties. While our structure is generic, so are the uncertainties, nonetheless, a few notable ones are: the wavepacket size of the external particles, the effective interaction volume at production and detection, the energy reconstruction model and the neutrino production profile. Furthermore, we estimate the experimental sensitivities for SD and PD parameterized by the uncertainty parameters, for reactor neutrinos and decay-at-rest neutrinos, using a traditional rate measuring method and a novel phase measuring method.

By virtue of the more and more precisely measured phenomenon of neutrino oscillation [1], the quantum coherence of neutrino mass eigenstates can routinely be observed on a macroscopic level. There are two common bases that are utilized to express the quantum state of neutrinos, viz. the mass eigenstates and the flavor eigenstates. While the former determines how neutrinos propagate, the charged-current interaction, for the production and detection of neutrinos, is characterized by the latter. Evolution of neutrinos is effectively encapsulated by the flavor transition probability (FTP), P να→ν β , which represents the probability that a neutrino produced as flavor ν α is detected as flavor ν β . The non-trivial mixing between these two bases, described by the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix [2][3][4], implies that P να→ν β is not diagonal. Consequently, after a neutrino is produced as flavor eigenstate, it propagates in a superposition of mass eigenstates. Such superposition in the Hilbert space describes quantum coherence, which could lead to observational interference patterns. In fact, the loss of coherence, decoherence, represents a transition from the quantum to classical level, for it describes the loss of interference pattern from the correlation between quantum states [5,6]. In addition, the smallness of the neutrino masses and of their difference admits that quantum coherence can be observed macroscopically through the FTP. Moreover, with increasing precision of neutrino oscillation experiments, the degree of the quantum correlation may be measured more accurately, therefore, decoherence effects warrant further investigation.
As all observable effects of mixed quantum states, neutrino coherence is expected to be lost at some stage. Theories for calculating quantum decoherence in neutrino oscillation (or neutrino decoherence) have been widely discussed in the literature. These theories include the degree of wavepacket (WP) separation calculated through quantum mechanics (QM) [7][8][9][10][11] and quantum field theory (QFT) [9,[12][13][14][15][16]. Another approach deals with the effect of losing information to an open quantum system from Liouville dynamics calculated with density matrices (e.g. through the Lindblad equation) [17][18][19][20][21][22][23][24][25], or through the Wigner quasiprobability distribution [21,26,27]. There is also literature comparing one approach with another, for instance, QM vs. QFT approach for WP separation in [9], Lindblad equation vs. the WP format in [19], and Lindblad equation vs. Wigner quasi-probability distribution in [21]. Among these theories, QFT is able to describe the situation on the most fundamental level by considering neutrino oscillation as the propagator of a full process described by a Feynman diagram. However, the open quantum system method is more tailored for quantum decoherence effect in a generic way by considering a system of interest in an environment. In this sense, decoherence in a system reflects loosing information to the environment, and is calculated by tracing out states of the environment entangled to the system. Nonetheless, regardless of the way one chooses to formulate the decoherence effect, it will result in additional terms to the coherent interference patterns caused by quantum correlation. For neutrino oscillation, the decoherence effect appears as a complex function Ψ jk in the FTP as where ψ jk is the coherent phase, usually estimated as ∆m 2 jk L 0 /(2E 0 ) for some traveling distance L 0 and energy E 0 . The decoherence term Ψ jk would, in general, erase the interference pattern, hence, |Ψ jk | ≤ 1. However, since it could also be complex, is might also cause a phase shift w.r.t. ψ jk .
In this work, we introduce the concept of the open quantum system method to the QFT calculations by considering the propagator describing neutrino oscillation as the system of interest, and everything else in the diagram as the environment, which is to be integrated out. Furthermore, since neutrino oscillation is considered as a phenomenon resulting from the coherence of kinematics between mass eigenstates, the states of the environment we integrate out are in the phase space (PS). In particular, if a state is described by creation and annihilation operators represented in the coordinate and/or the momentum space, we call it the "Fock-PS"; on the other hand, if a state is described by occupation numbers on a PS forming a Wigner quasi-probability distribution [26], we call it the "Wigner-PS". Notably, since the Fock space representations for mass basis and flavor basis are unitarily inequivalent with each other [28,29], at least one of them must be unphysical, and the debate on which of them is unphysical is still on-going, e.g. [30][31][32][33][34]. As both representations approximately agree with each other in the relativistic limit, we choose to build the Fock-PS for mass states here, with flavor states represented by a superposition of mass states. Nonetheless, one can also build a flavor based Fock-PS for the layer structure, for it does not specify the representation we choose.
Concretely, we start from following the QFT description for calculating neutrino oscillation in [12], which already includes integrating out the momentum space of the external particles, so we complete the picture by also integrating out the space-time components of the PS (although the PS does not include a temporal dimension, the kinematic of the states is time-dependent, therefore when we say "PS variables", a temporal component is also included). Furthermore, we focus on the structure of the PS while calculating the FTP, which will be called the "layer structure", composed of three layers. Vertical-wise, the layer structure includes three layers, namely the microscopic layer, the physical layer and the measurement layer; and horizontalwise, each layer represents a PS composed of space-time variables and momentum variables, which are able to determine the kinematics of the neutrinos fully. The microscopic layer relates to the theories in the literature [7-9, 12-15, 17-21, 26, 27], which can be represented by either the Fock-PS or the Wigner-PS on which Fock states and the Wigner quasi-probability distribution are represented, respectively. With the layer structure we are able to account for the decoherence effect from information loss to the environment, from a microscopic to macroscopic level. Different from the existing literature in which one also calculates neutrino decoherence on a macroscopic level [15,35], we focus on classifying and understanding decoherence in a generic picture. At the end, neutrino decoherence for continuously emitted neutrinos is mainly parameterized by four uncertainties appearing on the layer structure. Those are the coordinate/momentum uncertainty on the microscopic layer (quantum effects quantified by σ x /σ p ), and that on physical layer (macroscopic effects such as energy resolution or neutrino production profile σ L /σ E ), providing an interface between microscopic mechanisms and the macroscopic experiments. As for non-continuously emitted neutrinos, we would simply have an additional temporal uncertainty on the physical layer, σ T .
As for the phenomenology part of neutrino decoherence, roughly speaking, measurements of neutrinos produced in both long and short baseline experiments and in the atmosphere are best fitted with neutrinos considered as fully coherent, see for instance [36] for an updated global fit; as for neutrinos produced outside the Earth, such as solar and supernova neutrinos, it is best described as fully incoherent. There are also many discussions on neutrino decoherence phenomenons, such as general neutrino decoherence for reactor experiments in [37][38][39]; gravitation fluctuation or cosmological effects causing atmospheric neutrino decoherence in [18,[40][41][42]; matter effect responsible for accelerator or atmospheric neutrino decoherence in [22,[43][44][45][46]. Our structure carries the potential of including all the mechanisms above and more, since both the QFT approach and the open quantum system concept is included. Hence, we do not go into the details of these theories but give a generic picture on what measurable parameters it could reflect on. Furthermore, we also include decoherence signatures caused by classical uncertainties or the ignorance of the observer, such as the production profile (for instance, the exact shape of the neutrino source) of the neutrino and the energy reconstruction model.
Notably, we introduce the phase wash-out (PWO) effect, which is an averaging effect with respect to some phase structure, washing-out the oscillation signatures. In this paper, we will show that decoherence signatures of neutrino oscillation can all be described by some PWO effects, and the distinction between difference decoherence parameters comes from different dependence on the phase structure(s). In other words, while quantum coherence reflects on the oscillation signature, decoherence can be described by some wash-out of such signature, and the PS (e.g. the traveling distance and the energy of neutrinos) dependence of each decoherence effect comes from the formalism of the phase structure being washed-out. Therefore, the PWO effect arising from neutrino decoherence results in a damping and/or phase shift signature to the oscillation. We further analyze the damping/phase shift signatures w.r.t. both classical and QM uncertainty parameters for reactor/decay-at-rest (DAR) neutrinos. In particular, for damping signatures which are expected in all literature mentioned above, we estimate the sensitivity of the parameters by directly analysing the neutrino count rate in a conventional way. However, the phase shift signals are estimated to be not as suitable for such method, so instead, we evaluate the possibility of measuring the distance dependence of the oscillation phase for the phase shift signals, considering that a moving detector is possible.
The paper is organised as follows. In Sec. 2, we introduce the layer structure and calculate neutrino FTP for neutrinos propagating in vacuum throughout the layers, as a demonstration for decoherence com-ing from only the production and detection site. In Sec. 3, with the help of the layer structure, we classify neutrino decoherence into "state decoherence" and "phase decoherence", where the former is related to the WP separation mechanism and the latter to the information loss mechanism. We also show how both types of decoherence are the result of the PWO effect. In Sec. 4, we discuss two kinds of phenomenological analyses for the decoherence parameters, namely, the "rate measuring method" (RMM) and "phase measuring method" (PMM). In particular, we estimate the sensitivity to the parameters for reactor and DAR neutrinos. We summarize our results and conclude with some remarks in Sec. 5. Technical details are delegated to appendices. Before we discuss the physics in detail, we provide a glossary of terms and definitions needed in this paper.

Phase space (PS) variables
Three coordinate variables, three momentum variables (composing the six-dimension PS) and a temporal variable.
Layer structure Composed of three layers (layer 1-3) from the microscopic Hilbert space to the macroscopic measurement space, where all spaces are represented by PS variables. The structure is illustrated in Sec. 2, and its value of providing a simple and generic picture of decoherence effects is shown in Sec. 3.

Weighting function
Localized distributions that characterise uncertainties included in the LMO. Examples of uncertainties from layer 1 and layer 2 for neutrino oscillation are summerized in Sec. 3

.2 and Sec. 3.3, respectively.
Phase wash-out (PWO) effect An averaging effect which washes out oscillation signatures by introducing a damping term and a phase shift term. Mathematical formalism and properties are given in Appendix A.
Uncertainty parameters (σ n ) Widths of the weighting functions w.r.t. some PS variable n which parameterize decoherence signatures. Some analysing methods, as well as its sensitivity estimation of three relevant uncertainty parameters are shown in Sec. 4 for neutrino oscillation experiments.

State decoherence (SD)
Decoherence by the separation of superposition states on the physical layer, which is equivalent to a PWO effect on the Wigner-PS under a factorisation condition (see Appendix C) and is dominated by uncertainties on layer 1 (see Sec. 3.2).

Phase decoherence (PD)
Decoherence by the PWO effect on the physical layer dominated by the macroscopic uncertainties on layer 2 (see Sec. 3.3).

The Layer Structure
In Fig. 1, we show the layer structure for calculating the expectation value of some observable, which is composed of three layers of PS, from microscopic to macroscopic, including the "microscopic layer" (layer 1), the "physical layer" (layer 2) and the "measurement layer" (layer 3). As illustrated in the introduction, this structure combines QFT and the concept of having an open quantum system while including statistical effects for actual measurements. Quantum effects, such as the superposition of states in the Hilbert space are described on layer 1, and it is impossible for both coordinate and momentum uncertainty to be zero due to the uncertainty principle. These uncertainties are parameterized as σ p and σ x which are in general independent of each other; the former is for uncertainties from external states on the mass-shell, which are described as WPs in momentum space; the latter is for uncertainties around the vertices in coordinate space, i.e. how non point-like the effective vertices are in a simplified effective diagram with only the external states and the neutrino propagator such as Fig. 2. On the other hand, layer 2 describes our ignorance towards the system in the classical regime, where the probability is summed (integrated) over, instead of the amplitude. The uncertainties on this layer include those transmitted from the first layer, and additional ones, which are usually, but not necessarily, macroscopic. These additional (macroscopic) uncertainties contain the energy uncertainty σ E , such as the energy resolution and the energy reconstruction model, on top of the coordinate uncertainty σ L , e.g. the neutrino production profile. More examples contributing to σ p , σ x , σ E and σ L will be discussed in the following sections when its corresponding layer is introduced for calculating the FTP of neutrinos. Consider the double slit experiment observed by taking a photo of the interference pattern as an analogy: the two slits represent uncertainties on the first layer, while the resolution of the camera taking the photo is on the second layer. The former creates freedom for a superposition state while the latter is a classical effect. For the case of neutrino oscillation, although the uncertainties on layer 2 are macroscopic, we still observe quantum coherence, due to the smallness of neutrino mass splitting. In the following, we will discuss the layer structure more formally, including aspects in Fig. 1, such as the representation of the PS and the layer-moving-operators (LMOs), LMO, connecting each layer. In particular, an important remark is that, as we will show later, the uncertainty parameters on each layer are the width of a weighting function in the LMO, which carries information of the environment entangled to the system (e.g. the propagating neutrino). The uncertainty parameters therefore characterize neutrino decoherence in experiments. Figure 1: Illustration of the layer structure, and the notation of each phase space variable deciding the kinematics of states: t,t, T, T 0 are the temporal variables; x,x, L, L 0 are the spatial variables; and p,p, P, P 0 are the momentum variables; E 0 is the energy and Ω represents the solid angle. The layers are linked by the layer-moving operator in Eq. (2), and the uncertainties are discussed in the text.
Each layer in our structure could have different representations for the PS, for instance, the first layer could be represented by either the Fock-PS or the Wigner-PS; for the representation of the second layer, we take the expectation values of the PS variables on the first layer assuming massless neutrinos (more on this in Sec. 2.3), which will be called the relativistic-PS; as for the third layer, we simply represent the PS variables in terms of the expectation values of the relativistic-PS, which should coincide with the actual measurement values. Otherwise, it will have a dependence on the mass of the neutrino, and the states will collapse to a certain neutrino mass state, giving us no oscillation. Hence, we call such representation the measurement-PS. Each PS is composed of its own temporal, coordinate and momentum space, while the energy is implied by the momentum variables through the dispersion relation, and the notations are given in Fig. 1. The layer structure could generically be applied for the calculation of the expectation value for any measurement, but we only focus on the calculation of the FTP for neutrino oscillation in this work.
As illustrated in Fig. 1, the layers are connected by the LMOs, LMO i , moving a quantity B i (x i , p i ), such as the FTP or the transition amplitude, from layer i to layer i + 1. This is done by integrating out the PS variables, x i , p i , of layer i, while considering additional uncertainties by including the weighting function W i , such that In particular, corresponding to the open quantum system concept, B 1 would be the system of interest, and W 1 includes the environment entangled with it. Furthermore, corresponding to our notation of phase space variables in Fig. 1, x 1 = (t, x), p 1 = p for the Fock-PS; x 1 = (t,x), p 1 =p for the Wigner-PS, representing the PS for the occupation number of quasi-probability distributions [26]; x 2 = (T, L), p 2 = P for the relativistic-PS, and x 3 = (T 0 , L 0 ), p 3 = P 0 for the measurement-PS. Each of these PS will be explored one by one in the following subsections for the neutrino case. Also, B i on the first, second and third layer represents the system of interest, the observable and the measured value, respectively. Moreover, W i is a probability density function (PDF) defined as In fact, the normalization of W i does not matter for now, as we will show later in Sec. 3.1 that the FTP will automatically be normalized. Nonetheless, we define the weighting function as a PDF simply for the convenience to observe the width, relating to the definition of "width" described in Appendix A. Moreover, x i+1 and p i+1 are the next layer variables, usually defined as (functions of) the expectation values of x i and p i . Hence, the width of W 1 (x, p; T, L, P) gives microscopic quantum uncertainties of t, x and p, as σ t , σ x , and σ p , respectively. Equivalently, that of W 2 (T, L, P; T 0 , L 0 , P 0 ) gives macroscopic statistical uncertainties for T , L and P as σ T , σ L and σ P , respectively. In fact, if the weighting function can be written as , the layer moving operator is an act of convolution between W i and B i , see Appendix A. In addition, LMO i refers to calculating the expectation value of the quantity B i , and the layer structure is mathematically fibre bundles [47]. In other words, looking from upper layers to lower layers, each PS point (x i+1 , p i+1 ) can be expanded into a whole PS composed of (x i , p i ) on the lower layer. Also, on each layer, operations mapping one state to the other could be made, depending on what we wish to observe. Viewing the LMOs as vertical operators, such operators can be referred to as horizontal operators such that the state remains on the same layer.
indicating that the uncertainty principle between x i , p i remains fulfilled on each layer. This is because LMO i along with e ix i p i means to first project everything onto the x i or the p i space, and then integrate over that space, while the projection process secures the uncertainty principle. In fact, this is exactly the case for the position-space representation of the wavefunction for some considered particle. In this case, B 1 (x, p) = e ixp∆ (p), where∆(p) is the propagator in momentum space. We will demonstrate this explicitly for the case of neutrinos in Sec. 2.1. As a matter of fact, if B i = exp(iη(x i , p i )), for some phase structure η(x i , p i ), the LMO meets the condition of giving rise to a phase wash-out (PWO) effect described in Appendix A. The PWO effect is an averaging effect over the phase structure caused by the non-trivial width of the (normalized) weighting function, resulting in an additional suppression term Φ as where |Φ(x i+1 , p i+1 )| ≤ 1, ∀ (x i+1 , p i+1 ). Only when the weighting function is symmetric with respect to the phase structure would Φ be a real function (see again Appendix A). Furthermore, when there is a substructure of B, i.e. B i = ν B νi , then the summation rule is simply Finally, the determination of the measurement expectation value of the FTP (P 3 ), is by doing the statistical averaging (LMO 2 ) over the FTP on the physical layer (P 2 ); P 2 , however, can either be calculated by squaring the transition amplitude (A 2 ) on layer 2, or directly by moving up (LMO1) the quasi-probability distribution (P1) from the Wigner-PS. Here, A 2 is calculated by integrating over all the quantum-mechanical configurations of the environment (LMO 1 ), moving the system of interest in the Fock-PS (A 1 ) up to the physical layer; and LMO1 performs an effective statistical averaging over an effective FTP, the quasiprobability distribution, for the quantum-mechanical superposition effect in the Wigner-PS [26,48]. That is, In the following subsections, we will introduce each layer and its part in calculating the expectation value for the measurement of the FTP for neutrino oscillation in vacuum.

Microscopic Layer (Layer 1): QFT Transition Amplitude
In this subsection, we calculate the transition amplitude on the Fock-PS with the QFT approach. Although QM can also describe neutrino coherence and decoherence on the Fock-PS, it could not answer a number of Figure 2: A simplified Feynman diagram where the neutrino propagating a macroscopic distance is treated as a propagator of a full diagram, and the kinematics of the external particles are described with wavepackets. At the production/detection vertex site, the diameter of the shaded blue/green areas represent the uncertainties of external wavepackets projected onto the coordinate space, and their mean value is labeled as x 1 /x 2 . The inner circles with solid lines at both sites are the additional coordinate uncertainties from the blob vertices for the internal states regardless of the external particles. In other words, it represents the uncertainties of x 1 /x 2 by g P (x 1 )/g D (x 2 ) in Eq. (8). Hence the total uncertainty on the coordinate space at this layer would be the diameter of the dashed-lined circles.

Production Detection
questions while the QFT approach can, see for example [9,12]. Moreover, for the purpose of investigating the quantum decoherence effect and its implication for fundamental physics, it is necessary to use the QFT framework, even for the scenario of vacuum propagation, since the weighting functions on the first layer could originate from uncertainties of the interactions around the vertices. Moreover, in order to have the weighting functions which are determined by the states entangled to the neutrinos explicitly, the phase space of this layer will be p/x, the momentum/space-time coordinate of the neutrino given by its entangled states. As for the weighting functions, we treat the external particles as WPs (also referred to as the Jacob-Sachs model [49] in [12]) through Eq. (7), resulting in a microscopic uncertainty, σ p , represented in the momentum space. Hence, σ p includes information such as the life-time of the external particles [50] for neutrinos produced by decaying particles, or the mean free path of processes before the production of neutrinos [51]. In addition, regardless of the uncertainties of the external particles, the microscopic space and time uncertainties of the interaction around the vertices are taken into account by another microscopic uncertainty, σ x . This parameter depends on the internal states and the scattering/collision process. In principle, one can also write the first layer Fock-PS directly in terms of the neutrinos, as in Appendix B. However, in this case we can only obtain an effective weighting function in either the energy-momentum space or the space-time coordinate space. We calculate the transition amplitude for neutrinos in the first layer with the Fock-PS representation, by applying the S-matrix method following [12], where the traveling neutrino is treated as an internal propagator in a diagram of a full process including the production and detection process as illustrated in Fig. 2. The kinematics of the neutrino's initial and final state are written in the form of WPs in the momentum space, which will eventually be combined as a weighting function with width σ p . Hence, without loss of generality we can write for the initial/final state at the production site (|P i/f ) and the initial/final state at the detection site 3 , for each h = {q, k, q , k }. Additionally, the internal states (excluding the neutrino propagator) of the process are included by the distributions g P (x 1 ) and g D (x 2 ) which represent space-time uncertainties around the vertex at the production and detection site, respectively. Note that since these states are not restricted on the mass-shell, such uncertainties include four degrees of freedom, namely, a temporal uncertainty and three spatial ones, while the on-shell WPs only have three degrees of freedom. However, such uncertainties are usually not explicitly referred to explicitly in the literature, since in terms of WP separation, we will show that it can be combined with σ p as an effective value, and in terms of a localization term (such as in [13]), it is microscopic compared the scale of the experiment. Nonetheless, since the external particles are better known and are, in principle, observable, one may still be able to extract the contribution of such uncertainties upon measurement. We can readily calculate the transition amplitude of a neutrino that is produced as flavor α while detected as flavor β and propagates in mass eigenstates of mass m j , and integrate over x 1 and x 2 for the sake of completion of our structure as Here M P (q, k) and M D (q , k ) are the plane-wave amplitudes determined by particles involved in the production and detection process, respectively. With the goal of leaving only the neutrino momentum (p = q−k = k − q ) and traveling distance and time (x = x 2 − x 1 ) determined by the entangled states unintegrated, we derive the form of the layer structure as demonstrated explicitly in Appendix B. Here, F (p; P) represents the effective PDF from the WPs of the external particles, and G x (x; X) is that from the vertices. Hence, the layer-moving-operator is where F (p; P )G x (x; X) is the weighting function, and the width of these uncertainties, σ p and σ x , are the observational parameters. From Fig. 18 in Appendix B we see how these parameters are related to each of the original distributions in Eq. (8). Here, the notation G x (x; X) means that X is the expectation value of x for the PDF G x (x), as well as other functions in this paper. In general, the first layer transition amplitude takes the form of Eq. (116), where all configuration of the internal energy of the neutrino propagator has to be included. However, since the measurement is done macroscopically, we can consider the propagating neutrino to be on the mass-shell, hence the first layer transition amplitude for neutrino oscillation is Note that although we impose here the on-shell approximation by hand, it will still be on-shell on the second layer even if we do not. This is due to the fact that the variables on the second layer are macroscopic while those on the first are microscopic. Hence, the dynamics becomes classical and we will obtain the energymomentum dispersion relation. Additionally, note that e ipx in A 1,α→β ensures the uncertainty principle all the way up to the measurement layer, as illustrated previously. Finally, we expand the neutrino energy around p = P j , the saddle point of the total weighting function on the momentum space, then keep terms up to the first order, i.e. E j (p) E j + v j (p − P j ), where E j ≡ E j (P j ) and v j ≡ ∂E j (p)/∂p| p=P j = P j /E j . Hence, corresponding to Eq. (121), we obtain where L and P are the second layer PS variables, such that P j = P j (P) (the explicit relation for the latter will be discussed in Sec. 2.3). In principle, it is possible to do a full analysis of weighting functions from the neutrino decoherence phenomenon, and inspect the mechanism behind it. Nonetheless, this would require very high precision measurements and a wide spectrum in both coordinate and momentum space, and an analysis which is beyond the scope of this paper. Here, we investigate the weighting function through the width of the PDFs, σ x /σ p for G(x, L j )/F j (p, P j ) in Eq. (13) below. If the weighting functions are symmetric, thenΦ j (L, P) is real, according to Appendix A. Otherwise, an additional parameter for the phase ofΦ j (L, P) would be needed. Additionally, note that σ x and σ p are not the total momentum and coordinate uncertainty of the system. In fact, the total coordinate uncertainty is calculated in Appendix B.2, which turns out to be the convolution between the coordinate distributions at the production and detection site. The total coordinate distribution for each production/detection site is the convolution between g P (x 1 )/g D (x 2 ) and the momentum distributions of the external states projected onto the coordinate space (F tot P/D ). In other words, the total coordinate uncertainty is (g P * F tot P ) * (g D * F tot D ), where * denotes convolution of two functions as noted in Table 1. This is illustrated at the vertices in Fig. 2: while the external particles already cause some coordinate uncertainties by conjugating the momentum uncertainties onto the coordinate space (the individual blue and green circles), the internal process would give rise to spatial uncertainties on top of that (the inner solid circle lines) resulting in a total coordinate uncertainty (the outer dashed lines) larger than that of the individual uncertainties.
In particular, for the case of Gaussian WPs, the formalism of F j (p, P) is derived widely in literature, e.g. [9,12,13]. However, in the following we take a generic Gaussian form for the weighting function on the Fock-PS for simplicity, and also integrate out x 0 in Eq. (9) for later convenience. Therefore, from Appendix B.2, where while L j = L − v j T and ∆ = 1 + 4σ 2 x σ 2 p . The relation between σ x and σ p with the original input functions in Eq. (8), namely, f P i , f P f , f Di , f Df , g P and g D is calculated in the Appendix B and summarized in Fig. 18. At the end, σ x includes all the individual uncertainties, temporal and coordinate, production and detection in a convolutional type, hence, the larger one of which will dominate. On the other hand, σ p merges the momentum uncertainties of the production and detection in a product type referred to in Table 1; however, at each site, the total momentum uncertainty combines momentum uncertainties of individual external states in a convolutional way, therefore, the largest one among them would dominate. This relation coincides with the one from Ref. [13] assuming Gaussian distributions. Note that the normalisation of the weighting functions is irrelevant here, since we will later show that the FTP will be automatically normalised on the measurement layer with the definition in Eq. (50). Hence, after inserting the Gaussian distributions, we easily find that Thus, when σ x 1/(2σ p ), i.e. when the width of the inner solid line in Fig. 2 is much smaller than that of the blue/green circle's width, σ x can be neglected, and vice versa.

Microscopic Layer (Layer 1): Quasi-Transition Probability
In this section we calculate the quasi-probability distribution (or Wigner function in e.g. [48]) of the FTP on the Wigner-PS. Such distributions bridge QM (or QFT in this case) to statistical probability distributions, such that the expectation value is calculated by direct PS integration. Moreover, in Section 3.2 we will illustrate that it is insightful and useful to look at state decoherence as a PWO effect from the Wigner-PS perspective. Although we do not obtain the quasi-probability distribution P1 from the Wigner transformation directly, we end up with the same formalism by doing a change of variables such that Eq. (6) is fulfilled. This means we find P1 such that where A 2,j/k are given by Eq. (13). Therefore, for Eq. (16) implies that in which the right-hand side involves mixing of the two PS, (x, p) and (x , p ), while that on the left-hand side does not. Therefore, P1 ,jk is the quasi-probability distribution representing the occupation number of having both the jth and the kth mass eigenstate simultaneously in the Wigner-PS. The equation above is achieved when whereW G jk (x,p) andW F jk (x,p) take the form of the Wigner quasi-probability distribution as follows: In this case LMO1 = d 3x d 3p is simply the integration over the Wigner-PS, and the FTP on this layer is the quasi-probability distributionW G jkW F jk , which includes both W 1 and B 1 in Eq. (2). In particular, if we assume all weighting functions to be Gaussian on the Fock-PS given in Eq. (14), the quasi-probability distributions becomeW . When we move the FTP up to the physical layer by integrating over the Wigner-PS, there will be a PWO effect suppressing the plane wave term on the physical layer. Moreover, the PWO effect is determined by the width of the weighting function relative to the wavelength of the phase structure, which is P j − P k forx and L j − L k forp in this case. In fact, the wider the weighting function is relative to the wavelength of the phase structure, the more suppression will the PWO effect cause. We can also view this as how many periods (e.g. 2π/(P j − P k ) is one period inx) there are determined by the phase structure within some width of the weighting function. We call the number of periods within some area on the PS the phase density. Hence, the higher the phase density is, the more damping we will get from the PWO effect. Fig. 3 and Fig. 4 are plotted as an example to illustrate the quasi-probability distribution on the Wigner-PS before integrating out the PS, which would lead to a PWO effect on the physical layer after the integration. In particular, Fig. 4 demonstrates the projection of 3D plots such as Fig. 3 onto a 2D contour plots,  where the areas within both the shaded and non-shaded circles represent contour lines of one standard deviation of a Gaussian distribution with Eq. (21). For both plots, the more times we see positive and negative values (red and blue circles) alter, the higher the phase densities are within one standard deviation range of the weighting function, and the larger the PWO effect will be. In particular, relative to the left plot which is taken for some time T left and some width σ p = σ x = σ left for the weighting functions, the middle plot has the same width but a larger propagating time, i.e. T middle > T left . Furthermore, the right plot has width as either σ p = σ x = σ right = 0.325 σ left or 1.376 σ left , and T middle = T right . By comparing the left and middle plot, we see that on this layer, as time evolves, although the FTP of the two mass eigenstates separate, the width of the overlapping FTP does not. However, the phase density does increase, hence, so does the PWO effect. On the other hand, since the width of the FTP is σ 2 This indicates that the quasi-probability is localised and the uncertainty principle holds automatically. Hence, there are two solutions of σ right , for some value of σx = σp < 1/2. Also, the decrease in σx and σp would also reduce the phase density, relaxing the PWO effect. Note that on the Wigner-PS, when the quasi-probability distribution is localised to a point, it describes a classical monochromatic field [48].

Physical Layer (Layer 2): Transition Probability
In this section, we move onto the physical layer, where experimental uncertainties are considered in addition. Since we consider isotropic emission of neutrinos for simplicity, the additional uncertainties introduced through the weighting function on this layer are the (macroscopic) energy uncertainty σ E and the (macroscopic) coordinate uncertainty σ L . The former include the energy uncertainties, such as the energy resolution of the experiment (σ E ∝ 1/ √ E 0 ) and the energy reconstruction models (for example, see [52,53]); whereas the latter include any uncertainties on the propagation distance. For instance: the core size and the distribution of multiple reactors for reactor neutrinos; the length of the decay pipe and the velocity of the parent particles before they decay into neutrinos for accelerator neutrinos; the uncertainty in the altitude of where the neutrinos are produced in the atmosphere for atmospheric neutrinos; or even more exotic effects that would introduce uncertainties to how long a neutrino propagates before it gets detected, such as space time fluctuations.
Before taking the additional uncertainties on layer 2 into consideration, we need to first move the FTP from layer 1 to this layer, which can either be transferred from the first layer by the Fock-PS via where the FTP is or by the Wigner-PS via Both approaches give us the same result as a useful consistency check. By taking the Gaussian distribution on the first layer in Eq. (14), we arrive at by either Eq. (15) or Eq. (21), and it is also useful to rewrite the widths in terms of σx = ∆ σ x and σp = σ p /∆, where ∆ = 1 + 4σ 2 x σ 2 p . So far we have not specified how the second layer momentum variable P in the relativistic-PS is related to the expectation value of momentum P j for each mass eigenstate. The relation is obtained below by comparing the massless case with the small mass case. For a certain energy E, in the massless case, we have E = |P|; but for nonzero neutrino mass, some of the energy, δE j , will be consumed by the mass, and this freedom is limited by the uncertainty of the energy on the first layer. In other words, this picture is equivalent to having some P j (E ) instead of E j (p) while deriving Eq. (12), with E being the mean of energies E of the neutrinos given by the external particles in the first layer. Therefore, we can write If we expand P j with respect to m j and keep only the leading order, we get where ξ p = P/E = P/|P|, so that | ξ p | 2 = | ξ j | 2 = 1. This allows the identification of δE j , by substituting to lowest order in m j . The actual energy which the mass eigenstate carries is decided by the dispersion relation: Therefore, we can see when ξ p ξ j = 1, i.e. the case where all mass eigenstates as well as the massless case are co-linear with each other, we have equal energy of mass states; and when ξ p ξ j = 0, we have equal momentum modulus instead. To have exact equal momentum, we need ξ j = 0. However, according to e.g. [54,55], neither of this should be the case due to Lorentz invariance. Additionally, it is also useful to derive the group velocity With the approximated relation from above, Eq. (25) becomes where the phase structure, the momentum weighting function and the coordinate weighting function are respectively. Also, we have taken ξ j = ξ k ≡ ξ and η = ξ ξ p = ξ ξ L as the alignment factor. Note that in Sec. 3 we will show that terms such as the second term in the brackets of Eq. (33) will be cancel out by normalization. Furthermore, for experiments with continuously emitted neutrinos within a sufficiently long period of time, we should integrate out T , since there is no temporal information. In this scenario, ψ jk would be replaced by ψ jk and D x would be replaced by D x in Eq. (31), where Note that ψ jk being independent of ξ j and ξ k implies that "equal energy", "equal momentum" or anything in between, will lead to the same phase structure if we have no temporal information. In fact, when σx ∼ 0, we will arrive at the standard decoherence formula in [13], namely where for freely propagating neutrinos.
Up to now, we have not yet considered the weighting functions on the second layer, but only how the weighting functions on the first layer are transmitted onto the second layer. But before going into that, we simplify our discussion by changing the variables {T, L, P} to {T, L = |L|, E = |P|, Ω L , Ω P }, where Ω L /Ω P are the solid angles for L/P. Therefore, the layer-moving-operator from the physical layer to the measurement layer is Then by assuming that neutrinos are isotropically emitted and that we have no temporal information, i.e. we consider only uncertainties of E and L, P 2,jk is moved to the third layer as For a counting experiment, the transformation from the second layer to the third layer performs the convolution between the FTP on the second layer which includes uncertainties from the first layer, and the weighting function H X (X; X 0 ), which is the PDF of the true value X, for some measured value X 0 . For example, the measured rate at L 0 does not only have contributions from neutrinos actually propagating the distance L 0 but is a sum of all possible contributions within a time window given by the uncertainty of time, which is taken as infinity for the case of continuous emission of neutrinos for a sufficiently long period of time. In particular, with the commutative and associative properties of convolution, the coordinate uncertainty, H L (L; L 0 ), comes from the convolution of the spatial PDF of the production process and the detection process. Therefore, the largest uncertainty would dominate, which is usually the production PDF, i.e. the source profile or the PDF of neutrino production. Since Eq. (40) also results in a PWO effect, Fig. 5 is plotted in the same way as Fig. 4 to see the phase density of the weighting function. The difference between these two plots is that on the relativistic-PS, the separation of P 2,jj (L, P) and P 2,kk (L, P) affects the width of P 2,jk (L, P) (black dotted circle), unlike the case of the Wigner-PS in Fig. 4. In Fig. 5 we illustrate the time-dependent case for P 2,jk on the physical layer for an example, the phase structure is the same as Eq. (32), as well as the spatial uncertainty with Eq. (34). However, since the energy range of interest is much larger than the central value of D p which is suppressed by the neutrino mass, in Eq. (33) we consider the second layer weighting function H E (E; E 0 ) for the energy uncertainty. Hence, the contour lines in Fig. 5 are for one standard deviation of We also show two cases of the alignment η between ξ p and ξ = ξ j = ξ k in Eq. (27), to see how it affects the time-dependent phase structure. In addition, the heavier mass eigenstate P 2,jj (the yellow shaded area) is a bit tilted at large T , due to the energy dependence of the group velocity in Eq. (30). Nonetheless, this tilting is negligible in reality, since the mass splitting is a lot smaller compared to the uncertainty of the energy, σ E . We will show that Fig. 4 and Fig. 5 correspond to state decoherence and phase decoherence, respectively, and illustrate the time-independent version of the plots in Sec. 3.

Measurement Layer (Layer 3): Transition Probability
In this section, we finally reach the measurement layer where experimental data is collected. On top of the uncertainties arising from for the previous layers, i.e. the phase space uncertainties (PSUs), the final data also has to take into consideration the count uncertainties (CUs). While the PSUs are uncertainties of the PS variables (for instance, σ x , σ p , σ L and σ E discussed previously), the CUs are uncertainties of the neutrino flux, for example, the statistical uncertainties and the background uncertainties. The PSUs and the CUs will be treated differently in our likelihood/χ 2 analysis in Sec. 4: the CUs will be treated as the conventional "uncertainties" in the analysis, whereas the PSUs are included in the theoretical prediction. The total count rate for energy-distance binned data is where the CUs are included in Φ 0 (E 0 ), which is the neutrino flux at L 0 = 0 taking into account the decay rate of the production process; D(E 0 ) is the detection rate including the cross section involved in the detection process. On the other hand, the PSUs are included in P 3,α→β (L 0 , E 0 ), the FTP on the measurement layer, and will be discussed in the following paragraph. Below, we demonstrate effects of PSUs on the FTP by integrating over the approximated time-ignorant FTP on the second layer, and take Gaussian PDFs for both H L (L; L 0 ) and H E (E; E 0 ). However, since only the integration over a Gaussian H L , can be calculated analytically, the integration over E will be done numerically, and the FTP on the measurement layer, P 3 (L 0 , E 0 ) is plotted in Fig. 6 and Fig. 7. In Eq. (44), the total spatial uncertainty (width of H L ) for a vacuum propagating neutrino is σ 2 Here σ S describes the production profile at the source and σ D represents the spatial resolution of the detector, hence, σ L is normally dominated by σ S , as mentioned previously, and In the second line of Eq. (44), we can see that the two last terms in the exponent are locality terms, indicating that the more local uncertainty (σx, 1/σp and σ L ) we have, microscopic or macroscopic, the more the FTP will be smeared out, and since σx and 1/σp σ L the macroscopic one dominates as we can see in the third line. As for the first term in the second exponent, the coherence length is modified by L coh kj → (L coh kj ) 2 + 4σ 2 L , but since L coh kj is inversely proportional to the mass splitting of neutrinos, it is usually much larger than σ L for ground-based experiments. In this case, the microscopic uncertainty would dominate for this term. However, this would not be the case for neutrinos produced in continuously emitting celestial objects for a long period of time, see e.g. [56]. For instance, the size of the Sun's core would be much larger than the coherence length in the three neutrino paradigm. Nonetheless, the coherence length might be stretched out for neutrinos produced in an extreme environment, such as supernovae [57,58], or by going to ultrahigh energy for long oscillation length (∝ 1/E 0 ) [59]. Additionally, since the integration dEH E (E; E 0 )P semi 3,jk (L 0 , E) satisfies the factorization condition in Appendix C, we can write Thus, the effects of H L barely depend on L 0 , unless σ L accumulates w.r.t. L 0 , since it would then only affect the locality term. The FTP of neutrinos detected at a longer distance is plotted in Fig. 6 and the left plot in Fig. 7, with oscillation parameters taken from NuFit 5.1 global fit results [36]. From these figures we can see that the FTP is barely sensitive to σ L compared to σ p and σ E , and σ E would not only cause damping to the oscillation, but also a phase shift, which is more pronounced at low energies. On the other hand, the right plot of Fig. 7 is plotted for a very short traveling distance, and the sensitivity for σ L is much larger than that for σ p and σ E . Although the effect is very small here compared to the long distance case, it benefits from larger statistics. Moreover, since the phase-density is higher for lower energies, due to a smaller oscillation length, the PWO effect is more significant and the coherence length is shortened, hence neutrino decoherence effects are more enhanced for lower energies, as we can see in Fig. 6 and Eq. (44).

Formalism
In this section we show how neutrino decoherence can be further classified into two categories, depending on either coherence is lost by the separation of the mass eigenstates, or by statistical averaging. Moreover, in the language of the layer structure, we will see that both categories of decoherence effects result from PWO effects, therefore, we formulate the effects of neutrino decoherence with a damping term (φ jk ) and a phase shift term (β jk ), which can be written as for   where ψ jk is the phase structure on the second layer in Eq. (32). Also, σ = { σ x , σ p , σ L , σ E , σ T } represents the set of parameters that describe the weighting functions, such as the width and asymmetry parameters for the respective variables. In general X 3 = {T 0 , L 0 , E 0 , Ω L0 , Ω P 0 } are the third layer's temporal variable, spatial variable, energy variable, spatial solid angle and momentum solid angle. As for isotropic neutrinos X 3 = {T 0 , L 0 , E 0 }, and on top of that, if the neutrinos are continuously emitted for a sufficiently long period of time, X 3 = {L 0 , E 0 }. In this paper, we will only consider the two latter cases for simplicity. In particular, if all the weighting PDFs are Gaussian distributions, then the damping term can be parameterized as Moreover, since P 3,jk is on the measurement layer, it is necessary to properly define the operational FTP.
Here, we define the FTP by each P 3,jk (X 3 ), as the ratio between the total count with and without oscillation, namely where Γ 2,jk (X 2 ; X 3 ) ∝ P 2,jk (X 2 ; X 3 ) is the un-normalized FTP on the second layer. Here X 2 are the second layer variables analogous to X 3 . We note that the so-defined P 3,jk (X 3 ) is not affected by the scale of the weighting functions on each layer, but is only dependent on their widths and shapes.
One might be concerned that the definition in Eq. (50) is not justified since we measure the FTP by the total P 3,α→β instead of each P 3,jk . However, theoretically speaking, it is possible to measure P 3,jk if we measure the denominator by measuring the exact neutrino mass, and the numerator from a well controlled oscillation experiment. Note that there would be no oscillation for a mass measuring experiment, for one would know exactly which mass eigenstate the neutrino propagates on. Therefore, a mass measuring experiment is only useful for determining the normalization of an oscillation experiment, whereas the numerator in Eq. (50) can be found in oscillation experiments for neutrino mixing where only two mass eigenstates (with eigenvalues m j and m k ) are allowed/sensitive. Nonetheless, with the smallness of the mass splitting, the shape of the weighting functions on each layer can be approximated as being independent of j, k. Therefore is the production rate for flavor α, and σ det β (X 2 ; X 3 ) is the detection cross section for flavor β. In this case, the difference between the weighting functions only comes from the different group velocities v j , which is only relevant when j = k, hence, the FTP would become This expression and its condition for the uncertainties to be independent of j, k coincide with [9]. Additionally, the definition of Eq. (50) automatically normalizes, without any further condition, since it implies P 3,jj = 1 for any j, then and similarly for β P 3,α→β = 1. In fact, for 0 ≤ U βj U * βj ≤ 1 and 0 ≤ P 3,jj ≤ 1 ∀j, if and only if P 3,jj = 1 will the normalization condition be satisfied.
Furthermore, the definition of Eq. (50) also serves the purpose of analyzing the decoherence effect, by writing it as State decoherence term Phase decoherence term × × e iψ jk P 3,jk = Figure 8: Demonstration of how we separate the FTP defined in Eq. (50) into two terms: the state decoherence term and the phase decoherence term in Eq. (55) and Eq. (56) in terms of probability. The blue and red shaded circles represent two different mass eigenstates on the physical layer, and while the state decoherence term represents the separation of the two mass eigenstates, the phase decoherence term demonstrates a phase wash-out effect. where and Here, we have written for simplicity. Both functions in Eq. (55) and Eq. (56) are the decoherence terms, which are both unitary for the fully coherent case, and with modulus ≤ 1 in general. In particular, S 3,jk represents the probability of the two mass states overlapping with each other, while Φ 3,jk gives us the PWO effect introduced in previous sections and in Appendix A. Therefore, we call the former part "the state decoherence (SD) term", where S 3jk (X 3 ) = 1 indicates that two mass eigenstates are fully overlapping. The latter part will be called "the phase decoherence (PD) term", where Φ 3,jk (X 3 ) = 1 represents the case where there is no PWO effect on the physical layer. The trick to separate the decoherence term into SD and PD in Eq. (54) is demonstrated in Fig. 8, where we insert 1 = dX 2 H|Γ 2,jk | e −iψ jk / dX 2 H|Γ 2,jk | e −iψ jk into Eq. (50). The blue and red shaded circles represent two different mass eigenstates on the physical layer, and the total FTP would be the purple area of the numerator in the phase decoherence term normalized by the purple area in the denominator of the state decoherence term. Finally, we see that SD is the probability of the two mass eigenstates being in superposition with each other on the relativistic-PS, while PD indicates our ignorance to the system on the relativistic-PS by averaging over all possibilities. In particular, corresponding to Eq. (44), the term with the coherence length describing the WP separation represents SD, and the localization term shows PD. Therefore, we can already expect that the SD term will be dominated by the (microscopic) uncertainties on the first layer, while the PD term will be dominated by the (macroscopic) uncertainties on the second layer.

State Decoherence
In this section, we will present how the SD term can be further analyzed and approximated. Eventually, we will show that although the SD represents state separation on the physical layer (2nd layer) it is equivalent to a PWO effect on the Wigner-PS (1st layer) under most conditions. Therefore, the dominating uncertainties for SD are those on the first layer, namely σ x and σ p .
• σ x (coordinate uncertainty on the Fock-PS layer 1): This uncertainty originates from the off-shell intrinsic uncertainties related to the finite space-time extension of the vertices, which depends on the neutrino interaction at the production and detection site. In fact, the Fourier transformation of the weighting function g P /g D results in an effective form factor as a function of neutrino momentum, which can be seen by integrating out x 1 /x 2 in Eq. (8). Hence, depending on the interaction, σ x could be the size of the charge radius of the proton or neutron, which is at O(0.1 − 1) fm according to [60], or at most the inter-atomic distance at O(0.1 − 1) nm, which is discussed in [61]. In Sec. 2.1, we included uncertainty sources as the unrelated spatial and temporal uncertainties (since they are not restricted to the mass-shell) for both the production and detection site in the calculation of the transition amplitude. In Appendix B, we have shown that σ x collects the sources of uncertainties in a convolution way, hence, it is dominated by the largest uncertainty source.
• σ p (momentum uncertainty on the Fock-PS layer 1): This uncertainty arises from the WP description of the external states on a quantum mechanical level, which also appears when we calculate the transition amplitude. For instance, σ p would be related to the mean free path of interactions before the external particles interact with the neutrinos or the life-time of the parent particles, see e.g. [13,50]. Moreover, since the external particles are on the mass-shell, the energy uncertainties are correlated to the momentum uncertainties, and unlike σ x , σ p collect the sources of uncertainties, including the production's and detection's momentum and energy uncertainties, in a product way, hence the smallest one among them would dominate. However, the total momentum uncertainty at each site is the convolution of the initial and final state WP, hence, the smallest one among them would dominate (see Appendix B).
Note that σ p and σ x and independent of each other and can both be represented in either coordinate space or momentum space. The difference between these two uncertainties is whether they describe the uncertainties from the external states or the effective vertices (internal states). Moreover, from Sec. 2.3 and Sec. 2.4, we have seen that it is more effective to consider the width of the quasi-probability distribution: σx = σ x /∆ and σp = σ p /∆, where ∆ = 1 + 4σ 2 x σ 2 p . • σx (coordinate uncertainty on the Wigner-PS layer 1): By Eq. (66), the damping term in Eq. (49) for a time-independent Gaussian distributed PDF is However, as we will see in the next section, this structure is exactly the same as that for σ L , which is usually macroscopic, and can therefore be neglected.
• σp (momentum uncertainty on the Wigner-PS layer 1): By Eq. (36), the damping term in Eq. (49) for a time-independent Gaussian distributed PDF is and the time-dependent one can be taken care of by replacing L 0 with T 0 . Therefore, the decoherence effect resulting from σp is should be searched for at long distance. Such effects have been widely studied, and explored, e.g. for reactor neutrinos in [37,38], which excludes (2σp) −1 < 2.08 × 10 −4 nm at 90 % CL. In fact (2σp) −1 represents the total coordinate uncertainty combining the production and detection region. Therefore, in the following, we sometimes take σp around 0.1 MeV as a "reasonable" value, since it corresponds to (2σp) −1 ∼ 10 −3 nm, which is within the range of the proton/neutron charge radius and the inter-atomic distance (O(0.1-1) nm), see [61] for an estimation for reactor neutrinos.
For the derivation of SD, we start by writing down the SD term following Eq. (55) as by replacing Γ 2,jk with where S 2,jk = S 2,jk (X 2 ), Φ 2,jk = Φ 2,jk (X 2 ) and ψ jk = ψ jk (X 2 ). Moreover, we have ψ jj (X 2 ) = 0 for any j according to Eq. (31), and Γ1 ,jk (x,p; X 2 ) ∝ P1 ,jk (x,p; X 2 ) is the un-normalized transition probability distribution on the Wigner-PS, such that Γ 2,jk (X 2 ) = d 3x d 3p Γ1 ,jk (x,p; X 2 ). Therefore, the remaining term becomes indicating that S 2jj (X 2 ) = 1, since ψ jj = η jj = 0 ∀j. Hence, we only need to consider thex andp dependent terms forD jk . In fact, Eq. (61) has the same formalism as Eq. (96) (i.e. it is a PWO effect) with the phase structureη jk averaged over within the normalizedD jk region. For instance, considering Gaussian distributed weighting functions (Eq. (21)), then and the phase structure is in general where the relation between T, L, E andL jk ,P jk , P j , P k is given in Sec. 3. In the following, we again consider all distributions as isotropic, and simplify our discussion to a one-dimensional scenario. From Appendix C, we find that for σ S and ∆ X 2 being the width of S 2,jk (X 2 ) and Y 2,jk (X 2 ; X 3 ) ≡ H(X 2 ; X 3 )Φ 2,jk (X 2 ) respectively, the first approximation in Eq. (64) can be made in the limit of X 3 ∆ X 2 and σ S ∆ X 2 . The picture of this factorisation condition is to neglect the uncertainty of X 3 on the second layer for the term that can be taken out of the integral, i.e. X 2 X 3 , in S 2,jk for some ∆ X 2 , but at the same time σ S can not be neglected, in order to have SD. Additionally, if the total uncertainty ∆ X 2 is dominated by the physical layer uncertainty H, which is exactly the case for macroscopic measurements, then Y 2,jk (X 2 ) Y 2,jj (X 2 ) Y 2,kk (X 2 ) and we arrive at the second approximation in the following equation: The evaluation of the comparison between the width sizes of S 2,jk , Φ 2,jk and H can be done by taking Gaussian distributions for weighting function on each layer. Precisely speaking, for both S 2,jk (with width σ S,X 2 ) and Φ 2,jk (with width σ Φ,X 2 ), we use Eq. (14); as for H(X 2 ) (with width σ H,X 2 ), we simply write it as a Gaussian distribution around X 3 . In this case, we have: and Therefore the width of Φ jk w.r.t. L and T is σ Φ,L = (2σp) −1 and σ Φ, ] −1 and σ S,L → ∞ for the time-dependent case, since S 2,jk does not depend on L. As for the energy uncertainties, we have and State decoherence Phase decoherence With the 1/E dependence, the exponent in Eq. (68) → 1 as E → ∞, so σ S,E diverges, and the width σ Φ,E would also be stretched out. Therefore, σ H,L σ Φ,L , σ S,L and σ H,E σ Φ,E , σ S,E . Hence, the macroscopic uncertainty σ H,L , σ H,E will dominate over the transferred microscopic ones of Φ 2,jk for Y 2,jk (X 2 ; X 3 ) in Eq. (64). Thus, ∆ L σ H,L and ∆ E σ H,E . Accordingly, the factorization condition for Eq. (64) is fulfilled for the L and E part if the measurement values L 0 and E 0 are also much larger than σ H,L and σ H,E , which is exactly the case for neutrino experiments which have the resolution to measure neutrino oscillation.
As for the temporal part, we will discuss two scenarios: σ H,T σ Φ,T , σ S,T and σ H,T → ∞. If it is neither of these two cases, one should integrate over T in advance while taking H T into account. In the former case, the factorization condition for T is also satisfied, and H(X 2 ) dominates over Φ 2,jk (X 2 ) for all X 2 = {L, E, T }, then S 3,jk (X 3 ) S 2,jk (X 3 ). Therefore, we can directly obtain the observational time-dependent SD weighting function D jk (x,p; T 0 , L 0 , E 0 ) and phase structure η jk by replacing X 2 with X 3 . Specifically, if all quantum uncertainties are Gaussian distributed, then D jk =D jk (x,p; T 0 , L 0 , E 0 ) and η jk =η jk (p,x; T 0 , E 0 ). In Fig. 4, we illustrate such time-dependent PWO effect which increases as the WP of two states separate with time. In fact, this would be translated into WP separation on the physical layer in the sense of Eq. (59) and Fig. 8. On the other hand, if we have no temporal information during the detection process, i.e. σ H,T → ∞, then we should integrate out X 2 = T in Eq. (59) first, then look at the SD term in the same form of Eqs. (59)-(64) but with X 3 = {L 0 , E 0 }, X 2 = {L, E} and Γ 2,jk (X 2 ) replaced by dT Γ 2,jk (L, E, T ). In this case, after the approximation given in Eqs. (26)- (29) and including terms up to O(m 2 ), the Gaussian example in Eq. (14), leads to the left plot in Fig. 9, and the time-independent SD weighting function is In the left panel of Fig. 9, where the PWO effect for the above phase structure and weighting functions is shown, we take η = 1, σx = 0 and amplify ∆m 2 jk to one order smaller than E 0 for illustration purpose. Nonetheless, when ∆m 2 jk E 0 , the integration overx would only have negligible contribution to the S 3,jk , and the SD term becomes where D jk,σp is a product of three Gaussian distributions. Moreover, for E 0 m 2 j /E 0 , the dominate one would be the distribution with width σp/2 and centred at E 0 , while the phase structure appears to be η jk −∆m 2 jkp L 0 /(2E 2 0 ). This can be seen from the upper row in the left panel of Fig. 9, where the phase averaging mainly results from the integration overp. Hence, the SD is mainly decided by the Wigner distribution w.r.t.p, which we will call Dp-induced decoherence. In fact, the resulting SD term, S 3,jk , from such PWO effect not only agrees with the standard decoherence formula in Eq. (37), but also shows that such L 0 and E 0 dependence is a consequence of the phase structure. Moreover, although we demonstrated the case where all weighting functions are Gaussian distributed, the phase structure also applies to arbitrary distributions if the saddle point approximation is adopted in Eq. (20). Therefore, as the phase structure implies a Fourier transformation fromp to α jk = ∆m 2 jk L 0 /(2E 2 0 ), we plot the damping term and the phase shift term from SD for some typical distributions in Fig. 10. In particular, while the one Gaussian case represents a typical statistical distribution for a single process, the two-Gaussian case considers neutrinos produced simultaneously by two different processes with slightly different expectation values for momentum. The other two distributions are more suitable for describing an H L -induced PD effect which will be introduced in the next subsection, where more discussions on this plot will be given along with the H E -induced PD effect. Nonetheless, since both decoherence effects are described by Fourier transformation, only with different space mappings, we demonstrate both effects in the same plot. In fact, this shows that different sources of decoherence effect are distinguishable by their (L 0 , E 0 ) dependence, which origins from the difference in their phase structures.
To sum up, while the SD represents the separation of two mass eigenstates on the physical layer, which is quantified by the overlapping area between them (Fig. 8), it is not the case if we move down to the Wigner-PS. From Fig. 4, we can see that the width of P1 ,jk does not get smaller as the two mass eigenstates, P1 ,jj and P1 ,kk separate. Nevertheless, the phase density increases as the two mass eigenstates depart from each other, and the PWO effect is stronger, which would give the same result as calculating the amount of overlap of two mass eigenstates on the physical layer. In addition, from the colorful background in Fig. 9, we see that the phase structure would vary with PS variables on the third layer, T 0 , L 0 and E 0 , while the width of the overlapping weighting function does not. For the time-independent case on the left plot, where σ H,T → ∞, it is straightforward to see that there is no dependence on T 0 ; yet, for the time-dependent case in Fig. 4, it turns out that there is no dependence on L 0 . This is because in the difference L j − L k the common factor L cancels, which can be seen from either the phase structure in Eq. (21), or directly from Eq. (66).

Phase Decoherence
As we will show in this section, the dominating uncertainties deciding the effect of PD are the macroscopic ones summarized below: • σ L (coordinate uncertainty on the relativistic-PS layer 2): This uncertainty mainly comes from the macroscopic spatial uncertainty of the full process, which is dominated by the uncertainty on  Figure 10: Decoherence damping terms (second row) and phase shift terms (third row) for different shaped weighting functions (first row) of H L -induced or Dp-induced decoherence. While both types of decoherence are described by Fourier transformation, the former transfers from L space to the α jk = ∆m 2 jk /(2E 0 ) corresponding to Eq. (76); and the later fromp to α jk = ∆m 2 jk L 0 /(2E 2 0 ) relates to Eq. (71). In particular, the same coloured lines represent weighting functions with the same widths. the production profile of the neutrino source for the vacuum propagation case. Unlike σ x , σ L does not enter the Feynman diagram (Fig. 2) which calculates the transition amplitude. Instead, it accounts for uncertainty in the traveling distance of the transition probability on the second layer. For instance, σ L could be dominated by the reactor core size (∼ 3 − 5 m) for reactor neutrinos and the distance mesons/muons travel before they decay into neutrinos in accelerator experiments. For Gaussian distributed PDF, the damping term in Eq. (49) is to a good approximation due to the smallness of mass splitting when we integrate out E in Eq. (44). In other words, it is nearly independent of the traveling distance. Hence, it is favourable to search for such effects close to the neutrino source for the sake of higher statistics. In Fig. 10, we show the effects of PD for different production profiles: The one Gaussian PDF can be used for neutrinos produced at rest such as reactor neutrinos and DAR neutrinos. The two Gaussian PDFs are shown for the scenarios with multiple sources/detectors. The box PDF should be adopted when constraints put forth by the experimental setup dominates. For instance, the pipe/rod volume of some accelerator/reactor producing neutrinos, which cuts the production profile to an idealized box shape. Finally, the exponential decaying PDF is for neutrinos produced by decaying particles decaying at flight, for instance, accelerator neutrinos and atmospheric neutrinos. However, if the propagation process also contributes to σ L , then σ L will accumulate over distance, such as through matter effects [46], or some exotic effects [18,23,24,59]. In this case, we can write σ 2 L ∝ L 0 , to agree with the dependence on the traveling distance in the damping term calculated by the Lindblad equation, which the literature mentioned above adopts.
• σ E (energy uncertainty on the relativistic-PS layer 2): This uncertainty is mainly governed by the energy resolution and reconstruction model of the experiment. Typically, for neutrinos detected with photomultiplier tubes, the energy resolution is given by Additionally, there will also be different contributions due to the energy reconstruction model, for instance, the degree of quasi-elastic scattering in [52,53], leading to a tail in H E . Moreover, due to the 1/E dependence in ψ jk , which makes σ E asymmetric w.r.t. the phase structure, there will be a phase shift even for a Gaussian distributed PDF, as we have plotted numerically in Fig. 11. Then from such numerical results, we can extract an E 0 and L 0 behavior. For example, we can see from Fig. 6 that the effect of σ E also increases with L 0 in a comparable way as σp.
Similar to state decoherence, phase decoherence also describes a PWO effect, which can be seen directly from its definition: Here is a real and normalized PDF, which is dominated by the second layer macroscopic weighting function H X 2 for X 2 = L, E. When σ H,T σ S,T , σ Φ,T , Y 2,jk is completely dominated by the macroscopic weighting functions, in this case the phase structure will be On the other hand, if we have no temporal information, we should again integrate over T before looking at the decoherence effect. The phase structure up to O(m 2 ) for the time-dependent scenario is Eq. (32), while the one for the time-independent one is Eq. (35). From Eq. (44) we can see that after L is integrated out, the factorization condition in Appendix C is satisfied, such that we can factorize Eq. (73) as for the time-independent case, where both H L (L; L 0 ) and H E (E; E 0 ) are normalized PDFs. Therefore, we can treat the (macroscopic) coordinate and energy uncertainties on the second layer as separate PWO effects, naming the former H L -induced PD, and the latter H E -induced PD. In Fig. 10 and Fig. 11 these two PWO effects are demonstrated, respectively, where the weighting functions, H L and H E are taken as PDFs with the same width for the same colored lines. In other words, dX 2 H X 2 = 1 and max{H X 2 (X 2 )} = 1/(σ X 2 √ 2π), according to our definition of "width" in Appendix A. Therefore, the only parameter for the two figures presenting the PWO effect on the second layer is the width σ X 2 , for some L 0 or E 0 , which we take σ L = 1 m for the blue and yellow line, σ L = 2 m for the red line in Fig. 10, at distance L 0 = 10 m; also, MeV for the red line in Fig. 11, at energy E 0 = 10 MeV. In particular, analogously to the Dp-induced SD, the H L -induced PD takes the form of a Fourier transformation from L to α jk = ∆m 2 jk /(2E 0 ), as plotted in Fig. 10. For instance, the Gaussian PDF transforms into a Gaussian distribution in the ∆m 2 jk /E 0 space, the box PDF transforms into a sinc function, and the PDF for exponential decay (for neutrinos produced by decaying charged leptons) is transformed into a Lorentzian function. As for the phase shift term, from Appendix C we know that only the asymmetric functions (the yellow-lined two Gaussian PDF and the exponential decaying PDF) have non-zero and non-π phase shift. For the symmetric ones, while there is no phase shift for a single Gaussian PDF, the phase ranges would jump from 0 to π (still no imaginary part) for the box PDF and the symmetric two-Gaussian PDF, due to negative values of the function from the Fourier transformation. However, owing to the smallness of neutrino mass splitting, α jk is typically small for both Dp-induced and H L -induced decoherence effects. Hence, the range of interest would lie in a small range around α jk → 0, where the phase is zero for asymmetric cases. On the other hand, regardless of how "symmetric" the weighting PDF spectrum is, it is not symmetric w.r.t. 1/E. Therefore, there is always a non-trivial H E -induced phase shift, even for the Gaussian distributed PDF. Moreover, it is also clear that the larger the width of the PDF is, the more significant the decoherence effect (both the damping term and the phase shift term) becomes. In the end, if we find some phase structure dependence as the damping term and/or the phase shift term in Fig. 10 and Fig. 11, we should be able to reconstruct the production profile and cross-check the energy reconstruction model.
Comparing PD with SD, both come from PWO effect, but with different phase structures (η jk /ψ jk for time-dependent SD/PD and η jk /ψ jk for time-dependent SD/PD) and distributions (D jk / H T H L H E for time-dependent SD/PD and D jk Dp / H L H E for time-dependent SD/PD). In particular, for SD both the phase structure and the distribution varies with L 0 , E 0 , while in the PD case, only the mean of the distribution does. Such structure is illustrated in the right plot of Fig. 9, where the phase structure in the background only depends on the second layer PS variables, L and E, but not on the third layer PS variables, L 0 and E 0 . This also implies that the second and third layer share the same coherent phase structure, but not with the first layer. The third layer PS only decides where the weighting functions are centred, and when it appears at a higher phase density area (lower energy or larger distance), there will be a stronger PWO effect. Illustrating the PWO effect for the time-independent case, Fig. 9 is plotted in the same ways as Fig. 4 and Fig. 5, while also showing the phase structures on the Wigner-PS and the relativistic-PS, respectively, in the colourful background. The final observational effects on the measurement layer of SD, H L -induced PD and H E -induced PD are further plotted in Fig. 10 and Fig. 11.
The advantage of putting SD and PD effect in terms of PWO effect is that one can estimate the decoherence from S 3,jk and Φ 3,jk numerically, and the only input would be the weighting functions. Therefore, the steps to estimate the SD/PD effect numerically are simply: 1) Determine the weighting functions on the Wigner-PS and relativistic-PS; 2) Multiply it with the corresponding phase structure derived in this section; 3) Integrate out the respective PS. Such integration are certain to converge, since the weighting functions are localised and distributed around the next level PS variables. For instance, even for the simplest phase structure -the time-independent case on layer 2 -the PD terms can only be evaluated numerically as we did for Fig. 11. Therefore, in principle, by analysing the waveform and spectrum of neutrino oscillation, we should be able to reconstruct the weighting functions once we identify its corresponding phase structure.

Phenomenology of Neutrino Decoherence
From the previous section, we saw that neutrino decoherence effects can be classified into SD and PD, which are both a consequence of the PWO effect. Therefore, both would result in damping terms and phase shift terms, where the latter is nontrivial only when the weighting function is not symmetric w.r.t. the phase structure. In particular, SD is dominated by the (microscopic) uncertainties on the first layer (σ x and σ p ), while PD is governed by the (macroscopic) uncertainties (σ T , σ L and σ E ) on the second layer. Nonetheless, we can only observe σp for the first layer uncertainties since σx σ L . In addition, we only consider the time-independent case in this section, since current experiments shown in Fig. 12 continuously emit neutrinos for a sufficient long period of time, therefore, σ T → ∞. Additionally, we also parameterize the first layer uncertainties as σx and σp, then only σp would be an valid observational parameter, while σx would be eaten by σ L . In this section we estimate how far we are experimentally from having a 90% CL sensitivity for observing damping and/or phase shifting signatures in Eq. (47) from the Dp-induced SD, H L -induced PD and H E -induced PD.
For the damping signatures, we consider all weighting functions distributed as single Gaussians. Hence the damping term is then parameterised by σp, σ L and σ E in where the first two parts are taken from Eq. (58) and Eq. (72), and the latter can only be found numerically as we showed in Fig. 11. Furthermore, we estimate the sensitivity towards the three uncertainty parameters through a χ 2 -analysis for the traditional rate measuring method (RMM) to look for unexpected disappearance or appearance signals caused by SD and/or PD. Moreover, since neutrino decoherence is more enhanced at low energy, from Fig. 12 we find that for current experiments, reactor neutrinos should have the best sensitivity. Hence, in the section below, as a benchmark experiment to estimate how far we are from detecting neutrino decoherence through the damping term, we choose the RENO experiment. This is because it has less uncertainty compared to the Double Chooz experiment and a less complicated structure for the distribution of reactors and detectors compared with the Daya Bay experiment, which would affect the damping signature through σ L in a non-trivial way, see Appendix D. In addition, different dependence on L 0 and E 0 in the damping signatures in neutrino oscillation is discussed in e.g. [39,62]. The RMM, on the other hand, is a lot less sensitive to the phase shift terms, since a shift in the phase would barely vary the shape of the neutrino spectrum, as we will further illustrate in Fig. 17. Hence, for the phase shifting signatures, we introduce a method to measure the oscillation phase directly, which could be done by moving the detector around an expected oscillation minimum, namely the phase measuring method (PMM). Such method not only has the purpose of measuring asymmetries in the weighting functions, it is also a cleaner way to measure neutrino oscillation signatures as we will further discuss in Sec. 4.2. As for the theory input, only the H E -induced PD will have a non-trivial contribution for Gaussian distributed weighting functions. We thus consider a two-Gaussian distributed Dp to introduce an asymmetry for the quantum uncertainties, and ignore H L -induced PD, since it is has only negligible contribution for groundbased experiments. At last, after introducing the PMM and evaluating theoretic inputs and parametrizations (an asymmetry parameter "a" for the quantum uncertainties and the width σ E for the energy uncertainty), we estimate the statistically and systematic uncertainties required for having a 90% CL sensitivity in the parameter space, considering a πDAR neutrino source.

Rate Measuring Method
In order to observe the damping term φ jk which, in general, results from any source of decoherence, we analyze the total count rate of neutrinos for ground-based neutrino experiments. This is done by fitting our theory for SD through σp, as well as PD through σ L and σ E to oscillation data, defining the χ 2 -function: Here R i could either be the observed rate or the ratio of the detected rate of the near and far detectors; b j are the pull parameters, which include the oscillation parameters and the experimental uncertainties of the rate, and U i represents the statistical uncertainty of the bin.
With the purpose of determining what experiments are more sensitive for each decoherence parameter, we plot Fig. 12, in which the red, blue and yellow lines represent contours of |P (σ n = 0)−P (σ n = 0)| = 10 −6 for the solid lines and 10 −4 for the dashed lines, for n = L, E andp respectively, giving us a hint of what experiments to look at for a certain σ n . For simplicity, P represents the FTP on the third layer forν e →ν e in this section. We can see from the figure that for all decoherence effects, the influence would be larger at lower energies, because the oscillation structure is denser along thep and E axes, so the PWO effect is enhanced. Therefore, reactor neutrinos having the lowest energies for ground-based experiments would be the best candidate. Additionally, for vacuum oscillation, the decoherence effect by σ L is small and does not depend on L 0 , hence, it is suitable for experiments near the source where the statistics are high, whereas σp and σ E are more pronounced at larger distance. Moreover, a more realistic version of Fig. 12 for reactor neutrinos is plotted in Fig. 13 for the each σ n , where the contour lines represent the rate difference between decoherent and coherent fluxes, |Φ(σ n = 0) − Φ(σ n = 0)|, considering the energy spectrum of neutrinos for RENO and also the diffusion over distance by where L bm /N (E 0 ) is the average distance/spectrum of the near detector for the RENO experiment. We take the RENO experiment as a benchmark, and see how much more statistic is needed to have enough sensitivity for each of the three parameters. The sensitivity we get from fitting the decoherence parameters with current RENO far-to-near ratio data in [64] by Eq. (78) is given in three left plots of Fig. 14. Corresponding with the formalism in [64][65][66] , is the observed far-to-near ratio of IBD candidates in the i-th energy bin after background subtraction, taken from the supplemental material in [64]; and the theoretical input is where f, , τ, b F , b N are the pull terms for the systematic uncertainties, namely the uncorrelated reactorflux systematic uncertainty, the uncorrelated detection, the timing veto systematic uncertainty, and the background uncertainties for near and far detectors, respectively, given in [65]. Also, τ , the uncorrelated energy-scale systematic uncertainty, is inserted by scaling the energy, i.e. E 0 → (1 + τ )E 0 , and N MC i is the near-to-far ratio without oscillation traced back from the Monte Carlo data given in [64]; Furthermore, L near l /L far l is the distance to the l-th near/far detector, and sin 2 (2θ 13 ), ∆m 2 ee = cos 2 θ 12 ∆m 2 31 + sin 2 θ 12 ∆m 2 32 are the coherent oscillation parameters. Finally, we plot the three left plots in Fig. 14  Additionally, experiments (see [63] for a review) with their corresponding baseline and typical neutrino energies are labeled on the plot, for accelerator neutrinos (blue), decay-at-rest neutrinos (pink) and reactor neutrinos (green). χ 2 -function with the Python package "iminuit" over α = (sin 2 (2θ 13 ), ∆m 2 ee , f, , τ, b F , b N ) and β = ( α, σ n ) for ∆χ 2 (σ n ) = min The oscillation parameters are effectively marginalized by first minimizing over all the parameter, β, to get the best-fit values s 0 = 0.087, m 0 = 2.66 × 10 −3 eV 2 with errors σ s = 0.023, σ m = 0.12 × 10 −3 eV 2 for sin 2 (2θ 13 ) and ∆m 2 ee , respectively. Then we add two more pull terms to the χ 2 -function, The fitting results are shown in the three left plots of Fig. 14, from which we take the 90% CL limit as our benchmark points. They read σ bm E = 0.12 √ E 0 MeV, σ bm L = 548 m, and σ bm p = 1.6 MeV. In particular, σ bm p = 1.6 MeV is consistent with the analysis for the RENO experiment in [37,38]; it is obvious that σ bm L is still far from being a realistic value, which should be a few meters; nonetheless, σ bm E seems to be close to the energy resolution which is σ E / √ E 0 = 0.08 E 0 (MeV) + 0.3 [65]. Furthermore, in order to estimate how many times more statistics we need to reach the required sensitivity for a reasonable σ n , we assume that the statistical uncertainty will be increased with some value √ N i for each energy bin. On top of that, if uncertainties of the pull parameters are small enough and the signal count is much larger than the background count, we may consider only having the statistical uncertainties left in χ 2 . In this case Therefore, λ represents the enhancement of statistics needed to have a 90% CL signal for some value σ n . Moreover, for λ = 1, Eq. (83) means to sum over all the energy bins in Fig. 13 for some σ bm n . Hence, with the benchmark values we get from the three left plots of Fig. 14, and the flux difference from Fig. 13 inserted into Eq. (83), we arrive at the right plot in Fig. 14, where we see how much more statistics we need (λ in Eq. (83)) compared to RENO to gain a 90% CL sensitivity. This increment (λ) could be achieved by lowering the energy threshold for neutrino detection, increasing the reactor power or simply waiting for more data to be collected over time. Unsurprisingly, we are more sensitive to σ E and σp at larger distances, despite that the statistics drops by 1/L 2 0 , while σ L prefers a shorter propagation distance. Hence, the ranges of L 0 are chosen accordingly. Furthermore, since L 0 cannot not be smaller then σ L , we leave the triangular area on the upper left blank, and the edge of that area means that the detection is exactly beside the source.

Phase Measuring Method
Decoherence effects not only include damping signatures, but would also cause a shift to the coherent phase ψ jk on the measurement layer by β jk ( σ) defined in Eq. (47), when there is an asymmetry in the weighting function. In this section, we propose an illustration of a, in principle realistic, method which is extensively more sensitive to the phase shift terms compared to the damping terms in contrary to the RMM, by measuring where an oscillation extremum occurs on the third layer. Specifically, for some neutrino energy  Figure 14: Taking RENO as a benchmark experiment for reactor neutrinos, the left three plots present the constraints on different decoherence parameters from a fit to RENO data. The right plot takes the obtained 90% CL limit on σ n from the left plots. It shows the contour lines for evaluating how-many-times statistics compared to the current RENO data are needed to achieve a 90% CL sensitivity for some decoherence parameter and baseline, which is λ in Eq. (83). The while area in the middle figure is cut out since it would indicate that the detector is inside the reactor core.
E 0 , the goal is to search for the deviation in distance ( * L min -L osc ) cased by the phase shift. For the case of two neutrino mixing, we measure L min = L jk min in and see how it differs from L jk osc = 4πnE 0 /∆m 2 jk due to an asymmetry in the decoherence effect. Moreover, we can search for where an extremum occurs, and the signal is thus concentrated at a single L min , instead of having an entire distribution. Henceforward, an advantage of such method is that since we search for an extremum, the non-oscillating part of the event rate would only influence the signal trivially. Therefore, as long as other factors, such as the production rate, detection rate, background, etc. do not have an extremum within the range in (L 0 , E 0 ) of interest, it would barely contribute to the signal.
At the first oscillation minimum in a two-neutrino oscillation case, we look for The approximation above is done by taking β(L jk min , E 0 ; σ n ) β(L jk osc , E 0 ; σ n ) as we have checked that higher orders in the expansion of the left-hand side around the right-hand side can be neglected. In fact, due to the large difference between the atmospheric and solar mass splitting, when we search around L 13 osc , the total L min is L 13 min for the three neutrino mixing paradigm. Therefore, complications arising from the interference between different mass splittings, such as that from the damping terms, can also be negligible. Nonetheless, we still consider a full three neutrino-mixing scenario in our simulation below, and find L min numerically even for the fully coherent case. Additionally, in the same fashion, it is also possible to find a certain E min for some L 0 according to Eq. (84). However, it is usually not possible to look for effects by σ E through ∆E min or σ L through ∆L min , since the former is usually much larger than the latter. Hence, while the uncertainties give rise to a phase shift, it is likely to lower the sensitivity of the corresponding variable even more. Also, since σ L is too small for ground-based neutrino sources even in terms of ∆E min , we do not consider it in this section. For the purpose of this paper, we focus on the phase shift caused by the decoherehnce effect. Nonetheless, the PMM simply measures the (effective) neutrino oscillation phase exclusively. Hence, this method would also include measurements such as the neutrino mass splitting, the CP phase and the mass hierarchy or even the existence of an additional sterile neutrino. Therefore, in principle, a global analysis of all relevant experiments including decoherence effects would need to be performed. Fortunately, as we will soon show, while ∆L min from the errors in the mass splitting (which would also indicate the mass hierarchy) scales with E 0 , that from the asymmetry of the intrinsic quantum uncertainties would saturate to a constant value when E 0 is above ∼ 5 − 10 MeV. Hence, there will be a distinctive dependence on E 0 between these contributions in ∆L min . Moreover, we will also show that the phase shift term from quantum decoherence effect is insensitive to traditional measurements of the neutrino spectrum, while other fundamental oscillation parameters are determined by these measurements with increasing precision. Therefore, for simplicity and illustration propose, we fix the neutrino mass at values determined by global analyses [36], assume a three flavor oscillation with normal mass ordering and take δ CP = π in this section.
Operationally, the way to find L min is to scan over L 0 (e.g. by moving the detector) around where we expect to observe the first local minimum for some neutrino energy E 0 . In particular, we consider counting neutrinos within some position bin ∆L bin , i.e.
where N (E 0 ) is the number of neutrinos produced times the detection rate, which is independent of the traveling distance L 0 . Here, we have assumed that we do not lose or gain neutrinos during its propagation. Nonetheless, even if we do take such consideration into account, it would only affect the signal of L min trivially, as long as it does not create a bump or dip for a certain L 0 within our range of interest. Next, in order to find L min , we look at when the normalized (numerical) derivative of N i , is zero. We have F i plotted as the red dots in Fig. 15, and the vertical dashed lines represent the position bins, L i ± ∆L bin /2. HereN i = (N i+1 (E 0 ; σ n ) + N i (E 0 ; σ n ))/2 is the normalization factor which would eliminate the correlated uncertainties between the position bins, similar to the purpose of having near-far detectors. Note that the bin size is required to be L osc , such that a "local" minimum would be observed. Therefore, L min is where the blue line connecting all dots (F (L 0 )) intersect with the black F (L 0 ) = 0 line with uncertainty labeled as red horizontal bars in Fig. 15. As for the determination of such uncertainty, we adopt the following steps: 1. Propagate the uncertainty of the count numbers N i (∆ sys /∆ stat for the systematic/statistic uncertianties) to the uncertainties of F i by the relation in Eq. (87) as the blue error bars in the left plot of Fig. 15. In particular, only the uncorrelated uncertainties remain due to the normalization factorN .
2. Connect (or fit) the error bars of F i (i.e. F i ± ∆ 2 sys + ∆ 2 stat ) and draw an uncertainty band as we have demonstrated in the left plot of Fig. 15. 3. The uncertainty of L min is then the intersection between the uncertainty band and F (L 0 ) = 0 labeled as red error bars in the left plot of Fig. 15. This uncertainty is therefore determined by the relation between the signal (F i ) and its uncertainty, as we will discuss in the following. Note that it is possible for such intersection to be infinite when the error of the maximal |F i | surpasses its value. See the right plot in Fig. 15 for instance, where the sensitivity of ∆L min goes to infinity when the statistics is too low. Moreover, the uncertainties depend on the chosen bin size, demonstrated in the middle plot of Fig. 15, in fact, while ∆L bin L osc , increasing the position bin size would lower the uncertainty. The reason of this is two-fold: 1) the statistics for one bin would increase, reducing the statistical uncertainty, 2) F i being enhanced by the bin size reduces the uncertainty of L min . As a matter of fact, for an oscillator such as sin(L/L osc ) (the role of the FTP), its derivative made discrete by bins (the role of F i ) is The approximation is valid when (L i ± ∆L bin )/L osc ∼ 2nπ for some integer n, which is well justified since we only search around the oscillation minimum. Here, we can see that for a fixed bin size, the signal is smaller at higher energies, hence the uncertainties would be larger. For instance, for ∆L bin = O(10) m and ∆ sys = O(1)%, the energy range with finite sensitivity lies within a few MeV. Furthermore, although α 0 in the expansion N i = i α i L i 0 is canceled out for the signal, it still contains uncorrelated uncertainties. Thus, without α 0 α i , for some i = 0, there will be a significant increase in the uncertainty of L min . In fact, this is what we have for observing neutrinos around the maximum oscillation value or if we look for disappearing neutrinos. Therefore, rare event measurements of appearance channels at minimum oscillation phase would be the better option for our method.
By taking a 50 m (around the detector size of Hyper-K [67] and the DUNE far detector [68]) bin size, only reactor neutrinos and DAR neutrino are in the energy range which leads to an appripriate oscillation length. In particular, the monochromatic neutrinos from πDAR are most suitable for the PMM with just one measurement, due to the following reasons: • The monochromatic neutrinos are produced sharply around 30 MeV, which is suitable for a 50 m bin size as we have demonstrated in Fig. 15.
• It provides a detectable appearance channel by producing ν µ which could oscillate into ν e . On the other hand, reactors only produceν e , hence its appearance channels are not detectable sinceν µ will be below the Cerenkov threshold in the sub MeV range.
• Since we consider a fixed E 0 , the monochromatic feature automatically satisfies the condition without wasting any neutrinos spread out in the spectrum. Hence, statistics-wise, on top of the bright spallation source, it would be better than having µDAR neutrinos if we only consider measurements of L min at a single E 0 .
• The systematic uncertainty would also be strongly reduced for πDAR neutrinos. First of all, the timing structure of DAR experiments [69][70][71] would enable identification between πDAR neutrinos and µDAR neutrinos. In fact,ν µ →ν e from µDAR would suffer from an intrinsic uncertainty since theν µ andν e production are indistinguishable [69][70][71]. Secondly, the energy reconstruction would be highly accurate for monochromatic neutrinos. Discussions on this topic can be found in [72].
Similar to what we did for RMM, we consider a benchmark experiment, and ask how far we are to having enough sensitivity for some decoherence parameters. In particular, we take numbers from the existing JSNS experiment [69], i.e. 1.114 × 10 23 proton-on-targets for 1 MeV power within 3 year, from which 64% would contribute to a πDAR process, producing monochromatic ν µ which would oscillate into ν e and be detected by a 17-ton gadolinium loaded liquid scintillator. Hence, for Eq. (86), we obtain Here, we adopted the cross-section for quasielastic scattering of ν e on proton from [73] as 7.5 × 10 −41 cm −2 at 30 MeV, and assume that the detector is moved to the oscillation minimum (the actual JSNS detector is placed 24 m from the source). In Fig. 15 we adjust the equation above by moving L 0 around its first oscillation minimum, then increase it λ times. In addition, similar to other DAR channels [69,74], the systematic uncertainties should be dominated by intrinsic uncertainties, i.e. theν e produced by µDAR, which take up approximately 3% of total amount of neutrinos produced at 30 MeV. Furthermore, one could also identify whether a neutrino comes from πDAR from the timing structure, for instance, in the JSNS setup,ν e from µDAR takes up only < 10% of the early time bin which is dominated by ν e from πDAR [69]. Hence, we take various systematic uncertainties in the range of 0.1-2% in Fig. 15. From the middle and right plot of Fig. 15, we find the sensitivity for some ∆L min by first estimating the uncertainty of L osc (i.e. when ∆L min =0) for some systematic and statistical uncertainty (from ∆ sys and λ), then further identify what values of ∆L min would be rejected by such data at 90% CL. Theoretical estimates for decoherence effects which lead to such ∆L min will be shown in the following paragraph. The phase shift from decoherence effect for ground-based neutrinos would mainly come from the asymmetry of quantum uncertainties decided by the weighting function Dp and the classical (statistical) energy uncertainty with weighting function H E . In particular, we consider H E being dominated by the energy resolution (i.e. H E is Gaussian distributed) and Dp as a two-Gaussian distribution generically formalised as where the width is σp = (1 + rs)σ p , according to the definition in Appendix A. This formalism represents scenarios such as neutrino produced or detected with two types of interactions simultaneously, with different probabilities and widths (r, s) and have slightly different expectation values for E 0 (E 0 ± dE 0 , in particular). Moreover, with the phase structure given in Eq. (71), the decoherence term is simply the Fourier transformation of D(p; E 0 ) fromp to α jk = ∆m 2 jk L 0 /(2E 2 0 ), and the phase shift is Furthermore, the fact that we search around the first minimum (ψ = 2π) and β jk L osc min implies that at high energies, where This can be seen in Fig. 16, where lines having the same a merge to one constant value at higher energies which is independent of both E 0 and σp. Such property is not generic for all sources of decoherence effect, in fact, only αp ,jk from the phase structure η jk = iαp ,jkp , cancels out the energy dependence with L jk osc in Eq. (85) exactly. For instance, in Plot D of Fig. 17, ∆L 13 min increases with energy only because σ E does as well. In fact, if σ E is not energy dependent, it would approach zero at large E 0 . Moreover, from Eq. (92) we can see that when s = 1, i.e. the two bumps have the same width, σp would have no role in the phase shift. In addition, the variance of H E (∆ E , weighting function on the second layer with width σ E ) must be larger or equal to that of Dp (∆p). In fact, when they are equal to one another, the energy would be measured to a quantum level. Hence if one keeps on lowering ∆ E , ∆p would be forced to lower accordingly and the uncertainty ofx would increase in order to fulfil the uncertainty principle. In this case, if we consider σ E = 0.08 √ E 0 , and scale Dp by scaling σp and dE 0 simultaneously to fit the constrain ∆p = ∆ E , we find a change from Plot B to Plot A and C in Fig. 16. Plot B, on the other hand, assumes that σ E is large enough (in this case, σ E ≥ 0.17 √ E 0 ) such the the quantum uncertainties are un-squeezed. Furthermore, we can see that the dashed lines are more influenced by the constraint from ∆ E than the solid lines, since they either have a larger dE 0 or σp, both indicating a larger ∆p. From plot A, we can see that if σ E is small, it would squeeze Dp and lower the phase shift; on the other hand, if σ E is large, then the blue area covers all the lines and there will not be enough sensitivity. Therefore, while there is still a little space out of the sensitivity line, ∆E min it is also not a suitable approach to measure a Dp-induced phase shift.
Finally, we estimate the sensitivity for the benchmarks in Fig. 16 for the PMM in the right plot of Fig. 17. Furthermore, in the left plot, we demonstrate how the RMM is not as sensitive to the phase shift term compared to the damping term. The blue band is the range of Dp-induced state decoherence which is not constrained by the combined analysis of reactor experiments from [38], i.e. the upper edge of the band represents Wp as a Gaussian with width σp = 0.47 MeV. On the other hand, while having O(100) m of ∆L min for the PMM, the colored lines (with the same parameter as those in Fig. 16) do not vary the FTP to an extend that is close to the limit set by the combined analysis (not to mention for just one single experiment). Moreover, while the RMM significantly depends on how the neutrino spectrum would be without oscillation, the phase shift, which slightly shifts the FTP, does not change the shape of the spectrum as the damping term does, hence, it can be easily compensated by non-oscillation related models. On the contrary, for the PMM, the signal is amplified by the oscillation length and is nearly independent of non-oscillation related models. The main disadvantage is the lack of statistics since we aim at searching for appearing flavors at the oscillation minimum. Nonetheless, from the right plot in Fig. 17, we see that  Figure 17: The left plot shows how the phase shift term is not sensitive to the rate measuring method compared to the damping term: The upper edge of the blue band represents the transition probability for σp at its upper limit given in [38] for a Gaussian distributed Dp, and the lower edge is the fully coherent case. The colour code of the lines corresponds to the same decoherence parameters given in Fig. 16, which are within the parameter space in the right plot. The right plot shows values of λ in Eq. (89) required to achieve a 90% CL sensitivity for the decoherence parameter space by the colour bar. The red line labeled "J-PARC+DUNE" gives the required λ by assuming a J-PARC-like source combined with a DUNE-like detector; and similarly, the blue line labeled "ESS (proposed)" considers one year of data taking of the ESS source and the water Cherenkov detector proposed in [76]. Both cases are assumed to have the baseline L 0 at the first oscillation minimum at 30 MeV (31829 m).
with the increment mainly by the detector size, the statistics would be enough for a 90% CL sensitivity for a range of decoherence asymmetry parameters of the quantum and classical uncertainties. Specifically, compared to the 17 T detector mass and a cross section of 7.5 × 10 −41 cm 2 of JSNS, the DUNE detector would have an increased detector mass of 40 kT, and the liquid argon material of the detector also enhances the cross section to 2.5 × 10 −40 cm 2 at 30 MeV [75], hence λ 7.8 × 10 3 in this case (red lines). As for the ESS setup proposed in [76], while using a water Cherenkov detector implies a lower cross section (3 × 10 −42 cm 2 at 30 MeV [77]), the detector mass would be increased to 538 kT, and the spallation source is also brighter by having 2.7 × 10 23 POT per year.

Conclusion
Owing to the increasing precision of neutrino oscillation experiments, neutrino decoherence effects may become approachable in future experiments, opening a new window to probe new physics. In this work we introduce the "layer structure" (illustrated in Fig. 1), which includes the concept of an open quantum system and classical statics while having QFT as the fundamental theory. This structure is particularly useful for understanding mechanisms behind decoherence signatures in neutrino oscillation experiments. For instance, quantum uncertainties such as coordinate uncertainties around the vertices and the lifetime of particles entangled with the system are parameterised as σx and σp. These two uncertainty parameters are the width of weighting functions w.r.t. the coordinate and momentum variables in the Wigner phase space (on layer 1), respectively. On the other hand, classical uncertainties caused by a lack of knowledge would also contribute to decoherence signatures, such as the production profile of neutrinos, energy resolution and errors of the energy reconstruction model. The former cause dominates the uncertainty parameter σ L while the latter two are included in σ E . These two parameters are the width of weighting functions w.r.t. the coordinate and energy variables in the relativistic phase space (on layer 2), respectively. We have shown that decoherence effects from all these uncertainty parameters come from phase wash-out effects, which are determined by a phase structure and some distribution. For each uncertainty parameter, there is a certain phase structure and some localized distribution with width as the corresponding parameter, resulting in a phase wash-out effect suppressing and/or causing a phase shift in the oscillation signature. The phase structure also characterises dependence on the traveling distance (L 0 ) and energy spectrum (E 0 ) for each uncertainty parameter, hence, enables us to identify the mechanisms behind decoherence signatures by analysing these parameters in the neutrino detection profile and/or spectrum. The phase structures are given in Eq. (63) and Eq. (32) (Eq. (69) and Eq. (35)) for uncertainties on the Wigner phase space and the relativistic phase space, respectively, for the time dependent (independent) case. Furthermore, we have classified neutrino decoherence in terms of its mechanism as state decoherence and phase decoherence. The former represents the separation of superposition (mass) states, and is dominated by quantum uncertainties; while the latter indicate averaging effect due to the information loss, and is mainly decided by macroscopic classical uncertainties.
In particular, we calculate the case of Gaussian distributed weighting functions and estimate how much more statistics we need for certain σp, σ L and σ E , to be sensitive to them at 90% CL in Fig. 14, by taking the RENO experiment as a benchmark. We find from Fig. 12 that when the current far detector is located e.g. 14 (1.4457, 0.01) km away from the source † , we need approximately 26 (2, 82) times more statistics to reach sensitivity to meaningful values σp = 0.1 MeV (σ E = 0.08 √ E 0 √ MeV, σ L = 3 m). Furthermore, we propose a novel method, the phase measuring method, to measure the asymmetry of weighting functions by searching an oscillation minimum. Particularly, we estimate the sensitivity of this method for a two-Gaussian distributed, Dp-induced, quantum mechanical uncertainty as well as the statistical uncertainty from the energy resolution in Fig. 17. While the energy resolution ranges typically from 1 − 10%/ √ E 0 ( √ MeV) for neutrino detectors, the asymmetry parameter a, could be caused by the quantum effect of having a superposition of different processes. For instance, having simultaneously quasi-elastic scatterings and inelastic scatterings for neutrinos scattering on nucleons, or by nuclear effects such as the Fermi motion [78]. In fact, quantitative estimation of the asymmetry parameter would need further investigation. To sum up, while the four uncertainty parameters σ x , σ p , σ L and σ E in our structure can be determined by some theoretical mechanisms, such as the wave packet size of the external particles, the type of collisions, matter effect, exotic effects like space-time fluctuation, etc; it could also be potentially measured experimentally through rate or phase measuring methods. Our considerations presented here provide the theoretical background for such analyses and can be applied to any experiment. Experimental improvements are necessary, for instance via better energy resolution or larger event numbers, or by other detection techniques made possible by e.g. developments in coherent elastic neutrino-nucleus scattering.
Decoherence effect in neutrino oscillation also has many potential beyond ground-based experiments. For instance, atmospheric neutrinos might be a promising possibility to see decoherence effects, especially the one mediated by σ L . The production profile of σ L for atmospheric neutrinos would be of O(10) km and asymmetric. Aside from the production profile, matter effects may also be included. Moreover, since the atmospheric neutrinos have a broad energy spectrum and can be detected at different zenith angles, we should be able to do a wide range tomography on the (L 0 , E 0 ) space to analysis decoherence effects. Furthermore, there could be contributions to the uncertainties that are not directly measurable by other approaches, such as off-shell mediators contributing to σ x . Therefore, a better understanding of the measurable uncertainties by other approaches would increase the sensitivity of probing new physics through neutrino decoherence. In addition to aiming at searching for mechanisms causing neutrino decoherence, one can also make use of the decoherence effect to investigate other aspects, for instance, the search of sterile neutrinos or the measurement of CP violation, by designing and/or engineering these parameters, such as controlling σ L by the distribution of neutrino sources/detector, and σx, σp by manipulating squeezed states used in quantum optics.
where L is the central value of Γ(x), such that Γ even (x + L) is even. In fact, according to the layer structure presented in the main text, L would also be the next level PS variable in our structure corresponding to x. Hence, similar to property 2, |Φ| ≤ 1 and β is non-zero only when Γ(x; L) is symmetric w.r.t. L. In fact, this is why we call Φ the damping term and β the phase shift term in this paper. Moreover, the wider Γ(x) is relative to the wavelength for the phase structure η(x), the smaller will Φ become.
• Property 4: For two distributions f (x) and g(x) with width σ f and σ g , respectively, the width of (f * g)(L), σ f * g is larger than either σ f or σ g , where " * " represents the convolution of two distributions. Whenever two function are related with the form there is a convolution between these two functions. This usually occurs when there are multiple sources of uncertainties taken into consideration, for instance, the total uncertainties of the PS variables from both the initial state and the final state (Eq. (108)), the production site and the detection site (Eq. (115)), or the external process and the internal process (Eq. (129)). The width of an arbitrary localized function f (x) is defined here as where f (x) is the normalized function of f (x), and f (x) = f (x)/max{|f (x)|}, such that its global maximum is unitary. Also, 1/2 √ π is inserted such that width of a Gaussian distributed function would have the width at one standard deviation and the other distributions are then defined accordingly. On the other hand, the width of the product of two function, σ f g , will be smaller than the individual widths of the functions σ f and σ g , since for |f (x)| ≤ 1 and |g (x)| ≤ 1. Moreover, by the convolution theorem, we can see that comparing to the trivial case where g is a delta function, and we have f = FT −1 [FT (f )], the width of FT (f ) FT (g) would decrease when the width of g is no longer zero, and hence σ f * g would increase. For example, if f (x) and g(x) are Gaussian distributions, then σ 2 f * g = σ 2 f + σ 2 g .
• Property 5: Convolution of a complex function, h(x) = f (x) e ipx , and a real function, g(x), is (h * g) (y) = e ip yf (p − p)g(p ) ≡ e ipY (y) I 1 (y)I 2 (p), wheref = FT [f ],g = FT [g] and the width of I(y) ∈ R is the same as that of (f * g)(y).
By the convolution theorem, Then by doing an inverse Fourier transformation from p to y on Eq. (102), we arrive at Eq. (101). Furthermore, the width of the product of two functions,f (p − p)g(p ) is independent of the parallel  Table 1: Properties of the width of the function h in terms of the original real function(s) f (and g). Here NC/PC means that σ h is negatively/positively correlated to σ f (and σ g ).
Function (width notation), Type Width relation Gaussian case H(y) (σ H ), Convolution type I PC , σ H > {σ f , σ g } σ H = σ f * g I(p) (σ I ), Convolution type II NC σ I = 1/σ f g production/detection vertices and q, k, q , k are the momenta of the initial and final states of the production and detection sites. Therefore, x and p represent the traveling distance and the momentum of the neutrino decided by the external particles and the position of the vertices, i.e. the states entangled to the neutrino. The process to reach Eq. (9) includes a series of Fourier transformations and convolutions, which is illustrated in Fig. 18, from which we can clearly see how each of the uncertainties carried by each of the external states and the vertices affect the weighting functions with the help of Table 1 and Table 2. Below, we will show the derivation from Eq. (8) to Eq. (9), and also the relations in Fig. 18. We first include all the uncertainties following Eq. (8): For all h = {q, k, q , k }, since the external states are on the mass-shell, h 0 = E h (h) = h 2 − m 2 h . Furthermore, if the wavepackets are sharply peaked at the expectation value h , then by the saddle point approximation, we can write E h (h) Then, for the production site, by doing a change of variables: {q, k} → {p, k}, where p = k − q, the integration over k performs the convolution between the initial state and the final state including the plane wave amplitudes M P j (q, k), i.e.
where we write M P j (p, k) = M P j (k − p, k)| p 0 =Eq(k−p)−E k (k), k 0 =E k (k) for convenience. For f P i and f P j being Gaussian functions with width σ q and σ k respectively, and applying the saddle point approximation, the momentum uncertainties from the external states at the production site are where P = q − k . The other terms are F P,j (p) F Dj (p) F′ P (ξ P (p) − E j (p)) F′ P ((y 0 2 − x 0 2 )v qk )F ′ P (ξ P (p) − E j (p)) F′ P ((y 0 2 − x 0 2 )v qk ) g 1 (x 1 ) g 2 (x 2 ) F j (p; P) G(x; X)

FT Type
Product Type Convolution Type I Convolution Type II Figure 18: The final (width of the) weighting functions (G(x; X) and F j (p, P)) in terms of the wavepacket (size) of each of the external particles, and the spatial uncertainty (size) at the vertices. Referring to Table 1 and Table 2, this diagram is useful for finding how the widths are related, and how they contribute to the width of the weighting functions, which are σ x and σ p in the main text.
whereF P is the Fourier transformation of F P . After analogous calculations for the detection site, Eq. (112) finally takes the form of the layer moving operator, with weighting function F j (p; P)G(x, X): A 2,j = d 3 p d 4 x F j (p; P)G(x; X) A 1,j (x, p), where x = x 2 − x 1 , F j (p, P) = F P j (p)F Dj (p), and the first layer transition amplitude is A 1,j (x, p) = e ipx dp 0 ν e −ip 0 This appears as the collection of configuration of all the energetically allowed states for some p given by the external particles. Moreover, since the neutrino propagates a macroscopic distance, it can be approximated as traveling on the mass-shell. Hence, p 0 ν = E j (p) ≡ p 2 + m 2 j , then we can replace A 1,j (x, p) → e −itE j (p)+ixp and F j (p) → F j (p)F P (ξ(p) − E j (p))F D (ξ (p) − E j (p)), in Eq. (114). Therefore, either from the calculation we have shown, or more efficiently from Fig. 18, the width of F j , σ p is Next, to move A 1,j onto the next layer, we first integrate out the coordinate space d 4 x, for X = (T, L), and have A 2,j = d 3 p e −iT E j (p)+iLp F j (p; P)G(p).
Then as long as one of the functions F P , F D ,F P ,F D orG is sharply peaked, we can apply the saddle point approximation at P j , such that d dp F j (p; P)G(p) then the approximation gives E j (p) E j + v j (p − P j ), resulting in the final form of the second layer transition amplitude as A 2,j = e −iE j T +iP j LΦ j (L j , P j ), where L j = L − v j T . In fact, with such approximation, and taking G(x; X) as the function after the integration of x (a Fourier transformation to the momentum space) is: wherem j = m j v j σ 2 t /(σ 2 t v 2 j + σ 2 x ). Therefore, Eq. (114) can be written as e −iE j T +v j P j T dp 3 dx 3 e ixp G(x; L j )F j (p; P j ), where G has width σ 2 x = σ 2 t v 2 j + σ 2 x and is centred at L j = L − v j T ; F j has width σ p and centred at P j = ∆ P j , where ∆ = 1 + 4σ 2 x σ 2 p , such that the saddle point from Eq. (120) is at P j . In general, if all the input distributions are Gaussian distributed and with the saddle point approximation, x and p both linear dependent, F j and G will also be Gaussian distributed with some width σ p and σ x , respectively. In this case, we obtain the formalism in Eq. (15).

B.2 Neutrinos represented directly
Instead of representing the first layer phase space for the neutrinos by its entangled states, we represent it by the neutrinos directly. Therefore, we need to leave y = y 2 − y 2 and p ν non-integrated. Nonetheless, as we will see later, with this representation, we cannot have uncertainties representing the external particles and the internal vertex explicitly in the weighting function, but as an effective one with either just space-time or energy-momentum uncertainties. Therefore, we still apply the other representation in the main text for the sake of investigating neutrino decoherence in terms of these two sorts of uncertainties. Nonetheless, the representation presented in this subsection could also be used in our structure, resulting in the same effects. In particular, with this representation we can derive the total spatial uncertainty mentioned in Sec. 2.1, since the width of the function of y would indicate the total coordinate uncertainty, when it is the only non-integrated variable. Hence following Eq. (112), × d 4 x 1 g P (x 1 )F P (y 1 − x 1 )e iP P (y 1 −x 1 ) d 4 x 2 g D (x 2 )F D (y 2 − x 2 )e −iP D (y 2 −x 2 ) , where ∆(p ν ) is the neutrino propagator in the momentum space. Here we write the wavepackets at the production and detection site in coordinate space as: F P j (y 1 − x 1 ) e iP P (y 1 −x 1 ) d 3 pF P j (p)e −i(y 0 1 −x 0 1 )ξ P (p)+i(y 1 −x 1 )p , where P P and P D are the saddle point of F P j (p) and F Dj (p ), respectively. If there is negligible energy loss during the neutrino propagation, then P P = P D P . Next, after the integration over x 1 and x 2 , Eq. (126) becomes d 4 y d 4 y 1 I P (y 1 ; P P )I D (y 1 − y; P D )e −[x P (y 1 )−x D (y 1 −y)] e iP y∆ (y), where y = y 1 − y 2 , and∆ (y) is the Fourier transformation of the propagator, i.e. the two point function with distance y of the neutrino. Here I P and I D represent the total coordinate uncertainties for the production and detection site respectively, which are the convolutions between the coordinate uncertainties of the external states and the vertices, i.e.
I P (y 1 ; P )e iP [y 1 −x P (y 1 )] = d 4 x 1 g P (x 1 )F P j (y 1 − x 1 )F P (y 0 1 − x 0 1 )e iP (y 1 −x 1 ) (130) This can be related to Property 5 in Appendix A. Finally, the large bracket in Eq. (129) represents the total coordinate uncertainty, which turns out to be the width of the convolution function (I P * I D )(y), or ((g P * F P j F P ) * (g D * F Dj F D ))(y). Moreover, with the association and commutation property for convolution, we can rewrite the total coordinate width as (G * F tot P j * F tot Dj )(y), where G is in Eq. (115),F tot P j =F P j F P andF tot Dj =F Dj F D . X 3 X 3 − Λ X 3 + Λ Figure 19: The cumulative distribution function (CDF) of Y (X 2 ; X 3 = 50) as a Gaussian PDF with width = 2 AU (arbitrary units of X 2 ), centered at 50 AU: N (X 2 ) for the red line; a delta function centered at 50 AU for the blue line; and N (X 2 ) cos(X 2 )/N (X 2 ) cos(X 2 /2) for the yellow/purple line; Λ is the cutoff value in Eq. (134).

D Phase Decoherence for Discrete Neutrino Sources
In this section, we show the formalism of phase decoherence effect including a damping term and a phase shift term for neutrino detection coming from multiple sources. In other words, we formulate the case where the weighting function on layer 2 for the coordinate uncertainty is composed of multiple delta functions discretely scattered. We start with the simple case where there are only two point-like sources located at and thus x (1) indicating that the damping term φ ≤ 1 as excepted since a + b = 1. In the case where the detector is placed far from all the sources, x 1 , x 2 ∆x, then x eff = x 2 , and c = a + b cos(∆x). Similarly, if there are three point like neutrino sources, Φ jk = e −iα jk L 3 a e iα jk x 1 + b e iα jk x 2 + c e iα jk x 3 ≡ φ (2) e iα jk (x (2) eff −L 3 ) , the damping term (φ (2) ) and the phase term (x (2) eff ) are obtained by replacing eff − x 3 , a → φ (1) and b → c in Eq. (139), and so on for more point-like sources.