††institutetext: Max-Planck-Institut für Kernphysik, Saupfercheckweg 1, 69117
Heidelberg, Germany

# Inferring astrophysical neutrino sources from the Glashow resonance

Guo-yuan Huang  Manfred Lindner  and Nele Volmer<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We infer the ultrahigh energy neutrino source by using the Glashow resonance
candidate event recently identified by the IceCube Observatory. For the
calculation of the cross section for the Glashow resonance, we incorporate
both the atomic Doppler broadening effect and initial state radiation
$\overline{\nu}_{e}e^{-}\to W^{-}\gamma$, which correct the original cross
section considerably. Using available experimental information, we have set a
generic constraint on the $\overline{\nu}_{e}$ fraction of astrophysical
neutrinos, which excludes the $\mu$-damped ${\rm p}\gamma$ source around
$2\sigma$ confidence level under the assumption that neutrino production is
dominated by the $\Delta$-resonance. While a weak preference has been found
for the pp source, next-generation measurements will be able to distinguish
between ideal pp and p$\gamma$ sources with a high significance assuming an
optimistic single power-law neutrino spectrum. The inclusion of multi-pion
production at very high energies for the neutrino source can weaken the
discrimination power. In this case additional multimessenger information is
needed to distinguish between pp and p$\gamma$ sources.

††arxiv: 2303.13706

##  1 Introduction

The IceCube Observatory has successfully established the observation of
ultrahigh energy (UHE) neutrino flux below a few PeV energies IceCube:2021rpz
; IceCube:2013low ; IceCube:2013cdw ; IceCube:2018cha ; IceCube:2020abv ;
IceCube:2020wum ; IceCube:2022der ; Halzen:2022pez . However, it remains a
mystery as to where those neutrinos come from. One of the most popular
mechanisms rests on accelerated cosmic rays colliding with ambient targets
around the source Gaisser:1994yf ; Bhattacharjee:1999mup ; Beatty:2009zz ;
gaisser_engel_resconi_2016 ; anchordoqui2019ultra ; tjus2020closing . There is
a variety of source models for UHE neutrinos Murase:2019tjj ; Murase:2022feu ;
Troitsky:2021nvu ; Xing:2011zza  which can usually be classified into the
p$\gamma$ and pp types depending on whether the target particle is a photon or
a proton.

For both p$\gamma$ and pp sources, after traveling an astronomical distance
the fluxes of three neutrino flavors strongly mix with each other due to
neutrino oscillations, which ends up with a nearly democratic flavor
composition
$\phi^{\oplus}_{\nu_{e}}+\phi^{\oplus}_{\overline{\nu}_{e}}:\phi^{\oplus}_{\nu_{\mu}}+\phi^{\oplus}_{\overline{\nu}_{\mu}}:\phi^{\oplus}_{\nu_{\tau}}+\phi^{\oplus}_{\overline{\nu}_{\tau}}\approx
1:1:1$ at Earth 444Throughout this work, we use the superscript ‘$\oplus$’ to
denote the quantity at Earth and ‘S’ to denote that at source.. It is unlikely
to disentangle those two sources by traditional flavor ratio measurements
Mena:2014sja ; Chen:2014gxa ; Palomares-Ruiz:2015mka ; Aartsen:2015ivb ;
Palladino:2015zua ; Arguelles:2015dca ; Bustamante:2015waa ; Aartsen:2015knd ;
Brdar:2016thq ; DAmico:2017dwq ; Pagliaroli:2015rca ; Rasmussen:2017ert ;
Brdar:2018tce ; Bustamante:2019sdb ; Palladino:2019pid ; Stachurska:2019srh ;
Song:2020nfh . The difference between those two sources lies in the
composition of neutrinos and antineutrinos. For the p$\gamma$ neutrino source,
cosmic rays collide with photons to produce charged pions (mostly $\pi^{+}$
below a certain energy threshold) followed by the decays
$\pi^{+}\to\mu^{+}\nu_{\mu}$ and $\mu^{+}\to
e^{+}\overline{\nu}_{\mu}\nu_{e}$, which results in more neutrino flux than
antineutrino flux, i.e., $\phi^{\rm S}_{\nu}:\phi^{\rm
S}_{\overline{\nu}}=2:1$. In comparison, the pp source will give rise to
nearly equal fractions of $\pi^{+}$ and $\pi^{-}$, which leads to $\phi^{\rm
S}_{\nu}:\phi^{\rm S}_{\overline{\nu}}=1:1$.

The key to distinguishing those two sources is by measuring the
$\overline{\nu}_{e}$ fraction
$f_{\overline{\nu}_{e}}\equiv\phi_{\overline{\nu}_{e}}/(\phi_{\overline{\nu}_{e}}+\phi_{{\nu}_{e}})$,
thanks to the Standard Model process $\overline{\nu}_{e}e^{-}\to W^{-}\to X$
predicted by S. L. Glashow Glashow:1960zz . Due to the resonance enhancement,
the cross section of $\overline{\nu}_{e}e^{-}$ scattering around
$E_{\nu}\approx 6.3~{}{\rm PeV}$ is larger than that of the deep inelastic
scattering (DIS) by more than two orders of magnitude. This promises us an
excellent channel to differentiate between the ideal pp (with
$f^{\oplus}_{\overline{\nu}_{e}}\approx 0.5$) and p$\gamma$ (with
$f^{\oplus}_{\overline{\nu}_{e}}\approx 0.23$) sources, as continuously
anticipated in previous works Berezinsky:1977sf ; Brown:1981ns ;
Anchordoqui:2004eb ; Hummer:2010ai ; Xing:2011zm ; Bhattacharya:2011qu ;
Bhattacharya:2012fh ; Barger:2012mz ; Barger:2014iua ; Palladino:2015uoa ;
Shoemaker:2015qul ; Anchordoqui:2016ewn ; Kistler:2016ask ; Biehl:2016psj ;
Sahu:2016qet ; Huang:2019hgs ; Zhou:2020oym . It is worthwhile to emphasize
that the above argument holds only for ideal pp and p$\gamma$ cases. Other
possible neutrino sources such as neutron and charm decays can give rise to a
$\overline{\nu}_{e}$ fraction different from the typical pion decays.
Moreover, if the multi-pion channel in the p$\gamma$ scattering dominates the
neutrino production when the collision energy is very high, more
multimessenger information about the source in addition to the $\nu_{e}$
fraction is necessary to disentangle the source degeneracy. In practice, the
overall diffuse neutrino flux might be contributed by different types of
sources, and the $\overline{\nu}_{e}$ fraction
$f^{\oplus}_{\overline{\nu}_{e}}$ can take any reasonable values in between.
In this regard, we shall first treat the $\overline{\nu}_{e}$ fraction as a
free parameter to be determined by the experimental probes regardless of the
model assumptions. The experimental information about this parameter can then
be used for theoretical interpretations under specific source assumptions.

Excitingly, with its unprecedented detection volume, the IceCube Observatory
has collected one candidate event with an energy deposition $E_{\rm
dep}=6.05\pm 0.72~{}{\rm PeV}$ in the sample of partially contained events
IceCube:2021rpz . The probability that this event stems from the Glashow
resonance (GR) is high, around $99\,\%$ by using the best-fit neutrino flux
taken from ref. IceCube:2015fuw . In this work, a timely quantitative
assessment is carried out to infer the $\overline{\nu}_{e}$ fraction by taking
$f^{\oplus}_{\overline{\nu}_{e}}$ as a free parameter and to explain the level
that we can differentiate between p$\gamma$ and pp sources. We have included
both the radiation of initial photons Gauld:2019pgt ; Garcia:2020jwr ;
Alikhanov:2009vcg  and the Doppler broadening effect Glashow:2014  while
calculating the GR events. Using the updated cross section, we investigate
both the results for the current GR candidate in IceCube as well as the
prospects of next-generation experiments.

##  2 A full treatment of Glashow resonance

As more and more UHE neutrino data have been accumulated, it becomes
increasingly important to take into account the subleading effects for the
theoretical evaluation of the GR. There are mainly two effects that should be
emphasized: (i) the initial state radiation (ISR) Gauld:2019pgt ;
Garcia:2020jwr ; (ii) the Doppler broadening effect Glashow:2014 . At the
leading level, the cross section for the process $\overline{\nu}_{e}e^{-}\to
W^{-}\to X$ reads IceCube:2021rpz

$\sigma^{(0)}(s)=24\pi{\rm\Gamma}^{2}_{W}{\rm
Br}_{W^{-}\to\overline{\nu}_{e}e^{-}}\frac{s/M_{W}^{2}}{(s-M_{W}^{2})^{2}+{\rm\Gamma}_{W}^{2}M_{W}^{2}}\;,$
(1)

where $M_{W}\approx 80.433\,\text{GeV}$ is the mass of the $W$ boson,
${\rm\Gamma}_{W}\approx 2.09~{}{\rm GeV}$ is the total decay width and ${\rm
Br}_{W^{-}\to\overline{\nu}_{e}e^{-}}\approx 10.7\,\%$ is the branching ratio
of the channel $W^{-}\to\overline{\nu}_{e}e^{-}$. The ISR and the Doppler
broadening effect are found to considerably modify the above picture and
should be included for completion.

Let us start with the ISR. This effect becomes increasingly notable when the
center-of-mass (COM) energy is much higher than the mass of the initial
charged lepton, for which the collinear emission of photons is significant.
For instance, in the Large Electron-Positron Collider (LEP), the ISR should be
taken into account when analyzing the $Z$ boson peak ALEPH:2005ab . For UHE
neutrino telescopes like IceCube, the ISR cross section near the GR will
receive a large enhancement factor of $\ln(M_{W}/m_{e})\approx 12$ on top of
the fine structure constant $\alpha$.

The ISR can be consistently included by using the structure function approach
in analogy with the DIS off hadrons. The modified cross section will be
Garcia:2020jwr

$\sigma(E_{\nu})=\int\mathrm{d}x\,{\rm\Gamma}_{e/e}(x,Q^{2})\sigma^{(0)}(x,Q^{2},E_{\nu})\;,$
(2)

where $Q$ represents the energy scale, $x$ is the longitudinal momentum
fraction of the electron after the photon radiation, $\sigma^{(0)}$ is the
cross section without the initial-state photon, and ${\rm\Gamma}_{e/e}$ is the
structure function of the electron. We take the structure function from ref.
Cacciari:1992  which includes soft photons resummed to all orders and hard
photons up to $\mathcal{O}(\alpha^{3})$.

Figure 1: Cross section for the Glashow resonance process
$\overline{\nu}_{e}+e^{-}\rightarrow W^{-}\rightarrow X$ with and without the
initial state radiation and Doppler broadening effect. The black curve shows
the cross section without initial state radiation and Doppler broadening, the
blue dotted one includes initial state radiation and the orange dotted one
includes Doppler broadening. The red curve is the cross section with both
Doppler broadening and initial state radiation effects, and the tabulated
result of this curve is given in our supplementary material. Both the
broadening and the radiative return are visible. For the Glashow resonance
curves we averaged over the electrons in H2O for the target.

The second effect of interest is the Doppler broadening due to the motion of
atomic electrons Glashow:2014 . The velocity of atomic electrons $\beta$ is
typically of the order $\mathcal{O}(\alpha\,c)$. A simple estimation shows
that this velocity will shift the COM energy square from $s=2E_{\nu}m_{e}$ to
$2E_{\nu}m_{e}(1-\beta\cos\theta)$, where $\theta$ is the angle between the
electron velocity and the incoming neutrino in the laboratory frame. This
broadens the COM energy by around $0.6~{}{\rm GeV}$ in comparison to the $W$
decay width ${\rm\Gamma}_{W}=~{}2.09~{}{\rm GeV}$. Non-relativistic electrons
in the atom have the four-momentum
$(m_{e}+{\left|\mathbold{k}\right|^{2}}/({2m_{e}}),\mathbold{k})$, where
$\left|\mathbold{k}\right|\approx m_{e}\beta$. By integrating over the
electron wave function, one can arrive at the total cross section Glashow:2014

$\sigma(E_{\nu})=\frac{1}{4\pi}\int\mathrm{d}\phi\int\mathrm{d}\beta
F(\beta)\int\mathrm{d}x^{\prime}\sigma^{(0)}[E_{\nu}(1-\beta x^{\prime})]\;,$
(3)

where $\phi$ represents the azimuth angle, $F(\beta)$ is the velocity
distribution of electrons and $x^{\prime}=\cos(\theta)$. Since the calculation
framework was already outlined in ref. Glashow:2014 , we give more details
about the updated calculation in appendix A.

Those two effects can be combined, and their joint result is shown as the red
curve in fig. 1 for the ${\rm H_{2}O}$ target, along with the cross sections
without (solid black curve) or modified by only one (blue and orange curves)
of those effects. In comparison, the charged-current (CC) and neutral-current
(NC) interactions are depicted as dashed and dotted black curves,
respectively. Some remarks on the results are given below.

  * •

The ISR will reduce the peak at the resonance energy $E_{\nu}\approx
6.3~{}{\rm PeV}$ by almost $20\,\%$. Furthermore, the cross section above the
resonance energy is enhanced by a factor of more than two. This is due to the
radiative return phenomenon, for which the photon in the process
$\overline{\nu}_{e}e^{-}\to W^{-}\gamma$ carries away some energy such that
the $W$ production will be made on shell even if $\sqrt{s}>M_{W}$.

  * •

The Doppler broadening effect for the ${\rm H_{2}O}$ target is small compared
to the ISR in the logarithmic scale. To see the detailed impact we also show
the result in a flat scale as fig. 4. The resonance peak is reduced slightly,
while the width is broadened due to the motion of atomic electrons.

  * •

The combined result of the ISR and the Doppler broadening is obtained with a
convolution, which reduces the peak by around $30\,\%$. However, we should
note that those effects will be partly smeared by the finite energy resolution
of the IceCube detector. We have checked that the eventual effect can decrease
the events within the energy window near the GR by almost $10\%$.

With the full GR cross section, we are able to calculate the event rate in
IceCube and compare it to both experimental data available now and those from
future experiments.

##  3 Analysis framework

In order to constrain the $\overline{\nu}_{e}$ fraction in the total diffuse
neutrino flux, we calculate the likelihood by fitting models with different
values of
$f_{\overline{\nu}_{e}}=\phi_{\overline{\nu}_{e}}/(\phi_{\overline{\nu}_{e}}+\phi_{{\nu}_{e}})$
to the available IceCube data. The reason why we use $f_{\overline{\nu}_{e}}$
to measure the $\overline{\nu}_{e}$ fraction is that it almost solely
determines the spectrum of single cascade event topology at PeV energies in
the IceCube detector.

The observed GR candidate in IceCube belongs to the PeV energy partially
contained events (PEPEs), in comparison to the high energy starting events
(HESEs) where the shower is fully contained inside the fiducial volume. Even
though the PEPE effective volume is nearly twice the volume of HESE at PeV
energies only one event with an energy deposition $E_{\rm dep}=6.05\pm
0.72~{}{\rm PeV}$ has been observed within the energy window $4~{}{\rm
PeV}<E_{\rm dep}<10~{}{\rm PeV}$. For HESE three PeV events have been
collected IceCube:2014stg , nicknamed Bert, Ernie and Big Bird. However, all
of them have energies below $3~{}{\rm PeV}$, which are most likely contributed
by the DIS. Even though the GR has not significantly arisen in the HESE
sample, HESE is useful to fix the normalization and shape of UHE neutrino flux
which are crucial for our extraction of the $\overline{\nu}_{e}$ fraction.

In ref. IceCube:2020wum , the IceCube Collaboration has analyzed the overall
UHE neutrino flux with HESEs collected over 7.5 years, assuming a flavor ratio
$\phi^{\oplus}_{\nu_{e}}+\phi^{\oplus}_{\overline{\nu}_{e}}:\phi^{\oplus}_{\nu_{\mu}}+\phi^{\oplus}_{\overline{\nu}_{\mu}}:\phi^{\oplus}_{\nu_{\tau}}+\phi^{\oplus}_{\overline{\nu}_{\tau}}=1:1:1$.
During our analysis we will use the HESE results including uncertainties from
ref. IceCube:2020wum  to set the spectrum of neutrino flux and use PEPE to
extract the $\overline{\nu}_{e}$ fraction $f^{\oplus}_{\overline{\nu}_{e}}$.
Note that a more thorough analysis would assume a completely free flavor
ratio. However, on the one hand, the latest IceCube HESE fit available has
fixed the flavor ratio IceCube:2020wum . On the other hand, ideal pp and
p$\gamma$ astrophysical models reasonably prefer such a democratic ratio after
neutrino oscillations over an astronomical distance.

For demonstration, we choose two benchmark flux models in our analysis: (i)
the unbroken single power-law model; (ii) the single power-law model with an
exponential energy cutoff. The former one reads

$\displaystyle\frac{\mathrm{d}{\rm\Phi}_{6\nu}}{\mathrm{d}E_{\nu}}={\rm\Phi}_{0}\left(\frac{E_{\nu}}{100~{}{\rm
TeV}}\right)^{-\gamma}10^{-18}~{}{\rm GeV^{-1}cm^{-2}s^{-1}sr^{-1}}\;,$ (4)

which represents models consistent with the Fermi acceleration mechanism and
extends to infinite energies. In practice, the reachable energy of
astrophysical accelerators always features a cutoff due to the Hillas
criterion Hillas:1984 . For the cutoff model, the flux in eq. (4) will be
multiplied by a suppression factor $\exp\left({-E_{\nu}/E_{\rm
cutoff}}\right)$. To confine the flux parameters, we construct a likelihood
based on the results in ref. IceCube:2020wum :

$\displaystyle-2\ln{\mathcal{L}_{6\nu}}=\frac{({\rm\Phi_{0}}-{\rm\Phi}^{\rm
bf}_{0})^{2}}{\sigma({\rm\Phi_{0}})^{2}}+\frac{(\gamma-\gamma^{\rm
bf})^{2}}{\sigma(\gamma)^{2}}\;,$ (5)

with the best-fit values ${\rm\Phi}^{\rm bf}_{0}=6.37$ and $\gamma^{\rm
bf}=2.87$, as well as the $1\sigma$ errors $\sigma({\rm\Phi_{0}})=1.54$ and
$\sigma(\gamma)=0.2$. For the cutoff model we further derive the likelihood
for $E_{\rm cutoff}$ from fig. VI.9 of ref. IceCube:2020wum  where the test-
statistic has been marginalized. Note that in this case we have ignored
possible correlations among ${\rm\Phi}_{0}$, $\gamma$ and $E_{\rm cutoff}$,
which are not provided. Nevertheless, such a choice will be more conservative
because less information is utilized in our analysis.

Figure 2: The likelihood (in blue) or posterior (in brown) of the
$\overline{\nu}_{e}$ fraction $f^{\oplus}_{\overline{\nu}_{e}}$ inferred from
the Glashow resonance event in IceCube with $4.6$ years of data taking. The
upper panel assumes a single power-law flux model with central values and
uncertainties from ref. IceCube:2020wum , while the lower one has incorporated
an exponential cutoff $E_{\rm cutoff}$ in the neutrino spectrum. The expected
$\overline{\nu}_{e}$ fractions of three representative ultrahigh energy
neutrino source models, including the ideal pp
($f^{\oplus}_{\overline{\nu}_{e}}\approx 0.5$), the ideal p$\gamma$
($f^{\oplus}_{\overline{\nu}_{e}}\approx 0.23$) and the ideal $\mu$-damped
p$\gamma$ ($f^{\oplus}_{\overline{\nu}_{e}}\approx 0$) sources, are indicated
by the dashed vertical lines. When the p$\gamma$ collision energy at the
source is very high, we may expect deviations from those ideal source models.
By assuming an equal mixture of single-pion and multi-pion production at the
source, we find the expected value $f^{\oplus}_{\overline{\nu}_{e}}\approx
0.36$ at Earth, shown as the dotted vertical lines. If multi-pion production
is more dominant than this assumption, the expected
$f^{\oplus}_{\overline{\nu}_{e}}$ should move to even larger values. The
sensitivity of the future IceCube-Gen2 project with an effective exposure of
ten (fifty) years is shown as the dashed (dotted) blue curves, assuming that
the pp source is dominant with ${\rm\Phi}_{0}=6.37$, $\gamma=2.7$ and $E_{\rm
cutoff}=5~{}{\rm PeV}$.

After the prior knowledge of $\\{{\rm\Phi}_{0},\gamma(,E_{\rm cutoff})\\}$ has
been established by HESE, we continue with fitting
$f^{\oplus}_{\overline{\nu}_{e}}$ to PEPE. The task is to calculate the
likelihood $\mathcal{L}_{\overline{\nu}_{e}}(f^{\oplus}_{\overline{\nu}_{e}})$
with the GR candidate we have. The joint likelihood can then be obtained with
$\mathcal{L}_{\rm
tot}=\mathcal{L}_{6\nu}\times\mathcal{L}_{\overline{\nu}_{e}}$ for the
parameter set ${\rm\Theta}\equiv\\{{\rm\Phi}_{0},\gamma(,E_{\rm
cutoff}),f^{\oplus}_{\overline{\nu}_{e}}\\}$. In the frame of extended
likelihood analysis of unbinned data Cowan:1998ji , the likelihood is
calculated with

$\displaystyle\mathcal{L}_{\overline{\nu}_{e}}=$
$\displaystyle\prod_{i=1}^{n}\left[\mu_{\rm DIS}{P}_{\rm
DIS}(\scalebox{0.8}{\raisebox{1.29167pt}{\char
35\relax}}{i}|{\rm\Theta})+\mu_{\rm GR}{P}_{\rm
GR}(\scalebox{0.8}{\raisebox{1.29167pt}{\char
35\relax}}{i}|{\rm\Theta})\right]\times\frac{1}{n!}\mathrm{e}^{-(\mu_{\rm
DIS}+\mu_{\rm GR})}\;,$ (6)

where $\mu_{\rm DIS}$ and $\mu_{\rm GR}$ are the expected event numbers within
the energy window $E_{\rm dep}\in[4,10]~{}{\rm PeV}$ for the DIS and the GR,
respectively, and $\scalebox{0.8}{\raisebox{1.29167pt}{\char 35\relax}}i$
represents in general all possible GR candidates. Moreover, ${P}_{\rm
DIS/GR}(\scalebox{0.8}{\raisebox{1.29167pt}{\char 35\relax}}{i}|{\rm\Theta})$
is the normalized probability to have an event at
$\scalebox{0.8}{\raisebox{1.29167pt}{\char 35\relax}}{i}$’s energy for the
given model parameter set ${\rm\Theta}$. Since there is only one GR candidate
so far we have $n=1$ in eq. (6). The event numbers can be obtained by
integrating the flux and cross sections with the detector configuration.

##  4 Main results

With the framework above, we can compute the total likelihood
$\mathcal{L}_{\rm tot}$ as a function of the parameter set
$\\{{\rm\Phi}_{0},\gamma(,E_{\rm cutoff}),f^{\oplus}_{\overline{\nu}_{e}}\\}$.
The likelihood can then be used for either frequentist or Bayesian
interpretations. For the frequentist interpretation, we obtain the likelihood
maximum $\mathcal{L}^{\rm max}_{\rm tot}(f^{\oplus}_{\overline{\nu}_{e}})$ by
marginalizing over the other parameters. For the Bayesian interpretation, we
need to derive the posterior distribution of $f^{\oplus}_{\overline{\nu}_{e}}$
by integrating over the likelihood and priors. We choose flat priors on
${\rm\Phi}_{0}$, $\gamma$, $f^{\oplus}_{\overline{\nu}_{e}}$ and $\ln{E_{\rm
cutoff}}$ for illustration.

Our main results are given in fig. 2, which shows the likelihood function (in
blue) or posterior distribution (in brown) of the $\overline{\nu}_{e}$
fraction $f^{\oplus}_{\overline{\nu}_{e}}$ inferred from the IceCube 4.6-year
data. Uncertainties from neutrino flux parameters have been systematically
included and marginalized when we constrain $f^{\oplus}_{\overline{\nu}_{e}}$.
The upper and lower panels stand for the assumptions of an unbroken single
power-law flux model and a single power-law model with a varying exponential
energy cutoff, respectively IceCube:2020wum . For blue curves, the horizontal
lines with $-2\ln{\mathcal{L}}=1$ and $4$ roughly set the $1\sigma$ and
$2\sigma$ confidence levels. For brown regions, the $1\sigma$ and $2\sigma$
credible intervals have been covered from dark to light colors.

We find that for all cases, the $\mu$-damped p$\gamma$ source with
$f^{\oplus}_{\overline{\nu}_{e}}\approx 0$ (single-pion production via the
$\Delta$-resonance for the ideal scenario) is excluded by around $2\sigma$
level. The current IceCube 4.6-year data weakly favor the pp source but are
not able to exclude the ideal p$\gamma$ source considerably (only at $1\sigma$
or so); see the dashed vertical lines. While interpreting the above results,
one must keep in mind that neutrinos may not only be produced by the ideal
$\Delta$-resonance of the p$\gamma$ scattering, but also by other possible
effects that can dominate at high energies Baerwald_2011 ; Huemmer_2010 , such
as multi-pion production, higher resonances, and the direct (t-channel)
production of pions. Note that the above considerations do not affect our
model-independent results of $f^{\oplus}_{\overline{\nu}_{e}}$ extracted from
experimental data. For those cases, the theoretically expected value of
$f^{\oplus}_{\overline{\nu}_{e}}$ for the $p\gamma$ source will shift towards
larger values. The actual magnitude of the deviation depends on the details of
the $\pi^{+}$ and $\pi^{-}$ mixture at the source. For demonstration, we
assume that the single-pion and multi-pion channels have the same production
rate at the source and draw the expected value
$f^{\oplus}_{\overline{\nu}_{e}}\approx 0.36$ as the dotted vertical lines in
fig. 2. If the multi-pion channel contributes more, this vertical line should
move even further to the right. On the other hand, for the pp source the
multi-pion contribution does not change the expected value of
$f^{\oplus}_{\overline{\nu}_{e}}$.

Figure 3: The log-likelihood of the energy cutoff $E_{\rm cutoff}$. The dashed
curve is taken from fig. VI.9 of ref. IceCube:2020wum , while the solid curve
is derived from the Glashow resonance candidate event by marginalizing over
the other model parameters.

Last but not least we should emphasize that the GR event can also constrain
the possible energy cutoff $E_{\rm cutoff}$ in the neutrino spectrum. The
original best-fit value of $E_{\rm cutoff}$ without GR is around $5~{}{\rm
PeV}$ in ref. IceCube:2020wum , with a $2\sigma$ lower boundary at $0.5~{}{\rm
PeV}$. The presence of the GR candidate event will push the $2\sigma$ lower
boundary to $2.2~{}{\rm PeV}$, as illustrated in fig. 3.

##  5 Outlook

Using the recent GR candidate event identified by IceCube, we have performed
an analysis to infer the $\overline{\nu}_{e}$ content in UHE astrophysical
neutrinos. We treat the $\overline{\nu}_{e}$ fraction as a model-independent
free parameter and have set a generic constraint on it by including the
uncertainties in the UHE neutrino flux. From the candidate event measured so
far, we find a weak preference for the pp source under the ideal assumption.
The situation will be greatly improved by the upcoming next-generation
neutrino telescopes. If the neutrino production in the p$\gamma$ source is
dominated by other channels at higher energies such as multi-pion production,
we would need additional information from the multi-wavelength observations of
the source to distinguish between pp and p$\gamma$ sources.

In the future, there are many projects such as IceCube-Gen2 IceCube:2014gqr ;
IceCube-Gen2:2020qha , Baikal-GVD Baikal-GVD:2018isr , KM3NeT KM3Net:2016zxf ,
P-ONE P-ONE:2020ljt , TAMBO Romero-Wolf:2020pzh , TRIDENT Ye:2022vbk  and so
on, which will provide very valuable sensitivities to PeV astrophysical
neutrinos Huang:2021mki ; Coleman:2022abf ; Ackermann:2022rqc ; Valera:2022wmu
. We take IceCube-Gen2 for demonstration by rescaling the current IceCube
target mass by ten times, and perform a count analysis in the energy window of
$[4,10]~{}{\rm PeV}$. The sensitivity for ten (fifty) years of effective
exposure is shown as the dashed (dotted) curves in fig. 2. Because the flux
parameters $\\{{\rm\Phi}_{0},\gamma(,E_{\rm cutoff})\\}$ can be very precisely
determined in the future IceCube-Gen2:2020qha , we choose a reasonably
optimistic spectrum as ${\rm\Phi}_{0}=6.37$, $\gamma=2.7$ and $E_{\rm
cutoff}=5~{}{\rm PeV}$ in making the forecast; see fig. 16 of ref. IceCube-
Gen2:2020qha  for example. It is worth noting that the tau neutrino telescope
TAMBO can be sensitive to the Glashow resonance event by searching for the
tau-induced showers from $W$ decays. Unfortunately, with a dedicated
simulation close to the TAMBO setup, we find that the event number is too
small, e.g., only two events with an optimistic flux and ten years of
exposure, compared to the DIS background of $\mathcal{O}(100)$. We will
elaborate on the related analysis in a future work.

Assuming the pp type as the true source, i.e.,
$f^{\oplus}_{\overline{\nu}_{e}}=0.5$, we expect eleven GR events in IceCube-
Gen2 with ten years of exposure for the optimistic single power-law model. If
we take an exponential cutoff $E_{\rm cutoff}=5~{}{\rm PeV}$ in the spectrum,
the event expectation would be reduced to three. The expected number of events
is still diverse due to low statistics of events at PeV energies. For the
single power-law model, IceCube-Gen2 with ten years of exposure can already
differentiate ideal pp from p$\gamma$ sources with a $2\sigma$ confidence
level. However, if there is an exponential cutoff at $5$ PeV, an effective
exposure of fifty years would be required to reach the $2\sigma$ level. Those
results can also be applied to other telescopes by adjusting the effective
exposure. By measuring the spectrum precisely in the future, one may go beyond
the assumptions of single power-law flux model (with cutoff) and take the
spectrum with a general energy dependence.

The hybrid cascade and early muon reconstruction in IceCube can already
greatly improve the angular resolution of the GR shower. In case of the
increased statistics, GR events detected in future experiments can also be
used to produce a map of the sky and identify associated PeVatrons
LHAASO:2021cbz ; LHAASO_nature ; Sudoh:2022sdk . Our main point is that
knowledge about neutrino sources will be significantly improved by those
upcoming facilities with large statistics, which also guarantees a robust
frontier for possible new physics studies Bustamante:2020niz ; Jezo:2014kla ;
Babu:2019vff ; Dey:2020fbx ; Babu:2022fje ; Xu:2022svm ; Arguelles:2022xxa ;
Huang:2022pce ; Huang:2022ebg ; Heighton:2023qpg .

##  Appendix A Appendix: Details of the Doppler broadening effect

Figure 4: Cross section for the Glashow resonance process
$\bar{\nu}_{e}+e^{-}\rightarrow W^{-}\rightarrow X$ with and without Doppler
broadening and assuming ice (${\rm H}_{2}{\rm O}$) as the target. The black
curve represents the cross section without Doppler broadening and for the
orange curve Doppler broadening is included.

We follow the procedure outlined in ref. Glashow:2014  to include the Doppler
broadening effect of atomic electrons. By integrating over angular variables
in eq. (3), we arrive at

$\displaystyle\sigma(E_{\nu})=\frac{6\pi{\rm\Gamma}_{W}^{2}{\rm
Br}_{W^{-}\rightarrow\overline{\nu}_{e}e^{-}}}{M_{W}m_{e}E_{\nu}}\int$
$\displaystyle
d\beta\frac{F(\beta)}{\beta}\left\\{\frac{1}{2M_{W}}\left[\ln(y_{h}^{2}+1)-\ln(y_{l}^{2}+1)\right]\right.$
$\displaystyle\left.+\frac{1}{{\rm\Gamma}_{W}}\left[\arctan(y_{h})-\arctan(y_{l})\right]\right\\}$
(7)

where

$y_{h}=\frac{2m_{e}E_{\nu}(1+\beta)+m_{e}^{2}-M_{W}^{2}}{{\rm\Gamma}_{W}M_{W}}\text{\hskip
14.22636pt and \hskip
14.22636pt}y_{l}=\frac{2m_{e}E_{\nu}(1-\beta)+m_{e}^{2}-M_{W}^{2}}{{\rm\Gamma}_{W}M_{W}}\,.$
(8)

Now the problem is attributed to the integration over the averaged electron
velocity distribution $F(\beta)$. In terms of the wave function of an electron
with quantum numbers $n$ and $l$, the distribution reads

$f_{n\,l}(\beta)=m_{e}\int\mathrm{d}{\rm\Omega}_{k}k^{2}|{\rm\Psi}_{n\,l}(k)|^{2}\text{\hskip
14.22636pt with \hskip 14.22636pt}{\rm\Psi}_{n\,l}(\mathbf{k})\propto
Y_{l\,m}^{*}({\rm\Omega}_{k})\int_{0}^{\infty}\mathrm{d}r\,r^{n+1}e^{-\mu
r}j_{l}(kr)\;,$ (9)

where $k=m_{e}\beta$ and $\mu_{n\,l}=\xi_{n\,l}/a_{0}$. Here, $a_{0}$ denotes
the Bohr radius, $\xi_{n\,l}=Z_{\rm eff}/n=(Z-\sigma_{n\,l})/n$, and
$\sigma_{n\,l}$ accounts for the screening of the nuclear charge by the other
electrons in the atom.

After the integration, we can get the velocity distribution for atoms up to
$Z=26$ Glashow:2014 :

$\displaystyle f_{1s}(k)=$
$\displaystyle\frac{32}{\pi}\frac{\mu_{1s}^{5}k^{2}}{(\mu_{1s}^{2}+k^{2})^{4}}\;,$
(10) $\displaystyle f_{2s}(k)=$
$\displaystyle\frac{32}{3\pi}\frac{\mu_{2s}^{5}(3\mu_{2s}^{2}k-k^{3})^{2}}{(\mu_{2s}^{2}+k^{2})^{6}}\;,$
(11) $\displaystyle f_{2p}(k)=$
$\displaystyle\frac{512}{3\pi}\frac{\mu_{2p}^{7}k^{4}}{(\mu_{2p}^{2}+k^{2})^{6}}\;,$
(12) $\displaystyle f_{3s}(k)=$
$\displaystyle\frac{1024}{5\pi}\frac{\mu_{3s}^{7}(\mu^{3}_{3s}k-\mu_{3s}k^{3})^{2}}{(\mu_{3s}^{2}+k^{2})^{8}}\;,$
(13) $\displaystyle f_{3p}(k)=$
$\displaystyle\frac{1024}{45\pi}\frac{\mu_{3p}^{7}(5\mu^{2}_{3p}k^{2}-k^{4})^{2}}{(\mu_{3p}^{2}+k^{2})^{8}}\;,$
(14) $\displaystyle f_{3d}(k)=$
$\displaystyle\frac{4096}{5\pi}\frac{\mu_{3d}^{9}k^{6}}{(\mu_{3d}^{2}+k^{2})^{8}}\;,$
(15) $\displaystyle f_{4s}(k)=$
$\displaystyle\frac{512}{35\pi}\frac{\mu_{4s}^{9}(5\mu^{4}_{4s}k-10\mu^{2}_{4s}k^{3}+k^{5})^{2}}{(\mu_{4s}^{2}+k^{2})^{10}}\;.$
(16)

Note that we have checked the expressions in ref. Glashow:2014  and corrected
possible discrepancies as in our eqs. (11) and (16).

We take the ice molecule ${\rm H}_{2}{\rm O}$ as an example. For oxygen,
$\mu_{1s}=7.6579$, $\mu_{2s}=2.2458$ and $\mu_{2p}=2.2266$ Clementi:1963 , and
for hydrogen $\mu_{1s}=1$. We weigh the distribution functions by averaging
over the electron numbers:

$F_{\rm ice}(\beta)=\frac{2F_{\rm H}(\beta)+8F_{\rm O}(\beta)}{10}\;.$ (17)

Using eq. (A) together with eq. (17) we get the Doppler broadened cross
section for ice as the target, which is depicted in fig. 4. The effect reduces
the peak by about $12\,\%$. Even though the total cross section integrated
over the initial neutrino energy is barely altered, the broadening effect will
make a difference when a non-uniform neutrino spectrum is considered.

###### Acknowledgements.

GYH would like to thank Xiao-Jun Bi, Xun-Jie Xu, Shoushan Zhang and Shun Zhou
for inspiring comments and discussions. GYH is supported in part by the
Alexander von Humboldt Foundation.

## References

  * (1) IceCube, M. G. Aartsen et al., “Detection of a particle shower at the Glashow resonance with IceCube,” Nature 591 (2021) no. 7849, 220–224, arXiv:2110.15051. [Erratum: Nature 592, E11 (2021)].
  * (2) IceCube, M. G. Aartsen et al., “Evidence for High-Energy Extraterrestrial Neutrinos at the IceCube Detector,” Science 342 (2013) 1242856, arXiv:1311.5238.
  * (3) IceCube, M. G. Aartsen et al., “First observation of PeV-energy neutrinos with IceCube,” Phys. Rev. Lett. 111 (2013) 021103, arXiv:1304.5356.
  * (4) IceCube, M. G. Aartsen et al., “Neutrino emission from the direction of the blazar TXS 0506+056 prior to the IceCube-170922A alert,” Science 361 (2018) no. 6398, 147–151, arXiv:1807.08794.
  * (5) IceCube, R. Abbasi et al., “Measurement of Astrophysical Tau Neutrinos in IceCube’s High-Energy Starting Events,” arXiv:2011.03561.
  * (6) IceCube, R. Abbasi et al., “The IceCube high-energy starting event sample: Description and flux characterization with 7.5 years of data,” Phys. Rev. D 104 (2021) 022002, arXiv:2011.03545.
  * (7) IceCube, R. Abbasi et al., “Evidence for neutrino emission from the nearby active galaxy NGC 1068,” Science 378 (2022) no. 6619, 538–543, arXiv:2211.09972.
  * (8) F. Halzen and A. Kheirandish, Chapter 5: IceCube and High-Energy Cosmic Neutrinos, pp. 107–235.  2023\.  arXiv:2202.00694.
  * (9) T. K. Gaisser, F. Halzen, and T. Stanev, “Particle astrophysics with high-energy neutrinos,” Phys. Rept. 258 (1995) 173–236, arXiv:hep-ph/9410384. [Erratum: Phys.Rept. 271, 355–356 (1996)].
  * (10) P. Bhattacharjee and G. Sigl, “Origin and propagation of extremely high-energy cosmic rays,” Phys. Rept. 327 (2000) 109–247, arXiv:astro-ph/9811011.
  * (11) J. J. Beatty and S. Westerhoff, “The Highest-Energy Cosmic Rays,” Ann. Rev. Nucl. Part. Sci. 59 (2009) 319–345.
  * (12) T. K. Gaisser, R. Engel, and E. Resconi, Cosmic Rays and Particle Physics.  Cambridge University Press, 2 ed., 2016.
  * (13) L. A. Anchordoqui, “Ultra-high-energy cosmic rays,” Physics Reports 801 (2019) 1–93.
  * (14) J. B. Tjus and L. Merten, “Closing in on the origin of Galactic cosmic rays using multimessenger information,” Physics Reports 872 (2020) 1–98.
  * (15) K. Murase and I. Bartos, “High-Energy Multimessenger Transient Astrophysics,” Ann. Rev. Nucl. Part. Sci. 69 (2019) 477–506, arXiv:1907.12506.
  * (16) K. Murase and F. W. Stecker, “High-Energy Neutrinos from Active Galactic Nuclei,” arXiv:2202.03381.
  * (17) S. Troitsky, “Constraints on models of the origin of high-energy astrophysical neutrinos,” Usp. Fiz. Nauk 191 (2021) no. 12, 1333–1360, arXiv:2112.09611.
  * (18) Z.-z. Xing and S. Zhou, Neutrinos in particle physics, astronomy and cosmology.  Zhejiang University Press, 2011.
  * (19) O. Mena, S. Palomares-Ruiz, and A. C. Vincent, “Flavor Composition of the High-Energy Neutrino Events in IceCube,” Phys. Rev. Lett. 113 (2014) 091103, arXiv:1404.0017.
  * (20) C.-Y. Chen, P. S. Bhupal Dev, and A. Soni, “Two-component flux explanation for the high energy neutrino events at IceCube,” Phys. Rev. D92 (2015) no. 7, 073001, arXiv:1411.5658.
  * (21) S. Palomares-Ruiz, A. C. Vincent, and O. Mena, “Spectral analysis of the high-energy IceCube neutrinos,” Phys. Rev. D91 (2015) no. 10, 103008, arXiv:1502.02649.
  * (22) IceCube, M. G. Aartsen et al., “Flavor Ratio of Astrophysical Neutrinos above 35 TeV in IceCube,” Phys. Rev. Lett. 114 (2015) no. 17, 171102, arXiv:1502.03376.
  * (23) A. Palladino, G. Pagliaroli, F. L. Villante, and F. Vissani, “What is the Flavor of the Cosmic Neutrinos Seen by IceCube?,” Phys. Rev. Lett. 114 (2015) no. 17, 171101, arXiv:1502.02923.
  * (24) C. A. Argüelles, T. Katori, and J. Salvado, “New Physics in Astrophysical Neutrino Flavor,” Phys. Rev. Lett. 115 (2015) 161303, arXiv:1506.02043.
  * (25) M. Bustamante, J. F. Beacom, and W. Winter, “Theoretically palatable flavor combinations of astrophysical neutrinos,” Phys. Rev. Lett. 115 (2015) no. 16, 161302, arXiv:1506.02645.
  * (26) IceCube, M. G. Aartsen et al., “A combined maximum-likelihood analysis of the high-energy astrophysical neutrino flux measured with IceCube,” Astrophys. J. 809 (2015) no. 1, 98, arXiv:1507.03991.
  * (27) V. Brdar, J. Kopp, and X.-P. Wang, “Sterile Neutrinos and Flavor Ratios in IceCube,” JCAP 01 (2017) 026, arXiv:1611.04598.
  * (28) G. D’Amico, “Flavor and energy inference for the high-energy IceCube neutrinos,” Astropart. Phys. 101 (2018) 8–16, arXiv:1712.04979.
  * (29) G. Pagliaroli, A. Palladino, F. L. Villante, and F. Vissani, “Testing nonradiative neutrino decay scenarios with IceCube data,” Phys. Rev. D 92 (2015) no. 11, 113008, arXiv:1506.02624.
  * (30) R. W. Rasmussen, L. Lechner, M. Ackermann, M. Kowalski, and W. Winter, “Astrophysical neutrinos flavored with Beyond the Standard Model physics,” Phys. Rev. D 96 (2017) no. 8, 083018, arXiv:1707.07684.
  * (31) V. Brdar and R. S. L. Hansen, “IceCube Flavor Ratios with Identified Astrophysical Sources: Towards Improving New Physics Testability,” JCAP 02 (2019) 023, arXiv:1812.05541.
  * (32) M. Bustamante and M. Ahlers, “Inferring the flavor of high-energy astrophysical neutrinos at their sources,” Phys. Rev. Lett. 122 (2019) no. 24, 241101, arXiv:1901.10087.
  * (33) A. Palladino, “The flavor composition of astrophysical neutrinos after 8 years of IceCube: an indication of neutron decay scenario?,” Eur. Phys. J. C79 (2019) no. 6, 500, arXiv:1902.08630.
  * (34) IceCube, J. Stachurska, “IceCube High Energy Starting Events at 7.5 Years – New Measurements of Flux and Flavor,” EPJ Web Conf. 207 (2019) 02005, arXiv:1905.04237.
  * (35) N. Song, S. W. Li, C. A. Argüelles, M. Bustamante, and A. C. Vincent, “The Future of High-Energy Astrophysical Neutrino Flavor Measurements,” JCAP 04 (2021) 054, arXiv:2012.12893.
  * (36) S. L. Glashow, “Resonant Scattering of Antineutrinos,” Phys. Rev. 118 (1960) 316–317.
  * (37) V. S. Berezinsky and A. Z. Gazizov, “Cosmic neutrino and the possibility of Searching for W bosons with masses 30-100 GeV in underwater experiments,” JETP Lett. 25 (1977) 254–256.
  * (38) R. W. Brown and F. W. Stecker, “Cosmic Ray Neutrino Tests for Heavier Weak Bosons and Cosmic Antimatter,” Phys. Rev. D 26 (1982) 373.
  * (39) L. A. Anchordoqui, H. Goldberg, F. Halzen, and T. J. Weiler, “Neutrinos as a diagnostic of high energy astrophysical processes,” Phys. Lett. B621 (2005) 18–21, arXiv:hep-ph/0410003.
  * (40) S. Hummer, M. Maltoni, W. Winter, and C. Yaguna, “Energy dependent neutrino flavor ratios from cosmic accelerators on the Hillas plot,” Astropart. Phys. 34 (2010) 205–224, arXiv:1007.0006.
  * (41) Z.-z. Xing and S. Zhou, “The Glashow resonance as a discriminator of UHE cosmic neutrinos originating from p-gamma and p-p collisions,” Phys. Rev. D84 (2011) 033006, arXiv:1105.4114.
  * (42) A. Bhattacharya, R. Gandhi, W. Rodejohann, and A. Watanabe, “The Glashow resonance at IceCube: signatures, event rates and $pp$ vs. $p\gamma$ interactions,” JCAP 1110 (2011) 017, arXiv:1108.3163.
  * (43) A. Bhattacharya, R. Gandhi, W. Rodejohann, and A. Watanabe, “On the interpretation of IceCube cascade events in terms of the Glashow resonance,” arXiv:1209.2422.
  * (44) V. Barger, J. Learned, and S. Pakvasa, “IceCube PeV Cascade Events Initiated by Electron-Antineutrinos at Glashow Resonance,” Phys. Rev. D87 (2013) no. 3, 037302, arXiv:1207.4571.
  * (45) V. Barger, L. Fu, J. G. Learned, D. Marfatia, S. Pakvasa, and T. J. Weiler, “Glashow resonance as a window into cosmic neutrino sources,” Phys. Rev. D90 (2014) 121301, arXiv:1407.3255.
  * (46) A. Palladino, G. Pagliaroli, F. L. Villante, and F. Vissani, “Double pulses and cascades above 2 PeV in IceCube,” Eur. Phys. J. C76 (2016) no. 2, 52, arXiv:1510.05921.
  * (47) I. M. Shoemaker and K. Murase, “Probing BSM Neutrino Physics with Flavor and Spectral Distortions: Prospects for Future High-Energy Neutrino Telescopes,” Phys. Rev. D93 (2016) no. 8, 085004, arXiv:1512.07228.
  * (48) L. A. Anchordoqui, M. M. Block, L. Durand, P. Ha, J. F. Soriano, and T. J. Weiler, “Evidence for a break in the spectrum of astrophysical neutrinos,” Phys. Rev. D95 (2017) no. 8, 083009, arXiv:1611.07905.
  * (49) M. D. Kistler and R. Laha, “Multi-PeV Signals from a New Astrophysical Neutrino Flux Beyond the Glashow Resonance,” Phys. Rev. Lett. 120 (2018) no. 24, 241105, arXiv:1605.08781.
  * (50) D. Biehl, A. Fedynitch, A. Palladino, T. J. Weiler, and W. Winter, “Astrophysical Neutrino Production Diagnostics with the Glashow Resonance,” JCAP 1701 (2017) 033, arXiv:1611.07983.
  * (51) S. Sahu and B. Zhang, “On the non-detection of Glashow resonance in IceCube,” JHEAp 18 (2018) 1–4, arXiv:1612.09043.
  * (52) G.-y. Huang and Q. Liu, “Hunting the Glashow Resonance with PeV Neutrino Telescopes,” JCAP 03 (2020) 005, arXiv:1912.02976.
  * (53) S. Zhou, “Cosmic Flavor Hexagon for Ultrahigh-energy Neutrinos and Antineutrinos at Neutrino Telescopes,” arXiv:2006.06181.
  * (54) IceCube, M. G. Aartsen et al., “The IceCube Neutrino Observatory - Contributions to ICRC 2015 Part II: Atmospheric and Astrophysical Diffuse Neutrino Searches of All Flavors,” in 34th International Cosmic Ray Conference.  10, 2015.  arXiv:1510.05223.
  * (55) R. Gauld, “Precise predictions for multi-TeV and PeV energy neutrino scattering rates,” Phys. Rev. D100 (2019) no. 9, 091301, arXiv:1905.03792.
  * (56) A. Garcia, R. Gauld, A. Heijboer, and J. Rojo, “Complete predictions for high-energy neutrino propagation in matter,” JCAP 09 (2020) 025, arXiv:2004.04756.
  * (57) I. Alikhanov, “On parton distributions in a photon gas,” Eur. Phys. J. C 65 (2010) 269–273, arXiv:0812.0937.
  * (58) A. Loewy, S. Nussinov, and S. L. Glashow, “The Effect of Doppler Broadening on the $6.3\ PeV$ $W^{-}$ Resonance in $\bar{\nu}_{e}e^{-}$ Collisions,” arXiv:1407.4415.
  * (59) ALEPH, DELPHI, L3, OPAL, SLD, LEP Electroweak Working Group, SLD Electroweak Group, SLD Heavy Flavour Group, S. Schael et al., “Precision electroweak measurements on the $Z$ resonance,” Phys. Rept. 427 (2006) 257–454, arXiv:hep-ex/0509008.
  * (60) M. Cacciari, A. Deandrea, G. Montagna, and O. Nicrosini, “QED Structure Functions: A Systematic Approach,” Europhysics Letters (EPL) 17 (1992) no. 2, 123–128.
  * (61) IceCube, M. G. Aartsen et al., “Observation of High-Energy Astrophysical Neutrinos in Three Years of IceCube Data,” Phys. Rev. Lett. 113 (2014) 101101, arXiv:1405.5303.
  * (62) A. M. Hillas, “The Origin of Ultra-High-Energy Cosmic Rays,” Annual Review of Astronomy and Astrophysics 22 (1984) no. 1, 425–444.
  * (63) G. Cowan, Statistical data analysis.  Clarendon Press, 1998.
  * (64) P. Baerwald, S. Hümmer, and W. Winter, “Magnetic field and flavor effects on the gamma-ray burst neutrino flux,” Phys. Rev. D 83 (2011) no. 6, , arXiv:1009.4010.
  * (65) S. Hümmer, M. Rüger, F. Spanier, and W. Winter, “Simplified models for photohadronic interactions in cosmic accelerators,” Astrophys. J. 721 (2010) no. 1, 630–652, arXiv:1002.1310.
  * (66) IceCube, M. G. Aartsen et al., “IceCube-Gen2: A Vision for the Future of Neutrino Astronomy in Antarctica,” arXiv:1412.5106.
  * (67) IceCube-Gen2, M. G. Aartsen et al., “IceCube-Gen2: the window to the extreme Universe,” J. Phys. G 48 (2021) no. 6, 060501, arXiv:2008.04323.
  * (68) Baikal-GVD, A. D. Avrorin et al., “Baikal-GVD: status and prospects,” EPJ Web Conf. 191 (2018) 01006, arXiv:1808.10353.
  * (69) KM3Net, S. Adrian-Martinez et al., “Letter of intent for KM3NeT 2.0,” J. Phys. G 43 (2016) no. 8, 084001, arXiv:1601.07459.
  * (70) P-ONE, M. Agostini et al., “The Pacific Ocean Neutrino Experiment,” Nature Astron. 4 (2020) no. 10, 913–915, arXiv:2005.09493.
  * (71) A. Romero-Wolf et al., “An Andean Deep-Valley Detector for High-Energy Tau Neutrinos,” in Latin American Strategy Forum for Research Infrastructure.  2, 2020.  arXiv:2002.06475.
  * (72) Z. P. Ye et al., “Proposal for a neutrino telescope in South China Sea,” arXiv:2207.04519.
  * (73) G.-y. Huang, S. Jana, M. Lindner, and W. Rodejohann, “Probing new physics at future tau neutrino telescopes,” JCAP 02 (2022) no. 02, 038, arXiv:2112.09476.
  * (74) A. Coleman et al., “Ultra high energy cosmic rays The intersection of the Cosmic and Energy Frontiers,” Astropart. Phys. 147 (2023) 102794, arXiv:2205.05845.
  * (75) M. Ackermann et al., “High-energy and ultra-high-energy neutrinos: A Snowmass white paper,” JHEAp 36 (2022) 55–110, arXiv:2203.08096.
  * (76) V. B. Valera, M. Bustamante, and C. Glaser, “Near-future discovery of the diffuse flux of ultra-high-energy cosmic neutrinos,” arXiv:2210.03756.
  * (77) LHAASO*†, LHAASO, Z. Cao et al., “Peta–electron volt gamma-ray emission from the Crab Nebula,” Science 373 (2021) no. 6553, 425–430, arXiv:2111.06545.
  * (78) Z. Cao et al., “Ultrahigh-energy photons up to 1.4 petaelectronvolts from 12 $\gamma$-ray Galactic sources,” Nature 594 (2021) no. 7861, 33–36.
  * (79) T. Sudoh and J. F. Beacom, “Where are Milky Way’s hadronic PeVatrons?,” Phys. Rev. D 107 (2023) no. 4, 043002, arXiv:2209.03970.
  * (80) M. Bustamante, “New limits on neutrino decay from the Glashow resonance of high-energy cosmic neutrinos,” arXiv:2004.06844.
  * (81) T. Ježo, M. Klasen, F. Lyonnet, F. Montanet, I. Schienbein, and M. Tartare, “Can new heavy gauge bosons be observed in ultra-high energy cosmic neutrino events?,” Phys. Rev. D89 (2014) no. 7, 077702, arXiv:1401.6012.
  * (82) K. S. Babu, P. S. Dev, S. Jana, and Y. Sui, “Zee-Burst: A New Probe of Neutrino Non-Standard Interactions at IceCube,” arXiv:1908.02779.
  * (83) U. K. Dey, N. Nath, and S. Sadhukhan, “Charged Higgs effects in IceCube: PeV events and NSIs,” JHEP 09 (2021) 113, arXiv:2010.05797.
  * (84) K. S. Babu, P. S. B. Dev, and S. Jana, “Probing neutrino mass models through resonances at neutrino telescopes,” Int. J. Mod. Phys. A 37 (2022) no. 11n12, 2230003, arXiv:2202.06975.
  * (85) D.-H. Xu and S.-J. Rong, “Connect the Lorentz Violation to the Glashow Resonance Event,” arXiv:2211.05478.
  * (86) C. A. Argüelles et al., “Snowmass White Paper: Beyond the Standard Model effects on Neutrino Flavor,” in 2022 Snowmass Summer Study.  3, 2022.  arXiv:2203.10811.
  * (87) G.-y. Huang, S. Jana, M. Lindner, and W. Rodejohann, “Probing Heavy Sterile Neutrinos at Ultrahigh Energy Neutrino Telescopes via the Dipole Portal,” arXiv:2204.10347.
  * (88) G.-y. Huang, “Double and multiple bangs at tau neutrino telescopes,” Eur. Phys. J. C 82 (2022) no. 12, 1089, arXiv:2207.02222.
  * (89) R. Heighton, L. Heurtier, and M. Spannowsky, “Hunting for Neutral Leptons with Ultra-High-Energy Cosmic Rays,” arXiv:2303.11352.
  * (90) E. Clementi and D. L. Raimondi, “Atomic Screening Constants from SCF Functions,” jcp 38 (1963) no. 11, 2686–2689.