Cosmic tau neutrino detection via Cherenkov signals from air showers from Earth-emerging taus

TL;DR

This paper models Earth-emerging tau neutrinos and their air shower signals to evaluate the sensitivity of space-based detectors like POEMMA for high-energy neutrino detection above 10 PeV.

Contribution

It provides a detailed flux and spectrum calculation of Earth-emerging tau leptons using a layered Earth model and assesses detector sensitivity considering systematic uncertainties.

Findings

Quantifies systematic uncertainties from cross sections and energy loss models.

Estimates neutrino sensitivity for space-based detectors above 10 PeV.

Simulates Cherenkov signals from tau-induced air showers.

Abstract

We perform a new, detailed calculation of the flux and energy spectrum of Earth-emerging $\tau$-leptons generated from the interactions of tau neutrinos and antineutrinos in the Earth. A layered model of the Earth is used to describe the variable density profile of the Earth. Different assumptions regarding the neutrino charged- and neutral-current cross sections as well as the $\tau$-lepton energy loss models are used to quantify the systematic uncertainty from these on the results. A baseline simulation is then used to generate the optical Cherenkov signal from upward-moving extensive air showers generated by the $\tau$-lepton decay in the atmosphere, applicable to a range of space-based instruments. We use this simulation to determine the neutrino sensitivity for $E_\nu>$ 10 PeV for a space-based experiment with performance similar to that for the Probe of Extreme MultiMessenger Astrophysics (POEMMA) mission currently under study.

Figures

Tables

Table 1. Table 1: Density ρ 𝜌 \rho as a function of vertical depth d = R E − r 𝑑 subscript 𝑅 𝐸 𝑟 d=R_{E}-r below the Earth’s surface for β tr ≤ 50 ∘ subscript 𝛽 tr superscript 50 \beta_{\rm tr}\leq 50^{\circ} based on the PREM parameterization Dziewonski and Anderson ( 1981 ) of the Earth’s density.

$\rho \left[g/{\mathrm{cm}}^3\right]$	Depth $d=R_E-r$
1.02	$d\le 3.0\mathrm{km}$
2.6
2.9
3.4
3.8
4.4
4.8

Table 2. Table 2: For a given Δ ​ α = α H − α Δ 𝛼 subscript 𝛼 𝐻 𝛼 \Delta\alpha=\alpha_{H}-\alpha as measured from altitude h = 33 , 525 ℎ 33 525 h=33,\ 525 and 1000 km, the viewing angle relative to the horizon at Earth β v subscript 𝛽 𝑣 \beta_{v} , all in degrees.

$\Delta \alpha$	${\beta }_v\left(33\mathrm{km}\right)$	${\beta }_v\left(525\mathrm{km}\right)$	${\beta }_v\left(1000\mathrm{km}\right)$
1	3.6	7.0	8.2
2	5.2	10.0	11.7
3	6.6	12.3	14.5
4	7.9	14.4	16.9
5	9.1	16.2	19.0
6	10.3	18.0	21.0
7	11.4	19.6	22.8
8	12.6	21.2	24.6
9	13.6	22.6	26.3
10	14.7	24.1	27.9
15	20.0	30.8	35.4
20	25.2	37.0	42.2

Table 3. Table 3: Comparison of the cap ⟨ A cap ​ Ω ⟩ geo subscript delimited-⟨⟩ subscript 𝐴 cap Ω geo \langle A_{\rm cap}\Omega\rangle_{\rm geo} and zone ⟨ A zone ​ Ω ⟩ geo subscript delimited-⟨⟩ subscript 𝐴 zone Ω geo \langle A_{\rm zone}\Omega\rangle_{\rm geo} geometric apertures for several altitudes h ℎ h when θ d = 1.5 ∘ subscript 𝜃 𝑑 superscript 1.5 \theta_{d}=1.5^{\circ} and Δ ​ α = 7 ∘ Δ 𝛼 superscript 7 \Delta\alpha=7^{\circ} . The final column is the ratio ⟨ A zone ​ Ω ⟩ geo / ⟨ A cap ​ Ω ⟩ geo subscript delimited-⟨⟩ subscript 𝐴 zone Ω geo subscript delimited-⟨⟩ subscript 𝐴 cap Ω geo \langle A_{\rm zone}\Omega\rangle_{\rm geo}/\langle A_{\rm cap}\Omega\rangle_{\rm geo} .

Altitude $h$ [km]	Cap [km2 sr]	Zone [km2 sr]	Zone/Cap
3	5.2	4.5	0.87
4	7.9	6.8	0.85
33	178	124	0.70
525	8,480	4,072	0.48
1000	18,857	8,538	0.45

Table 4. Table 4: The tau exit probability for E ν τ = 10 7 subscript 𝐸 subscript 𝜈 𝜏 superscript 10 7 E_{\nu_{\tau}}=10^{7} , 10 8 superscript 10 8 10^{8} , 1 0 9 superscript 0 9 0^{9} and 10 10 superscript 10 10 10^{10} GeV assuming the ALLM structure function extrapolation for the photonuclear energy loss parameter, as a function of β tr subscript 𝛽 tr \beta_{\rm tr} . The standard model cross section for neutrino-nucleon interactions is assumed.

${\beta }_{\mathrm{tr}}$ [∘]	${10}^7$ GeV	${10}^8$ GeV	${10}^9$ GeV	${10}^{10}$ GeV
1	2.89e-05	7.41e-04	6.85e-03	2.54e-02
3	2.48e-05	5.55e-04	3.27e-03	5.15e-03
5	2.07e-05	3.92e-04	1.49e-03	1.23e-03
7	2.01e-05	2.50e-04	5.79e-04	3.36e-04
10	1.29e-05	1.48e-04	2.20e-04	1.38e-04
12	1.49e-05	1.06e-04	1.35e-04	8.20e-05
15	9.30e-06	6.42e-05	6.81e-05	5.48e-05
17	8.90e-06	4.99e-05	5.40e-05	3.77e-05
20	7.90e-06	3.63e-05	3.31e-05	2.59e-05
25	4.70e-06	1.57e-05	1.43e-05	1.31e-05
30	2.44e-06	6.73e-06	5.92e-06	5.85e-06
35	1.44e-06	3.13e-06	2.84e-06	2.67e-06
40	9.00e-07	1.95e-06	1.83e-06	1.59e-06

Table 5. Table 5: The cumulative distribution function displayed in Fig. 26 for E ν τ = 10 7 subscript 𝐸 subscript 𝜈 𝜏 superscript 10 7 E_{\nu_{\tau}}=10^{7} GeV and β tr = 1 ∘ subscript 𝛽 tr superscript 1 \beta_{\rm tr}=1^{\circ} , 5 ∘ superscript 5 5^{\circ} , 10 ∘ superscript 10 10^{\circ} and 20 ∘ superscript 20 20^{\circ} , as a function of z τ = E τ / E ν τ subscript 𝑧 𝜏 subscript 𝐸 𝜏 subscript 𝐸 subscript 𝜈 𝜏 z_{\tau}=E_{\tau}/E_{\nu_{\tau}} .

$z_{\tau }$	$1^{\circ }$	$5^{\circ }$	${10}^{\circ }$	${20}^{\circ }$
1.41e-01	3.45e-03	0.00e+00	0.00e+00	2.53e-02
1.78e-01	3.45e-03	4.83e-03	4.65e-02	5.06e-02
2.24e-01	6.90e-03	9.66e-03	6.98e-02	6.33e-02
2.82e-01	1.03e-02	1.93e-02	1.09e-01	8.86e-02
3.55e-01	2.41e-02	3.38e-02	1.32e-01	1.77e-01
4.47e-01	4.13e-02	5.31e-02	1.63e-01	2.78e-01
5.62e-01	1.41e-01	1.59e-01	2.48e-01	3.92e-01
7.08e-01	2.96e-01	3.62e-01	3.80e-01	5.31e-01
8.91e-01	9.96e-01	1.00e+00	1.00e+00	9.99e-01

Table 6. Table 6: The cumulative distribution function displayed in Fig. 27 for E ν τ = 10 8 subscript 𝐸 subscript 𝜈 𝜏 superscript 10 8 E_{\nu_{\tau}}=10^{8} GeV and β tr = 1 ∘ subscript 𝛽 tr superscript 1 \beta_{\rm tr}=1^{\circ} , 5 ∘ superscript 5 5^{\circ} , 10 ∘ superscript 10 10^{\circ} and 20 ∘ superscript 20 20^{\circ} , as a function of z τ = E τ / E ν τ subscript 𝑧 𝜏 subscript 𝐸 𝜏 subscript 𝐸 subscript 𝜈 𝜏 z_{\tau}=E_{\tau}/E_{\nu_{\tau}} .

$z_{\tau }$	$1^{\circ }$	$5^{\circ }$	${10}^{\circ }$	${20}^{\circ }$
1.41e-02	0.00e+00	0.00e+00	1.35e-03	0.00e+00
1.78e-02	1.35e-04	0.00e+00	2.70e-03	2.75e-03
2.24e-02	1.35e-04	5.10e-04	4.05e-03	8.26e-03
2.82e-02	1.35e-04	7.65e-04	1.08e-02	2.21e-02
3.55e-02	1.35e-04	3.32e-03	1.55e-02	3.86e-02
4.47e-02	2.70e-04	5.62e-03	2.23e-02	6.61e-02
5.62e-02	8.10e-04	9.70e-03	3.18e-02	1.07e-01
7.08e-02	2.16e-03	1.53e-02	4.26e-02	1.40e-01
8.91e-02	3.64e-03	2.35e-02	6.42e-02	1.74e-01
1.12e-01	6.07e-03	3.42e-02	9.53e-02	2.42e-01
1.41e-01	9.58e-03	5.28e-02	1.25e-01	3.06e-01
1.78e-01	1.55e-02	8.21e-02	1.69e-01	3.77e-01
2.24e-01	3.01e-02	1.11e-01	2.26e-01	4.71e-01
2.82e-01	4.98e-02	1.53e-01	2.97e-01	5.48e-01
3.55e-01	9.11e-02	2.18e-01	3.71e-01	6.34e-01
4.47e-01	1.71e-01	3.14e-01	4.71e-01	7.08e-01
5.62e-01	3.22e-01	4.57e-01	5.92e-01	8.10e-01
7.08e-01	5.86e-01	6.77e-01	7.78e-01	8.87e-01
8.91e-01	1.00e+00	9.99e-01	9.98e-01	9.97e-01

Table 7. Table 7: The cumulative distribution function displayed in Fig. 28 for E ν τ = 10 9 subscript 𝐸 subscript 𝜈 𝜏 superscript 10 9 E_{\nu_{\tau}}=10^{9} GeV and β tr = 1 ∘ subscript 𝛽 tr superscript 1 \beta_{\rm tr}=1^{\circ} , 5 ∘ superscript 5 5^{\circ} , 10 ∘ superscript 10 10^{\circ} and 20 ∘ superscript 20 20^{\circ} , as a function of z τ = E τ / E ν τ subscript 𝑧 𝜏 subscript 𝐸 𝜏 subscript 𝐸 subscript 𝜈 𝜏 z_{\tau}=E_{\tau}/E_{\nu_{\tau}} .

$z_{\tau }$	$1^{\circ }$	$5^{\circ }$	${10}^{\circ }$	${20}^{\circ }$
1.41e-03	0.00e+00	3.36e-05	5.37e-04	2.60e-03
2.24e-03	8.74e-06	3.83e-04	2.60e-03	1.46e-02
3.55e-03	4.23e-05	1.26e-03	8.33e-03	4.06e-02
5.62e-03	1.24e-04	3.43e-03	2.17e-02	9.14e-02
8.91e-03	3.05e-04	8.29e-03	4.87e-02	1.85e-01
1.41e-02	7.71e-04	1.82e-02	9.69e-02	3.39e-01
2.24e-02	2.24e-03	3.73e-02	1.73e-01	5.06e-01
3.55e-02	6.95e-03	6.89e-02	2.79e-01	6.64e-01
5.62e-02	2.28e-02	1.22e-01	4.00e-01	8.05e-01
8.91e-02	6.79e-02	2.05e-01	5.27e-01	8.94e-01
1.41e-01	1.66e-01	3.28e-01	6.54e-01	9.47e-01
2.24e-01	3.26e-01	4.86e-01	7.73e-01	9.76e-01
3.55e-01	5.36e-01	6.68e-01	8.74e-01	9.88e-01
5.62e-01	7.71e-01	8.53e-01	9.52e-01	9.97e-01
8.91e-01	1.00e+00	1.00e+00	1.00e+00	9.98e-01

Table 8. Table 8: The cumulative distribution function displayed in Fig. 29 for E ν τ = 10 10 subscript 𝐸 subscript 𝜈 𝜏 superscript 10 10 E_{\nu_{\tau}}=10^{10} GeV and β tr = 1 ∘ subscript 𝛽 tr superscript 1 \beta_{\rm tr}=1^{\circ} , 5 ∘ superscript 5 5^{\circ} , 10 ∘ superscript 10 10^{\circ} and 20 ∘ superscript 20 20^{\circ} , as a function of z τ = E τ / E ν τ subscript 𝑧 𝜏 subscript 𝐸 𝜏 subscript 𝐸 subscript 𝜈 𝜏 z_{\tau}=E_{\tau}/E_{\nu_{\tau}} .

$z_{\tau }$	$1^{\circ }$	$5^{\circ }$	${10}^{\circ }$	${20}^{\circ }$
1.41e-04	7.86e-07	2.45e-05	5.88e-04	3.09e-03
2.24e-04	6.69e-06	3.51e-04	4.04e-03	1.82e-02
3.55e-04	1.77e-05	1.60e-03	1.24e-02	4.53e-02
5.62e-04	5.46e-05	4.41e-03	3.13e-02	1.15e-01
8.91e-04	1.36e-04	1.13e-02	7.64e-02	2.24e-01
1.41e-03	3.43e-04	2.47e-02	1.51e-01	3.77e-01
2.24e-03	9.75e-04	5.02e-02	2.67e-01	5.68e-01
3.55e-03	2.90e-03	9.29e-02	4.14e-01	7.40e-01
5.62e-03	9.70e-03	1.57e-01	5.74e-01	8.77e-01
8.91e-03	2.99e-02	2.42e-01	7.28e-01	9.50e-01
1.41e-02	7.46e-02	3.41e-01	8.38e-01	9.83e-01
2.24e-02	1.49e-01	4.45e-01	9.07e-01	9.92e-01
3.55e-02	2.50e-01	5.49e-01	9.47e-01	9.96e-01
5.62e-02	3.67e-01	6.47e-01	9.70e-01	9.97e-01
8.91e-02	4.89e-01	7.34e-01	9.81e-01	9.97e-01
1.41e-01	6.10e-01	8.12e-01	9.90e-01	9.97e-01
2.24e-01	7.23e-01	8.77e-01	9.95e-01	9.97e-01
3.55e-01	8.27e-01	9.33e-01	9.98e-01	9.97e-01
5.62e-01	9.20e-01	9.75e-01	9.99e-01	9.97e-01
8.91e-01	1.00e+00	1.00e+00	1.00e+00	9.97e-01

Equations

$$c=2R_E\sin {\beta }_{\mathrm{tr}},$$

$$X\left({\beta }_{\mathrm{tr}}\right)={\int }_0^c\rho \left(s^{\prime },{\beta }_{\mathrm{tr}}\right)d{s}^{\prime }$$

$$\Biggl{\langle}\frac{dE_{\tau}}{dX}\Biggr{\rangle}=-(a_{\tau}+b_{\tau}E_{\tau})\ .$$

$$b_{\tau }=\sum _i{b}_{\tau }^i=b_{\tau }^{\mathrm{brem}}+b_{\tau }^{\mathrm{pair}}+b_{\tau }^{\mathrm{nuc}},$$

$$b_{\tau }^i(E)=\frac{N}{A}\int dyy\frac{d{\sigma }_{\tau A}^i(y,E)}{dy}$$

$$b_{\tau }^i(E)$$

$$\frac{d{n}_{\tau \to {\nu }_{\tau }}}{dy}=g_0(y)-g_1(y)$$

$$g_0(y)$$

$$g_1(y)$$

$$P_{\mathrm{surv}}$$

$$P_{\mathrm{decay}}$$

$$s(a,{\beta }_{tr})$$

$${\rho }_{\gamma }(a,{\beta }_{\mathrm{tr}},E_{\mathrm{shr}})$$

$$⟨A\Omega ⟩\left(E_{{\nu }_{\tau }}\right)={\int }_S{\int }_{\Delta {\Omega }_{\mathrm{tr}}}{P}_{\mathrm{obs}}\hat{r}\cdot \hat{n}dSd{\Omega }_{\mathrm{tr}},$$

$$G$$

$$⟨A\Omega {⟩}_{\mathrm{geo}}$$

$$P_{\mathrm{obs}}$$

$$\times \left[\int d{s}^{\prime }{p}_{\mathrm{decay}}(s^{\prime })p_{\det }\left(E_{\tau },{\theta }_E,{\beta }_{\mathrm{tr}},s^{\prime }\right)\right]d{E}_{\tau },$$

$$p_{\det }=H\left[{\theta }_{\mathrm{Ch}}-\theta \right]H\left[s_{\mathrm{win}}-s\right]H\left[N_{\mathrm{PE}}-N_{\mathrm{PE}}^{\min }\right],$$

$$H\left(x\right)=\left\{\begin{array}[]{ll}0&\mbox{if $x<0$};\\ 1&\mbox{if $x\geq 0$}.\end{array}\right.$$

$$N_{\mathrm{PE}}={\rho }_{\gamma }(a,{\beta }_{\mathrm{tr}},E_{\mathrm{shr}})\times 2.5{\mathrm{m}}^2\times 0.2.$$

$$F_{\mathrm{sens}}\left(E_{{\nu }_{\tau }}\right)=\frac{2.44\times N_{\nu }}{\ln (10)\times E_{{\nu }_{\tau }}\times ⟨A\Omega ⟩\left(E_{{\nu }_{\tau }}\right)\times t_{\mathrm{obs}}},$$

$$F_{\nu}(E_{\nu})=F_{0}\times\Biggl{(}\frac{\rm GeV}{E_{\nu}}\Biggr{)}^{\gamma}$$

$$N_{\mathrm{evts}}^{{\nu }_{\tau }}={\int }_{10^{-0.5}{E}_{\nu }}^{10^{0.5}{E}_{\nu }}dE\frac{F_{\nu }(E)}{N_{\nu }}⟨A\Omega ⟩(E)t_{\mathrm{obs}}=2.44.$$

$$N_{\mathrm{evts}}\left(>E_{\tau }^{\min }\right)$$

$$\sin {\alpha }_H$$

$${\theta }_E^H$$

$$\frac{\cos {\beta }_v}{R_E+h}$$

$$v\sin \alpha$$

$$v^2=(R_E+h{)}^2+R_E^2-2R_E(R_E+h)\cos {\theta }_E,$$

$$\cos {\theta }_v=\frac{\left(R_E+h\right)\cos {\theta }_E-R_E}{\sqrt{{\left(R_E+h\right)}^2+R_E^2-2\left(R_E+h\right)R_E\cos {\theta }_E}}.$$

$$\alpha +{\beta }_v+{\theta }_E=90^{\circ }.$$

$$G$$

Full text

I Introduction

The measurement of the spectrum of the very-high energy (VHE: E_{\nu}\mathrel{\hbox{\hbox to0.0pt{\lower 3.22916pt\hbox{\sim}\hss} \kern-3.00003pt\raise 1.72218pt\hbox{>}}}1 PeV) neutrino and antineutrino (hereafter denoted collectively as neutrinos) component of the cosmic radiation and its angular distribution on the sky provides a unique probe of high-energy astrophysical phenomena. A by-product of cosmic ray acceleration, astrophysical neutrinos can reveal the environments of sources of ultra-high energy cosmic rays (UHECR) Gaisser et al. (1995); Ahlers and Halzen (2017). The results from IceCube (Aartsen et al., 2013, 2014, 2015) demonstrate the existence of an extra-solar system astrophysical neutrino flux with energies from above 10 TeV to potentially as high as 10 PeV. A neutrino event detected by IceCube that appears correlated at $3\sigma$ with gamma-flaring data from an active galactic nuclei source has been reported Aartsen et al. (2018a). Gamma ray bursts Murase (2007), newborn pulsars Fang et al. (2014), active galactic nuclei Murase (2017), galactic clusters with central sources Murase et al. (2008); Fang and Murase (2018) and UHECR photodisintegration within cosmic ray sources Unger et al. (2015) are among candidates sources for the diffuse astrophysical flux of neutrinos. Astrophysical neutrinos are key to the multi-messenger approach to understand sources of cosmic radiation Spurio (2018).

At the highest energies, neutrinos are anticipated from UHECR that attenuate through interactions with photons in transit from sources Berezinsky and Zatsepin (1969); Hill and Schramm (1985); Engel et al. (2001); Anchordoqui et al. (2007); Kotera et al. (2010); Decerprit and Allard (2011); Roulet et al. (2013) with the Greisen-Zatsepin-Kuzmin (GZK) cutoff Greisen (1966); Zatsepin and Kuzmin (1966) being a signal of this process. The details of the cosmogenic neutrino spectrum would provide invaluable information to the cosmic ray acceleration process, source distribution, source evolution, and the UHECR nuclear composition (Kotera et al., 2010). While the existence of the cosmogenic neutrinos is implied by the baryonic component in cosmic rays, the detection of these neutrinos has remained elusive nul . It is one of the most important measurements in astroparticle physics.

The weak neutrino interaction cross section makes neutrinos a critical and unique component of multi-messenger astronomy and astrophysics. The neutrino horizon extends far beyond that for UHECR. Measurements of cosmogenic neutrinos have the potential to provide information about sources much farther than those responsible for the observed flux of UHECR. Measurements of astrophysical neutrinos probe the environments of cosmic accelerators.

While the weak interactions of neutrinos give a benefit in their transit to the Earth, their detection requires very large target volumes. Direct detection of neutrino interactions in, e.g., the IceCube instrumented volume, are also augmented by muons from $\nu_{\mu}$ charged current interactions that produce a muon outside the detector, thereby increasing the effective detection volume. Tau neutrino production of $\tau$-leptons can also increase the effective volume at higher energies due to the Lorentz-boosted lifetime. Over cosmological distances, the $1:2$ ratio of $\nu_{e}+\bar{\nu}_{e}:\nu_{\mu}+\bar{\nu}_{\mu}$ produced in sources of UHECR or in their transit to the Earth yields a nearly equal flux of electron, muon and tau neutrinos Learned and Pakvasa (1995). This leads to characteristic upward-going tau neutrino induced signals in underground detectors and ground-based detection of air showers Halzen and Saltzberg (1998); Bertou et al. (2002); Feng et al. (2002); Lachaud et al. (2002); Hou and Huang (2002); Tseng et al. (2003); Aramo et al. (2005); Dutta et al. (2005); Asaoka and Sasaki (2013); Fargion (2002); Bottai and Giurgola (2003); Fargion et al. (2004); Fargion (2003); Palomares-Ruiz et al. (2006); Abreu et al. (2012); Zas (2018); Fang et al. (2018); Alvarez-Muñiz et al. (2018a); Ahnen et al. (2018); Otte (2019); Aab:2019auo . The potential to use sub-orbital and space-based measurements of extensive air showers (EAS) induced from neutrino interactions either in the Earth (Domokos and Kovesi-Domokos, 1997, 1998) or Earth’s atmosphere has been recognized as a way to achieve even larger neutrino target masses, greater than $10^{13}$ metric tons for the atmosphere (Stecker et al., 2004), for example. In particular, the signals from upward-moving EAS that come from tau neutrino interactions within the Earth Fargion (2002); Bottai and Giurgola (2003); Fargion et al. (2004); Fargion (2003); Palomares-Ruiz et al. (2006) provide a path to measure the astrophysical and the cosmogenic neutrino fluxes above $\sim 1$ PeV, with a huge neutrino target mass (Krizmanic and Mitchell, 2011; CHANT, ).

Indirect detection techniques for measuring the characteristics of upward-going EASs from VHE neutrinos include: 1) the detection of the beamed optical Cherenkov radiation from EAS particles, and 2) the detection of the coherent radio radiation from the EAS. The latter comes from electric fields induced by separation of positive and negative EAS particles in the Earth’s magnetic field (geomagnetic radiation) and from a time-varying net charge in showers initiated by VHE neutrinos in dense media such as ice (the Askaryan effect Askar’yan (1962)). Measurements of either type of signal, optical Cherenkov or radio, can be leveraged to determine the energy of the neutrino primary while having excellent angular resolution in the incident neutrino direction. The long-duration balloon flights of ANITA have demonstrated the capability of using a sub-orbital instrument for detecting radio signals from downward-going and horizontally propagating UHECRs Hoover et al. (2010); Schoorlemmer et al. (2016); Gorham et al. (2016). Searches for radio signals from upward-going, neutrino-induced cascades in the Antarctic ice have led to limits on the diffuse flux of UHE neutrinos at energies above $10^{10}$ GeV Gorham et al. (2018, 2019), culminating in stringent constraints Barwick et al. (2006) on “top-down” models of UHECR and UHE neutrino production via phenomena associated with relics from the Big Bang, from phase transitions or from physics near the Grand Unified scale in the early universe (e.g. Refs. Yoshida et al. (1997); Berezinsky and Vilenkin (2000); Sigl et al. (1998)).

Recently, the ANITA Collaboration reported the detection of two anomalous events with radio signal characteristics that appear to be compatible with upward-moving EASs. The inferred EAS energies and the projected path lengths through the Earth imply that these events were likely not initiated by $\nu_{\tau}$ interactions in the Earth Gorham et al. (2016, 2018); Romero-Wolf et al. (2018a). Nevertheless, these events have fueled speculation that they may signify physics beyond the Standard Model, such as Earth-interacting sterile neutrinos Cherry and Shoemaker (2018); Huang (2018), the interactions or production of supersymmetric particles Collins et al. (2019); Connolly et al. (2018); Fox et al. (2018), neutrino-induced supersymmetric sphaleron transitions Anchordoqui and Antoniadis (2019), or the decay of superheavy dark matter Anchordoqui et al. (2018).

Several sub-orbital (EUSO-SPB2, Adams et al. (2017)), space-based (CHANT, CHANT ; POEMMA, Olinto et al. (2018)), and ground-based (Trinity, Otte (2019)) instruments have been proposed to search for EASs from Earth-skimming tau neutrinos via optical Cherenkov radiation, which would allow for sensitivity to neutrino energies above $\sim 10$ PeV, even though the duty cycle for optical Cherenkov is $\sim 20$% (compared to $\sim 100$% for radio). POEMMA, the Probe Of Extreme Multi-Messenger Astrophysics, is a NASA probe-class mission under study that will detect UHECRs and VHE neutrinos from space Olinto et al. (2018). In this paper, we report on the development and the results from an initial, end-to-end calculation of the optical Cherenkov signal generated from upward-moving EAS from the decay of $\tau$-leptons sourced by tau neutrino interactions in the Earth. Many of the results presented here are not specific to the POEMMA concept. To perform a baseline simulation to calculate the neutrino sensitivity for a space-based experiment, we use the 525 km orbiting altitude, a 2.5 m2 detector collecting area, a photo-detection efficiency of 20%, and other detector characteristics of POEMMA.

In the next section, we give a brief overview of satellite and balloon-based neutrino detection starting with the geometry of the viewed surface of the Earth. Appendix A has many of the detailed relations between the various angles that describe the instrument viewing and the emerging showers. The treatment of tau neutrino interactions to produce $\tau$-leptons and the subsequent tau energy loss and/or decay in the Earth to produce a lower energy neutrino, are discussed in Sec. II. We show the flux independent and flux dependent $\nu_{\tau}\to\tau$ transmission results.

In Sec. III, we review $\tau$-lepton decay probabilities in the atmosphere and the geometry of decays as a function of altitude. We also describe how air showers in the atmosphere are modeled to get the photon density to arrive at a detector like POEMMA at 525 km altitude.

Sec. IV includes results that incorporate detection by a POEMMA-like optical Cherenkov detector as an application of our transmission results. We evaluate the aperture and flux independent sensitivity from our simulations. We also estimate the number of events that could be detected for two flux models for the isotropic cosmogenic flux Kotera et al. (2010).

We summarize our conclusions and discuss uncertainties and potential improvements to the simulation in Sec. V. In Appendix A, we collect further details of the geometry of the Earth, as seen at altitude. We display tables for the cumulative probability functions for $\nu_{\tau}\to\tau$ as a function of outgoing $\tau$-lepton energy fraction of the incident neutrino energy, for selected incident neutrino energies and angles, in Appendix B. Finally, in Appendix C we outline in more detail inputs to the Monte Carlo simulation of the detection probabilities as a function of energy and angle.

II Neutrino and tau propagation in Earth

The Earth-skimming technique for detecting tau neutrinos relies on using the Earth as a neutrino converter. The feasibility of this technique depends on the details of tau neutrino propagation along trajectories through the Earth that are determined by the detection geometry, shown in Fig. 1. We begin with a summary of our notation in Sec. II.1. Along the chord $c$ through the Earth, a neutrino interacts to produce a tau, then the tau loses energy and may decay with a neutrino in the final state. This regeneration process can occur multiple times, depending on the energy and angle. In Secs. II.2, II.3 and II.4, we describe how we incorporate these elements in our simulation. Sec. II.5 shows results for tau exit probabilities, and tau energy distributions, for a fixed tau neutrino incident energy and angle. Transmission functions, the ratios of outgoing taus to incoming tau neutrinos, as functions energy and angle are shown in Sec. II.6. The results presented in Sec. II are related to tau neutrino and tau propagation in the Earth. Sec. III is focused on tau neutrino induced air showers from tau decay and their detection.

II.1 Geometry of Neutrino Trajectories and Column Depth Through the Earth

For neutrino energies above a PeV, the neutrino and antineutrino interaction cross sections are essentially equal, and their interaction length is shorter than the column depth along the diameter of the Earth Gandhi et al. (1996, 1998). Consequently, neutrino detectors at altitude $h$ above the Earth need to point near the limb. The limb is at a viewing angle $\alpha_{H}$ away from the nadir in our notation. A given point on the Earth’s surface in view of the detector is described by the angle $\theta_{E}$, as shown in Fig. 1. It corresponds to a detector viewing angle $\alpha$ away from the satellite nadir. The line of sight from the point to the detector is at an elevation angle $\beta_{v}$ and angle $\theta_{v}$ from the local zenith. The distance $v$ is the distance along the line of sight. Relations between angles and distances are listed in Appendix A.

In order for the detector to receive Cherenkov radiation from a neutrino-induced EAS, the $\tau$-lepton must emerge with a trajectory with an elevation angle $\beta_{\rm tr}$ ($\theta_{\rm tr}$ from the local zenith) that is within a factor of the Cherenkov angle $\theta_{\rm Ch}$ from the line of sight with the detector, namely, $\beta_{\rm tr}=\beta_{v}\pm\theta_{\rm Ch}$. Hence, at any given location within the field of view of the detector, only those neutrinos with specific trajectories through the Earth will be detectable. These trajectories determine the column depth of material neutrinos must traverse as they propagate through the Earth. Similar arguments can be made for instruments that make use of the radio detection technique.

The Cherenkov angle in air at sea level is $\theta=1.5^{\circ}$, so as Fig. 1 shows, $\beta_{\rm tr}\simeq\beta_{v}$ and $\theta_{\rm tr}\simeq\theta_{v}$. The geometry will be discussed in more detail in Sec. IV.1.

Assuming a spherical Earth, the chord length $c$ of the path through the Earth for a neutrino trajectory with elevation angle $\beta_{\rm tr}$ is given by

[TABLE]

where $R_{E}=6371$ km is the radius of the Earth. The column depth through the Earth for a given trajectory with $\beta_{\rm tr}$ is given by

[TABLE]

where $s$ is the distance along the path traversed through the Earth and $\rho\left(s,\beta_{\rm tr}\right)$ is the density of the Earth along the path length at $s$ for the trajectory. We distinguish between column depth and chord length because the neutrino interaction length and the energy losses for the produced $\tau$-lepton depend on the column depth whereas the probability of the $\tau$-lepton decaying prior to emerging through the surface depends on the remaining distance to the surface.

For the density of the Earth, we use a multi-shell model based on the average radial profile given by the Preliminary Reference Earth Model (PREM) parameterization in Ref. Dziewonski and Anderson (1981). For neutrino trajectories with angles $\beta_{\rm tr}\leq 50^{\circ}$, we find that a model consisting of seven shells of constant density is a good approximation to the PREM column depth, as shown in Fig. 2 with solid and dashed lines, respectively. The density parameters are listed in Table 1 as a function of vertical depth $d=R_{E}-r$, where $r$ is the radial distance from the center of the Earth. To evaluate the column depth, $r$ is related to $s$ and $\beta_{\rm tr}$ by the law of cosines.

II.2 Neutrino interactions

The first step in the evaluation of skimming tau neutrinos is their interactions via charged current or neutral current weak interactions. The high energy neutrino cross section depends on the small-$x$ behavior of the structure functions of the nucleon targets. There is an extensive literature on the high energy extrapolations of the neutrino-nucleon cross section Gandhi et al. (1996, 1998); Reno (2005); Jeong and Reno (2010); Connolly et al. (2011); Cooper-Sarkar et al. (2011); Arguelles et al. (2015). We use as a standard the leading order neutrino cross section evaluated with next-to-leading order (NLO) parton distribution functions (PDFs) for free protons provided by the nCTEQ group (nCTEQ15-1) Kovarik et al. (2016), adjusted for an isoscalar target. It has been shown that the NLO matrix element squared changes the cross section by less than 5% compared to using the leading order matrix element squared with NLO PDFs, a much smaller correction than the uncertainty associated with the choice of small-$x$ extrapolations of the PDFs Jeong and Reno (2010).

Alternate neutrino cross sections can be used to assess uncertainties. Two neutrino cross sections are evaluated with small-$x$ extrapolations based on the Abramowicz et al. (ALLM) extrapolation Abramowicz et al. (1991); Abramowicz and Levy (1997) of the electromagnetic structure function $F_{2}(x,Q^{2})$ and the Block et al. (BDHM) extrapolation Block et al. (2013a, b, 2014) of $F_{2}$. As discussed below, these two expressions for the electromagnetic structure function will be used for the photonuclear energy loss of the $\tau$-lepton . Details of the correspondence between tau energy loss and the tau neutrino cross section appear in ref. Jeong et al. (2017).

With the cross section and differential energy distributions, we use the column depth determined from the Earth density shells to find the interaction point of the neutrino. High densities are favorable for neutrino conversions to taus, however, tau energy loss in dense materials promotes tau decay.

For reference, in Fig. 3 we show the neutrino charged current cross sections using our standard nCTEQ15-1 Kovarik et al. (2016) cross section. and the two alternate cross sections. Fig. 4 shows the ratio of the column depth for fixed angles, $X(\beta_{\rm tr})$, to the nCTEQ15-1 neutrino charged current interaction length as a function of neutrino energy. This shows the importance of attenuation at large angles.

Our Monte Carlo program propagates the neutrino along the chord with repeated interactions of neutral currents as applicable, until a tau is produced (or not) over the neutrino trajectory. The outgoing lepton energy is determined by the differential weak interaction cross section $d\sigma/dy$ where $y=(E_{\nu_{\tau}}-E_{\ell})/E_{\nu_{\tau}}$ where $\ell=\nu_{\tau}$ or $\ell=\tau$ for neutral current and charged current interactions, respectively.

II.3 Tau propagation

Taus produced by charged current interactions of tau neutrinos with nucleons then propagate along the trajectory through the Earth. The lifetime, including the $\gamma$-factor, influences the distance traveled, as does the electromagnetic energy loss. For reference, the time dilated decay length of a tau of energy $E_{\tau}$ is $\gamma c\tau\simeq 5\ {\rm km}[E_{\tau}/(10^{8}\ {\rm GeV})]$ given $c\tau=87.11$ $\mu$m and $m_{\tau}=1.776$ GeV/$c^{2}$.

Electromagnetic energy loss by charged leptons $\ell$ in transit through materials comes from ionization, bremsstrahlung, electron-positron pair production and photonuclear interactions Lohmann et al. (1985); Lipari and Stanev (1991); Dutta et al. (2001, 2005); Armesto et al. (2008); Bigas et al. (2008); Koehne et al. (2013); Alvarez-Muñiz et al. (2018b). The average energy loss per unit column depth $X$ for a charged lepton $\ell$, here with $\ell=\tau$, is written

[TABLE]

At high energies, the ionization loss characterized by $a_{\tau}$ is small. The energy loss parameter $b_{\tau}$ is written as

[TABLE]

where the $b_{\tau}^{i}$ are each calculated in terms of the differential cross section for tau electromagnetic scattering

[TABLE]

weighted by the inelasticity parameter $y=(E_{\tau}-E_{\tau}^{\prime})/E_{\tau}$ for incident lepton energy $E$ and outgoing energy $E^{\prime}$. The $b_{\tau}^{i}$ that characterize $\langle dE/dX\rangle$ depend weakly on energy. It can be shown that $b_{\ell}^{\rm brem}$ scales with lepton mass as $\sim 1/m_{\ell}^{2}$, while for pair production and photonuclear interaction, $b_{\ell}^{\rm pair,nuc}\sim 1/m_{\ell}$ Tannenbaum (1991); Reno et al. (2005). For muons, the three $\beta_{\mu}^{i}$ are similar in scale Lohmann et al. (1985), however, for tau energy loss, our interest here, only pair production and photonuclear interactions are important.

We can use the $b_{\tau}$ in eq. (3), written in terms of the column depth, to understand qualitative features of tau propagation in materials. The density along chord $c$ accounts for most of the differences in materials of interest here (7 density shells), however, there are some material dependent corrections. For pair production, dominated by coherent scattering on the nucleus, the differential cross section scales as the nuclear charge $Z^{2}$, normalized by $A$ in the definition of $b_{\tau}^{\rm pair}$. For photonuclear interactions, the scaling is with $Z/A$, so to a good approximation $b^{\rm nuc}_{\tau}$ is independent of material. We use two sets of parameters for the calculation below: one set of parameters for water, and a second set evaluated with standard rock $Z$ and $A$ to use for rock and the other density shells. Of course, for tau propagation, the density of each shell is included in the column depth for that shell.

We remark that electromagnetic energy loss of taus by photonuclear interactions, at high energies, is tied to the small $x$, and the $Q\sim 1$ GeV regime of the electromagnetic structure function $F_{2}(x,Q)$. As with the high energy neutrino cross section, small-$x$ extrapolations of $F_{2}$, as yet unconstrained by experiment, are required. Different extrapolations give different values for $b_{\tau}^{\rm nuc}$. In Fig. 5, we show how two extrapolations of $F_{2}$ translate to $\beta_{\tau}^{\rm nuc}$ to characterize theoretical uncertainties in the electromagnetic energy loss. We show the Abramowicz et al. (ALLM) form Abramowicz et al. (1991); Abramowicz and Levy (1997), and the Block, Durand, Ha and McKay (BDHM) parameterization Block et al. (2013b, a, 2014). The ALLM extrapolation of $F_{2}$ gives a more strongly energy dependent $b_{\tau}^{\rm nuc}$ than the BDHM parameterization, however, in both cases, the growth with energy is only (approximately) logarithmic in tau energy. The values of $b_{\tau}^{\rm pair}$ for rock and water are also shown in Fig. 5. A more complete discussion of connection between tau energy loss and the small-$x$ extrapolations of $F_{2}$ appears in, e.g., ref. Jeong et al. (2017).

It has been known for some time that stochastic energy loss effects are important, especially for muons Lipari and Stanev (1991). At high energies, stochastic losses from bremsstrahlung and the photonuclear interaction dominate muon energy loss in IceCube Aartsen et al. (2016). For $\tau$-leptons , $\beta_{\tau}^{\rm brem}\ll\beta_{\tau}^{\rm nuc},\beta_{\tau}^{\rm pair}$. As Fig. 5 shows, $b_{\tau}^{\rm nuc}>b_{\tau}^{\rm pair}$, so we propagate $\tau$-leptons stochastically. We account for stochastic effects using the procedure outlined in Refs. Lipari and Stanev (1991); Antonioli et al. (1997); Dutta et al. (2001); Jeong et al. (2017), in a one-dimensional approximation. Since we are at high energies, a one-dimensional approach is appropriate. The procedure involves rewriting

[TABLE]

We account for energy loss with a continuous term for $y\leq y_{\rm cut}$ and simulate tau energy loss stochastically for $y_{\rm cut}<y\leq 1$. We take $y_{\rm cut}=10^{-3}$, a value that reliably includes stochastic energy loss for muons Antonioli et al. (1997); Chirkin and Rhode (2004), and by extension, $\tau$-leptons Dutta et al. (2001). It has been shown that an effective $b_{\tau}$ that describes the average energy of emerging taus from stochastic propagation is typically larger than what one obtains from eq. (5) Dutta et al. (2005).

II.4 Tau decays and regeneration

Over large column depths, tau neutrino regeneration from $\nu_{\tau}\to\tau$ charged current conversion and $\tau\to\nu_{\tau}$ decays is an important effect Halzen and Saltzberg (1998); Beacom et al. (2002); Fargion (2004); Alvarez-Muñiz et al. (2018b). The tau neutrino from a tau decay has approximately one-third of the tau energy. This neutrino then has the potential to interact with the remaining column depth, producing a tau. Each step in the regeneration process reduces the tau energy.

One can make an approximate model of the energy distribution of neutrinos from tau decays for the two- and three-body decays $\tau\to\nu_{\tau}\ell\nu_{\ell},\ \tau\to\pi\nu_{\tau},\ \tau\to\rho\nu_{\tau}$ and $\tau\to a_{1}\nu_{\tau}$, accounting for $\sim 85\%$ of the tau decay width Barr et al. (1988); Pasquali and Reno (1999). The energy distribution of the tau neutrino can be evaluated directly where the form of the distribution is Gaisser et al. (2016):

[TABLE]

for $y\equiv E_{\nu_{\tau}}/E_{\tau}$, given that the taus produced by neutrino interactions are left-handed and antineutrinos produce right-handed $\tau^{+}$. The detailed formulas for $g_{0}$ and $g_{1}$ appear in, e.g., Ref. Gaisser et al. (2016) and the appendix of Ref. Bhattacharya et al. (2016). For three-body leptonic decays,

[TABLE]

In fact, the expression for the three-body purely leptonic decay is a good representation of the energy distribution of tau neutrinos from tau decays evaluated in PYTHIA 8 Sjöstrand et al. (2015). In Fig. 6, we show with solid line histogram the sum of the decay channels from PYTHIA. The solid curve shows just the decay $\tau\to\nu_{\tau}\ell\nu_{\ell}$, normalized to one. Because of the correspondence shown in Fig. 6 between PYTHIA and the analytic leptonic decay distribution, we approximate the full decay distribution of the tau by the semileptonic distribution. For reference, $\langle y\rangle=\langle E_{\nu_{\tau}}/E_{\tau}\rangle\simeq 0.4$ in Fig. 6. Depolarization can in principle occur through multiple scattering of the $\tau$’s, suppressing the term $g_{1}(y)$ in Eq. (7). We neglect this effect here. The spectrum of regenerated taus is not very sensitive to the details of the energy distribution.

To model regeneration in the propagation of neutrinos and taus in the Earth for angles relative to the horizon of $\beta_{\rm tr}\leq 25^{\circ}$, we account for up to five tau decays to neutrinos, where at each step, the tau neutrino is propagated to determine whether or not it reinteracts to produce a tau. This means we simulate up to six neutrino charged current interactions. Any number of neutrino neutral current interactions are included.

The number of interactions becomes progressively more important as a function of increasing energy and angle $\beta_{\rm tr}$. Alvarez-Muniz et al. Alvarez-Muñiz et al. (2018b) have evaluated the number of charged current events for fixed initial tau neutrino energy. They find that the average number of charged current interactions for $E_{\nu}=10^{10}-10^{12}$ GeV is less than $\sim 3$ for $\beta_{\rm tr}\leq 40^{\circ}$. We checked that 5 regeneration steps are sufficient for our purposes here: the tau exit probabilities are not significantly changed by including the last regeneration step for angles less than $20^{\circ}$. For angles larger than those accessible by a POEMMA detector, the final regeneration step does not have a large impact on the exit probability, even at high energy. The exit probability for $E_{\nu_{\tau}}=10^{10}$ GeV changes by 2% for $\beta_{\rm tr}=25^{\circ}$ and by 7% for $\beta_{\rm tr}=35^{\circ}$.

II.5 Flux independent $\nu_{\tau}\to\tau$ results

To illustrate the effects of regeneration, we show the tau exit probability for fixed incident tau neutrino energies: $10^{7}$, $10^{8}$, $10^{9}$ and $10^{10}$ GeV, as a function of angle $\beta_{\rm tr}$. The dashed lines Fig. 7 show the exit probability neglecting regeneration, while the solid lines include regeneration. As has recently been emphasized by Alvarez-Muniz et al. Alvarez-Muñiz et al. (2018b), regeneration is important for all but the smallest values of $\beta_{\rm tr}$.

The details of the tau exit probabilities are difficult to model simply, but the case without regeneration follows from a few factors. The probability for the tau to exit the Earth depends on the attenuation of the neutrino flux and the ratio of the tau range to the neutrino interaction length. This approximation works well for low energies (at or below $E_{\nu_{\tau}}\sim 10^{8}$ GeV) where the range is essentially the time dilated decay length. For higher energies, tau energy loss reduces the range relative to the decay length, and regeneration enhances the large angle exit probabilities, albeit with lower energy taus. Further details of the enhancement of exit probabilities with regeneration are discussed in Appendix B. For a range of 12.5 km, $\rho=\rho_{w}$ for $\beta_{\rm tr}\leq 14^{\circ}$, and for 22 km, trajectory is solely in water for $\beta_{\rm tr}\leq 8^{\circ}$.

The BDHM energy loss model yields somewhat higher exit probabilities at high energies. For $E_{\nu_{\tau}}=10^{8}$ GeV, the enhancement is only a few percent. At this energy, the tau range is only slightly smaller than the decay length. The enhancement of the BDHM evaluation relative to the ALLM electromagnetic energy loss model is of order 25% for $E_{\nu_{\tau}}=10^{9}$ GeV and of order 50% for $10^{10}$ GeV. These enhancements are only due to a smaller $b_{\tau}^{\rm nuc}$ relative to that of the ALLM energy loss model. The neutrino cross section is the same in these comparisons, set to the standard model evaluation with the nCTEQ15-1 PDFs.

The energy distributions of directly produced taus and regenerated taus are combined in the curves in Fig. 8 for $E=10^{9}$ and $10^{10}$ GeV. Plotted are the relative probabilities, the probability for a given bin normalized by the total exit probability shown in Fig. 7, for taus to emerge as a function of energy fraction $z_{\tau}=E_{\tau}/E_{\nu_{\tau}}$. For low angles, the contributions are mainly from direct neutrino production of taus that emerge, while at larger angles, the lower tau energies reflect the fact that regenerated taus dominate the number of taus that exit the Earth. The tau energy fraction upon exit depends on both the energy of the initial tau neutrino and the angle. For example, for $\beta_{\rm tr}=20^{\circ}$, the typical tau energy fraction is a few per cent of the tau neutrino energy for $E_{\nu}=10^{9}$ GeV, but only a few tenths of a percent for $E_{\nu}=10^{10}$ GeV, even though the overall tau exit probabilities for these two energies differ by less than a factor of 2. Tables of tau exit probabilities and the relative probability for an exiting tau with energy $E_{\tau}$ given the incident neutrino energy $E_{\nu_{\tau}}$ and angle $\beta_{\rm tr}$ are used in the Monte Carlo evaluation of the aperture and sensitivity independent of neutrino flux. Representative tables are included in Appendix B.

II.6 Flux dependent $\nu_{\tau}\to\tau$ results

The exit probabilities $P_{\tau}^{\rm exit}$ shown in Fig. 7 are for fixed neutrino energies, integrated over all exiting tau energies. We are also interested in the spectrum of exiting taus given an incident spectrum of neutrinos. As Fig. 8 illustrates, the taus exiting with energy $E_{\tau}$ come from neutrinos with higher energies. We can also evaluate the ratio of the flux of emerging taus to an incident isotropic tau neutrino flux. To begin, we show the ratio $F_{\tau}(E)/F_{\nu}(E)$, the ratio of out-going taus to incident neutrinos, each at the same energy $E$. We denote the ratio by the terminology “transmission function” to emphasize that it includes not just the conversion of tau neutrinos to $\tau-$leptons, but also requires that the $\tau$-leptons exit the Earth before decaying. We use the cosmogenic neutrino fluxes of Kotera et al. in Ref. Kotera et al. (2010) to exhibit the energy and angular dependence of the tau transmission functions.

The range of cosmogenic neutrino fluxes in Ref. Kotera et al. (2010) follow from reasonable inputs to the evolution of the source emissivity, maximum acceleration energy, chemical composition and galactic to extragalactic transition model. Injection spectra and overall normalizations for the cosmic rays were adjusted to best fit the Auger data, then used to predict the associated neutrino spectra Kotera et al. (2010). The neutrino fluxed from six combinations of these inputs are shown in Fig. 9. Plotted is the sum over all flavors of neutrinos plus antineutrinos from all sources, with the differential flux $F(\nu)$ scaled by the neutrino energy squared. The tau neutrino plus antineutrino flux is one third of the flux shown in Fig. 9.

In Fig. 9, curves 1, 5 and 6 use a uniform source emissivity (no evolution of the sources) up to redshift $z=8$. Curves 2 and 3 use a star formation rate (SFR1) according to Hopkins and Beacom Hopkins and Beacom (2006). Curve 3 has an adjusted gamma ray burst evolution (GRB2) following Le and Dermer Le and Dermer (2007). Curve 4 in Fig. 9 uses the evolution of Faranoff-Riley type II galaxies (FRII) of Wall et al. in Ref. Wall et al. (2005). A mixed cosmic ray composition is used in curves 1 and 2, pure protons in 3 and 4, an iron rich composition in curve 5 and pure iron in curve 6. The maximum proton acceleration energy is $E_{p,{\rm max}}=10^{11}$ GeV in curves 1, 2 and 6, $10^{10}$ GeV for curve 5, and $10^{12.5}$ GeV in curves 3 and 4. For mixed composition $E_{Z,{\rm max}}=ZE_{p,{\rm max}}$. Curves 3 and 4 use a dip model for the transition from galactic to extragalactic sources.

For the results presented here, we use two of these representative cosmogenic neutrino fluxes, those labeled by 1 and 4, to evaluate transmission functions that depend on angle. Flux 4 is excluded by Auger in the energy range of $E_{\nu}\sim 4\times 10^{8}-4\times 10^{9}$ GeV Zas (2018); Aab:2019auo . We use it here to illustrate the effect of a harder high energy cosmogenic spectrum than that of flux 1. Fig. 10 shows the transmission ratios for fluxes 1 and 4, with the ALLM energy loss model and the standard model neutrino cross section. The dashed histograms show the transmission functions without regeneration, while the solid histograms include the effects of regeneration. Regeneration effects increase as $\beta_{\rm tr}$ increases, as expected.

The transmission functions for the two fluxes show similar features. The increase in transmission function for $E=10^{7}-10^{8}$ GeV is largely due to the $\gamma$-factor in the decay length and the increase in the cross section with energy. Some regeneration effects are evident already at $\beta_{\rm tr}=5^{\circ}$. At higher energies, the flattening in the ratio occurs because the tau range saturates. Eventually, attenuation of the neutrino flux cuts off the transmission ratio. For $E_{\tau}$ larger than a few $\times 10^{8}$ GeV, the transmission ratio is primarily from a single charged-current interaction. Multiple interactions in regeneration feed down the neutrino flux to a lower emerging tau flux. The harder neutrino spectrum of flux 4 at high energies yields a larger ratio of the tau flux to the incident neutrino flux 4 compared to the ratio for flux 1 for energies above a few times $10^{7}$ GeV to a few times $10^{8}$ GeV, and for small angles $\sim 1^{\circ}-5^{\circ}$, above $E\sim 10^{9}$ GeV. This can be seen in a comparison of the upper and lower panels of Fig. 10.

In Fig. 11, we show $EF_{\tau}(E)$ rather than the transmission function for flux 1 to illustrate the difference in the energy behavior of exiting $\tau$-leptons compared to incident tau neutrinos. The figure comes from using the ALLM energy loss model, again for fixed angles $\beta_{\rm tr}$ relative to the horizon. The much larger incident isotropic tau neutrino flux is scaled by a factor of 1/10.

The energy loss model makes some difference in the predictions. In Fig. 12, the ALLM model results are shown with the solid histograms while the dashed histograms are results using the BDHM model for tau electromagnetic energy loss, both with standard model (SM) neutrino-nucleon cross section. The parameter $b_{\tau}^{\rm nuc}(E)$ evaluated using BDHM is smaller than for ALLM, so tau energy loss at high energies is smaller for BDHM than ALLM evaluations. This effect accounts for the difference at high energies. We note, however, that we use stochastic energy loss rather than $\langle dE_{\tau}/dX\rangle=-b_{\tau}E$ for the tau energy loss to better model the exiting tau energy after transport through the column depth $X$.

Below $E_{\tau}=10^{8}$ GeV, there is little difference in the exiting tau fluxes for a fixed incident neutrino flux because the main feature is that taus are produced in the final few kilometers before exiting the Earth. The predicted tau fluxes with the two energy loss models differ by $\sim\pm 15\%$ for energies below $E_{\tau}=10^{8}$ GeV. Between $E_{\tau}=10^{8}-10^{9}$ GeV, the ratio of the flux prediction from BDHM energy loss grows to a factor of $\sim 1.7$ relative to the ALLM energy loss prediction, increasing further, by more than a factor of 2, as the tau energy increases.

By changing the cross section for neutrino interactions, the variation in the predictions at high energies is wider than if our default SM neutrino-nucleon cross section is used, as shown with the error band in Fig. 12. The error band shows the minimum and maximum exiting tau flux where we also use BDHM extrapolations for both the energy loss and the neutrino cross section, and ALLM extrapolations for both the energy loss and the neutrino cross section.

III Tau air showers

III.1 Tau decays in the atmosphere

The signals that could be detected by sub-orbital or space-based instruments come from air showers produced by tau decays in the atmosphere. The characteristics of the air shower discussed below in Sec. III.2 depend on the tau emergence angle and the altitude at which the tau decays and the shower begins. For a given emergence angle, $\beta_{\rm tr}$, the probability $P_{\rm decay}$ that a tau with energy $E_{\tau}$ will decay at an altitude, $a$, is given by

[TABLE]

where $s\left(a,\beta_{\rm tr}\right)$ is its path length through the atmosphere as a function of $a$ and $\beta_{\rm tr}$ and $p_{\rm decay}(s)=\exp(-s/\gamma c\tau)/(\gamma c\tau)$. Tau energy loss in the atmosphere will be small, so we neglect it here. Fig. 13 shows the path length as a function of altitude, derived from geometry and $\beta_{\rm tr}$ to be

[TABLE]

for selected emergence angles $\beta_{\rm tr}$ between $1^{\circ}-20^{\circ}$. Fig. 14 shows the tau survival probability as a function of altitude for two tau energies: $10^{8}$ GeV (upper figure) and $10^{10}$ GeV (lower figure). Note that at higher energies, the survival probability changes slowly with altitude, so the $y$-axis in the lower figure is linear, whereas it is logarithmic in the upper figure.

In addition to the emergence angle and the altitude, the characteristics of the air shower also depend on the amount of energy that goes into the shower, $E_{\rm shr}$. The $\tau$-lepton is massive enough to produce quarks through its decay, giving rise to a variety of possible final states with substantial variations in the amounts of energy being given to daughter particles that will produce extensive air showers (EASs), namely, electrons and hadrons. Full-scale air shower simulations that would include Monte Carlo simulations of $\tau$-lepton decays is beyond the scope of this paper. Instead, for the purposes of modeling the $\tau$-lepton air showers, we scale the energy of the shower as a fraction of the tau energy. To that end, we performed Monte Carlo simulations of tau decays in PYTHIA 8 and determined that the mean fraction of the tau energy that is given to showering daughters is approximately $\sim 50$% (see Fig. 30 in Appendix C), so for all of the results shown below, $E_{\rm shr}=0.5\,E_{\tau}$.

III.2 Air shower modeling

Once the emergent $\tau$-lepton energy is determined, the resultant EAS needs to be generated based on the location of the $\tau$-lepton decay in the atmosphere. We neglect the decay channel $\tau\rightarrow\mu\nu_{\mu}\nu_{\tau}$, with a branching fraction of 17%. The muonic channel, while interesting, has different shower characteristics than the hadronic and electromagnetic induced extensive air showers modeled in this paper (Stanev and Vankov, 1989).

The EAS modeling philosophy we employ uses the established technique of using a parametric model to describe the average EAS development based on shower universality arguments Hillas (1982a, b); Nerling et al. (2006). We develop the shower and generate the Cherenkov light in a modeled atmosphere and take into account the dominant attenuating atmospheric processes when propagating the signal to the detector. The unique nature of modeling the optical Cherenkov signal induced by upward-moving EASs and measuring the signal using a space-based instrument motivates this philosophy as a well-defined initial simulation. This traditional modeling approach also offers a method to assess the use of simulation packages such as CORSIKA Heck et al. (2012) that were developed for ground-based measurements but require significant modifications to adapt their use for upward-moving rather than down-going EAS modeling Ahnen et al. (2018); Otte (2019).

A detailed atmospheric model is required to define the EAS development, the beamed Cherenkov light emission and the Cherenkov light attenuation based on the optical depth between the EAS and observation point. We employ a static, baseline model for the definition of the atmospheric index of refraction (NAir) as a function of altitude based on that given by Hillas Hillas (1982a), which provides NAir as a function of temperature and atmospheric overburden, g/cm2. We use the model of Shibata (Gaisser, 1990) to define the overburden.

The Cherenkov light attenuation includes the effects of Rayleigh scattering Sokolsky (1989), Mie (aerosol) scattering and ozone absorption. A model for calculating the wavelength dependence of Mie scattering uses the data presented by Elterman in Ref. (Elterman, 1968), which also defines the atmospheric aerosol profile. The Earth’s ozone layer efficiently attenuates optical signals at shorter wavelengths (\lambda\mathrel{\hbox{\hbox to0.0pt{\lower 3.22916pt\hbox{\sim}\hss} \kern-3.00003pt\raise 1.72218pt\hbox{<}}}330 nm). An ozone attenuation model (Krizmanic, 1999) is used with an altitude dependent profile derived from Total Ozone Mapping Spectrometer (TOMS) measurements (McPeters et al., 1996). This composite parametric atmospheric model easily accommodates the underlying spherical geometry.

For the results presented in this paper, the EAS are modeled using the Greisen parameterization Hillas (1982a) with the optical Cherenkov signal calculated over the 200 - 900 nm wavelength band Fernow (1986). The EAS and Cherenkov signal development is generated in 100 meter linear increments as a function of Earth-emergence angle ($\beta_{\rm tr}$) and altitude $a$ of the $\tau$-lepton decay using the inherent spherical geometry. The 3-dimensional nature of the EAS development is modeled using Hillas’s parameterization of the angular and energy spectra of the EAS-electron distributions as a function of shower age (Hillas, 1982b), with the Cherenkov angles and energy thresholds based on the local index of refraction. As noted above, the EAS development simulation results rely on shower universality arguments. CONEX simulations Bergmann et al. (2007) have been used to test our approach for upward EAS modeling, and the comparison shows very good agreement up to an altitude of 15 km, above which CONEX modestly enhances the particle yields.

Fig. 15 shows the Cherenkov profile (upper figure), in photons/m2, for a 100 PeV EAS initiated at sea level with an Earth-emergence angle of $\beta_{\rm tr}=10^{\circ}$ and the attenuated Cherenkov light spectrum (lower figure), both calculated at 525 km altitude. The photon density is approximately constant within the Cherenkov ring (horn-like structure in the profile), but for the Cherenkov light that falls outside the ring, the flux drops off as a power law with radius. Note that at 525 km altitude, the Cherenkov ring spans $\sim$ 100 km and the photon density is quite substantial for a 100 PeV EAS. The Cherenkov ring is very pronounced in the profile due to fact that for modest $\beta_{\rm tr}$ the EAS completes its development at low altitudes, $<5$ km in this case where the Cherenkov emission angle varies only by $\sim 0.2^{\circ}$. The top panel in Fig. 15 also shows a profile function ANe that describes the Cherenkov spatial profile as a flat top with an exponential falloff. This function is used to describe the Cherenkov signal intensity in the $\tau$-lepton EAS simulation for the results presented in this paper. Details of this implementation are in Appendix C.

The intensity and spectrum of the Cherenkov light delivered to a specific altitude is a function of the $\tau$-lepton energy, Earth-emergence angle $\beta_{\rm tr}$, and attenuation in the atmosphere, which can be severe for small values $\beta_{\rm tr}$ mainly due to the effects of the low-altitude aerosols. The interplay between the EAS development and Cherenkov light attenuation is shown in Fig. 16 where the simulated intensities and Cherenkov spectra are shown for 100 PeV EAS with $\beta_{\rm tr}=5^{\circ}$ but with the EAS initiated at different altitudes. The EAS energy 100 PeV is chosen for comparison purposes. At the lowest altitudes, aerosol absorption decimates the Cherenkov intensity and pushes the spectrum towards the longest wavelengths. However, the exponential nature of both the aerosol layer ($\sim$ 1 km scale height) and atmosphere itself ($\sim$ 8 km scale height) leads to signals with higher Cherenkov intensities and spectra peaked at lower wavelengths fairly quickly as a function of EAS starting altitude, due to the nature of the upward-moving $\tau$-lepton EASs.

Eventually the atmosphere becomes too rarefied for complete EAS development at altitudes \mathrel{\hbox{\hbox to0.0pt{\lower 3.22916pt\hbox{\sim}\hss} \kern-3.00003pt\raise 1.72218pt\hbox{>}}} 17 km. The Cherenkov angle becomes significantly reduced at higher altitudes and the Cherenkov threshold energy is increased due to the index of refraction of air (Nair) approaching unity. These combine to lead to a significant reduction in the Cherenkov intensity for EAS that have a large fraction of their development above $\sim$ 17 km altitudes. The energy scale for a $\tau$-lepton to have a non-negligible decay probability at these high altitudes depends on $\beta_{\rm tr}$, but for $\beta_{\rm tr}=5^{\circ}$ and $E_{\tau}=2.5$ EeV, 25% of the $\tau$-leptons will decay above and altitude of 17 km.

The fact that UHE $\tau$-leptons can decay at altitudes comparable to that used by balloon-borne experiments is an interesting phenomena. Initial studies indicate that instruments on scientific balloons at altitudes $a\simeq 33$ km could be in the electromagnetic part of the EAS itself for $\tau$-lepton energies above $\sim$ EeV for Earth emergence angles below $\beta_{\rm tr}\sim 10^{\circ}$. For $\beta_{\rm tr}=1^{\circ}$ and $E_{\tau}=10^{10}$ GeV, at an altitude of $33$ km, the probability for a $\tau$-lepton decay above that altitude is 0.33, while for $\beta_{\rm tr}=10^{\circ}$ for the same $\tau$-lepton energy, the decay probability above $a=33$ km is 0.70. This is illustrated in Fig. 17, which shows the fraction of $\tau$-lepton decays at an altitude larger than 33 km as a function of $\tau$-lepton energy and $\beta_{\rm tr}=1^{\circ}-10^{\circ}$, with colored bands marking increments of 0.1 in the tau decay fraction. For small $\beta_{\rm tr}$, the line of sight distance $v$ is large, so except for the highest energies, almost all of the $\tau$-leptons decay below $a=33$ km. On the other hand, for larger $\beta_{\rm tr}$, $v$ is shorter, so there is the possibility for more $\tau$-leptons to decay at higher altitudes. The accurate calculation of the Cherenkov signal for this case requires a 3-dimensional particle-cascade simulation, e.g., CORSIKA Heck et al. (2012) or Cosmos Roh et al. (2013) but with modifications for the modeling of the Cherenkov signal for upward-moving EASs. However, our simulation approach is valid for balloon altitudes ($\sim$ 33 km) for below an EeV, where the $\tau$-lepton decays below $\sim$ 5 km.

For space-based observations using an $\sim 0.1^{\circ}$ focal plane pixel field-of-view (FoV), a 1-dimensional treatment of the EAS signal is sufficient. This can be understood by considering the relevant distance scales. Assuming the EAS width is defined by a Molière radius value 8.83 g/cm2 for air at STP PDG (a), near sea level the EAS radius is $\sim$ 100 m. From the view of the EAS from 525 km altitude, the 100 m radius is well contained in a single 0.1∘ pixel, even for nadir viewing. For viewing a highly inclined EAS originating near the Earth’s limb, the distance to shower maximum is $>1000$ km (assuming a 525 km orbit) for the Earth-emergence angles ($\beta_{\rm tr}$) with reasonable $\tau$-lepton exit probabilities. This distance scale includes those $>1$ EeV $\tau$-leptons that can decay at altitudes $\sim$ 20 km. While the the EAS radius will widen to $\sim$ 1 km at an altitude of 20 km ($\sim$ 10% atmospheric pressure), the width of the visible portion of the EAS is still well contained in a $0.1^{\circ}$ pixel. In contrast, for observations on balloon-borne experiments ($\sim$ 33 km altitude) or on a mountain-top, such as Trinity ($\sim$ 3 km altitude) the width of an $\tau$-lepton EAS can be large compared the pixel FoV and a 3-dimensional EAS cascade development model is more appropriate.

Thus for the calculation of the Cherenkov signal intensity, spatial extent, and spectrum for low-Earth orbits, we use a parametric model based on our EAS 3-dim Cherenkov approach, which is much more computationally efficient when sampled in a Monte Carlo. The Cherenkov intensities and angles as functions of $\beta_{\rm tr}$ and EAS decay altitude, are tabulated for a fixed, 100 PeV EAS energy in a library format. A profile function fit is used, shown in Fig. 15, to describe the beamed Cherenkov “flattop” signal within the Cherenkov cone, ignoring the horns. As discussed in Appendix C, we scale the intensity as a function of $\tau$-lepton energy and use a mathematical function to account for the increase in the effective Cherenkov acceptance angle for bright signals that place portions of the power-law part of the Cherenkov profile (outside the Cherenkov ring) above the detection threshold of the instrument. This models the increase in acceptance solid angle for brighter EASs.

The Cherenkov angle $\theta_{\rm Ch}^{0}$ as a function of starting altitude, for a 100 PeV showers, is shown by the upper figure of Fig. 18, based on an evaluation of 3-dimensional EAS Cherenkov simulations. Showers that start at lower altitudes have a Cherenkov angle between $\sim 1.2^{\circ}-1.3^{\circ}$. The Cherenkov angle reduces with altitude due to the reduction of the atmospheric index of refraction. The detection of the air shower depends on the photon density at the detector, which in turn depends on the altitude of the detector, the altitude of the start of the air shower, and the Earth-emergence angle. For our evaluation of the sensitivity of instruments with POEMMA performance, we consider a detector at an altitude of $h=525$ km. The photon density within the Cherenkov cone that arrives at such a detector is shown in the lower panel of Fig. 18.

The photon density for $E_{\rm shr}=100$ PeV = $10^{8}$ GeV is the starting point for the photon density at other energies. We approximate the photon density at POEMMA as a function of energy to be

[TABLE]

We discuss below how large photon densities effectively increase the Cherenkov signal acceptance angle, or solid angle, in our evaluation of the POEMMA tau neutrino sensitivity.

We have not included geomagnetic effects in our EAS modeling. The angular spread of electrons in an EAS is due to Coulomb scattering. At the higher elevations, the rarified atmosphere leads to an longer Coulomb scattering length and thus geomagnetic bending can lead to an appreciable enhancement of the angular distributions of the shower Cocconi (1954a, b). We find that the effective Cherenkov acceptance angle enhancement may be as much as a factor of $\sim 2.2$ larger than what we use in our model. The fact that this is a high altitude effect means that the enhancement is for showers with energies higher than $\sim$EeV energies, most of which are already well above detection thresholds. As indicated in Fig. 12, the geomagnetic enhancement has the most impact for \beta_{\rm tr}\mathrel{\hbox{\hbox to0.0pt{\lower 3.22916pt\hbox{\sim}\hss} \kern-3.00003pt\raise 1.72218pt\hbox{<}}}10^{\circ}. Thus, our sensitivity calculations above the EeV scale would be only modestly corrected by including magnetic field effects. Without geomagnetic effects, our sensitivity results are somewhat conservative.

IV Aperture, sensitivity and event rates for POEMMA

IV.1 Aperture

A given instrument’s capability to detect Earth-skimming tau neutrinos is determined by its aperture (or acceptance). Following Motloch et al. (2014), the detector aperture is given by

[TABLE]

where $S$ is the area of the observation region on the Earth, $\hat{r}\cdot\hat{n}=\cos\theta_{\rm tr}$ is the cosine of the angle between the $\tau$-lepton trajectory and the local zenith (see Fig. 1), and $\Omega_{\rm tr}$ is the solid angle of the particle trajectories, of interest in the observable solid angle $\Delta\Omega_{\rm tr}$ around the line of sight. Space-based and balloon borne detectors point near the Earth’s limb for viewing air showers from Earth-skimming tau neutrinos. The observation region is determined by the altitude of the detector and the minimum nadir angle at which the detector can view air showers. For the calculation for POEMMA, we assume a configuration where POEMMA views the Earth $\Delta\alpha\simeq 7^{\circ}$ below the limb. At an orbit altitude of $h=525$ km, this translates to $\theta_{\rm tr}^{\rm min}=90^{\circ}-\beta_{\rm tr}^{\rm max}\simeq 70^{\circ}$. Accounting for just the area of the zone observed on the surface of the Earth over the full $2\pi$ azimuth, the effective area $G$ and geometric aperture $\langle A\Omega\rangle_{\rm geo}$ are

[TABLE]

where the last approximation is for a fixed maximum detection angle $\theta_{d}$ away from the viewing angle, here, the Cherenkov angle. The geometric aperture of this zone is $\langle A\Omega\rangle_{\rm geo}=4,072$ km2sr for the detection angle $\theta_{d}=1.5^{\circ}$, and $\Delta\alpha=7^{\circ}$ at altitude $h=525$ km. The full surface cap under the detector, with $\alpha=0\to\alpha_{H}$, has a geometric aperture approximately twice the zone geometric aperture. For reference, we include in Appendix A the details of the viewing geometry at $h=525$ km as for POEMMA, and for detectors at $h=33$ km and 1000 km.

The geometric aperture is modified by $P_{\rm obs}$, the probability that a tau neutrino with energy $E_{\nu_{\tau}}$ produces a shower that would be detectable. The observation probability $P_{\rm obs}$ is given by

[TABLE]

as also discussed in the context of ANITA in ref. Romero-Wolf et al. (2018b). In Eq. (IV.1), $p_{\rm exit}$ is the differential probability that a $\tau$-lepton of energy $E_{\tau}$ emerges from the Earth given a parent tau neutrino energy of $E_{\nu_{\tau}}$ and an Earth-emergence angle of $\beta_{\rm tr}$, as described in Sec. II.3. Our default energy loss for determining $p_{\rm exit}$ uses the ALLM electromagnetic structure function for photonuclear energy loss. The differential decay probability $p_{\rm decay}\left(s\right)$ is for a tau to decay a distance $s$ from the Earth along its trajectory, discussed in Sec. II.4. The quantity $p_{\rm det}$ is the probability that the emerging tau produces an air shower that would be detected by space-based detector. It depends on the shower energy ($E_{\tau}/2$ here), the position angles of the point of emergence on the Earth $\theta_{E}$ (related to $\theta_{v}$) and $\phi_{E}$, and $\beta_{\rm tr}=\pi/2-\theta_{\rm tr}$ and $\phi_{\rm tr}$. Fig. 19 shows the geometry, with an exaggerated difference between $\theta_{\rm tr}$ and $\theta_{v}$ to show the distinction.

For the detection probability, we approximate $p_{\rm det}$,

[TABLE]

in terms of the Heaviside function $H(x)$:

[TABLE]

The angle $\theta$ is the angle between the $\tau$-lepton trajectory and the line of sight to the detector labeled by $v$ in Fig. 19.

The $\tau$-lepton that decays a distance $s$ from its exit point must be within an “observing cone” of the detector. A two-dimensional projection of the observing cone is shown by the dashed red lines in Fig. 19. To be observed, the tau must decay before it passes outside of the observing cone. The maximum path length for detection of the $\tau$-lepton emerging from the Earth is $s_{\rm win}$, labeled for one of the $\tau$-lepton trajectories in Fig. 19. The value of $s_{\rm win}$ depends on $\theta_{\rm tr}$.

The signal in an instrument is given by the number of photoelectrons which is evaluated from the number density of photons in the Cherenkov cone, multiplied by the area of the detector times the quantum efficiency of the photo-detector for Cherenkov photons. For POEMMA, we assume $A=2.5$ m2 for the effective optical collecting area and 0.2 for the quantum efficiency, based on the average Cherenkov-spectra-weighted photon detection efficiency (PDE) of an typical silicon photomultier (SiPM) Ham :

[TABLE]

For the results shown here, we take $N_{\rm PE}^{\rm min}=10$ with $E_{\rm shr}=0.5\,E_{\tau}$. This choice for $N_{\rm PE}^{\rm min}$ follows from considerations of the night-time air glow which could give false signals of neutrino events, an estimate of the temporal width of the Cherenkov signal based on a geometrical calculation Che , and the largest viewing angles away from the EAS trajectory that leads to measurable signals based on our POEMMA performance model. Note that we have not included the effects of the point-spread-function of POEMMA optics. We assume the Cherenkov signal is effectively delivering into a single 0.084∘ pixel. A model of the air glow background 314 - 900 nm band Mac is used based on VLT/UVES measurements Hanuschik (2003); Cosby et al. (2006) and the van Rhijn enhancement van Rhijn (1921); Roach and Meinel (1955); Best (1965). This model yields an background-spectrum-weighted average PDE of 0.1 using the same SiPM performance for the Cherenkov signal $\langle{\rm PDE}\rangle$. With a collecting area of $A=2.5$ m2 for POEMMA and a 60 ns coincidence window for neutrino events with stereo viewing, the false positive rate due to air glow background in the 314-900 nm band is effective eliminated for N_{\rm PE}^{\rm min}\mathrel{\hbox{\hbox to0.0pt{\lower 3.22916pt\hbox{\sim}\hss} \kern-3.00003pt\raise 1.72218pt\hbox{>}}}10.

We use an effective Cherenkov angle $\theta_{\rm Ch}$ that depends on $\beta_{\rm tr}$, altitude and number of photons. The results we show below in Figs. 20, 21 and 22 have the full integration using Eq. (17) to determine the observation probability. To integrate Eq. (13), we have developed a code that performs the integration via Monte Carlo using importance sampling. More details of the numerical evaluation are discussed in Appendix C. A reasonable approximation overall for the sensitivity is to take $\beta_{\rm tr}=\beta_{v}$ in $P_{\rm obs}$ and integrate $d\Omega_{\rm tr}$ independently to a maximum angle equal to the effective Cherenkov angle. In our consideration of variations to energy loss and neutrino cross section inputs, we use this approximation.

The tau neutrino aperture for POEMMA is shown in Fig. 20. The solid black curve shows the aperture for a $360^{\circ}$ configuration (POEMMA360) with $N_{\rm PE}^{\rm min}=10$. For a configuration with $\Delta\phi=30^{\circ}$ (POEMMA30), the aperture is reduced by a factor of $12$. This is shown with the dashed black curve.

The above calculations do not account for the loss of aperture due to cloud coverage. We do not model the effects of clouds here. As an approximate worst-case scenario of dense, optically opaque clouds over the entire field of view and below an altitude of 5 km, we can reevaluate the effective aperture. Mathematically, this involves multiplying eq. (17) by another Heaviside function: $H\left[a_{\rm decay}-5\mbox{ km}\right]$. The resulting aperture curves are plotted as red lines in the Fig. 20, with the solid red curve for $\Delta\phi=360^{\circ}$ and dashed curve for $\Delta\phi=30^{\circ}$.

IV.2 Sensitivity

The tau neutrino aperture as a function of neutrino energy permits us to evaluate the sensitivity for POEMMA at $h=525$ km altitude to an isotropic tau neutrino flux. The sensitivity over a decade in energy for $N_{\nu}=3$ flavors is given by

[TABLE]

where the factor of $2.44$ events arises from the unified confidence upper limit, (i.e., the upper edge of the two-sided interval for which the lower limit is [math]) at the $90\%$ confidence level Feldman and Cousins (1998). The unified confidence upper limit includes all hypothetical Poisson means for which $n=0$ observed events would be a reasonable realization (i.e., $n=0$ is within the $90$% acceptance interval of observed numbers of events) when drawing from a given Poisson distribution within the unified confidence interval. As such, for signals that are expected to fluctuate about their true values, our use of the unified confidence interval ensures that possible realizations in that observed number of events will be “covered” to the desired confidence level, in this case $90$% (i.e., “coverage probability” of $90$%).111Note that the value of $2.3$ that is often used in the literature excludes values in the interval $\left[2.3,2.44\right]$ for which $n=0$ is a reasonable realization to within $90$% and hence, does not fully cover the $90$% confidence region. In this case, the coverage probability would in fact be less than $90$%. For more in depth discussions, we refer the reader to Refs. Feldman and Cousins (1998); Tanabashi et al. (2018). For the results shown here, we take $t_{\rm obs}=0.2\times\mbox{5 years}$ assuming a twenty percent duty cycle over five years. The assumption for the twenty percent duty cycle is motivated by the relatively large N_{\rm PE}\mathrel{\hbox{\hbox to0.0pt{\lower 3.22916pt\hbox{\sim}\hss} \kern-3.00003pt\raise 1.72218pt\hbox{>}}}10 threshold needed to eliminate the effects of the large air glow background in the 314-900 nm range, e.g., some modest amount of moonlight can be tolerated.

The resulting three-flavor sensitivity curves $E^{2}F_{\rm sens}$ are plotted as black lines in Fig. 21, the dashed curve for $\Delta\phi=360^{\circ}$ and solid curve for $\Delta\phi=30^{\circ}$. The closed circular markers come from evaluating an integral flux scaling like $E_{\nu}^{-\gamma}$ for $\gamma=2$ that yields 2.44 events per neutrino flavor for a given decade of energy centered (on the $\log_{10}$ scale) at the energy of the marker for $\Delta\phi=360^{\circ}$. Thus, we find the normalization $F_{0}$ of

[TABLE]

from

[TABLE]

Plotted in Fig. 21 with the closed circular markers are the values of $E^{2}F_{0}$. There is some variation in the sensitivity to the spectral index $\gamma$, on the order of $\sim\pm 20\%$ for $\gamma=1.5-2.5$.

A range of cosmogenic fluxes from Ref. Kotera et al. (2010) are shown in Fig. 21, with the top of the shaded region bounded by the prediction labeled flux 4 in sec. II.F. We also show 90% CL upper limits for Auger Aab:2019auo , IceCube Aartsen et al. (2018b), and ANITA Gorham et al. (2018), and the predicted sensitivities for ARIANNA Barwick et al. (2015), ARA-37 Allison et al. (2012) and GRAND10k Fang et al. (2018); Alvarez-Muñiz et al. (2018a).

Fig. 22 shows our POEMMA sensitivity calculation as compared to some source models, summed over sources to get a diffuse neutrino flux. Shown with the red band are the all flavor neutrino predictions from newborn pulsar sources in ref. Fang et al. (2014). The dark blue curve labeled AGN is a prediction from active galactic nuclei sources Murase (2017). Neutrinos from galactic clusters with central sources Murase et al. (2008); Fang and Murase (2018) are shown with the light blue curve. The yellow curve shows a prediction from late flares and late prompt emission from gamma ray bursts Murase (2007). Finally, the curve labeled UFA shows a neutrino flux prediction that comes from UHECR photodisintegration within a source from Unger et al. (2015).

The comparison of the sensitivities shows that an azimuthal coverage of $360^{\circ}$ would be required for a 5 year sensitivity to be competitive with other detectors. The POEMMA360 sensitivity in 5 years would probe cosmogenic fluxes in the upper range of predictions in Ref. Kotera et al. (2010) and diffuse astrophysical fluxes from a range of models, e.g., pulsar models of Ref. Fang et al. (2014). The POEMMA360 sensitivity illustrates the benefits of full azimuthal coverage and demonstrates the potential of using the optical Cherenkov signal from upward-moving $\tau$-lepton EAS induced from tau neutrino interactions in the Earth.

IV.3 Flux dependent event rates

In addition to computing the energy-dependent aperture and sensitivity, we also calculate the event rate for a given flux of tau neutrinos above a specified energy. We follow the same procedure as for the flux independent results, however, with a factor of the tau flux and an integration over the energy of the tau (twice the shower energy),

[TABLE]

The tau flux is determined from the transmission functions discussed in sec. II.F.

For the POEMMA $360^{\circ}$ configuration and an observation time of five years with a duty cycle of $\sim 20$%, we find the expected number of events above $E_{\nu_{\tau},{\rm min}}=10^{7}$ GeV is $\sim 0.2$ events in the flux 1 scenario, representing the prediction with a uniform distribution of sources with no evolution, a mixed composition of UHECRs and maximum cosmic ray energy of $E_{\rm max}=100$ EeV. For the flux 4 scenario that has a source evolution following that of Fanaroff-Riley Type II active galactic nuclei and a pure proton UHECR composition with $E_{\rm max}=3160$ EeV, the number of events in the five-year, 20% duty cycle time frame is $\sim 13$ events. A restriction to decays above an altitude of $a=5$ km reduces the number of events, for example, to $\sim 3$ event for Flux 4 in the same time period.

V Discussion

How the tau neutrino and $\tau$-lepton interactions are modeled affects the sensitivity of POEMMA, as does the density model. We have examined some of these features using a simplified evaluation of the aperture and sensitivity for POEMMA, setting $\beta_{\rm tr}=\beta_{v}$ so $H[s_{\rm win}-s_{d}]=1$ in eq. (17). Numerically, the simplified evaluation gives results very close to the full Monte Carlo evaluation, so we used this simplification to study variations in the predictions due to these effects.

The relative benefits of observing upward-going air showers over land and water have been discussed by Palomares-Ruiz et al. (2006). They argue that the very high energy shower rate is significantly enhanced over water compared to over rock. Tau energy loss in water is less than in rock because of the different densities, but the density of rock favors neutrino interactions. In the results shown thus far, the sensitivity is evaluated assuming the final density shell of the Earth is water, according to the PREM model. We can do the same evaluation of the sensitivity assuming the final density shell of the Earth is standard rock. A similar evaluation has been performed in Ref. Alvarez-Muñiz et al. (2018b). We find that most of the energy range to which POEMMA is sensitive is not high enough for the onset of an enhancement of the over-water event rate, in qualitative agreement with the results presented in Ref. Alvarez-Muñiz et al. (2018b). The water versus land effect can be understood by considering the distance scales as a function of energy.

The angles $\beta_{\rm tr}=1^{\circ}-20^{\circ}$ correspond to a range of chord lengths in the final density shell, for example, the whole trajectory of 222 km for $\beta_{\rm tr}=1^{\circ}$ to a final 10.3 km in the outer shell for $\beta_{\rm tr}=20^{\circ}$. For $E_{\tau}=10^{8}$ GeV, the time dilated decay length is 5 km. At energies E_{\nu}\mathrel{\hbox{\hbox to0.0pt{\lower 3.22916pt\hbox{\sim}\hss} \kern-3.00003pt\raise 1.72218pt\hbox{<}}}10^{8} GeV, the produced tau’s lifetime, not energy loss, determines the tau range in the final density shell. To first approximation, there is a benefit to a rock target in the last density shell rather than a water target, since the column depth for neutrino interactions is $X\sim\rho\gamma c\tau$ and the exit probability is $P^{\rm exit}_{\tau}\sim X/\lambda_{\nu}$ for the interaction length (in units of column depth) $\lambda_{\nu}=(N_{A}\sigma_{\nu N})^{-1}$, so $P^{\rm exit}_{\tau}$ is larger for higher density $\rho$.

As the neutrino energy increases, the time-dilated lifetime increases and tau electromagnetic energy loss becomes important. Fig. 5 shows that the tau energy loss parameter is smaller for water than for rock for the $e^{+}e^{-}$ pair production process. In addition, the column depth in water is smaller than in rock for the final shell, so $\Delta E_{\tau}\sim b_{\tau}\Delta X$ is smaller for water (by a factor of $\rho_{\rm water}/\rho_{\rm rock}$) than for rock, allowing more taus to emerge from water than from rock. An evaluation of the sensitivity for an Earth model with the outer shell density set to $\rho_{\rm rock}$, in the approximation that $\beta_{\rm tr}=\beta_{v}$, is shown in Fig. 23 by the dashed black curve, compared to an evaluation with the same approximation with the final density shell with $\rho_{\rm water}$, as in Table I (solid black curve). For showers from taus that emerge from rock, the sensitivity is lower (better) by a factor of $1/2.5$ for $E_{\nu_{\tau}}=10^{7}$ GeV. For $E_{\nu_{\tau}}=10^{8}$ GeV, the reduction is a factor of 1/2. When $E_{\nu_{\tau}}=10^{10}$ GeV, the sensitivities are equal, whether the showers occur over rock or water. For POEMMA detector thresholds, energies above $E_{\nu_{\tau}}=10^{10}$ GeV show only a modest (up to $\sim 30\%$) advantage for observations over water compared to rock.

The different extrapolations of the electromagnetic structure function needed for photonuclear energy loss also shift the sensitivity curve. Our default choice is the ALLM extrapolation of $F_{2}(x,Q^{2})$ in the evaluation of $b_{\tau}^{\rm nuc}$ with a neutrino nucleon cross section that relies on the nCTEQ-1 PDFs ($\sigma_{\rm SM}$), which also extrapolate structure functions. The blue solid curve in Fig. 23 shows the sensitivity with the BDHM extrapolation of $F_{2}(x,Q^{2})$, keeping the neutrino nucleon cross section evaluated with nCTEQ-1 PDFs. The blue shaded region is indicative of the uncertainty associated with the $F_{2}$ extrapolation, keeping the neutrino-nucleon cross section fixed.

Again, at low energy where tau energy loss is not very important, the ALLM and BDHM energy loss evaluations yield nearly identical results. As the energy increases, the smaller value of $b_{\tau}^{\rm nuc}$ with the BDHM extrapolation means less tau energy loss, a larger aperture and a lower sensitivity. The sensitivity curve from the BDHM evaluation is a factor of two lower than the ALLM evaluation at the highest energy shown in Fig. 23.

Changing the neutrino cross section to evaluations using ALLM (BDHM) structure function extrapolation instead of using nCTEQ15-1 in the next-to-leading order QCD calculation, along with the respective $b_{\tau}^{\rm nuc}$, gives results shown with the black (blue) dot-dashed lines. At the highest energy, the change is modest. For $E_{\nu}=10^{9}$ GeV, the lowest curve is a factor of 1.75 below the default curve (the black curve) with ALLM used for energy loss and the standard model neutrino cross section.

To what degree do the POEMMA360 detection characteristics limit the sensitivity? The sensitivity for POEMMA viewing for angles within $\Delta\alpha\sim 7^{\circ}$ below the horizon is shown in Figs. 21 and 22. Increasing the viewing to $\Delta\alpha=9^{\circ}$ only marginally improves the sensitivity. For $\Delta\alpha=15^{\circ}$, corresponding to \beta_{\rm tr}\mathrel{\hbox{\hbox to0.0pt{\lower 3.22916pt\hbox{\sim}\hss} \kern-3.00003pt\raise 1.72218pt\hbox{<}}}31^{\circ}, the sensitivity at $E_{\nu}=10^{7}$ GeV is $\sim 10^{-7}$ for $E_{\nu}^{2}dN/dE_{\nu}$ in Fig. 21. The sensitivity for POEMMA360 would reduce by a factor of 0.42 for $E_{\nu}=2\times 10^{7}$ GeV, and only by a factor of 0.87 for $E_{\nu}=10^{8}$ GeV. The larger values of $\beta_{\rm tr}$ favor the lower neutrino energies relative to higher energies because the lower energy neutrino flux is less attenuated, so a relatively higher fraction of taus emerge to shower at angles that are better detected at low energies than for small $\beta_{\rm tr}$ where the distance $v$ from the point of emergence and the detector is much larger.

The threshold for the number of photoelectrons detected is an important feature of the sensitivity. If $N_{\rm PE}^{\rm min}=5$, or alternatively, if the detection area times quantum efficiency were a factor of two higher, we find that the sensitivity is a factor of 0.10 of the black curve for POEMMA360 in Fig. 21 for $E_{\nu}=2\times 10^{7}$ GeV, a factor of 0.43 lower for $E_{\nu}=10^{8}$ GeV and lower by a factor of $0.60$ for $E_{\nu}=10^{9}$ GeV, for example. The photoelectron threshold is set for POEMMA to achieve a negligible false positive neutrino rate from the night-time air glow background. Lowering the photoelectron threshold may be possible with a restricted wavelength range, however, further work is needed to assess whether or not the sensitivity would be lowered.

The results presented here for the detection of skimming tau neutrinos via upward-going air showers from tau decays, with POEMMA detection as a specific example, come from a Monte Carlo evaluation of neutrino interaction and tau energy loss, then a one-dimensional model of the extensive air shower from the tau decay. A broader program to simulate signals of skimming tau neutrinos is underway to address some of the approximations used here, and to compare with approximations in other approaches. For example, in Ref. Alvarez-Muñiz et al. (2018b), tau energy loss is treated as continuous, whereas we use stochastic energy loss here.

An evaluation of timing of the extensive air showers shows that for detectors more than 16 km away from an extensive air shower with a radius of 100 m, our one-dimensional approximation for the air shower modeling is reliable. A full three-dimensional shower simulation will be useful for showers at larger elevation angles for balloon detection.

The air shower modeling uses the 83% of non-muonic tau decays, with half of the tau energy going into the shower. Incorporation of the energy distribution of $E_{\rm shr}$ is planned. An additional consideration is whether or not air showers associated with high energy muons are detectable Stanev and Vankov (1989). The particle yield from HE muon showers is suppressed compared to that of a similar energy electron (or photon) EAS. The length of a muonic air shower has a much longer slant depth due to this and the long lifetime of the muon. From a practical standpoint, the propagation of the muonic shower could yield an EAS signature at much high altitudes, even for $\tau$-lepton decays that initiate near the Earth’s surface.

To summarize, we have presented a new calculation of the flux and spectrum of Earth emerging $\tau$-lepton from an isotropic flux of cosmic neutrinos, then applied our results to a space-based experiment with a performance modeled on the POEMMA mission. We have illustrated many features of neutrino and tau propagation and tau shower detection. Optical Cherenkov signals from upward-going air showers show promise for detecting cosmic neutrinos, especially below 100 PeV. We find that a POEMMA-like instrument requires a $360^{\circ}$ azimuthal optical Cherenkov coverage to be competitive with other detectors current or planned. While our focus has been on an isotropic diffuse flux, the calculational tools developed here can be applied to searches for individual neutrino sources. For a detector like POEMMA, the ability to quickly reorient the detectors will permit tracking of target-of-opportunity neutrino sources. An assessment of POEMMA’s sensitivity to these sources appears in Ref. Venters et al. (2019).

Acknowledgements.

We thank Y. Akaike, L. Anchordoqui, D. Bergman, C. Guépin, S. Mackovjak, A. Neronov, A. Olinto, P. Privitera, A. Romero-Wolf, F. Sarazin and E. Zas for discussions, and the POEMMA collaboration for motivating this analysis. Our work was supported in part by a US Department of Energy grant DE-SC-0010113, NASA grant 17-APRA17-0066 and NASA award NNX17AJ82G.

Appendix A Geometry of Earth at altitude

The detection of tau neutrino induced tau air showers relies on the Earth as a neutrino converter. In this appendix, we show the angle relations needed to evaluate the geometric aperture in terms of angles labeled in Fig. 1. Air shower detectors at altitude $h$ above the Earth will point near the limb, which is at a viewing angle $\alpha_{H}$ away from the nadir. For POEMMA, we take $h=525$ km. Emergence angles at the limb are related by

[TABLE]

while more generally,

[TABLE]

where $v$ is the path length in the atmosphere,

[TABLE]

and $\theta_{E}$ is the polar angle of the given position on the Earth. These relations are used to find the difference in the viewing angle $\alpha$ relative to the angle to the limb, $\Delta\alpha\equiv\alpha_{H}-\alpha$. For $h=33,\ 525,\ 1000$ km, $\alpha_{H}=84.2^{\circ},\ 67.5^{\circ},\ 59.8^{\circ}$, respectively. Table 2 shows $\Delta\alpha=\alpha_{H}-\alpha$ as a function of $\beta_{E}$ for these three altitudes. The planned POEMMA Cherenkov detector will have a detection viewing angle range of $\Delta\alpha\simeq 7^{\circ}$, when pointed near the limb, so the angular coverage for $h=525$ km is for $\beta_{v}\sim 0\to 20^{\circ}$. In fact, very close to the limb, the signal will be overcome by background, but for $\Delta\alpha>1^{\circ}$, the backgrounds are significantly reduced.

The viewing angle relative to the local zenith, $\theta_{v}$, is given by

[TABLE]

The angles $\theta_{v}=90^{\circ}-\beta_{v}$, $\alpha$ and $\theta_{E}$ are related by

[TABLE]

Absent detector field of view considerations and in the approximation that $\theta_{\rm tr}\simeq\theta_{v}$, the accessible flux from the area of a cap of the spherical Earth below a detector at altitude $h$ comes from a surface area $S$ from $\theta_{E}^{min}=0^{\circ}$ to $\theta_{E}^{max}=\theta_{E}^{H}$. The effective area for air showers from taus emerging from the Earth with angle $\theta_{v}$ relative to the normal to the cap surface is Motloch et al. (2014)

[TABLE]

The effective area of the cap below the detector can be written as

[TABLE]

For a field of view characterized by $\Delta\alpha=7^{\circ}$, the effective area is reduced. The accessible area is a band around the cap. For $\Delta\alpha=7^{\circ}$ and $h=525$ km, $\theta_{E}^{min}=9.9^{\circ}$. We refer to this as the effective area of the “zone,” smaller than the “cap” described by Eq. (31).

The geometric aperture (geometry factor) $\langle A\Omega\rangle_{\rm geo}$ defined in eq. (IV.1), however, depends on $\hat{r}\cdot\hat{n}=\cos\theta_{\rm tr}$. In terms of the angles $\delta$ and $\phi_{\delta}$ that the $\tau$-lepton trajectory makes with respect to the line of sight, and the angle $\theta_{v}$ the line of sight makes with respect to the local zenith,

[TABLE]

Then for $d\Omega_{\rm tr}=\sin\delta\,d\delta\,d\phi_{\delta}$ with the full $\phi_{\delta}$ integral over $2\pi$ and $\delta=0\to\theta_{d}$, the geometric aperture is

[TABLE]

Table 3 compares the cap and zone geometric apertures for several altitudes when $\theta_{d}=1.5^{\circ}$ and $\Delta\alpha=7^{\circ}$. For $h=525$ km, the ratio of apertures for the zone and cap is 0.48 with these assumptions. Fig. 24 shows the geometric aperture of the cap and of the zone (for $\Delta\alpha=7^{\circ}$) as a function of detector altitude for $\theta_{d}=1.5^{\circ}$.

Appendix B Tau exit probability and energy distributions

In this section, we include tables for the $\tau$-lepton exit probabilities for a given tau neutrino energy and tables for the exiting tau energy given a fixed tau neutrino energy and $\beta_{\rm tr}$. Table 4 lists the exit probabilities for fixed energies $E_{\nu_{\tau}}=10^{7}$, $10^{8}$, $10^{9}$ and $10^{10}$ GeV when the ALLM extrapolation is assumed for $b_{\tau}^{\rm nuc}$. The standard model neutrino cross section, as described in sec. II.B is assumed. For reference, we show in Fig. 25 the ratio of the exit probabilities in Table 4 to the exit probabilities without regeneration, namely, assuming a single $\nu_{\tau}\to\tau$ conversion. The probabilities by themselves do not reflect the shift in energy from the multiple interactions when regeneration is important.

Given an exit probability for a given neutrino energy, the outgoing tau energy distribution depends on elevation angle $\beta_{\rm tr}$, as shown, for a few energies, in Fig. 8. In our evaluation of the aperture and sensitivity, instead of the distributions like those in Fig. 8, we use the cumulative distribution functions,

[TABLE]

Figs. 26, 27, 28, and 29 show the cumulative distribution function as a function of the scaled energy $z_{\tau}=E_{\tau}/E_{\nu_{\tau}}$. Tables 5-8 list numerical values for four tau neutrino energies and $\beta_{\rm tr}=1^{\circ}$, $5^{\circ}$, $10^{\circ}$ and $20^{\circ}$.

The average exiting tau energy decreases with an increase in $\beta_{\rm tr}$. Multiple interactions and energy loss are responsible for the shift to a lower exiting tau energy as the incident neutrino energy increases. This means that even though regeneration significantly enhances the exit probability, the energy of the exiting tau is lower, so regeneration does not necessarily translate to a better sensitivity. For example, with the POEMMA360 detection characteristics modeled here, the correction to the sensitivity due to regeneration is at most a $\sim 20\%$ effect.

Appendix C Detection probability

In this appendix, we give more details for our evaluation of the detection probability in our Monte Carlo computer program. The probability of detecting the $\tau$-lepton shower depends on the shower energy, the altitude of the tau decay and the detection angle. In the full Monte Carlo simulation calculating the probability that the emerging tau produces a detectable air shower, $p_{\rm det}$, we include the following requirements:

The trajectory of the parent neutrino, assumed to be collinear with the emerging tau lepton, must be appropriately aligned with the line of sight between the detector and the point of emergence. We require the angle between the trajectory of the neutrino and the line of sight of the detector, $\theta=\lvert\theta_{v}-\theta_{\rm tr}\rvert$ in Fig. 14, to be less than the effective Cherenkov angle, $\theta_{\rm Ch}$. 2. 2.

The tau must decay before it leaves the observation window, the three-dimensional zone that is visible to the detector denoted by the red lines in Figure 19. 3. 3.

The shower from the tau neutrino must be able to trigger the detector, namely, the number of photoelectrons in the detector generated by light from the shower, $N_{\rm PE}$, must be greater than a threshold value, taken to be $10$ for POEMMA.

For the purposes of calculation, we model each requirement using a Heaviside function.

The detection window is determined by the Cherenkov angle, discussed in Sec. III.2. The Cherenkov angle for $E_{\rm shr}=10^{8}$ GeV is approximated by the results shown in the lower plot of Fig. 18. For showers $N_{\rm PE}\gg N_{\rm PE}^{\rm min}$, the effective Cherenkov angle is larger than what is shown in the upper panel of Fig. 18. The tails of the Cherenkov photon density in the upper panel of Fig. 18 show that a Cherenkov angle based on the plateau of the photon density underestimates the width of the Cherenkov signal. We use an effective Cherenkov angle

[TABLE]

This comes from assuming the one-dimensional profile of the shower is approximately Gaussian and scaling by the half-width at $f=N_{\rm PE}/N_{\rm PE}^{\rm min}$ times the maximum.

The observation window requires the pathlength of the $\tau$-lepton from exit point to decay $s_{d}$ to be less than $s_{\rm win}$. To compute $s_{\rm win}$, we consider three cases, where the zenith angle $\theta_{\rm tr}$ of the trajectory of the tau lepton is less than, greater than, or equal to the zenith angle $\theta_{v}$ of the detector line of sight (black line segment labelled “v” in Fig. 19). Considering all three cases, $s_{\rm win}$ is given by

[TABLE]

where $v$ is defined in eq. (27). The quantity $\alpha_{\rm min}$ is the minimum nadir angle of the detector viewing zone, and $\xi$ is the angle opposite the detector’s line of sight, given by

[TABLE]

The value of $\alpha$ can be found using the law of cosines:

[TABLE]

As discussed in Sec. II.5, the number of photoelectrons detected depends on the photon number density, the elevation angle $\beta_{\rm tr}$ and the altitude of the decay, related to $s_{d}$. In Fig. 30, we show the frequency of a given $E_{\rm shr}/E_{\tau}$ from a PYTHIA8 simulation of tau decays. Without the electron decay channel, the average energy of the shower is $\sim 0.6E_{\tau}$, but including the electron channel lowers the average to $\sim 0.5E_{\tau}$. The results in this paper use the approximation $E_{\rm shr}=E_{\tau}/2$. We take $N_{\rm PE}^{\rm min}=10$.

We integrate eq. (13) via Monte Carlo Integration using importance sampling (also known as the inverse transform method) PDG (b). In this sampling method, random variables are drawn from cumulative distribution functions (CDFs) constructed from selected probability density functions (PDFs). Ideally, the selected PDFs would be as similar to the functions being integrated as possible, though normalized, in order to minimize sample variance. In general, for a function, $f\left(x\right)$ being integrated and samples drawn from the selected PDF, $p\left(x\right)$, the Monte Carlo estimator is given by

[TABLE]

where $N$ is the number of samples and $X_{i}$ is the $i^{\rm th}$ drawn random variable. The formula for the Monte Carlo estimator can be verified by taking the expectation of $F_{N}$ over $p\left(x\right)$ over the interval of integration.

For the full Monte Carlo integration, the integrand of eq. (13) includes $\hat{r}\cdot\hat{n}=\cos\theta_{\rm tr}$. Putting $\hat{r}$ and $\hat{n}$ in a frame in which $v$ points in the $\hat{z}$ direction and taking the dot-product, we find that $\cos\theta_{\rm tr}=\cos\theta_{v}\cos\delta-\sin\theta_{v}\sin\delta\sin\phi_{p}$, where $\theta_{v}$ is the local zenith angle of the line-of-sight between the spot on the ground and the detector, $\delta$ is the angle between the particle’s trajectory and the line-of-sight to the detector, $\phi_{p}$ is the azimuthal angle of the particle in the frame in which $v$ points in the $\hat{z}$ direction. Then, eq. (13) becomes

[TABLE]

where $\theta_{E}$ is the zenith angle of the position on the surface of the Earth, $R_{E}$ is the radius of the Earth, and $P_{\rm obs}$ is the observation probability given by eq. (IV.1). For the full Monte Carlo integration, the chosen PDF is

[TABLE]

where $p_{\rm decay}\left(s\right)$ is the probability that the $\tau$ decays after traveling a path length $s$. Then, the Monte Carlo estimator is given by

[TABLE]

where $\mathcal{N}$ is a factor that includes the normalization of $p\left(\theta_{E},\delta,\phi_{p},s\right)$ and the factor $2\pi R^{2}_{E}$.

The POEMMA sensitivity is determined with the full Monte Carlo. A simpler formalism in which $\theta_{\rm tr}\rightarrow\theta_{v}$ and the integration is performed over $d\Omega_{\rm tr}$ yields sensitivities that are reasonably close to the full calculation.