SM

    Standard Model
CPV

    CP violating
CPC

    CP conserving
EH

    Euler-Heisenberg
BSM

    beyond the Standard Model
COM

    center of mass
QED

    Quantum Electrodynamics
EFT

    Effective Field Theory
ALP

    axion-like particle
EOM

    Equation of Motion
RHS

    right hand side
LHS

    left hand side
FP

    Fabry-Perot
EM

    Electromagnetism
SRF

    superconducting radio frequency
SNR

    signal to noise ratio
EP

    Equivalence Principle
CL

    confidence level
DM

    dark matter
ULDM

    ultra-light dark matter
VEV

    vacuum expectation value
EWSB

    Electroweak Symmetry Breaking

Time Dependent CP-even and CP-odd Signatures of Scalar Ultra-light Dark Matter
in Neutrino Oscillations

Marta Losada1, Yosef Nir2, Gilad Perez2, Inbar Savoray2 and Yogev Shpilman2

1New York University Abu Dhabi, PO Box 129188, Saadiyat Island, Abu Dhabi,
United Arab Emirates
2Department of Particle Physics and Astrophysics,
Weizmann Institute of Science, Rehovot 7610001, Israel

<EMAIL_ADDRESS>yosef.nir,gilad.perez,inbar.savoray,
<EMAIL_ADDRESS>
Scalar ultra-light dark matter (ULDM) interacting with neutrinos can induce,
under certain conditions, time-dependent modifications to neutrino oscillation
probabilities. The limit in which the ULDM perturbation can be treated as
constant throughout the neutrino propagation time has been addressed by
several previous works. We complement these by systematically analyzing the
opposite limit – accounting for the temporal-variations of the ULDM potential
by solving time-dependent Schrödinger equations. In particular, we study a
novel two-generations-like CP violating (CPV) signature unique to rapidly
oscillating ULDM. We derive the leading order, time-dependent, corrections to
the oscillation probabilities, both for CP conserving (CPC) and CPV couplings,
and explain how they can be measured in current and future experiments.

##  1 Introduction

Scalar ULDM with mass $m_{\phi}\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip
1.0pt$\sim$}\hss}\raise 1.0pt\hbox{$<$}}{\rm eV}$, may leave signatures in
neutrino oscillation experiments, that cannot be interpreted within the
Standard Model (SM) supplemented by neutrino masses and lepton mixing. The
scalar can be treated as a classical bosonic field that oscillates with time,

$\displaystyle\phi=\phi_{0}\sin(m_{\phi}(t-t_{0}))\,,$ (1)

with an arbitrary initial time $t_{0}$ and amplitude

$\displaystyle\phi_{0}=\frac{\sqrt{2\rho_{\phi}^{\oplus}}}{m_{\phi}}\sim 2\
{\rm GeV}\ \left(\frac{10^{-12}\ {\rm eV}}{m_{\phi}}\right),$ (2)

where $\rho_{\phi}^{\oplus}$ is the local ULDM density. Consider the effective
mass and $\phi$-Yukawa terms for the neutrinos, arising from dimension-five
and dimension-six terms in the Lagrangian, respectively. In the basis in which
the ULDM-independent mass matrix is diagonal, they read

$\displaystyle{\cal
L}_{m_{\nu}}=m_{i}\nu_{i}^{T}\nu_{i}+\hat{y}_{ij}\phi\nu_{i}^{T}\nu_{j}+{\rm
h.c}\,.$ (3)

Treating $\phi$ as a classical field, it modifies the neutrino mass matrix as

$\displaystyle(\hat{m}_{\nu})_{ij}=m_{i}\delta_{ij}+\hat{y}_{ij}\phi\,,$ (4)

thus inducing a time-dependent component to the neutrino propagation
Hamiltonian.

Some of the implications of the ULDM contribution to the neutrino mass matrix
have been studied in previous works [1, 2, 3, 4, 5, 6, 7, 8]. Most previous
analyses focused on ULDM candidates with typical oscillation times,
$\tau_{\phi}\equiv 2\pi/m_{\phi}$, much larger than the neutrino propagation
time, $\tau_{d}$, but smaller than the total run-time of the experiment,
$\tau_{e}$ [4, 1, 5, 6]. In this limit, the effective neutrino masses and
mixing do not change during the propagation from the source to the detector,
but do change between neutrino events over the duration of the experiment.
Depending on the value of $m_{\phi}$, the experiment time resolution
$\tau_{r}$, and the way one analyses the data, the ULDM effect can be either
time-resolved, or time-averaged. In the former case, the temporal modulation
could be directly observed in some spectral analysis of the data [5, 6]. In
the latter case, the oscillatory contributions are averaged out, and one must
use indirect methods for recovering evidence of the time dependence, e.g. via
the distortion of the dynamical range of the neutrino oscillation parameters
or of relations between oscillation amplitudes at different $L/E$ frequencies
[1, 7, 5, 6] ($L$ is the source to detector distance and $E$ is the neutrino
energy).

However, the experimental signatures of heavier ULDM, with
$\tau_{\phi}\lesssim\tau_{d}$, require a dedicated analysis. In the limit of
rapid oscillations, one must solve time-dependent equations of motion to
correctly account for the ULDM effects on neutrino oscillation probabilities,
and the picture of slowly varying oscillation parameters is no longer valid
[5, 9]. In addition, thus far, the majority of previous works discussing ULDM-
modified neutrino oscillations have only considered CPC couplings. CPV effects
in neutrino oscillations induced by ULDM-neutrino interactions have been
analyzed by authors of this manuscript in Ref. [6], for the case of slow
oscillations. Furthermore, past works have either studied specific flavor
structures for the ULDM couplings [2, 7, 6], or studied the “model-agnostic”
temporal-modulations of the mixing angles, mass differences and the CPV phase
of the PMNS matrix [1, 5, 6]. To the best of our knowledge, the relations
between the ULDM couplings and these variations of oscillation parameters have
not been explicitly discussed, and could be non trivial (both mathematically
and phenomenologically). Therefore, translating the conclusions of the
existing analyses to generic ULDM models could be challenging.

In this work, we fill in those gaps, and provide a more complete description
of the phenomenology of ULDM-neutrino couplings. We are interested in mapping
those couplings onto their signatures in neutrino oscillation experiments, for
a wide range of ULDM masses, covering both the slow oscillations and the rapid
oscillations limits. We present a generic prescription for estimating the
first-order, time-dependent, modifications to neutrino oscillations from ULDM
couplings to neutrinos, and the resulting bounds on them from current and
future experiments. While we study both CPC and CPV phenomena, we give special
attention to the implications of the CPV ULDM couplings. In particular, we
describe new CPV effects generated by them, that can be observable even in
two-generations neutrino systems, and cannot be addressed by analyses assuming
a constant neutrino Hamiltonian through propagation.

The paper is organized as follows. In Sec. 2, we describe the new CPV effect
in the two neutrino generations picture - which we refer to as the “2-$\nu$
CPV”. In Sec. 3, we derive the full analytical expressions for the leading
order modifications to neutrino oscillation probabilities induced by the ULDM,
both for CPC and CPV phenomena. In Sec. 4, we describe the method for
detecting such time-dependent modifications in neutrino oscillations data. In
Sec. 5, we show the current bounds, and the future-projected sensitivities,
for the ULDM parameter space. We conclude in Sec. 6.

##  2 $2-\nu$ CP violation

We are interested in the possibility of observing CPV in neutrino oscillations
induced by the time-dependent vacuum expectation value (VEV) of $\phi$. In the
usual two-generations picture, one can show that by field redefinitions, the
mass matrix $\hat{m}_{\nu}$ can always be brought to a form in which it has no
physical phases affecting neutrino oscillations (see e.g. [10, 11]), and thus
CPV cannot be observed in these processes. This can be seen by the simple
following argument. In the two-neutrino picture, the Equation of Motion (EOM)
is given by

$\displaystyle i\partial_{t}\left(\begin{matrix}\nu_{\alpha}\\\
\nu_{\beta}\end{matrix}\right)$ $\displaystyle={\cal
H}\left(t\right)\left(\begin{matrix}\nu_{\alpha}\\\
\nu_{\beta}\end{matrix}\right)=\left(\begin{matrix}{\cal H}_{11}&{\cal
H}_{12}e^{i\varphi}\\\ {\cal H}_{12}e^{-i\varphi}&{\cal
H}_{22}\end{matrix}\right)\left(\begin{matrix}\nu_{\alpha}\\\
\nu_{\beta}\end{matrix}\right)\,.$ (5)

Since the Hamiltonian ${\cal H}$ is Hermitian, ${\cal H}_{ij}$ and $\varphi$
above must be real. Let $\nu_{i}(t)$ be the time-evolved state, according to
Eq. (5), which was $\nu_{i}$ at $t=0$. If $\nu_{\alpha,\beta}$ are interaction
eigenstates, then the initial conditions at $t=0$ are zero for one of them,
and one for the other, and the probability for detecting the states
$\nu_{\alpha,\beta}$ at some later time $t$ is given by solving the EOM and
calculating

$\displaystyle
P_{\nu_{i}\rightarrow\nu_{f}}=|\langle\nu_{f}|\nu_{i}(t)\rangle|^{2}\,.$ (6)

While ${\cal H}$ carries a CPV Majorana phase $\varphi$, it can be absorbed in
the redefinition of the field $\nu^{\prime}_{\beta}=e^{i\varphi}\nu_{\beta}$,
which makes the EOM real

$\displaystyle i\partial_{t}\left(\begin{matrix}\nu_{\alpha}\\\
\nu^{\prime}_{\beta}\end{matrix}\right)$ $\displaystyle={\cal
H}^{\prime}\left(t\right)\left(\begin{matrix}\nu_{\alpha}\\\
\nu^{\prime}_{\beta}\end{matrix}\right)=\left(\begin{matrix}{\cal
H}_{11}&{\cal H}_{12}\\\ {\cal H}_{12}&{\cal
H}_{22}\end{matrix}\right)\left(\begin{matrix}\nu_{\alpha}\\\
\nu^{\prime}_{\beta}\end{matrix}\right)\,.$ (7)

In the context of neutrino oscillations, this field redefinition has no
physical significance, since $\nu^{\prime}_{\beta}$ is still a flavor
eigenstate, and the oscillation probabilities for any $\nu^{\prime}$ and $\nu$
states that are related by arbitrary phases
$\nu^{\prime}_{\alpha,\beta}=e^{i\varphi_{\alpha,\beta}}\nu_{\alpha,\beta}$
are the same, since

$\displaystyle
P_{\nu^{\prime}_{i}\rightarrow\nu^{\prime}_{f}}=|\langle\nu^{\prime}_{f}|\nu^{\prime}_{i}(t)\rangle|^{2}=|\langle\nu_{f}|e^{i(\varphi_{i}-\varphi_{f})}|\nu_{i}(t)\rangle|^{2}=P_{\nu_{i}\rightarrow\nu_{f}}\,.$
(8)

For anti-neutrinos, which are the CP conjugates of $\nu$, i.e.
${\bar{\nu}}=i\sigma_{2}\nu^{*}$, the EOM is given by

$\displaystyle i\partial_{t}\left(\begin{matrix}\overline{\nu}_{\alpha}\\\
\overline{\nu}_{\beta}\end{matrix}\right)$ $\displaystyle={\cal
H}^{*}\left(t\right)\left(\begin{matrix}\overline{\nu}_{\alpha}\\\
\overline{\nu}_{\beta}\end{matrix}\right)\,,$ (9)

and thus we may again redefine our fields as
$\overline{\nu^{\prime}}_{\beta}=e^{-i\varphi}\overline{\nu}_{\beta}$ and
obtain

$\displaystyle i\partial_{t}\left(\begin{matrix}\overline{\nu}_{\alpha}\\\
\overline{\nu^{\prime}}_{\beta}\end{matrix}\right)$ $\displaystyle={\cal
H}^{\prime}\left(t\right)\left(\begin{matrix}\overline{\nu}_{\alpha}\\\
\overline{\nu^{\prime}}_{\beta}\end{matrix}\right).$ (10)

We then learn that, since ${\cal H^{\prime}}$ does not depend on $\varphi$,
anti-neutrinos propagate like neutrinos.

In the presence of the oscillating ULDM background, $\phi(t)$, the CPV phase,
$\varphi$, could become time-dependent. Then, the field redefinition
$\nu^{\prime}_{\beta}=e^{i\varphi(t)}\nu_{\beta}$ yields

$\displaystyle i\partial_{t}\left(\begin{matrix}\nu_{\alpha}\\\
\nu^{\prime}_{\beta}\end{matrix}\right)$
$\displaystyle=\left(\begin{matrix}{\cal H}_{11}&{\cal H}_{12}\\\ {\cal
H}_{12}&{\cal
H}_{22}-\dot{\varphi}\end{matrix}\right)\left(\begin{matrix}\nu_{\alpha}\\\
\nu^{\prime}_{\beta}\end{matrix}\right)\,,$ (11)

for neutrinos, and

$\displaystyle i\partial_{t}\left(\begin{matrix}\overline{\nu}_{\alpha}\\\
\overline{\nu^{\prime}}_{\beta}\end{matrix}\right)$
$\displaystyle=\left(\begin{matrix}{\cal H}_{11}&{\cal H}_{12}\\\ {\cal
H}_{12}&{\cal
H}_{22}+\dot{\varphi}\end{matrix}\right)\left(\begin{matrix}\overline{\nu}_{\alpha}\\\
\overline{\nu^{\prime}}_{\beta}\end{matrix}\right)\,,$ (12)

for anti-neutrinos. Note that the oscillation probabilities for $\nu$ and
$\nu^{\prime}$ are still identical, since equation (8) still holds even when
promoting $\varphi$ to be time-dependent. Thus, neutrinos and anti-neutrinos
would evolve differently in time, as a result of the time-dependent Majorana
phase.

While the instantaneous phase of the fields is not physically observable, its
derivative – which can be thought of as an instantaneous frequency – shifts
the energy of the system, and is thus observable. Indeed, the Hamiltonians in
Eq. (11) and Eq. (12) are very similar to those resulting from charged current
matter effects, identifying $-\dot{\varphi}$ with the matter potential $V_{C}$
[12]. Note that in our case $-\dot{\varphi}$ is a periodic function of time,
which will ultimately result in a time periodic CP asymmetry.

For concreteness, let us explicitly examine this effect in a two generations
system. Consider the following Yukawa matrix $\hat{y}$ describing the ULDM
coupling to the unperturbed mass eigenstates in two generations,

$\displaystyle\hat{y}$ $\displaystyle=iy_{I}\left(\begin{matrix}0&1\\\
1&0\end{matrix}\right),$ (13)

with real $y_{I}$. Given the Yukawa interaction in Eq. (3) is Majorana-like,
the most general form of $\hat{y}$ is $\hat{y}=a_{\mu}\sigma_{\mu}$, with
$\mu=0,1,3\,,$ $\sigma_{0}=\mathbf{1}_{2}$, $\sigma_{1,3}$ are the
corresponding symmetric Pauli matrices, and the $a_{\mu}$’s are complex
numbers. Note that within the two generation framework, at any given time the
unperturbed instantaneous PMNS matrix can be brought to be real, and so it is
straightforward to show that the form of $\hat{y}$ written above uniquely
leads to two-generation CPV. In the interaction basis, in the ultra-
relativistic limit, the Hamiltonian can be written as

$\displaystyle{\cal H}$
$\displaystyle=\frac{1}{2E}\hat{m}_{\nu}^{\dagger}\hat{m}_{\nu}\sim\frac{\Delta
m^{2}}{4E}\left(\begin{matrix}-\cos 2\theta&\sin 2\theta\\\ \sin 2\theta&\cos
2\theta\end{matrix}\right)+\frac{\phi_{0}}{2E}\left(\begin{matrix}0&-iy_{I}\Delta
m\\\ iy_{I}\Delta
m&0\end{matrix}\right)\sin\left(m_{\phi}\left(t-t_{0}\right)\right)+{\cal
O}\left(\frac{\phi_{0}^{2}}{E}\right),$ (14)

where $\theta$ is the mixing angle between the interaction eigenstates in the
unperturbed mass eigenstates, and $\Delta m=m_{2}-m_{1}$, $\Delta
m^{2}=m_{2}^{2}-m_{1}^{2}$. To leading order in the ULDM perturbation, we
identify

$\displaystyle\dot{\varphi}\approx-\frac{2m_{\phi}\phi_{0}y_{I}}{\left(m_{1}+m_{2}\right)\sin{2\theta}}\cos\left(m_{\phi}\left(t-t_{0}\right)\right)\,.$
(15)

Following the prescription presented in the next section, the leading order
difference between the survival probability of neutrinos $P_{\alpha\alpha}$
and the survival probability of anti-neutrinos $P_{\bar{\alpha}\bar{\alpha}}$
is given by

$\displaystyle\Delta P^{(1)}_{\alpha\alpha}$
$\displaystyle=\frac{2\phi_{0}\sin{4\theta}y_{I}\left(m_{1}-m_{2}\right)}{E\left(\left(\frac{\Delta
m^{2}}{2E}\right)^{2}-m_{\phi}^{2}\right)}\times$ $\displaystyle\times\Big{[}$
$\displaystyle\cos\left(m_{\phi}\left(\left(t-L-
t_{0}\right)+\frac{L}{2}\right)\right)\left(\frac{\Delta
m^{2}}{4E}\sin{\frac{\Delta
m^{2}L}{2E}}\sin{\frac{m_{\phi}L}{2}}-m_{\phi}\sin^{2}{\frac{\Delta
m^{2}L}{4E}}\cos{\frac{m_{\phi}L}{2}}\right)\Big{]}\,.$ (16)

In the following sections, we discuss the implications of this “2-$\nu$ CPV”
effect, and how it can be searched for in neutrino oscillations experiments.

##  3 Analytic approximation for small perturbations

###  3.1 Generic Prescription

The neutrino Schrödinger EOMs can be written as

$\displaystyle
i\partial_{t}\nu\left(t\right)=\frac{1}{2E}{\hat{m}}^{\dagger}\hat{m}=\frac{1}{2E}\left(m^{\dagger}m+\left(m^{\dagger}y+y^{\dagger}m\right)\phi\left(t\right)+y^{\dagger}y\phi^{2}\left(t\right)\right)\,.~{}$
(17)

The solutions of systems with time dependent Hamiltonians
$\mathcal{H}(t)=\mathcal{H}_{0}+V(t)$ can be found using the Dyson series. A
neutrino state at any final time point $t_{f}$, is evolved from an initial
neutrino state at time $t_{i}$ as
$\ket{\nu\left(t_{f}\right)}=\mathcal{U}(t_{f},t_{i})\ket{\nu\left(t_{i}\right)}$
with

$\displaystyle\mathcal{U}(t_{f},t_{i})=\sum_{n=0}^{\infty}\mathcal{U}^{(n)}(t_{f},t_{i})\equiv\sum_{n=0}^{\infty}\left(-i\right)^{n}\int\limits_{t_{f}\geq
t_{1}...\geq t_{n}\geq t_{i}}dt_{1}...dt_{n}e^{-i{\cal
H}_{0}\left(t-t_{1}\right)}V(t_{1})e^{-i{\cal
H}_{0}\left(t_{1}-t_{2}\right)}...V(t_{n})e^{-i{\cal
H}_{0}\left(t_{n}-t_{i}\right)}\,.$ (18)

The probability for a neutrino produced in a state $\nu_{\alpha}$ at time
$t_{i}$ to be measured as neutrino state $\nu_{\beta}$ at a later time $t_{f}$
is

$\displaystyle P_{\alpha\beta}\left(t_{f}\right)$
$\displaystyle=\left|\left\langle\nu_{\beta}|\mathcal{U}\left(t_{f},t_{i}\right)|\nu_{\alpha}\right\rangle\right|^{2}\equiv\left|\mathcal{U}_{\beta\alpha}\left(t_{f},t_{i}\right)\right|^{2}\,.$
(19)

In our case, we identify ${\cal H}_{0}$ with the ULDM-independent contribution
to Eq. (17), and denote $V(t)=V^{(1)}(t)+V^{(2)}\left(t\right)$ with

$\displaystyle V^{(1)}(t)$
$\displaystyle\equiv\frac{\phi_{0}}{2E}\tilde{m}\sin(m_{\phi}(t-t_{0}))\,,$
(20) $\displaystyle V^{(2)}(t)$
$\displaystyle\equiv\frac{\phi^{2}_{0}}{2E}y^{\dagger}y\sin^{2}(m_{\phi}(t-t_{0}))\,,$
(21)

where we define $\tilde{m}\equiv
m^{\dagger}_{\nu}\hat{y}+{\hat{y}}^{\dagger}m_{\nu}$ .

Then, the ULDM-independent probabilities $P^{(0)}_{\alpha\beta}$ are given by

$\displaystyle
P^{(0)}_{\alpha\beta}=\left|\mathcal{U}^{(0)}_{\beta\alpha}\left(t_{f},t_{i}\right)\right|^{2}\,,$
(22)

where

$\displaystyle\mathcal{U}^{(0)}\left(t_{f},t_{i}\right)$
$\displaystyle=e^{-i{\cal H}_{0}\left(t_{f}-t_{i}\right)}\,.$ (23)

As expected
$\mathcal{U}^{(0)}\left(t_{f},t_{i}\right)=\mathcal{U}^{(0)}\left(t_{f}-t_{i}\right)$.
From this point onwards, we assume $t_{f}-t_{i}$ corresponds to the
propagation time over a distance $L$, i.e $t_{f}-t_{i}\approx L$. Thus
$\mathcal{U}^{(0)}$ and $P^{(0)}$ are functions of $L$ alone, and do not vary
with time. For the sake of brevity, moving forward we keep the $L$ dependence
implicit, and only explicitly note time dependencies. In addition, in the
following we assume $\alpha$ and $\beta$ are interaction eigenstates.

The leading order perturbation to the transition/survival probabilities is
given by (no summation is implied on $\alpha,\beta$)

$\displaystyle P^{(1)}_{\alpha\beta}(t)$
$\displaystyle=2\text{Re}\left(\mathcal{U}_{\beta\alpha}^{(0)}{\mathcal{U}_{\beta\alpha}^{(1)}\left(t\right)}^{*}\right)\,,$
(24)

where

$\displaystyle\mathcal{U}^{(1)}\left(t\right)$
$\displaystyle=-i\int\limits_{t-L}^{t}dt_{1}e^{-i{\cal
H}_{0}\left(t-t_{1}\right)}V^{(1)}\left(t_{1}\right)e^{-i{\cal
H}_{0}\left(t_{1}-t+L\right)}=-i\int\limits_{0}^{L}dt_{1}e^{-i{\cal
H}_{0}\left(L-t_{1}\right)}V^{(1)}\left(t_{1}+t-L\right)e^{-i{\cal
H}_{0}t_{1}}\,.$ (25)

Since $V^{(1)}$ given in Eq. (20) is sinusoidal, the integral is easily solved
analytically, and we obtain

$\displaystyle P^{(1)}_{\alpha\beta}\left(t\right)$
$\displaystyle=\frac{2\phi_{0}}{E}\text{Im}\left(U_{\beta k}U^{*}_{\alpha
k}U^{*}_{\beta i}U_{\alpha j}{\tilde{m}_{ij}}^{*}e^{i\frac{\left(\Delta
m^{2}_{jk}+\Delta
m^{2}_{ik}\right)L}{4E}}{\kappa^{*}\left(t\right)_{ij}}\right)\,,$ (26)

where the $i,j,k$ indices correspond to the eigenstates of the ULDM-
independent Hamiltonian $\mathcal{H}_{0}$, and $U$ is the ULDM-independent
PMNS matrix, namely - the transition matrix between the interaction
eigenstates and the eigenstates of $\mathcal{H}_{0}$. Specifying the form of
the ULDM couplings $\hat{y}_{ij}$ in the basis in which $\mathcal{H}_{0}$ is
diagonal,

$\displaystyle\tilde{m}_{ij}=\left(m^{\dagger}_{\nu}\hat{y}+{\hat{y}}^{\dagger}m_{\nu}\right)_{ij}=m_{i}\hat{y}_{ij}+{\hat{y}_{ij}}^{*}m_{j}=\text{Re}\left(\hat{y}_{ij}\right)\left(m_{i}+m_{j}\right)+i\text{Im}\left(\hat{y}_{ij}\right)\left(m_{i}-m_{j}\right)\,.$
(27)

Lastly, $\kappa_{ij}\left(t\right)$ is complex, with phases that are only
dynamical (CPC). Its real and imaginary parts are given by

$\displaystyle\text{Re}\left(\kappa_{ij}\right)$
$\displaystyle=\frac{\sin\left(m_{\phi}\left(t-\frac{L}{2}-t_{0}\right)\right)}{\left(\frac{\Delta
m^{2}_{ij}}{2E}\right)^{2}-m_{\phi}^{2}}\left(m_{\phi}\cos{\frac{\Delta
m^{2}_{ij}L}{4E}}\sin{\frac{m_{\phi}L}{2}}-\frac{\Delta
m^{2}_{ij}}{2E}\sin{\frac{\Delta
m^{2}_{ij}L}{4E}}\cos{\frac{m_{\phi}L}{2}}\right)\,,$ (28)
$\displaystyle\text{Im}\left(\kappa_{ij}\right)$
$\displaystyle=\frac{\cos\left(m_{\phi}\left(t-\frac{L}{2}-t_{0}\right)\right)}{\left(\frac{\Delta
m^{2}_{ij}}{2E}\right)^{2}-m_{\phi}^{2}}\left(\frac{\Delta
m^{2}_{ij}}{2E}\cos{\frac{\Delta
m^{2}_{ij}L}{4E}}\sin{\frac{m_{\phi}L}{2}}-m_{\phi}\sin{\frac{\Delta
m^{2}_{ij}L}{4E}}\cos{\frac{m_{\phi}L}{2}}\right)\,.$ (29)

From the temporal dependence of $\kappa$, and thus also of the neutrino
probabilities, we see that simply adding up data points gathered over times
larger than $m_{\phi}^{-1}$ would suppress the effect due to averaging over
different neutrino arrival times $t$. We therefore require time-stamping the
events with a minimal temporal resolution corresponding to $\sim 1/m_{\phi}$.
While most previous analyses assumed the ULDM field to be constant during the
propagation time $L$, we are interested in extending these results into the
$m_{\phi}L\gtrsim 1$ region. For long baseline neutrino experiments, this is
equivalent to

$m_{\phi}\gtrsim 10^{-12}\text{ eV},$ (30)

or

$\tau_{\phi}\equiv\frac{2\pi}{m_{\phi}}\lesssim 1\text{ ms}.$ (31)

This time scale is much shorter than the typical time between events in any
neutrino experiment. In order to study how to probe these time variations, we
will use the unbinned Rayleigh periodogram, as will be explained in Sec. 4.

Note that in the limit $m_{\phi}L\rightarrow 0\,,2Em_{\phi}/\Delta
m_{ij}^{2}\rightarrow 0\,$, we have
$\text{Im}\left(\kappa_{ij}\right)/\text{Re}\left(\kappa_{ij}\right)\rightarrow
0$ ($\text{Re}\left(\kappa_{ij}\right)$ remains finite as
$V\left(t\right)\propto\sin{m_{\phi}\left(t-t_{0}\right)}$). Therefore,
previous analyses that assumed the ULDM potential to be constant within the
propagation time, have effectively only taken into account effects that are
associated with the real part of $\kappa$, which in that limit reduces to

$\displaystyle\text{Re}\left(\kappa_{ij}\right)$
$\displaystyle\approx-\frac{\sin\left(m_{\phi}\left(t-t_{0}\right)\right)}{\frac{\Delta
m^{2}_{ij}}{2E}}\sin{\frac{\Delta m^{2}_{ij}L}{4E}}\,.$ (32)

In this limit, for ULDM couplings that are off-diagonal in the unperturbed
mass basis, the result can be interpreted as following from the first order
correction to the diagonalizing matrix of the neutrino Hamiltonian in the
interaction basis. Namely,

$\displaystyle U_{\alpha k}\rightarrow\tilde{U}_{\alpha
k}\left(t\right)=\left(UT^{\dagger}\left(t\right)\right)_{\alpha k}\,,$ (33)

where

$\displaystyle
T\left(t\right)_{ij}=\delta_{ij}+\left(1-\delta_{ij}\right)\frac{\phi_{0}\tilde{m}_{ij}}{\Delta
m^{2}_{ij}}\sin\left(m_{\phi}\left(t-t_{0}\right)\right)\,.$ (34)

Note that in this limit, the new diagonalizing matrix $\tilde{U}$ is treated
as constant within propagation, but changes its value for each measurement at
time $t$.

To understand what happens when going beyond the constant potential limit, let
us rewrite the EOM in the interaction basis as

$\displaystyle
i\partial_{t}\nu^{I}\left(t\right)=UT^{\dagger}(t)\mathcal{H}_{D}T(t)U^{\dagger}\nu^{I}\left(t\right)\approx
UT^{\dagger}(t)\mathcal{H}_{0}T(t)U^{\dagger}\nu^{I}\left(t\right)\,,$ (35)

where $\nu^{I}$ is given in the interaction basis, and $\mathcal{H}_{D}$ is
the instantaneously diagonal Hamiltonian, which to linear order in the ULDM
perturbation is just $\mathcal{H}_{0}$ . To solve the EOM, one may define

$\displaystyle\tilde{\nu}\left(t\right)={\tilde{U}\left(t\right)}^{\dagger}\nu^{I}\left(t\right)\,,$
(36)

and obtain

$\displaystyle
i\partial_{t}\tilde{\nu}\left(t\right)=\left(\mathcal{H}_{0}-i{\tilde{U}\left(t\right)}^{\dagger}\partial_{t}\tilde{U}\left(t\right)\right)\tilde{\nu}\left(t\right)=\left(\mathcal{H}_{0}-iT\left(t\right)\partial_{t}{T\left(t\right)}^{\dagger}\right)\tilde{\nu}\left(t\right)\,.$
(37)

As the EOM is not diagonal due to the contributions coming from the time-
derivative of the Hamiltonian, one should in principle diagonlize the new
effective Hamiltonian by another time-dependent matrix, the derivative of
which will again contribute to the Hamiltonian which should again be
diagonalized, and so on. Keeping only the linear terms in the ULDM
perturbation, this procedure will yield an infinite series in
$\partial^{(n)}_{t}V^{(1)}_{ij}/\left(\Delta m^{2}_{ij}\right)^{n}$, where
$i\neq j$ and $V^{(1)}$ is defined in Eq. (20). The result of this series is
the effective diagonalizing matrix

$\displaystyle\tilde{U}\left(t\right)=U\tilde{T}^{\dagger}\left(t\right)\,,$
(38)

with

$\displaystyle\tilde{T}_{ij}\left(t\right)=\delta_{ij}+\left(1-\delta_{ij}\right)\frac{\phi_{0}\tilde{m}_{ij}}{\Delta
m^{2}_{ij}}\sum_{n=0}^{\infty}i^{n}\left(\frac{2E}{\Delta
m_{ij}^{2}}\right)^{n}\partial^{(n)}_{t}\sin\left(m_{\phi}\left(t-t_{0}\right)\right)\,.$
(39)

Since our potential is harmonic, we can easily sum all the contributions up to
infinity, and obtain

$\displaystyle\tilde{T}_{ij}\left(t\right)=\delta_{ij}+\left(1-\delta_{ij}\right)\frac{\phi_{0}\tilde{m}_{ij}}{\Delta
m^{2}_{ij}}\frac{1}{1-\left(\frac{2Em_{\phi}}{\Delta
m_{ij}^{2}}\right)^{2}}\left(\sin\left(m_{\phi}\left(t-t_{0}\right)\right)+i\frac{2Em_{\phi}}{\Delta
m_{ij}^{2}}\cos\left(m_{\phi}\left(t-t_{0}\right)\right)\right)\,.$ (40)

Therefore, the neutrino flavor probabilities are given by taking the zeroth
and first order terms in $\phi_{0}\tilde{m}_{ij}/\Delta m^{2}_{ij}$ of

$\displaystyle P_{\alpha\beta}\left(t\right)\approx\tilde{U}_{\beta
k}\left(t\right)\tilde{U}^{*}_{\alpha k}\left(t-L\right)\tilde{U}^{*}_{\beta
i}\left(t\right)\tilde{U}_{\alpha i}\left(t-L\right)e^{-i\frac{\Delta
m^{2}_{ki}L}{2E}}\,.$ (41)

In the next subsections, we describe the full ULDM effects in three
generations, first considering CPC observables in subsection 3.2, and then CPV
phenomena in subsection 3.3. Note, however, that the expressions above are
completely generic for any number of generations, provided the indices are
appropriately summed. The probabilities for the two generations case are
presented in Appendix A.

###  3.2 CPC effects

In the standard case with three neutrino generations, the neutrino
probabilities in three generations are given by

$\displaystyle P^{(0)}_{\alpha\beta}$ $\displaystyle=U_{\beta k}U^{*}_{\alpha
k}U^{*}_{\beta i}U_{\alpha i}e^{-i\frac{\Delta m^{2}_{ki}L}{2E}}\,.$ (42)

The CPC part of this probability is given by averaging over neutrinos and
anti-neutrinos

$\displaystyle\Sigma P^{(0)}_{\alpha\beta}$
$\displaystyle\equiv\frac{1}{2}\left(P^{(0)}_{\alpha\beta}+P^{(0)}_{\bar{\alpha}\bar{\beta}}\right)=\text{Re}\left(U_{\beta
k}U^{*}_{\alpha k}U^{*}_{\beta i}U_{\alpha i}\right)\cos\left(\frac{\Delta
m^{2}_{ki}L}{2E}\right)\,.$ (43)

We now discuss the first order corrections from ULDM to CPC neutrino
probabilities. From Eq. (26), we find CPC contributions to be associated with
the real parts of the complex matrices $U$ and $\hat{y}$. Let us denote these
contributions by $\Sigma P^{(1)}_{\alpha\beta}$:

$\displaystyle\Sigma
P^{(1)}_{\alpha\beta}\equiv\frac{1}{2}\left(P^{(1)}_{\alpha\beta}+P^{(1)}_{\bar{\alpha}\bar{\beta}}\right)$
$\displaystyle=\frac{2\phi_{0}}{E}\text{Re}\left(U_{\beta k}U^{*}_{\alpha
k}U^{*}_{\beta i}U_{\alpha
j}{\tilde{m}_{ij}}^{*}\right)\text{Im}\left(e^{i\frac{\left(\Delta
m^{2}_{jk}+\Delta
m^{2}_{ik}\right)L}{4E}}{\kappa^{*}\left(t\right)_{ij}}\right)\,.$ (44)

In three generations, for diagonal couplings $\hat{y}_{ii}$, both
$\tilde{m}_{ii}$ and $\kappa_{ii}$ are real, and thus for a specific $i$,
specific $j\neq k\neq i$ and any $\alpha,\beta$

$\displaystyle\Sigma P^{(1)}_{\alpha\beta}\left(t\right)$
$\displaystyle=-\sin\left(m_{\phi}\left(t-\frac{L}{2}-t_{0}\right)\right)\sin\left(\frac{m_{\phi}L}{2}\right)\frac{2\tilde{m}_{ii}\phi_{0}}{Em_{\phi}}\times$
$\displaystyle\times\left(\text{Re}\left(U_{\beta k}U^{*}_{\alpha
k}U^{*}_{\beta i}U_{\alpha i}\right)\sin{\left(\frac{\Delta
m^{2}_{ik}L}{2E}\right)}+\text{Re}\left(U_{\beta j}U^{*}_{\alpha
j}U^{*}_{\beta i}U_{\alpha i}\right)\sin{\left(\frac{\Delta
m^{2}_{ij}L}{2E}\right)}\right)\,.$ (45)

This result is nothing but an effective, time-dependent, shift of the neutrino
mass differences $\Delta m^{2}_{ij}$

$\displaystyle\Delta m^{2}_{ij}\rightarrow\Delta
m^{2}_{ij}+\frac{2m_{i}\phi_{0}\hat{y}_{ii}}{L}\int^{t}_{t-L}\sin\left(m_{\phi}\left(t_{1}-t_{0}\right)\right)dt_{1}\,.$
(46)

Effects of this sort have already been quite thoroughly discussed in previous
works such as [1, 2, 7, 6].

For off-diagonal couplings $\hat{y}_{ij}$, with $i\neq j\neq k$ transition
probabilities ($\alpha\neq\beta$) are modified by

$\displaystyle\Sigma
P^{(1)}_{\alpha\beta}=\frac{4\phi_{0}}{E}\text{Im}\left(\kappa^{*}\left(t\right)_{ij}\right)\,$
$\displaystyle\text{Re}\left(U_{\beta k}U^{*}_{\alpha k}\left(U^{*}_{\beta
i}U_{\alpha j}{\tilde{m}_{ij}}^{*}-U^{*}_{\beta j}U_{\alpha
i}{\tilde{m}_{ij}}\right)\right)\sin{\frac{\Delta
m^{2}_{ik}L}{4E}}\sin{\frac{\Delta m^{2}_{kj}L}{4E}}+$
$\displaystyle\frac{4\phi_{0}}{E}\text{Re}\left(\kappa^{*}\left(t\right)_{ij}\right)\Bigg{[}$
$\displaystyle\text{Re}\left(U_{\beta j}U^{*}_{\alpha j}\left(U^{*}_{\beta
i}U_{\alpha j}{\tilde{m}_{ij}}^{*}+U^{*}_{\beta j}U_{\alpha
i}{\tilde{m}_{ij}}\right)\right)\sin{\frac{\Delta
m^{2}_{kj}L}{4E}}\cos{\frac{\Delta m^{2}_{ki}L}{4E}}$ $\displaystyle+$
$\displaystyle\,\text{Re}\left(U_{\beta i}U^{*}_{\alpha i}\left(U^{*}_{\beta
i}U_{\alpha j}{\tilde{m}_{ij}}^{*}+U^{*}_{\beta j}U_{\alpha
i}{\tilde{m}_{ij}}\right)\right)\sin{\frac{\Delta
m^{2}_{ki}L}{4E}}\cos{\frac{\Delta m^{2}_{kj}L}{4E}}\Bigg{]}\,,$ (47)

and survival probabilities ($\alpha=\beta$) by

$\displaystyle\Sigma P^{(1)}_{\alpha\alpha}$
$\displaystyle=\frac{4\phi_{0}}{E}\text{Re}\left(\kappa\left(t\right)_{ij}\right)\text{Re}\left(U^{*}_{\alpha
i}U_{\alpha j}{\tilde{m}_{ij}}^{*}\right)\times$
$\displaystyle\times\Bigg{[}\left(2|U_{\alpha
j}|^{2}-1\right)\sin{\frac{\Delta m^{2}_{kj}L}{4E}}\cos{\frac{\Delta
m^{2}_{ki}L}{4E}}+\left(2|U_{\alpha i}|^{2}-1\right)\sin{\frac{\Delta
m^{2}_{ki}L}{4E}}\cos{\frac{\Delta m^{2}_{kj}L}{4E}}\Bigg{]}\,.$ (48)

As mentioned, previous analyses considered the constant ULDM limit, in which
$\kappa_{ij}$ is strictly real, and given by Eq. (32). In this case, one may
interpret the first-order corrections as simple modifications to the PMNS
matrix $U\rightarrow\tilde{U}$ as in Eq. (33) and Eq. (34). However, moving
away from this limit, we may no longer use the constant $\tilde{U}$
description. This leads to the different energy dependence of the result
associated with $\text{Im}\left(\kappa\right)$ and with the non-zeroth order
terms in $2Em_{\phi}/\left(\Delta m^{2}_{ij}\right)$ appearing in
$\text{Re}\left(\kappa\right)$ .

###  3.3 CPV effects

CPV is realized in vacuum in the three neutrino picture. The difference
between the neutrino and antineutrino transition probability is given by

$\displaystyle\Delta P_{\alpha\beta}^{\text{standard}}$
$\displaystyle=4\sum_{i>j}\text{Im}\left(U_{\alpha i}^{*}U_{\beta i}U_{\alpha
j}U_{\beta j}^{*}\right)\sin{\frac{\Delta
m_{ij}^{2}L}{2E}}=16J\sin{\frac{\Delta m_{31}^{2}L}{4E}}\sin{\frac{\Delta
m_{21}^{2}L}{4E}}\sin{\frac{\Delta
m_{32}^{2}L}{4E}}\sum_{\gamma}\varepsilon_{\alpha\beta\gamma}\,,$ (49)

where the $\varepsilon_{\alpha\beta\gamma}$ is the three-dimensional Levi-
Civita tensor and $J$ is the Jarlskog invariant:

$J=\cos\theta_{13}\sin 2\theta_{12}\sin 2\theta_{13}\sin
2\theta_{23}\sin\delta_{\text{CP}}.$ (50)

Notice that CPV vanishes if we set any $\Delta m^{2}_{ij}$ to be zero, which
will result in an effective two-neutrino picture. We therefore name this
source of CPV the “$3-\nu$ CPV”.

Let us now discuss CPV in the presence of ULDM interactions with neutrinos.
From Eq. (26), we find the CPV part of neutrino probabilities, i.e. the
difference between neutrinos and anti-neutrinos probabilities $\Delta
P^{(1)}_{\alpha\beta}$ to be

$\displaystyle\Delta P^{(1)}_{\alpha\beta}\equiv
P^{(1)}_{\alpha\beta}-P^{(1)}_{\bar{\alpha}\bar{\beta}}$
$\displaystyle=-\frac{4\phi_{0}}{E}\sum_{i,j,k}\text{Im}\left(U^{*}_{\beta
k}U_{\alpha k}U_{\beta i}U^{*}_{\alpha
j}{\tilde{m}_{ij}}\right)\text{Re}\left(e^{i\frac{\left(\Delta
m^{2}_{jk}+\Delta m^{2}_{ik}\right)L}{4E}}{\kappa_{ij}}^{*}\right)\,.$ (51)

Previous works have assumed the neutrino Hamiltonian to be constant throughout
propagation, and thus the ULDM effects were thought of as slow temporal
modulations of the standard neutrino oscillation parameters. In this limit, as
mentioned, $\kappa_{ij}$ is strictly real, and CPV can be interpreted as a
slow variation of the standard CPV result [6]. Explicitly, in three
generations, the contribution of the real part of $\kappa$ for $i\neq j\neq k$
and $\alpha\neq\beta$, for a specific $\hat{y}_{ij}=\hat{y}_{ji}\neq 0$ is

$\displaystyle\Delta P^{(1)}_{\alpha\beta}$
$\displaystyle=-8\frac{\phi_{0}}{E}\text{Re}\left(\kappa_{ij}\right)\sin{\frac{\Delta
m^{2}_{jk}L}{4E}}\sin{\frac{\Delta
m^{2}_{ik}L}{4E}}\text{Im}\left(U^{*}_{\beta k}U_{\alpha k}\left(U_{\beta
i}U^{*}_{\alpha j}\tilde{m}_{ij}+U_{\beta j}U^{*}_{\alpha
i}\tilde{m}^{*}_{ij}\right)\right)\,,$ (52)

which yields in the $m_{\phi}L\rightarrow 0,2Em_{\phi}/\Delta
m_{ij}^{2}\rightarrow 0\,,V\left(t\right)\propto\sin{m_{\phi}t}$ limit

$\displaystyle\Delta P^{(1)}_{\alpha\beta}\approx$
$\displaystyle-16\frac{\phi_{0}}{\Delta
m^{2}_{ij}}\sin\left(m_{\phi}\left(t-\frac{L}{2}-t_{0}\right)\right)\cos{\frac{m_{\phi}L}{2}}\times$
$\displaystyle\times\sin{\frac{\Delta m^{2}_{ij}L}{4E}}\sin{\frac{\Delta
m^{2}_{jk}L}{4E}}\sin{\frac{\Delta
m^{2}_{ik}L}{4E}}\text{Im}\left(U^{*}_{\beta k}U_{\alpha k}\left(U_{\beta
i}U^{*}_{\alpha j}\tilde{m}_{ij}+U_{\beta j}U^{*}_{\alpha
i}\tilde{m}^{*}_{ij}\right)\right)\,.$ (53)

Since we keep only first-order terms in the ULDM couplings, which here we
assume to be off-diagonal in the ULDM-independent mass basis, one may
interpret this result as the effect of the first-order correction to the
standard Jarlskog invariant.

For diagonal couplings $\hat{y}_{ii}$, both $\tilde{m}_{ii}$ and $\kappa_{ii}$
are real, and thus

$\displaystyle\Delta P^{(1)}_{\alpha\beta}\left(t\right)$
$\displaystyle=-\frac{4\tilde{m}_{ii}\phi_{0}}{Em_{\phi}}\sin\left(m_{\phi}\left(t-\frac{L}{2}-t_{0}\right)\right)\sin\left(\frac{m_{\phi}L}{2}\right)\times$
$\displaystyle\times\text{Im}\left(U_{\beta j}U^{*}_{\alpha j}U^{*}_{\beta
i}U_{\alpha i}\right)\sin{\left(\frac{\left(\Delta m^{2}_{ik}+\Delta
m^{2}_{ij}\right)L}{4E}\right)}\sin{\left(\frac{\Delta
m^{2}_{kj}L}{4E}\right)}\,.$ (54)

This CPV term is of course proportional to $J$, as only the real part of
diagonal couplings affects the neutrino Hamiltonian. This term is simply a
correction to the mass differences of the neutrinos coming from the ULDM
coupling, as shown in Eq. (46).

For off-diagonal couplings, stepping away from the $m_{\phi}\rightarrow 0$
limit, we should also consider the contribution of the imaginary part of
$\kappa$

$\displaystyle\Delta P^{(1)}_{\alpha\beta}=$ $\displaystyle
4\text{Im}\left(\kappa_{ij}\right)\sum_{i,j,k}\text{Im}\left(U^{*}_{\beta
k}U_{\alpha k}U_{\beta i}U^{*}_{\alpha
j}{\tilde{m}_{ij}}\right)\sin\left(\frac{\left(\Delta m^{2}_{jk}+\Delta
m^{2}_{ik}\right)L}{4E}\right)\,.$ (55)

For $i\neq j\neq k$ and $\alpha\neq\beta$, for a specific
$\hat{y}_{ij}=\hat{y}_{ji}\neq 0$, it becomes

$\displaystyle\Delta
P^{(1)}_{\alpha\beta}=8\text{Im}\left(\kappa_{ij}\right)\Big{[}$
$\displaystyle\text{Im}\left(U^{*}_{\beta i}U_{\alpha i}\left(U_{\beta
j}U^{*}_{\alpha i}{\tilde{m}^{*}_{ij}}-U_{\beta i}U^{*}_{\alpha
j}{\tilde{m}_{ij}}\right)\right)\sin{\frac{\Delta
m^{2}_{ik}L}{4E}}\cos{\frac{\Delta m^{2}_{jk}L}{4E}}$ $\displaystyle+$
$\displaystyle\text{Im}\left(U^{*}_{\beta j}U_{\alpha j}\left(U_{\beta
j}U^{*}_{\alpha i}{\tilde{m}^{*}_{ij}}-U_{\beta i}U^{*}_{\alpha
j}{\tilde{m}_{ij}}\right)\right)\cos{\frac{\Delta
m^{2}_{ik}L}{4E}}\sin{\frac{\Delta m^{2}_{jk}L}{4E}}\Big{]}\,.$ (56)

This expression is fundamentally different from the “3-$\nu$ CPV” effect, as
it “contains” the 2-$\nu$ CPV effect we discussed in Section 2 (see Eq. (16)).
This can be seen in three distinct ways, all pointing to cases in which there
would be no CPV in the constant ULDM potential limit.

First, note that Eq. (3.3) vanishes for $U^{*}_{\beta k}U_{\alpha k}=0$ ,
namely, when the projection of the mass eigenstate not coupled to the ULDM on
either the outgoing or incoming interaction eigenstates vanishes. This case
corresponds to an effective $J=0$. However, taking this limit for Eq. (56) ,
we obtain

$\displaystyle\Delta P^{(1)}_{\alpha\beta}=$ $\displaystyle
8\text{Im}\left(\kappa_{ij}\right)\text{Im}\left(U^{*}_{\beta i}U_{\alpha
i}\left(U_{\beta j}U^{*}_{\alpha i}{\tilde{m}^{*}_{ij}}-U_{\beta
i}U^{*}_{\alpha j}{\tilde{m}_{ij}}\right)\right)\sin{\frac{\Delta
m^{2}_{ij}L}{4E}}\,.$ (57)

This expression does not vanish generically. Specifically, if say only
$U_{\alpha k}=0$ , then unitarity implies $|U_{\alpha i}|^{2}+|U_{\alpha
j}|^{2}=1$ , and we may denote $|U_{\alpha
i}|^{2}\equiv\sin^{2}{\theta_{\alpha}}$ and obtain

$\displaystyle\Delta P^{(1)}_{\alpha\beta}=$ $\displaystyle
16\cos{2\theta_{\alpha}}\text{Im}\left(\kappa_{ij}\right)\text{Im}\left(U_{\beta
i}U^{*}_{\beta j}{\tilde{m}_{ij}}\right)\sin{\frac{\Delta
m^{2}_{ij}L}{4E}}\,.$ (58)

This effect is indeed the result of the two-generations-like effect, and thus
does not vanish.

Second, as explained above, we expect the standard CPV effect in three
generations to vanish when any of the mass differences $\Delta m^{2}_{ij}$
vanishes. Indeed, the standard CPV in Eq. (49) vanishes at this limit, as well
as the constant-ULDM CPV in Eq. (3.3) when considering a specific off-diagonal
$\hat{y}_{ij}$ , and setting $\Delta m^{2}_{jk}=0$, for $k\neq i\neq j$.
However, at this limit, Eq. (56) yields a non-zero result, and reads

$\displaystyle\Delta P^{(1)}_{\alpha\beta}=$ $\displaystyle
8\text{Im}\left(\kappa_{ij}\right)\text{Im}\left(U^{*}_{\beta i}U_{\alpha
i}\left(U_{\beta j}U^{*}_{\alpha i}{\tilde{m}^{*}_{ij}}-U_{\beta
i}U^{*}_{\alpha j}{\tilde{m}_{ij}}\right)\right)\sin{\frac{\Delta
m^{2}_{ij}L}{4E}}\,,$ (59)

exactly as in the $U^{*}_{\beta k}U_{\alpha k}=0$ case (however here we do not
assume a specific structure for $U$). Given the known hierarchy in the
neutrino mass differences, $\Delta m^{2}_{12}\ll\Delta m^{2}_{31},\Delta
m^{2}_{32}$, we then expect that the new 2-$\nu$ CPV effect resulting from
$\hat{y}_{31}$ or $\hat{y}_{32}$ would dominate over the constant ULDM CPV
roughly for $m_{\phi}\gtrsim\Delta m^{2}_{12}/\left(2E\right)$. Note that if
the neutrinos coupled to the ULDM are degenerate in the unperturbed system,
i.e. $\hat{y}_{ij}\neq 0$ for $\Delta m^{2}_{ij}L/\left(4E\right)\rightarrow
0$, the constant ULDM limit $m_{\phi}L\rightarrow 0$ yields a non-zero CPV, as
$\text{Re}\left(\kappa_{ij}\right)\rightarrow\sin\left(m_{\phi}\left(t-t_{0}\right)\right)L/2$.
This corresponds to a leading order effective mass difference induced by the
ULDM potential. In this limit, we obtain
$\text{Im}\left(\kappa_{ij}\right)\propto L^{3}m_{\phi}\Delta
m^{2}_{ij}/\left(4E\right)$ , and thus the 2-$\nu$ CPV effect is expected to
be sub-dominant to the constant ULDM effect.

Finally, note that standard CPV does not affect survival probabilities.
Accordingly, the constant-ULDM CPV in Eq. (3.3) is also zero for
$\alpha=\beta$. The new CPV effect we found in Eq. (55) for $\alpha=\beta$
does not vanish generically, and is given for a specific
$\hat{y}_{ij}=\hat{y}_{ji}\neq 0$ with $i\neq j\neq k$ by

$\displaystyle\Delta P^{(1)}_{\alpha\alpha}=$
$\displaystyle-8\text{Im}\left(\kappa_{ij}\right)\text{Im}\left(\tilde{m}_{ij}U_{\alpha
i}U^{*}_{\alpha j}\right)\times$ $\displaystyle\times$
$\displaystyle\left(\left(2|U_{\alpha j}|^{2}+1\right)\sin\left(\frac{\Delta
m^{2}_{jk}L}{4E}\right)\cos\left(\frac{\Delta
m^{2}_{ik}L}{4E}\right)+\left(2|U_{\alpha
i}|^{2}+1\right)\cos\left(\frac{\Delta
m^{2}_{jk}L}{4E}\right)\sin\left(\frac{\Delta
m^{2}_{ik}L}{4E}\right)\right)\,.$ (60)

Therefore, even experiments only measuring survival probabilities would be
sensitive to CPV ULDM-neutrino couplings. This is a completely new prediction
resulting from our analysis, considering the time variations of the
Hamiltonian within the neutrino propagation time.

##  4 Experimental implications - Rayleigh Periodogram

In Section 3, we introduced analytical expressions for the first-order
modifications to neutrino probabilities resulting from ULDM-neutrino
interactions. We are interested in devising a method for detecting these
modifications in neutrino oscillations experiments. As the effects we found
are time dependent, characterized by simple oscillatory functions with periods
$\tau_{\phi}=2\pi/m_{\phi}$, we are essentially interested in measuring the
spectral power of the neutrino probabilities at an angular frequency
corresponding to $m_{\phi}$.

Consider the following time scales relevant for neutrino oscillation
experiments:

  1. 1.

$\tau_{e}$: The running time of the experiment.

  2. 2.

$\tau_{s}=\frac{\tau_{e}}{N_{\nu}}$: The average spacing between events, where
$N_{\nu}$ is the total number of measured events.

  3. 3.

$\tau_{r}$: The resolution of the event timing. While in practice the
effective resolution could be set by different sources of temporal uncertainty
(e.g. beam spread), we will refer to it as the “clock” resolution.

We notice the following hierarchy between the time scales:

$\tau_{r}\ll\tau_{s}\ll\tau_{e}.$ (61)

Since the frequencies we are interested in are determined by $\tau_{\phi}$,
one would naively assume that the experimental sensitivity to ULDM masses
$m_{\phi}\gtrsim 2/\tau_{s}$ would be suppressed, as the statistical
uncertainty on the neutrino probability binned over the corresponding
$\tau_{\phi}$ would be quite large, due to a small number of events expected
to occur within that period. However, as we will now show, the true limiting
factor, in terms of the experimental sensitivity, is rather the effective
clock resolution $\tau_{r}$, which is much smaller than $\tau_{s}$.

To understand this, first recall that the events are inherently binned over
times $\tau_{r}$, while $\tau_{s}$ is not associated with some instrumental
resolution, and binning over it is completely artificial, and is done post-
measurement. Let us then use the full data, given as the number of neutrinos
of a certain flavor detected at a sequence of $N_{t}=\tau_{e}/\tau_{r}$ clock
ticks, each with short duration $\tau_{r}$. Since $\tau_{r}\ll\tau_{s}$, we
ignore the possibility that two neutrinos might be detected during $\tau_{r}$,
and define the following power spectrum:

$\displaystyle
z\left(f\right)=\frac{2}{N_{\nu}}\left(\left[\sum_{n=1}^{N_{t}}h\left(t_{n}\right)\cos\left(2\pi
ft_{n}\right)\right]^{2}+\left[\sum_{n=1}^{N_{t}}h\left(t_{n}\right)\sin\left(2\pi
ft_{n}\right)\right]^{2}\right)\,,$ (62)

with

$h\left(t_{i}\right)=\left\\{\begin{matrix}1&\text{event detected,}\\\
0&\text{no event detected.}\end{matrix}\right.$ (63)

This definition gives us the well known Rayleigh power spectrum

$z\left(f\right)=\frac{2}{N_{\nu}}\left(\left[\sum_{n=1}^{N_{\nu}}\cos\left(2\pi
ft_{n}\right)\right]^{2}+\left[\sum_{n=1}^{N_{\nu}}\sin\left(2\pi
ft_{n}\right)\right]^{2}\right).$ (64)

Notice that the summation now is over the neutrino events instead of over time
bins. For a list of events detected in a time series $\\{t_{n}\\}$, we
calculate the power spectrum $z$, which is an analytic function of the probed
frequencies $f$. For $\\{t_{n}\\}$ uniformly random distributed between
$0<t_{n}<\tau_{e}$, the sum of sines and sum of cosines in Eq. (64) are each
normally distributed with mean value of zero and standard deviation of
$\sqrt{N_{\nu}/2}$. Thus $z(f)$, being their appropriately normalized sum of
squares, is Chi-squared distributed with two degrees of freedom, which is
simply an exponential distribution of $z$. This means, under the background
only hypothesis, that if we consider a specific frequency $f^{\prime}$ and
calculate its power $z\left(f^{\prime}\right)$, its probability to be higher
than some value $Z$ is independent of the frequency $f^{\prime}$ (assuming
$f^{\prime}\gtrsim 1/\tau_{e}$ and no other spectral noise sources), and given
by

$p\left(z>Z\right)=\int\limits_{Z}^{\infty}\frac{1}{2}e^{-z/2}dz=e^{-Z/2}.$
(65)

Let us estimate the signal in case of a temporal modulation of the neutrino
probability. The probability of detecting a neutrino of some flavor at time
$t_{i}$ in an energy bin $[E_{j,\min},E_{j,\max}]$ is given by

$P\left(h\left(t_{i},E_{j}\right)=1\right)=\int\limits_{E_{j,\min}}^{E_{j,\max}}\int\limits_{t_{i}-\tau_{r}/2}^{t_{i}+\tau_{r}/2}dtdEF\left(E\right)P\left(E,t\right),$
(66)

where $F$ is the unoscillated neutrino flux, i.e. the flux that would have
reached the detector had all the produced neutrinos oscillated into the flavor
in question

$\displaystyle F\left(E\right)$
$\displaystyle=\frac{dN_{\text{unoscillated}}}{dEdt},$ (67)

and $P\left(E,t\right)$ is the oscillation probability given by

$\displaystyle P\left(E,t\right)$
$\displaystyle=P^{\left(0\right)}\left(E\right)+P^{\left(1\right)}\left(E,t\right)\,,$
(68)

with $P^{(0)}$, $P^{(1)}$ given in Eqs. (22) and (24). The total number of
events collected from different energy bins over the entire duration of the
experiment $\tau_{e}$, assuming the ULDM modification to the probability is
very small, is then

$\displaystyle N_{\nu}$
$\displaystyle\approx\int\limits_{E_{\min}}^{E_{\max}}\int\limits_{0}^{\tau_{e}}dtdEF\left(E\right)P^{\left(0\right)}\left(E\right)=\tau_{e}\int\limits_{E_{\min}}^{E_{\max}}dEF\left(E\right)P^{\left(0\right)}\left(E\right)\,.$
(69)

Similarly, one may define the corresponding number of events associated with
the small ULDM probability

$\displaystyle N_{\nu,s}^{\left(1\right)}$
$\displaystyle=\tau_{e}\int\limits_{E_{\min}}^{E_{\max}}dEF\left(E\right)P_{s}^{\left(1\right)}\left(E\right)\,,$
(70) $\displaystyle N_{\nu,c}^{\left(1\right)}$
$\displaystyle=\tau_{e}\int\limits_{E_{\min}}^{E_{\max}}dEF\left(E\right)P_{c}^{\left(1\right)}\left(E\right)\,,$
(71)

where $P_{s}$ and $P_{c}$ can be read off from

$P^{\left(1\right)}\left(E,t\right)=P_{s}^{\left(1\right)}(E)\sin\left(m_{\phi}t\right)+P_{c}^{\left(1\right)}(E)\cos\left(m_{\phi}t\right)\,.$
(72)

To identify $P_{s}^{\left(1\right)}(E)$ and $P_{c}^{\left(1\right)}(E)$, note
that the imaginary and real components of $\kappa_{ij}$ determine the time
dependence of $P^{(1)}$, as can be seen in Eqs. (28) and (29). Since
$\text{Re}(\kappa)$ and $\text{Im}(\kappa)$ are orthogonal in time, but have
the same frequency, $P^{(1)}$ can be written as

$\displaystyle
P^{(1)}=a\sin(m_{\phi}(t-L/2-t_{0}))+b\cos(m_{\phi}(t-L/2-t_{0}))\,,$ (73)

and thus

$\displaystyle P_{s}^{\left(1\right)}$
$\displaystyle=a\cos\left(m_{\phi}\left(L/2+t_{0}\right)\right)+b\sin\left(m_{\phi}\left(L/2+t_{0}\right)\right)\,,$
(74) $\displaystyle P_{c}^{\left(1\right)}$
$\displaystyle=-a\sin\left(m_{\phi}\left(L/2+t_{0}\right)\right)+b\cos\left(m_{\phi}\left(L/2+t_{0}\right)\right)\,.$
(75)

We may now calculate the expected value of $z(f)$ at
$f=m_{\phi}/\left(2\pi\right)$ in the presence of a ULDM, with mass
$m_{\phi}$, interacting with neutrinos111 We do not specify the range and
resolution of frequencies at which a peak in the Rayleigh spectrum is searched
for, as there is no clear recipe for determining them. One can scan
frequencies as small as $\sim\tau_{e}^{-1}$, and as large as $\tau_{r}^{-1}$,
since the effect is suppressed for frequencies higher than that. Typically,
the resolution of the scan would be approximately $\sim\tau_{e}^{-1}$, as this
is the width of the features in $z(f)$. However, the authors of [13], for
example, chose to oversample with a higher resolution. This comes at the cost
of a more significant look-elsewhere effect [14], which reduces the confidence
level of the signal. Since this effect is mainly determined by the number of
frequencies one scans, a number which depends neither on the experiment nor on
the ULDM parameters, we did not include it in the analyses in this paper.. One
can show, see Appendix B, that for $m_{\phi}\lesssim 2/\tau_{r}$

$\displaystyle\left\langle\left[\sum_{i=1}^{N_{\nu}}\sin\left(2\pi
t_{i}f\right)\right]^{2}\right\rangle$
$\displaystyle=\frac{N_{\nu}}{2}+\frac{\left(N_{s}^{\left(1\right)}\right)^{2}}{4}\,,$
(76) $\displaystyle\left\langle\left[\sum_{i=1}^{N_{\nu}}\cos\left(2\pi
t_{i}f\right)\right]^{2}\right\rangle$
$\displaystyle=\frac{N_{\nu}}{2}+\frac{\left(N_{c}^{\left(1\right)}\right)^{2}}{4}\,,$
(77)

and thus

$\displaystyle\left\langle z\left(f\right)\right\rangle$
$\displaystyle=\frac{2}{N_{\nu}}\left[\left\langle\left[\sum_{i=1}^{N_{\nu}}\sin\left(2\pi
t_{i}f\right)\right]^{2}\right\rangle+\left\langle\left[\sum_{i=1}^{N_{\nu}}\cos\left(2\pi
t_{i}f\right)\right]^{2}\right\rangle\right]=2+\frac{\left(N_{s}^{\left(1\right)}\right)^{2}+\left(N_{c}^{\left(1\right)}\right)^{2}}{2N_{\nu}}\,.$
(78)

Note a few key points regarding the above result:

  * •

The expected sensitivity to a ULDM coupling $\hat{y}$ improves with the
square-root of the (unoscillated) neutrino flux of the experiment $\sqrt{F}$,
and decreases as the square-root of the expected vacuum oscillation
probability $\sqrt{P^{(0)}}$.

  * •

Assuming the ULDM perturbation is coherent throughout the experimental time ,
the sensitivity to $\hat{y}$ improves as $\sqrt{\tau_{e}}$. For a coherence
time $\tau_{\text{c}}$ satisfying $L,\tau_{r}<\tau_{\text{c}}<\tau_{e}$, the
sensitivity to $\hat{y}$ would follow
$\sqrt{\tau_{e}}\left(\frac{\tau_{c}}{\tau_{e}}\right)^{\frac{1}{4}}$. While
our results were derived assuming $\tau_{c}\gtrsim\tau_{r},L$, it is possible
to derive the appropriate expressions for $\tau_{c}<L,\tau_{r}$ by
statistically averaging over the ULDM phase (either in determining $P^{(1)}$
from the Dyson integral over $L$ in Eq. (25), or in the integral over
$\tau_{r}$ in Eq. (66)), which might significantly suppress the linear
contributions to the neutrino probabilities.

  * •

As stated, since the distribution of $z\left(f\right)$ under the null
hypothesis is frequency independent, the expected sensitivity of the
experiment for different ULDM masses is completely determined by the mass
dependence of $P^{(1)}$, assuming $m_{\phi}\lesssim 2/\tau_{r}$ and
$\tau_{c}\gtrsim\tau_{e}$ (or mass-independent coherence times). If the ULDM
coherence time is determined by $\tau_{c}=2\pi/\left(m_{\phi}\beta^{2}\right)$
with the virial velocity $\beta\approx 10^{-3}$ , then for
$m_{\phi}>2\pi/\left(\tau_{e}\beta^{2}\right)$ one can think of $P^{(1)}$ as
being effectively suppressed by an additional factor of
$\left(\frac{2\pi}{\tau_{e}m_{\phi}\beta^{2}}\right)^{\frac{1}{4}}$. For
$m_{\phi}>2/\tau_{r}$, the integration over time $\tau_{r}$ in Eq. (66) would
yield an effective additional $2/\left(m_{\phi}\tau_{r}\right)$ suppression to
$P^{(1)}$. Also note that our result coincides with the result obtained by
previous analyses assuming $\tau_{s}<\tau_{\phi}$ if one sets
$\tau_{r}=\tau_{s}$. Therefore, compared to previous analyses, the Rayleigh
method would yield bounds on the couplings of ULDM of masses
$2/\tau_{r}>m_{\phi}>2/\tau_{s}$ that are stronger by $m_{\phi}\tau_{s}$, and
bounds that are stronger by $\tau_{s}/\tau_{r}$ for $m_{\phi}>2/\tau_{r}$.

  * •

While the ULDM phase (denoted by $m_{\phi}t_{0}$) affects the sum of cosines
and the sum of sines in Eq (77) and Eq. (76), respectively, it does not affect
$z(f)$, as they are added in quadrature. Namely, the calculation of $z(f)$
does not require a knowledge of the ULDM phase, and the phase of the cosines
and sines used for it may be chosen arbitrarily. Consequently, in order to
calculate the expected bounds, one could treat the time-independent factors
associated with $\text{Re}\left(\kappa_{ij}\right)$ and
$\text{Im}\left(\kappa_{ij}\right)$ separately, as they will be added in
quadrature. Of course, scanning over the phase and examining the sum of
cosines and sum of sines independently would allow recovering the ULDM phase
information.

  * •

While in our calculation we used the simple modulation function $h(t_{i})$ in
Eq. (63), one may use other choices for $h$ in order to isolate or amplify
specific properties of the signal. Specifically, assigning
$h\left(t_{i}\right)=h\left(E\right)$ would effectively set a different,
energy-dependent, weight for the neutrino events, which could be helpful, for
example, in recovering energy-dependent highly oscillatory components of the
signal. Another useful choice could be made to isolate either the CPC or the
CPV probabilities in an experiment that allows for detecting both neutrinos
and anti-neutrinos, as

$h_{\pm}\left(t_{i}\right)=\left\\{\begin{matrix}1/n_{\nu}&\nu\text{ event
detected,}\\\ \pm 1/n_{\bar{\nu}}&\overline{\nu}\text{ event detected,}\\\
0&\text{no event detected\,,}\end{matrix}\right.$ (79)

where $n_{\nu}$ is the total number of neutrinos detected and $n_{\bar{\nu}}$
is the total number of anti-neutrinos detected. Generally speaking, for any
choice of a modulation function $h(\phi)$ which is time-independent (let it be
a function of the neutrino energy, flavor, lepton number etc.), we may define

$\displaystyle\tilde{r}\left(f\right)=\left(\left[\sum_{i=1}^{N_{\nu}}h\left({\phi}_{i}\right)\sin\left(2\pi
t_{i}f\right)\right]^{2}+\left[\sum_{i=1}^{N_{\nu}}h\left(\phi_{i}\right)\cos\left(2\pi
t_{i}f\right)\right]^{2}\right)\,.$ (80)

Then, the modified Rayleigh score

$\displaystyle\tilde{z}\left(f\right)=\frac{2}{\tilde{N}^{0}_{\nu}}\tilde{r}\left(f\right)\,,$
(81)

is $\chi_{2}^{2}$ distributed, and in the presence of a signal its expectation
value is

$\displaystyle\langle\tilde{z}\left(f\right)\rangle=2+\frac{\left(\tilde{N}_{s}^{\left(1\right)}\right)^{2}+\left(\tilde{N}_{c}^{\left(1\right)}\right)^{2}}{2\tilde{N}^{0}_{\nu}}\,,$
(82)

with

$\displaystyle\tilde{N}^{0}_{\nu}$ $\displaystyle=\left\langle\sum_{\rm
bkg}h^{2}\right\rangle=\tau_{e}\int\limits_{E_{\min}}^{E_{\max}}\int_{\phi}\frac{dN}{dEdtd\phi}P^{\left(0\right)}h^{2}d\phi
dE\,,$ (83) $\displaystyle\tilde{N}_{\nu,s}^{\left(1\right)}$
$\displaystyle=\left\langle\sum_{\rm
s}h\right\rangle=\tau_{e}\int\limits_{E_{\min}}^{E_{\max}}\int_{\phi}\frac{dN}{dEdtd\phi}P_{s}^{\left(1\right)}hd\phi
dE\,,$ (84) $\displaystyle\tilde{N}_{\nu,c}^{\left(1\right)}$
$\displaystyle=\left\langle\sum_{\rm
c}h\right\rangle=\tau_{e}\int\limits_{E_{\min}}^{E_{\max}}\int_{\phi}\frac{dN}{dEdtd\phi}P_{c}^{\left(1\right)}hd\phi
dE\,.$ (85)

Then, concretely, for $h_{\pm}$ one would obtain

$\displaystyle\langle\tilde{z}\left(f\right)\rangle_{\pm}\approx
2+\frac{\left(\epsilon_{s}^{\pm}\right)^{2}+\left(\epsilon_{c}^{\pm}\right)^{2}}{2\left(\frac{1}{n_{\nu}}+\frac{1}{n_{\bar{\nu}}}\right)}\,,$
(86)

with

$\displaystyle\epsilon_{s,c}^{+}$
$\displaystyle=\frac{\int\limits_{E_{\min}}^{E_{\max}}dEF\left(E\right)2\Sigma
P_{s,c}^{\left(1\right)}\left(E\right)}{\int\limits_{E_{\min}}^{E_{\max}}dEF\left(E\right)P^{\left(0\right)}\left(E\right)}\,,$
(87) $\displaystyle\epsilon_{s,c}^{-}$
$\displaystyle=\frac{\int\limits_{E_{\min}}^{E_{\max}}dEF\left(E\right)\Delta
P_{s,c}^{\left(1\right)}\left(E\right)}{\int\limits_{E_{\min}}^{E_{\max}}dEF\left(E\right)P^{\left(0\right)}\left(E\right)}\,,$
(88)

with $\Sigma P^{\left(1\right)}$ and $\Delta P^{\left(1\right)}$ defined in
Eq. (44) and Eq. (51), respectively.

##  5 Current and future-projected bounds

In this section we study the sensitivity of various neutrino oscillation
experiments to the ULDM modulation amplitude using the Rayleigh periodogram.
We derive bounds on the neutrino-ULDM couplings from existing experiments
which have already looked for time modulations in neutrino survival and
transition probabilities, and projected sensitivities for future experiments.
We do not use Monte-Carlo simulations, but use the scaling that we derived for
the expectation value of the magnitude of the Rayleigh spectrum peak at the
modulation frequency in Eq. (78), with the parameters of the experiment and of
the model. When calculating the height of the peak, we assumed that the ULDM
modulation amplitude is predominantly affected by one entry of the $\hat{y}$
matrix at a time. We also provide a separate analysis for the real and
imaginary components of the off-diagonal entries of $\hat{y}$, assuming
$\delta_{CP}=0$. Note that in our calculations we assumed a normal hierarchy,
and that $m_{1}=0$. The latter assumption can be relaxed by re-identifying
$\text{Im}\left(y_{1i}\right)\rightarrow\left(m_{i}-m_{1}\right)\text{Im}\left(y_{1i}\right)/m_{i}$
,
$\text{Re}\left(y_{1i}\right)\rightarrow\left(m_{i}+m_{1}\right)\text{Re}\left(y_{1i}\right)/m_{i}$
for the off-diagonal couplings, and $y_{22}\rightarrow
y_{22}-m_{1}y_{11}/m_{2}\,,y_{33}\rightarrow y_{33}-m_{1}y_{11}/m_{3}$ for the
diagonal ones. For the projected sensitivities, we find the curve in parameter
space $\hat{y}\left(m_{\phi}\right)$ for which the peak is expected to match
the 95% percentile of the $\chi^{2}_{2}$ distribution, characterizing the
background (corresponding to a CL=0.95). The corresponding bounds in terms of
the ULDM parameters are shown in Fig. 1.

###  5.1 Current bounds

####  5.1.1 Daya-Bay

The Daya Bay experiment collects $\bar{\nu}_{e}$ events from nearby reactors
which are approximately $800$ meters away from the detector. The collaboration
has searched for time-dependent modulations of the survival probability
$P_{\bar{e}\bar{e}}$, for frequencies between $5.9\times 10^{-5}$/sidereal
hour to 0.5/sidereal hour, and found no significant evidence [15]. We plot in
Fig. 1 the ULDM couplings corresponding to an expected sensitivity at CL=0.95,
assuming 621 days of data taking, with a total of 800 unoscillated events per
day [16].

####  5.1.2 Solar experiments - SNO

The Super-K and SNO experiments collect solar $\nu_{e}$ with event rates of
$15$/day [17] and $10$/day [18], respectively. Due to the MSW effect,
neutrinos with energy larger than $\sim 1$ MeV leave the sun as the $\nu_{2}$
mass states. Therefore, their survival probability does not oscillate with
propagation, and is given by
$P^{0}_{2e}=|U_{e2}|^{2}\approx\sin^{2}{\theta_{12}}$. This means that up to
day/night matter effects, the solar neutrino flux should be constant along the
year.

To calculate the first-order correction coming from the ULDM, we assume that
$m_{\phi}\tau_{\odot}\ll 1$, where $\tau_{\odot}$ is the longest time-scale
for neutrino effects in the sun, such that the neutrino leaves the sun as an
eigenstate of the instantaneous Hamiltonian,
$\mathcal{H}(t_{i})\approx\mathcal{H}_{0}+\tilde{m}\sin\left(m_{\phi}\left(t_{i}-t_{0}\right)\right)$,
denoted by $\nu_{m}\left(t_{i}\right)$. We then turn to our previously
obtained expressions for $P_{\alpha\beta}$, noting that
$\nu_{\alpha}=\nu_{m}\left(t_{i}\right)$ with $m=2$. Namely, the PMNS-like
matrix $U_{\alpha i}$ should be replaced by the projection of the
instantaneous mass eigenstate at time $t_{i}$ onto the ULDM-independent mass
eigenstates - $U_{\alpha i}\rightarrow T_{2i}\left(t_{i}\right)$ where
$T(t_{i})$ is the diagonalizing matrix of $\mathcal{H}(t_{i})$ in the ULDM-
independent mass basis. The deviation of $T(t_{i})$ from the identity yields a
first order correction in the ULDM coupling to $\mathcal{U}^{(0)}$. When
calculating $\mathcal{U}^{(1)}$, one may set $U_{\alpha i}=\delta_{2i}$, since
it is already linear in the ULDM potential. For off-diagonal ULDM couplings
$\hat{y}_{i2}$ with $i\neq 2$, the leading ULDM correction in our formalism is
given by

$\displaystyle P^{(1)}_{2e}$
$\displaystyle=2\frac{\phi_{0}}{E}\text{Re}\left(U^{*}_{e2}U_{ei}\tilde{m}_{i2}\right)\left(-\frac{E\sin\left(m_{\phi}\left(t-L-
t_{0}\right)\right)}{\Delta m^{2}_{i2}}\cos{\frac{\Delta
m^{2}_{i2}L}{2E}}-\cos{\frac{\Delta
m^{2}_{i2}L}{4E}}\text{Im}\left(\kappa_{i2}\right)+\sin{\frac{\Delta
m^{2}_{i2}L}{4E}}\text{Re}\left(\kappa_{i2}\right)\right)$
$\displaystyle-2\frac{\phi_{0}}{E}\text{Im}\left(U^{*}_{e2}U_{ei}\tilde{m}_{i2}\right)\left(\frac{E\sin\left(m_{\phi}\left(t-L-
t_{0}\right)\right)}{\Delta m^{2}_{i2}}\sin{\frac{\Delta
m^{2}_{i2}L}{2E}}+\cos{\frac{\Delta
m^{2}_{i2}L}{4E}}\text{Re}\left(\kappa_{i2}\right)+\sin{\frac{\Delta
m^{2}_{i2}L}{4E}}\text{Im}\left(\kappa_{i2}\right)\right)\,,$ (89)

whereas for the diagonal coupling $\hat{y}_{22}$ and for all $\hat{y}_{ij}$
with $j\neq 2$ , there are no first-order contributions to the solar
probability.

The SNO collaboration looked for a time variation in the solar neutrino flux
with periods between $10$ minutes to $1$ day, and found a null result with
confidence level of $90\%$ [13]. The search was expected to be sensitive (have
a 90% probability of making a CL=0.99 discovery) to sinusoidal modulations at
these frequencies with an amplitude of 0.12 or greater, compared to the time-
averaged flux. We may then interpret this result as a bound on the ULDM
parameters, assuming $m_{\phi}$ is equal to the modulation frequency. In the
$m_{\phi}\ll\frac{\Delta m^{2}_{i2}}{2E}$ limit, which applies throughout the
frequency range of the analysis of [13], the correction to the solar
probability in Eq. (89) becomes

$\displaystyle P^{(1)}_{2e}$ $\displaystyle\approx\frac{2\phi_{0}}{\Delta
m^{2}_{i2}}\text{Re}\left(U_{e2}U^{*}_{ei}\tilde{m}_{i2}^{*}\right)\Big{(}-\sin{\left(m_{\phi}\left(t-t_{0}\right)\right)}+\frac{2Em_{\phi}}{\Delta
m^{2}_{i2}}\sin{\frac{\Delta m^{2}_{i2}L}{2E}}\cos\left(m_{\phi}\left(t-L-
t_{0}\right)\right)\Big{)}$ $\displaystyle+4\frac{Em_{\phi}\phi_{0}}{{\Delta
m^{2}_{i2}}^{2}}\text{Im}\left(U_{e2}U^{*}_{ei}\tilde{m}_{i2}^{*}\right)\Big{(}-\cos\left(m_{\phi}\left(t-t_{0}\right)\right)+\cos{\frac{\Delta
m^{2}_{i2}L}{2E}}\cos\left(m_{\phi}\left(t-L-t_{0}\right)\right)\Big{)}\,.$
(90)

###  5.2 Projected sensitivities

####  5.2.1 Reactor experiments

In addition to current bounds from the Daya-Bay experiment, we study the
projected sensitivities of ULDM parameters from other reactor neutrino
experiments. We consider JUNO that measures the survival probability of
$\bar{\nu}_{e}$ that travel $\sim 53$ km. The energy resolution and baseline
length of JUNO enable it to be sensitive to both neutrino mass squared
differences $\Delta m^{2}_{21}$, $\Delta m^{2}_{31}$ and thus to determine the
neutrino mass hierarchy. It has an expected flux of $83$ unoscillated events
per day [19], and its projected sensitivity to ULDM parameters is presented in
Fig. 1.

We did not consider the KamLAND experiment in our analysis, as due to its long
baseline of $\sim 180$ km has a relatively low flux of $1-2$ events per day.
As opposed to accelerator experiments which fire a narrow beam of neutrinos,
reactor experiments fire neutrinos uniformly in every direction and therefore
their flux decays with $L^{-2}$. It also collects events from different
reactors with different distances, which considerably complicates the
analysis.

####  5.2.2 Accelerator experiments

Accelerator neutrino experiments fire a beam of protons towards a target that
produces charged pions. These can be manipulated into a narrow beam and decay
mostly into $\nu_{\mu}$ or $\bar{\nu}_{\mu}$, depending if the positive or
negative pions are focused. These experiments therefore measure $P_{\mu\mu}$
and $P_{\mu e}$ (or $P_{\bar{\mu}\bar{\mu}}$ and $P_{\bar{\mu}\bar{e}}$)
simultaneously. We studied the DUNE, ESS (far and near detectors), and Hyper-K
experiments, out of which the latter has produced overall slightly better
sensitivities to the ULDM parameters, and only these results are presented in
Fig. 1. Hyper-K is expected to collect approximately $2000$ to $4000$ signal
events in each mode of neutrinos that travel $295$ km, over a run of 10 years
[20].

Figure 1:  Current bounds on the ULDM parameter space from Daya-Bay and SNO
experiments, alongside the projected sensitivities of JUNO and Hyper-K. To
generate these plots we assumed normal hierarchy, $m_{1}=0$, and
$\delta_{\text{CP}}=0$.

###  5.3 Asymptotic behavior

To qualitatively understand the various plots in Fig. 1, let us consider the
asymptotic behavior following from the analytical expressions we found in Sec.
3. Note that we only account for the first order corrections to the neutrino
oscillation probability. Therefore, near the resonance $m_{\phi}=\Delta
m^{2}/2E$ [9] the result may be inaccurate, as higher order terms would yield
considerable contributions.

Let us first discuss the behavior of the CPC ULDM effects (corresponding to
the plots presenting the real parts of the couplings, assuming
$\delta_{\rm{CP}}=0$). For couplings that are diagonal in the unperturbed mass
basis, the correction to the transition or survival probability in Eq. (3.2)
is proportional to

$\displaystyle\Sigma P^{\left(1\right)}$
$\displaystyle\propto\sin\left(\frac{m_{\phi}L}{2}\right)\frac{2y_{ii}\phi_{0}}{Em_{\phi}}\propto\begin{cases}\frac{\text{Re}\left(y_{ii}\right)}{m_{\phi}}&m_{\phi}L\ll
1\,,\\\
\sin\left(\frac{m_{\phi}L}{2}\right)\frac{\text{Re}\left(y_{ii}\right)}{m_{\phi}^{2}}&m_{\phi}L\gg
1\,.\end{cases}$ (91)

Therefore, in an experiment sensitive to a certain value of $P^{(1)}$, the
corresponding diagonal coupling $y_{ii}$ that the experiment will be sensitive
to is proportional to

$\displaystyle\text{Re}\left(y_{ii}\right)$
$\displaystyle\propto\begin{cases}m_{\phi}&m_{\phi}L\ll 1\,,\\\
m_{\phi}^{2}/\sin\left(m_{\phi}L/2\right)&m_{\phi}L\gg 1\,.\end{cases}$ (92)

We can see the transition between the two limits in the two bottom plots of
Fig. 1.

As expected, when the couplings are off-diagonal, the probability also depends
on the neutrino oscillation frequency, and we obtain the following relation

$\displaystyle\text{Re}\left(y_{ij}\right)$
$\displaystyle\propto\begin{cases}m_{\phi}&L^{-1}\gg m_{\phi}\,,\\\
m_{\phi}/\sqrt{A^{2}\cos^{2}\left(m_{\phi}L/2\right)+B^{2}\sin^{2}\left(m_{\phi}L/2\right)}&L^{-1}\ll
m_{\phi}\ll\Delta m_{ij}^{2}/2E\,,\\\
m_{\phi}^{2}/\sqrt{A^{2}\cos^{2}\left(m_{\phi}L/2\right)+B^{2}\sin^{2}\left(m_{\phi}L/2\right)}&L^{-1},\Delta
m_{ij}^{2}/2E\ll m_{\phi}\,,\end{cases}$ (93)

where $A$ and $B$ are determined by the coefficients of $\text{Re}(\kappa)$
and $\text{Im}(\kappa)$ in $P^{(1)}$, which depend on PMNS elements. If $A$
and $B$ are of similar magnitude, as is the case when $\text{Re}(\kappa)$ and
$\text{Im}(\kappa)$ yield similar contributions to the probability, we expect
the oscillations in $m_{\phi}$ to be suppressed. However, if one is much
larger than the other, there are approximate poles of “bad sensitivity” around
the nodes of the bigger term, at which the bound follows the inverse of the
smaller term (for example, if $A^{2}/B^{2}\gg 1$, then around
$m_{\phi}L/2=n\pi$ the bound is proportional to $1/A$, while around
$m_{\phi}L/2=\pi/2+n\pi$ the bound is proportional to $1/B$). This is the case
when either $\text{Re}(\kappa)$ or $\text{Im}(\kappa)$ dominate the
probability. In the case of diagonal entries of $\hat{y}$, only
$\text{Re}(\kappa)$ affects the probability, and therefore the bound
oscillates rapidly for larger $m_{\phi}$. For off-diagonal $y$ entries, the
modification to the CPC survival probability is affected only by
$\text{Re}(\kappa)$, while the transition probability is affected also by
$\text{Im}(\kappa)$. Therefore, the sensitivity for off-diagonal CPC couplings
from measurements of survival probabilities will oscillate rapidly for larger
$m_{\phi}$, while these oscillations are suppressed for transition
probabilities.

The two extreme limits of Eq. (93) are similar to the diagonal couplings case,
while we may also obtain an in-between behavior if $\Delta
m_{ij}^{2}L/\left(2E\right)>1$. This region may be probed due to the neutrino
mass hierarchy $\Delta m^{2}_{21}\ll\Delta m^{2}_{31},\Delta m^{2}_{32}$. The
JUNO experiment is designed in a way such that $\Delta
m^{2}_{21}L/\left(2E\right)={\cal O}(1)$, and $\Delta
m^{2}_{31}L/\left(2E\right)\gg 1$. This allows for the in-between region where
$L^{-1}\ll m_{\phi}\ll\Delta m^{2}_{31}/\left(4E\right)$. We can indeed see
this region in the plots of $y_{13},y_{23}$, but not for $y_{12}$.

The behavior of the CPV couplings is similar to that of the CPC couplings,
with one exception: the 2-$\nu$ CPV effect. The 2-$\nu$ CPV is the only CPV
effect in survival probabilities, and its mass dependence follows that of
$\text{Im}\left(\kappa_{ij}\right)$ in Eq. (29), and thus we expect that the
bounds from survival probabilities follow

$\displaystyle\text{Im}\left(y_{ij}\right)$
$\displaystyle\propto\begin{cases}\text{const}\left(m_{\phi}\right)&L^{-1}\gg
m_{\phi}\,,\\\
m_{\phi}/\sqrt{A^{2}\cos^{2}\left(m_{\phi}L/2\right)+B^{2}\sin^{2}\left(m_{\phi}L/2\right)}&L^{-1}\ll
m_{\phi}\ll\Delta m_{ij}^{2}/2E\,,\\\
m_{\phi}^{2}/\sqrt{A^{2}\cos^{2}\left(m_{\phi}L/2\right)+B^{2}\sin^{2}\left(m_{\phi}L/2\right)}&L^{-1},\Delta
m_{ij}^{2}/2E\ll m_{\phi}\,.\end{cases}$ (94)

The main difference from the CPC bounds is in the low mass region, where the
bound flattens. This is a result of the fact that in the two generations case,
as we discussed in Sec 2, the only observable CPV effect comes from the time-
derivative of the ULDM-induced phase. The 2-$\nu$ CPV effect occurs also in
transition probabilities, but overcomes the 2-$\nu$ CPV effect in the low mass
region if $\Delta m^{2}_{21}/E\ll m_{\phi}\ll\Delta m^{2}_{31}/E,\Delta
m^{2}_{32}/E$ (this effectively sets $\Delta m^{2}_{21}=0$ which gives the
two-neutrino picture). We can see this in-between behavior of the transition
probability in Hyper-K for $y_{13}$ and $y_{23}$, but not for $y_{12}$. Since
the 2-$\nu$ CPV effect takes over at larger masses (for all off-diagonal
couplings), the CPV transition probabilities would inherit the behavior of the
CPV survival probabilities, which are solely affected by
$\imaginary\left(\kappa\right)$, and thus they would significantly oscillate
as a function of $m_{\phi}$.

##  6 Conclusions

The inclusion of a scalar ULDM field that couples to neutrinos provides new
rich phenomenology. In this work, we have presented a generic analysis of the
impact of time varying ULDM interactions with neutrinos. For the case when the
ULDM rapidly oscillates during neutrino propagation, we have determined the
corrections to the survival and transition probabilities for both CPC and CPV
quantities. Some of the novel effects that appear are the following:

  * •

A new CPV effect is present even in the two neutrino-like case. We have also
explicitly checked that in the three generation case this new effect is also
present. Moreover, for three generations of neutrinos, the CPV effect does not
disappear in the limit of one of the $\Delta m_{ij}^{2}\rightarrow 0$, and
both survival and transition probabilities are modified by the new CPV effect.

  * •

The phenomena we are interested in presents oscillatory behavior and thus we
can look at the spectral power of the neutrino probabilities for angular
frequencies defined by $m_{\phi}$ to determine the parameters in the model.

Interestingly, we have identified that several various neutrino oscillation
experiments can be sensitive to the fast variations of the ULDM field. Given
this we have determined the bounds on CPC and CPV couplings from the neutrino
survival probability measurements at Daya Bay and SNO. Furthermore, we have
determined the sensitivity to measure the couplings at future reactor (JUNO)
and long baseline neutrino oscillation experiments (DUNE, ESS and Hyper-K).

## Acknowledgements

YN is the Amos de-Shalit chair of theoretical physics, and is supported by
grants from the Israel Science Foundation (grant number 1124/20), the United
States-Israel Binational Science Foundation (BSF), Jerusalem, Israel (grant
number 2018257), by the Minerva Foundation (with funding from the Federal
Ministry for Education and Research), and by the Yeda-Sela (YeS) Center for
Basic Research. The work of GP is supported by grants from BSF-NSF, Friedrich
Wilhelm Bessel research award, GIF, ISF, Minerva, SABRA - Yeda-Sela - WRC
Program, the Estate of Emile Mimran, and The Maurice and Vivienne Wohl
Endowment. IS is supported by a fellowship from the Ariane de Rothschild Women
Doctoral Program.

##  Appendix A Two Generations

For diagonal ULDM couplings, $\hat{y}_{ii}$, and $k\neq i$, the survival
probability is modified by

$\displaystyle P^{(1)}_{\alpha\alpha}\left(t\right)$
$\displaystyle=-\sin\left(m_{\phi}\left(t-\frac{L}{2}-t_{0}\right)\right)\sin\left(\frac{m_{\phi}L}{2}\right)\frac{\sin^{2}{2\theta}\hat{y}_{ii}m_{i}\phi_{0}}{Em_{\phi}}\sin{\left(\frac{\Delta
m^{2}_{ik}L}{2E}\right)}\,,$ (95)

where $\theta$ is the mixing angle between the unperturbed mass eigenstates
and interaction eigenstates. For non-diagonal ULDM couplings $\hat{y}_{ij}$,
with $i\neq j$

$\displaystyle P^{(1)}_{\alpha\alpha}$
$\displaystyle=\frac{\phi_{0}\sin{4\theta}}{E\left(\left(\frac{\Delta
m^{2}_{ij}}{2E}\right)^{2}-m_{\phi}^{2}\right)}\times$
$\displaystyle\times\Big{[}$
$\displaystyle\text{Im}\left(\hat{y}_{12}\right)\left(m_{1}-m_{2}\right)\cos\left(m_{\phi}\left(\left(t-L-
t_{0}\right)+\frac{L}{2}\right)\right)\left(\frac{\Delta
m^{2}_{21}}{4E}\sin{\frac{\Delta
m^{2}_{21}L}{2E}}\sin{\frac{m_{\phi}L}{2}}-m_{\phi}\sin^{2}{\frac{\Delta
m^{2}_{21}L}{4E}}\cos{\frac{m_{\phi}L}{2}}\right)$ $\displaystyle+$
$\displaystyle\text{Re}\left(\hat{y}_{12}\right)\left(m_{1}+m_{2}\right)\sin\left(m_{\phi}\left(\left(t-L-
t_{0}\right)+\frac{L}{2}\right)\right)\left(\frac{m_{\phi}}{2}\sin{\frac{\Delta
m^{2}_{21}L}{2E}}\sin{\frac{m_{\phi}L}{2}}-\frac{\Delta
m^{2}_{21}}{2E}\sin^{2}{\frac{\Delta
m^{2}_{21}L}{4E}}\cos{\frac{m_{\phi}L}{2}}\right)\Big{]}\,.$ (96)

##  Appendix B Mean value of Rayleigh power at the modulation frequency

Recall the time scales of the problem:

  1. 1.

$\tau_{e}$: The running time of the experiment.

  2. 2.

$\tau_{s}=\frac{N_{\nu}}{\tau_{e}}$: The average spacing between events, where
$N_{\nu}$ is the total number of measured events.

  3. 3.

$\tau_{r}$: The resolution of the clock that times the events.

Therefore, during the experiment running time, there are
$N_{t}=\frac{\tau_{e}}{\tau_{r}}$ clock ticks, each with small duration
$\tau_{r}$. At each clock tick we define the function:

$h\left(t_{i}\right)=\left\\{\begin{matrix}1&\text{event detected,}\\\
0&\text{no event detected.}\end{matrix}\right.$ (97)

We calculate the expectation value of the Rayleigh power spectrum, namely

$\displaystyle\left\langle z\left(f\right)\right\rangle$
$\displaystyle=\frac{2}{N_{\nu}}\left(\left\langle\left|\sum_{i=1}^{N_{t}}h\left(t_{i}\right)\cos\left(2\pi
t_{i}f\right)\right|^{2}\right\rangle+\left\langle\left|\sum_{i=1}^{N_{t}}h\left(t_{i}\right)\sin\left(2\pi
t_{i}f\right)\right|^{2}\right\rangle\right),$ (98)

were $f$ is the modulation frequency of the probability modulation,
corresponding to $m_{\phi}/2\pi$. We calculate the expectation value of a
generic function $F\left(h\left(t_{i}\right),t_{i}\right)$ in the sense that

$\left\langle
F\left(h\left(t_{i}\right),t_{i}\right)\right\rangle=F\left(1,t_{i}\right)\cdot
P\left(h\left(t_{i}\right)=1\right)+F\left(0,t_{i}\right)\cdot
P\left(h\left(t_{i}\right)=0\right),$ (99)

where the probability that an event is detected is given by

$\displaystyle P\left(h\left(t_{i}\right)=1\right)$
$\displaystyle=\frac{N_{\nu}}{N_{t}}\left[1+\epsilon\sin\left(2\pi
t_{i}f\right)\right]=N_{\nu}\frac{\tau_{\text{r}}}{\tau_{\text{e}}}\left[1+\epsilon\sin\left(2\pi
t_{i}f\right)\right],$ (100)

where $\epsilon$ is the amplitude of the modulation. We assumed for simplicity
that the phase of the modulation is zero, such that it appears as a sine, but
the final result holds for a generic phase. Let us start with calculating the
expectation value of the sine term. Notice that

$\left\langle\left|\sum_{i=1}^{N_{t}}h\left(t_{i}\right)\sin\left(2\pi
t_{i}f\right)\right|^{2}\right\rangle=\text{Var}\left(\sum_{i=1}^{N_{t}}h\left(t_{i}\right)\sin\left(2\pi
t_{i}f\right)\right)+\left\langle\left|\sum_{i=1}^{N_{t}}h\left(t_{i}\right)\sin\left(2\pi
t_{i}f\right)\right|\right\rangle^{2}\,.$ (101)

Calculating the variance term first,

$\displaystyle\text{Var}\left[h\left(t_{i}\right)\sin\left(2\pi
t_{i}f\right)\right]$ $\displaystyle=\left\langle
h^{2}\left(t_{i}\right)\sin^{2}\left(2\pi
t_{i}f\right)\right\rangle-\left\langle h\left(t_{i}\right)\sin\left(2\pi
t_{i}f\right)\right\rangle^{2}=\left(P\left(t_{i}\right)-P^{2}\left(t_{i}\right)\right)\sin^{2}\left(2\pi
t_{i}f\right)$
$\displaystyle\approx\frac{N_{\nu}}{N_{t}}\left[1+\epsilon\sin\left(2\pi
t_{i}f\right)\right]\sin^{2}\left(2\pi t_{i}f\right))\,.$ (102)

In the last approximation we assume $\frac{N_{\nu}}{N_{t}}\ll 1$, which is
equivalent to the statement that the resolution of the clock is much better
than the typical time between events. The variance of sum, is simply the sum
of variances, since in our case the measurements at different times are
uncorrelated, thus

$\displaystyle\text{Var}\left(\sum_{i=1}^{N_{t}}h\left(t_{i}\right)\sin\left(2\pi
t_{i}f\right)\right)$
$\displaystyle=\sum_{i=1}^{N_{t}}\frac{N_{\nu}}{N}\left[1+\epsilon\sin\left(2\pi
t_{i}f\right)\right]\sin^{2}\left(2\pi t_{i}f\right)\,.$ (103)

Changing variables:

$\displaystyle t_{n}$ $\displaystyle\rightarrow n\tau_{\text{res}}\,,$ (104)
$\displaystyle 0$ $\displaystyle<n\leq N_{t}\,,$ (105)

we obtain

$\displaystyle\text{Var}\left(\sum_{i=1}^{N_{t}}h\left(t_{i}\right)\sin\left(2\pi
t_{i}f\right)\right)$
$\displaystyle=\frac{N_{\nu}}{N_{t}}\left[\sum_{n=1}^{N_{t}}\sin^{2}\left(2\pi
n\tau_{\text{res}}f\right)+\epsilon\sum_{n=1}^{N_{t}}\sin^{3}\left(2\pi
n\tau_{\text{res}}f\right)\right]\,.$ (106)

The first term yields

$\displaystyle\sum_{n=1}^{N_{t}}\sin^{2}\left(2\pi n\tau_{\text{res}}f\right)$
$\displaystyle\approx\int\limits_{0}^{N_{t}}dz\sin^{2}\left(2\pi
z\tau_{\text{res}}f\right)=\frac{N_{t}}{2}\left[1-\text{sinc}\left(4\pi
N_{t}\tau_{\text{res}}f\right)\right]\approx\frac{N_{t}}{2}\,,$ (107)

where we use the approximation

$N_{t}\tau_{\text{res}}f=\tau_{\exp}f\gg 1.~{}$ (108)

The second term gives

$\displaystyle\sum_{n=1}^{N_{t}}\sin^{3}\left(2\pi n\tau_{\text{res}}f\right)$
$\displaystyle\approx\int\limits_{0}^{N_{t}}dz\sin^{3}\left(2\pi
z\tau_{\text{res}}f\right)=\frac{2N_{t}\left(2+\cos\left(2fN_{t}\pi\tau_{\text{res}}\right)\right)\sin^{3}\left(fN_{t}\pi\tau_{\text{res}}\right)}{3}\text{sinc$\left(fN_{t}\pi\tau_{\text{res}}\right)$}\,,$
(109)

and may therefore be neglected under the approximation in Eq. (108).
The second term in Eq. (101) is

$\displaystyle\left\langle\left|\sum_{i=1}^{N_{t}}h\left(t_{i}\right)\sin\left(2\pi
t_{i}f\right)\right|\right\rangle^{2}$
$\displaystyle=\left|\sum_{i=1}^{N_{t}}P(t_{i})\sin\left(2\pi
t_{i}f\right)\right|^{2}=\left[\sum_{i=1}^{N_{t}}\frac{N_{\nu}}{N_{t}}\left[1+\epsilon\sin\left(2\pi
t_{i}f\right)\right]\sin\left(2\pi
t_{i}f\right)\right]^{2}=\frac{\epsilon^{2}}{4}N_{\nu}^{2}.$ (110)

Finally we get

$\displaystyle\left\langle\left|\sum_{i=1}^{N_{t}}h\left(t_{i}\right)\sin\left(2\pi
t_{i}f\right)\right|^{2}\right\rangle$
$\displaystyle=\frac{N_{\nu}}{2}+\frac{\epsilon^{2}}{4}N_{\nu}^{2}\,,$ (111)
$\displaystyle\left\langle\left|\sum_{i=1}^{N_{t}}h\left(t_{i}\right)\cos\left(2\pi
t_{i}f\right)\right|^{2}\right\rangle$ $\displaystyle=\frac{N_{\nu}}{2}\,,$
(112)

and thus

$\displaystyle\left\langle z\left(f\right)\right\rangle$
$\displaystyle=2+\frac{\epsilon^{2}}{2}N_{\nu}\,.$ (113)

## References

  * [1] G. Krnjaic, P. A. N. Machado, and L. Necib, Distorted neutrino oscillations from time varying cosmic fields, Phys. Rev. D 97 (2018), no. 7 075017, [1705.06740].
  * [2] V. Brdar, J. Kopp, J. Liu, P. Prass, and X.-P. Wang, Fuzzy dark matter and nonstandard neutrino interactions, Phys. Rev. D 97 (2018), no. 4 043001, [1705.09455].
  * [3] F. Capozzi, I. M. Shoemaker, and L. Vecchi, Neutrino Oscillations in Dark Backgrounds, JCAP 07 (2018) 004, [1804.05117].
  * [4] A. Berlin, Neutrino Oscillations as a Probe of Light Scalar Dark Matter, Phys. Rev. Lett. 117 (2016), no. 23 231801, [1608.01307].
  * [5] A. Dev, P. A. N. Machado, and P. Martínez-Miravé, Signatures of ultralight dark matter in neutrino oscillation experiments, JHEP 01 (2021) 094, [2007.03590].
  * [6] M. Losada, Y. Nir, G. Perez, and Y. Shpilman, Probing scalar dark matter oscillations with neutrino oscillations, JHEP 04 (2022) 030, [2107.10865].
  * [7] E. J. Chun, Neutrino Transition in Dark Matter, 2112.05057.
  * [8] G.-y. Huang and N. Nath, Neutrino meets ultralight dark matter: 0$\nu$$\beta$$\beta$ decay and cosmology, JCAP 05 (2022), no. 05 034, [2111.08732].
  * [9] M. Losada, Y. Nir, G. Perez, I. Savoray, and Y. Shpilman, Parametric resonance in neutrino oscillations induced by ultra-light dark matter and implications for KamLAND and JUNO, 2205.09769.
  * [10] S. Bilenky, J. Hošek, and S. Petcov, On the oscillations of neutrinos with dirac and majorana masses, Phys. Lett. B 94 (1980), no. 4 495–498.
  * [11] C. Giunti, No Effect of Majorana Phases in Neutrino Oscillations, Phys. Lett. B 686 (2010) 41–43, [1001.0760].
  * [12] L. Wolfenstein, Neutrino oscillations in matter, Phys. Rev. D 17 (1978) 2369–2374.
  * [13] SNO Collaboration, B. Aharmim et al., Searches for High Frequency Variations in the 8B Solar Neutrino Flux at the Sudbury Neutrino Observatory, Astrophys. J. 710 (2010) 540–548, [0910.2433].
  * [14] E. Gross and O. Vitells, Trial factors for the look elsewhere effect in high energy physics, Eur. Phys. J. C 70 (2010) 525–530, [1005.1891].
  * [15] Daya Bay Collaboration, D. Adey et al., Search for a time-varying electron antineutrino signal at Daya Bay, Phys. Rev. D 98 (2018), no. 9 092013, [1809.04660].
  * [16] Daya Bay Collaboration, F. P. An et al., Measurement of electron antineutrino oscillation based on 1230 days of operation of the Daya Bay experiment, Phys. Rev. D 95 (2017), no. 7 072006, [1610.04802].
  * [17] Super-Kamiokande Collaboration, J. Yoo et al., A Search for periodic modulations of the solar neutrino flux in Super-Kamiokande I, Phys. Rev. D 68 (2003) 092002, [hep-ex/0307070].
  * [18] SNO Collaboration, N. Tolich, Sudbury Neutrino Observatory: Latest results and future prospects, Nucl. Phys. B Proc. Suppl. 217 (2011) 107–112.
  * [19] A. Giaz, Status and perspectives of the JUNO experiment, in Prospects in Neutrino Physics, pp. 53–60, 4, 2018.  1804.03575.
  * [20] Hyper-Kamiokande Working Group Collaboration, E. Kearns et al., Hyper-Kamiokande Physics Opportunities, in Community Summer Study 2013: Snowmass on the Mississippi, 9, 2013.  1309.0184.

  *[ULDM]: ultra-light dark matter
  *[CPV]: CP violating
  *[CPC]: CP conserving
  *[SM]: Standard Model
  *[VEV]: vacuum expectation value
  *[EOM]: Equation of Motion
  *[EOMs]: Equation of Motion