# Constraints on the Cosmic Neutrino Background from NGC 1068

Jack Franklin<EMAIL_ADDRESS>Ivan Martinez-Soler
<EMAIL_ADDRESS>Yuber F. Perez-Gonzalez  yuber.f.perez-
<EMAIL_ADDRESS>Jessica Turner<EMAIL_ADDRESS>Institute
for Particle Physics Phenomenology, Durham University, South Road DH1 3LE,
Durham, U.K.

###### Abstract

We use recent evidence of TeV neutrino events from NGC 1068, detected by the
IceCube experiment, to constrain the overdensity of relic neutrinos locally
and globally. Since these high-energy neutrinos have travelled long distances
through a sea of relic neutrinos, they could have undergone scattering,
altering their observed flux on Earth. Considering only Standard Model
interactions, we constrain the relic overdensity to be $\eta\leq 3.85\times
10^{8}~{}(5.39\times 10^{11})$ at the 95$\%$ confidence level for
overdensities with a radius of 14 Mpc (10 kpc), assuming the sum of neutrino
masses saturates the cosmological bound, $\sum_{i}m_{i}=0.13$ eV. We
demonstrate that this limit improves with larger neutrino masses and how it
depends on the scale of the overdensity region.

††preprint: IPPP/24/12

The $\Lambda$CDM standard cosmology model predicts a neutrino background akin
to the observed cosmic microwave background. The discovery of such relic
neutrinos would mark a significant breakthrough for cosmology and particle
physics, offering a path to unveil fundamental neutrino properties such as
their masses and whether they are Dirac or Majorana particles Long  _et al._
(2014). However, detecting these relic neutrinos experimentally has proven to
be remarkably challenging. This is primarily because these neutrinos have low
momenta, approximately $\sim 0.6~{}\mu$eV, resulting in only a few possible
interaction channels to detect these neutrinos. Nonetheless, there are
numerous proposals for relic neutrino detection including capture on
radioactive nuclei Weinberg (1962); Baracchini  _et al._ (2018), observing the
annihilation of ultra-high energy cosmic ray neutrinos with the relic neutrino
background using the $Z$-resonance Eberle  _et al._ (2004), elastic scattering
of the relic neutrino wind on a test mass Stodolsky (1975); Opher (1974); Duda
_et al._ (2001); Domcke and Spinrath (2017); Shvartsman  _et al._ (1982);
Smith and Lewin (1983); Cabibbo and Maiani (1982), looking for alterations in
atomic deexcitation spectra Yoshimura  _et al._ (2015) and resonant neutrino
capture in accelerator experiments Bauer and Shergold (2021).

These ideas, which may become experimentally feasible in the future, are
hindered by anticipated minuscule rates of detected neutrinos; see Ref. Bauer
and Shergold (2023a) for a comprehensive overview of experimental
sensitivities. This necessitates either large quantities of detector material
or the experimental detection of extremely subtle effects, rendering detection
exceptionally challenging. However, there is a caveat that could potentially
alter the anticipated rates in a future facility: while the $\Lambda$CDM model
predicts a neutrino density of approximately $\sim 56~{}\nu/{\rm cm^{3}}$ for
each mass state (and similarly for antineutrinos) Giunti and Kim (2007), there
could be an enhancement of the neutrino number density due to gravitational
clustering. Such an overdensity seems to be $\sim 20\%$ above the average
predicted by the $\Lambda$CDM model. Nonetheless, it could be much larger
depending on whether there are additional beyond-the-Standard Model (BSM)
interactions impacting neutrinos Ringwald and Wong (2004); de Salas  _et al._
(2017); Holm  _et al._ (2023); Zimmer  _et al._ (2023); Elbers  _et al._
(2023); Elbers (2023). Thus, it may be prudent to first experimentally
constrain potential overdensities before attempting a direct detection of the
Cosmic Neutrino Background (C$\nu$B).

Constraints on overdensities stem from both theoretical and experimental
sides. On the one hand, clustering large amounts of neutrinos would be
prohibited by the Pauli exclusion principle Tremaine and Gunn (1979); Bauer
and Shergold (2021). Since the standard C$\nu$B is expected to be close to the
Pauli limit, any large overdensity would indicate not only the presence of
additional interactions, but also that the relic neutrinos would possess a
larger temperature to the one expected in standard cosmology. On the other
hand, the KATRIN experiment recently has managed to constrain the local
overdensity through a direct search for events consistent with the detection
of the C$\nu$B, resulting in a limit of $\sim 10^{11}$ times the average
density predicted by cosmology Aker  _et al._ (2022). Another possible way to
constrain relic neutrino overdensities consists in determining the effects
from scattering of high energy particles with the C$\nu$B, such as detecting
C$\nu$B neutrinos boosted by the Diffuse Supernova Neutrino Background Das
_et al._ (2022), cosmic rays Císcar-Monsalvatje  _et al._ (2024), detecting
absorption features in the cosmogenic neutrinos spectrum Brdar  _et al._
(2022), or considering the gravitational effect on Solar System objects Tsai
_et al._ (2024).

A recent analysis of the high-energy events detected by IceCube has unveiled
the origins of some of these astrophysical neutrinos. The most significant
among the candidate sources is the active galaxy NGC 1068 Abbasi  _et al._
(2022); Aartsen  _et al._ (2020), with a global significance of $4.2\sigma$.
The discovery of these neutrino point sources presents an opportunity for
studying scenarios Beyond the Standard Model by exploring deviations in the
expected neutrino flux from the source, see e.g. Ng and Beacom (2014); Esteban
_et al._ (2021); Carloni  _et al._ (2022); Rink and Sen (2024); Cline and Puel
(2023); Döring and Vogl (2023). Crucially, these high-energy neutrinos
traverse extensive distances through the C$\nu$B on their way to Earth.
Overdensities of these relic neutrinos, which lie along the path traversed by
the astrophysical neutrinos, can potentially modify the astrophysical neutrino
fluxes through the self-interactions of neutrinos, even within the Standard
Model (SM). In this Letter, we seek to constrain the overabundance of relic
neutrinos by analysing neutrino events from NGC 1068 using publicly available
data from IceCube.

Relic Overdensity — Any deviation from the expected average number density of
relic neutrinos, $n_{0}\approx 56\,\mathrm{cm}^{-3}$, can be parameterised by
$\eta=n/n_{0}$, where $n$ is the true number density. Local relic
overdensities are predicted by standard $\Lambda$CDM cosmology since neutrinos
will cluster due to gravitational effects de Salas  _et al._ (2017); Holm  _et
al._ (2023). However, these effects are predicted to be small, with $\eta$
ranging between 1.2 and 20 Aker  _et al._ (2022). On the other hand, it is
possible that some exotic scenarios could modify this significantly, producing
larger overabundances over different scales Ringwald and Wong (2004); Tsai
_et al._ (2024); Smirnov and Xu (2022); Fardon  _et al._ (2004); Dvali and
Funcke (2016). If there is a significant overabundance of relic neutrinos,
interactions between relic and astrophysical neutrinos could impact the
neutrino flux from NGC 1068 observed at IceCube.

To take into account the effect of scattering with relic neutrinos of the flux
propagating to the Earth, it is necessary to solve a transport equation for
the flux of neutrinos with mass state $i$,

$\begin{split}\frac{\partial\Phi_{i}(r,E)}{\partial
r}=&-\Phi_{i}(r,E)\sum_{j}n_{j}\sigma_{ij}(E)\\\
&+\sum_{j,k,l}n_{k}\int_{E}^{\infty}dE^{\prime}\Phi_{j}(r,E^{\prime})\frac{d\sigma_{jk\rightarrow
il}(E^{\prime},E)}{dE}\,,\end{split}$ (1)

where $\Phi_{i}$ denotes the combined flux of neutrinos and anti-neutrinos
with mass state $i$, $r$ is the distance that the neutrinos have travelled
from NGC 1068, $n_{j}$ is the number density of mass states $j$ and $\sigma$
is the neutrino-neutrino cross-section. In the second term of Eq. (11), the
$j$ state is the incoming neutrino with energy $E^{\prime}$, the $k$ state is
the relic neutrino, $i$ is the outgoing neutrino with the mass state of
interest, and $l$ is the other neutrino state produced in the interaction. The
first term of Eq. (11) accounts for the decrease in neutrino flux caused by
scattering off of the C$\nu$B, while the second term accounts for both down-
scattering (energy loss of the neutrinos) and up-scattering of the relic
neutrino. The centre-of-mass energy of the interaction of the astrophysical
neutrino, with energy $E$, with the relic neutrino of mass $m_{j}$, assumed to
be at rest111Note that large overdensities due to some BSM interactions might
also imply that the relic neutrinos could be relativistic today, enhancing the
center-of-mass energy. However, we refrain from considering such a scenario to
keep our discussion as model-independent as possible., is

$\displaystyle\sqrt{s_{j}}$ $\displaystyle=\sqrt{2m_{j}E}$ $\displaystyle\sim
1.41~{}{\rm MeV}\left(\frac{m_{j}}{0.1~{}{\rm
eV}}\right)^{1/2}\left(\frac{E}{10~{}{\rm TeV}}\right)^{1/2}\,,$ (2)

and is always sufficiently small that the only important processes are
$\nu\nu\rightarrow\nu\nu$, $\nu\bar{\nu}\rightarrow\nu\bar{\nu}$, and
$\nu\bar{\nu}\rightarrow e^{+}e^{-}$.

For the neutrinos originating from NGC 1068, we assume an initial power-law
(PL) flux, where the flux is parameterised in terms of a normalisation
$\Phi_{0}$, at a reference energy $E_{0}$, and a spectral index $\gamma$ such
that $\Phi=\Phi_{0}(E/E_{0})^{-\gamma}$; throughout this Letter we take
$E_{0}=1$ TeV. We also assume that the neutrino flavours are produced with the
initial ratio of 1:2:0 for $\nu_{e}$:$\nu_{\mu}$:$\nu_{\tau}$, which
corresponds to pion decays Abdullahi and Denton (2020). To find the final flux
at Earth, we solve the transport equation numerically over the distance that
the neutrinos travel through the overdense region, $d_{\mathrm{eff}}$. For a
density that varies in value over spatial coordinates, for example, with
gravitational clustering de Salas  _et al._ (2017), $d_{\mathrm{eff}}$ is the
radius of a constant density profile that would produce the same number of
interactions. The effect of neutrino self-interactions over the astrophysical
neutrino flux can be seen in Fig. 1, where the initial flux is plotted
alongside the final flux for several different relic neutrino densities. The
main effect is a suppression of the flux at high energies. Regeneration
processes lead to an increased flux at lower energies; however, the influence
at lower energies remains minimal due to the rapid decrease of the flux with
energy.

Figure 1: Muon neutrino fluxes from NGC 1068 at IceCube under different
scenarios. The black line is the initial flux, taken to be a power law. The
parameters of the initial flux are the best-fit values from our TS analysis,
$\gamma=3.37$ and $\Phi_{0}=2.66\times
10^{-14}\mathrm{cm}^{-2}\mathrm{s}^{-1}\mathrm{GeV}^{-1}$. Given this initial
flux, the other lines are the solutions to the transport equation, with the
relevant neutrino number density indicated. The number density is assumed to
be constant over the distance between NGC 1068 and Earth, and $m_{1}$ (the
lightest neutrino mass) is $0.07$ eV.

Standard neutrino self-interactions — For the sake of simplicity, we limit our
discussion to SM neutrino interactions and assume that the number density is
the same for all mass states and for neutrinos/antineutrinos. We require both
total and differential SM cross-sections to solve the transport equations in
Eq. (11), which we will briefly outline next.

The total cross-section for the SM interactions between the neutrino flux from
NGC 1068 and relic neutrinos can be split into two contributions: one from the
production of neutrino final states only and one from the production of
electron-positron pairs:

$\sigma_{ij}=\sigma^{\nu}_{ij}+\sigma^{e}_{ij}\,,$ (3)

where $i$ is the incoming neutrino mass state, and $j$ is the relic neutrino
mass state, as before. From this, we find that

$\sigma^{\nu}_{ij}=\frac{G_{F}^{2}(7\delta_{ij}+2)s_{j}}{3\pi}\,,$ (4)

where $G_{F}$ is the Fermi constant. This cross-section has been summed over
the contributions from relic neutrinos and antineutrinos. For the cross-
section with $e^{+}e^{-}$ final states, the interaction must be calculated in
the mass basis, such that

$\sigma^{e}_{ij}=\sum_{\alpha}\left\lvert U_{ei}\right\rvert^{2}\left\lvert
U_{ej}\right\rvert^{2}\sigma^{e}_{\alpha}\,,$ (5)

where $\alpha$ is the flavour of the initial state neutrinos in the
interaction. The production of $e^{+}e^{-}$ only occurs with neutrino-
antineutrino annihilation, so we can consider the initial flavours to be
identical. The total cross-sections for this process are

$\displaystyle\sigma^{e}_{\\{\mu,\tau\\}}$
$\displaystyle=\frac{G_{F}^{2}\left((2\kappa-1)m_{e}^{2}+(\kappa+1)s_{j}\right)}{12\pi}\Theta(s_{j}-4m_{e}^{2})$
(6a) $\displaystyle\sigma^{e}_{e}$
$\displaystyle=\frac{G_{F}^{2}(s_{j}-3m_{e}^{2})(s_{j}-m_{e}^{2})}{3\pi
s_{j}}\Theta(s_{j}-4m_{e}^{2})$ (6b)

where $\kappa=4\sin^{2}\theta_{W}(2\sin^{2}\theta_{W}-1)$, with $\theta_{W}$
the weak mixing angle, $m_{e}$ the electron mass, and $\Theta(x)$ the
Heaviside step function. Unlike in the case of the total cross-sections, when
calculating the differential cross-section in the second term on the RHS of
Eq. (11), we are interested in the kinematics of the final state of the
interaction. In particular, we wish to obtain the differential cross-section
in terms of the energy of the outgoing $i$ mass state (anti)neutrino. We find
that the differential cross-section for neutrino pair production, taking all
processes into account, is

$\frac{d\sigma_{jk\rightarrow
il}(E^{\prime},E)}{dE}=\frac{G_{F}^{2}m_{k}}{2\pi}\left(A_{ijkl}+B_{ijkl}\left(\frac{E}{E^{\prime}}\right)^{2}\right)\,,$
(7)

where we have defined

$\displaystyle A_{ijkl}$
$\displaystyle=(\delta_{ij}\delta_{kl}+\delta_{ik}\delta_{jl})^{2}\,,$ (8a)
$\displaystyle B_{ijkl}$
$\displaystyle=(\delta_{jk}\delta_{il}+\delta_{ij}\delta_{kl})^{2}\,.$ (8b)

We provide details of this calculation in the Supplementary Material for the
interested reader.

Analysis — To search for C$\nu$B overdensities within the IceCube data, we
perform a maximum likelihood test Braun  _et al._ (2008); Abbasi  _et al._
(2022). The likelihood function for $N$ events, given a neutrino overdensity
$\eta$, is given by

$\mathcal{L}(\eta,n_{s}|\mathbf{x},N)=\prod_{i=1}^{N}\left(\frac{n_{s}}{N}f_{S}(\mathbf{x}_{i}|\eta)+\left(1-\frac{n_{s}}{N}\right)f_{B}(\mathbf{x}_{i})\right)\,,$
(9)

where $\mathbf{x}_{i}$ are the observables of the event $i$, $n_{s}$ is the
number of events associated with the signal, and $f_{S}$ and $f_{B}$ are the
signal and background probability distribution functions (pdf), respectively.
When performing the analysis, we consider the events from 2012-2018 taken from
the public release IceCube Collaboration (2021), and select those within a 15∘
radius of NGC 1068.

We compare the hypothesis of a neutrino overdensity $\eta$ with the standard
scenario by taking the log of the likelihood ratios

$TS=-2\Delta\log\mathcal{L}=-2\log\left(\frac{\mathcal{L}(\eta,ns|\mathbf{x}_{i},N)}{\mathcal{L}_{0}}\right)\,,$
(10)

for each value of the neutrino overdensity, we minimise the test statistic
($TS$) by marginalising over $n_{s}$ and $\gamma$. The likelihood of the null
hypothesis ($\mathcal{L}_{0}$) corresponds to the scenario where the C$\nu$B
density follows the prediction of the $\Lambda$CDM model and the neutrino flux
originated from NGC 1068 is described by $n_{s}=29.6$ and $\gamma=3.37$. These
values are obtained through a $TS$ analysis comparing the likelihoods between
the power-law model and the scenario where all the data corresponds to the
background.

The signal and background pdfs depend on the observables IceCube uses to
reconstruct the astrophysical events. These include the reconstructed energy
($\hat{E}_{\mu}$) and direction ($\hat{\bm{d_{i}}}$), which can be decomposed
into declination ($\hat{\delta}_{i}$) and right ascension
($\hat{\alpha}_{i}$), along with their uncertainties ($\hat{\sigma}_{i}$). As
the predominant backgrounds are atmospheric and diffuse astrophysical
neutrinos, the background pdf is independent of the right ascension. We can
thus construct the background pdf as

$f_{B}(\hat{E}_{\mu,i},\hat{\mathbf{d_{i}}},\hat{\sigma}_{i})=\frac{1}{2\pi}f_{B}(\hat{E}_{\mu,i},\sin\hat{\delta}_{i})\,.$
(11)

We neglect the dependence of $f_{B}$ on $\hat{\sigma}_{i}$ as it will be
approximately equal to the dependence of $f_{S}$ on $\hat{\sigma}_{i}$, and
thus will cancel in the $TS$ as can be seen in Eq. (10) Abbasi  _et al._
(2022). To estimate the background pdf, we use a data scrambling method where
we calculate the pdf of all events in the public data release IceCube
Collaboration (2021) that were detected with declination $\hat{\delta}_{i}$ in
a window of 15∘ above and below NGC 1068 Aartsen  _et al._ (2017). Doing so
effectively integrates over right ascension. Since the number of background
events far outweighs the number of signal events, this will approximate the
true pdf of the background events well Aartsen  _et al._ (2017).

The signal pdf can be decomposed into the product of the directional and the
energy pdfs Abbasi  _et al._ (2022)

$\begin{split}f_{S}(\hat{E}_{\mu,i},\hat{\bm{d}_{i}},\hat{\sigma}_{i}|\sin\delta_{\mathrm{src}},\bm{\theta})\approx&\frac{1}{2\pi\hat{\psi}_{i}}f_{S}(\hat{\psi}_{i}|\hat{E}_{\mu,i},\sigma_{i},\bm{\theta})\\\
&\times
f_{S}(\hat{E}_{\mu,i}|\sin\delta_{\mathrm{src}},\bm{\theta})\,,\end{split}$
(12)

where $\delta_{\mathrm{src}}$ is the declination of the assumed source and
$\bm{\theta}$ are the model parameters. For a PL flux, there is only one
relevant parameter, $\gamma$, whereas for the C$\nu$B scattered flux, there
are two parameters $\bm{\theta}=\\{\gamma,\eta\\}$. We approximate the
directional pdf,
$f_{S}(\hat{\psi}_{i}|\hat{E}_{\mu,i},\sigma_{i},\bm{\theta})$, as a Rayleigh
distribution, as was done in Braun  _et al._ (2008). This can be improved by
using better models of the detector response, for example, considering the
effect of energy on the spatial reconstruction Abbasi  _et al._ (2022). The
energy pdf in Eq. (12) is calculated from the following integral Bellenghi
_et al._ (2023)

$\begin{split}&f_{S}(\hat{E}_{\mu,i}|\sin\delta_{\mathrm{src}},\bm{\theta})=\\\
&\int
dE_{\nu}\,f(E_{\nu}|\sin\delta_{\mathrm{src}},\bm{\theta})f(\hat{E}_{\mu,i}|E_{\nu},\sin\delta_{\mathrm{src}})\,,\end{split}$
(13)

where $f(\hat{E}_{\mu}|E_{\nu},\sin\delta_{\mathrm{src}})$ is the pdf of
reconstructed muon energies given an initial neutrino energy $E_{\nu}$ coming
from a source with declination $\delta_{\mathrm{src}}$. This is taken from the
detector response functions. The neutrino energy pdf,
$f(E_{\nu}|\sin\delta_{\mathrm{src}},\bm{\theta})$, is calculated from the
flux and effective area. The detector response functions and the effective
area are taken from the public release from IceCube IceCube Collaboration
(2021).

Applying the maximum likelihood method outlined above, we utilised IceCube
data to assess the significance of the NGC 1068 source against a background-
only hypothesis to ensure the validity of our code, and we found a good
agreement with the official results. Details regarding the comparison of the
NGC 1068 source parameters obtained through our code with the IceCube official
results in Ref. Abbasi  _et al._ (2022) can be found in the Supplementary
Material.

Figure 2: The t-statistic $-2\Delta\log{\mathcal{L}}$ as a function of the
relic neutrino overabundance. The effect of the value of $m_{1}$ is also
demonstrated by taking different limits for the mass scale, as explained in
the text. Here we take $d_{\rm eff}=14$ Mpc for the radius of the overdense
region.

Results — In Fig. 2, we show the test statistic as a function of the C$\nu$B
overabundance parameter $\eta$. Since the SM cross-sections are proportional
to the centre-of-mass energy $s$, and because we assume the relic neutrinos to
be at rest, the values of $\eta$ that are probed will depend on the neutrino
masses. As the absolute scale of the neutrino masses is not known, this
analysis was repeated with different assumptions on the mass of the lightest
neutrino, assuming normal ordering (NO) and using the mixing best-fit
parameters from the NuFit global analysis Esteban  _et al._ (2020).

We consider three scenarios for the value of the lightest neutrino mass.
First, we employ the constraints on the sum of the neutrino masses coming from
cosmological measurements, i.e. $\sum_{i}m_{i}<0.13$ eV Abbott  _et al._
(2022), and the best-fit values of the mass splittings from neutrino
oscillation experiments to obtain a mass for the lightest neutrino of
$m_{1}=0.0342$ eV. Direct searches for neutrino masses, such as the one
carried away by KATRIN experiment, set bound on the effective electron
neutrino mass of $\left\lvert
U_{ei}\right\rvert^{2}m_{i}^{2}<0.8\,\mathrm{eV}^{2}$ The KATRIN Collaboration
(2022), which is then translated to the value of the lightest neutrino mass of
$m_{1}=0.8$ eV. Finally, we also consider a case where the lightest neutrino
mass is small, but $\nu_{1}$ is still non-relativistic today, i.e.
$m_{1}=0.001$ eV.

The strongest constraints on $\eta$ come from the larger values of $m_{1}$,
resulting from the increased centre-of-mass energy of the scattering
processes, which leads to stronger interactions between the neutrinos. As seen
in Fig. 2, there are some values of $\eta$ where the test statistic is
negative, indicating a potential preference for this value. However, the value
is small ($-2\Delta\log\mathcal{L}>-0.36$) and so well below being
statistically significant. From further inspection of the contributions to the
test statistic from each event, we find that the reduced statistics of higher
energy muon events allow for this preference. The limiting factor on the
strength of the constraint, i.e. the plateau observed at higher values of
$\eta$, is due to the limited strength of our analysis of NGC 1068 as a point
source. At these values, the signal pdf of the scattered neutrino model goes
to zero for all events as no signal events are predicted to be detected; this,
in turn, means that the model likelihood tends towards a pure background. To
push the exclusion bounds to higher confidence levels, we would require a
larger likelihood value for the PL point-source analysis of NGC 1068.

Extending this analysis to a range of values of $m_{1}$ gives the results
shown in Fig. 3, where the 90$\%$ and 95$\%$ C.L. are plotted for different
effective distances of the overabundance. A value of $d_{\mathrm{eff}}=14$ Mpc
implies that the relic density overabundance exists between Earth and NGC
1068, which could be the case if either there is some overdensity due to
clustering on large scales. Also shown are the bounds for more conservative
effective distances, corresponding to a local overdensity in the Milky Way or
possibly in NGC 1068. We compare the bounds from this analysis to those
produced by the KATRIN collaboration in Aker  _et al._ (2022), which limit
$\eta$ to less than $1.94\times 10^{11}\mathrm{cm}^{-3}$ for Dirac neutrinos
at the 90$\%$ C.L. Unlike the KATRIN bounds, our constraints improve as
$m_{1}$ increases due to the increase in the centre-of-mass energy. On the
other hand, for small $m_{1}$, the mass squared differences dominate in
setting the mass scale, which results in an asymptotic limit.

Figure 3:  The 90$\%$ confidence limits on a global relic overabundance as a
function of the mass of the lightest neutrino for the overdense region radius
of $d_{\rm eff}=$ 14 Mpc (full), 100 kpc (dashed), and 10 kpc (dot-dashed).
$d_{\mathrm{eff}}=14$ Mpc corresponds to an overdensity that extends between
NGC 1068 and Earth, but the limit also holds for larger scales. We assume
normal ordering and take the mass splittings from NuFIT 2022 Esteban  _et al._
(2020). The solid blue line corresponds to the limit on overdensities from the
Pauli exclusion principle with an escape velocity of
$\beta_{\mathrm{esc}}=1.8\times 10^{-3}$ Bauer and Shergold (2023b); for the
dashed blue line the same principle is used but we take the limit
$\beta_{\mathrm{esc}}\rightarrow 0$. The light purple line comes from the
upscattering of relic neutrinos by cosmic rays from Ref. Císcar-Monsalvatje
_et al._ (2024). Finally, the grey line illustrates values inconsistent with
the measured Milky Way’s mass. The green region indicates the KATRIN bound,
while the orange highlights constraints on the neutrino masses from Cosmology.

It is also possible that other BSM scenarios could be constrained using a
similar analysis. We performed a log-likelihood analysis on neutrino decays
arising from a coupling to a massless Majoron. However, we found that current
public data cannot constrain this scenario to a statistically significant
degree as $\left\lvert TS\right\rvert<0.2$ for all values of the coupling.
This may change with improved reconstruction and more data Valera  _et al._
(2023).

Final Thoughts — While there is a long way to go before experiments can begin
to probe overdensities close to the Pauli exclusion limit, which is
$\mathcal{O}(1)$ for a global overdensity Bauer and Shergold (2023b), the
results from this analysis offer a significant improvement on the direct
probing of the C$\nu$B without any assumption of BSM physics. For any
overabundance that spans between NGC 1068 and the Milky Way, this analysis
places a bounds on $\eta$ that reaches $3.48\times 10^{8}$ ($3.85\times
10^{8}$) at the 90$\%$ (95$\%$) CL. These bounds will improve as more data
from astrophysical sources are gathered. For example, the most recent result
from NGC 1068 Abbasi  _et al._ (2022), where the number of signal events has
increased to $n_{s}=79$, will enhance the exclusion overdensity to nearly
$5\sigma$. Additionally, as highlighted in Carloni  _et al._ (2022); Blanco
_et al._ (2023), other possible point sources could be added to this dataset,
thereby enhancing the potential to constrain the relic neutrino density
further. This outcome is expected to strengthen with upcoming experiments,
such as IceCube Gen2 Aartsen  _et al._ (2021), which is projected to increase
the statistics by a factor of eight. Additionally, background reconstruction
improvements will increase the significance of the limit on the neutrino
overdensities.

###  .1 Acknowledgements

###### Acknowledgements.

We want to thank Jack Shergold for useful discussions and Carlos A. Argüelles
and Martin Bauer for reading the draft version of this paper and for their
helpful comments. The UK Science and Technology Facilities Council (STFC) has
funded this work under grant ST/T001011/1. This project has received
funding/support from the European Union’s Horizon 2020 research and innovation
programme under the Marie Skłodowska-Curie grant agreement No 860881-HIDDeN.
This work has made use of the Hamilton HPC Service of Durham University. This
work used the DiRAC@Durham facility managed by the Institute for Computational
Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk), which is
part of the National eInfrastructure and funded by BEIS capital funding via
STFC capital grants ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham
University and STFC operations grant ST/R000832/1.

## References

  * Long  _et al._ (2014) A. J. Long, C. Lunardini, and E. Sabancilar, JCAP 08, 038 (2014), arXiv:1405.7654 [hep-ph] .
  * Weinberg (1962) S. Weinberg, Phys. Rev. 128, 1457 (1962).
  * Baracchini  _et al._ (2018) E. Baracchini _et al._ (PTOLEMY), (2018), arXiv:1808.01892 [physics.ins-det] .
  * Eberle  _et al._ (2004) B. Eberle, A. Ringwald, L. Song, and T. J. Weiler, Phys. Rev. D 70, 023007 (2004), arXiv:hep-ph/0401203 .
  * Stodolsky (1975) L. Stodolsky, Phys. Rev. Lett. 34, 110 (1975), [Erratum: Phys.Rev.Lett. 34, 508 (1975)].
  * Opher (1974) R. Opher, Astron. Astrophys. 37, 135 (1974).
  * Duda  _et al._ (2001) G. Duda, G. Gelmini, and S. Nussinov, Phys. Rev. D 64, 122001 (2001), arXiv:hep-ph/0107027 .
  * Domcke and Spinrath (2017) V. Domcke and M. Spinrath, JCAP 06, 055 (2017), arXiv:1703.08629 [astro-ph.CO] .
  * Shvartsman  _et al._ (1982) B. F. Shvartsman, V. B. Braginsky, S. S. Gershtein, Y. B. Zeldovich, and M. Y. Khlopov, JETP Lett. 36, 277 (1982).
  * Smith and Lewin (1983) P. F. Smith and J. D. Lewin, Phys. Lett. B 127, 185 (1983).
  * Cabibbo and Maiani (1982) N. Cabibbo and L. Maiani, Phys. Lett. B 114, 115 (1982).
  * Yoshimura  _et al._ (2015) M. Yoshimura, N. Sasao, and M. Tanaka, Phys. Rev. D 91, 063516 (2015), arXiv:1409.3648 [hep-ph] .
  * Bauer and Shergold (2021) M. Bauer and J. D. Shergold, Phys. Rev. D 104, 083039 (2021), arXiv:2104.12784 [hep-ph] .
  * Bauer and Shergold (2023a) M. Bauer and J. D. Shergold, JCAP 01, 003 (2023a), arXiv:2207.12413 [hep-ph] .
  * Giunti and Kim (2007) C. Giunti and C. W. Kim,  _Fundamentals of Neutrino Physics and Astrophysics_ (2007).
  * Ringwald and Wong (2004) A. Ringwald and Y. Y. Y. Wong, JCAP 12, 005 (2004), arXiv:hep-ph/0408241 .
  * de Salas  _et al._ (2017) P. F. de Salas, S. Gariazzo, J. Lesgourgues, and S. Pastor, JCAP 09, 034 (2017), arXiv:1706.09850 [astro-ph.CO] .
  * Holm  _et al._ (2023) E. B. Holm, I. M. Oldengott, and S. Zentarra, Phys. Lett. B 844, 138073 (2023), arXiv:2305.13379 [hep-ph] .
  * Zimmer  _et al._ (2023) F. Zimmer, C. A. Correa, and S. Ando, JCAP 11, 038 (2023), arXiv:2306.16444 [astro-ph.CO] .
  * Elbers  _et al._ (2023) W. Elbers, C. S. Frenk, A. Jenkins, B. Li, S. Pascoli, J. Jasche, G. Lavaux, and V. Springel, JCAP 10, 010 (2023), arXiv:2307.03191 [astro-ph.CO] .
  * Elbers (2023) W. H. Elbers,  _Neutrinos from horizon to sub-galactic scales_ , Ph.D. thesis, Durham U. (2023).
  * Tremaine and Gunn (1979) S. Tremaine and J. E. Gunn, Phys. Rev. Lett. 42, 407 (1979).
  * Aker  _et al._ (2022) M. Aker _et al._ (KATRIN), Phys. Rev. Lett. 129, 011806 (2022), arXiv:2202.04587 [nucl-ex] .
  * Das  _et al._ (2022) A. Das, Y. F. Perez-Gonzalez, and M. Sen, Phys. Rev. D 106, 095042 (2022), arXiv:2204.11885 [hep-ph] .
  * Císcar-Monsalvatje  _et al._ (2024) M. Císcar-Monsalvatje, G. Herrera, and I. M. Shoemaker, (2024), arXiv:2402.00985 [hep-ph] .
  * Brdar  _et al._ (2022) V. Brdar, P. S. B. Dev, R. Plestid, and A. Soni, Phys. Lett. B 833, 137358 (2022), arXiv:2207.02860 [hep-ph] .
  * Tsai  _et al._ (2024) Y.-D. Tsai, J. Eby, J. Arakawa, D. Farnocchia, and M. S. Safronova, JCAP 02, 029 (2024), arXiv:2210.03749 [hep-ph] .
  * Abbasi  _et al._ (2022) R. Abbasi, M. Ackermann, J. Adams, J. A. Aguilar, M. Ahlers, and M. Ahrens (IceCube Collaboration), Science 378, 538 (2022), https://www.science.org/doi/pdf/10.1126/science.abg3395 .
  * Aartsen  _et al._ (2020) M. G. Aartsen, M. Ackermann, J. Adams, J. A. Aguilar, M. Ahlers, M. Ahrens, C. Alispach, and et al, Phys. Rev. Lett. 124, 051103 (2020).
  * Ng and Beacom (2014) K. C. Y. Ng and J. F. Beacom, Phys. Rev. D 90, 065035 (2014).
  * Esteban  _et al._ (2021) I. Esteban, S. Pandey, V. Brdar, and J. F. Beacom, Phys. Rev. D 104, 123014 (2021).
  * Carloni  _et al._ (2022) K. Carloni, C. A. Arguelles, I. Martinez-Soler, K. S. Babu, and P. S. B. Dev (2022) arXiv:2212.00737 [astro-ph.HE] .
  * Rink and Sen (2024) T. Rink and M. Sen, Phys. Lett. B 851, 138558 (2024), arXiv:2211.16520 [hep-ph] .
  * Cline and Puel (2023) J. M. Cline and M. Puel, JCAP 06, 004 (2023), arXiv:2301.08756 [hep-ph] .
  * Döring and Vogl (2023) C. Döring and S. Vogl, (2023), arXiv:2304.08533 [hep-ph] .
  * Smirnov and Xu (2022) A. Y. Smirnov and X.-J. Xu, JHEP 08, 170 (2022), arXiv:2201.00939 [hep-ph] .
  * Fardon  _et al._ (2004) R. Fardon, A. E. Nelson, and N. Weiner, JCAP 10, 005 (2004), arXiv:astro-ph/0309800 .
  * Dvali and Funcke (2016) G. Dvali and L. Funcke, Phys. Rev. D 93, 113002 (2016), arXiv:1602.03191 [hep-ph] .
  * Abdullahi and Denton (2020) A. Abdullahi and P. B. Denton, Phys. Rev. D 102, 023018 (2020).
  * Braun  _et al._ (2008) J. Braun, J. Dumm, F. De Palma, C. Finley, A. Karle, and T. Montaruli, Astroparticle Physics 29, 299 (2008).
  * IceCube Collaboration (2021) IceCube Collaboration, “Icecube data for neutrino point-source searches years 2008-2018,”  (2021).
  * Aartsen  _et al._ (2017) M. G. Aartsen _et al._ (IceCube), Astrophys. J. 835, 151 (2017), arXiv:1609.04981 [astro-ph.HE] .
  * Bellenghi  _et al._ (2023) C. Bellenghi _et al._ (IceCube), PoS ICRC2023, 1061 (2023), arXiv:2308.12733 [astro-ph.HE] .
  * Esteban  _et al._ (2020) I. Esteban, M. Gonzalez-Garcia, M. Maltoni, T. Schwetz, and A. Zhou, Journal of High Energy Physics (2020), 10.1007/JHEP09(2020)178.
  * Abbott  _et al._ (2022) T. M. C. Abbott, M. Aguena, A. Alarcon, S. Allam, O. Alves, A. Amon, and et al (DES Collaboration), Phys. Rev. D 105, 023520 (2022).
  * The KATRIN Collaboration (2022) The KATRIN Collaboration, Nature Physics 18 (2022), 10.1038/s41567-021-01463-1.
  * Bauer and Shergold (2023b) M. Bauer and J. D. Shergold, Journal of Cosmology and Astroparticle Physics 2023, 003 (2023b).
  * Valera  _et al._ (2023) V. B. Valera, D. Fiorillo, I. Esteban, and M. Bustamante, PoS ICRC2023, 1066 (2023).
  * Blanco  _et al._ (2023) C. Blanco, D. Hooper, T. Linden, and E. Pinetti, (2023), arXiv:2307.03259 [astro-ph.HE] .
  * Aartsen  _et al._ (2021) M. G. Aartsen _et al._ (IceCube-Gen2), J. Phys. G 48, 060501 (2021), arXiv:2008.04323 [astro-ph.HE] .

Constraints on the Cosmic Neutrino Background from NGC1068 — Supplementary
Material

##  Appendix A I. SM INTERACTIONS

(a)  (b)  (c)  (d)  (e)

Figure 1: Feynman diagrams for the relevant interactions between neutrinos and
relic neutrinos. (a) and (b) are the $t$\- and $u$-channel diagrams,
respectively, for elastic scattering of incident neutrinos off of relic
neutrinos. (c) is for lepton pair production where $l=\nu$ or $l=e^{-}$. (d)
is the $t$-channel $e^{+}e^{-}$ pair production diagram. (e) is the
$t$-channel elastic scattering of neutrinos on relic antineutrinos.

###  A.1 A. Total cross-sections

We assume that all particles in the incoming flux are neutrinos rather than
antineutrinos. This is valid as the flux observed at IceCube is a combination
of $\nu_{\mu}$ and $\bar{\nu}_{\mu}$ fluxes, and the total cross-section is
invariant under swapping $\nu\leftrightarrow\bar{\nu}$.

We start with the $\nu\nu\rightarrow\nu\nu$ scattering. Since the $Z$-boson
mediates this process, it is convenient to work in the neutrino mass basis.
The total cross-section is:

$\sigma_{ij}=\frac{G_{F}^{2}s_{j}(3\delta_{ij}+1)}{2\pi}\,,$ (1)

where $G_{F}$ is Fermi’s constant and the $3\delta_{ij}$ occurs because of the
interference between $t$ and $u$-channel diagrams, as shown in Fig. 1a and
Fig. 1b respectively, which occurs when $i=j$. In the case of $\nu\bar{\nu}$
scattering, we separate the cross-section calculation into two categories -
$\nu\bar{\nu}$ and $e^{+}e^{-}$ production. The first of these follows
similarly to the previous case, in particular when $i\neq j$:

$\sigma_{ij}=\frac{G_{F}^{2}s_{i}}{6\pi}\,,$ (2)

while for $i=j$, we have to account for the annihilation and production of new
neutrino pairs. Combined, this gives:

$\sigma_{ii}=\frac{G_{F}^{2}s_{j}}{\pi}\left(\frac{2}{3}+2\times\frac{1}{6}\right)=\frac{G_{F}^{2}s_{j}}{\pi}\,,$
(3)

where the first term arises from production on a $\nu_{i}\bar{\nu}_{i}$ pair,
which receives an enhancement from the additional $t$-channel diagram as shown
in Fig. 1e. The second term is from the production of a $k$ state mass pair,
where $i\neq k$, of which there are two possibilities.

For the electron pair production process, the cross-section must be calculated
on a weak basis. We write the mass-basis cross-section in terms of the weak
basis cross-section using the PMNS matrix:

$\sigma_{ij}=\sum_{\alpha,\beta}\left\lvert
U_{i\alpha}\right\rvert^{2}\left\lvert
U_{j\beta}\right\rvert^{2}\sigma_{\alpha\beta}\,.$ (4)

However, we know that the flavours must be identical to produce an electron
pair. Focusing first on $\alpha\neq e$ we find that

$\begin{split}\sigma_{\alpha\alpha}=\frac{G_{F}^{2}}{12\pi}\left((2\kappa-1)m_{e}^{2}+(\kappa+1)s_{j}\right)\Theta(s_{j}-4m_{e}^{2})\,,\end{split}$
(5)

where $\kappa=4\sin^{2}\theta_{W}(2\sin^{2}\theta_{W}-1)$, with $\theta_{W}$
is the weak mixing angle, $m_{e}$ is the mass of the electron, and $\Theta(x)$
is the Heaviside step function. In the case of $\alpha=e$ there is an
additional $t$-channel, W boson mediated, interaction which leads to an
enhancement in the cross-section:

$\sigma_{ee}=\frac{G_{F}^{2}}{3\pi
s_{j}}(s_{j}-3m_{e}^{2})(s_{j}-m_{e}^{2})\Theta(s_{j}-4m_{e}^{2})\,.$ (6)

###  A.2 B. Differential cross-sections

For the process $\nu_{j}\nu_{k}\rightarrow\nu_{i}\nu_{l}$ we find that:

$\frac{d\sigma_{jk\rightarrow
il}}{dt}=\frac{G_{F}^{2}}{4\pi}(\delta_{ij}\delta_{kl}+\delta_{ik}\delta_{jl})^{2}\,,$
(7)

where the Mandelstam variables are $t=-2m_{k}(E^{\prime}-E)$ and $u=-2m_{k}E$.
The process $\nu_{j}\bar{\nu}_{k}\rightarrow\nu_{i}\bar{\nu}_{l}$ follows
similarly, with the Mandelstam variables remaining the same:

$\frac{d\sigma_{jk\rightarrow
il}}{dt}=\frac{G_{F}^{2}}{4\pi}\frac{u^{2}}{s_{k}^{2}}(\delta_{jk}\delta_{il}+\delta_{ij}\delta_{kl})^{2}\,.$
(8)

Finally, we also need to account for up-scattered relic antineutrinos, i.e.
the process $\nu_{j}\bar{\nu}_{k}\rightarrow\bar{\nu}_{i}\nu_{l}$. This
differs from the previous two cases as the Mandelstam variables $u$ and $t$
are swapped. The differential cross-section for this process is then:

$\frac{d\sigma_{jk\rightarrow
il}}{dt}=\frac{G_{F}^{2}}{4\pi}\frac{t^{2}}{s_{k}^{2}}(\delta_{jk}\delta_{il}+\delta_{ij}\delta_{kl})^{2}\,.$
(9)

Combining these different processes gives the final differential cross-
section:

$\frac{d\sigma_{jk\rightarrow
il}(E^{\prime},E)}{dE}=\frac{G_{F}^{2}m_{k}}{2\pi}\left((\delta_{ij}\delta_{kl}+\delta_{ik}\delta_{jl})^{2}+(\delta_{jk}\delta_{il}+\delta_{ij}\delta_{kl})^{2}\left(\frac{E}{E^{\prime}}\right)^{2}\right)\,,$
(10)

##  Appendix B II. NUMERICALLY SOLVING THE NEUTRINO TRANSPORT EQUATION

To produce the flux of muon neutrinos at Earth, it is necessary to solve the
following equation for the three mass states $i=1,2,3$:

$\frac{\partial\Phi_{i}(r,E)}{\partial
r}=-\Phi_{i}(r,E)\sum_{j}n_{j}\sigma_{ij}(E)+\sum_{j,k,l}n_{k}\int_{E}^{\infty}dE^{\prime}\
\Phi_{j}(r,E^{\prime})\frac{d\sigma_{jk\rightarrow il}}{dE}(E^{\prime},E)\,.$
(11)

We use a similar method to Esteban  _et al._ (2021) to discretise the fluxes
into bins of energy, and then solve over distance using an implicit finite
difference method. Doing so provides a numerically stable solution to the
transport equation.

###  B.1 A. Upper And Lower Bounds On Neutrino Energy

Before describing the details of the numerical method, we will explain the
bounds of the neutrino energy integral in the second term of Eq. 11. The lower
bound, somewhat trivially, is the lowest energy of neutrinos that can be
detected by IceCube, which we take to be $E_{\mathrm{min}}=100\,$GeV. We can
do this because, as seen from the integral in (11), the flux at a specific
energy depends only on the flux at higher energies. This results in the
neutrinos “flowing down” in energy. As such, there is no dependence on
undetectable neutrinos.

Physically, the upper bound on the neutrino energy is infinity; however, when
solving this integral numerically, we need a finite upper bound which will
approximate the integral well. To find a value for the finite upper bound, it
is useful to look at the integral of the initial power-law flux from this
lower bound up to some energy cutoff $E$:

$I(E)=\int_{100\,\mathrm{GeV}}^{E}\Phi(0,E^{\prime})dE^{\prime}=\frac{\Phi_{0}}{1-\gamma}\frac{1}{E_{0}^{-\gamma}}\left(E^{-\gamma+1}-\left(100\,\mathrm{GeV}\right)^{-\gamma+1}\right)\,,$
(12)

where the initial flux is a power law:

$\Phi(0,E)=\Phi_{0}\left(\frac{E}{E_{0}}\right)^{-\gamma}\,,$ (13)

and we have assumed that $\gamma>1$ to ensure that the integral converges when
$E\rightarrow\infty$. The fraction of the total flux above the cutoff is then
given by:

$f(E)=1-\frac{I(E)}{I(\infty)}=\left(\frac{E}{10^{2}\,\mathrm{GeV}}\right)^{-\gamma+1}\,.$
(14)

Since the total flux above our cutoff energy is always greater than or equal
to the scattered flux above the cutoff energy, this fraction is then an upper
bound on the error in the approximation of the integral in (11). If we want
the fraction of the total flux above our cutoff to be smaller than some
$\epsilon$, we can find the corresponding cutoff $E_{\mathrm{max}}$ from
rearranging (14). This gives the relation:

$\log_{10}\left(\frac{E_{\mathrm{max}}}{1\,\mathrm{GeV}}\right)=2+\frac{\log_{10}(\epsilon)}{\gamma-1}\,.$
(15)

For example, if we want to limit the flux ignored to $\epsilon=10^{-6}$ we
require $E_{\mathrm{max}}=10^{8}$ GeV assuming $\gamma=2$. In practice, we
take a fixed value of $E_{\mathrm{max}}=10^{10}$ GeV, which satisfies
$\epsilon\leq 10^{-6}$ for $\gamma\geq 1.5$.

###  B.2 B. Discretisation

To discretise the transport equation, we first rewrite it in a simpler
notation:

$\frac{1}{n_{\nu}}\frac{\partial\Phi_{i}(r,E)}{\partial
r}=-\Phi_{i}(r,E)K_{i}(E)+\sum_{j}\int_{E}^{\infty}dE^{\prime}\
\Phi_{j}(r,E^{\prime})J_{ji}(E,E^{\prime})$ (16)

where the function $K_{i}$ contains the cross-sections for incoming neutrino
with mass $m_{i}$, and the kernel function $J_{ji}$ contains the differential
cross-sections for the incoming neutrino with mass $j$ and outgoing neutrino
with mass state $i$. The sum over $k$ and $l$ has been moved inside of
$J_{j}i$. To discretise the flux over energy, we follow a similar procedure to
Esteban  _et al._ (2021) by integrating the differential equation over some
energy bin $(E_{n-1/2},E_{n+1/2})$. This improves the numerical stability of
the solution when there are discontinuities, such as those in the cross-
sections which arise from $e^{+}e^{-}$ production (in our implementation, we
space the bins logarithmically and take the total number of bins as 300). This
results in:

$\frac{\Delta E_{n}}{n_{\nu}}\frac{\partial\Phi_{i}^{n}(r)}{\partial
r}=-\Phi_{i}^{n}(r)K_{i}^{n}+\sum_{j}\sum_{m}\Phi_{j}^{m}J_{ji}^{mn}\,,$ (17)

where we have defined:

$\displaystyle\Phi_{i}^{n}(r)$ $\displaystyle=\frac{1}{\Delta
E_{n}}\int_{E_{n-1/2}}^{E_{n+1/2}}dE\ \Phi(r,E)$ $\displaystyle K_{i}^{n}$
$\displaystyle=\int_{E_{n-1/2}}^{E_{n+1/2}}dE\ K_{i}(r,E)$ $\displaystyle
J_{ji}^{mn}$
$\displaystyle=\int_{E_{n-1/2}}^{E_{n+1/2}}dE\int_{E_{m-1/2}}^{E_{m+1/2}}dE^{\prime}\
J_{ji}(E,E^{\prime}),$ $\displaystyle\text{ if }m>n$ $\displaystyle
J_{ji}^{mn}$
$\displaystyle=\int_{E_{n-1/2}}^{E_{n+1/2}}dE\int_{E}^{E_{m+1/2}}dE^{\prime}\
J_{ji}(E,E^{\prime}),$ $\displaystyle\text{ if }m=n$ $\displaystyle
J_{ji}^{mn}$ $\displaystyle=0,$ $\displaystyle\text{ if }m<n$
$\displaystyle\Delta E_{m}$ $\displaystyle=E_{m+1/2}-E_{m-1/2}\,.$

There are three possible values of $J_{ji}^{mn}$, each corresponding to
different cases. First, we have that the $m$th energy bin is higher than the
$n$th (or equivalently $m>n$). In this case we can integrate over both energy
limits independently. In the second case however, we have $m=n$ which implies
that the energy bins are the same. Since the initial energy $E^{\prime}$ must
be greater than the final energy $E$, the lower energy limit in the second
integral is $E$ rather than $E_{m-1/2}$. Finally, we have the case of $m<n$
which is not possible as the energy must decrease, resulting in a value of
zero. We calculate the integrals of $K_{i}^{n}$ and $J_{ji}^{mn}$
analytically, which reduces the time needed to solve the equation.

We now discretise the distance $r$, following a half-step implicit scheme.
This amounts to the substitutions:

$\displaystyle\frac{\partial\Phi_{i}^{n}(r)}{\partial r}$
$\displaystyle\rightarrow\frac{\Phi_{i}^{n}(r_{a+1})-\Phi_{i}^{n}(r_{a})}{\Delta
r}$ (18) $\displaystyle\Phi_{i}^{n}(r)$
$\displaystyle\rightarrow\frac{1}{2}\left(\Phi_{i}^{n}(r_{a+1})+\Phi_{i}^{n}(r_{a})\right)$
(19)

Since $K$ and $J$ are not functions of $r$, this discretisation does not
affect them. Performing these substitutions in Eq. (17) gives:

$\frac{\Delta
E_{n}}{n_{\nu}}\frac{\Phi_{i}^{n}(r_{a+1})-\Phi_{i}^{n}(r_{a})}{\Delta
r}=-\frac{1}{2}\left(\Phi_{i}^{n}(r_{a+1})+\Phi_{i}^{n}(r_{a})\right)K_{i}^{n}+\sum_{j}\sum_{m}\frac{1}{2}\left(\Phi_{i}^{m}(r_{a+1})+\Phi_{i}^{m}(r_{a})\right)J_{ji}^{mn}\,.$
(20)

We now have a fully discretised form of Eq. (11), which we need to solve
starting from our initial power-law flux.

(a)  (b)

Figure 2: Contour plots in the number of signal events $n_{s}$ and the
spectral index $\gamma$ of the flux from NGC 1068. Left: The 68 and 95$\%$
contours from our power law analysis, neglecting the effect of scattering with
the C$\nu$B. Also shown are the best-fit point from our analysis (black
cross), which we compare to the best fit from SkyLLH Bellenghi  _et al._
(2023) using the full 2008-2018 public dataset (red cross) and the most recent
IceCube analysis Abbasi  _et al._ (2022) (blue cross). Right: A comparison of
the contours and best-fit parameters for the $\Lambda$CDM case and the case of
a large overdensity.

###  B.3 C. Solving The Discretised Equation

Eq. (20) can be reinterpreted as a matrix equation with vectors over energy
bins. The flux vector at $r=r_{a}$ is then $\bm{\Phi}_{i}^{a}$, whilst $K_{i}$
and $J_{ji}$ are now matrices over energy bins. Following this procedure, and
then rearranging the equation, we arrive at a solution which is defined
iteratively:

$\left(\frac{\bm{\Delta}_{E}}{n_{\nu}\Delta_{r}}+\frac{1}{2}K_{i}\right)\bm{\Phi}_{i}^{a+1}-\sum_{j}\frac{1}{2}J_{ji}\bm{\Phi}_{j}^{a+1}=\left(\frac{\bm{\Delta}_{E}}{n_{\nu}\Delta_{r}}-\frac{1}{2}K_{i}\right)\bm{\Phi}_{i}^{a}+\sum_{j}\frac{1}{2}J_{ji}\bm{\Phi}_{j}^{a}\,,$
(21)

where $\bm{\Delta}_{E}=\text{diag}(\Delta E_{0},\Delta E_{1},\dots,\Delta
E_{N})$ and $\Delta_{r}=\Delta r$. Since the RHS is just a vector, we can
solve the simultaneous set of equations for $i,j=1,2,3$ down to a linear
equation with just one unknown vector $\bm{\Phi}_{1}^{a+1}$; this is then
solved for and substituted back into the relations to find all the unknown
vectors. To solve each step, we use the Eigen linear algebra library eigenweb.
Since all the matrices are upper triangular, solving these equations is very
efficient.

One final step is to introduce a lower bound value for the flux in a bin,
below which it is set to zero. Doing so reduces noise and improves the
stability of the solution dramatically. The cutoff value was found by trial
and error not to lose any useful information about the flux. We found that a
value of $\Phi_{\mathrm{min}}=10^{-30}\,\Phi_{0}$ was sufficient, as it
provided a stable solution without truncating the flux at too low an energy.
Since the output of the numerical solver was used to calculate a pdf, the
normalisation $\Phi_{0}$ does not matter, so we set it to be 1.

##  Appendix C III. ANALYSIS RESULTS

In Fig. 2a, we compare the best-fit results from our analysis and those from
SkyLLH Bellenghi  _et al._ (2023) and the most recent analysis of NGC 1068
from IceCube Abbasi  _et al._ (2022). The SkyLLH analysis best-fit parameters
are within 2$\sigma$ of our best-fit points. This is expected as both analyses
use similar methods and data from the IceCube public release however we only
use the events from the full 86 string configuration whereas SkyLLH uses the
full 10 year dataset. In Fig. 2b, we also show how the contour for the initial
flux parameters is affected by the inclusion of relic neutrino scattering
effects. We have specifically chosen a neutrino density value corresponding to
a negative value of $TS$. The number of signal events is almost identical for
this value, $n=10^{8}\,\mathrm{cm}^{-3}$, but the spectral index has become
small. On the other hand, the 1$\sigma$ range almost contains the best-fit
parameter for the original best-fit, which shows that the spectral index may
play a sub-dominant role in the final likelihood. When investigating the
origin of the slight preference for particular values of the relic density, we
found that the lack of statistics at higher muon energies allowed for
ambiguities in the analysis. We believe these values will become constrained
with improved statistics at higher energies or at least have non-negative $TS$
values.


* Raissi et al. [2017] M. Raissi, P. Perdikaris, and G. E. Karniadakis.  Inferring solutions of differential equations using noisy multi-fidelity data.  _J. Comput. Phys._ , 335:736–746, 2017.
  * Hornik et al. [1989] K. Hornik, M. Stinchcombe, and W. Halbert.  Multilayer feedforward networks are universal approximators.  _Neural Netw._ , 2:359–366, 1989.
  * Zhang et al. [2019] D. Zhang, L. Guo, and G. E. Karniadakis.  Learning in modal space: solving time-dependent stochastic PDEs using physics-informed neural networks.  _SIAM J. Sci. Comput._ , 42:A639–A665, 2019.
  * Michoski et al. [2019] C. Michoski, M. Milosavljevic, T. Oliver, and D. Hatch.  Solving differential equations using deep neural networks.  _Neurocomputing_ , 399:193–212, 2019.
  * Mao et al. [2020] Z. Mao, A. D. Jagtap, and G. E. Karniadakis.  Physics-informed neural networks for high-speed flows.  _Comput. Method. Appl. M._ , 360:112789, 2020.
  * Duraisamy et al. [2019] K. Duraisamy, G. Iaccarino, and H. Xiao.  Turbulence modeling in the age of data.  _Annu. Rev. Fluid Mech._ , 51:357–377, 2019.
  * Xie et al. [2018] Y. Xie, E. Franz, M. Chu, and N. Thuerey.  tempogan: A temporally coherent, volumetric gan for super-resolution fluid flow.  _ACM Trans. Graph._ , 37(4):1–15, 2018.
  * Chu and Thuerey [2017] M. Chu and N. Thuerey.  Data-driven synthesis of smoke flows with cnn-based feature descriptors.  _ACM Trans. Graph._ , 36(4):1–14, 2017.
  * Anderson et al. [1996] J. Anderson, I. Kevrekidis, and R. Rico-Martinez.  A comparison of recurrent training algorithms for time series analysis and system identification.  _Comput. Chem. Eng._ , 20:S751–S756, 1996.
  * Crutchfield and McNamara [1987] James P Crutchfield and Bruce S McNamara.  Equations of motion from a data series.  _Complex Syst._ , 1(417-452):121, 1987.
  * Daniels and Nemenman [2015] Bryan C Daniels and Ilya Nemenman.  Automated adaptive inference of phenomenological dynamical models.  _Nat. Commun._ , 6(1):1–8, 2015.
  * Wang et al. [2017] J. Wang, J. Wu, and H. Xiao.  Physics-informed machine learning approach for reconstructing reynolds stress modeling discrepancies based on dns data.  _Phys. Rev. Fluids_ , 2(3):034603, 2017.
  * Hammond et al. [2022] J. Hammond, F. Montomoli, M. Pietropaoli, R. D. Sandberg, and V. M.  Machine learning for the development of data-driven turbulence closures in coolant systems.  _J. Turbomach._ , 144(8):081003, 2022.
  * Xu et al. [2022] X. Xu, F. Waschkowski, A. S. Ooi, and R. D. Sandberg.  Towards robust and accurate reynolds-averaged closures for natural convection via multi-objective cfd-driven machine learning.  _Int. J. Heat Mass Transf._ , 187:122557, 2022.
  * Mohan et al. [2020a] A. T. Mohan, N. Lubbers, D. Livescu, and M. Chertkov.  Embedding hard physical constraints in convolutional neural networks for 3D turbulence.  In _International Conference on Learning Representations_ , 2020a.
  * Yang et al. [2019] X. Yang, S. Zafar, J. Wang, and H. Xiao.  Predictive large-eddy-simulation wall modeling via physics-informed neural networks.  _Phys. Rev. Fluids_ , 4:034602, 2019.
  * Raissi et al. [2020] M. Raissi, A. Yazdani, and G. E. Karniadakis.  Hidden fluid mechanics: Learning velocity and pressure fields from flow visualizations.  _Science_ , 367(6481):1026–1030, 2020.
  * Belbute-Peres et al. [2020] F. Belbute-Peres, T. Economon, and Z. Kolter.  Combining differentiable pde solvers and graph neural networks for fluid flow prediction.  In _International Conference on Machine Learning_ , pages 2402–2411, 2020.
  * Lye et al. [2020] K. Lye, S. Mishra, and D. Ray.  Deep learning observables in computational fluid dynamics.  _J. Comput. Phys._ , 410:109339, 2020.
  * White et al. [2019] Cristina White, Daniela Ushizima, and Charbel Farhat.  Neural networks predict fluid dynamics solutions from tiny datasets.  _arXiv preprint arXiv:1902.00091_ , 2019.
  * Mohan et al. [2020b] Arvind T Mohan, Nicholas Lubbers, Daniel Livescu, and Michael Chertkov.  Embedding hard physical constraints in neural network coarse-graining of 3d turbulence.  _arXiv preprint arXiv:2002.00021_ , 2020b.
  * Kingma and Ba [2014] D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization.  _arXiv:1412.6980_ , 2014.
  * He et al. [2016] K. He, X. Zhang, S. Ren, and J. Sun.  Deep residual learning for image recognition.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognitionn_ , pages 770–778, 2016.
  * Chen et al. [2018] R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. Duvenaud.  Neural ordinary differential equations.  In _Conference on Neural Information Processing Systems_ , pages 6571–6583, 2018.
  * Pontryagin [2018] Lev Semenovich Pontryagin.  _Mathematical theory of optimal processes_.  Routledge, 2018.
  * Forest and Ruth [1990] E. Forest and R. D. Ruth.  Fourth-order symplectic integration.  _Physica D_ , 43:105–117, 1990.
  * Yoshida [1990] H. Yoshida.  Construction of higher order symplectic integrators.  _Phys. Lett. A_ , 150:262–268, 1990.
  * Candy and Rozmus [1991] J. Candy and W. Rozmus.  A symplectic integration algorithm for separable Hamiltonian functions.  _J. Comput. Phys._ , 92:230–256, 1991.
  * Tao [2016] Molei Tao.  Explicit symplectic approximation of nonseparable hamiltonians: Algorithm and long time performance.  _Physical Review E_ , 94(4):043303, 2016.
  * Wu et al. [2015] J. Z. Wu, H. Y. Ma, and M. D. Zhou.  _Vortical Flows_.  Springer, 2015.
  * Evans [2010] L. C. Evans.  _Partial Differential Equations_.  American Mathematical Society, 2 edition, 2010.
  * Roe [1981] P. L. Roe.  Approximate riemann solvers, parameter vectors and difference schemes.  _J. Comput. Phys._ , 43:357–372, 1981.
  * Helmholtz [1858] H. Helmholtz.  Uber integrale der hydrodynamischen Gleichungen welche den Wirbel-bewegungen ensprechen.  _J. Reine Angew. Math_ , 55:25–55, 1858.
  * Yang and Pullin [2010] Y. Yang and D. I. Pullin.  On Lagrangian and vortex-surface fields for flows with Taylor–Green and Kida–Pelz initial conditions.  _J. Fluid Mech._ , 661:446–481, 2010.
  * Xiong and Yang [2017] S. Xiong and Y. Yang.  The boundary-constraint method for constructing vortex-surface fields.  _J. Comput. Phys._ , 339:31–45, 2017.
  * Hao et al. [2019] J. Hao, S. Xiong, and Y. Yang.  Tracking vortex surfaces frozen in the virtual velocity in non-ideal flows.  _J. Fluid Mech._ , 863:513–544, 2019.
  * Cottet and Koumoutsakos [2000] G.H. Cottet and P.D. Koumoutsakos.  _Vortex Methods: Theory and Practice_.  Cambridge University Press, 2000.
  * Redmon et al. [2016] J. Redmon, S. Divvala, R. Girshick, and A. Farhadi.  You only look once: unified, real-time object detection.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognitionn_ , pages 779–788, 2016.
  * Lin et al. [2017] T. Lin, P. Goyal, R. Girshick, K. He, and P. Dollár.  Focal loss for dense object detection.  _IEEE Trans. Vis. Comput. Graph._ , pages 2980–2988, 2017.
  * Paszke et al. [2019] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, et al.  Pytorch: An imperative style, high-performance deep learning library.  In _Advances in Neural Information Processing Systems_ , pages 8026–8037, 2019.
  * Glorot and Bengio [2010] X. Glorot and Y. Bengio.  Understanding the difficulty of training deep feedforward neural networks.  In _Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics_ , pages 249–256, 2010.
  * Qu et al. [2019] Z. Qu, X. Zhang, M. Gao, C. Jiang, and B. Chen.  Efficient and conservative fluids using bidirectional mapping.  _ACM Trans. Graph._ , 38:1–12, 2019.
  * Sod [1978] G. A. Sod.  A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws.  _J. Comput. Phys._ , 27:1–31, 1978.
  * Lu et al. [2019] L. Lu, X. Meng, Z. Mao, and G. E. Karniadakis.  DeepXDE: A deep learning library for solving differential equations.  _SIAM Rev. Soc. Ind. Appl. Math._ , 63:208–228, 2019.
  * Rogallo [1981] R. S. Rogallo.  Numerical experiments in homogeneous turbulence.  In _Technical Report TM81315, NASA_ , 1981.
  * Xiong and Yang [2019] S. Xiong and Y. Yang.  Construction of knotted vortex tubes with the writhe-dependent helicity.  _Phys. Fluids_ , 31:047101, 2019.
  * Xiong and Yang [2020] S. Xiong and Y. Yang.  Effects of twist on the evolution of knotted magnetic flux tubes.  _J. Fluid Mech._ , 895:A28, 2020.
  * Yang et al. [2010] Y. Yang, D. I. Pullin, and I. Bermejo-Moreno.  Multi-scale geometric analysis of Lagrangian structures in isotropic turbulence.  _J. Fluid Mech._ , 654:233–270, 2010.
  * Yang and Pullin [2011] Y. Yang and D. I. Pullin.  Geometric study of Lagrangian and Eulerian structures in turbulent channel flow.  _J. Fluid Mech._ , 674:67–92, 2011.
  * Zhao et al. [2016] Y. Zhao, Y. Yang, and S. Chen.  Evolution of material surfaces in the temporal transition in channel flow.  _J. Fluid Mech._ , 793:840–876, 2016.
  * Zheng et al. [2016] W. Zheng, Y. Yang, and S. Chen.  Evolutionary geometry of Lagrangian structures in a transitional boundary layer.  _Phys. Fluids_ , 28:035110, 2016.
  * Zheng et al. [2019] W. Zheng, S. Ruan, Y. Yang, L. He, and S. Chen.  Image-based modelling of the skin-friction coefficient in compressible boundary-layer transition.  _J. Fluid. Mech._ , 875:1175–1203, 2019.
  * Lundgren [1982] T. S. Lundgren.  Strained spiral vortex model for turbulent fine structure.  _Phys. Fluids_ , 25:2193–2203, 1982.
  * Lundgren [1993] T. S. Lundgren.  A small-scale turbulence model.  _Phys. Fluids A_ , 5:1472, 1993.