# Neutrino Oscillations as a Gravitational Wave Detector?

Dominik<EMAIL_ADDRESS>, Sara
<EMAIL_ADDRESS>, Heinrich Päs1333heinrich.paes@tu-
dortmund.de , Mustafa<EMAIL_ADDRESS>
1 Fakultät für Physik, Technische Universität Dortmund, Germany

###### Abstract

Gravitational Waves (GWs) can alter the neutrino propagation distance and thus
affect neutrino oscillations. This can result in a complete disappearance of
the oscillatory behavior that competes with other sources of neutrino
decoherence. We develop a set of criteria that determines under which
conditions neutrino oscillations are sensitive to this effect, and discuss
three concrete scenarios for neutrinos from astrophysical sources. We find
that neutrino oscillations may probe so far unexplored regions of the GW
parameter space.

##  1 Introduction

Gravitational wave (GW) astronomy provides a new window for our exploration of
the universe. This includes the observation of previously inaccessible
phenomena, such as the merger of black holes [1] and allows for the
investigation and detection of neutron star mergers [2], collisions between
neutron stars and black holes [3], and potentially the explosion of supernovae
in the future [4]. Moreover, GWs can be employed in the search for exotic
objects such as cosmic strings [5] or primordial black holes [6], and thus
provide clues to physics beyond the Standard Model (SM). In principle GWs in
the sub $\mathrm{nHz}$ range may even provide novel possibilities to constrain
inflationary models [7]. Recent achievements like NANOGrav’s detection of
evidence for a GW background using pulsar timing arrays (PTAs) [8] prove that
we are at the beginning of an exciting new era of GW and multimessenger
astronomy. Further efforts are underway to expand the capabilities of existing
gravitational wave detectors such as LIGO and Virgo, while plans for next
generation detectors like the Einstein Telescope [9] and LISA (Laser
Interferometer Space Antenna) [10] promise an even greater sensitivity and a
broader range of observations in the near future.

While much attention has been focused on the detection and characterization of
GWs using conventional methods, such as laser interferometry [11], there is a
growing recognition of the potential role that neutrinos could play in this
area of research [12, 13, 14, 15, 16, 17, 18, 19, 20, 21]. As has been already
pointed out in references [22, 23, 24, 25, 26, 27, 28, 29], gravitational
fields modify the time evolution phase of quantum particles. Since neutrinos
are able to traverse vast distances through astrophysical environments without
being deflected from their trajectory, they could offer unique insights into
the GW landscape. One promising research direction is to study how GWs affect
quantum interference phenomena like neutrino flavor oscillations.

So far, most works on this subject consider the influence of a stochastic GW
background (SGWB) as a source of quantum decoherence in neutrino oscillations,
see for example references [12, 13, 14, 15, 16, 17, 18, 19]. In this work we
argue that not only stochastic perturbations of the neutrino evolution phase
can cause decoherence in the neutrino system, but that also GW signals from
individual sources will at least in principle have a similar effect that may
be searched for. We discuss in detail the challenges and prospects of such an
approach and investigate the importance of the interplay of GW induced
decoherence with other sources of decoherence, as described by references [30,
31, 32, 33, 34, 35, 36, 37, 38, 39, 40]. These decoherence effects arise, for
example, from the separation of neutrino wave packets, potentially macroscopic
neutrino production regions and finite energy binning and result in a similar
damping of the oscillation behavior of neutrino flavor transitions. Therefore
it is important to account for these effects in a coherent statistical
treatment in order to determine if an experiment would, in principle, be
sensitive to the GW induced effect, a step that has not yet been undertaken in
previous studies.

This paper is organized as follows. In section 2 we develop a simple yet
robust heuristic model accounting for the effects of the curved spacetime
induced by GWs on neutrino propagation by considering deviations from the
static trajectory approximation. A time average of the oscillation probability
over the running time of a neutrino experiment is employed in order to account
for a varying travel distance and to predict the expected neutrino oscillation
pattern. In section 3 we develop several criteria that determine whether
neutrino oscillations can be sensitive to GWs, and constrain the amplitudes
and frequencies of GWs that can be probed within this framework. In section 4
we study concrete scenarios for the impact of GWs on neutrino oscillations,
starting with baselines smaller than the diameter of the solar system in
section 4.1. Baselines exceeding the diameter of the solar system using the
example of the Vela pulsar are considered in section 4.2, and in section 4.3
we investigate neutrino oscillations in presence of a SGWB. A summary and
conclusions are provided in section 5. Technical details regarding the
calculation and the statistical method employed for the analyses can be found
in the appendices, where we also comment on further potential scenarios.

##  2 Modified Oscillation Probability in the Presence of Gravitational Waves

To investigate how the propagation of neutrinos through a curved spacetime
induced by GWs affects their flavor oscillations, we start with the standard
neutrino oscillation probability

$\displaystyle\hat{P}_{ab}(E,L)=\sum_{j}|U_{aj}|^{2}|U_{bj}|^{2}+2\sum_{j<k}\mathrm{Re}\left(U_{aj}^{\ast}U_{bj}U_{ak}U_{bk}^{\ast}\exp\left[-2\pi
i\frac{L}{L_{jk}^{\mathrm{osc}}}-\mathcal{D}_{jk}(E,L)\right]\right)\,,$ (1)

where $\mathcal{D}_{jk}$ is a positive damping term describing the decoherence
due to wave packet separation and $L_{jk}^{\mathrm{osc}}=4\pi E/\Delta
m_{jk}^{2}$ is the oscillation length. The matrix $U$ describes the neutrino
mixing, e.g. for three neutrino flavors $U$ corresponds to the
Pontecorvo–Maki–Nakagawa–Sakata (PMNS) matrix, $\Delta
m_{jk}^{2}=m_{j}^{2}-m_{k}^{2}$ are the squared mass splittings of the
neutrino mass eigenstates and $E$ is the average energy of the neutrino
system. Depending on the shape of the neutrino wave packets the decoherence
function $\mathcal{D}_{jk}$ can differ. A common assumption for the neutrino
wave packets is that they are of Gaussian shape implying

$\displaystyle\mathcal{D}_{jk}(E,L)=\left(\frac{L}{L_{jk}^{\mathrm{coh}}}\right)^{2}\quad\text{with}\quad
L_{jk}^{\mathrm{coh}}=2\sqrt{2}\frac{\sigma_{x}}{|\Delta
v_{jk}|}\quad\text{and}\quad\Delta v_{jk}=\frac{\Delta m_{jk}^{2}}{2E^{2}}\,,$
(2)

with the position space wave packet width $\sigma_{x}$ that we take to be
equal for all mass eigenstates for simplicity.

We now describe how GWs influence neutrino oscillations by modifying the
length $L$ in the standard probability, cf. equation (1), and replacing it
with the propagation distance of neutrinos in curved spacetime. For this, we
consider the average geodesic $\gamma$ that the neutrino wave packets follow
which now are slightly perturbed due to the presence of the GW,

$\displaystyle\gamma_{0}^{\mu}(\tau)$
$\displaystyle=\underbrace{\frac{X_{D}^{\mu}-X_{P}^{\mu}}{\tau_{D}-\tau_{P}}}_{u^{\mu}}(\tau-\tau_{P})+X_{P}^{\mu}\quad\longrightarrow\quad\gamma^{\mu}(\tau)=\gamma_{0}^{\mu}(\tau)+\Delta\gamma^{\mu}(\tau)\,,$
(3)

where $X_{P(D)}^{\mu}$ denote the coordinates of neutrino production
(detection), $\tau$ is the proper time parameterizing the neutrino geodesic
and hence $\tau_{P(D)}$ denote its corresponding value at production
(detection), respectively. To obtain the first order correction $\Delta\gamma$
to the flat spacetime geodesic $\gamma_{0}$, we solve the geodesic equation
for the metric555In this work we use the west coast convention, i.e.
$\eta=\mathrm{diag}(+1,-1,-1,-1)$. $g_{\mu\nu}=\eta_{\mu\nu}+h_{\mu\nu}$ with
the pertubation to flat space time $h_{\mu\nu}$, describing the GW, with
boundary values666Especially we assume that this problem is well posed, i.e.
that geodesics connecting those points are unique. $X_{P}^{\mu}$ and
$X_{D}^{\mu}$. For the step by step calculation see appendix A.

The spatial length $L$ that neutrino wave packets have traveled between
production and detection coordinate time $t_{P}$ and $t_{D}$ is given by

$\displaystyle L(t_{P},t_{D})$
$\displaystyle=\int_{\tau_{P}}^{\tau_{D}}\sqrt{-g_{jk}(\gamma(\tau))\dot{\gamma}^{j}(\tau)\dot{\gamma}^{k}(\tau)}\;\mathrm{d}\tau$
(4) $\displaystyle\approx
L_{0}-\frac{u^{0}}{2}\int_{\tau_{P}}^{\tau_{D}}h_{\parallel}(\gamma_{0}(t))\;\mathrm{d}\tau=:L_{0}+\Delta
L\,,$ (5)

which is then substituted into the standard oscillation probability. Here
$h_{\parallel}$ is the metric perturbation projected onto the zeroth order
neutrino trajectory $\gamma_{0}$. Moreover, we use the ultrarelativistic
approximation for the neutrino velocity and neglect terms of
$\mathcal{O}(h^{2})$ in order to express the detection time as $t_{D}\approx
L_{0}+t_{P}\equiv L_{0}+t$. For the details of the calculation see Appendix B.
In order to evaluate the proper time intergration in equation (5), we use the
fact that any GW can be decomposed into plus and cross polarized plane waves,
i.e.

$\displaystyle h_{\mu\nu}(x)$
$\displaystyle=\sum_{r\in\\{+,\times\\}}\int\mathrm{d}^{3}\vec{k}\;\psi_{r}(\vec{k})A_{\mu\nu}^{r}\cos(\omega
t-\vec{k}\vec{x}+\phi^{r}(\vec{k}))\,,$ (6)

where $A_{\mu\nu}^{r}$ are normalized polarization tensors,
$\psi_{r}(\vec{k})$ are real, positive momentum space wave packets and
$\phi^{r}(\vec{k})$ denote the phase shifts of the corresponding modes.
Exploiting this decomposition, equation (5) becomes

$\displaystyle\Delta L(t)$
$\displaystyle\approx-\frac{1}{2}\sum_{r\in\\{+,\times\\}}\int\mathrm{d}^{3}\vec{k}\;\psi_{r}(\vec{k})\frac{A_{\parallel}^{r}(\varphi,\theta)}{\tilde{\omega}}\left\\{\sin(\tilde{\omega}L_{0})\cos(\omega
t+\phi^{r})+\left[\cos(\tilde{\omega}L_{0})-1\right]\sin(\omega
t+\phi^{r})\right\\}\,,$ (7)

where $\theta$ and $\varphi$ are the angles defining the relative orientation
of the neutrino trajectory $\gamma_{0}$ and the GW propagation direction. In
terms of these angles we define the reduced GW frequency, $\tilde{\omega}$,
and the projected polarization tensors, $A_{\parallel}^{r}(\varphi,\theta)$,
as

$\displaystyle\tilde{\omega}:=\omega(1-\cos(\theta))\,,$ $\displaystyle
A_{\parallel}^{+}(\varphi,\theta):=\sin^{2}(\theta)\cos(2\varphi)\,,$
$\displaystyle
A_{\parallel}^{\times}(\varphi,\theta):=\sin^{2}(\theta)\sin(2\varphi)\,.$ (8)

Depending on the experimental setup, the angles $\theta$ and $\varphi$ could
change with time. This can weaken the overall effect since it is only maximal
for constant perpendicular neutrino and GW directions, i.e. $\theta=\pi/2$,
see Appendix C for a detailed discussion.

Next, we need to average over the duration of the data taking period $T$,
since neither production nor detection time are usually analyzed in neutrino
experiments. Thus, the oscillation probability is given by

$\displaystyle P_{ab}(E,L_{0})$
$\displaystyle=\frac{1}{T}\int_{0}^{T}\hat{P}_{ab}(E,L_{0}+\Delta
L(t))\;\mathrm{d}t\equiv\langle\hat{P}_{ab}(E,L_{0}+\Delta L)\rangle_{T}\,.$
(9)

This gives rise to an additional damping of the oscillating terms in the
oscillation probability $\hat{P}$ which can be interpreted as a new source of
decoherence.

Finally, we need to reconsider the physical meaning of our parameter $L_{0}$.
In the comparison of the prediction of the neutrino oscillation pattern with
and without the GW effect, $L_{0}$ plays the role of the physical distance of
neutrino source and detector measured at a certain time $t_{0}$. If a GW of
any kind is present in the system this will replace the original value of the
parameter $L_{0}$ by $\tilde{L}_{0}=L_{0}+\Delta L(t_{0})$. Therefore
$\tilde{L}_{0}$ is the actual, physical reference length entering the neutrino
flavor transition probability implying that equation (9) needs to be evaluated
at $L_{0}=\tilde{L}_{0}-\Delta L(t_{0})$, corresponding to the subtraction of
a constant off-set. This subtlety is especially important if the GW does not
oscillate over a full period during the time of data taking, $T$, since then
the parameter $L_{0}$ is unphysical and the actual observable length is
$\tilde{L}_{0}$. Hence, only time dependent changes in $L(t)$ can be measured
using the method described above. From now on, we eliminate the unphysical
parameter $L_{0}$ in favor of the physical reference length $\tilde{L}_{0}$
from our description and afterwards rename $\tilde{L}_{0}$ to $L_{0}$.

In order to gain a better understanding of the impact of GWs on the neutrino
oscillation probability, we plot the $\nu_{\mu}\to\nu_{\mu}$ survival
probability in figure 1. For illustration purposes, we choose a baseline of
$L_{0}=${10}^{10}\text{\,}\mathrm{km}$$, neutrino energies on the order of
$\mathcal{O}(10\text{--}100\,\mathrm{MeV})$, standard neutrino mass and mixing
parameters for normal mass ordering [41] and an infinite wave packet width,
i.e. we ignore decoherence from wave packet separation. Furthermore, we
consider a simple, plus polarized, plane GW with variable strain
$\psi_{+}\equiv h$ and frequency
$f\sim$3\text{\times}{10}^{-7}\text{\,}\mathrm{Hz}$$. It is obvious that for
small strains, the modified probability coincides with the standard
expectation while for larger strains, we expect a significant deviation due to
the averaging procedure just described. The non-trivial, oscillating
dependence on the strain can be explained by considering the mechanism of
averaging the flavor transition probability over time. As the strain
increases, so does the amplitude of $\Delta L$ and the averaging proceeds over
many neutrino oscillation cycles resulting in the stable decoherence limit for
large strains. For strains from the intermediate region the time integral only
produces averaging over a partial oscillation cycle leading to an average
value varying with the strain.

Figure 1: Muon neutrino survival probability, $P_{\mu\mu}$, for different
background GWs of different strains $h$. Left: $P_{\mu\mu}$ with (blue) and
without (light blue) the influence of a GW for $E=10\,\mathrm{MeV}$,
$L_{0}=${10}^{10}\text{\,}\mathrm{km}$$. Right: $P_{\mu\mu}$ for different
fixed GW strains $h\in\\{0,10^{-12},10^{-8}\\}$ and variable energy at
$L_{0}=${10}^{10}\text{\,}\mathrm{km}$$. The black oscillation curve ($h=0$)
almost coincides with the blue curve ($h=10^{-12}$) and is therefore hidden.

##  3 General Conditions for Observability of Gravitational Waves Induced
Decoherence

In the following, we discuss the conditions that must be fulfilled for the new
decoherence effect to be observable in experiments:

  1. (i)

Neutrino oscillations proceed sufficiently rapid such that for given GW
parameters the deviation $\Delta L$ is comparable to at least one of the
oscillation lengths.

  2. (ii)

Decoherence due to wave packet separation occurs later than GW induced
decoherence, i.e. at larger baselines or smaller energies.

  3. (iii)

The energy resolution or binning, $\Delta E$, of the neutrino source and
detector are sufficiently accurate such that oscillations are not averaged
out, i.e. neutrino oscillations have to be sufficiently slow for the
experiment to resolve them.

  4. (iv)

The spatial extent of the region of neutrino origin must be much smaller than
the oscillation length of the neutrinos to be resolved.

  5. (v)

Neutrino source and detector must not be fixed such that their physical
distance can be influenced by the GW.

Unfortunately, requirements one and three compete with each other: On the one
hand neutrino oscillations have to be fast enough such that the variation of
the deviation $\Delta L$ over the run time of the experiment is sufficiently
large in order for the time average in equation (9) to damp the oscillations.
On the other hand the oscillations must not be too fast since otherwise the
finite energy bin width will damp the oscillations, and the effect due to the
GW is unresolvable. In the following we further quantify these loosely
formulated criteria.

#### Criterion (i) and (ii)

First we examine the maximal deviation $\Delta L$ from the baseline $L_{0}$.
Since we want to explore the usefulness of neutrino oscillations as probes for
coherent GW signals, we now restrict our discussion to plane GWs777That means
GW signals with sufficiently peaked momentum distributions, i.e.
$\psi_{r}(\vec{k})\approx a_{r}\delta^{(3)}(\vec{k}-\vec{k}_{0})$. In this
case the distance correction $\Delta L$ from equation (7) becomes

$\displaystyle\Delta
L(t)\approx-\frac{1}{2\tilde{\omega}}\sum_{r\in\\{+,\times\\}}a_{r}A_{\parallel}^{r}(\varphi,\theta)\left\\{\sin(\tilde{\omega}L_{0})\cos(\omega
t+\phi^{r})+\left[\cos(\tilde{\omega}L_{0})-1\right]\sin(\omega
t+\phi^{r})\right\\}\,,$ (10)

with mean frequency $\omega$ and polarization amplitudes $a_{r}$. It is
evident that for a fixed length $L_{0}$ the amplitudes appearing in this
expression oscillate with respect to the reduced frequency $\tilde{\omega}$.
This makes it more difficult to set bounds on the high frequency region
$\tilde{\omega}\geq L_{0}^{-1}$ which is why, for the following discussion, we
focus on the low frequency region $\tilde{\omega}\ll L_{0}^{-1}$. In this
region, we can express the maximal baseline deviation as

$\displaystyle\max_{\omega t\in[0,2\pi)}|\Delta L(t)|=\frac{hL_{0}}{2}\,,$
(11)

with

$\displaystyle h:\approx\sin^{2}(\theta)\max_{\omega
t\in[0,2\pi)}|a_{+}\cos(2\varphi)\cos(\omega
t+\phi^{+})+a_{\times}\sin(2\varphi)\cos(\omega t+\phi^{\times})|\,.$ (12)

For a discernible effect, it is also crucial that the maximum deviation
$\Delta L$ is comparable to the oscillation length $L^{\text{osc}}_{jk}$ of
neutrinos:

$r=\frac{\max_{t\in\tau}|\Delta L(t)|}{L^{\text{osc}}_{jk}}=\frac{\Delta
m_{jk}^{2}}{4\pi E}\frac{hL_{0}}{2}\,,$ (13)

where we conservatively choose the ratio $r\approx 0.1$. The reason for this
is that the $\nu_{j}$–$\nu_{k}$ interference term needs to be sufficiently
damped by the averaging procedure. For each scenario we consider the
oscillation term with the largest mass difference $\Delta m_{jk}^{2}$, since
we are interested in the sensitivity of our framework. However, the decision
on which oscillation dominates must be made based on the experiment under
consideration. From this we can derive a lower limit on the baseline
$L_{\text{min}}$ such that the effect becomes visible, i.e.

$L_{\text{min}}(E)=8\pi r\frac{E}{h\Delta m_{jk}^{2}}<L_{0}\,.$ (14)

On the other hand, it is essential for the baseline $L$ to be smaller than the
coherence length $L_{jk}^{\text{coh}}$ of neutrinos to ensure that the effect
of wave packet decoherence does not dominate and thus hide the desired effect.
This results in an upper limit implying that the baseline $L_{0}$ has to vary
within the following window to enable the observation of an effect:

$\displaystyle L_{\text{min}}(E)$ $\displaystyle<L_{0}<L^{\text{coh}}_{jk}$
(15) $\displaystyle\Leftrightarrow 8\pi r\frac{E}{h\Delta m_{jk}^{2}}$
$\displaystyle<L_{0}<\frac{4\sqrt{2}E^{2}}{\Delta m_{jk}^{2}}\sigma_{x}\,.$
(16)

#### Criterion (iii) and (iv)

The energy resolution $\Delta E$ must be sufficiently high to resolve the
neutrino oscillations by the experiment. This implies that the frequency of
the oscillation cannot be too large as too many wave cycles are contained
within an energy bin, thus averaging out the effect. In the case of a
sinusoidal oscillation of the form $\sin{\left(\omega t\right)}$, this would
mean that the period $\Delta t$, corresponding to the bin width, should be
smaller than the inverse frequency: $\Delta t\lesssim\omega^{-1}$. This
consideration can be applied to our case, leading to an energy resolution
estimate for the resolution of the $jk$-oscillation:

$\displaystyle\left|\Delta\left(\frac{L_{0}}{E}\right)\right|\lesssim\frac{2}{\Delta
m_{jk}^{2}}\quad\Leftrightarrow\quad\Delta E\lesssim\frac{2}{\Delta
m_{jk}^{2}}\frac{E^{2}}{L_{0}}\,,$ (17)

In addition to energy resolution, position resolution is also relevant. It is
important that we know the origin of the neutrinos to a certain degree, hence
the width of the production region $\Sigma_{P}$ must be smaller than the
oscillation length $L^{\text{osc}}_{jk}$ to prevent additional averaging
hiding the gravitationally induced effect:

$\Sigma_{P}\ll L^{\text{osc}}_{jk}\,.$ (18)

#### Criterion (v)

In the measurement of GW effects, it is crucial to account for the absence of
strong forces maintaining a constant physical distance between two reference
points. For instance, in the case of LIGO, this is achieved through freely
suspended mirrors [42]. In our context, this demands the forces between the
detector and the source to be negligible compared to the force of displacement
induced by the GW. Astrophysical neutrino sources do not pose a problem in
this regard, as detector and source are freely falling in the local
gravitational fields. One possibility would be to detect neutrinos from
natural sources like the sun or, assuming they exist, decaying dark matter
(DM) particles within the gravitational wells of celestial bodies. At even
larger distances pulsars, blazars or supernovae are possible candidates for
neutrino sources.

#### Selection of Benchmark Scenarios

We can finally use these conditions to find interesting parameter
configurations for baselines $L_{0}$, the mass squared difference $\Delta
m_{jk}^{2}$, the neutrino wave packet width $\sigma_{x}$ and energy $E$. For
SM neutrinos we know the mass squared differences to be $\Delta
m_{21}^{2}=$7\text{\times}{10}^{-5}\text{\,}{\mathrm{eV}}^{2}$$ [41] and
$\Delta m_{32}^{2}=$2.437\text{\times}{10}^{-3}\text{\,}{\mathrm{eV}}^{2}$$
[41], respectively, while typical wave packet widths are in the $\mathrm{nm}$
range or below [43]. From the condition $\tilde{\omega}>L_{0}^{-1}$, cf.
criterion (i), we can now infer interesting baseline configurations depending
on which GW frequencies we want to be sensitive to.

We are particularly interested in frequency ranges of GWs that are not yet
explored by other existing or planned experiments. EPTA, IPTA and NANOgrav can
cover frequency ranges from roughly ${10}^{-8}\text{\,}\mathrm{Hz}$ to
${10}^{-6}\text{\,}\mathrm{Hz}$ [44, 45] and LIGO is sensitive in the range of
$10\text{\,}\mathrm{Hz}$ to ${10}^{4}\text{\,}\mathrm{Hz}$ [46]. The planned
experiment LISA and proposed eLISA would be sensitive to frequencies from
${10}^{-5}\text{\,}\mathrm{Hz}$ to $1\text{\,}\mathrm{Hz}$ [10]. From this,
three frequency ranges with $f=\omega/(2\pi)$ can be identified that are not
yet covered by currently running or approved experiments:

$\displaystyle f_{\text{GW},1}$
$\displaystyle\lesssim${10}^{-9}\text{\,}\mathrm{Hz}$\,,$ (19)
$\displaystyle${10}^{-6}\text{\,}\mathrm{Hz}$\lesssim f_{\text{GW},2}$
$\displaystyle\lesssim${10}^{-5}\text{\,}\mathrm{Hz}$\,,$ (20)
$\displaystyle$1\text{\,}\mathrm{Hz}$\lesssim f_{\text{GW},3}$
$\displaystyle\lesssim$10\text{\,}\mathrm{Hz}$\,.$ (21)

In order to detect an effect, the baseline $L_{0}$ should be smaller than
$f^{-1}$. Therefore, we suggest three baselines $L_{0}$:

$\displaystyle L_{0,1}$ $\displaystyle=${10}^{5}\text{\,}\mathrm{km}$\,,$ (22)
$\displaystyle L_{0,2}$ $\displaystyle=${10}^{8}\text{\,}\mathrm{km}$\,,$ (23)
$\displaystyle L_{0,3}$ $\displaystyle=${10}^{10}\text{\,}\mathrm{km}$\,.$
(24)

All three values are comparable to distances within our solar system, leading
us to the experimental challenge discussed above, cf. criterion (v). An
additional difficulty arises from the duration of coherent GW signals
$T_{\text{GW}}$ compared to the data taking period of the neutrino experiment
$T_{\text{exp}}$. If $T_{\text{GW}}\ll T_{\text{exp}}$ the averaging effect
will be washed out and the sensitivity decreases. Therefore we expect neutrino
based GW detection to be most sensitive to astrophysical events emitting GW
signals of durations of $\mathcal{O}\left($10$\,\text{yr}\right)$ or greater,
like supermassive blackhole binary (SMBHB) [47] mergers or extreme-mass-ratio
inspirals (EMRI) [48]. It is worth noting that SMBHB mergers typically last
longer than 10 years while emitting GWs with a constant frequency [49].

Next we want to estimate the sensitivity to the strain of a given GW signal
for the parameter configurations discussed above. To this end we show the
energy and strain ranges for which we expect an observable effect at the three
baselines from equation (22)–(24) of interest in figure 2.

Figure 2: Sensitivity Range determined from the requirement that neutrino
oscillations proceed sufficiently rapid (criterion (i)) for fixed baselines
$L_{0}$ with neutrino mass differences $\Delta m_{jk}^{2}$ consistent with the
neutrino oscillation data and neutrino wave packet widths in the $\mathrm{nm}$
range. The red area marks the range where neutrino oscillations can be
sensitive on GWs in each case, while the blue line marks the upper limit on
the energy and the orange line the lower limit as a function of the strain
$h$.

Here, the orange line represents the lower limit on the energy $E$ as a
function of the strain $h$ for a wave packet width of
$\sigma_{x}=$1\text{\,}\mathrm{nm}$$ fulfilling the condition in equation
(16). The blue line shows the upper limit on the energy $E$ as a function of
the strain $h$ for the respective baselines and a mass splitting $\Delta
m_{jk}^{2}=${10}^{-3}\text{\,}{\mathrm{eV}}^{2}$$, which results from the
condition in equation (14). The intersection of the two lines marks the lowest
strain $h_{\text{min}}$ that is observable for the scenarios just discussed.
To clarify which ranges of strain can realistically be accessed,
$h_{\text{min}}$ and the corresponding energy $E$ are shown as a function of
the baseline $L_{0}$ in figure 3.

Figure 3: Estimation of the minimally observable strain $h_{\text{min}}$ in
the case of SM neutrinos with a neutrino baseline $L_{0}$, energy $E$ and
neutrino wave packet width $\sigma_{x}=$1\text{\,}\mathrm{nm}$$.

It is evident that the bounds for the strain are several orders of magnitude
weaker compared to LIGO and LISA, thus favoring very large length $L_{0}$ and
energy $E$ combinations. As a consequence, this method is particularly
sensitive to SMBHB and EMRI below the frequency range probed by LISA. In the
following, we will study more closely how criteria (i)–(v) can be met in
concrete scenarios for the potentially interesting baselines $L_{0,1}$,
$L_{0,2}$ and $L_{0,3}$ identified above and perform a precise statistical
analysis demonstrating that our estimations for the minimum strains
$h_{\text{min}}$ are indeed accurate.

##  4 Parameter Constraints for Gravitational Waves: A Toy Analysis

In order to estimate the actual sensitivity of neutrino oscillations on the
presence of a GW, we perform a toy analysis for each scenario discussed in the
preceding section taking into account all deviations from the standard
oscillation pattern. To this end, we model several hypothetical experimental
set-ups and estimate the strain and frequency regions they would be sensitive
to. Each of these hypothetical experiments is assumed to count neutrino event
rates in finite, equal width energy bins and the corresponding neutrino signal
is expected to originate from a source at some distance $L_{0}$. We discuss
two different classes of sources. First, we assume a peaked neutrino energy
spectrum, like for example found for neutrinos from fixed target experiments
or from decaying particles at rest, approximated by a Gaussian energy spectrum
with width $\sigma_{E}$ centered at some mean energy $\bar{E}$. Second, we
assume a power law spectrum for neutrinos from astrophysical environments
where charged particles are accelerated in jets and interact with each other
or the environment to produce neutrinos.

While the power law dependence is a rather accurate model of a natural
neutrino flux, the Gaussian flux is more like a crude approximation to flux
shapes exhibiting a peak at a certain energy, $\bar{E}$. We still adopt this
approximation for simplicity because the sensitivity to the GW parameters does
not strongly depend on the shape of the flux spectrum as long as there are
sufficiently high event counts in the chosen energy range
$E\in[E_{\mathrm{min}},E_{\mathrm{max}}]$. For the same reason we generically
choose $\bar{E}=(E_{\mathrm{max}}+E_{\mathrm{min}})/2$ and
$\sigma_{E}=0.1(E_{\mathrm{max}}-E_{\mathrm{min}})$. The number of neutrino
events is then predicted according to the probability density,

$\displaystyle\rho_{b}(E)$
$\displaystyle=\sum_{a}\varphi_{a}(E)P_{ab}(E,L_{0})\,,$ (25)

where $\varphi_{a}$ is the neutrino flavor and energy distribution at the
source and $P_{ab}$ is the flavor transition probability of the given model.

We generate a toy data set, assuming the expected number of events for each
energy and flavor bin to follow the flat spacetime prediction with fixed
length $L_{0}$. On this data set, we perform a likelihood ratio test where we
evaluate the negative logarithmic likelihood ratio (NLLR) on each point of a
region in our parameter space. We only consider “plus” $(+)$ polarized GWs
perpendicular to the neutrino trajectory since both possible polarizations are
expected to cause a similar effect. Furthermore, we assume a data taking
period of $T_{\mathrm{exp}}=20\,\mathrm{yr}$ for each scenario for the details
of the time average, cf. equation(9). For GW frequencies $f\gtrsim
T_{\mathrm{exp}}^{-1}\approx$1.6\text{\,}\mathrm{nHz}$$ the effect of GW
induced decoherence is insensitive to the phase $\phi^{+}$ of the GW since we
are averaging over multiple cycles of the wave period. Hence, we set it to
zero if not stated otherwise.

At frequencies below $T_{\mathrm{exp}}^{-1}$ the decoherence effect caused by
the GW weakens since $\Delta L$ only changes slightly during the data taking
period. Thus, we mostly restrict the frequency range to values above this
threshold. It is only for neutrinos from galactic pulsars where we also
consider lower frequencies and discuss the implications of the slow variation
of $\Delta L$.

The remaining parameter space is two dimensional and consists of the possible
strain, $h$, and frequency values, $f$. Further details of the setup and
statistical analysis can be found in appendix D. The precise statistical
analysis for SM neutrinos with baselines $L_{0}$ corresponding to distances in
the solar system can be found in section 4.1. The analysis of astrophysical
neutrino sources beyond the solar system can be found in section 4.2. Finally,
we consider the influence of a SGWB in section 4.3.

###  4.1 Neutrinos from within the Solar System

To gain an intuition for the physics involved, we first discuss the simplest
yet rather unrealistic scenario of either free falling artificial neutrino
sources placed, for example, in a stable Lagrange point within the solar
system or neutrinos orignating from dark matter decaying or annihilating
inside the gravitational potential wells of celestial bodies like the moon.
While we will find that neutrinos from farther sources turn out to be more
promising for GW detection, the following discussion will help us to identify
the challenges to obtain an observable effect. We do not consider solar
neutrinos where oscillations are already impeded by strong matter effects.
Therefore, we assume an initially pure muon neutrino flux modelled by a
Gaussian shaped energy spectrum with a total event count of $N=${10}^{5}$$ for
an experimental running time of 20 years. We operate within realistically
achievable parameter configurations, selecting a neutrino wave packet width in
the $\mathrm{nm}$ range [50], and conduct a statistical analysis using the
parameters in table 1 to verify the estimations from section 3. This enables
us to determine the strains and frequencies of GWs to which corresponding
experiments would be sensitive to.

Table 1: Wavelengths $\lambda_{\text{GW}}$, baselines $L_{0}$, energy ranges $E$ and bin widths $\Delta E$, where effects would be observable for the three chosen GW frequency ranges $f_{\text{GW}}$ and neutrino wave packet widths of $\sigma_{x}=$1\text{\,}\mathrm{nm}$$. The oscillation lengths for the different energies are given in table 2. Scenario | $f_{\text{GW}}\,/\,$\mathrm{Hz}$$ | $\lambda_{\text{GW}}\,/\,$\mathrm{km}$$ | $L_{0}\,/\,$\mathrm{km}$$ | $E$ | $\Delta E$
---|---|---|---|---|---
1 | $1$–$10$ | ${10}^{7}$–${10}^{8}$ | $\thicksim${10}^{5}$$ | $[1,10]$\mathrm{keV}$$ | $\thicksim$0.45\text{\,}\mathrm{keV}$$
2 | $[1,10]$\mathrm{MeV}$$ | $\thicksim$0.45\text{\,}\mathrm{MeV}$$
3 | $[1,10]$\mathrm{GeV}$$ | $\thicksim$0.45\text{\,}\mathrm{GeV}$$
4 | ${10}^{-6}$–${10}^{-5}$ | ${10}^{10}$–${10}^{11}$ | $\thicksim${10}^{8}$$ | $[1,10]$\mathrm{keV}$$ | $\thicksim$0.45\text{\,}\mathrm{keV}$$
5 | $[1,10]$\mathrm{MeV}$$ | $\thicksim$0.45\text{\,}\mathrm{MeV}$$
6 | $[1,10]$\mathrm{GeV}$$ | $\thicksim$0.45\text{\,}\mathrm{GeV}$$
7 | ${10}^{-9}$ | ${10}^{17}$ | $\thicksim${10}^{10}$$ | $[1,10]$\mathrm{keV}$$ | $\thicksim$0.45\text{\,}\mathrm{keV}$$
8 | $[1,10]$\mathrm{MeV}$$ | $\thicksim$0.45\text{\,}\mathrm{MeV}$$
9 | $[1,10]$\mathrm{GeV}$$ | $\thicksim$0.45\text{\,}\mathrm{GeV}$$

Figure 4 depicts the results of our analysis for the parameters from table 1
with sensitivity curves for the experiments taken out of reference [51]. LIGO
and LIGOa give limits at 90 % CL [52] as well as Virgo [52], and LISA [53].
EPTA, IPTA, and eLISA give limits at 95 % CL [54, 55, 56]. In figure 4 we show
the upper limit of the strain $h$ as a function of the frequency $f$ at
$95\,\%$ CL assuming the absence of the GW induced decoherence in the
experimental data set.

Figure 4: Upper bounds on the GW strain $h$ as a function of the frequency $f$
at 95 % confidence level and wave packet widths
$\sigma_{x}=$1\text{\,}\mathrm{nm}$$ for the scenarios from table 1. We also
show the sensitivity curves for EPTA [54], IPTA [55], eLISA [56], LISA [53],
LIGO [52], LIGOa [52] and Virgo [52].

It is evident that only three parameter configurations corresponding to the
scenarios 2, 6, and 9 from table 1 yield meaningful constraints on the strain
$h$ and frequency $f$ of the GW. This is because the condition from equation
(16) is not satisfied for the remaining six parameter configurations, and
thus, the GW decoherence effect is obfuscated by wavepacket decoherence.
Another significant factor at play is the energy binning, which according to
equation (17), for example, should be in the $\mathrm{meV}$ range for the
baseline $L_{0,2}$, which is way smaller than what experiments can resolve.

The sensitivity curves in figure 4 can be divided into two regions: a plateau
occuring at frequencies $f\ll L_{0}^{-1}$ and an oscillatory region at higher
frequencies $f\gtrsim L_{0}^{-1}$. This behavior is understandable when
considering equation (10). In the oscillatory region, values of
$\tilde{\omega}L_{0}=2\pi k$ with $k\in\mathbb{Z}$ can occur, resulting in
$\sin{\left(\tilde{\omega}L_{0}\right)}=\cos{\left(\tilde{\omega}L_{0}\right)}-1=0$,
and overall no deviation $\Delta L=0$ from the length $L_{0}$ occurs. This
behavior can be explained by the fact that at these frequencies, an equal
number of wave troughs and wave crests influence the neutrino propagation,
leading to no net effect. Hence, this region is not suitable for investigating
the parameters of the GW. The most interesting region is the plateau where
$fL_{0}\ll 1$ holds. Here,
$\sin{\left(\tilde{\omega}L_{0}\right)}/\tilde{\omega}L_{0}$ approaches one,
resulting in a frequency independent amplitude of $\Delta L(t)$ shaping this
plateau. Hence, the frequency at which the contour enters the plateau region
depends on the baseline $L_{0}$. Furthermore, at the lowest end of the
frequency region one can see that the contours transition from a plateau into
an oscillating (intermediate) phase. This oscillation is due to the decreasing
number of GW cycles $n:=\lfloor fT_{\mathrm{exp}}\rfloor$ over which we
average the modified flavor transition probability in equation (9). As
$n\rightarrow\infty$, i.e. $f\gg T_{\mathrm{exp}}^{-1}$, the contribution of
the remaining incomplete period $T_{\mathrm{exp}}-nf$ is negligible compared
to the $n$ complete periods. Therefore, the time average corresponds to the
average over a single full period. For $n=\mathcal{O}(1)$, i.e. $f\rightarrow
T_{\mathrm{exp}}^{-1}$, the suppression of the incomplete cycle becomes weaker
resulting in periodically changing weaker or stronger deviations from the flat
spacetime predictions, ultimately resulting in the oscillating contour at the
lowest frequencies.

According to figure 4, it is possible to achieve constraints on the strain $h$
as low as $h\approx$3.2\text{\times}{10}^{-7}$$ with scenario 9 from table 1,
i.e. SM neutrinos with realistic wave packet widths, and baselines smaller
than the diameter of the solar system. These results confirm our prediction
from figure 3. Despite these rather weak constraints on the strain $h$, it is
interesting to note that with this method so far unexplored frequency ranges
can be probed. An attempt to improve on these sensitivities would require
larger wave packet widths $\sigma_{x}$, an approach discussed in Appendix E.
Nevertheless, the consideration of neutrinos originating from distances within
the solar system poses several challenges. For instance, deploying a detector
or the source in a Lagrange point or launching them into the interstellar
medium is technologically challenging, to put it mildly. Moreover, even if
dark matter exists inside of celestial bodies and decays or annihilates into
neutrinos, the corresponding neutrino production regions would be rather large
and hence positional averaging would already average out the oscillation
pattern and thereby obfuscating the GW induced decoherence effect, c.f.
criterion (iv) in section 3. Furthermore, matter effects for the propagation
out of the core of the celestial bodies have been ignored and finally
reasonable wave packet widths yield such weak constraints that this approach
is not promising to be investigated further. An alternative option would be a
scenario where a hypothetical heavy, sterile neutrino flavor mixes with the
active neutrinos. As a consequence of larger mass splitting $\Delta
m_{4k}^{2}$ might improve the sensitivity towards smaller strains. We
investigate these models in Appendix F, there we also discuss further
scenarios. However, due to strong, mass dependent constraints on the sterile-
active mixing angles that suppress this oscillation channel, this approach
does not provide an advantage over the light neutrino case. Therefore, we turn
to the more promising scenario considering neutrinos from outside the solar
system.

###  4.2 Neutrinos from Galactic Pulsars

After pointing out that the prospects to probe GWs at distances within our
solar system are rather bleak, we now consider astrophysical neutrino sources
outside the solar system. We assume a neutrino flux equivalent to a power-law
spectrum $\varphi\propto E^{-2}$ and an initial
$(\nu_{e},\nu_{\mu},\nu_{\tau})=(1,2,0)$ flavor composition as the commonly
considered scenario [58, 57]. We note however, that also different initial
flavor ratios are possible depending on the properties and surroundings of the
source [59, 60, 61, 62, 63, 64, 65].

At these length scales, one challenge lies in the requirement that the
neutrino production region needs to be sufficiently small compared to the
oscillation length $L^{\text{osc}}_{jk}$, cf. equation (18), in order to
detect the decoherence effect caused by GWs. The ratio
$L^{\text{osc}}_{jk}/\Sigma_{P}$ scales linearly with the energy $E$.
Consequently, the higher the energy of the neutrinos, the larger the neutrino
production region can be without significantly averaging out the oscillations.
Table 2 gives the oscillation lengths, $L_{jk}^{\mathrm{osc}}$, of neutrinos
at various energies, $E$, and thereby provides an intuition about how accurate
the region of origin would need to be known.

Table 2: Values for the oscillation length $L^{\text{osc}}_{jk}$ for different neutrino energies $E$ and mass splitting of $\Delta m_{jk}^{2}\sim${10}^{-3}\text{\,}{\mathrm{eV}}^{2}$$. $E$ | $L^{\text{osc}}_{jk}\,/\,$\mathrm{km}$$
---|---
$\mathrm{keV}$ | $2.4\text{\times}{10}^{-3}$
$\mathrm{MeV}$ | $2.4$
$\mathrm{GeV}$ | $2.4\text{\times}{10}^{3}$
$\mathrm{TeV}$ | $2.4\text{\times}{10}^{6}$
$\mathrm{PeV}$ | $2.4\text{\times}{10}^{9}$
$\mathrm{EeV}$ | $2.4\text{\times}{10}^{12}$

A possible astrophysical source for neutrinos are supernovae [66], typically
emitting $\mathrm{MeV}$ neutrinos [67]. However, since we would need to know
the location of the neutrino source within a $\mathrm{km}$ range to identify
decoherence from GWs no sensitivity is expected for such sources, see table 2.

Another astrophysical source of neutrinos are blazars [68], i.e. active
galactic nuclei (AGNs) with jets pointing in our direction. Neutrinos from
such a source are expected to be in the energy range from $\mathrm{PeV}$ to
$\mathrm{EeV}$ [69]. According to table 2, for such high-energy neutrinos the
spatial size of the neutrino source can be as large as approximately
${10}^{14}\text{\,}\mathrm{km}$. While this may be possible for neutrinos
produced in the vicinity of a blazar, these sources suffer from an expected
low flux of high energy neutrinos reaching the Earth, resulting in a non-
significant neutrino count, cf. reference [70], compared to the required event
count of $N=${10}^{4}$$–${10}^{5}$ in our set-up.

Thus, the most promising astrophysical sources of neutrinos are the jets of
pulsars [71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81] from our own Galaxy, like
for example the Vela ($L_{0}=294^{+76}_{-50}\,\mathrm{pc}$ [82]), Crab
($L_{0}=1.9^{+0.22}_{-0.18}\,\mathrm{kpc}$ [83]) or Cas A ($L_{0}=3.5\pm
0.2\,\mathrm{kpc}$ [84]) pulsars. Due to the relatively high expected neutrino
fluxes of about $\mathcal{O}(30\,\mathrm{km}^{-2}\mathrm{yr}^{-1})$, event
counts of $N\sim${10}^{4}$$ seem to be realistic at future
$100\,\mathrm{km}^{2}$-neutrino telescopes and a data taking period of about
$T_{\mathrm{exp}}\sim 20\,\mathrm{yr}$.

Very high energy neutrinos emerging from pulsars may be generated in two
regions: The vicinity of the pulsar’s magnetosphere or in the surrounding
pulsar wind nebula. In both cases protons are accelerated to high energies and
interact with photons or other protons [71, 72, 73, 74, 75, 76, 77, 78, 79,
80, 81, 85, 86] producing pions which in turn decay into neutrinos. So far it
is unkown, which of these regions yields the dominant neutrino flux. While the
size of the magnetosphere is of the order of several $\mathrm{km}$, depending
on the specific properties of the pulsar [85], and thus fulfills the criterion
$\Sigma_{P}\ll L^{\mathrm{osc}}_{jk}$, the pulsar wind nebula is much larger.
It is, however, uncertain how large the actual size of the neutrino production
region within the nebula really is. In the following, we assume that the width
of the production region is sufficiently small compared to the oscillation
length, cf. table 2. Even if this assumption does not apply, pulsars remain
suitable sources when a potential additional oscillation damping originating
from a necessary averaging over the neutrino production region is included in
the analysis. Furthermore, we assume neutrino energies in the energy range
$E\in\left[$100\text{\,}\mathrm{TeV}$,$2\text{\,}\mathrm{PeV}$\right]$, as
predicted e.g. in reference [80], and a linear energy binning of $\Delta
E=$95\text{\,}\mathrm{TeV}$$. With this binning, we conservatively mimic the
logarithmic choice of binning used in a recent IceCube analysis of the
astrophysical neutrino flux, cf. reference [87]. In the following, we
concentrate on neutrinos from the well studied Vela Pulsar, located in the
Vela supernova remnant.

Figure 5: As figure 4, for neutrinos from the Vela Pulsar for different GW
phases $\phi^{+}\in\\{0,\pi/4,\pi/2\\}$ and wave packet width
$\sigma_{x}=$1\text{\,}\mathrm{nm}$$.

Figure 5 depicts the the upper limit on the strain $h$ as a function of the
frequency $f$ at $95\,\%$ CL resulting from a statistical analysis for the
mentioned parameter configurations. As discussed in the last section, we
expect the best sensitivity at $f\ll L_{0}^{-1}$. At astrophysical distances,
like $L_{0}\sim 1000\,\mathrm{lyr}$, a new issue arises: In order to achieve a
maximum oscillation damping in the time averaged oscillation probability, the
averaging interval needs to be of the order of the inverse frequency of the
GW, i.e. $T_{\mathrm{exp}}\gtrsim T_{\mathrm{GW}}:=f^{-1}$. Since we are
mostly interested in $T_{\mathrm{GW}}\gg L_{0}$ to eliminate a phase
dependence completely, we would also need $T_{\mathrm{exp}}\gtrsim
T_{\mathrm{GW}}\gg L_{0}$, corresponding to a $1000\,\mathrm{yr}$ period of
data taking or more. This is of course impossible to achieve. Nevertheless, a
decoherence effect may still be observed given that the neutrino oscillation
length is much smaller than the period of the GW, see table 2, albeit with a
smaller sensitivity that results from the fact that only a small portion of
the GW strain effects the experiment. As a consequence, the contours do not
enter the plateau phase as they would for $T_{\exp}\gg L_{0}$. Instead the
physical length888As discussed in section 2 we need to renormalize $L(t)$ with
respect to some time $t_{0}$ which we take to be zero for simplicity.,

$\displaystyle L(t)$ $\displaystyle=L_{0}+\Delta L(t)-\Delta
L(0)\stackrel{{\scriptstyle T_{\mathrm{GW}}\gg L_{0}\gg
T_{\mathrm{exp}}}}{{\approx}}L_{0}\left(1+\frac{h}{2}\sin(\phi^{+})(2\pi
ft)+\mathcal{O}\left(\left(2\pi ft\right)^{2}\right)\right)\,,$ (26)

approaches the same value, $L_{0}$, as the one used in the flat spacetime
prediction as $f\to 0$. Obviously, due to the lack of averaging over multiple
periods of the GW, this method becomes dependent on the initial phase of the
GW, $\phi^{+}$. This behavior can indeed be seen in the contours shown in
figure 5 for 20 years of data taking where those upper limits on $h$ for phase
shifts closer to $\pi/2$ are stronger than those corresponding to phases
closer to zero, cf. equation (26).

Thus, the consideration of neutrinos from pulsars enables us in principle to
explore frequency domains below $f\approx${10}^{-12}\text{\,}\mathrm{Hz}$$,
where we expect signals from SMBHBs and EMRIs. However, for SMBHBs, only
strains of $h\lesssim 10^{-15}$ [88] and for EMRIs $h\lesssim 10^{-21}$ [51]
are expected, which is significantly below the range of $h\approx${10}^{-7}$$
that can be probed realistically. Besides these expected sources of GW signals
in the $\mathrm{nHz}$ band neutrino oscillation driven GW detectors would in
principle also be sensitive to more exotic models or so-far unknown effects
that produce GWs at even lower frequencies below $f_{\mathrm{GW}}\lesssim
10^{-13}\,\mathrm{Hz}$. Therefore, this method constitutes a complementary
approach to the study of anisotropies of the cosmic microwave background [19,
89, 90] in order to probe such low frequency GW signals.

Moreover, at least for approximately plane GW signals, the method described
above now depends on the phase of the wave and results in a weaker sensitivity
for experimental running times smaller than the equivalent distance between
neutrino source and detector. The only way to improve this weak sensitivity is
either by taking neutrino data over an unrealistically long time period or by
employing various measurements over a very long period.

A more realistic way out of this dilemma exists for another source, the SGWB,
a random superposition of low frequency plane wave signals of different phases
and directions. Since each neutrino trajectory is perturbed by a slightly
different superposition of such waves the effect might become visible even if
each signal taken on its own would not lead to a visible effect in the
oscillation pattern. We study this effect in more detail in the next section.

A further, common problem with neutrinos from astrophysical sources is the
accurate determination of the distance $L_{0}$ which is subject to rather
large uncertainties compared to the baseline variation caused by the GW. This
might also impair the sensitivity to the influence of GW induced decoherence.

###  4.3 Probing the Stochastic Gravitational Wave Background with Standard
Model Neutrinos from Galactic Pulsars

Until now, we considered the interplay of neutrino oscillations and plane GW
signals of well defined origin. Another interesting phenomenon can be studied
in the $\mathrm{nHz}$ frequency band, i.e. the SGWB. In parallel to the recent
discovery of the SGWB by PTA collaborations, already several works about
possible consequences for neutrino oscillations have appeared [16, 19]. While
in these studies the interplay of GW induced decoherence and other sources of
oscillation damping was already discussed qualitatively, we now adopt a more
quantitative approach by analyzing the scenario considered in reference [19]
statistically.

To this end, we revisit equation (1) and incorporate not only the wave packet
separation effect into the decoherence term, but also the damping term induced
by the SGWB,

$\mathcal{D}_{jk}(E,L_{0})=\left(\frac{L_{0}}{L^{\text{coh}}_{jk}}\right)^{2}+\Gamma_{jk}^{\text{SGWB}}(E)\,.$
(27)

This approach, while being different to the one employed so far, roots in very
similar considerations of the modification of the neutrino propagation due to
the GW in the density matrix formalism, cf. reference [16] for further
details. For the following analysis, we do not average over the running time
of the neutrino experiment anymore as we have shown that this effect is
negligible at the low GW frequencies of the SGWB. Therefore, we only take into
account the stochastic perturbations of the neutrino trajectories caused by an
ensemble of GWs through which the neutrino system propagates.

The final expression for $\Gamma_{jk}^{\text{SGWB}}$ according to reference
[19] is given by

$\Gamma_{jk}^{\text{SGWB}}=\frac{3}{64\left(\gamma-1\right)}\left(\frac{|A_{*}|}{f_{\text{yr}}L^{\text{osc}}_{jk}}\right)^{2}\left(\frac{f_{\text{min}}}{f_{\text{yr}}}\right)^{1-\gamma}\,,$
(28)

where the dependence of the SGWB strain $h_{c}$ on the frequency $f$ is
modelled via a power-law spectrum of the form

$h_{c}(f)=A_{*}\left(\frac{f}{f_{\text{yr}}}\right)^{\frac{3-\gamma}{2}}\,.$
(29)

Here $f_{\text{yr}}=1\text{yr}^{-1}\approx$31.7\text{\,}\mathrm{nHz}$$ is the
reference frequency, $\gamma$ is the spectral index, $f_{\text{min}}$ is the
lowest frequency of the GW spectrum and $A_{*}$ is the amplitude.

We concentrate on the case of an unspecified source of the SGWB, meaning that
$A_{\ast}$ as well as the spectral index $\gamma$ are free parameters for
several values of the minimum frequency $f_{\mathrm{min}}$. In the full
framework accounting for neutrino wave packet separation as well as finite
energy binning, we expect the best results to occur for very high energy
neutrinos. This corresponds to those scenarios considered in reference [19]
where $f_{\mathrm{min}}\in[10^{-18},10^{-15}]\,$\mathrm{Hz}$$. For our
analysis, we again assume the same setting as in section 4.2, i.e. neutrinos
emmited from the Vela Pulsar in the energy range
$E\in[100\,\mathrm{TeV},2\,\mathrm{PeV}]$ with a power law flux spectrum
$\varphi(E)\propto E^{-2}$. Furthermore, we assume a pure $\nu_{\mathrm{e}}$
flux as in reference [19] and wave packet widths of
$\sigma_{x}\sim$1\text{\,}\mathrm{nm}$$.

Figure 6: Shown are the limits at 95 % CL on the GW ampltiude $A_{\ast}$ and
the spectral index $\gamma$ for different values of minimum frequencies,
$f_{\text{min}}\in\\{10^{-18},10^{-17},10^{-16},10^{-15}\\}$. Parameter
configurations to the upper right of these lines are excluded if the effect of
the SGWB is not observed in neutrino oscillation experiments at baselines
$L_{0}\sim 1000\,\mathrm{lyr}$, energies
$E\in[100\,\mathrm{TeV},2\,\mathrm{PeV}]$, bin widths $\Delta E\lesssim
90\mathrm{TeV}$ and neutrino wave packet widths $\sigma_{x}\sim
1\,\mathrm{nm}$.

The results of our analysis are visualized in figure 6. Shown are the $95\,\%$
confidence level contours of the NLLR function in the $\gamma$–$A_{\ast}$
plane for different values of minimum frequencies, $f_{\text{min}}$,
respectively. The parameter regions to the upper right of these contours are
excluded for a data set generated assuming the absence of a GW. For the
corresponding configurations, i.e. higher strain and bigger spectral index
than indicated by the contour, neutrino oscillations are damped compared to
the prediction in flat spacetime. Therefore, our results are in very good
agreement with those from reference [19], demonstrating that neutrinos in the
$\mathrm{PeV}$ energy range can indeed provide new insights to properties of
the SGWB, even when taking wave packet decoherence and finite energy
resolution of neutrino experiments into account.

For much larger cut-off frequencies like
$f_{\mathrm{min}}\in[10^{-12},10^{-9}]\,$\mathrm{Hz}$$ on the other hand,
neutrinos with energies in the $\mathrm{MeV}$ range would be needed in order
to observe a significant effect from the SGWB. A possible source for
astrophysical neutrinos of this energy would be supernova neutrinos where
typical wave packet widths are $\sigma_{x}\sim$0.1\text{\,}\mathrm{pm}$$ [38].
In this case, as discussed in reference [19], decoherence from wave packet
separation has already occured and no effect from the SGWB can be expected
unless coherence is restored, e.g. by matter effects in the so called catch-up
mechanism [38]. Even if coherence would be restored by any effect, it is
questionable, if the neutrino oscillations could be resolved, since that would
require an energy resolution or binning of

$\displaystyle\Delta E\sim 4.17\times
10^{-12}\left(\frac{10^{-5}\,\mathrm{eV}^{2}}{\Delta
m_{jk}^{2}}\right)\left(\frac{E}{\mathrm{MeV}}\right)^{2}\left(\frac{\mathrm{lyr}}{L}\right)\,\mathrm{MeV}\,,$
(30)

cf. equation (17). In addition and as already discussed in section 4.2,
uncertainties in the baseline of astrophysical neutrinos can also affect the
constraints on GW parameters.

##  5 Conclusions

In this paper we have studied the influence of a GW signal on neutrino
oscillations for a variety of baselines and parameter sets. We have developed
a heuristic yet robust model describing the propagation of neutrinos in a
spacetime curved by GWs. After pointing out how GWs influence neutrino flavor
oscillations by altering the propagation length of neutrinos, we have
systematically established conditions under which this effect is observable.
Subsequently, we have performed an NLLR test comparing this effect to the
absence of GW induced decoherence in a toy neutrino experiment that allowed us
to determine the conditions under which an effect becomes observable. For
neutrinos with propagation lengths smaller than the diameter of the solar
system, we obtained a sensitivity of $h\approx${10}^{-7}$$ with realistic wave
packet widths of the order $1\text{\,}\mathrm{nm}$. For larger wave packet
widths we have shown that respective experiments might even reach
sensitivities as low as $h\approx${10}^{-13}$$. While performing experiments
within the solar system poses significant challenges, such as placing a
neutrino source or detector within a Lagrange point or detecting neutrinos
from dark matter decaying in the core of the moon or planets, the situation
differs for astrophysical neutrino sources outside the solar system. In
particular, we have investigated neutrinos originating from the Vela Pulsar
and find a possible sensitivity on GWs ranging from $h\approx${10}^{-7}$$ to
$h\approx${10}^{-2}$$ in the frequency range of $f\approx${10}^{-15}$$ to
$f\approx${10}^{-11}$$, which is currently inaccessible by any other
experiment. Only significantly weaker bounds from nucleosynthesis exist in
this region. For such low frequencies however, the exclusion limits in this
case depend on the phase of the GW which impairs the sensitivity to small
strains.

Furthermore, we have examined the influence of the SGWB on neutrinos
originating from a pulsar, that does not feature the problem of a phase
dependent sensitivity, since the neutrino trajectories are altered by a
randomly distributed ensemble of GW signals of all phases and directions. We
improved on previous investigations of this scenario that have only
qualitatively discussed the effect of wave packet decoherence. When
considering wave packet decoherence and the finite energy resolution and
binning of realistic experiments, we show explicitly that $\mathrm{MeV}$
neutrinos are not sensitive to the SGWB, while $\mathrm{PeV}$ neutrinos still
provide an intriguing possibility to study the SGWB.

In summary, we have discussed how neutrino oscillations are affected by GW
signals. Neutrino oscillations provide a promising new avenue to probe the
SGWB, while the possibility to detect plane GW signals using neutrino flavor
oscillations suffers from relatively low sensitivities especially compared to
laser interferometry experiments, which is mostly due to the interplay of the
GW induced decoherence effect and other sources of oscillation averaging. Even
so, neutrino oscillations may provide a new window in so far unexplored
regions of the GW parameter space.

## References

  * [1] Y. Eroshenko and V. Stasenko, Symmetry 15 (2023) no.3, 637 doi:10.3390/sym15030637 [arXiv:2302.05167 [astro-ph.CO]].
  * [2] J. Casalderrey-Solana, D. Mateos and M. Sanchez-Garitaonandia, [arXiv:2210.03171 [hep-th]].
  * [3] S. Rahvar, Int. J. Mod. Phys. D 32 (2023) no.16, 2350103 doi:10.1142/S0218271823501031 [arXiv:2304.00348 [gr-qc]].
  * [4] J. Powell, A. Iess, M. Llorens-Monteagudo, M. Obergaulinger, B. Müller, A. Torres-Forné, E. Cuoco and J. A. Font, [arXiv:2311.18221 [astro-ph.HE]].
  * [5] A. Eichhorn, R. R. Lino dos Santos and J. L. Miqueleto, Phys. Rev. D 109 (2024) no.2, 2 doi:10.1103/PhysRevD.109.026013 [arXiv:2306.17718 [gr-qc]].
  * [6] E. Bagui et al. [LISA Cosmology Working Group], [arXiv:2310.19857 [astro-ph.CO]].
  * [7] M. C. Guzzetti, N. Bartolo, M. Liguori and S. Matarrese, Riv. Nuovo Cim. 39 (2016) no.9, 399-495 doi:10.1393/ncr/i2016-10127-1 [arXiv:1605.01615 [astro-ph.CO]].
  * [8] G. Agazie et al. [NANOGrav], Astrophys. J. Lett. 951 (2023) no.1, L8 doi:10.3847/2041-8213/acdac6 [arXiv:2306.16213 [astro-ph.HE]].
  * [9] M. Branchesi, M. Maggiore, D. Alonso, C. Badger, B. Banerjee, F. Beirnaert, E. Belgacem, S. Bhagwat, G. Boileau and S. Borhanian, et al. JCAP 07 (2023), 068 doi:10.1088/1475-7516/2023/07/068 [arXiv:2303.15923 [gr-qc]].
  * [10] J. Baker, J. Bellovary, P. L. Bender, E. Berti, R. Caldwell, J. Camp, J. W. Conklin, N. Cornish, C. Cutler and R. DeRosa, et al. [arXiv:1907.06482 [astro-ph.IM]].
  * [11] B. P. Abbott et al. [LIGO Scientific], Rept. Prog. Phys. 72 (2009), 076901 doi:10.1088/0034-4885/72/7/076901 [arXiv:0711.3041 [gr-qc]].
  * [12] M. Dvornikov, Phys. Rev. D 100 (2019) no.9, 096014 doi:10.1103/PhysRevD.100.096014 [arXiv:1906.06167 [hep-ph]].
  * [13] G. Koutsoumbas and D. Metaxas, Gen. Rel. Grav. 52 (2020) no.10, 102 doi:10.1007/s10714-020-02758-z [arXiv:1909.02735 [hep-ph]].
  * [14] M. Dvornikov, J. Phys. Conf. Ser. 1435 (2020) no.1, 012005 doi:10.1088/1742-6596/1435/1/012005 [arXiv:1910.01415 [hep-ph]].
  * [15] M. Dvornikov, JCAP 12 (2020), 022 doi:10.1088/1475-7516/2020/12/022 [arXiv:2009.02195 [astro-ph.HE]].
  * [16] M. Dvornikov, Phys. Rev. D 104 (2021) no.4, 043018 doi:10.1103/PhysRevD.104.043018 [arXiv:2103.15464 [hep-ph]].
  * [17] M. Dvornikov, Moscow Univ. Phys. Bull. 77 (2022) no.2, 439-441 doi:10.3103/S002713492202031X [arXiv:2207.03937 [hep-ph]].
  * [18] G. Lambiase, L. Mastrototaro and L. Visinelli, JCAP 01 (2023), 011 doi:10.1088/1475-7516/2023/01/011 [arXiv:2207.08067 [hep-ph]].
  * [19] G. Lambiase, L. Mastrototaro and L. Visinelli, Phys. Rev. D 108 (2023) no.12, 123028 doi:10.1103/PhysRevD.108.123028 [arXiv:2306.16977 [astro-ph.HE]].
  * [20] A. Domi, T. Eberl, M. J. Fahn, K. Giesel, L. Hennig, U. Katz, R. Kemper and M. Kobler, [arXiv:2403.03106 [gr-qc]].
  * [21] S. Mandal, to be published.
  * [22] L. Stodolsky, Gen. Rel. Grav. 11 (1979), 391-405 doi:10.1007/BF00759302
  * [23] D. Piriz, M. Roy and J. Wudka, Phys. Rev. D 54 (1996), 1587-1599 doi:10.1103/PhysRevD.54.1587 [arXiv:hep-ph/9604403 [hep-ph]].
  * [24] D. V. Ahluwalia and C. Burgard, Gen. Rel. Grav. 28 (1996), 1161-1170 doi:10.1007/BF03218936 [arXiv:gr-qc/9603008 [gr-qc]].
  * [25] D. V. Ahluwalia and C. Burgard, [arXiv:gr-qc/9606031 [gr-qc]].
  * [26] C. Y. Cardall and G. M. Fuller, Phys. Rev. D 55 (1997), 7960-7966 doi:10.1103/PhysRevD.55.7960 [arXiv:hep-ph/9610494 [hep-ph]].
  * [27] N. Fornengo, C. Giunti, C. W. Kim and J. Song, Phys. Rev. D 56 (1997), 1895-1902 doi:10.1103/PhysRevD.56.1895 [arXiv:hep-ph/9611231 [hep-ph]].
  * [28] G. Lambiase, G. Papini, R. Punzi and G. Scarpetta, Phys. Rev. D 71 (2005), 073011 doi:10.1103/PhysRevD.71.073011 [arXiv:gr-qc/0503027 [gr-qc]].
  * [29] L. Visinelli, Gen. Rel. Grav. 47 (2015) no.5, 62 doi:10.1007/s10714-015-1899-z [arXiv:1410.1523 [gr-qc]].
  * [30] B. Kayser, Phys. Rev. D 24 (1981), 110 doi:10.1103/PhysRevD.24.110
  * [31] P. H. Frampton and P. Vogel, Phys. Rept. 82 (1982), 339-388 doi:10.1016/0370-1573(82)90159-4
  * [32] C. Giunti, C. W. Kim and U. W. Lee, Phys. Rev. D 44 (1991), 3635-3640 doi:10.1103/PhysRevD.44.3635
  * [33] C. Giunti, C. W. Kim, J. A. Lee and U. W. Lee, Phys. Rev. D 48 (1993), 4310-4317 doi:10.1103/PhysRevD.48.4310 [arXiv:hep-ph/9305276 [hep-ph]].
  * [34] C. Giunti, C. W. Kim and U. W. Lee, Phys. Lett. B 421 (1998), 237-244 doi:10.1016/S0370-2693(98)00014-8 [arXiv:hep-ph/9709494 [hep-ph]].
  * [35] K. Kiers and N. Weiss, Phys. Rev. D 57 (1998), 3091-3105 doi:10.1103/PhysRevD.57.3091 [arXiv:hep-ph/9710289 [hep-ph]].
  * [36] E. K. Akhmedov and A. Y. Smirnov, Phys. Atom. Nucl. 72 (2009), 1363-1381 doi:10.1134/S1063778809080122 [arXiv:0905.1903 [hep-ph]].
  * [37] E. Akhmedov, J. Kopp and M. Lindner, [arXiv:1405.7275 [hep-ph]].
  * [38] J. Kersten and A. Y. Smirnov, Eur. Phys. J. C 76 (2016) no.6, 339 doi:10.1140/epjc/s10052-016-4187-5 [arXiv:1512.09068 [hep-ph]].
  * [39] E. Akhmedov, J. Kopp and M. Lindner, JCAP 09 (2017), 017 doi:10.1088/1475-7516/2017/09/017 [arXiv:1702.08338 [hep-ph]].
  * [40] D. V. Naumov and V. A. Naumov, Phys. Part. Nucl. 51 (2020) no.1, 1-106 doi:10.1134/S1063779620010050
  * [41] R. L. Workman et al. [Particle Data Group], PTEP 2022 (2022), 083C01 doi:10.1093/ptep/ptac097
  * [42] B. P. Abbott et al. [LIGO Scientific and Virgo], Phys. Rev. Lett. 116 (2016) no.6, 061102 doi:10.1103/PhysRevLett.116.061102 [arXiv:1602.03837 [gr-qc]].
  * [43] C. A. Argüelles, T. Bertólez-Martínez and J. Salvado, Phys. Rev. D 107 (2023) no.3, 036004 doi:10.1103/PhysRevD.107.036004 [arXiv:2201.05108 [hep-ph]].
  * [44] J. Antoniadis et al. [EPTA], Astron. Astrophys. 678 (2023), A48 doi:10.1051/0004-6361/202346841 [arXiv:2306.16224 [astro-ph.HE]].
  * [45] G. Agazie et al. [NANOGrav], Astrophys. J. Lett. 951 (2023) no.1, L9 doi:10.3847/2041-8213/acda9a [arXiv:2306.16217 [astro-ph.HE]].
  * [46] B. P. Abbott, R. Abbott, T. D. Abbott, C. Adams, R. X. Adhikari, R. A. Anderson, S. B. Anderson, K. Arai, M. A. Arain and S. M. Aston, et al. Phys. Rev. D 93 (2016) no.11, 112004 doi:10.1103/PhysRevD.93.112004 [arXiv:1604.00439 [astro-ph.IM]].
  * [47] A. Sesana, Class. Quant. Grav. 30 (2013), 244009 doi:10.1088/0264-9381/30/24/244009 [arXiv:1307.4086 [astro-ph.CO]].
  * [48] C. Q. Ye, H. M. Fan, A. Torres-Orjuela, J. d. Zhang and Y. M. Hu, [arXiv:2310.03520 [gr-qc]].
  * [49] S. Chen, A. Sesana and C. J. Conselice, Mon. Not. Roy. Astron. Soc. 488 (2019) no.1, 401-418 doi:10.1093/mnras/stz1722 [arXiv:1810.04184 [astro-ph.GA]].
  * [50] F. P. An et al. [Daya Bay], Eur. Phys. J. C 77 (2017) no.9, 606 doi:10.1140/epjc/s10052-017-4970-y [arXiv:1608.01661 [hep-ex]].
  * [51] C. J. Moore, R. H. Cole and C. P. L. Berry, Class. Quant. Grav. 32 (2015) no.1, 015014 doi:10.1088/0264-9381/32/1/015014 [arXiv:1408.0740 [gr-qc]].
  * [52] J. Abadie et al. [LIGO Scientific and VIRGO], Phys. Rev. D 81 (2010), 102001 doi:10.1103/PhysRevD.81.102001 [arXiv:1002.1036 [gr-qc]].
  * [53] B. S. Sathyaprakash and B. F. Schutz, Living Rev. Rel. 12 (2009), 2 doi:10.12942/lrr-2009-2 [arXiv:0903.0338 [gr-qc]].
  * [54] M. Kramer et al. [EPTA], Class. Quant. Grav. 30 (2013), 224009 doi:10.1088/0264-9381/30/22/224009
  * [55] J. Antoniadis, Z. Arzoumanian, S. Babak, M. Bailes, A. S. B. Nielsen, P. T. Baker, C. G. Bassa, B. Becsy, A. Berthereau and M. Bonetti, et al. Mon. Not. Roy. Astron. Soc. 510 (2022) no.4, 4873-4887 doi:10.1093/mnras/stab3418 [arXiv:2201.03980 [astro-ph.HE]].
  * [56] P. Amaro-Seoane, S. Aoudia, S. Babak, P. Binetruy, E. Berti, A. Bohe, C. Caprini, M. Colpi, N. J. Cornish and K. Danzmann, et al. GW Notes 6 (2013), 4-110 [arXiv:1201.3621 [astro-ph.CO]].
  * [57] J. F. Beacom, N. F. Bell, D. Hooper, S. Pakvasa and T. J. Weiler, Phys. Rev. D 68 (2003), 093005 [erratum: Phys. Rev. D 72 (2005), 019901] doi:10.1103/PhysRevD.68.093005 [arXiv:hep-ph/0307025 [hep-ph]].
  * [58] H. Athar, C. S. Kim and J. Lee, Mod. Phys. Lett. A 21 (2006), 1049-1066 doi:10.1142/S021773230602038X [arXiv:hep-ph/0505017 [hep-ph]].
  * [59] J. P. Rachen and P. Meszaros, Phys. Rev. D 58 (1998), 123005 doi:10.1103/PhysRevD.58.123005 [arXiv:astro-ph/9802280 [astro-ph]].
  * [60] L. A. Anchordoqui, H. Goldberg, F. Halzen and T. J. Weiler, Phys. Lett. B 593 (2004), 42 doi:10.1016/j.physletb.2004.04.054 [arXiv:astro-ph/0311002 [astro-ph]].
  * [61] T. Kashti and E. Waxman, Phys. Rev. Lett. 95 (2005), 181101 doi:10.1103/PhysRevLett.95.181101 [arXiv:astro-ph/0507599 [astro-ph]].
  * [62] M. Kachelrieß and R. Tomàs, Phys. Rev. D 74 (2006), 063009 doi:10.1103/PhysRevD.74.063009
  * [63] P. Lipari, M. Lusignoli and D. Meloni, Phys. Rev. D 75 (2007), 123005 doi:10.1103/PhysRevD.75.123005 [arXiv:0704.0718 [astro-ph]].
  * [64] S. R. Klein, R. E. Mikkelsen and J. Becker Tjus, Astrophys. J. 779 (2013), 106 doi:10.1088/0004-637X/779/2/106 [arXiv:1208.2056 [astro-ph.HE]].
  * [65] M. Bustamante, J. F. Beacom and W. Winter, Phys. Rev. Lett. 115 (2015) no.16, 161302 doi:10.1103/PhysRevLett.115.161302 [arXiv:1506.02645 [astro-ph.HE]].
  * [66] H. T. Janka, doi:10.1007/978-3-319-21846-5_4 [arXiv:1702.08713 [astro-ph.HE]].
  * [67] A. M. Suliga, doi:10.1007/978-981-15-8818-1_129-1 [arXiv:2207.09632 [astro-ph.HE]].
  * [68] F. Oikonomou, PoS ICRC2021 (2022), 030 doi:10.22323/1.395.0030 [arXiv:2201.05623 [astro-ph.HE]].
  * [69] M. Boettcher, M. Fu, T. Govenor, Q. King and P. Roustazadeh, Acta Phys. Polon. Supp. 15 (2022) no.3, 8 doi:10.5506/APhysPolBSupp.15.3-A8 [arXiv:2204.12242 [astro-ph.HE]].
  * [70] F. Oikonomou, K. Murase, P. Padovani, E. Resconi and P. Mészáros, Mon. Not. Roy. Astron. Soc. 489 (2019) no.3, 4347-4366 doi:10.1093/mnras/stz2246 [arXiv:1906.05302 [astro-ph.HE]].
  * [71] R. J. Protheroe, W. Bednarek and Q. Luo, Astropart. Phys. 9 (1998), 1-14 doi:10.1016/S0927-6505(98)00014-0 [arXiv:astro-ph/9703045 [astro-ph]].
  * [72] R. Gandhi, Nucl. Phys. B Proc. Suppl. 91 (2001), 453-461 doi:10.1016/S0920-5632(00)00975-0 [arXiv:hep-ph/0011176 [hep-ph]].
  * [73] J. H. Beall and W. Bednarek, Astrophys. J. 569 (2002), 343-348 doi:10.1086/339276 [arXiv:astro-ph/0108447 [astro-ph]].
  * [74] D. Guetta and E. Amato, Astropart. Phys. 19 (2003), 403-407 doi:10.1016/S0927-6505(02)00221-9 [arXiv:astro-ph/0209537 [astro-ph]].
  * [75] D. Guetta and J. Granot, Phys. Rev. Lett. 90 (2003), 201103 doi:10.1103/PhysRevLett.90.201103 [arXiv:astro-ph/0212045 [astro-ph]].
  * [76] S. Razzaque, P. Meszaros and E. Waxman, Phys. Rev. Lett. 90 (2003), 241103 doi:10.1103/PhysRevLett.90.241103 [arXiv:astro-ph/0212536 [astro-ph]].
  * [77] C. D. Dermer and A. Atoyan, Phys. Rev. Lett. 91 (2003), 071102 doi:10.1103/PhysRevLett.91.071102 [arXiv:astro-ph/0301030 [astro-ph]].
  * [78] W. Bednarek, Astron. Astrophys. 407 (2003), 1-6 doi:10.1051/0004-6361:20030929 [arXiv:astro-ph/0305430 [astro-ph]].
  * [79] B. Link and F. Burgio, Phys. Rev. Lett. 94 (2005), 181101 doi:10.1103/PhysRevLett.94.181101 [arXiv:astro-ph/0412520 [astro-ph]].
  * [80] B. Link and F. Burgio, Mon. Not. Roy. Astron. Soc. 371 (2006), 375-379 doi:10.1111/j.1365-2966.2006.10665.x [arXiv:astro-ph/0604379 [astro-ph]].
  * [81] K. Fang, JCAP 06 (2015), 004 doi:10.1088/1475-7516/2015/06/004 [arXiv:1411.2174 [astro-ph.HE]].
  * [82] P. A. Caraveo, A. De Luca, R. P. Mignani and G. F. Bignami, Astrophys. J. 561 (2001), 930 doi:10.1086/323377 [arXiv:astro-ph/0107282 [astro-ph]].
  * [83] R. Lin, M. H. van Kerkwijk, F. Kirsten, U. L. Pen and A. T. Deller, Astrophys. J. 952 (2023) no.2, 161 doi:10.3847/1538-4357/acdc98 [arXiv:2306.01617 [astro-ph.HE]].
  * [84] G. G. Pavlov, D. Sanwal and M. A. Teter, IAU Symp. 218 (2004), 239 [arXiv:astro-ph/0311526 [astro-ph]].
  * [85] C. Guépin, B. Cerutti and K. Kotera, Astron. Astrophys. 635 (2020), A138 doi:10.1051/0004-6361/201936816 [arXiv:1910.11387 [astro-ph.HE]].
  * [86] A. Asthana, A. A. Mushtukov, A. A. Dobrynina and I. S. Ognev, Mon. Not. Roy. Astron. Soc. 522 (2023) no.3, 3405-3411 doi:10.1093/mnras/stad1118 [arXiv:2304.04520 [astro-ph.HE]].
  * [87] R. Abbasi, M. Ackermann, J. Adams, S. K. Agarwalla, J. A. Aguilar, M. Ahlers, J. M. Alameddine, N. M. Amin, K. Andeen and G. Anton, et al. [arXiv:2402.18026 [astro-ph.HE]].
  * [88] R. M. Shannon, V. Ravi, W. A. Coles, G. Hobbs, M. J. Keith, R. N. Manchester, J. S. B. Wyithe, M. Bailes, N. D. R. Bhat and S. Burke-Spolaor, et al. Science 342 (2013) no.6156, 334-337 doi:10.1126/science.1238012 [arXiv:1310.4569 [astro-ph.CO]].
  * [89] L. M. Krauss and M. J. White, Phys. Rev. Lett. 69 (1992), 869-872 doi:10.1103/PhysRevLett.69.869 [arXiv:hep-ph/9205212 [hep-ph]].
  * [90] K. S. Thorne, [arXiv:gr-qc/9506086 [gr-qc]].
  * [91] P. D. Bolton, F. F. Deppisch and P. S. Bhupal Dev, JHEP 03 (2020), 170 doi:10.1007/JHEP03(2020)170 [arXiv:1912.03058 [hep-ph]].
  * [92] J. M. Berryman, Phys. Rev. D 100 (2019) no.2, 023540 doi:10.1103/PhysRevD.100.023540 [arXiv:1905.03254 [hep-ph]].
  * [93] I. Alekseev et al. [DANSS], Phys. Lett. B 787 (2018), 56-63 doi:10.1016/j.physletb.2018.10.038 [arXiv:1804.04046 [hep-ex]].
  * [94] M. Cirelli, G. Marandella, A. Strumia and F. Vissani, Nucl. Phys. B 708 (2005), 215-267 doi:10.1016/j.nuclphysb.2004.11.056 [arXiv:hep-ph/0403158 [hep-ph]].
  * [95] M. Dentler, Á. Hernández-Cabezudo, J. Kopp, P. A. N. Machado, M. Maltoni, I. Martinez-Soler and T. Schwetz, JHEP 08 (2018), 010 doi:10.1007/JHEP08(2018)010 [arXiv:1803.10661 [hep-ph]].
  * [96] Y. J. Ko et al. [NEOS], Phys. Rev. Lett. 118 (2017) no.12, 121802 doi:10.1103/PhysRevLett.118.121802 [arXiv:1610.05134 [hep-ex]].
  * [97] K. H. Hiddemann, H. Daniel and O. Schwentker, J. Phys. G 21 (1995), 639-650 doi:10.1088/0954-3899/21/5/008
  * [98] C. Kraus, A. Singer, K. Valerius and C. Weinheimer, Eur. Phys. J. C 73 (2013) no.2, 2323 doi:10.1140/epjc/s10052-013-2323-z [arXiv:1210.4194 [hep-ex]].
  * [99] A. I. Belesev, A. I. Berlev, E. V. Geraskin, A. A. Golubev, N. A. Likhovid, A. A. Nozik, V. S. Pantuev, V. I. Parfenov and A. K. Skasyrskaya, J. Phys. G 41 (2014), 015001 doi:10.1088/0954-3899/41/1/015001 [arXiv:1307.5687 [hep-ex]].
  * [100] J. N. Abdurashitov, A. I. Belesev, V. G. Chernov, E. V. Geraskin, A. A. Golubev, P. V. Grigorieva, G. A. Koroteev, N. A. Likhovid, A. A. Nozik and V. S. Pantuev, et al. Pisma Zh. Eksp. Teor. Fiz. 105 (2017) no.12, 723-724 doi:10.1134/S0021364017120013 [arXiv:1703.10779 [hep-ex]].
  * [101] E. Holzschuh, W. Kundig, L. Palermo, H. Stussi and P. Wenk, Phys. Lett. B 451 (1999), 247-255 doi:10.1016/S0370-2693(99)00200-2
  * [102] E. Holzschuh, L. Palermo, H. Stussi and P. Wenk, Phys. Lett. B 482 (2000), 1-9 doi:10.1016/S0370-2693(00)00476-7
  * [103] A. V. Derbin, A. I. Egorov, S. V. Bakhlanov and V. N. Muratova, JETP Lett. 66 (1997), 88-92 doi:10.1134/1.567508
  * [104] K. Schreckenbach, G. Colvin and F. Von Feilitzsch, Phys. Lett. B 129 (1983), 265-268 doi:10.1016/0370-2693(83)90858-4
  * [105] S. Friedrich, G. B. Kim, C. Bray, R. Cantor, J. Dilling, S. Fretwell, J. A. Hall, A. Lennarz, V. Lordi and P. Machule, et al. Phys. Rev. Lett. 126 (2021) no.2, 021803 doi:10.1103/PhysRevLett.126.021803 [arXiv:2010.09603 [nucl-ex]].
  * [106] A. Aguilar-Arevalo et al. [PIENU], Phys. Rev. D 97 (2018) no.7, 072012 doi:10.1103/PhysRevD.97.072012 [arXiv:1712.03275 [hep-ex]].
  * [107] D. A. Bryman and R. Shrock, Phys. Rev. D 100 (2019) no.5, 053006 doi:10.1103/PhysRevD.100.053006 [arXiv:1904.06787 [hep-ph]].
  * [108] D. A. Bryman and R. Shrock, Phys. Rev. D 100 (2019), 073011 doi:10.1103/PhysRevD.100.073011 [arXiv:1909.11198 [hep-ph]].
  * [109] G. Bernardi, G. Carugno, J. Chauveau, F. Dicarlo, M. Dris, J. Dumarchez, M. Ferro-Luzzi, J. M. Levy, D. Lukas and J. M. Perreau, et al. Phys. Lett. B 203 (1988), 332-334 doi:10.1016/0370-2693(88)90563-1
  * [110] S. A. Baranov, Y. A. Batusov, A. A. Borisov, S. A. Bunyatov, V. Y. Valuev, A. S. Vovenko, V. N. Goryachev, M. M. Kirsanov, D. Kish and O. L. Klimov, et al. Phys. Lett. B 302 (1993), 336-340 doi:10.1016/0370-2693(93)90405-7
  * [111] E. Cortina Gil et al. [NA62], Phys. Lett. B 807 (2020), 135599 doi:10.1016/j.physletb.2020.135599 [arXiv:2005.09575 [hep-ex]].
  * [112] K. Abe et al. [T2K], Phys. Rev. D 100 (2019) no.5, 052006 doi:10.1103/PhysRevD.100.052006 [arXiv:1902.07598 [hep-ex]].
  * [113] K. Abe et al. [Super-Kamiokande], Phys. Rev. D 91 (2015), 052019 doi:10.1103/PhysRevD.91.052019 [arXiv:1410.2008 [hep-ex]].
  * [114] S. Kovalenko, Z. Lu and I. Schmidt, Phys. Rev. D 80 (2009), 073014 doi:10.1103/PhysRevD.80.073014 [arXiv:0907.2533 [hep-ph]].
  * [115] A. Atre, T. Han, S. Pascoli and B. Zhang, JHEP 05 (2009), 030 doi:10.1088/1126-6708/2009/05/030 [arXiv:0901.3589 [hep-ph]].
  * [116] J. C. Helo, S. Kovalenko and I. Schmidt, Nucl. Phys. B 853 (2011), 80-104 doi:10.1016/j.nuclphysb.2011.07.020 [arXiv:1005.1607 [hep-ph]].
  * [117] A. Abada, V. De Romeri, M. Lucente, A. M. Teixeira and T. Toma, JHEP 02 (2018), 169 doi:10.1007/JHEP02(2018)169 [arXiv:1712.03984 [hep-ph]].
  * [118] F. Bergsma et al. [CHARM], Phys. Lett. B 166 (1986), 473-478 doi:10.1016/0370-2693(86)91601-1
  * [119] P. Vilain et al. [CHARM II], Phys. Lett. B 343 (1995), 453-458 doi:10.1016/0370-2693(94)01422-9
  * [120] R. Barouki, G. Marocco and S. Sarkar, SciPost Phys. 13 (2022), 118 doi:10.21468/SciPostPhys.13.5.118 [arXiv:2208.00416 [hep-ph]].
  * [121] M. Ablikim et al. [BESIII], Phys. Rev. D 99 (2019) no.11, 112002 doi:10.1103/PhysRevD.99.112002 [arXiv:1902.02450 [hep-ex]].
  * [122] J. Badier et al. [NA3], Z. Phys. C 31 (1986), 21 doi:10.1007/BF01559588
  * [123] D. Liventsev et al. [Belle], Phys. Rev. D 87 (2013) no.7, 071102 [erratum: Phys. Rev. D 95 (2017) no.9, 099903] doi:10.1103/PhysRevD.87.071102 [arXiv:1301.1105 [hep-ex]].
  * [124] P. Abreu et al. [DELPHI], Z. Phys. C 74 (1997), 57-71 [erratum: Z. Phys. C 75 (1997), 580] doi:10.1007/s002880050370
  * [125] O. Adriani et al. [L3], Phys. Lett. B 295 (1992), 371-382 doi:10.1016/0370-2693(92)91579-X
  * [126] P. Achard et al. [L3], Phys. Lett. B 517 (2001), 67-74 doi:10.1016/S0370-2693(01)00993-5 [arXiv:hep-ex/0107014 [hep-ex]].
  * [127] G. Aad et al. [ATLAS], JHEP 07 (2015), 162 doi:10.1007/JHEP07(2015)162 [arXiv:1506.06020 [hep-ex]].
  * [128] G. Aad et al. [ATLAS], JHEP 10 (2019), 265 doi:10.1007/JHEP10(2019)265 [arXiv:1905.09787 [hep-ex]].
  * [129] A. M. Sirunyan et al. [CMS], JHEP 01 (2019), 122 doi:10.1007/JHEP01(2019)122 [arXiv:1806.10905 [hep-ex]].
  * [130] A. M. Sirunyan et al. [CMS], Phys. Rev. Lett. 120 (2018) no.22, 221801 doi:10.1103/PhysRevLett.120.221801 [arXiv:1802.02965 [hep-ex]].
  * [131] A. Das, P. S. B. Dev and C. S. Kim, Phys. Rev. D 95 (2017) no.11, 115013 doi:10.1103/PhysRevD.95.115013 [arXiv:1704.00880 [hep-ph]].
  * [132] M. Blennow, P. Coloma, E. Fernandez-Martinez, J. Hernandez-Garcia and J. Lopez-Pavon, JHEP 04 (2017), 153 doi:10.1007/JHEP04(2017)153 [arXiv:1609.08637 [hep-ph]].
  * [133] F. del Aguila, J. de Blas and M. Perez-Victoria, Phys. Rev. D 78 (2008), 013010 doi:10.1103/PhysRevD.78.013010 [arXiv:0803.4008 [hep-ph]].
  * [134] J. de Blas, EPJ Web Conf. 60 (2013), 19008 doi:10.1051/epjconf/20136019008 [arXiv:1307.6173 [hep-ph]].
  * [135] S. Antusch and O. Fischer, JHEP 10 (2014), 094 doi:10.1007/JHEP10(2014)094 [arXiv:1407.6607 [hep-ph]].
  * [136] S. Bridle, J. Elvin-Poole, J. Evans, S. Fernandez, P. Guzowski and S. Soldner-Rembold, Phys. Lett. B 764 (2017), 322-327 doi:10.1016/j.physletb.2016.11.050 [arXiv:1607.00032 [astro-ph.CO]].
  * [137] A. C. Vincent, E. F. Martinez, P. Hernández, M. Lattanzi and O. Mena, JCAP 04 (2015), 006 doi:10.1088/1475-7516/2015/04/006 [arXiv:1408.1956 [astro-ph.CO]].
  * [138] A. Boyarsky, O. Ruchayskiy and M. Shaposhnikov, Ann. Rev. Nucl. Part. Sci. 59 (2009), 191-214 doi:10.1146/annurev.nucl.010909.083654 [arXiv:0901.0011 [hep-ph]].
  * [139] O. Ruchayskiy and A. Ivashko, JCAP 10 (2012), 014 doi:10.1088/1475-7516/2012/10/014 [arXiv:1202.2841 [hep-ph]].
  * [140] K. C. Y. Ng, B. M. Roach, K. Perez, J. F. Beacom, S. Horiuchi, R. Krivonos and D. R. Wik, Phys. Rev. D 99 (2019), 083005 doi:10.1103/PhysRevD.99.083005 [arXiv:1901.01262 [astro-ph.HE]].
  * [141] B. M. Roach, K. C. Y. Ng, K. Perez, J. F. Beacom, S. Horiuchi, R. Krivonos and D. R. Wik, Phys. Rev. D 101 (2020) no.10, 103011 doi:10.1103/PhysRevD.101.103011 [arXiv:1908.09037 [astro-ph.HE]].
  * [142] X. Shi and G. Sigl, Phys. Lett. B 323 (1994), 360-366 [erratum: Phys. Lett. B 324 (1994), 516-516] doi:10.1016/0370-2693(94)91232-7 [arXiv:hep-ph/9312247 [hep-ph]].
  * [143] G. G. Raffelt and S. Zhou, Phys. Rev. D 83 (2011), 093014 doi:10.1103/PhysRevD.83.093014 [arXiv:1102.5124 [hep-ph]].
  * [144] S. Zhou, Int. J. Mod. Phys. A 30 (2015) no.13, 1530033 doi:10.1142/S0217751X15300331 [arXiv:1504.02729 [hep-ph]].
  * [145] H. Nunokawa, J. T. Peltoniemi, A. Rossi and J. W. F. Valle, Phys. Rev. D 56 (1997), 1704-1713 doi:10.1103/PhysRevD.56.1704 [arXiv:hep-ph/9702372 [hep-ph]].
  * [146] M. Galeazzi, F. Fontanelli, F. Gatti and S. Vitale, Phys. Rev. Lett. 86 (2001), 1978-1981 doi:10.1103/PhysRevLett.86.1978
  * [147] F. P. Calaprice and D. J. Millener, Phys. Rev. C 27 (1983), 1175-1181 doi:10.1103/PhysRevC.27.1175
  * [148] A. I. Derbin, A. V. Chernyi, L. A. Popeko, V. N. Muratova, G. A. Shishkina and S. I. Bakhlanov, JETP Lett. 57 (1993), 768-772
  * [149] C. Hagner, M. Altmann, F. von Feilitzsch, L. Oberauer, Y. Declais and E. Kajfasz, Phys. Rev. D 52 (1995), 1343-1352 doi:10.1103/PhysRevD.52.1343
  * [150] G. Bellini et al. [Borexino], Phys. Rev. D 88 (2013) no.7, 072010 doi:10.1103/PhysRevD.88.072010 [arXiv:1311.5347 [hep-ex]].

##  Appendix A Geodesics in Gravitational Wave Spacetimes: The Boundary Value
Problem

Within our heuristic approach, we replace the baseline with a modified length,
as it is subject to the effects of GW as can be seen in equation (4). In order
to calculate this spatial length, we first rewrite equation (4) by using the
on-shell condition

$\displaystyle 1$
$\displaystyle=g_{\mu\nu}\dot{\gamma}^{\mu}\dot{\gamma}^{\nu}$ (31)
$\displaystyle\Leftrightarrow 1$
$\displaystyle=(\dot{\gamma}^{0})^{2}+g_{jk}\dot{\gamma}^{j}\dot{\gamma}^{k}$
(32) $\displaystyle\Leftrightarrow-g_{jk}\dot{\gamma}^{j}\dot{\gamma}^{k}$
$\displaystyle=(\dot{\gamma}^{0})^{2}-1\,.$ (33)

Here we furthermore employed the transverse traceless (TT) gauge implying
$g_{j0}\equiv 0$ and $g_{00}\equiv 1$ in our convention. Substituting this
result into equation (4) yields

$L=\int_{\tau_{P}}^{\tau_{D}}\left(\left(\dot{\gamma}^{0}\right)^{2}-1\right)^{\frac{1}{2}}\text{d}\tau\,.$
(34)

To calculate this modification, it is initially necessary to compute the
neutrino geodesic in our GW spacetime. To this end, we solve the geodesic
equation for given boundary conditions (BC), i.e.

$\displaystyle\ddot{\gamma}^{\mu}(\tau)$
$\displaystyle=-\Gamma^{\mu}_{\alpha\beta}(\gamma(\tau))\dot{\gamma}^{\alpha}(\tau)\dot{\gamma}^{\beta}(\tau)\,,$
(35) $\displaystyle\gamma^{\alpha}(\tau_{P})$
$\displaystyle=X_{P}^{\alpha}\,,$ (36)
$\displaystyle\gamma^{\alpha}(\tau_{D})$ $\displaystyle=X_{D}^{\alpha}\,.$
(37)

Since our metric is close to $\eta=\mathrm{diag(1,-1,-1,-1)}$, we adopt the
ansatz

$\displaystyle\gamma^{\alpha}(\tau)=\gamma_{0}^{\alpha}(\tau)+\epsilon\Delta\gamma^{\alpha}(\tau)\,.$
(38)

In general the boundary value problem (BVP) of the geodesic equation is not
well posed for arbitrarily far separated production and detection points since
there may be more than one shortest curve connecting these points. In our
specific scenario, however, we assume that geodesics between two points are
inherently unique because our spacetime closely resembles Minkowski spacetime.

Upon substituting this ansatz into the geodesic equation, we obtain, keeping
only terms up to first order in $\epsilon$,

$\ddot{\gamma}_{0}^{\mu}(\tau)+\epsilon\Delta\ddot{\gamma}^{\mu}(\tau)\approx-\frac{\epsilon}{2}\eta^{\mu\nu}\left(\partial_{\alpha}h_{\nu\beta}+\partial_{\beta}h_{\nu\alpha}-\partial_{\nu}h_{\alpha\beta}\right)|_{\gamma_{0}(\tau)}\dot{\gamma}_{0}^{\alpha}\dot{\gamma}_{0}^{\beta}\,.$
(39)

By matching terms of the same order in $\epsilon$ we arrive at the following
set of equations:

$\displaystyle\ddot{\gamma}_{0}^{\mu}(\tau)$ $\displaystyle=0\,,$ (40)
$\displaystyle\Delta\ddot{\gamma}^{\mu}(\tau)$
$\displaystyle=-\frac{\eta^{\mu\nu}}{2}\left(\partial_{\alpha}h_{\nu\beta}+\partial_{\beta}h_{\nu\alpha}-\partial_{\nu}h_{\alpha\beta}\right)|_{\gamma_{0}(\tau)}\dot{\gamma}_{0}^{\alpha}\dot{\gamma}_{0}^{\beta}\,.$
(41)

Equation (40) is easily solved by $\gamma_{0}^{\mu}(\lambda)=a^{\mu}+\lambda
u^{\mu}$. Using our BC from above (at zeroth order in $\epsilon$) on the flat
spacetime part of the geodesic, the constants of integration can be determined
to yield

$\displaystyle\gamma_{0}^{\mu}(\tau)=X_{P}^{\mu}+\left(\tau-\tau_{P}\right)u^{\mu}\,,$
(42)

where $u^{\mu}=(X_{D}^{\mu}-X_{P}^{\mu})/(\tau_{D}-\tau_{P})$. For the
geodesic correction $\Delta\gamma$, we hence obtain the differential equation

$\Delta\ddot{\gamma}^{\mu}(\tau)=-\frac{\eta^{\mu\nu}}{2}\left(2\partial_{\alpha}h_{\nu\beta}-\partial_{\nu}h_{\alpha\beta}\right)|_{\gamma_{0}(\tau)}u^{\alpha}u^{\beta}=:\Delta
a^{\mu}(\tau)\,,$ (43)

and BC, $\Delta\gamma(\tau_{P})=\Delta\gamma(\tau_{D})=0$, at first order in
$\epsilon$. Integrating $\Delta\ddot{\gamma}^{\mu}(\tau)$ leads to the
correction to the four velocity, i.e.

$\Delta\dot{\gamma}^{\mu}(\tau)=\Delta
u_{P}^{\mu}+\int_{\tau_{P}}^{\tau}\Delta
a^{\mu}(\tau^{\prime})d\tau^{\prime}\,.$ (44)

Integrating once again finally results in

$\displaystyle\Delta\gamma^{\mu}(\tau)=\Delta\gamma_{P}^{\mu}+\Delta
u_{P}^{\mu}(\tau_{D}-\tau_{P})+\int_{\tau_{P}}^{\tau}\int_{\tau_{P}}^{\tau^{\prime}}\Delta
a^{\mu}(\tau^{\prime\prime})\;\mathrm{d}\tau^{\prime\prime}\mathrm{d}\tau^{\prime}\,.$
(45)

By using our BC from equation (36) and (37), we are able to determine both
constants of integration, i.e. $\Delta\gamma_{P}^{\mu}$ as well as $\Delta
u_{P}^{\mu}$:

$\displaystyle 0$
$\displaystyle\stackrel{{\scriptstyle!}}{{=}}\Delta\gamma^{\mu}(\tau_{P})=\Delta\gamma_{P}^{\mu}\,,$
(46) $\displaystyle 0$
$\displaystyle\stackrel{{\scriptstyle!}}{{=}}\Delta\gamma^{\mu}(\tau_{D})=\Delta
u_{P}^{\mu}\left(\tau_{D}-\tau_{P}\right)\
+\int_{\tau_{P}}^{\tau_{D}}\text{d}\tau^{\prime}\int_{\tau_{P}}^{\tau^{\prime}}\text{d}\tau^{\prime\prime}\Delta
a^{\mu}(\tau^{\prime\prime})\,,$ (47) $\displaystyle\Leftrightarrow\Delta
u_{P}^{\mu}$
$\displaystyle=-\Delta\bar{u}^{\mu}=\frac{-1}{\tau_{D}-\tau_{P}}\int_{\tau_{P}}^{\tau_{D}}\text{d}\tau^{\prime}\int_{\tau_{P}}^{\tau^{\prime}}\text{d}\tau^{\prime\prime}\Delta
a^{\mu}(\tau^{\prime\prime})\,.$ (48)

Here we also defined $\Delta\bar{u}^{\mu}$, i.e. the average velocity
correction between $X_{P}$ and $X_{D}$. In total, the solution is given by

$\displaystyle\gamma^{\mu}(\tau)$
$\displaystyle=u^{\mu}(\tau-\tau_{P})+X_{P}^{\mu}-\frac{\tau-\tau_{P}}{\tau_{D}-\tau_{P}}\int_{\tau_{P}}^{\tau_{D}}\int_{\tau_{P}}^{\tau^{\prime}}\Delta
a^{\mu}(\tau^{\prime\prime})\;\mathrm{d}\tau^{\prime\prime}\mathrm{d}\tau^{\prime}+\int_{\tau_{P}}^{\tau}\int_{\tau_{P}}^{\tau^{\prime}}\Delta
a^{\mu}(\tau^{\prime\prime})\;\mathrm{d}\tau^{\prime\prime}\mathrm{d}\tau^{\prime}\,,$
(49) $\displaystyle\Delta a^{\mu}(\tau)$
$\displaystyle=-\frac{\eta^{\mu\nu}}{2}\left(2\partial_{\alpha}h_{\nu\beta}-\partial_{\nu}h_{\alpha\beta}\right)|_{\gamma_{0}(\tau)}u^{\alpha}u^{\beta}\,.$
(50)

Since this curve solves the geodesic equation to first order in epsilon and we
have parametrized it with respect to the particles proper time, the on-shell
condition

$\displaystyle 1\equiv
g_{\mu\nu}(\gamma(\tau))\dot{\gamma}^{\mu}(\tau)\dot{\gamma}^{\nu}(\tau)\,,$
(51)

enforced by the BC, also holds along the whole trajectory. In our specific
scenario the on-shell relation can be approximated as

$\displaystyle 1\approx(\eta_{\mu\nu}+\epsilon
h_{\mu\nu}(\gamma_{0}(\tau)))u^{\mu}u^{\nu}+\underbrace{2\eta_{\mu\nu}\Delta\epsilon\dot{\gamma}^{\mu}(\tau)u^{\nu}}_{(\ast)}\,.$
(52)

The last term, $(\ast)$, expands to

$\displaystyle 2\eta_{\mu\nu}\Delta\epsilon\dot{\gamma}^{\mu}(\tau)u^{\nu}$
$\displaystyle=-\frac{2\epsilon}{\tau_{D}-\tau_{P}}\int_{\tau_{P}}^{\tau_{D}}\int_{\tau_{P}}^{\tau}\eta_{\mu\nu}\Delta
a^{\mu}(\tau^{\prime})u^{\nu}\;\mathrm{d}\tau^{\prime}\mathrm{d}\tau+2\epsilon\int_{\tau_{P}}^{\tau}\eta_{\mu\nu}\Delta
a^{\mu}(\tau^{\prime})u^{\nu}\;\mathrm{d}\tau^{\prime}\,,$ (53)

where

$\displaystyle 2\eta_{\mu\nu}\Delta a^{\mu}(\tau^{\prime})u^{\nu}$
$\displaystyle=-2\eta_{\mu\nu}\frac{\eta^{\mu\sigma}}{2}\left(2\partial_{\alpha}h_{\sigma\beta}-\partial_{\sigma}h_{\alpha\beta}\right)|_{\gamma_{0}(\tau)}u^{\alpha}u^{\beta}u^{\nu}$
(54)
$\displaystyle=-\left(2\partial_{\alpha}h_{\nu\beta}-\partial_{\nu}h_{\alpha\beta}\right)|_{\gamma_{0}(\tau)}u^{\alpha}u^{\beta}u^{\nu}$
(55)
$\displaystyle=-\partial_{\alpha}h_{\nu\beta}(\gamma_{0}(\tau))u^{\alpha}u^{\beta}u^{\nu}$
(56)
$\displaystyle=-\frac{\mathrm{d}}{\mathrm{d}\tau}\left(h_{\nu\beta}(\gamma_{0}(\tau))u^{\beta}u^{\nu}\right)\,.$
(57)

Plugging this result back into $(\ast)$ and using the fundamental theorem of
calculus yields

$\displaystyle\begin{split}2\eta_{\mu\nu}\epsilon\Delta\dot{\gamma}^{\mu}(\tau)u^{\nu}&=\frac{\epsilon}{\tau_{D}-\tau_{P}}\int_{\tau_{P}}^{\tau_{D}}\left[h_{\alpha\beta}(\gamma_{0}(\tau))u^{\alpha}u^{\beta}-h_{\alpha\beta}(X_{P})u^{\alpha}u^{\beta}\right]\;\mathrm{d}\tau\\\
&\quad-\epsilon\left[h_{\alpha\beta}(\gamma_{0}(\tau))u^{\alpha}u^{\beta}-h_{\alpha\beta}(X_{P})u^{\alpha}u^{\beta}\right]\end{split}$
(58)
$\displaystyle=\frac{\epsilon}{\tau_{D}-\tau_{P}}\int_{\tau_{P}}^{\tau_{D}}h_{\alpha\beta}(\gamma_{0}(\tau))u^{\alpha}u^{\beta}\;\mathrm{d}\tau-\epsilon
h_{\alpha\beta}(\gamma_{0}(\tau))u^{\alpha}u^{\beta}\,.$ (59)

Finally, substituting $(\ast)$ back into the on-shell condition yields

$\displaystyle 1$ $\displaystyle\approx(\eta_{\mu\nu}+\epsilon
h_{\mu\nu}(\gamma_{0}(\tau)))u^{\mu}u^{\nu}+\frac{\epsilon}{\tau_{D}-\tau_{P}}\int_{\tau_{P}}^{\tau_{D}}h_{\alpha\beta}(\gamma_{0}(\tau))u^{\alpha}u^{\beta}\;\mathrm{d}\tau-\epsilon
h_{\alpha\beta}(\gamma_{0}(\tau))u^{\alpha}u^{\beta}$ (60)
$\displaystyle=\eta_{\mu\nu}u^{\mu}u^{\nu}+\frac{\epsilon}{\tau_{D}-\tau_{P}}\int_{\tau_{P}}^{\tau_{D}}h_{\mu\nu}(\gamma_{0}(\tau))u^{\mu}u^{\nu}\;\mathrm{d}\tau\,.$
(61)

In the next section, we will use this result and the ultra relativistic
approximation in order to obtain the correction to the spatial distance
traveled by a neutrino in a GW spacetime.

##  Appendix B Modified Neutrino Path Length for general Gravitational Waves

Using equation (34) and the $\mathcal{O}(\epsilon)$ solution to the geodesic
equation found in the last section, we are now able to compute the corrected
physical length. First we show that $L(t_{P},t_{D})\approx t_{D}-t_{P}$ in the
ultra relativistic limit and subsequently, we use this fact and the on-shell
condition (61) to derive a relation between the physical, spatial distance
between neutrino production and detection and the coordinate distance,

$\displaystyle L_{0}:=||\vec{X}_{D}-\vec{X}_{P}||\,.$ (62)

The latter distance may be identified with the flat spacetime distance, as the
physical distance reduces to the coordinate distance for $h\equiv 0$. Indeed
we will see that in most cases the physical distance oscillates around this
flat spacetime value justifying that to some degree it still has physical
significance.

First we prove that $L(t_{P},t_{D})\approx t_{D}-t_{P}$ by evaluating equation
(34). Since we have already rewritten $L$ in terms of $\dot{\gamma^{0}}$ using
TT gauge this task is trivial:

$\displaystyle
L(t_{P},t_{D})=\int_{\tau_{P}=\tau(t_{P})}^{\tau_{D}=\tau(t_{D})}\left(\left(\dot{\gamma}^{0}\right)^{2}-1\right)^{\frac{1}{2}}\text{d}\tau\stackrel{{\scriptstyle\dot{\gamma}^{0}\gg
1}}{{\approx}}\int_{\tau_{P}}^{\tau_{D}}\dot{\gamma}^{0}\text{d}\tau=t_{D}-t_{P}\,.$
(63)

The other important ingredient is the deviation of the flow of proper time
from the flat spacetime case quantified in equation (61),

$\displaystyle 1$
$\displaystyle=\eta_{\mu\nu}u^{\mu}u^{\nu}+\frac{\epsilon}{\tau_{D}-\tau_{P}}\int_{\tau_{P}}^{\tau_{D}}h_{\mu\nu}(\gamma_{0}(\tau))u^{\mu}u^{\nu}\;\mathrm{d}\tau$
(64) $\displaystyle\Leftrightarrow(u^{0})^{2}$
$\displaystyle\approx\vec{u}^{2}-\frac{\epsilon}{\tau_{D}-\tau_{P}}\int_{\tau_{P}}^{\tau_{D}}h_{\mu\nu}(\gamma_{0}(\tau))u^{\mu}u^{\nu}\;\mathrm{d}\tau$
(65)
$\displaystyle\Leftrightarrow\frac{(t_{D}-t_{P})^{2}}{(\tau_{D}-\tau_{P})^{2}}$
$\displaystyle\approx\frac{||\vec{X}_{D}-\vec{X}_{P}||^{2}}{(\tau_{D}-\tau_{P})^{2}}-\frac{\epsilon}{\tau_{D}-\tau_{P}}\int_{\tau_{P}}^{\tau_{D}}h_{\mu\nu}(\gamma_{0}(\tau))u^{\mu}u^{\nu}\;\mathrm{d}\tau$
(66) $\displaystyle\Leftrightarrow t_{D}-t_{P}$
$\displaystyle\approx\left[L_{0}^{2}-\epsilon(\tau_{D}-\tau_{P})\int_{\tau_{P}}^{\tau_{D}}h_{\mu\nu}(\gamma_{0}(\tau))u^{\mu}u^{\nu}\;\mathrm{d}\tau\right]^{\frac{1}{2}}$
(67) $\displaystyle\Leftrightarrow t_{D}-t_{P}$ $\displaystyle\approx
L_{0}\left[1-\epsilon(\tau_{D}-\tau_{P})\frac{u^{j}}{L_{0}}\frac{u^{k}}{L_{0}}\int_{\tau_{P}}^{\tau_{D}}h_{jk}(\gamma_{0}(\tau))\;\mathrm{d}\tau\right]^{\frac{1}{2}}$
(68) $\displaystyle\Leftrightarrow t_{D}-t_{P}$ $\displaystyle\approx
L_{0}-\frac{\epsilon}{2}L_{0}(\tau_{D}-\tau_{P})\frac{u^{j}}{L_{0}}\frac{u^{k}}{L_{0}}\int_{\tau_{P}}^{\tau_{D}}h_{jk}(\gamma_{0}(\tau))\;\mathrm{d}\tau\,.$
(69)

Here, we used that $h$ is given in the TT gauge in order to eliminate all
terms involving $h_{0\mu}$. Next, we exploit $L\approx t_{D}-t_{P}$ and
$u^{j}=(X_{D}^{j}-X_{P}^{j})/(\tau_{D}-\tau_{P})$ and define the unit vectors
$\hat{l}^{j}:=(X_{D}^{j}-X_{P}^{j})/L_{0}$ pointing into the direction of
neutrino flow. Then we can write the physical length as

$\displaystyle L(t_{P},t_{D})\approx
L_{0}-\frac{\epsilon}{2}\frac{L_{0}}{\tau_{D}-\tau_{P}}\hat{l}^{j}\hat{l}^{k}\int_{\tau_{P}}^{\tau_{D}}h_{jk}(\gamma_{0}(\tau))\;\mathrm{d}\tau\,.$
(70)

Furthermore neglecting terms of $\mathcal{O}(\epsilon^{2})$, we can write
$\frac{L_{0}}{\tau_{D}-\tau_{P}}\approx u^{0}$. Thus the most general
expression without further specifying $h_{\mu\nu}$ reads

$\displaystyle L(t_{P},t_{D})$ $\displaystyle\approx
L_{0}-\frac{\epsilon}{2}u^{0}\int_{\tau_{P}}^{\tau_{D}}h_{\parallel}(\gamma_{0}(\tau))\;\mathrm{d}\tau\,,$
(71) $\displaystyle h_{\parallel}(\gamma_{0}(\tau))$
$\displaystyle:=\hat{l}^{j}\hat{l}^{k}h_{jk}(\gamma_{0}(\tau))\,.$ (72)

We can make this formula more explicit by employing the Fourier decomposition
of $h_{\mu\nu}$ in plane wave modes which is given by

$\displaystyle h_{\mu\nu}(x)$
$\displaystyle=\sum_{r\in\\{+,\times\\}}\int\mathrm{d}^{3}\vec{k}\;\mathrm{Re}\left[\Psi_{r}(\vec{k})A_{\mu\nu}^{r}(\hat{k})\exp(-ikx)\right]\,,$
(73)

where $\Psi_{r}$ are the complex, scalar weight function (wave packets) of the
corresponding modes and $A_{\mu\nu}^{r}(\hat{k})$ are the normalized
polarization tensors depending only on the direction $\hat{k}$ of the
respective mode. Next we briefly discuss the polarization tensors and our
choice of coordinates before we carry on the calculation of the integrated GW
tensor along the neutrino trajectory.

A reasonable choice of coordinate system to describe the neutrino oscillation
experiment of interest would be to define the coordinate axes such that the
neutrino source resides at $\vec{X}_{P}\equiv 0$ and the detector remains at
$\vec{X}_{D}\equiv L_{0}\vec{e}_{z}$. This way we also assume that source and
detector are not bound to each other by any forces since the physical length
can vary freely as the GW passes through the system. Furthermore, we define
the propagation direction of the GW with respect to this coordinate system. In
order to obtain the parallel polarization components

$\displaystyle A_{\parallel}^{r}(\hat{k})$
$\displaystyle:=\hat{l}^{j}\hat{l}^{k}A_{lk}^{r}(\hat{k})\,,$ (74)

which we need for our length expression, we simply rotate from the coordinate
system, $\tilde{\Sigma}$, where $\hat{k}=\vec{e}_{z}$ to our system, $\Sigma$,
using one rotation around the $y$ axis with angle $\theta$ followed by a
rotation around the $z$ axis with angle $\varphi$. This yields

$\displaystyle A_{\parallel}^{r}(\hat{k})$
$\displaystyle=\hat{l}^{j}\hat{l}^{k}A_{jk}^{r}(\hat{k})$ (75)
$\displaystyle=R^{m}{}_{j}\hat{l}^{j}R^{n}{}_{k}\hat{l}^{k}\tilde{A}_{mn}^{r}(R^{T}\hat{k})$
(76) $\displaystyle=R^{m}{}_{3}R^{n}{}_{3}\tilde{A}_{mn}^{r}(\hat{e}_{z})$
(77) $\displaystyle=\begin{pmatrix}\cos(\varphi)\sin(\theta)\\\
\sin(\varphi)\sin(\theta)\\\
\cos(\theta)\end{pmatrix}^{T}\begin{pmatrix}A_{11}^{r}&A_{12}^{r}&0\\\
A_{12}^{r}&-A_{11}^{r}&0\\\ 0&0&0\\\
\end{pmatrix}\begin{pmatrix}\cos(\varphi)\sin(\theta)\\\
\sin(\varphi)\sin(\theta)\\\ \cos(\theta)\end{pmatrix}\,,$ (78)

where $A^{+}_{11}=1$, $A^{+}_{12}=0$ and $A^{\times}_{11}=0$,
$A^{\times}_{12}=1$. Multiplying everything out, we get

$\displaystyle A_{\parallel}^{+}(\hat{k})$
$\displaystyle=\sin^{2}(\theta)(\cos^{2}(\varphi)-\sin^{2}(\varphi))=\sin^{2}(\theta)\cos(2\varphi)\,,$
(79) $\displaystyle A_{\parallel}^{\times}(\hat{k})$
$\displaystyle=2\sin^{2}(\theta)\cos(\varphi)\sin(\varphi)=\sin^{2}(\theta)\sin(2\varphi)\,.$
(80)

The angle $\theta$ describes the angle between the zeroth order neutrino
trajectory and the GW propagation direction while $\varphi$ quantifies the
direction of the GW perpendicular to the neutrino trajectory.

Next, we decompose the complex wave functions $\Psi_{r}$ into their absolute
values and phases, i.e.

$\displaystyle\Psi_{r}(\vec{k})=:\psi_{r}(\vec{k})\exp(-i\phi^{r}(\vec{k}))\,,$
(81)

leading us to the final expression for $h_{\parallel}$:

$\displaystyle h_{\parallel}(x)$
$\displaystyle=\sum_{r\in\\{+,\times\\}}\int\mathrm{d}^{3}\vec{k}\;\psi_{r}(\vec{k})A_{\parallel}^{r}(\varphi,\theta)\cos(kx+\phi^{r}(\vec{k}))\,.$
(82)

Using this result, we can perform the integration along $\gamma_{0}$
explicitly and obtain a formula for $\Delta L$ in terms of the momentum
distribution of $h$:

$\displaystyle\Delta L(t_{P},t_{D})$
$\displaystyle=-\frac{\epsilon}{2}u^{0}\int_{\tau_{P}}^{\tau_{D}}h_{\parallel}(\gamma_{0}(\tau))\;\mathrm{d}\tau$
(83)
$\displaystyle=-\frac{\epsilon}{2}u^{0}\sum_{r\in\\{+,\times\\}}\int\mathrm{d}^{3}\vec{k}\;\psi_{r}(\vec{k})A_{\parallel}^{r}(\varphi,\theta)\int_{\tau_{P}}^{\tau_{D}}\cos(k\gamma^{0}+\phi^{r}(\vec{k}))\;\mathrm{d}\tau$
(84)
$\displaystyle=-\frac{\epsilon}{2}\sum_{r\in\\{+,\times\\}}\int\mathrm{d}^{3}\vec{k}\;\frac{u^{0}}{ku}\psi_{r}(\vec{k})A_{\parallel}^{r}(\varphi,\theta)\int_{k\gamma_{0}(\tau_{P})+\phi^{r}(\vec{k})}^{k\gamma_{0}(\tau_{D})+\phi^{r}(\vec{k})}\cos(\xi)\;\mathrm{d}\xi$
(85)
$\displaystyle\approx-\frac{\epsilon}{2}\sum_{r\in\\{+,\times\\}}\int\mathrm{d}^{3}\vec{k}\;\frac{\psi_{r}(\vec{k})A_{\parallel}^{r}(\varphi,\theta)}{\omega(1-\cos(\theta))}\int_{k\gamma_{0}(\tau_{P})+\phi^{r}(\vec{k})}^{k\gamma_{0}(\tau_{D})+\phi^{r}(\vec{k})}\cos(\xi)\;\mathrm{d}\xi$
(86)
$\displaystyle=-\frac{\epsilon}{2}\sum_{r\in\\{+,\times\\}}\int\mathrm{d}^{3}\vec{k}\;\frac{\psi_{r}(\vec{k})A_{\parallel}^{r}(\varphi,\theta)}{\omega(1-\cos(\theta))}\left(\sin(kX_{D}+\phi^{r}(\vec{k}))-\sin(kX_{P}+\phi^{r}(\vec{k}))\right)\,,$
(87)

where $\omega=||\vec{k}||$ and $\vec{u}\vec{k}\approx
u^{0}\omega\cos(\theta)$.

Again making use of our choice of coordinates, the ultrarelativistic limit and
neglecting $\mathcal{O}(\epsilon^{2})$ terms, we can write $X_{P}=(t,\vec{0})$
and $X_{D}=(t+L_{0},L_{0}\vec{e}_{z})$ for some time $t$ parameterizing the
trajectories of neutrino source and detector. Thus, from now on we will also
parametrize the distance deviation in terms of $t$, i.e.

$\displaystyle\Delta
L(t)\approx-\frac{\epsilon}{2}\sum_{r\in\\{+,\times\\}}\int\mathrm{d}^{3}\vec{k}\;\frac{\psi_{r}(\vec{k})A_{\parallel}^{r}(\varphi,\theta)}{\omega(1-\cos(\theta))}$
$\displaystyle\left(\sin(\omega(t+L_{0})-\omega
L_{0}\cos(\theta)+\phi^{r}(\vec{k}))-\sin(\omega t+\phi^{r}(\vec{k}))\right)$
(88)
$\displaystyle\begin{split}=-\frac{\epsilon}{2}\sum_{r\in\\{+,\times\\}}\int\mathrm{d}^{3}\vec{k}\;\frac{\psi_{r}(\vec{k})A_{\parallel}^{r}(\varphi,\theta)}{\omega(1-\cos(\theta))}&\left(\sin\left\\{\omega(1-\cos(\theta))L_{0}\right\\}\cos(\omega
t+\phi^{r}(\vec{k}))\right.\\\
&\quad+\left.\left[\cos\left\\{\omega(1-\cos(\theta))L_{0}\right\\}-1\right]\sin(\omega
t+\phi^{r}(\vec{k}))\right)\,.\end{split}$ (89)

Lastly, we want to consider the case of a sharply peaked momentum space wave
packet, i.e.

$\displaystyle\psi_{r}(\vec{k})$
$\displaystyle=a_{r}\delta_{\varepsilon}^{(3)}(\vec{k}-\vec{K})\,,$ (90)

where $\delta_{\varepsilon}^{(3)}$ represents a sequence of functions
converging to the delta distribution under the integral sign for
$\varepsilon\to 0$. In the case where all other functions than $\psi_{r}$ are
sufficiently weakly momentum dependent on the support of
$\delta_{\varepsilon}$ we can pull them out of the integral and evaluate them
at the peak momentum $\vec{K}$ which yields

$\Delta
L(t)=-\frac{\epsilon}{2}\sum_{r\in\\{+,\times\\}}\frac{a_{r}A_{\parallel}^{r}(\varphi,\theta)}{\omega_{K}(1-\cos(\theta))}\left(\sin\left\\{\omega_{K}(1-\cos(\theta))L_{0}\right\\}\cos(\omega_{K}t+\phi^{r}(\vec{K}))\right.\\\
\quad+\left.\left[\cos\left\\{\omega_{K}(1-\cos(\theta))L_{0}\right\\}-1\right]\sin(\omega_{K}t+\phi^{r}(\vec{K}))\right)\,.$
(91)

##  Appendix C The influence of non-trivial Neutrino Trajectories on the
Gravitational Wave Effect

Especially, when considering neutrinos in the solar system, the orientation of
the GW and neutrino trajectories might vary non-trivially during the time of
data taking. Usually this includes the motion of the neutrino detector and
source relatively to each other and the orbit of their center of mass around
the sun. Consequently this also has an non-negligible impact on the effect of
the GW on the neutrino flavor transition probability since the angle between
the GW direction $\vec{k}$ and the zeroth order neutrino trajectory direction
$\vec{v}$ might become time dependent. Figure 7 illustrates the scenario just
discussed.

SunRotation axisProductionDetection$\vec{k}_{GW}$$\vec{k}_{\nu}$

$\omega t$

Figure 7: Illustration of the rotation of the detction and production region
around the sun. The GW background is shown by the dashed gray lines with a
resulting angle $\Theta(t)$ between the direction of propagation of the GW
$\vec{k}_{\mathrm{GW}}$ and the neutrino propagation $\vec{k}_{\nu}$.

In the following, we assume that the relative motion of neutrino source and
detector is negligible during the propagation time of a single neutrino. This
will later allow us to repeat most of the steps of the derivation of $\Delta
L$ shown in appendix B since during the integration over the neutrino path we
can neglect the newly introduced time dependence. On the time scale of the
neutrino oscillation experiment, however, the relative orientation of the
neutrino velocity and GW vector as well as the GW tensor components may change
significantly.

Assuming that the GW direction remains fixed in the sun’s rest frame,
$\tilde{\Sigma}$ described by coordinates $\tilde{x}^{\mu}$, enables us to
calculate the now time dependent GW parameters in the co-rotating center of
mass system (CMS) of the neutrino source-detector system, $\Sigma$, via a
coordinate change to the new coordinates $x^{\mu}(\tilde{x}^{\alpha})$. The
corresponding GW tensors are, of course, related by the usual transformation
law

$\displaystyle h_{\alpha\beta}(x)$
$\displaystyle=\frac{\partial\tilde{x}^{\mu}}{\partial
x^{\alpha}}\frac{\partial\tilde{x}^{\nu}}{\partial
x^{\beta}}\tilde{h}_{\mu\nu}(\tilde{x}(x))\,,$ (92)

where the GW tensor in the sun’s rest frame, $\tilde{h}_{\mu\nu}$, is still
described by

$\displaystyle\tilde{h}_{\mu\nu}(\tilde{x})$
$\displaystyle=\sum_{r\in\\{+,\times\\}}\int\mathrm{d}^{3}\vec{k}\;\psi_{r}(\vec{k})A_{\mu\nu}^{r}(\hat{k})\cos(k\tilde{x}+\phi^{r}(\vec{k}))\,,$
(93)

c.f. equation (73), where we take $\Psi_{r}=\psi_{r}\exp(-i\phi^{r})$ with
$\psi_{r}\geq 0$.

In order to find the GW tensor in the CMS of the neutrino experiment, our task
now is to derive an expression for the coordinate transformation
$\tilde{x}^{\mu}(x^{\alpha})$. To that end, assume for simplicity that
neutrino source and detector move at a constant coordinate distance $L_{0}$,
as described in appendix B. Now, let $X_{P}(t)$ and $X_{D}(t)$ be the
trajectories of corresponding neutrino production and detection events,
respectively. We define absolute and relative coordinates in $\tilde{\Sigma}$
by the relations

$\displaystyle\begin{rcases}\tilde{D}^{\mu}(t):=\frac{\tilde{X}_{P}^{\mu}(t)+\tilde{X}_{D}^{\mu}(t)}{2}\\\
\tilde{r}^{\mu}(t):=\frac{\tilde{X}_{D}^{\mu}(t)-\tilde{X}_{P}^{\mu}(t)}{2}\end{rcases}\Leftrightarrow\begin{cases}\tilde{X}_{P}^{\mu}(t):=\tilde{D}^{\mu}(t)-\tilde{r}^{\mu}(t)\\\
\tilde{X}_{D}^{\mu}(t):=\tilde{D}^{\mu}(t)+\tilde{r}^{\mu}(t)\\\
\end{cases}\,,$ (94)

where $||\vec{\tilde{r}}(t)||\equiv L_{0}/2$ and $\tilde{r}^{0}(t)\equiv
L_{0}/2$ at zeroth order in the GW and using the ultrarelativistic
approximation. Thus, at most $\tilde{r}$ performs a rotating motion in
$\tilde{\Sigma}$. We specify the coordinate system $\Sigma$ by demanding that

  * •

$t=\tilde{t}$

  * •

$\vec{D}(t)\equiv 0$

  * •

$\vec{r}(t)\equiv\frac{L_{0}}{2}\vec{e}_{z}$

leading to the transformation law

$\displaystyle\begin{rcases}t(\tilde{x})=t\\\
\vec{x}(\tilde{x})=R^{T}(\vec{\Phi})R^{T}(\vec{\Theta}(t))\left(\vec{\tilde{x}}-\vec{\tilde{D}}(\tilde{t})\right)\end{rcases}\Leftrightarrow\begin{cases}\tilde{t}(x)=t\\\
\vec{\tilde{x}}(x)=\vec{\tilde{D}}(t)+R(\vec{\Theta}(t))R(\hat{\Phi})\vec{x}\end{cases}\,,$
(95)

where $\Theta(t)$ is the axis (and angle) of rotation about which
$\vec{\tilde{r}}$ rotates in $\tilde{\Sigma}$ and $\vec{\Phi}$ is the rotation
axis (and angle) needed to rotate $L_{0}\cdot\hat{\tilde{e}}_{z}/2$ to
$\vec{\tilde{r}}(0)$. From now on, we only refer to this combination of
rotations as $O(t):=R(\vec{\Theta}(t))R(\vec{\Phi})$. The components of the
Jacobian of this coordinate transformation thus read

$\displaystyle\frac{\partial\tilde{x}^{0}}{\partial x^{0}}$
$\displaystyle=1\,,$ $\displaystyle\frac{\partial\tilde{x}^{0}}{\partial
x^{k}}=0\,,$ $\displaystyle\frac{\partial\tilde{x}^{j}}{\partial
x^{0}}=\dot{\tilde{D}}^{j}(t)+\dot{O}^{j}{}_{l}(t)x^{l}\,,$
$\displaystyle\frac{\partial\tilde{x}^{j}}{\partial x^{k}}=O^{j}{}_{k}(t)\,.$
(96)

And the GW tensor becomes

$\displaystyle\begin{split}h_{\alpha\beta}(\tilde{x})&=\left(O^{\mu}{}_{\alpha}(t)+\delta^{\mu}(t)\delta_{\alpha
0}\right)\left(O^{\nu}{}_{\beta}(t)+\delta^{\nu}(t)\delta_{\beta 0}\right)\\\
&\quad\times\sum_{r\in\\{+,\times\\}}\int\mathrm{d}^{3}\vec{k}\;\psi_{r}(\vec{k})A_{\mu\nu}^{r}(\hat{k})\cos\left(\omega
t-\vec{k}^{T}O^{T}(t)(\vec{x}-\vec{\tilde{D}}(t))+\phi^{r}(\vec{k})\right)\end{split}$
(97)
$\displaystyle\begin{split}&=\sum_{r\in\\{+,\times\\}}\int\mathrm{d}^{3}\vec{k}\;\psi_{r}(\vec{k})O^{\mu}{}_{\alpha}(t)O^{\nu}{}_{\beta}(t)A_{\mu\nu}^{r}(\hat{k})\cos\left(\omega
t-(O(t)\vec{k})^{T}\vec{x}+(O(t)\vec{k})^{T}\vec{\tilde{D}}(t)+\phi^{r}(\vec{k})\right)\\\
&\quad+\delta_{\alpha
0}\sum_{r\in\\{+,\times\\}}\int\mathrm{d}^{3}\vec{k}\;\psi_{r}(\vec{k})\delta^{\mu}(t)O^{\nu}{}_{\beta}(t)A_{\mu\nu}^{r}(\hat{k})\cos\left(\omega
t-(O(t)\vec{k})^{T}\vec{x}+(O(t)\vec{k})^{T}\vec{\tilde{D}}(t)+\phi^{r}(\vec{k})\right)\\\
&\quad+\delta_{\beta
0}\sum_{r\in\\{+,\times\\}}\int\mathrm{d}^{3}\vec{k}\;\psi_{r}(\vec{k})O^{\mu}{}_{\alpha}(t)\delta^{\nu}(t)A_{\mu\nu}^{r}(\hat{k})\cos\left(\omega
t-(O(t)\vec{k})^{T}\vec{x}+(O(t)\vec{k})^{T}\vec{\tilde{D}}(t)+\phi^{r}(\vec{k})\right)\\\
&\quad+\delta_{\alpha 0}\delta_{\beta
0}\sum_{r\in\\{+,\times\\}}\int\mathrm{d}^{3}\vec{k}\;\psi_{r}(\vec{k})\delta^{\mu}(t)\delta^{\nu}(t)A_{\mu\nu}^{r}(\hat{k})\cos\left(\omega
t-(O(t)\vec{k})^{T}\vec{x}+(O(t)\vec{k})^{T}\vec{\tilde{D}}(t)+\phi^{r}(\vec{k})\right)\,,\end{split}$
(98)

where we define $\delta^{j}=\dot{\tilde{D}}^{j}(t)+\dot{O}^{j}{}_{l}(t)x^{l}$
and $\delta^{0}=0$. To simplify this expression, we perform an orthogonal
transformation of the integration variable $\vec{q}=O(t)\vec{k}$ which we then
rename to $\vec{k}$ again. Furthermore, we exploit the transformation behavior
of the polarization tensor
$A_{\mu\nu}(R\hat{k})={R}^{\alpha}{}_{\mu}{R}^{\beta}{}_{\nu}A_{\alpha\beta}(\hat{k})$
and define $\chi^{r}(\vec{k},t):=\phi^{r}(\vec{k})+\vec{k}\vec{\tilde{D}}(t)$.
This yields

$\displaystyle\begin{split}h_{\alpha\beta}(\tilde{x})&=\sum_{r\in\\{+,\times\\}}\int\mathrm{d}^{3}\vec{k}\;\psi_{r}(O^{T}\vec{k})A_{\alpha\beta}^{r}(\hat{k})\cos\left(\omega
t-\vec{k}\cdot\vec{x}+\chi^{r}(\vec{k},t)\right)\\\ &\quad+\delta_{\alpha
0}\sum_{r\in\\{+,\times\\}}\int\mathrm{d}^{3}\vec{k}\;\psi_{r}(\vec{k})\delta^{\mu}(t){O^{T}}^{\rho}{}_{\mu}A_{\rho\beta}^{r}(\hat{k})\cos\left(\omega
t-\vec{k}\cdot\vec{x}+\chi^{r}(\vec{k},t)\right)\\\ &\quad+\delta_{\beta
0}\sum_{r\in\\{+,\times\\}}\int\mathrm{d}^{3}\vec{k}\;\psi_{r}(\vec{k})\delta^{\nu}(t){O^{T}}^{\sigma}{}_{\nu}A_{\alpha\sigma}^{r}(\hat{k})\cos\left(\omega
t-\vec{k}\cdot\vec{x}+\chi^{r}(\vec{k},t)\right)\\\ &\quad+\delta_{\alpha
0}\delta_{\beta
0}\sum_{r\in\\{+,\times\\}}\int\mathrm{d}^{3}\vec{k}\;\psi_{r}(\vec{k})\delta^{\mu}(t){O^{T}}^{\rho}{}_{\mu}\delta^{\nu}(t){O^{T}}^{\sigma}{}_{\nu}A_{\sigma\rho}^{r}(\hat{k})\cos\left(\omega
t-\vec{k}\cdot\vec{x}+\chi^{r}(\vec{k},t)\right)\,.\end{split}$ (99)

One can immediately notice that this coordinate system does not belong to the
TT gauge orbit since $h_{0\mu}\neq 0$ in general. Thus in the most general
case we get corrections to $\Delta L$ from the non-inertial coordinate
transformation we just discussed. However, when considering sufficiently slow
motion of the coordinate system during the propagation of a single neutrino,
i.e.

$\displaystyle(t_{D}-t_{P})\cdot|\dot{\delta}^{j}(t)|\ll 1\,,$ (100)

we can actually neglect these terms in our expression in $\Delta L$ as being
of higher than linear order in small quantities. Thus in this approximation
the TT gauge condition holds again and we can proceed to calculate $\Delta L$
as before in appendix B, yielding

$\displaystyle\Delta L(t)$
$\displaystyle\approx-\frac{\epsilon}{2}\sum_{r\in\\{+,\times\\}}\int\mathrm{d}^{3}\vec{k}\;\frac{\psi_{r}(O^{T}(t)\vec{k})A_{\parallel}^{r}(\hat{k})}{\omega-\vec{k}\vec{v}}\left(\sin(kX_{D}(t)+\chi^{r}(\vec{k},t))-\sin(kX_{P}(t)+\chi^{r}(\vec{k}))\right)\,.$
(101)

In order to obtain a better impression on what changes compared to the case
discussed in appendix B, we now apply the inverse integral transform as before
leading to

$\displaystyle\Delta L(t)$
$\displaystyle\approx-\frac{\epsilon}{2}\sum_{r\in\\{+,\times\\}}\int\mathrm{d}^{3}\vec{k}\;\frac{\psi_{r}(\vec{k})A_{\parallel}^{r}(O(t)\hat{k})}{\omega-\vec{k}O(t)\vec{v}}\left(\sin(kO(t)X_{D}(t)+\chi^{\prime
r}(\vec{k},t))-\sin(kO(t)X_{P}(t)+\chi^{\prime r}(\vec{k}))\right)\,.$ (102)

Therefore, we see that the relative orientation of the neutrino trajectory and
the GW direction now changes with time, as expected, implying a variable
influence of the effect. If for example $\vec{k}$ becomes (anti-)parallel to
$\vec{v}$ the polarization vector vanishes and the corresponding mode doesn’t
contribute anymore.

##  Appendix D Toy Analysis: Likelihoods

In this section, we discuss the statistical procedure used to obtain the
exclusion limits for the GW parameter space. First, we describe the setup of a
toy experiment roughly modeling a typical neutrino oscillation experiment.
Afterwards, we describe likelihood ratio tests and how we use them in order to
obtain the allowed parameter space of our model for a given toy data set
generated assuming no GW background.

###  D.1 Experimental Data Model

We attempt to model the unfolded experimental neutrino event counts for each
flavor. To keep things simple since this is meant to be a proof of principle
analysis rather than a complete experimental treatment, we neglect the
following effects:

  * •

Matter effects

  * •

Systematic uncertainties like detector efficiency, accuracy and acceptance

  * •

The procedure of unfolding the measured data

  * •

Theoretical uncertainties

Furthermore, we do not simulate multiple data sets by sampling random event
counts from the corresponding Poisson distribution. Instead, we take the
underlying theory prediction for the expected number of events as the outcome
of our experiment. The details of the statistical methods employed in this
work will be described in the next subsection.

Regarding the experimental setup, we take into account the finite energy
resolution and the binning of neutrino counts. Moreover, experiments never
directly measure the flavor transition probabilities, but only in convolution
with the energy flux spectrum. The form of this spectrum mostly depends on the
origin of neutrinos, e.g. an energy power-law for atmospheric and
astrophysical neutrinos.

Therefore, we need to compute the expected number of events from the combined
probability distribution function (pdf):

$\rho_{b}(E,L)=\sum_{a}\varphi_{a}(E)P_{ab}(E,L_{0})\,,$ (103)

where $\varphi_{a}$ is the normalized neutrino flux for an initial flavor $a$
and $P_{ab}$ is the transition probability of a neutrino from flavor $a$ to
flavor $b$. Since we aim to predict the event count for each energy bin, we
need the accumulated probability for the respective energy ranges. Thus, for
the $i$-th energy bin and flavor $b$ we arrive at the final discrete
probabilities:

$p_{ib}(L_{0}):=\int_{E_{i}}^{E_{i+1}}\rho_{b}(E,L_{0})\;\text{d}E\,,$ (104)

where $E_{i}$ and $E_{i+1}$ are the lower and upper bin edges of the $i$-th
bin as illustrated in figure 8.

$E$Events$\nu_{e}$$\nu_{\mu}$$\nu_{\tau}$$E_{1}$$E_{2}$First energy
bin$\nu_{e}$$\nu_{\mu}$$\nu_{\tau}$$E_{n}$$E_{n+1}$$n\text{-th energy bin}$
Figure 8: Illustration of the number of events for each flavor in the
different energy bins.

Finally, we want to turn towards modeling the neutrino flux spectrum. We will
especially discuss the neutrino flux expected from natural neutrino sources,
i.e. atmospheric and astrophysical neutrinos, and that expected from fixed
target experiments. As already mentioned, in the simplest scenario a naturally
produced flux spectrum can be approximated by a power law energy dependence,
i.e.

$\varphi_{a}(E)=\frac{c_{a}}{\phi}E^{-\gamma}\,,$ (105)

where $\gamma>1$, $c_{a}$ determines the flavor composition of the source flux
and $\phi$ is a normalization constant such that

$1=\sum_{a}\int_{E_{\text{min}}}^{E_{\text{max}}}\varphi_{a}(E)\;\text{d}E\,.$
(106)

Here, $E_{\text{min}}$ and $E_{\text{max}}$ are the minimum and maximum
energies of neutrinos emerging from the source, respectively. Hence, in total
we arrive at

$\varphi_{a}(E)=(\gamma-1)\frac{c_{a}}{\sum_{b}c_{b}}\frac{E^{-\gamma}}{\left(E_{\text{min}}^{1-\gamma}-E_{\text{max}}^{1-\gamma}\right)}\,.$
(107)

A typical flavor composition of a neutrino source based on pion decay would
for example be $\vec{c}=(c_{a})^{3}_{a=1}=(1,2,0)$.

For a fixed target experiment the precise shape of the flux depends on the
details of the neutrino production mechanism. In order to keep our
considerations simple, we do not specify a particular production process but
approximate the neutrino energy spectrum by a Gaussian shape in order to
account for the fact that artificially produced neutrino spectra usually peak
at some mean energy $\bar{E}$ and exhibit some width $\sigma_{E}$. We adopt
this rather crude approximation, since our effect is not sensitive to the
exact shape of the energy spectrum. It is rather the location of the energy
peak determining if an effect is visible. Therefore, in the currently
considered setting the neutrino flux spectrum reads

$\displaystyle\varphi_{a}(E)$
$\displaystyle=\frac{c_{a}}{\phi}\exp{\left(-\frac{1}{2}\left(\frac{E-\bar{E}_{a}}{\sigma_{E,a}^{2}}\right)^{2}\right)}\,,$
(108)

where the normalization constant is given by

$\displaystyle\phi$
$\displaystyle=\sum_{a}c_{a}\sigma_{E,a}\sqrt{\frac{\pi}{2}}\left(\mathrm{erf}\left(\frac{E_{\mathrm{max}}-\bar{E}_{a}}{\sqrt{2}\sigma_{E,a}}\right)-\mathrm{erf}\left(\frac{E_{\mathrm{min}}-\bar{E}_{a}}{\sqrt{2}\sigma_{E,a}}\right)\right)\,.$
(109)

###  D.2 Likelihood Ratio Tests

In order to tell which parameter configurations $(h,\lambda)$ or $(h,f)$, i.e.
GW strain and wave length or frequency, lead to observable differences in
comparison to the scenario where no GW is present, we use a likelihood ratio
test. This method is especially suitable to tell differences between nested
theoretical models, i.e. models belonging to the same class differing only in
their respective parameter configurations. In our case, the standard scenario
is represented by the parameter configuration $h=0$.

In general the test statistic for this test is the negative logarithmic
likelihood ratio (NLLR) given by

$\displaystyle\Lambda(\vec{X})$
$\displaystyle=-2\ln\left(\frac{\sup_{\vec{\vartheta}\in\Theta_{0}}\mathcal{L}(\vec{\vartheta}\;|\;\vec{X})}{\sup_{\vec{\vartheta}\in\Theta}\mathcal{L}(\vec{\vartheta}\;|\;\vec{X})}\right)\in[0,\infty)\,,$
(110)

where $\vec{X}$ is a set of observations, $\vec{\vartheta}$ is the vector of
parameters, $\Theta_{0}$ is the space of parameter configurations
corresponding to the null hypothesis ($H_{0}$) and $\Theta$ represents the
space of all possible parameter configurations. If Wilk’s theorem applies,
$\Lambda$ follows a $\chi^{2}$ distribution with
$d:=\mathrm{dim}(\Theta)-\mathrm{dim}(\Theta_{0})$ degrees of freedom (dof).

In our scenario, $\Theta_{0}$ is a single point in our parameter space, i.e.
$\vec{\vartheta}_{0}=(h_{0},\lambda_{0})$, while $\Theta$ is the whole space
of $(h,\lambda)$ pairs. Moreover, all conditions of Wilk’s theorem are
fulfilled:

  * •

$\Theta_{0}\subset\Theta$: Since $\Theta_{0}=\\{\theta_{0}\\}$ consists of a
single point $\theta_{0}$ in $\Theta=\mathbb{R}\times(0,\infty)$

  * •

The optimum parameter configuration lies in the interior of $\Theta$:
$\theta_{\mathrm{max}}=(0,\lambda)$ where $\lambda\in(0,\infty)$ is arbitrary.

Thus, assuming sufficiently high event rates, $\Lambda$ follows a $\chi^{2}$
distribution with $2$ dof. As we consider counting experiments, our likelihood
function is given by a product of Poisson distributions, i.e.

$\mathcal{L}(h,\lambda\;|\;\vec{n})=\prod_{b=1}^{n_{f}}\prod_{i=1}^{n}\text{Pois}(n_{ib},\eta_{ib}(h,\lambda))\,,$
(111)

where $\eta_{ib}$ denotes the theoretical expected value of events in bin $i$
and flavor $b$:

$\eta_{ib}(h,\lambda)=p_{ib}(h,\lambda)N_{\text{tot}}\,.$ (112)

Here $N_{\mathrm{tot}}$ is the total number of neutrino events.

Subsequently, the NLLR of our model reads

$\Lambda(\vec{n})=-2\ln{\left(\frac{\mathcal{L}(h_{0},\lambda_{0})}{\sup_{(h,\lambda)\in\mathbb{R}^{2}}\left(\mathcal{L}(h,\lambda)\right)}\right)}\,,$
(113)

expanding to

$\displaystyle\Lambda(\vec{n})=-2\sum_{i,b}\left[n_{ib}\left(\ln{(\eta_{ib}(h_{0},\lambda_{0}))}-\ln{(\eta_{ib}(h^{\ast},\lambda^{\ast}))}\right)-\left(\eta_{ib}(h_{0},\lambda_{0})-\eta_{ib}(h^{\ast},\lambda^{\ast})\right)\right]\,,$
(114)

where $h^{\ast}$ and $\lambda^{\ast}$ are the parameters maximizing the
likelihood on $\Theta$.

For $H_{0}$ we assume that $(h_{0},\lambda_{0})$, i.e. a certain combination
of GW parameters, are the underlying parameters of the given data set. We then
proceed as follows: We generate a data set representing the expected number of
neutrino events if no GW is present. Then a grid of parameter configurations
$(h_{0},\lambda_{0})$ is scanned and the corresponding NLLR values are
compared to the $\alpha$-quantile value $\Lambda_{\alpha}$ obtained from the
$\chi^{2}$ distribution with $2$ dof. If $\Lambda(\vec{n})>\Lambda_{\alpha}$
the null hypothesis is rejected meaning that the GW parameters would lead to
an observable effect at $(1-\alpha)\times 100\,\%$ confidence level.

##  Appendix E Neutrinos from within the Solar System with Large Wave Packet
Widths

We observe that SM neutrinos with feasible wave packet widths can only provide
weak constraints. To investigate the maximal limits achievable for SM
neutrinos, we now consider the smallest still measurable energies possible in
typical neutrino experiments, i.e. $E=\mathcal{O}(\mathrm{keV})$. To estimate
the lowest possible strain $h$ in this scenario, we rearrange equation (16)
and obtain:

$h>\frac{8\pi rE}{L_{0}\Delta m_{jk}^{2}}\,.$ (115)

We are able to determine the wave packet width we need for achieving an
observable effect with the right side of equation (16). Table 3 presents the
values obtained in this process.

Table 3: Values of the three chosen frequency ranges $f_{\text{GW}}$, wavelengths $\lambda_{\text{GW}}$, baselines $L_{0}$, wave packet widths $\sigma_{x}$ and estimated strain $h$ sensitivity for $\mathrm{keV}$ neutrinos. For all benchmark points we assume $E\in[1,10]\,\mathrm{keV}$ with a bin width of $\Delta E\sim 0.45\,\mathrm{keV}$. $f_{\text{GW}}\,/\,$\mathrm{Hz}$$ | $\lambda_{\text{GW}}\,/\,$\mathrm{km}$$ | $L_{0}\,/\,$\mathrm{km}$$ | $\sigma_{x}\,/\,$\mathrm{m}$$ | $h$
---|---|---|---|---
$1$–$10$ |  ${10}^{7}$–${10}^{8}$ | $\thicksim${10}^{5}$$ | $4\text{\times}{10}^{-3}$ | ${10}^{-8}$
${10}^{-6}$–${10}^{-5}$ |  ${10}^{10}$–${10}^{11}$ | $\thicksim${10}^{8}$$ | $4$ | ${10}^{-11}$
${10}^{-9}$ | ${10}^{17}$ | $\thicksim${10}^{10}$$ | $4\text{\times}{10}^{2}$ | ${10}^{-13}$

The results of the analysis conducted with these parameters is depicted in
figure 9.

Figure 9: As figure 4, for wave packet widths from table 3.

We find that significantly stronger constraints as low as
$h\approx${10}^{-13}$$ with the longest baseline $L_{0,3}$ and the wavepacket
width $\sigma_{x}=$4\text{\times}{10}^{2}\text{\,}\mathrm{m}$$ can be achieved
in this manner. Interestingly this even results in an overlap with the
sensitivity range of EPTA.

##  Appendix F Alternative Scenarios

In this appendix, we discuss the scenarios of adding new heavy neutrino
flavours, meson mixing and an alternativ experimental approach.

#### New Heavy Neutrino Flavor

As we argued in the section 4.1 and 4.2, the sensitivity of neutrino
oscillations to the GW strain $h$ might improve for larger squared mass
differences, $\Delta m_{jk}^{2}$. Hence, in the following we want to turn our
attention to exotic models involving heavy sterile neutrinos. For simplicity
we consider the addition of a fourth neutrino generation mostly mixing with
the electron neutrino. Apart from that, we consider here the same experimental
scenarios already discussed in section 4.1. In spite of the freedom of
choosing $\Delta m_{41}^{2}$, we need to take into account mass dependent
upper limits on the mxing matrix element $|U_{\text{e}4}|^{2}$ [91, 92, 93,
94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110,
111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125,
126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140,
141, 142, 143, 144, 145, 146, 147, 148, 149, 150]. The most stringent
constraints are collected in figure 10 plotted against the mass of the fourth
mass eigenstate.

Figure 10: Shown are the bounds for the squared mixing matrix element
$|U_{e4}|^{2}$ as a function of fourth mass eigenstate $m_{4}$. The red
crosses indicate the mass and mixing configurations chosen for the benchmark
points in table 4 [91].

It can be seen that the constraints on the mixing angle are less stringent for
small masses $\mathcal{O}($\mathrm{eV}$)$ and for very high masses of
$\mathcal{O}($100\text{\,}\mathrm{GeV}$)$. It, therefore, seems intriguing to
consider these high masses, such as $\Delta
m_{41}^{2}=$100\text{\,}{\mathrm{GeV}}^{2}$$; however, assuming a wave packet
width of $\sigma_{x}=$1\text{\,}\mathrm{nm}$$, the condition from equation
(16) that the coherence length $L_{41}^{\mathrm{coh}}$ must be larger than the
considered baseline is only met for energies above the thresholds
$E_{\text{thr,1}}=$132\text{\,}\mathrm{PeV}$$,
$E_{\text{thr,2}}=$42\text{\,}\mathrm{EeV}$$ and
$E_{\text{thr,3}}=$420\text{\,}\mathrm{EeV}$$ for $L_{0,1}$, $L_{0,2}$ and
$L_{0,3}$, respectively. These threshold energies are far beyond the energies
achievable at colliders, so we do not consider them. At slightly smaller
masses of a few $\mathrm{keV}$ to $\mathrm{GeV}$, the mixing matrix element
$|U_{\text{e}4}|^{2}$ is already strongly constrained. Therefore, the
additional terms in the flavor transition probability describing the
oscillation over the large mass gap would be suppressed by
$|U_{e4}|^{2}\lesssim 10^{-8}$. Consequently, the difference to the standard
three neutrino case is nearly unobservable resulting in a very similar
sensitivity to the GW strain $h$. In the light of these considerations, we
concentrate on mass configurations of the order of $1$–$100\,$\mathrm{eV}$$
highlighted by the red crosses in figure 10.

Table 4: Values of the three chosen GW frequency ranges $f_{\text{GW}}$, wavelengths $\lambda_{\text{GW}}$, baselines $L_{0}$, various energies $E$, bin widths $\Delta E$, mass splittings $\Delta m_{41}^{2}$, mixing matrix elements $|U_{e4}|^{2}$ and threshold energies $E_{\text{thr}}$ where effects would be observable for wave packet widths $\sigma_{x}=$1\text{\,}\mathrm{nm}$$. $f_{\text{GW}}\,/\,$\mathrm{Hz}$$ | $\lambda_{\text{GW}}\,/\,$\mathrm{km}$$ | $L_{0}\,/\,$\mathrm{km}$$ | $E$ | $\Delta E$ | $\Delta m_{41}^{2}\,/\,$\mathrm{eV}$^{2}$ | $\qquad|U_{e4}|^{2}$ | $E_{\text{thr}}\,/\,$\mathrm{GeV}$$
---|---|---|---|---|---|---|---
$1$–$10$ | ${10}^{4}$–${10}^{5}$ | $\thicksim${10}^{3}$$ | $[10,100]$\mathrm{MeV}$$ | $4.5\text{\,}\mathrm{MeV}$ | $>$1$$ | $<$7.4\text{\times}{10}^{-4}$$ | $0.01$
$[1,10]$\mathrm{GeV}$$ | $0.45\text{\,}\mathrm{GeV}$ | $>$1$$ | $<$7.4\text{\times}{10}^{-4}$$ | $4.20$
${10}^{-6}$–${10}^{-5}$ | ${10}^{10}$–${10}^{11}$ | $\thicksim${10}^{8}$$ | $[1,10]$\mathrm{GeV}$$ | $0.45\text{\,}\mathrm{GeV}$ | $>$1$$ | $<$7.4\text{\times}{10}^{-4}$$ | $4.20$
$[1,10]$\mathrm{TeV}$$ | $0.45\text{\,}\mathrm{TeV}$ | $>$1$$ | $<$7.4\text{\times}{10}^{-4}$$ | $4.20$
$[1,10]$\mathrm{TeV}$$ | $0.45\text{\,}\mathrm{TeV}$ | $>$5.6\text{\times}{10}^{3}$$ | $<${10}^{-6}$$ | $314.63$
${10}^{-9}$ | ${10}^{17}$ | $\thicksim${10}^{10}$$ | $[10,100]$\mathrm{GeV}$$ | $4.5\text{\,}\mathrm{GeV}$ | $>$1$$ | $<$7.4\text{\times}{10}^{-4}$$ | $42.04$
$[1,10]$\mathrm{TeV}$$ | $0.45\text{\,}\mathrm{TeV}$ | $>$1$$ | $<$7.4\text{\times}{10}^{-4}$$ | $42.04$
$[1,10]$\mathrm{TeV}$$ | $0.45\text{\,}\mathrm{TeV}$ | $>$5.6\text{\times}{10}^{3}$$ | $<${10}^{-6}$$ | $3146.34$

In table 4 we show our choice of parameter configurations for the three
interesting distances. There also the corresponding threshold energy
$E_{\text{thr}}$, changing with the baseline $L_{0}$ and mass difference
$\Delta m_{41}^{2}$, is shown. The energy binning for each scenario is denoted
by $\Delta E$. We consider energy ranges around and slightly above the lowest
energy allowed to ensure $L_{41}^{\mathrm{coh}}>L_{0}$, c.f. criterion (ii) in
section 3. Furthermore, condition (i) implies that by going to even higher
energies the oscillation length increases implying weaker constraints on the
GW strain, $h$, as can be seen in equation (115). For example, for a baseline
of $L_{0,1}=$1000\text{\,}\mathrm{km}$$ and $\Delta
m_{41}^{2}=$1\text{\,}\mathrm{eV}$$ a minimum energy
$E_{\text{thr}}=$0.01\text{\,}\mathrm{GeV}$$ is required which is why we
consider energies between $10\text{\,}\mathrm{MeV}$ and
$100\text{\,}\mathrm{MeV}$. The parameters for the baseline
$L_{0,2}=${10}^{8}\text{\,}\mathrm{km}$$ and
$L_{0,3}=${10}^{10}\text{\,}\mathrm{km}$$ follow the same reasoning. Due to
the small, additional mixing angle $\theta_{14}$ we need to increase the
number of events to $N=${10}^{9}$$ over 20 years in order to be able to
resolve the high frequency oscillation over the sterile-active mass gap.
Moreover, for the following analysis, we assume an initially pure electron
neutrino fixed-target flux and that conditions (iv) and (v) are fulfilled by
the neutrino source. Thus, the neutrino production region must fulfill
$\Sigma_{P}\ll 25\,\mathrm{m}$ for the lowest energies and $\Sigma_{P}\ll
442\,\mathrm{km}$ for the highest energies under consideration.

Firstly, we consider atmospheric-like distances of
$\mathcal{O}($1000\text{\,}\mathrm{km}$)$, the results of the respective
likelihood ratio test at 95% CL are depicted in figure 11.

Figure 11: As figure 4, for the scenarios from the first major row of table 4.

It can be observed that the contour line follows a similar shape to that of SM
neutrinos we discuss in section 4.1 and with this configuration, we can
achieve a minimal strain of $h\approx${10}^{-4}$$.

Figure 12: As figure 4, for the scenarios from the second major row of table
4.

The $95\,\%$ CL NLLR contours in the $f$–$h$ plane for baselines roughly of
the order of the Mars–Earth distance are depicted in figure 12. Besides the
typical behavior that the constraints for the strain $h$ weaken with
increasing neutrino energy, we also observe that the second $\mathrm{TeV}$
line, corresponding to the mass splitting $\Delta
m_{41}^{2}=$5.6\text{\times}{10}^{3}\text{\,}{\mathrm{eV}}^{2}$$, provides
weaker constraints due to the smaller mixing angle compared to the $\Delta
m_{41}^{2}=$1\text{\,}{\mathrm{eV}}^{2}$$ scenario. Hence, the best
constraints are achieved for a neutrino energy of $\mathrm{GeV}$, where it is
possible to probe strains of $h\approx${10}^{-7}$$.

Figure 13: As figure 4, for the scenarios from the last major row of table 4.

Lastly, we also consider baselines corresponding to the radius of the solar
system for heavy neutrinos. The results are shown in figure 13. Again, we
observe the same, typical behavior of the contours, with the difference that
this time the gap between the $\mathrm{TeV}$ contour lines is larger. This is
due to the longer baseline, as the energy resolution needs to increase in
order to probe the same oscillation lengths compared to smaller baselines.
With a neutrino energy of $\mathrm{GeV}$, neutrino experiments would be
sensitive to strains of $h\approx${10}^{-8}$$.

Overall, it can be concluded that heavy neutrinos do not significantly improve
the results for constraining GW parameters. The main reason for this is the