# Role of strong correlation and spin-orbit coupling in $\textrm{LuB}_{4}$: a
first principle study

Ismail Sk  Department of Physics, Bajkul Milani Mahavidyalaya, Purba
Medinipur,West Bengal, 721655, India  Department of Physics, Kazi Nazrul
University, Asansol, West Bengal, 713340, India  Joydeep Chatterjee
Department of Physics, Indian Institute of Technology Kharagpur, West Bengal,
721302, India.  A Taraphder  Department of Physics, Indian Institute of
Technology Kharagpur, West Bengal, 721302, India.  Nandan Pakhira 2

###### Abstract

The recent observation of magnetization plateaus in rare-earth metallic
tetraborides has drawn a lot of attention to this class of materials. In this
work, we investigate the electronic structure of one such canonical system
$\textrm{LuB}_{4}$, using first-principle density functional theory, together
with strong Coulomb correlation and spin-orbit coupling (SOC) effects. The
electronic band structures show that $\textrm{LuB}_{4}$ is a non-magnetic
correlated metal with completely filled $4f$ shell. The projected density of
states (DOS) shows a continuum at the Fermi level (FL), arising mainly from
hybridized Lu $d$ and B $p$ orbitals, along with some discrete peaks, well
separated from the continuum. These peaks arise mainly due to core-level Lu
$s$, $p$ and $4f$ atomic orbitals. Upon inclusion of SOC, the discrete peak
arising due to Lu $p$ is split into two peaks with $j=1/2$, $j=3/2$ while the
peak arising due to Lu $4f$ orbitals splits into two peaks with $j=5/2$ and
$j=7/2$. These peaks will give rise to multiplet structure in core level X-ray
photo-emission spectroscopy. Inclusion of strong correlation effects pushes
the Lu $4f$ peak away from the FL while the qualitative features remain
intact.

##  I Introduction

Compounds involving Boron like Boron Carbide (BC), $\textrm{MgB}_{2}$ article
, hexagonal Boron Nitride h-BN C7TC04300G , $\textrm{SmB}_{6}$ and various
metallic tetra-borides display exotic mechanical and electronic attributes:
very high tensile strength, high $\textrm{T}_{c}$ superconductivity, semi-
metallic behavior with topological properties, Kondo insulating states,
magnetization plateaus etc. PhysRevB.91.155151 ; PhysRevB.88.180405 ;
PhysRevX.3.011011 . These materials are under intense investigation due to
their fascinating properties like chemical stability, hardness, high melting
points, low work function and magnetic ground states with geometric
frustration. They have potential tehnological applications such as narrow-band
semiconductors CHEN2004L6 ; TAKEDA2004471 ; GAVILANO20001359 , thermoionic
emission devices. Compounds like $\textrm{CrB}_{4}$ and $\textrm{MnB}_{4}$
exhibit exceptional mechanical properties in the form of superior rigidity,
ideal tensile and shear strength 10.1063/1.3692777 .

The rare-earth metallic tetraborides with chemical formula $R\textrm{B}_{4}$
($R$ being the rare-earth atom) have drawn interest due to recent observations
of magnetization plateaus in $\textrm{TmB}_{4}$ PhysRevB.93.174408 ,
$\textrm{ErB}_{4}$ Michimura  and $\textrm{NdB}_{4}$, $\textrm{HoB}_{4}$ Brunt
. In these systems the rare-earth atoms are arranged in a geometrically
frustrated Shastry-Sutherland lattice and show wide range of magnetic ground
states SRIRAMSHASTRY19811069 . The $\textrm{CeB}_{4}$($\textrm{f}^{1}$) and
$\textrm{YbB}_{4}$($\textrm{f}^{13}$) do not show any order, while
$\textrm{PrB}_{4}$ is ferro-magnetically ordered article2 . On the other hand,
$\textrm{TmB}_{4}$ and $\textrm{ErB}_{4}$ order anti-ferromagnetically at low
temperatures ETOURNEAU1979531  in the presence of weak external magnetic
field. The role of spin-orbit coupling and strong Coulomb correlations on the
electronic properties of $\textrm{RB}_{4}$ are thought to be key to
understanding the novel and exotic properties observed in these systems.

The strong correlation effects in $\textrm{RB}_{4}$ arise due to the presence
of localized $4f$ orbitals in rare-earth ions. The trong correlation, present
in various $3d$, $4d$ transition metal compounds as well as actinides with
$5f$ orbitals, is central to the understanding of various novel and exotic
properties observed in these systems. The valence state(s) of rare-earth atoms
in $\textrm{RB}_{4}$ are found to be in di-, tri- and tetra-valent states in
various rare-earth tetraborides goryachev1975electronic . Recently, in
$\textrm{YbB}_{4}$ PMID:35667370  fractional valance state of Yb between
$\textrm{Yb}^{2+}$ and $\textrm{Yb}^{+3}$ has also been reported.

In an earlier work EPJB_IsmailSk ; dae_ismail , the role of strong spin-orbit
coupling effects on the electronic properties of various rare-earth
tetraborides in their non-magnetic ground states were investigated using first
principles electronic structure methods as implemented in the Quantum Espresso
(QE). One of the key features observed in various systems was the splitting of
the Boron-derived peak and rare-earth $p$ orbitals into $j=1/2$ and $j=3/2$
multiplets. The splitting $\Delta E$ was found to be proportional to $Z^{n}$,
where $Z$ is the atomic number of the rare-earth atoms and $n\approx 3.81$.
However, the strong Coulomb correlation effects were neglected, and more
importantly, the partial and total DOS had no contribution from $4f$ orbitals
due to the absence of contributions from $4f$ orbitals towards the pseudo-
potentials used in QE based calculations.

In the present work we study the interplay of strong spin-orbit coupling and
Coulomb correlation on the electronic structure of one such canonical system
$\textrm{LuB}_{4}$, using first principles electronic structure methods as
implemented in the Vienna ab initio simulation package (VASP) KRESSE199615 ;
PhysRevB.54.11169 . Most interestingly, the pseudo-potential used in VASP
involves contribution from $4f$ orbitals of the rare-earth atom (Lu) as well.
Also this system does not show any magnetic ordering due to closed shell
atomic structure of Lu ($\textrm{4f}^{14}$), giving rise to vanishing of
orbital and spin angular momentum. The organization of the rest of the paper
is as follows: in Section II, we describe the crystal structure of
$\textrm{LuB}_{4}$. In Section III we elaborate on the computational details
for electronic band structure calculations. In Section IV, we discuss our
results on the electronic band structure and density of states of
$\textrm{LuB}_{4}$. Finally, in Section V we conclude.

##  II Crystal structure

$\textrm{LuB}_{4}$ has a tetragonal symmetry and belongs to space group P4/mbm
article3 ; FISK19811189 . In Fig. 1 we summarize the crystal structure of
$\textrm{LuB}_{4}$ from different perspectives. Fig. 1(a) displays the full
tetragonal structure consisting of alternate layers of Lu and B ions stacked
along $c$-axis. Fig. 1(b)shows the top view of the crystal structure. There
are two distinct types of B atoms - (i) planar and (ii) octahedral. Boron
atoms form octahedra as well as 7-atom rings in the $a-b$ plane
https://doi.org/10.1002/adma.201604506 . Ring forming planar B atoms (shown in
blue) which are not part of octahedra also forms dimers and these dimers are
arranged in a regular pattern. It is important to mention that B atoms play a
crucial role in the electronic structure as they are in the $sp$-hybridized
state.

Figure 1: (Color online) Tetragonal crystal structure of $\textrm{LuB}_{4}$.
Panel (a) : Full structure consisting of different layers of rare-earth (Lu)
ions (red) and B (green) stacked along $c$-axis. Panel (b) : Top view of the B
(green and blue) sub lattice (in the $a-b$ plane) comprising of 7 atom ring
and a square formed by the Lu (red) atoms.

##  III Computational Details

First-principles electronic structure calculations were performed using
density functional theory (DFT) PhysRev.136.B864 ; PhysRev.140.A1133  as
implemented in VASP. All the calculations were performed with Perdew-Burke-
Ernzerhof (PBE) exchange-correlation (XC) functionals in the framework of
generalized gradient approximation (GGA) PhysRevLett.77.3865  and linearized
augmented plane wave (LAPW) PhysRevB.59.1758  basis. The experimental lattice
parameters and atomic positions of $\textrm{LuB}_{4}$ are taken from material
project website Jain2013 . The $\textrm{LuB}_{4}$ compound shows tetragonal
lattice structure consisting of 20 atoms in the unit cell. There are 4 Lu and
16 B atoms in each unit cell. For visualization and perodicity of the
structure, VESTA Momma:db5098  software were used. The constituent atoms and
their valence state are: $\textrm{Lu : [Xe]}4{f}^{14}5{d}^{1}6s^{2}$,
$\textrm{B}:2{s}^{2}2{p}^{1}$, respectively. The energy tolerance was set to
$10^{-7}$ eV for structural relaxations. We have chosen a cut off energy $520$
eV throughout the calculations. All the structures were then fully relaxed by
conjugate gradient (CG) algorithm when the difference in energy between
succeeding iteration on each atom was less than $10^{-2}$ eV/atom. The
standard Monkhorst- Pack PhysRevB.13.5188 (MP) $k$-grid was used for accurate
integration over the irreducible Brillouin zone. $k$-mesh of size $11\times
11\times 11$ was used for relaxation and GGA self-consistent calculations
whereas a denser $k$-meshe of size $13\times 13\times 13$ was used for the DOS
calculations. The chosen $k$-meshes were good enough for well-converged
calculations. In the GGA+U calculations, we have used $U=8.0$ eV and $J=1$ eV.
These are typical values for $4f$ systems as reported in many calculations
PhysRevB.77.035135 . The electronic band structure was plotted along the path
involving high symmetry points $\Gamma$=$(0,0,0)$, $X$ =
($\frac{\pi}{a},0,0)$, $M$=($\frac{\pi}{a}$,$\frac{\pi}{a}$,0),
$Z$=(0,0,$\frac{\pi}{c}$), $R$=($\frac{\pi}{a}$,0,$\frac{\pi}{c}$),
$A$=($\frac{\pi}{a}$,$\frac{\pi}{a}$, $\frac{\pi}{c}$). Calculated band
structures were plotted along the high symmetry directions $\Gamma-X-M-\Gamma$
and $Z-R-A-Z$.

##  IV Results and Discussions

In Table1 we compare the lattice constants obtained from structural relaxation
under GGA and GGA+U approximations to their observed experimental values. As
can be clearly observed the calculated values are very close to the
experimental values which validates the computed relaxed structures against
the observed lattice structure. The calculated Fermi energy is tabulated in
Table2. It is clear that the Fermi energy is not much affected by the
inclusion of either strong Coloumb correlation effect or spin-orbit coupling
effects in this system.

Table 1: Comparison of experimental and calculated (this work) lattice constants under various approximations | Experimental | GGA | GGA+U
---|---|---|---
Materials |  $a$ (Å) |  $c$ (Å) |  $a$ (Å) |  $c$ (Å) |  $a$ (Å) |  $c$ (Å)
$\textrm{LuB}_{4}$ | $7.02687$ | $3.96821$ | $7.03172$ | $3.97630$ | $6.99467$ | $3.94233$
Table 2: Calculated Fermi energy | GGA | GGA+SO | GGA+U | GGA+U+SO
---|---|---|---|---
$\textrm{E}_{f}(eV)$ | $5.7625$ | $5.7567$ | $5.7694$ | $5.7638$

###  IV.1 GGA and GGA+SO

From the relaxed structure we calculate the self-consistent field (SCF) and
subsequently the electronic band structure and density of states. Because of
the completely filled $\textrm{4f}^{14}$ configuration of the Lu atom,
$\textrm{LuB}_{4}$ is inherently non-magnetic. The band structure and DOS
plotted in Fig. 2 and Fig. 4 also confirm the non-magnetic metallic state of
$\textrm{LuB}_{4}$. The Fermi level, indicated by red dashed line in both the
plots, is set to zero. In Fig. 2(a) and 2(b) we show the band structure under
GGA and GGA+SO approximations, respectively. It is clearly observed from
Fig.2(a) that only four bands (indicated by numbers $1\cdots 4$) cross the
Fermi level. This is similar to the reported band structure in
$\textrm{TmB}_{4}$ pakhira2018electronic ; PRADHAN20175532 . Also as reported
earlier PhysRevB.95.205140  for $\textrm{TmB}_{4}$, the bands near the
symmetry point $M$ arises mainly due to hybridization between dimer B $2p$
orbital and Lu $4f$ orbital. In Fig. 2(b) we show the SOC effect on the
electronic band structure. A noticeable feature is the splitting of the
degenerate bands (labelled 5 and 6) at the $\Gamma$ point at an energy 1.5 eV.
The large splitting gap is due to large atomic spin-orbit coupling present in
Lu. Degeneracy-lifting is also present at various points along the $\Gamma-
X-M-\Gamma$ directions. Along the $\Gamma-Z$ direction the detailed features
of some of the bands (labelled 3 and 4) get modified in the presence of SOC.
Along the $Z-R-A-Z$ direction and especially at the $A$ point, there is a
large gap between the top and bottom bands. There are many flat bands arising
due to deep core level states far away from the Fermi level, not shown in the
figure. One additional band (orange, labelled 7) splits from band 3 and is
found to cross the FL along the $\Gamma$-Z path in the presence of spin-orbit
coupling. In Fig. 3(a) and (b) we have plotted the bands in a reduced energy
scale and only along the shortened $\Gamma-Z$ path. The additional band split
from band 3 in Fig. 3(a) is now clearly visible in Fig. 3(b).

Figure 2: (Color online)Panel : (a) Electronic band structure for
$\textrm{LuB}_{4}$ under GGA approximation. Panel : (b) Electronic band
structure under GGA+SO approximation. The Fermi level is set to zero. Figure
3: (Color online)Panel: (a) and (b) Electronic band structure for
$\textrm{LuB}_{4}$ under GGA and GGA+SO approximation plotted in reduced
energy window and along the shortened path $\Gamma$-$Z$ path. The band
(labeled 3) in the left panel splits into two bands (labeled 3 and 5) in the
right panel.

In Fig. 4 we plotted projected DOS (PDOS) as well as total DOS (TDOS) from
various orbitals at a given site under GGA and GGA+SO approximations,
respectively. In panel 4(a) and 4(b) are shown PDOS due to Lu and B,
respectively, in the absence of SOC. As can be observed in panel 4(a), a
discrete peak at -15.5 eV, well separated from the continuum DOS at Fermi
level, arises due to B-$s$ and B-$p$ orbitals. Seen in panel 4(b), discrete
large spectral peaks at -56.3 eV, -27.0 eV and -5.0 eV arise due to deep core
level Lu-$s$, Lu-$p$ and Lu-$4f$ orbitals, respectively. The DOS at FL, as
shown in the inset of panel 4(a), comes mainly from the hybridized B-$p$ and
Lu-$d$ orbitals, with minor contributions from B-$s$, Lu-$s$ and Lu-$4f$
orbitals.

In panel 4(c) and 4(d) we have shown the effect of SOC on PDOS due to Lu and
B, respectively. As depicted in panel 4(c), except minor redistribution, the
PDOS due to B largely remains unaffected. On the other hand, in panel 4(b),
the Lu-$p$ peak at -27.0 eV gets split into two peaks with $j=1/2$ and
$j=3/2$, respectively. This feature was also observed in our earlier study
EPJB_IsmailSk ; dae_ismail . Additionally, the discrete peak due to $4f$
orbital of Lu atom gets split into two peaks with $j=5/2$ and $j=7/2$,
respectively. The continuum density of states around the Fermi level gets
enhanced in the presence of SOC, shown in the inset of panel 4(d). The
enhancement is mainly due to Lu-$s$ and $d$ orbitals. However, the
contribution from Lu-$4f$ orbital decreases. The top panels show the TDOS in
each case. Due to the non-magnetic nature of $\textrm{LuB}_{4}$, TDOS for only
one spin species is shown.

Figure 4: (Color online) Panel (a) and (b) : PDOS of $\textrm{LuB}_{4}$ in the
absence of SOC. Contributions from different orbitals of B and Lu are shown.
Discrete peaks at -56.3 eV, -27 eV and -5 eV arises mainly due to Lu s, Lu p
and Lu f orbitals while the peak at -15.5 eV due to major B s with small
contribution from B p orbital. Panel (c) and (d) : PDOS in the presence of
spin-orbit coupling effect. The peak at -27 eV gets split into two peaks with
$j=1/2$ and $j=3/2$ and another peak due to Lu-$4f$ orbital get split into two
peaks with $j=5/2$ and $j=7/2$. Top left and right panel shows total dos from
all orbitals under both of these approximations.

###  IV.2 GGA+U and GGA+U+SO

The correlation effect in rare-earth tetra-borides is quite significant and
the electronic structure calculations should include these effects together
with the strong SOC. In Fig. 5(a) and 5(b) we have summarized the electronic
band structure of $\textrm{LuB}_{4}$ in the presence of U and U+SO,
respectively. As observed in Fig. 5(a), except an overall shift, the
qualitative features of the energy bands near the Fermi level remain the same
as in $U=0$. However, far away from the Fermi level, the qualitative features
of some of the bands (like band 5 and 6) are different from the non-
interacting case. It is important to mention that the effect of correlation,
treated within the DFT so far, is a static mean-field effect. However, treated
within the dynamical mean field theory (DMFT), it can drastically change the
electronic band structure through the time-dependent self-energy functional.
Inclusion of such effects is, for the time being, beyond the scope of this
study. There is significant degeneracy-lifting at the $\Gamma$ point between
band 5 and 6. Also, as in the earlier case, degeneracy lifting is visible at
various points along the $\Gamma-X-M-\Gamma$ direction. However, there is
degeneracy-lifting at an additional point ($Z$) at an energy $\sim 0.75$ eV. A
detailed analysis also shows Lu-$4f-2p$ hybridization near $Z$. Inclusion of
the SOC, as in $U=0$ case, gives rise to one additional Fermi level-crossing
band (orange colour, labelled 7) shown in Fig. 5(b).

As in the earlier case, in Fig. 6(a) and (b) we have plotted Fermi level-
crossing bands along the shortened $\Gamma$-Z path in a reduced energy window
under GGA+U and GGA+U+SO approximations, respectively. Due to the strong
correlation effect, one of the bands (labelled 3, Fig. 6(a)) now nearly
overlaps with another band (labelled 2). In the presence of strong SOC, as
shown in Fig. 6(b), an additional Fermi level-crossing band (labelled 5)
appears. However, at the Fermi level two bands (labelled 2 and 3) now overlap.
This overlap is due to interplay between strong correlation and strong SOC
effects.

In Fig. 7 we have summarized the projected as well as total density of states
for $\textrm{LuB}_{4}$ under GGA+U and GGA+U+SO approximations, respectively.
In panel 7(a) we have shown PDOS due to B orbitals. As in the earlier case
(GGA) a discrete peak at -15.5 eV, well separated from the continuum of
density of states at the Fermi level, arises mainly due to B-$s$ deep core
level states with minor contribution from B-$p$ orbital. The position of this
peak is unaffected by the correlation effects. In panel 7(b) we have shown
PDOS due to Lu orbitals. As in the earlier case three discrete peaks at -56
eV, -26.75 eV, and -8 eV arise mainly due to Lu-$s$, Lu-$p$, and Lu-$4f$ core
levels states. The position of the peaks due to Lu-$s$ and Lu-$p$ slightly
gets shifted. However, the position of the peak due to Lu-$4f$ orbital gets
significantly modified. This is largely expected as the correlation effect is
largest for localized $4f$ orbitals and the correlation effects on $p$, and
$s$ orbitals are only through hybridization with $4f$ orbitals. In panel 7(c)
and 7(d) we have shown the PDOS due to B and Lu orbitals under GGA+U+SO
approximation, respectively. As in the earlier case, inclusion of SOC splits
some of the discrete peaks while leaving other peaks largely unaffected. The
peak at -26.75 eV arising due to the Lu-$p$ orbital gets split into two peaks
with $j=1/2$ and $j=3/2$, and another peak at -8 eV arising due to Lu-$4f$
orbital gets split into two peaks with $j=5/2$ and $j=7/2$, respectively.
Discrete peaks at -56 eV arising due to deep core level Lu-$s$ orbitals and at
-15.5 eV arising due to B-$p$ and B-$s$ orbital remains unaffected.The DOS at
the Fermi level gets enhanced due to Lu $4f$, $p$ and $d$ orbitals while the
contribution from B orbitals largely remains the same.

Figure 5: (Color online)Panel : (a) and (b) Electronic band structure for
$\textrm{LuB}_{4}$ under GGA+U and GGA+U+SO approximations, respectively.
Figure 6: (Color online)Panel : (a) and (b) Electronic band structures for
$\textrm{LuB}_{4}$ under GGA+U and GGA+U+SO approximations plotted in reduced
energy window and along shortened $\Gamma$-Z path. Figure 7: (Color online)
Panel (a) and (b) : PDOS of $\textrm{LuB}_{4}$ under GGA+U approximation.
Contributions from different orbitals of B and Lu are shown. Discrete peaks at
-56 eV, -26.75 eV and -8 eV arises mainly due to Lu-$s$, Lu-$p$ and Lu-$4f$
orbitals while the peak at -15.5 eV arises due to major contributions from
B-$s$ with minor contribution from B-$p$ orbital. Panel (c) and (d): PDOS of
$\textrm{LuB}_{4}$ with GGA+U+SO consideration. The peak at -26.75 eV gets
split into two peaks with $j=1/2$ and $j=5/2$ and another peak for Lu-$4f$
orbital get split into two peaks with $j=5/2$ and $j=7/2$, respectively. The
density of states at the Fermi level arises mainly due to hybridized Lu-$d$
and Lu-$4f$ orbitals with minor contributions from other orbitals. Top left
and right panel shows TDOS from all the orbitals. Figure 8: (Color online)
TDOS of $\textrm{LuB}_{4}$ under GGA, GGA+U and GGA+U+SO approximations. Panel
(a) : Comparison between GGA and GGA+SO approximations. The discrete peaks due
to Lu-$p$ and Lu-$4f$ orbitals gets split into peaks characterized by $j=1/2$,
$j=3/2$ and $j=5/2$, $j=7/2$, respectively. Panel (b) : Comparison between GGA
and GGA+U approximations. The discrete peak due to Lu-$4f$ orbital gets
shifted away from the Fermi level by nearly 3 eV while the shifts of other
discrete peaks are less significant. Panel (c) : Comparison between GGA+U and
GGA+U+SO approximations. Overall features are similar to case (a) but the
position of various peaks are different due to strong correlation effects.

In Fig.8 we have summarized total density of states under various
approximations. In panel 8(a) we have shown comparison between GGA and GGA+SO
approximations. As mentioned earlier the discrete peak at -27 eV arising due
to Lu-$p$ and $s$ orbitals gets split into two peaks with $j=1/2$ and $j=3/2$
while the peak at -5 eV arising due to Lu-$4f$ orbital gets split into two
peaks with $j=5/2$ and $j=7/2$. The density of states at and around the Fermi
level enhances. However total integrated spectral weight under GGA from both
the spin channels and under GGA+SO approximations remain the same. In panel
7(b) we have shown comparison between GGA and GGA+U. The density of states at
the Fermi level does not change but the discrete peak at -5 eV shifts to -8 eV
due to correlation effects. The shape of the peak remains same. Further, as
shown in panel 8(c), upon turning on the SOC effect, the discrete peaks at -8
eV and at -27 eV gets split into peaks characterised by $j=5/2$, $j=7/2$ and
$j=1/2$, $j=3/2$ respectively.

The splitting of various peaks arising due to strong SOC effect will be
directly visible as detailed multiplet structure in core-level X-ray photo-
emission spectroscopy. However, various non-local screening effects not
included under GGA approximations can broaden various peaks. These non-local
screening effects can be included under various non-local approximations
beyond GGA like GW scheme and beyond the scope of the present study.

##  V Conclusions

To summarise, we have studied the electronic structure of $\textrm{LuB}_{4}$
under strong correlation and spin-orbit coupling. The electronic band
structure under GGA shows 4 Fermi level-crossing bands as reported earlier for
$\textrm{TmB}_{4}$. Inclusion of strong SOC under GGA+SO approximation shows
strong degeneracy-lifting at $\Gamma$ point between two bands (labelled 5 and
6). Also, the degeneracy-lifting is visible at various other points along the
$\Gamma-X-M-\Gamma$ direction. Most interestingly, in the presence of SOC, an
additional band appears to cross the Fermi level along the $\Gamma$-$Z$
direction. The projected density of states shows discrete peaks at -56 eV, -27
eV and -5 eV arising due to Lu-$s$, Lu-$p$ and Lu-$4f$ orbitals, respectively.
The density of states at the Fermi level is due to hybridized Lu-$p$, $s$ and
$4f$ orbitals. Inclusion of strong correlation effects through local $U$
causes the shifting of -5 eV peak further away from the Fermi level while the
other features remains nearly same. Upon inclusion of strong SOC effects the
peaks at -8 eV gets split into two peaks characterized by $j=5/2$, $j=7/2$
while the peak at -27 eV gets split into two peaks characterized by $j=1/2$,
$j=3/2$. Also, the density of states at and around the Fermi level shows
enhancement.

The present study will be extended to other rare-earth tetra-boride systems
with magnetically ordered ground states to understand the interplay of strong
correlation, strong SOC with magnetism in these fascinating systems. Also, the
dynamic correlation effects, not included under DFT approximations, are
extremely important and will be included in a future study using dynamical
mean field theory (DMFT) as implemented under LDA+DMFT approximation.

## Acknowledgements

This work is partially supported by WB-DSTBT research grant no. STBT-
11012(26)/31/2019-ST SEC. One of us (NP) would like to acknowledge hospitality
of IIT, Kharagpore where part of the manuscript was written.

## References

  * [1] J. Nagamatsu, N. Nakagawa, T. Muranaka, Y. Zenitani, and J. Akimitsu.  Nature, 410:63–4, (2001).
  * [2] K. Zhang, Y. Feng, F. Wang, Z. Yang, and J. Wang.  J. Mater. Chem. C, 5:11992, (2017).
  * [3] Tay-Rong Chang, T. Das, Peng-Jen Chen, and M. Neupane et al.  Phys. Rev. B, 91:155151, (2015).
  * [4] S. Wolgast, Çağlı yan Kurdak, and Kai . Sun et al.  Phys. Rev. B, 88:180405, (2013).
  * [5] X. Zhang, N. P. Butch, P. Syers, S. Ziemak, Richard L. Greene, and J. Paglione.  Phys. Rev. X, 3:011011, (2013).
  * [6] Chun-Hua Chen, Takashi Aizawa, Nobuo Iyi, Akira Sato, and Shigeki Otani.  Journal of Alloys and Compounds, 366:L6, (2004).
  * [7] M. Takeda, T. Fukuda, F. Domingo, and T. Miura.  Journal of Solid State Chemistry, 177:471, (2004).
  * [8] J.L Gavilano, B Ambrosini, H.R Ott, D.P Young, and Z Fisk.  Physica B: Condensed Matter, 284:1359, (2000).
  * [9] Huiyang Gou, Zhiping Li, Hui Niu, Faming Gao, Jingwu Zhang, Rodney C. Ewing, and Jie Lian.  Applied Physics Letters, 100:111907, (2012).
  * [10] Sai Swaroop Sunku, Tai Kong, Toshimitsu Ito, Paul C. Canfield, B. Sriram Shastry, Pinaki Sengupta, and Christos Panagopoulos.  Phys. Rev. B, 93:174408, (2016).
  * [11] Shinji Michimura, A. Shigekawa, F. Iga, M. Sera, Toshiro Takabatake, K. Ohoyama, and Yutaka Okabe.  Physica B: Condensed Matter, 378:596, (2006).
  * [12] D. Brunt, G. Balakrishnan, Daniel Mayoh, Martin Lees, D.I. Gorbunov, Navid Qureshi, and O. Petrenko.  Scientific Reports, 8, (2018).
  * [13] B. Sriram Shastry and Bill Sutherland.  Physica B+C, 108:1069, (1981).
  * [14] Hongchul Choi, Amel Laref, Ji Shim, S. Kwon, and B. Min.  Journal of Applied Physics, 105:07E107, (2009).
  * [15] J Etourneau, J.P Mercurio, A Berrada, P Hagenmuller, R Georges, R Bourezg, and J.C Gianduzzo.  Journal of the Less Common Metals, 67:531, (1979).
  * [16] Yu M Goryachev, BA Kovenskaya, and EM Dudnik.  Zh. Strukt. Khim.;(USSR), 16, 1975.
  * [17] F. Frontini, Blair W Lebert, KK Cho, MS Song, BK Cho, Christopher J Pollock, and Young-June Kim.  Journal of physics. Condensed matter : an Institute of Physics journal, 34, (2022).
  * [18] I. Sk and N. Pakhira.  The European Physical Journal B: Condensed Matter and Complex Systems, 96:1, (2023).
  * [19] I. Sk and N. Pakhira.  AIP Conference Proceedings, 2995:020099, 2024.
  * [20] G. Kresse and J. Furthmüller.  Computational Materials Science, 6:15, (1996).
  * [21] G. Kresse and J. Furthmüller.  Phys. Rev. B, 54:11169, (1996).
  * [22] J. Olsen, A. Waśkowska, L. Gerward, G. Vaitheeswaran, and K. Venkatakrishnan et al.  High Pressure Research, 31:3, (2011).
  * [23] Z. Fisk, M.B. Maple, D.C. Johnston, and L.D. Woolf.  Solid State Communications, 39:1189, (1981).
  * [24] Georgiy Akopov, Michael T. Yeung, and Richard B. Kaner.  Advanced Materials, 29:1604506, (2017).
  * [25] P. Hohenberg and W. Kohn.  Phys. Rev., 136:B864, (1964).
  * [26] W. Kohn and L. J. Sham.  Phys. Rev., 140:A1133, (1965).
  * [27] John P. Perdew, Kieron Burke, and Matthias Ernzerhof.  Phys. Rev. Lett., 77:3865, (1996).
  * [28] G. Kresse and D. Joubert.  Phys. Rev. B, 59:1758, (1999).
  * [29] A. Jain, Shyue Ping Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, and S. Cholia et al.  APL Materials, 1:011002, (2013).
  * [30] Koichi Momma and Fujio Izumi.  Journal of Applied Crystallography, 44:1272, (2011).
  * [31] Hendrik J. Monkhorst and James D. Pack.  Phys. Rev. B, 13:5188, (1976).
  * [32] Z. P. Yin and W. E. Pickett.  Phys. Rev. B, 77:035135, (2008).
  * [33] Nandan Pakhira, Jyoti Krishna, S. Nandy, T. Maitra, and A Taraphder.  Electronic structure of metallic tetra-boride $\textrm{TmB}_{\textrm{4}}$: An lda+dmft study, (2018).
  * [34] Subhasree Pradhan and A. Taraphder.  Materials Today: Proceedings, 4:5532, 2017.
  * [35] John Shin, Zack Schlesinger, and B. Sriram Shastry.  Phys. Rev. B, 95:205140, (2017).


$F_{i}(E_{0};\vec{\sigma}_{n})=\frac{1}{\bar{N}_{i}}\frac{N_{i+1}(E_{0};\vec{\sigma}_{n})-N_{i}(E_{0};\vec{\sigma}_{n})}{L_{i+1}-L_{i}},$
(87)

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}\pm\Delta L_{\rm bin}/2$. Here
$\bar{N}_{i}=(N_{i+1}(E_{0};\vec{\sigma}_{n})+N_{i}(E_{0};\vec{\sigma}_{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 $\ll L_{\rm osc}$, such
that a “local” minimum would be observed. Therefore, $L_{\rm 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. 1.

Propagate the uncertainty of the count numbers $N_{i}$ ($\Delta_{\rm
sys}/\Delta_{\rm 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 factor $\bar{N}$.

  2. 2.

Connect (or fit) the error bars of $F_{i}$ (i.e. $F_{i}\pm\sqrt{\Delta_{\rm
sys}^{2}+\Delta_{\rm stat}^{2}}$) and draw an uncertainty band as we have
demonstrated in the left plot of Fig. 15.

  3. 3.

The uncertainty of $L_{\rm 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 $\Delta L_{\rm 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 $\Delta L_{\rm bin}\ll L_{\rm 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_{\rm min}$. As a matter of fact, for an oscillator such as
$\sin(L/L_{\rm osc})$ (the role of the FTP), its derivative made discrete by
bins (the role of $F_{i}$) is

$\tilde{F}_{i}=\frac{1}{L_{\rm osc}}\int^{L_{i}+\Delta L_{\rm
bin}/2}_{L_{i}-\Delta L_{\rm bin}/2}dL\cos\left(\frac{L}{L_{\rm
osc}}\right)\simeq\frac{\Delta L_{\rm bin}}{L_{\rm osc}}.$ (88)

The approximation is valid when $(L_{i}\pm\Delta L_{\rm bin})/L_{\rm osc}\sim
2n\pi$ 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 $\Delta L_{\rm bin}=\mathcal{O}(10)$ m and
$\Delta_{\rm sys}=\mathcal{O}(1)\%$, the energy range with finite sensitivity
lies within a few MeV. Furthermore, although $\alpha_{0}$ in the expansion
$N_{i}=\sum_{i}\alpha_{i}L_{0}^{i}$ is canceled out for the signal, it still
contains uncorrelated uncertainties. Thus, without $\alpha_{0}\ll\alpha_{i}$,
for some $i\neq 0$, there will be a significant increase in the uncertainty of
$L_{\rm 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.

Figure 15:  The left plot shows how we determine $L_{\rm min}$ for neutrinos
with 30 MeV energy and its uncertainty (red error bars) from Eq. (87) (blue
error bars) considering two different systematic uncertainties labeled in the
plot. Here we consider the JSNS setup by taking Eq. (89) for $N(E_{0})$ in Eq.
(86), and each bin is separated by the black dashed lines, hence the
parameters which would influence the uncertainty of $L_{\rm min}$ are the
systematic uncertainty ($\Delta_{\rm sys}$), the statistical uncertainties
($\lambda$: increment w.r.t. the JSNS set up) and the bin size ($\Delta L_{\rm
bin}$). The middle and right plot show how these parameters would influence
the sensitivity for $\Delta L_{\rm min}$ through simulating the uncertainties
of $L_{\rm min}$. In particular, the middle plot corresponds to Eq. (88), and
the right plot shows when systematic/statistic uncertainties dominate over one
another. Specific steps to determine $L_{\rm min}$ and more discussions on the
plots are given in the text.

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 $\pi$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 $\nu_{\mu}$ which
could oscillate into $\nu_{e}$. On the other hand, reactors only produce
$\bar{\nu}_{e}$, hence its appearance channels are not detectable since
$\bar{\nu}_{\mu}$ 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 $\mu$DAR neutrinos if we only consider
measurements of $L_{\rm min}$ at a single $E_{0}$.

  * •

The systematic uncertainty would also be strongly reduced for $\pi$DAR
neutrinos. First of all, the timing structure of DAR experiments [69, 70, 71]
would enable identification between $\pi$DAR neutrinos and $\mu$DAR neutrinos.
In fact, $\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{e}$ from $\mu$DAR would suffer
from an intrinsic uncertainty since the $\bar{\nu}_{\mu}$ and $\bar{\nu}_{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\times 10^{23}$ proton-on-targets for 1 MeV power within 3 year, from
which $64\%$ would contribute to a $\pi$DAR process, producing monochromatic
$\nu_{\mu}$ which would oscillate into $\nu_{e}$ and be detected by a 17-ton
gadolinium loaded liquid scintillator. Hence, for Eq. (86), we obtain

$\frac{N(E_{0})}{4\pi L_{0}^{2}}\simeq
0.43\frac{1}{\rm{m}^{2}}\left(\frac{P}{1\,\rm MeV}\right)\left(\frac{T}{3\,\rm
yr}\right)\left(\frac{\rm M_{D}}{17\,\rm ton}\right)\left(\frac{31829\,\rm
m}{L_{0}}\right)^{2}.$ (89)

Here, we adopted the cross-section for quasielastic scattering of $\nu_{e}$ on
proton from [73] as $7.5\times 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
$\lambda$ times. In addition, similar to other DAR channels [69, 74], the
systematic uncertainties should be dominated by intrinsic uncertainties, i.e.
the $\bar{\nu}_{e}$ produced by $\mu$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 $\pi$DAR from the timing structure, for
instance, in the JSNS setup, $\bar{\nu}_{e}$ from $\mu$DAR takes up only
$<10\%$ of the early time bin which is dominated by $\nu_{e}$ from $\pi$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 $\Delta L_{\rm min}$ by first estimating the uncertainty of $L_{\rm
osc}$ (i.e. when $\Delta L_{\rm min}$=0) for some systematic and statistical
uncertainty (from $\Delta_{\rm sys}$ and $\lambda$), then further identify
what values of $\Delta L_{\rm min}$ would be rejected by such data at 90% CL.
Theoretical estimates for decoherence effects which lead to such $\Delta
L_{\rm min}$ will be shown in the following paragraph.

Figure 16:  Plot A (B, C) shows the variation of $E_{\rm min}^{13}$ ($L_{\rm
min}^{13}$) for a number of quantum uncertainty parameters for Eq. (90). The
labels for each line (as well as its colour and style) on all three plots (and
also in Fig. 17) correspond to the table on Plot A, while Plot A and C
consider a constrain by $\sigma_{E}=0.08\sqrt{E_{0}}$ MeV (blue area in Plot
A) in addition. Plot D shows the variation of $L_{\rm min}^{13}$ caused by
$\sigma_{E}$ from a Gaussian distributed energy resolution.

The phase shift from decoherence effect for ground-based neutrinos would
mainly come from the asymmetry of quantum uncertainties decided by the
weighting function $D_{\bar{p}}$ 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 $D_{\bar{p}}$ as a two-Gaussian distribution generically
formalised as

$D_{\bar{p}}(\bar{p}\equiv|\bar{\mathbb{p}}|;E_{0})=\frac{1}{2\sqrt{\pi}\sigma_{\bar{p}}^{\prime}(1+rs)}\left\\{\exp\left(\frac{-(\bar{p}-E_{0}+dE_{0})^{2}}{4\sigma_{\bar{p}}^{\prime
2}}\right)+r\exp\left(\frac{-(\bar{p}-E_{0}-dE_{0})^{2}}{4(s\sigma_{\bar{p}}^{\prime})^{2}}\right)\right\\},$
(90)

where the width is $\sigma_{\bar{p}}=(1+rs)\sigma_{\bar{p}}^{\prime}$,
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}\pm dE_{0}$, in
particular). Moreover, with the phase structure given in Eq. (71), the
decoherence term is simply the Fourier transformation of $D(\bar{p};E_{0})$
from $\bar{p}$ to $\alpha_{jk}=\Delta m_{jk}^{2}L_{0}/(2E_{0}^{2})$, and the
phase shift is

$\beta_{jk}(E_{0},\sigma_{\bar{p}})=\tan^{-1}\left(\frac{1-r\,s\,e^{-(s^{2}-1)\alpha_{\bar{p},jk}^{2}\sigma_{\bar{p}}^{\prime
2}}}{1+r\,s\,e^{-(s^{2}-1)\alpha_{\bar{p},jk}^{2}\sigma_{\bar{p}}^{\prime
2}}}\tan(\alpha_{\bar{p},jk}\,dE_{0})\right)\xrightarrow{\alpha_{\bar{p},jk}\ll
1}\,\frac{1-rs}{1+rs}\frac{\Delta m_{jk}^{2}L_{0}}{2E_{0}^{2}}\,dE_{0}.$ (91)

Furthermore, the fact that we search around the first minimum ($\psi=2\pi$)
and $\beta_{jk}\ll L_{\rm min}^{\rm osc}$ implies that $\Delta
m_{jk}^{2}L_{0}/(2E_{0})\simeq 2\pi$, hence

$\Delta L_{\rm min}^{jk}\simeq\frac{2\pi}{2.53\,\Delta m_{jk}^{2}}\,a,$ (92)

at high energies, where

$a=\frac{1-rs}{1-rs}\,dE_{0}.$ (93)

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
$\sigma_{\bar{p}}$. Such property is not generic for all sources of
decoherence effect, in fact, only $\alpha_{\bar{p},jk}$ from the phase
structure $\eta_{jk}=i\alpha_{\bar{p},jk}\bar{p}$, cancels out the energy
dependence with $L_{\rm osc}^{jk}$ in Eq. (85) exactly. For instance, in Plot
D of Fig. 17, $\Delta L_{\rm min}^{13}$ increases with energy only because
$\sigma_{E}$ does as well. In fact, if $\sigma_{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, $\sigma_{\bar{p}}$
would have no role in the phase shift. In addition, the variance of $H_{E}$
($\Delta_{E}$, weighting function on the second layer with width $\sigma_{E}$)
must be larger or equal to that of $D_{\bar{p}}$ ($\Delta_{\bar{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 $\Delta_{E}$, $\Delta_{\bar{p}}$
would be forced to lower accordingly and the uncertainty of $\bar{x}$ would
increase in order to fulfil the uncertainty principle. In this case, if we
consider $\sigma_{E}=0.08\sqrt{E_{0}}$, and scale $D_{\bar{p}}$ by scaling
$\sigma_{\bar{p}}$ and $dE_{0}$ simultaneously to fit the constrain
$\Delta_{\bar{p}}=\Delta_{E}$, we find a change from Plot B to Plot A and C in
Fig. 16. Plot B, on the other hand, assumes that $\sigma_{E}$ is large enough
(in this case, $\sigma_{E}\geq 0.17\sqrt{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 $\Delta_{E}$ than the solid lines,
since they either have a larger $dE_{0}$ or $\sigma_{\bar{p}}$, both
indicating a larger $\Delta_{\bar{p}}$. From plot A, we can see that if
$\sigma_{E}$ is small, it would squeeze $D_{\bar{p}}$ and lower the phase
shift; on the other hand, if $\sigma_{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, $\Delta E_{\rm min}$ it
is also not a suitable approach to measure a $D_{\bar{p}}$-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 $D_{\bar{p}}$-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
$W_{\bar{p}}$ as a Gaussian with width $\sigma_{\bar{p}}=0.47$ MeV. On the
other hand, while having $\mathcal{O}(100)$ m of $\Delta L_{\rm 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 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\times 10^{-41}$ cm2 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\times
10^{-40}$ cm2 at 30 MeV [75], hence $\lambda\simeq 7.8\times 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\times 10^{-42}$ cm2 at 30
MeV [77]), the detector mass would be increased to 538 kT, and the spallation
source is also brighter by having $2.7\times 10^{23}$ POT per year.

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 $\sigma_{\bar{p}}$ at its
upper limit given in [38] for a Gaussian distributed $D_{\bar{p}}$, 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
$\lambda$ 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 $\lambda$ 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).

##  5 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 $\sigma_{\bar{x}}$ and $\sigma_{\bar{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
$\sigma_{L}$ while the latter two are included in $\sigma_{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
$\sigma_{\bar{p}}$, $\sigma_{L}$ and $\sigma_{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†††The actual distance of RENO’s far detector is 1.4457
km., we need approximately 26 (2, 82) times more statistics to reach
sensitivity to meaningful values $\sigma_{\bar{p}}=0.1$ MeV
($\sigma_{E}=0.08\sqrt{E_{0}}\rm{\sqrt{MeV}}$, $\sigma_{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, $D_{\bar{p}}$-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\%/\sqrt{E_{0}}(\sqrt{\rm 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
$\sigma_{x}$, $\sigma_{p}$, $\sigma_{L}$ and $\sigma_{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 $\sigma_{L}$. The production profile of $\sigma_{L}$ for atmospheric
neutrinos would be of $\mathcal{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 $\sigma_{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 $\sigma_{L}$ by the distribution of neutrino
sources/detector, and $\sigma_{\bar{x}},\sigma_{\bar{p}}$ by manipulating
squeezed states used in quantum optics.

## Acknowledgments

We would like to thank Evgeny Akhmedov and Janina Hakenmüller for useful
discussions. TC acknowledges support by the IMPRS-PTFS.

## Appendix

##  Appendix A Fourier Transformation and Convolution Properties

We review some useful properties of Fourier transformation (FT) and
convolution in this section. These properties are particularly useful in our
structure, for the layer moving operators involve integrations of complex
functions, which can be parameterized as a plane wave term $\exp(i\eta(x,p))$
and a normalized-real probability density function (PDF) term. In particular,
we will show how the layer variables are connected with each other via FT
properties, how two sources of uncertainties are combined to an effective one
with convolutional properties, and demonstrate the phase washout effect.
Below, we will outline the properties in bold front followed by a more
detailed demonstration.

  * •

Property 1: The FT of an even function is real, while that of an odd function
is purely imaginary. For an even function $f(x)$,

$\begin{split}&\int^{\infty}_{-\infty}dx\,e^{-ipx}f(x)=\frac{1}{2}\int^{\infty}_{-\infty}dx\,e^{-ipx}\left\\{f(x)+f(-x)\right\\}\\\
&=\frac{1}{2}\int^{\infty}_{-\infty}dx\,e^{-ipx}f(x)-\frac{1}{2}\int^{\infty}_{-\infty}dx\,e^{ipx}f(x)=\int^{\infty}_{-\infty}dx\,f(x)\cos{(px)}\in\mathbb{R},\end{split}$
(94)

and similarly for odd functions being imaginary after Fourier transformation.

  * •

Property 2: For any probability density function (PDF), $W(x;L)$, the FT can
be written as

$P(L,p)=\int^{\infty}_{-\infty}dx\,e^{-ipx}W(x;L)=e^{-ipL}\tilde{W}(p)\equiv
e^{-i(pL-\beta(p))}|\tilde{W}(p)|,$ (95)

where the damping term $|\tilde{W}(p)|\leq\tilde{W}(0)=1$ and the phase shift
term $\beta(p)$ is non-zero and non-$\pi$ only when $W(x;L)$ is symmetric
w.r.t. $L$. If the shape of $W(x)$ is symmetric, $W(x+L)$ would be even for
$L=\int dx\,x\,W(x)$. Then, after shifting $x\rightarrow x+L$, we obtain Eq.
(95) with $\tilde{W}(p)\in\mathbb{R}$ and no phase shift, i.e. $\beta=0$. On
the other hand, if the shape of $W(x)$ is not symmetric, we can always write
it in terms of a an even function and an odd function, i.e. $W(x+L)=W_{\rm
even}(x+L)+W_{\rm odd}(x+L)$. Hence on top of the even part, which is treated
in the same manner as the symmetric case, the odd part would give rise to an
imaginary part in $\tilde{W}(p)$, or in terms of rotation coordinate, a phase
shift $\beta(p)\neq 0$ and $\pi$ relative to the oscillation phase, $pL$, on
the next layer. Moreover, $W(x)$ being a PDF indicates that $\int dx\,W(x)=1$,
and $W(x)\geq 0$, therefore, $|\tilde{W}(p)|\leq\int dx\,W(x)=\tilde{W}(0)=1$.
A list of examples is given in Table 10, showing how the asymmetry of the PDF
induces a non-zero and non-$\pi$ phase shift. Additionally, it is also clear
from the plots that, in most cases, the larger the width (labeled as
$\sigma_{n}$, for $n=\\{p,L,E\\}$) is for the PDF, the smaller will the width
of $\tilde{W}(p)$ be. Henceforward, since $|\tilde{W}(p)|\leq\tilde{W}(0)=1$,
the larger $\sigma_{n}$ is, the smaller $|\tilde{W}(p)|$ will be, for some
$p\neq 0$, and the more suppressed $P(L,p)$ will be. Another way to look at
this effect is that a wider width of the PDF indicates that there is a wider
range for $e^{ipx}$ to be averaged out upon the integration over $x$, namely,
the PWO effect.

  * •

Property 3: The PWO effect is the generic case of Eq. (95), for a complex
function $\Gamma(x;L)\equiv|\Gamma(x;L)|e^{i\eta(x)}$, and Eq. (95) is simply
when $\eta(x)$ is linear in $x$. The PWO effect is written as

$\frac{\int dx\,\Gamma(x;L)}{\int
dx\,|\Gamma(x;L)|}=e^{i(\eta(x)|_{x=L}-\beta)}\Phi(L),$ (96)

where $L$ is the central value of $\Gamma(x)$, such that $\Gamma_{\rm
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, $|\Phi|\leq 1$ and $\beta$
is non-zero only when $\Gamma(x;L)$ is symmetric w.r.t. $L$. In fact, this is
why we call $\Phi$ the damping term and $\beta$ the phase shift term in this
paper. Moreover, the wider $\Gamma(x)$ is relative to the wavelength for the
phase structure $\eta(x)$, the smaller will $\Phi$ become.

  * •

Property 4: For two distributions $f(x)$ and $g(x)$ with width $\sigma_{f}$
and $\sigma_{g}$, respectively, the width of $(f*g)(L)$, $\sigma_{f*g}$ is
larger than either $\sigma_{f}$ or $\sigma_{g}$, where “$*$” represents the
convolution of two distributions. Whenever two function are related with the
form

$\int dx\,f(x)g(x-L)\equiv(f*g)(L),$ (97)

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

$\sigma_{f}=\frac{1}{2\sqrt{\pi}}\int dx\,|f^{\prime}(x)|,$ (98)

where $f^{\prime}(x)$ is the normalized function of $f(x)$, and
$f^{\prime}(x)=f(x)/max\\{|f(x)|\\}$, such that its global maximum is unitary.
Also, $1/2\sqrt{\pi}$ 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, $\sigma_{fg}$, will be smaller than the
individual widths of the functions $\sigma_{f}$ and $\sigma_{g}$, since

$\displaystyle 4\pi\sigma_{fg}=\int dx\,|f^{\prime}(x)g^{\prime}(x)|\leq\int
dx\,|f^{\prime}(x)||g^{\prime}(x)|$ $\displaystyle\leq\int
dx\,|f^{\prime}(x)|=2\sqrt{\pi}\sigma_{f}\text{ and }\int
dx\,|g^{\prime}(x)|=2\sqrt{\pi}\sigma_{g},$ (99)

for $|f^{\prime}(x)|\leq 1$ and $|g^{\prime}(x)|\leq 1$. Moreover, by the
convolution theorem,

$f*g=\mathcal{FT}^{-1}[\mathcal{FT}(f)\,\mathcal{FT}(g)],$ (100)

we can see that comparing to the trivial case where $g$ is a delta function,
and we have $f=\mathcal{FT}^{-1}[\mathcal{FT}(f)]$, the width of
$\mathcal{FT}(f)\,\mathcal{FT}(g)$ would decrease when the width of $g$ is no
longer zero, and hence $\sigma_{f*g}$ would increase. For example, if $f(x)$
and $g(x)$ are Gaussian distributions, then
$\sigma_{f*g}^{2}=\sigma_{f}^{2}+\sigma_{g}^{2}$.

  * •

Property 5: Convolution of a complex function, $h(x)=f(x)\,e^{ipx}$, and a
real function, $g(x)$, is

$\left(h*g\right)(y)=\int
e^{ip^{\prime}y}\,\tilde{f}(p^{\prime}-p)\,\tilde{g}(p^{\prime})\equiv
e^{ipY(y)}I_{1}(y)I_{2}(p),$ (101)

where $\tilde{f}=\mathcal{FT}[f]$, $\tilde{g}=\mathcal{FT}[g]$ and the width
of $I(y)\in\mathbb{R}$ is the same as that of $(f*g)(y)$.

By the convolution theorem,

$\displaystyle\mathcal{FT}[h*g]=\mathcal{FT}[h]\mathcal{FT}[g]=\int
dxf(x)e^{-i(p^{\prime}-p)x}\int
dxg(x)e^{-ip^{\prime}x}=\tilde{f}(p^{\prime}-p)\tilde{g}(p^{\prime}).$ (102)

Then by doing an inverse Fourier transformation from $p^{\prime}$ to $y$ on
Eq. (102), we arrive at Eq. (101). Furthermore, the width of the product of
two functions, $\tilde{f}(p^{\prime}-p)\tilde{g}(p^{\prime})$ is independent
of the parallel shift from $p$, i.e. the width of
$\tilde{f}(p^{\prime}-p)\tilde{g}(p^{\prime})$ is the same as
$\tilde{f}(p^{\prime})\tilde{g}(p^{\prime})$, which is the case where the
convolution is between $f$ and $g$. For example, if $f(x)$ and $g(x)$ are
Gaussian distributed, i.e.

$f(x)=\exp\left[\frac{-(x-\mu_{f})^{2}}{4\sigma_{f}}\right],\quad
g(x)=\exp\left[\frac{-(x-\mu_{g})^{2}}{4\sigma_{g}}\right],$ (103)

then $\tilde{f}(p^{\prime})=\exp\left(-ip^{\prime}\mu_{f}-p^{\prime
2}\sigma_{f}^{2}\right)$,
$\tilde{g}(p^{\prime})=\exp\left(-ip^{\prime}\mu_{g}-p^{\prime
2}\sigma_{g}^{2}\right)$, and

$\tilde{f}(p^{\prime}-p)\tilde{g}(p^{\prime})=e^{-ip^{\prime}(\mu_{f}+\mu_{g})}\exp\left[-\left(\sigma_{f}^{2}+\sigma_{g}^{2}\right)\left(p^{\prime}-\frac{p}{\Delta}\right)^{2}-\sigma_{f}^{2}p^{2}\left(1-\frac{1}{\Delta}\right)\right],$
(104)

where $\Delta=(\sigma_{f}^{2}+\sigma_{g}^{2})/\sigma_{f}^{2}$. We can see that
the width of Eq. (104) w.r.t. $p^{\prime}$ is independent of $p$. Therefore,
according to Eq. (101), the convolution of $h$ and $g$ is then to do a Fourier
transformation from $p^{\prime}$ to $y-(\mu_{f}+\mu_{g})$, giving us

$(h*g)(y)=e^{i\frac{p}{\Delta}(y-\mu_{f}-\mu_{g})}\exp\left[-\frac{(y-\mu_{f}-\mu_{g})^{2}}{4(\sigma_{f}^{2}+\sigma_{g}^{2})}-\sigma_{f}^{2}p^{2}\left(1-\frac{1}{\Delta}\right)\right].$
(105)

Hence, the width w.r.t. $y$ is $\sigma_{f}^{2}+\sigma_{g}^{2}$, which is the
same the width of $(f*g)(y)$ shown in Property 4. Moreover, when
$\sigma_{f}=1/\sigma_{p}$ and $\sigma_{g}=2\sigma_{x}$, it follows
$\Delta=1+4\sigma_{x}^{2}\sigma_{p}^{2}$, which agrees with Eq. (15).

  * •

Property 6: In Table 1 & 2 we classify how the width would evolve after
Fourier transformation, product of functions, convolution of real functions
and convolution of complex functions as in “property 5”. By the properties
above, we summarize the width evolution of the first three types in the Table
1, specifying its relation with the original function(s) and give the example
of assuming all original functions are Gaussian distributed.

As for the last type, two functions are generated under such combination,
namely, $H(y)$ and $I(p)$ in

$\bigg{|}\int dxf(x)e^{ipx}g(x-y)\bigg{|}=H(y)I(p).$ (106)

By property 5, we see that the width of $H(y)$ and $I(p)$ has properties shown
in Table 2.

Type (notation for $\sigma_{h}$) | Function relation | Width relation | Gaussian case
---|---|---|---
FT type ($\tilde{\sigma}_{f}$) | $h=\mathcal{FT}(f)$ | NC | $\sigma_{h}=\frac{1}{2\sigma_{f}}$
Product type ($\sigma_{fg}$) | $h=f\times g$ | PC, $\sigma_{h}<\\{\sigma_{f},\sigma_{g}\\}$ | $\frac{1}{\sigma_{h}^{2}}=\frac{1}{\sigma_{f}^{2}}+\frac{1}{\sigma_{g}^{2}}$
Convolution type I ($\sigma_{f*g}$) | $h=f*g$ | PC, $\sigma_{h}>\\{\sigma_{f},\sigma_{g}\\}$ | $\sigma_{h}^{2}=\sigma_{f}^{2}+\sigma_{g}^{2}$
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 $\sigma_{h}$ is negatively/positively correlated to $\sigma_{f}$ (and $\sigma_{g}$). Function (width notation), Type | Width relation | Gaussian case
---|---|---
$H(y)$ $(\sigma_{H})$, Convolution type I | PC , $\sigma_{H}>\\{\sigma_{f},\sigma_{g}\\}$ | $\sigma_{H}=\sigma_{f*g}$
$I(p)$ $(\sigma_{I})$, Convolution type II | NC | $\sigma_{I}=1/\sigma_{fg}$
Table 2:  Properties of the width of functions $H$ and $I$ for convolution
with an addition complex phase in terms of their origin functions $f$ and $g$,
giving rise to an additional term $I(p)$.

##  Appendix B Calculation of the Neutrino Flavor Transition Amplitude

###  B.1 Neutrinos represented by Entanglement States

In this subsection, we work out the details in Eq. (9), where the phase space
representing layer 1 is composed by $x=x_{2}-x_{1}$ and
$\mathbb{p}=\mathbb{q}-\mathbb{k}=\mathbb{k}^{\prime}-\mathbb{q}^{\prime}$,
where $x_{1}/x_{2}$ are the space-time coordinates of the production/detection
vertices and $\mathbb{q},\mathbb{k},\mathbb{q}^{\prime},\mathbb{k}^{\prime}$
are the momenta of the initial and final states of the production and
detection sites. Therefore, $x$ and $\mathbb{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. (2.1) to Eq. (9), and also the
relations in Fig. 18. We first include all the uncertainties following Eq.
(2.1):

$\displaystyle
A_{2,j}(T,\mathbb{L},\mathbb{P})=\int[dq]f_{Pi}(\mathbb{q})\int[dk]\,f^{*}_{Pf}(\mathbb{k})\int[dq^{\prime}]\,f_{Di}(\mathbb{q}^{\prime})\int[dk^{\prime}]f^{*}_{Df}(\mathbb{k}^{\prime})$
$\displaystyle\times\int d^{4}x_{1}\,g_{P}(x_{1})\int
d^{4}x_{2}\,g_{D}(x_{2})\int
d^{4}\,y_{2}\,M_{Dj}(q^{\prime},k^{\prime})e^{-i(q^{\prime}-k^{\prime})(y_{2}-x_{2})}$
$\displaystyle\times\int\frac{d^{4}p_{\nu}}{(2\pi)^{4}}\,\frac{\mathrlap{\\!\not{\phantom{p_{\nu}}}}p_{\nu}+m_{j}}{p_{\nu}^{2}-m_{j}^{2}+i\epsilon}e^{-ip_{\nu}(y_{1}-y_{2})}\int
d^{4}y_{1}\,M_{Pj}(q,k)e^{-i(q-k)(y_{1}-x_{1})}.$ (107)

For all $h=\\{q,k,q^{\prime},k^{\prime}\\}$, since the external states are on
the mass-shell, $h^{0}=E_{h}(\mathbb{h})=\sqrt{\mathbb{h}^{2}-m_{h}^{2}}$.
Furthermore, if the wavepackets are sharply peaked at the expectation value
$\langle h\rangle$, then by the saddle point approximation, we can write
$E_{h}(\mathbb{h})\simeq E_{h}+\mathbb{v}_{h}(\mathbb{h}-\mathbb{h}_{0})$,
where $E_{h}=E_{h}(\langle h\rangle)$. Then, for the production site, by doing
a change of variables: $\\{q,k\\}\rightarrow\\{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_{Pj}(q,k)$, i.e.

$\displaystyle\int\frac{d^{3}k}{(2\pi)^{3}}f^{*}_{Pf}(\mathbb{k})f_{Pi}(\mathbb{k}-\mathbb{p})M^{\prime}_{Pj}(\mathbb{p},\mathbb{k})e^{-i(y_{1}^{0}-x_{1}^{0})(E_{q}(\mathbb{k}-\mathbb{p})-E_{k}(\mathbb{k}))}$
$\displaystyle=F_{Pj}(\mathbb{p})F_{P}^{\prime}(y_{1}^{0}-x_{1}^{0})e^{i\xi(\mathbb{p})(y_{1}^{0}-x_{1}^{0})},$
(108)

where we write
$M^{\prime}_{Pj}(\mathbb{p},\mathbb{k})=M_{Pj}(k-p,k)|_{p^{0}=E_{q}(\mathbb{k}-\mathbb{p})-E_{k}(\mathbb{k}),\,k^{0}=E_{k}(\mathbb{k})}$
for convenience. For $f_{Pi}$ and $f_{Pj}$ being Gaussian functions with width
$\sigma_{q}$ and $\sigma_{k}$ respectively, and applying the saddle point
approximation, the momentum uncertainties from the external states at the
production site are

$F_{Pj}=M^{\prime}_{Pj}(\langle\mathbb{p}\rangle,\langle\mathbb{k}\rangle)\exp\left[\frac{-(\mathbb{p}-\mathbb{P})^{2}}{4(\sigma_{q}^{2}+\sigma_{k}^{2})}\right],$
(109)

where $\mathbb{P}=\langle\mathbb{q}\rangle-\langle\mathbb{k}\rangle$. The
other terms are

$\displaystyle
F_{P}^{\prime}(y_{1}^{0}-x_{1}^{0})=\exp\left[-(y_{1}^{0}-x_{1}^{0})^{2}\sigma_{qk}^{2}\,\mathbb{v}_{qk}^{2}\right],$
(110)
$\displaystyle\xi_{P}(\mathbb{p})=E_{q}-E_{k}-\mathbb{v}_{q}\mathbb{q}_{0}+\mathbb{v}_{k}\mathbb{k}_{0}+\frac{\mathbb{v}_{qk}^{2}}{\Delta_{qk}}\mathbb{P}+\mathbb{p}\left(\mathbb{v}_{q}+\frac{\mathbb{v}_{qk}^{2}}{\Delta_{qk}}\right),$
(111)

where $\mathbb{v}_{qk}=\mathbb{v}_{q}-\mathbb{v}_{k}$,
$\Delta_{qk}=(\sigma_{q}^{2}+\sigma_{k}^{2})/\sigma_{k}^{2}$, and
$\sigma_{kq}$ refers to the notation in Table 1. Then analogously for the
detection site, with
$\mathbb{p}^{\prime}=\mathbb{k}^{\prime}-\mathbb{q}^{\prime}$ and also
$\mathbb{P}^{\prime}=\langle\mathbb{k}^{\prime}\rangle-\langle\mathbb{q}^{\prime}\rangle$,
Eq. (107) becomes:

$\displaystyle A_{2,j}=\int d^{3}p\,F_{Pj}(\mathbb{p})\int
d^{3}p^{\prime}\,F_{Dj}(\mathbb{p}^{\prime})\int d^{4}x_{1}\,g_{P}(x_{1})\int
d^{4}x_{2}\,g_{D}(x_{2})$ $\displaystyle\times\int
d^{4}y_{1}e^{-i(y_{1}^{0}-x_{1}^{0})\xi_{P}(\mathbb{p})+i(\mathbb{y}_{1}-\mathbb{x}_{1})\mathbb{p}}F^{\prime}_{P}(y_{1}^{0}-x_{1}^{0})\int
d^{4}y_{2}e^{-i(y_{2}^{0}-x_{2}^{0})\xi_{D}(\mathbb{p}^{\prime})+i(\mathbb{y}_{2}-\mathbb{x}_{2})\mathbb{p}^{\prime}}F^{\prime}_{D}(y_{2}^{0}-x_{2}^{0})$
$\displaystyle\times\int\frac{d^{4}p_{\nu}}{(2\pi)^{4}}\frac{\mathrlap{\\!\not{\phantom{p_{\nu}}}}p_{\nu}+m_{j}}{p_{\nu}^{2}-m_{j}^{2}+i\epsilon}\,e^{-ip_{\nu}(y_{1}-y_{2})}.$
(112)

The integration over $\int d^{3}y_{1}$ and $\int d^{3}y_{2}$ gives rise to
$\delta^{3}(\mathbb{p}-\mathbb{p}_{\nu})$ and
$\delta^{3}(\mathbb{p}^{\prime}-\mathbb{p}_{\nu})$, respectively, while the
integration over $\int dy_{1}^{0}$ returns

$\int
dy_{1}^{0}e^{-iy_{1}^{0}(\xi(\mathbb{p})-p_{\nu}^{0})}F_{P}^{\prime}(y_{1}^{0}-x_{1}^{0})=e^{-ix_{1}^{0}(\xi(\mathbb{p})-p_{\nu}^{0})}\tilde{F}^{\prime}_{P}(\xi(\mathbb{p})-p_{\nu}^{0}),$
(113)

where $\tilde{F}^{\prime}_{P}$ is the Fourier transformation of
$F^{\prime}_{P}$. After analogous calculations for the detection site, Eq.
(112) finally takes the form of the layer moving operator, with weighting
function $F_{j}(\mathbb{p};\mathbb{P})G(x,X)$:

$A_{2,j}=\int d^{3}p\,\int
d^{4}x\,F_{j}(\mathbb{p};\mathbb{P})G(x;X)\,A_{1,j}(x,\mathbb{p}),$ (114)

where $x=x_{2}-x_{1}$,
$F_{j}(\mathbb{p},\mathbb{P})=F_{Pj}(\mathbb{p})F_{Dj}(\mathbb{p})$,

$G(x)=\int d^{4}x_{2}g_{P}(x_{2}-x)g_{D}(x_{2}),$ (115)

and the first layer transition amplitude is

$A_{1,j}(x,\mathbb{p})=e^{i\mathbb{p}\mathbb{x}}\,\int
dp_{\nu}^{0}e^{-ip_{\nu}^{0}x^{0}}\tilde{F}^{\prime}_{P}(\xi_{P}(\mathbb{p})-p_{\nu}^{0})\tilde{F}^{\prime}_{D}(\xi_{D}(\mathbb{p})-p_{\nu}^{0})\frac{\mathrlap{\\!\not{\phantom{p_{\nu}}}}p_{\nu}+m_{j}}{p_{\nu}^{2}-m_{j}^{2}+i\epsilon}\bigg{|}_{\mathbb{p}_{\nu}=\mathbb{p}}.$
(116)

This appears as the collection of configuration of all the energetically
allowed states for some $\mathbb{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_{\nu}^{0}=E_{j}(\mathbb{p})\equiv\sqrt{\mathbb{p}^{2}+m_{j}^{2}}$, then we
can replace $A_{1,j}(x,\mathbb{p})\rightarrow
e^{-itE_{j}(\mathbb{p})+i\mathbb{x}\mathbb{p}}$ and

$F_{j}(\mathbb{p})\rightarrow
F_{j}(\mathbb{p})\tilde{F}^{\prime}_{P}(\xi(\mathbb{p})-E_{j}(\mathbb{p}))\tilde{F}^{\prime}_{D}(\xi^{\prime}(\mathbb{p})-E_{j}(\mathbb{p})),$
(117)

in Eq. (114). Therefore, either from the calculation we have shown, or more
efficiently from Fig. 18, the width of $F_{j}$, $\sigma_{p}$ is

$\frac{1}{\sigma_{p}^{2}}=\frac{1}{\sigma_{q}^{2}+\sigma_{k}^{2}}+\frac{1}{\sigma_{q^{\prime}}^{2}+\sigma_{k^{\prime}}^{2}}+\frac{|\mathbb{v}_{q}-\mathbb{v}_{j}+\mathbb{v}_{qk}^{2}/\Delta_{qk}|}{\sigma_{qk}^{2}\mathbb{v}_{qk}^{2}}+\frac{|\mathbb{v}_{q^{\prime}}-\mathbb{v}_{j}+\mathbb{v}_{q^{\prime}k^{\prime}}^{2}/\Delta_{q^{\prime}k^{\prime}}|}{\sigma_{q^{\prime}k^{\prime}}^{2}\mathbb{v}_{q^{\prime}k^{\prime}}^{2}}.$
(118)

Figure 18:  The final (width of the) weighting functions ($G(x;X)$ and
$F_{j}(\mathbb{p},\mathbb{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 $\sigma_{x}$ and $\sigma_{p}$ in the main text.

Next, to move $A_{1,j}$ onto the next layer, we first integrate out the
coordinate space $\int d^{4}x$, for $X=(T,\mathbb{L})$, and have

$A_{2,j}=\int
d^{3}p\,e^{-iTE_{j}(\mathbb{p})+i\mathbb{L}\mathbb{p}}F_{j}(\mathbb{p};\mathbb{P})\tilde{G}(\mathbb{p}).$
(119)

Then as long as one of the functions $F_{P}$, $F_{D}$,
$\tilde{F}_{P}^{\prime}$, $\tilde{F}_{D}^{\prime}$ or $\tilde{G}$ is sharply
peaked, we can apply the saddle point approximation at $\mathbb{P}_{j}$, such
that

$\frac{d}{d\mathbb{p}}F_{j}(\mathbb{p};\mathbb{P})\tilde{G}(\mathbb{p})\bigg{|}_{\mathbb{p}=\mathbb{P}_{j}}=0,$
(120)

then the approximation gives $E_{j}(\mathbb{p})\simeq
E_{j}+\mathbb{v}_{j}(\mathbb{p}-\mathbb{P}_{j})$, resulting in the final form
of the second layer transition amplitude as

$A_{2,j}=e^{-iE_{j}T+i\mathbb{P}_{j}\mathbb{L}}\hat{\Phi}_{j}(\mathbb{L}_{j},\mathbb{P}_{j}),$
(121)

where $\mathbb{L}_{j}=\mathbb{L}-\mathbb{v}_{j}T$. In fact, with such
approximation, and taking $G(x;X)$ as

$G_{x}(x,X)=\exp\left[-\frac{(t-T)^{2}}{4\sigma_{t}^{2}}-\frac{(\mathbb{x}-\mathbb{L})^{2}}{4\sigma_{\mathbb{x}}^{2}}\right],$
(122)

the function after the integration of $x$ (a Fourier transformation to the
momentum space) is:

$\displaystyle\int dt\int
d^{3}xe^{-it(E_{j}+\mathbb{v}_{j}\mathbb{p}-\mathbb{v}_{j}\mathbb{P}_{j})+i\mathbb{x}\mathbb{p}}G_{x}(x,X)$
(123)
$\displaystyle=e^{-i(E_{j}-\mathbb{v}_{j}\mathbb{P}_{j})T+i\mathbb{p}(\mathbb{L}-\mathbb{v}_{j}T)}\exp\left[-(\sigma_{t}^{2}\mathbb{v}_{j}^{2}+\sigma_{\mathbb{x}}^{2})(\mathbb{p}-\tilde{m}_{j})^{2}\right],$
(124)

where
$\tilde{m}_{j}=m_{j}\mathbb{v}_{j}\sigma_{t}^{2}/(\sigma_{t}^{2}\mathbb{v}_{j}^{2}+\sigma_{\mathbb{x}}^{2})$.
Therefore, Eq. (114) can be written as

$e^{-iE_{j}T+\mathbb{v}_{j}\mathbb{P}_{j}T}\int dp^{3}\int
dx^{3}\,e^{i\mathbb{x}\mathbb{p}}G(\mathbb{x};\mathbb{L}_{j})F_{j}(\mathbb{p};\mathbb{P}_{j}),$
(125)

where $G$ has width
$\sigma_{x}^{2}=\sigma_{t}^{2}\mathbb{v}_{j}^{2}+\sigma_{\mathbb{x}}^{2}$ and
is centred at $\mathbb{L}_{j}=\mathbb{L}-\mathbb{v}_{j}T$; $F_{j}$ has width
$\sigma_{p}$ and centred at $\mathbb{P}^{\prime}_{j}=\Delta\,\mathbb{P}_{j}$,
where $\Delta=1+4\sigma_{x}^{2}\sigma_{p}^{2}$, such that the saddle point
from Eq. (120) is at $\mathbb{P}_{j}$. In general, if all the input
distributions are Gaussian distributed and with the saddle point
approximation, $\mathbb{x}$ and $\mathbb{p}$ both linear dependent, $F_{j}$
and $G$ will also be Gaussian distributed with some width $\sigma_{p}$ and
$\sigma_{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_{\nu}$ 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),

$\displaystyle A_{2,j}(T,\mathbb{L},\mathbb{P})=\int d^{4}y_{1}\int
d^{4}y_{2}\int\frac{d^{4}p_{\nu}}{(2\pi)^{4}}\Delta(p_{\nu})e^{-ip_{\nu}(y_{1}-y_{2})}$
$\displaystyle\times\int
d^{4}x_{1}\,g_{P}(x_{1})\tilde{F}_{P}(y_{1}-x_{1})e^{iP_{P}(y_{1}-x_{1})}\int
d^{4}x_{2}\,g_{D}(x_{2})\tilde{F}_{D}(y_{2}-x_{2})e^{-iP_{D}(y_{2}-x_{2})},$
(126)

where $\Delta(p_{\nu})$ is the neutrino propagator in the momentum space. Here
we write the wavepackets at the production and detection site in coordinate
space as:

$\displaystyle\tilde{F}_{Pj}(y_{1}-x_{1})\,e^{iP_{P}(y_{1}-x_{1})}\simeq\int
d^{3}pF_{Pj}(\mathbb{p})e^{-i(y_{1}^{0}-x_{1}^{0})\xi_{P}(\mathbb{p})+i(\mathbb{y}_{1}-\mathbb{x}_{1})\mathbb{p}},$
(127)
$\displaystyle\tilde{F}_{Dj}(y_{2}-x_{2})\,e^{-iP_{D}(y_{2}-x_{2})}\simeq\int
d^{3}p^{\prime}F_{Dj}(\mathbb{p}^{\prime})e^{-i(y_{2}^{0}-x_{2}^{0})\xi_{D}(\mathbb{p}^{\prime})+i(\mathbb{y}_{2}-\mathbb{x}_{2})\mathbb{p}^{\prime}},$
(128)

where $P_{P}$ and $P_{D}$ are the saddle point of $F_{Pj}(\mathbb{p})$ and
$F_{Dj}(\mathbb{p}^{\prime})$, respectively. If there is negligible energy
loss during the neutrino propagation, then $P_{P}=P_{D}\simeq P$. Next, after
the integration over $x_{1}$ and $x_{2}$, Eq. (B.2) becomes

$\int d^{4}y\left\\{\int
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)]}\right\\}e^{iPy}\tilde{\Delta}(y),$
(129)

where $y=y_{1}-y_{2}$, and $\tilde{\Delta(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.

$\displaystyle I_{P}(y_{1};P)e^{iP[y_{1}-x_{P}(y_{1})]}=\int
d^{4}x_{1}\,g_{P}(x_{1})\tilde{F}_{Pj}(y_{1}-x_{1})F_{P}^{\prime}(y_{1}^{0}-x_{1}^{0})e^{iP(y_{1}-x_{1})}$
(130) $\displaystyle I_{D}(y_{2};P)e^{-iP[y_{2}-x_{D}(y_{2})]}=\int
d^{4}x_{2}\,g_{D}(x_{2})\tilde{F}_{Dj}(y_{2}-x_{2})F_{D}^{\prime}(y_{2}^{0}-x_{2}^{0})e^{-iP(y_{2}-x_{2})}.$
(131)

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}*\tilde{F}_{Pj}F_{P}^{\prime})*(g_{D}*\tilde{F}_{Dj}F_{D}^{\prime}))(y)$.
Moreover, with the association and commutation property for convolution, we
can rewrite the total coordinate width as $(G*\tilde{F}_{Pj}^{\rm
tot}*\tilde{F}_{Dj}^{\rm tot})(y)$, where $G$ is in Eq. (115),
$\tilde{F}_{Pj}^{\rm tot}=\tilde{F}_{Pj}F_{P}^{\prime}$ and
$\tilde{F}_{Dj}^{\rm tot}=\tilde{F}_{Dj}F_{D}^{\prime}$.

##  Appendix C Factorization Condition

In this appendix, we derive the conditions under which a function can be
factorized out of some integration, i.e.

$\int_{-\infty}^{\infty}dX_{2}S(X_{2})Y(X_{2};X_{3})\simeq
S(X_{3})\int_{-\infty}^{\infty}dX_{2}Y(X_{2};X_{3}),$ (132)

where $S(X_{2})$ and $Y(X_{2};X_{3})$ are both localized, i.e.
$S(X_{2})=Y(X_{2})=0$ as $X_{2}\rightarrow\pm\infty$. This is useful for
describing state decoherence as a phase wash-out effect on the Wigner-PS in
Eq. (64), and disentangling $\sigma_{E}$ and $\sigma_{L}$ in Eq. (76).
Intuitively, the condition where $S(X_{2})$ can be factorized out of the
integral as Eq. (132) is when the width of $S(X_{2})$ is much larger then
$Y(X_{2})$, since the product of the two functions would be dominated by the
function which is more localized. Nevertheless, in order to see if this
condition is sufficient and to have a more concrete idea, we derive Eq. (132)
as follows:

$\displaystyle LHS=$
$\displaystyle-\int_{-\infty}^{\infty}dX_{2}\frac{dS(X_{2})}{dX_{2}}\int_{-\infty}^{X_{2}}dX_{2}^{\prime}Y(X_{2}^{\prime};X_{3})$
(133) $\displaystyle\simeq$
$\displaystyle-\int_{X_{3}-\Lambda}^{\infty}dX_{2}\frac{dS(X_{2})}{dX_{2}}\int_{X_{3}-\Lambda}^{X_{2}}dX_{2}^{\prime}Y(X_{2}^{\prime};X_{3})$
(134) $\displaystyle\simeq$
$\displaystyle-\int^{\infty}_{X_{3}-\Lambda}dX_{2}\frac{dS(X_{2})}{dX_{2}}\int_{-\infty}^{\infty}dX_{2}^{\prime}Y(X_{2}^{\prime};X_{3})$
$\displaystyle+\int^{X_{3}+\Lambda}_{X_{3}-\Lambda}dX_{2}\frac{dS(X_{2})}{dX_{2}}\int_{X_{2}}^{X_{3}+\Lambda}dX_{2}^{\prime}Y(X_{2}^{\prime};X_{3})=RHS.$
(135)

Eq. (133) is achieved by doing integration by parts with the boundary terms
vanishing due to the localization property of $S(X_{2})$.

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; $\Lambda$ is the cutoff value in Eq. (134).

In Eq. (134) and Eq. (135), we use the localization property of $Y(X_{2})$ to
make a cut at $\Lambda$ w.r.t. $X_{3}$ such that the cumulative distribution
function
$\int_{-\infty}^{X_{3}-\Lambda}dX_{2}^{\prime}Y(X_{2}^{\prime};X_{3})=0$ and
$\int^{X_{3}+\Lambda}_{-\infty}dX_{2}^{\prime}Y(X_{2}^{\prime};X_{3})$
converges to a constant, so that
$\int^{\infty}_{X_{3}+\Lambda}dX_{2}^{\prime}Y(X_{2}^{\prime};X_{3})=0$, as
shown in Fig. 19. Finally, if $S(X_{2})$ varies slowly in the interval
$(X_{3}-\Lambda,X_{3}+\Lambda)$, i.e. the width of $S(X_{2})$($\sigma_{S}$),
is much larger than that of $Y(X_{2})$($\Lambda$), then the latter part in Eq.
(135) can be neglected and we arrive at the $RHS$ of Eq. (132), where
$S(X_{3}-\Lambda)\rightarrow S(X_{3})$ can be taken out of the integral if
$X_{3}\gg\Lambda$. Furthermore, $Y(X_{2};X_{3})$ could be any distribution as
long as it is localized, even if it includes a non-zero phase term, as we can
see from the orange and purple lines in Fig. 19.

##  Appendix 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 $x_{1}$ and $x_{2}$,
contributing neutrino flux $A$ and $B$, then the phase decoherence term in Eq.
(73) is simply

$\Phi_{jk}=e^{-i\alpha_{jk}L_{3}}\left(a\,e^{i\alpha_{jk}x_{1}}+b\,e^{i\alpha_{jk}x_{2}}\right)\equiv\phi^{(1)}\,e^{i\alpha_{jk}(x_{\rm
eff}^{(1)}-L_{3})},$ (136)

where $a=A/(A+B)$, $b=B/(A+B)$ and $0\leq c\leq 1$, are real. This requires

$A\sin\left(\alpha_{jk}\delta_{1}\right)+B\sin\left(\alpha_{jk}\delta_{2}\right)=0,$
(137)

where $x_{1}=x_{\rm eff}^{(1)}+\delta_{1}$ and $x_{2}=x_{\rm
eff}^{(1)}+\delta_{2}$. Hence, by solving Eq. (137) for $\Delta
x^{(1)}=x_{1}-x_{2}=\delta_{1}-\delta_{2}$, we have

$\delta_{1}\equiv f_{jk}(\Delta
x,\frac{a}{b})=\frac{1}{\alpha_{jk}}\tan^{-1}\left[\frac{-\sin(\alpha_{jk}\Delta
x)}{\frac{a}{b}+\cos(\alpha_{jk}\Delta x)}\right],$ (138)

and thus

$\displaystyle x_{\rm eff}^{(1)}=x_{1}+f_{jk}(\Delta x^{(1)},\frac{a}{b}),$
$\displaystyle\phi^{(1)}=a\cos\left[\alpha_{jk}f_{jk}(\Delta
x^{(1)},\frac{a}{b})\right]+b\cos\left[\alpha_{jk}(f_{jk}(\Delta
x^{(1)},\frac{a}{b})+\Delta x^{(1)})\right],$ (139)

indicating that the damping term $\phi\leq 1$ as excepted since $a+b=1$. In
the case where the detector is placed far from all the sources,
$x_{1},x_{2}\gg\Delta x$, then $x_{\rm eff}=x_{2}$, and $c=a+b\cos(\Delta x)$.
Similarly, if there are three point like neutrino sources,

$\Phi_{jk}=e^{-i\alpha_{jk}L_{3}}\left(a\,e^{i\alpha_{jk}x_{1}}+b\,e^{i\alpha_{jk}x_{2}}+c\,e^{i\alpha_{jk}x_{3}}\right)\equiv\phi^{(2)}\,e^{i\alpha_{jk}(x_{\rm
eff}^{(2)}-L_{3})},$ (140)

the damping term ($\phi^{(2)}$) and the phase term ($x_{\rm eff}^{(2)}$) are
obtained by replacing $x_{1}\rightarrow x_{\rm eff}^{(1)}$, $\Delta
x\rightarrow x_{\rm eff}^{(1)}-x_{3}$, $a\rightarrow\phi^{(1)}$ and
$b\rightarrow c$ in Eq. (D), and so on for more point-like sources.

## References

  * [1] M. Sajjad Athar et al., “Status and perspectives of neutrino physics,” Prog. Part. Nucl. Phys. 124 (2022) 103947, arXiv:2111.07586 [hep-ph].
  * [2] V. N. Gribov and B. Pontecorvo, “Neutrino astronomy and lepton charge,” Phys. Lett. B 28 (1969) 493.
  * [3] Z. Maki, M. Nakagawa, and S. Sakata, “Remarks on the unified model of elementary particles,” Prog. Theor. Phys. 28 (1962) 870–880.
  * [4] B. Pontecorvo, “Neutrino Experiments and the Problem of Conservation of Leptonic Charge,” Zh. Eksp. Teor. Fiz. 53 (1967) 1717–1725.
  * [5] M. Schlosshauer, “Decoherence, the Measurement Problem, and Interpretations of Quantum Mechanics,” Rev. Mod. Phys. 76 (2004) 1267–1305, arXiv:quant-ph/0312059.
  * [6] H. D. Zeh, “On the interpretation of measurement in quantum theory,” Foundations of Physics 1 (1970) 69–76.
  * [7] J. Rich, “Quantum mechanics of neutrino oscillations,” Phys. Rev. D 48 (Nov, 1993) 4318–4325. https://link.aps.org/doi/10.1103/PhysRevD.48.4318.
  * [8] C. Giunti and C. W. Kim, “Coherence of neutrino oscillations in the wave packet approach,” Phys. Rev. D 58 (1998) 017301, arXiv:hep-ph/9711363.
  * [9] E. K. Akhmedov and J. Kopp, “Neutrino Oscillations: Quantum Mechanics vs. Quantum Field Theory,” JHEP 04 (2010) 008, arXiv:1001.4815 [hep-ph]. [Erratum: JHEP 10, 052 (2013)].
  * [10] K. Kiers, S. Nussinov, and N. Weiss, “Coherence effects in neutrino oscillations,” Phys. Rev. D 53 (1996) 537–547, arXiv:hep-ph/9506271.
  * [11] B. Kayser, “On the Quantum Mechanics of Neutrino Oscillation,” Phys. Rev. D 24 (1981) 110.
  * [12] M. Beuthe, “Oscillations of neutrinos and mesons in quantum field theory,” Phys. Rept. 375 (2003) 105–218, arXiv:hep-ph/0109119.
  * [13] C. Giunti, “Neutrino wave packets in quantum field theory,” JHEP 11 (2002) 017, arXiv:hep-ph/0205014.
  * [14] E. K. Akhmedov and A. Y. Smirnov, “Neutrino oscillations: Entanglement, energy-momentum conservation and QFT,” Found. Phys. 41 (2011) 1279–1306, arXiv:1008.2077 [hep-ph].
  * [15] N. V. Naumov D.V., “Quantum Field Theory of Neutrino Oscillations,” Phys. Part. Nuclei 51 (2020) 1–106.
  * [16] W. Grimus, “Revisiting the quantum field theory of neutrino oscillations in vacuum,” J. Phys. G 47 no. 8, (2020) 085004, arXiv:1910.13776 [hep-ph].
  * [17] G. Lindblad, “On the generators of quantum dynamical semigroups,” Commun. Math. Phys. 48 (1976) 119.
  * [18] E. Lisi, A. Marrone, and D. Montanino, “Probing possible decoherence effects in atmospheric neutrino oscillations,” Phys. Rev. Lett. 85 (2000) 1166–1169, arXiv:hep-ph/0002053.
  * [19] R. S. L. Hansen and A. Y. Smirnov, “The Liouville equation for flavour evolution of neutrinos and neutrino wave packets,” JCAP 12 (2016) 019, arXiv:1610.00910 [hep-ph].
  * [20] F. Benatti and R. Floreanini, “Open system approach to neutrino oscillations,” JHEP 02 (2000) 032, arXiv:hep-ph/0002221.
  * [21] T. Stirner, G. Sigl, and G. Raffelt, “Liouville term for neutrinos: Flavor structure and wave interpretation,” JCAP 05 (2018) 016, arXiv:1803.04693 [hep-ph].
  * [22] J. a. A. B. Coelho and W. A. Mann, “Decoherence, matter effect, and neutrino hierarchy signature in long baseline experiments,” Phys. Rev. D 96 no. 9, (2017) 093009, arXiv:1708.05495 [hep-ph].
  * [23] A. L. G. Gomes, R. A. Gomes, and O. L. G. Peres, “Quantum decoherence and relaxation in neutrinos using long-baseline data,” arXiv:2001.09250 [hep-ph].
  * [24] Y. Farzan, T. Schwetz, and A. Y. Smirnov, “Reconciling results of LSND, MiniBooNE and other experiments with soft decoherence,” JHEP 07 (2008) 067, arXiv:0805.2098 [hep-ph].
  * [25] B. J. P. Jones, “Dynamical pion collapse and the coherence of conventional neutrino beams,” Phys. Rev. D 91 no. 5, (2015) 053002, arXiv:1412.2264 [hep-ph].
  * [26] E. P. Wigner, “On the quantum correction for thermodynamic equilibrium,” Phys. Rev. 40 (1932) 749.
  * [27] A. Vlasenko, G. M. Fuller, and V. Cirigliano, “Neutrino Quantum Kinetics,” Phys. Rev. D 89 no. 10, (2014) 105004, arXiv:1309.2628 [hep-ph].
  * [28] M. Blasone and G. Vitiello, “Quantum field theory of fermion mixing,” Annals Phys. 244 (1995) 283–311, arXiv:hep-ph/9501263. [Erratum: Annals Phys. 249, 363–364 (1996)].
  * [29] C. Giunti, C. W. Kim, and U. W. Lee, “Remarks on the weak states of neutrinos,” Phys. Rev. D 45 (Apr, 1992) 2414–2420. https://link.aps.org/doi/10.1103/PhysRevD.45.2414.
  * [30] C. Giunti, “Fock states of flavor neutrinos are unphysical,” Eur. Phys. J. C 39 (2005) 377–382, arXiv:hep-ph/0312256.
  * [31] A. Tureanu, “Comment on the Comment on the paper ”Can oscillating neutrino states be formulated universally?”,” arXiv:2005.02219 [hep-ph].
  * [32] B. d. S. L. Torres, T. R. Perche, A. G. S. Landulfo, and G. E. A. Matsas, “Neutrino flavor oscillations without flavor states,” Phys. Rev. D 102 no. 9, (2020) 093003, arXiv:2009.10165 [hep-ph].
  * [33] L. Smaldone and G. Vitiello, “Neutrino Mixing and Oscillations in Quantum Field Theory: A Comprehensive Introduction,” Universe 7 no. 12, (2021) 504, arXiv:2111.11809 [hep-th].
  * [34] M. Blasone and L. Smaldone, “A note on oscillating neutrino states in quantum field theory,” Mod. Phys. Lett. A 35 no. 38, (2020) 2050313, arXiv:2004.04739 [hep-ph].
  * [35] E. Akhmedov, “Neutrino oscillations in matter: from microscopic to macroscopic description,” JHEP 02 (2021) 107, arXiv:2010.07847 [hep-ph].
  * [36] I. Esteban, M. C. Gonzalez-Garcia, M. Maltoni, T. Schwetz, and A. Zhou, “The fate of hints: updated global analysis of three-flavor neutrino oscillations,” JHEP 09 (2020) 178, arXiv:2007.14792 [hep-ph].
  * [37] A. de Gouvea, V. de Romeri, and C. A. Ternes, “Probing neutrino quantum decoherence at reactor experiments,” JHEP 08 (2020) 018, arXiv:2005.03022 [hep-ph].
  * [38] A. de Gouvêa, V. De Romeri, and C. A. Ternes, “Combined analysis of neutrino decoherence at reactor experiments,” JHEP 06 (2021) 042, arXiv:2104.05806 [hep-ph].
  * [39] JUNO Collaboration, J. Wang et al., “Damping signatures at JUNO, a medium-baseline reactor neutrino oscillation experiment,” arXiv:2112.14450 [hep-ex].
  * [40] T. Stuttard and M. Jensen, “Neutrino decoherence from quantum gravitational stochastic perturbations,” Phys. Rev. D 102 no. 11, (2020) 115003, arXiv:2007.00068 [hep-ph].
  * [41] D. Boriero, D. J. Schwarz, and H. Velten, “Flavour composition and entropy increase of cosmological neutrinos after decoherence,” Universe 5 no. 10, (2019) 203, arXiv:1704.06139 [astro-ph.CO].
  * [42] A. De Gouvêa, I. Martinez-Soler, Y. F. Perez-Gonzalez, and M. Sen, “Fundamental physics with the diffuse supernova background neutrinos,” Phys. Rev. D 102 (2020) 123012, arXiv:2007.13748 [hep-ph].
  * [43] J. Carpio, E. Massoni, and A. M. Gago, “Revisiting quantum decoherence for neutrino oscillations in matter with constant density,” Phys. Rev. D 97 no. 11, (2018) 115017, arXiv:1711.03680 [hep-ph].
  * [44] J. A. Carpio, E. Massoni, and A. M. Gago, “Testing quantum decoherence at DUNE,” Phys. Rev. D 100 no. 1, (2019) 015035, arXiv:1811.07923 [hep-ph].
  * [45] G. Balieiro Gomes, D. V. Forero, M. M. Guzzo, P. C. De Holanda, and R. L. N. Oliveira, “Quantum Decoherence Effects in Neutrino Oscillations at DUNE,” Phys. Rev. D 100 no. 5, (2019) 055023, arXiv:1805.09818 [hep-ph].
  * [46] P. Coloma, J. Lopez-Pavon, I. Martinez-Soler, and H. Nunokawa, “Decoherence in Neutrino Propagation Through Matter, and Bounds from IceCube/DeepCore,” Eur. Phys. J. C 78 no. 8, (2018) 614, arXiv:1803.04438 [hep-ph].
  * [47] G. L. S. R. N. Sen, “Fiber bundles in quantum physics,” J. Math. Phys. 43 (2002) 1323–1339.
  * [48] W. B. Case, “Wigner functions and Weyl transforms for pedestrians,” Am.J.Phys. 76(10) (2008) 937.
  * [49] R. Jacob and R. G. Sachs, “Mass and Lifetime of Unstable Particles,” Phys. Rev. 121 (1961) 350.
  * [50] E. Akhmedov, D. Hernandez, and A. Smirnov, “Neutrino production coherence and oscillation experiments,” JHEP 04 (2012) 052, arXiv:1201.4128 [hep-ph].
  * [51] J. Kersten and A. Y. Smirnov, “Decoherence and oscillations of supernova neutrinos,” Eur. Phys. J. C 76 no. 6, (2016) 339, arXiv:1512.09068 [hep-ph].
  * [52] M. Martini, M. Ericson, and G. Chanfray, “Energy reconstruction effects in neutrino oscillation experiments and implications for the analysis,” Phys. Rev. D 87 no. 1, (2013) 013009, arXiv:1211.1523 [hep-ph].
  * [53] V. De Romeri, E. Fernandez-Martinez, and M. Sorel, “Neutrino oscillations at DUNE with improved energy reconstruction,” JHEP 09 (2016) 030, arXiv:1607.00293 [hep-ph].
  * [54] E. K. Akhmedov and A. Y. Smirnov, “Paradoxes of neutrino oscillations,” Phys. Atom. Nucl. 72 (2009) 1363–1381, arXiv:0905.1903 [hep-ph].
  * [55] C. Giunti and C. W. Kim, “Quantum mechanics of neutrino oscillations,” Found. Phys. Lett. 14 no. 3, (2001) 213–229, arXiv:hep-ph/0011074.
  * [56] Y. Farzan and A. Y. Smirnov, “Coherence and oscillations of cosmic neutrinos,” Nucl. Phys. B 805 (2008) 356–376, arXiv:0803.0495 [hep-ph].
  * [57] A. Esteban-Pretel, S. Pastor, R. Tomas, G. G. Raffelt, and G. Sigl, “Decoherence in supernova neutrino transformations suppressed by deleptonization,” Phys. Rev. D 76 (2007) 125018, arXiv:0706.2498 [astro-ph].
  * [58] E. Akhmedov, J. Kopp, and M. Lindner, “Collective neutrino oscillations and neutrino wave packets,” JCAP 09 (2017) 017, arXiv:1702.08338 [hep-ph].
  * [59] D. Hooper, D. Morgan, and E. Winstanley, “Probing quantum decoherence with high-energy neutrinos,” Phys. Lett. B 609 (2005) 206–211, arXiv:hep-ph/0410094.
  * [60] NuSTEC Collaboration, L. Alvarez-Ruso et al., “NuSTEC White Paper: Status and challenges of neutrino–nucleus scattering,” Prog. Part. Nucl. Phys. 100 (2018) 1–68, arXiv:1706.03621 [hep-ph].
  * [61] B. Kayser and J. Kopp, “Testing the Wave Packet Approach to Neutrino Oscillations in Future Experiments,” arXiv:1005.4081 [hep-ph].
  * [62] M. Blennow, T. Ohlsson, and W. Winter, “Damping signatures in future neutrino oscillation experiments,” JHEP 06 (2005) 049, arXiv:hep-ph/0502147.
  * [63] Particle Data Group Collaboration, P. A. Zyla et al., “Review of Particle Physics,” PTEP 2020 no. 8, (2020) 083C01.
  * [64] RENO Collaboration, G. Bak et al., “Measurement of Reactor Antineutrino Oscillation Amplitude and Frequency at RENO,” Phys. Rev. Lett. 121 no. 20, (2018) 201801, arXiv:1806.00248 [hep-ex].
  * [65] RENO Collaboration, S. H. Seo et al., “Spectral Measurement of the Electron Antineutrino Oscillation Amplitude and Frequency using 500 Live Days of RENO Data,” Phys. Rev. D 98 no. 1, (2018) 012002, arXiv:1610.04326 [hep-ex].
  * [66] RENO Collaboration, J. K. Ahn et al., “RENO: An Experiment for Neutrino Oscillation Parameter $\theta_{13}$ Using Reactor Neutrinos at Yonggwang,” arXiv:1003.1391 [hep-ex].
  * [67] K. Abe et al., “Letter of Intent: The Hyper-Kamiokande Experiment — Detector Design and Physics Potential —,” arXiv:1109.3262 [hep-ex].
  * [68] DUNE Collaboration, B. Abi et al., “Deep Underground Neutrino Experiment (DUNE), Far Detector Technical Design Report, Volume II: DUNE Physics,” arXiv:2002.03005 [hep-ex].
  * [69] S. Ajimura et al., “Technical Design Report (TDR): Searching for a Sterile Neutrino at J-PARC MLF (E56, JSNS2),” arXiv:1705.08629 [physics.ins-det].
  * [70] D. Baxter et al., “Coherent Elastic Neutrino-Nucleus Scattering at the European Spallation Source,” JHEP 02 (2020) 123, arXiv:1911.00762 [physics.ins-det].
  * [71] J. Alonso et al., “Expression of Interest for a Novel Search for CP Violation in the Neutrino Sector: DAEdALUS,” arXiv:1006.0260 [physics.ins-det].
  * [72] R. Harnik, K. J. Kelly, and P. A. N. Machado, “Prospects of Measuring Oscillated Decay-at-Rest Neutrinos at Long Baselines,” Phys. Rev. D 101 no. 3, (2020) 033008, arXiv:1911.05088 [hep-ph].
  * [73] A. Strumia and F. Vissani, “Precise quasielastic neutrino/nucleon cross-section,” Phys. Lett. B 564 (2003) 42–54, arXiv:astro-ph/0302055.
  * [74] S. Ajimura et al., “Proposal: JSNS2-II,” arXiv:2012.10807 [hep-ex].
  * [75] COHERENT Collaboration, D. Akimov et al., “First Measurement of Coherent Elastic Neutrino-Nucleus Scattering on Argon,” Phys. Rev. Lett. 126 no. 1, (2021) 012002, arXiv:2003.10630 [nucl-ex].
  * [76] ESSnuSB Collaboration, A. Alekou et al., “Updated physics performance of the ESSnuSB experiment: ESSnuSB collaboration,” Eur. Phys. J. C 81 no. 12, (2021) 1130, arXiv:2107.07585 [hep-ex].
  * [77] E. Kolbe, K. Langanke, and P. Vogel, “Estimates of weak and electromagnetic nuclear decay signatures for neutrino reactions in Super-Kamiokande,” Phys. Rev. D 66 (2002) 013007.
  * [78] A. Bodek and J. L. Ritchie, “Fermi Motion Effects in Deep Inelastic Lepton Scattering from Nuclear Targets,” Phys. Rev. D 23 (1981) 1070.