Tau energy loss and ultrahigh energy skimming tau neutrinos

TL;DR

This paper analyzes the propagation and energy loss of ultrahigh energy tau neutrinos and taus, quantifying uncertainties in flux predictions relevant for astrophysical neutrino detection experiments.

Contribution

It provides a comprehensive evaluation of tau energy loss models and their impact on tau flux predictions, including uncertainty analysis for various distances and energies.

Findings

Uncertainty in tau energy loss and neutrino cross section partially cancel, reducing overall flux uncertainty.

Theoretical uncertainty ranges from 30% to factors of 3.8 depending on energy and distance.

Monte Carlo simulations agree with analytic calculations for tau propagation.

Abstract

We consider propagation of high energy earth-skimming taus produced in interactions of astrophysical tau neutrinos. For astrophysical tau neutrinos we take generic power-law flux, $E^{-2}$ and the cosmogenic flux initiated by the protons. We calculate tau energy loss in several approaches, such as dipole models and the phenomenological approach in which parameterization of the $F_2$ is used. We evaluate the tau neutrino charged-current cross section using the same approaches for consistency. We find that uncertainty in the neutrino cross section and in the tau energy loss partially compensate giving very small theoretical uncertainty in the emerging tau flux for distances ranging from $2$ km to $100$ km and for the energy range between $10^6$ GeV and $10^{11}$ GeV, focusing on energies above $10^8$ GeV. When we consider uncertainties in the neutrino cross section, inelasticity in neutrino interactions and the tau energy loss, which are not correlated, i.e. they are not all calculated in the same approach, theoretical uncertainty ranges from about $30\%$ and $60 \%$ at $10^8$ GeV to about factors of 3.3 and 3.8 at $10^{11}$ GeV for the $E^{-2}$ flux and the cosmogenic flux, respectively, for the distance of 10 km rock. The spread in predictions significantly increases for much larger distances, e.g., $\sim 1,000$ km. Most of the uncertainty comes from the treatment of photonuclear interactions of the tau in transit through large distances. We also consider Monte Carlo calculation of the tau propagation and we find that the result for the emerging tau flux is in agreement with the result obtained using analytic approach. Our results are relevant to several experiments that are looking for skimming astrophysical taus, such as the Pierre Auger Observatory, HAWC and Ashra. We evaluate the aperture for the Auger and discuss briefly application to the the other two experiments.

Tables

Table 1. Table 1: Fit to β τ nuc superscript subscript 𝛽 𝜏 nuc \beta_{\tau}^{\rm nuc} in [ 10 − 6 superscript 10 6 10^{-6} cm 2 /g] using Eq. ( 15 ) for E = 10 8 − 10 12 𝐸 superscript 10 8 superscript 10 12 E=10^{8}-10^{12} GeV. For bremsstrahlung, β 0 brem = 7.9 × 10 − 9 superscript subscript 𝛽 0 brem 7.9 superscript 10 9 \beta_{0}^{\rm brem}=7.9\times 10^{-9} cm 2 /g and pair production, β 0 pair = 0.148 × 10 − 6 superscript subscript 𝛽 0 pair 0.148 superscript 10 6 \beta_{0}^{\rm pair}=0.148\times 10^{-6} cm 2 /g, essentially independent of energy for this energy range.

Model	${\beta }_0^{\mathrm{nuc}}$	${\beta }_1$
Soyez-Shadow	0.380	$3.21\times {10}^{-2}$
Soyez-ASW	0.471	$3.91\times {10}^{-2}$
Soyez-GG	0.383	$3.07\times {10}^{-2}$
AAMQS-GG	0.433	$3.46\times {10}^{-2}$
BDHM-Shadow	0.435	$4.04\times {10}^{-2}$

Table 2. Table 2: The average tau range [km] in rock as a function of initial tau energy using our tau transport Monte Carlo.

Energy	Shadow	Shadow	Shadow	GG
	BDHM	ALLM	Soyez	Soyez
${10}^8$	3.46	3.22	3.49	3.46
${10}^9$	12.7	10.1	13.1	12.8
${10}^{10}$	25.0	16.9	26.4	25.6
${10}^{11}$	35.9	21.5	38.4	37.2
${10}^{12}$	44.1	24.2	47.2	45.9

Table 3. Table 3: Parameters for the charged current neutrino cross sections from Eq. ( 24 ) and ( 25 ).

Model	E [GeV]	${\sigma }_0$	${\sigma }_1$	${\sigma }_2$
BDHM	${10}^6\le E\le {10}^9$	0.310	$-4.45\cdot {10}^{-2}$	1.63 $\cdot {10}^{-3}$
	$E>{10}^9$	1.32	$-0.141$	3.92 $\cdot {10}^{-3}$
Soyez	${10}^6\le E\le {10}^9$	0.269	$-3.80\cdot {10}^{-2}$	1.37 $\cdot {10}^{-3}$
	$E>{10}^9$	2.30	$-0.230$	5.91 $\cdot {10}^{-3}$
ALLM	${10}^6\le E\le {10}^9$	0.460	0.361	$-$
	$E>{10}^9$	1.05	0.321	$-$
JR Jeong and Reno (2010)	${10}^8\le E\le {10}^{10}$	0.718	0.348	$-$
	$E>{10}^{10}$	1.78	0.308	$-$

Table 4. Table 4: Parameterizations of ⟨ y ⟩ delimited-⟨⟩ 𝑦 \langle y\rangle using Eq. ( 26 ). The values of ⟨ y ⟩ delimited-⟨⟩ 𝑦 \langle y\rangle from these parameterizations have less than 1% error between E = 10 6 − 10 13 𝐸 superscript 10 6 superscript 10 13 E=10^{6}-10^{13} GeV.

Model	E [GeV]	$y_0$	$y_1$	$y_2$
BDHM	${10}^6\le E\le {10}^8$	0.909	$-5.95\cdot {10}^{-2}$	1.17 $\cdot {10}^{-3}$
		0.654	$-3.35\cdot {10}^{-2}$	5.01 $\cdot {10}^{-4}$
Soyez	${10}^6\le E\le {10}^8$	1.08	$-8.55\cdot {10}^{-2}$	2.07 $\cdot {10}^{-3}$
		0.478	$-2.05\cdot {10}^{-2}$	2.98 $\cdot {10}^{-4}$
ALLM	${10}^6\le E\le {10}^8$	1.17	$-9.99\cdot {10}^{-2}$	2.59 $\cdot {10}^{-3}$
		0.356	$-1.25\cdot {10}^{-2}$	2.27 $\cdot {10}^{-4}$

Table 5. Table 5: Path length and sagitta to correspond to the angle α 𝛼 \alpha relative to the horizon.

$\alpha$ [∘]	Sagitta [km]	Path length [km]
1	0.97	222
2	3.88	445
3	8.73	667
4	15.5	889
5	24.2	1,110
6	34.9	1,330

Table 6. Table 6: β τ nuc superscript subscript 𝛽 𝜏 nuc \beta_{\tau}^{\rm nuc} in [ 10 − 6 superscript 10 6 10^{-6} cm 2 /g] using Eq. ( 35 ) for E = 10 6 − 10 12 𝐸 superscript 10 6 superscript 10 12 E=10^{6}-10^{12} GeV.

Model	${\beta }_0^{\mathrm{nuc}}$	${\beta }_1$	${\beta }_2$
BDHM	0.425	4.04 $\times {10}^{-2}$	1.12 $\times {10}^{-3}$
Soyez	0.371	3.20 $\times {10}^{-2}$	9.54 $\times {10}^{-4}$
Soyez-ASW	0.461	3.90 $\times {10}^{-2}$	1.13 $\times {10}^{-3}$
ALLM	1.02	0.210	1.51$\times {10}^{-2}$

Table 7. Table 7: Couplings for neutrino charged current ( V = W 𝑉 𝑊 V=W ) and neutral current ( V = Z 𝑉 𝑍 V=Z ) scattering.

Process	$m_1$	$m_2$	$g_v^V$	$g_a^V$
$V=W$	$m_d,m_s,m_b$	$m_u,m_c,m_t$	1	1
$V=Z$	$m_d,m_s,m_b$	$m_d,m_s,m_b$	$\sqrt{2}\left(-\frac{1}{2}+\frac{2}{3}{\sin }^2{\theta }_W\right)$	$-\sqrt{2}/2$
$V=Z$	$m_u,m_c,m_t$	$m_u,m_c,m_t$	$\sqrt{2}\left(\frac{1}{2}-\frac{4}{3}{\sin }^2{\theta }_W\right)$	$\sqrt{2}/2$

Table 8. Table 8: Charged and neutral current weak interaction cross sections with the BDHM parameterization and Soyez dipole, with slow rescaling χ = x ​ ( 1 + ( m i + m j ) 2 / Q 2 ) 𝜒 𝑥 1 superscript subscript 𝑚 𝑖 subscript 𝑚 𝑗 2 superscript 𝑄 2 \chi=x(1+(m_{i}+m_{j})^{2}/Q^{2}) , on an isoscalar nucleon target.

Energy	${\sigma }_{\mathrm{CC}}^{\mathrm{BDHM}}$	${\sigma }_{\mathrm{NC}}^{\mathrm{BDHM}}$	${\sigma }_{\mathrm{CC}}^{\mathrm{Soyez}}$	${\sigma }_{\mathrm{NC}}^{\mathrm{Soyez}}$
$\left[\mathrm{GeV}\right]$	[cm2]	[cm2]	[cm2]	[cm2]
${10}^7$	$1.71\times {10}^{-33}$	$6.93\times {10}^{-34}$	$1.42\times {10}^{-33}$	$6.02\times {10}^{-34}$
${10}^8$	$4.26\times {10}^{-33}$	$1.78\times {10}^{-33}$	$3.37\times {10}^{-33}$	$1.47\times {10}^{-33}$
${10}^9$	$8.88\times {10}^{-33}$	$3.78\times {10}^{-33}$	$7.19\times {10}^{-33}$	$3.18\times {10}^{-33}$
${10}^{10}$	$1.63\times {10}^{-32}$	$7.03\times {10}^{-33}$	$1.41\times {10}^{-32}$	$6.35\times {10}^{-33}$
${10}^{11}$	$2.72\times {10}^{-32}$	$1.19\times {10}^{-32}$	$2.60\times {10}^{-32}$	$1.17\times {10}^{-32}$
${10}^{12}$	$4.24\times {10}^{-32}$	$1.86\times {10}^{-32}$	$4.63\times {10}^{-32}$	$2.05\times {10}^{-32}$

Equations

$$\Biggl{\langle}\frac{dE}{dX}\Biggr{\rangle}=-(\alpha+\beta E)\ ,$$

$$X(D)={\int }_0^D\rho (\ell )d\ell$$

$$\beta ={\beta }^{\mathrm{brem}}+{\beta }^{\mathrm{pair}}+{\beta }^{\mathrm{nuc}},$$

$$\beta (E)=\frac{N}{A}\int dyy\frac{d\sigma (y,E)}{dy}$$

$$\ell +N(p)\to {\ell }^{\prime }+X,q=\ell -{\ell }^{\prime }.$$

$$\frac{d^2{\sigma }^{\mathrm{nuc}}}{dxd{Q}^2}$$

$$F_L$$

$$F_2(x,Q^2)$$

$$S(A,x,Q^2)=\frac{F_2^A(x,Q^2)}{A{F}_2(x,Q^2)},$$

$$S(A,x)=\left\{\begin{array}[]{cl}A^{-0.1}&x<0.0014\\ A^{0.069\log_{10}x+0.097}&0.0014<x<0.04\\ 1&0.04<x\end{array}\right.$$

$${\sigma }_{\mathrm{dip}}^A(x,r)$$

$${\sigma }_{\mathrm{dip}}^A(x,r,b)$$

$${\rho }_A(z,\vec{b})$$

$$T_A(b)$$

$$\int d^2\vec{b}{T}_A(b)$$

$${\beta }_{\tau }^{\mathrm{fit}}(E)={\beta }_0+{\beta }_1\ln (E/E_0),$$

$$\Biggl{\langle}\frac{dE}{dX}\Biggr{\rangle}\simeq\frac{dE}{dX}$$

$$R(E_{\tau }^i)=\int dX{P}_{\mathrm{surv}}(E_{\tau }^i,X).$$

$$E_{\tau }(E_{\tau }^i,x)$$

$$E_{\tau }^i(E_{\tau },x)$$

$$P_{\rm surv}(E_{\tau}^{i},z)=\exp\Biggl{[}-\int_{0}^{z}\frac{dx}{c\tau_{\tau}E_{\tau}(E_{\tau}^{i},x)/m_{\tau}}\Biggr{]}\ .$$

$$\frac{dN}{d{E}_{\tau }}=K(E_{\tau }^i{)}^{-2}{P}_{\mathrm{surv}}(E_{\tau }^i,E_{\tau })\frac{E_{\tau }^i}{E_{\tau }}{e}^{{\beta }_{1}X},$$

$$F_2^{CC}$$

$$F_2^{NC}$$

$$\sigma (E)=[{\sigma }_0+{\sigma }_1\ln (E)+{\sigma }_2{\ln }^2(E)]\cdot 10^{-31}\mathrm{c}{\mathrm{m}}^2,$$

$$\sigma (E)={\sigma }_0{E}^{{\sigma }_1}\cdot 10^{-35}\mathrm{c}{\mathrm{m}}^2,$$

$$⟨y(E)⟩=y_0+y_1\ln (E)+y_2{\ln }^2(E).$$

$$F(i)=\frac{dN}{d{E}_i}$$

$$F_{\tau }(E_{\tau })$$

$$D=2R_{\oplus }\cos \theta .$$

$$\alpha =\frac{\pi }{2}-\theta .$$

$$F_{\tau }(E_{\tau },D)=\int d{E}_{\nu }{F}_{\nu }(E_{\nu },0)K(E_{\nu },D,E_{\tau }).$$

$$(R_{\oplus }+h{)}^2=R_{\oplus }^2+x^2+2R_{\oplus }x\cos \theta .$$

Full text

I Introduction

Astrophysical sources that produce observed high energy cosmic rays and photons should also be sources of neutrinos (see, e.g., Refs. Gaisser et al. (1995); Katz and Spiering (2012); Anchordoqui et al. (2014)). While individual neutrino sources have not yet been detected, the IceCube Neutrino Observatory has evidence of a diffuse flux of neutrinos that lies above the background of neutrinos produced in the atmosphere for neutrino energies above $\sim 10^{5}-10^{6}$ GeV Aartsen et al. (2014). Further characterization of the diffuse flux is underway at IceCube and at underwater observatories Fusco (2016).

Detection of neutrinos from astrophysical sources is expected to reveal details about the site of cosmic ray acceleration and about the propagation of cosmic rays through the thermal cosmic background. The feature that neutrinos can escape sources without interacting and are not affected by magnetic fields translates to a challenge in neutrino detection. The underground observatories in ice and water instrument a large volume of target material that is transparent to the light signals that neutrino weak interactions produce. At higher energies, the neutrino interaction probabilities increase, however, the expected flux of neutrinos decreases more quickly as a function of energy, requiring even larger target volumes.

A well-known technique to overcome the problem of large target volumes for ultrahigh energy neutrinos has been explored for detecting astrophysical tau neutrinos in air shower experiments. The idea is to use the Earth as a neutrino converter then to use an emerging tau decay into a shower to detect the ultrahigh energy tau neutrino flux Domokos and Kovesi-Domokos (1997); Fargion (2004, 2002); Fargion et al. (2004); Athar et al. (2000); Feng et al. (2002); Kusenko and Weiler (2002); Hou and Huang (2002); Tseng et al. (2003); Aramo et al. (2005); Palomares-Ruiz et al. (2006). Currently, the Pierre Auger Observatory has set limits on the tau neutrino flux in the energy range $2\times 10^{8}$ GeV up to $2\times 10^{10}$ GeV Abraham et al. (2009); Abreu et al. (2012) using this technique. The potential to use the volcano adjacent to HAWC as the neutrino converter has also been explored Vargas et al. (2016). In addition, there is a proposal for the Ashra detector that would also use a mountain as a neutrino converter, for example Mauna Kea Asaoka and Sasaki (2013).

Among the uncertainties in these tau neutrino flux limits or eventual observation are three elements: the uncertainty in the high energy neutrino cross section Gandhi et al. (1998); Reno (2005); Jeong and Reno (2010); Goncalves and Hepp (2011); Cooper-Sarkar et al. (2011); Fiore and Zoller (2012); Goncalves and Gratieri (2015); Albacete et al. (2015); Arguelles et al. (2015) required to convert the tau neutrino to the tau lepton in the Earth, the inelasticity for neutrino scattering with isoscalar nucleons (how much of the neutrino energy is transferred to the tau) and the tau electromagnetic loss Dutta et al. (2001); Chirkin and Rhode (2004); Dutta et al. (2005); Armesto et al. (2008); Bigas et al. (2008); Koehne et al. (2013) as the tau transits the remaining rock before it emerges to decay into an air shower.

Tau energy loss is through electromagnetic interactions via ionization, bremsstrahlung, pair production and photonuclear interactions, e.g., see Refs. Dutta et al. (2001); Chirkin and Rhode (2004); Dutta et al. (2005); Armesto et al. (2008); Lipari and Stanev (1991); Lohmann et al. (1985); Antonioli et al. (1997); Koehne et al. (2013); Bigas et al. (2008); Sokalski et al. (2001). At high energies, the ionization energy loss is negligible. For taus, the photonuclear and pair production processes are much larger than bremsstrahlung, and as the energy increases, the photonuclear interaction dominates. For muons, bremsstrahlung, pair production and photonuclear interactions give comparable contributions to the electromagnetic energy loss Lohmann et al. (1985). The uncertainties in the high energy behavior of the photonuclear interaction depend on the parameterizations and assumptions about the nucleon structure. The neutrino cross section uncertainties are also related to the nucleon structure at higher momentum transfers.

In this paper, we focus on two goals. The first is to evaluate the theoretical uncertainties in the photonuclear energy loss, neutrino cross sections and neutrino inelasticity in a way that correlates the assumptions about the nuclear structure for both scattering processes. We consider energies between $10^{6}-10^{12}$ GeV, for taus and their corresponding neutrinos that come from cosmic sources and focus more on energies above $10^{8}$ GeV. We use an input spectrum of neutrinos scaling with energy as $E^{-2}$ as a representative neutrino flux. In addition, we consider a cosmogenic neutrino flux which is the flux from the interactions of cosmic rays with the cosmic neutrino background (GZK neutrinos). Cosmic rays have been detected up to $3\times 10^{20}$ eV$=3\times 10^{11}$ GeV Bird et al. (1994), and the shape of the cosmogenic neutrino flux is deduced from those observations.

Our second goal is to provide a semi-analytic approximation to the stochastic evaluation of first, the neutrino interaction, then the tau energy loss in matter. Uncertainties from a continuous versus stochastic treatment of the tau energy loss are discussed here (see also Bigas et al. (2008)). We examine the impact of stochastic versus approximate treatments of the neutrino interaction itself as well. The semi-analytic treatment reliably allows a survey of the impact on the tau flux of multiple approaches to calculating the neutrino (differential) cross section and tau energy loss.

We begin the paper with a summary of electromagnetic energy loss of charged leptons, focusing on the photonuclear process and the range of predictions that come from different approaches to modeling the electromagnetic structure function $F_{2}$ beyond the range of experimental measurements. To exhibit the effectiveness of analytic approximations to the tau energy loss we show an example for an incident tau lepton flux given by power-law, $E^{-2}$. Section 3 connects $F_{2}$ from electromagnetic scattering of charged leptons to the high energy neutrino cross section. In Section 4, we show results for incident neutrino fluxes characterized by a power law and for a representative cosmogenic neutrino flux. We discuss applications to the Pierre Auger Observatory in Section 5. Our conclusions are in Section 6.

II Electromagnetic energy loss of charged leptons

II.1 Average energy loss

The propagation of charged leptons through materials can be described through the average energy loss through electromagnetic processes. The average energy loss is usually written as

[TABLE]

in terms of the column depth

[TABLE]

for density $\rho$ and distance $D$. The parameter $\alpha$ describes ionization energy loss from the Bethe-Block formula Rossi (1952), weakly dependent on energy, and at high energies, weakly dependent on lepton mass $m_{\ell}$. For $E=10^{6}-10^{9}$ GeV, $\alpha=3.1-3.6\times 10^{-3}$ GeV cm2/g. At high energies, above a muon energy $E\sim 10^{3}$ GeV and tau energy $E\sim 10^{4}$ GeV, one can ignore the energy independent ionization energy loss, which we will do here.

The three contributions from electromagnetic energy loss to $\beta$ come from bremsstrahlung, electron-positron pair production and the photonuclear (inelastic scattering) processes. They are written as

[TABLE]

where $\beta$ is calculated for each process by

[TABLE]

in terms of the inelasticity parameter $y=(E-E^{\prime})/E$ for incident energy $E$ and outgoing energy $E^{\prime}$. As discussed below, $\beta(E)$ depends on the charged lepton, so we label it with $\beta_{\ell}$ for $\ell=\mu$ or $\ell=\tau$.

We evaluate the bremsstrahlung energy loss following Petrukhin and Sestakov Petrukhin and Shestakov (1968). Pair production energy loss follows Ref. Kokoulin and Petrukhin (1971). As shown in Ref. Koehne et al. (2013) (see also Ref. Klein (1999)), suppression of pair production and bremsstrahlung by the Landau-Pomeranchuk-Migdal and Ter-Mikaelian effects from coherence effects and dielectric suppression do not affect muon (and tau) energy loss below $E_{\ell}=10^{12}$ GeV, so we do not discuss it here. Scaling of the energy loss dependence on the charged lepton mass is discussed in Refs. Tannenbaum (1991); Reno et al. (2005). For bremsstrahlung, $\beta_{\ell}\sim 1/m_{\ell}^{2}$, while for pair production, $\beta_{\ell}\sim 1/m_{\ell}$.

The photonuclear process also scales as $\beta_{\ell}\sim 1/m_{\ell}$. The photon exchange between charged lepton and nucleon probes the nucleon structure. The differential cross section in Eq. (4) depends on the nucleon structure functions $F_{1}(x,Q^{2})$ and $F_{2}(x,Q^{2})$ which depend on the lepton four momentum transfer squared $Q^{2}=-(\ell-\ell^{\prime})^{2}=-q^{2}$ and $x=Q^{2}/(2p\cdot q)$ for

[TABLE]

The differential distribution for the photonuclear interaction is

[TABLE]

The full $x$ and $Q^{2}$ dependence of $F_{2}(x,Q^{2})$ in the photonuclear electromagnetic interaction gives an energy dependence to $\beta_{\ell}^{\rm nuc}(E)$ Dutta et al. (2001).

The dependence on the structure functions, in particular at small values of $x$ and low $Q^{2}$, introduces an uncertainty in $\beta_{\ell}^{\rm nuc}$. One option is to use a direct parameterization of $F_{2}(x,Q^{2})$, and to use the leading order Callan-Gross relation $2xF_{1}=F_{2}$ to obtain $F_{1}(x,Q^{2})$. Corrections can be implemented with $R=F_{L}/(2xF_{1})$. The dependence of the photonuclear energy loss on a reasonable range of $R$ was shown to be small in Ref. Koehne et al. (2013). In the cases of direct parameterizations, we set $R=0$. Using structure functions based on parton model evaluations is not possible here, even with small-$x$ extrapolations, because of the low $Q^{2}$ required.

Photonuclear interactions, in general, consist of the soft component (non-perturbative) and the hard component (perturbative). Early evaluations (see, e.g., Ref. Lohmann et al. (1985)) of $\beta_{\ell}^{\rm nuc}$ relied on a $Q$ independent parameterization by Bezrukov and Bugaev in Ref. Bezrukov and Bugaev (1981). Bezrukov and Bugaev modeled the soft component of the structure functions with a modified generalized vector dominance model. In Ref. Dutta et al. (2001), the impact of the full range of $Q^{2}$ in the evaluation of $\beta_{\ell}^{\rm nuc}$ was discussed using the structure functions parameterization of Abramowicz et al. (ALLM) Abramowicz et al. (1991); Abramowicz and Levy (1997). A similar parameterization of Butkevich and Mikhailov Butkevich and Mikheyev (2002) based on the Capella et al. model (CKMT) Capella et al. (1994) has been used in Ref. Chirkin and Rhode (2004); Koehne et al. (2013). The Bezrukov and Bugaev result was augmented by mass corrections and a hard component that depends on $Q^{2}$ in Ref. Bugaev and Shlepin (2003); *Bugaev:2003sw, giving similar results to the ALLM parameterization. More recently, a new parameterization of $F_{2}$ by Block et al. (BDHM) Block et al. (2014), guided by unitarity considerations Block et al. (2013a, b), has been provided. It has a somewhat different behavior at low-$x$ than the ALLM parameterization. We show below that the difference in small-$x$ behavior of the BDHM and ALLM parameterizations of $F_{2}$ have implications for the high energy flux of taus emerging from the Earth.

A QCD-motivated approach to the photonuclear electromagnetic interactions has the photon, given its wave function, split into a $q\bar{q}$-pair which then interacts with the nucleon non-perturbatively for the case when their transverse separation is large and perturbatively when it is small. This so-called dipole model Nikolaev and Zakharov (1991); Mueller (1994), while based on perturbation theory, can be used as an approximate calculation even at low $Q^{2}$. The structure function therefore depends on the photon wave function (squared) which itself depends on the transverse separation $r$ and longitudinal momentum fraction $z$, and the dipole cross section $\sigma_{\rm dip}$, implicitly summed over quark flavors:

[TABLE]

Formulae for the wave functions squared are listed in the appendix. In the dipole model, the longitudinal structure function can be included explicitly as well. Different models are used to provide the dipole cross section. We use the dipole cross section of Soyez Soyez (2007). The Soyez dipole cross section is based on a functional form discussed by Iancu, Itakura and Munier Iancu et al. (2004) from an approximate solution to the non-linear Balitsky-Kovchegov (BK) equations Balitsky (1996); Kovchegov (1999). Dipole cross section parameters were fit to $F_{2}$ data using Eq. (8) for $0.045$ GeV${}^{2}\leq Q^{2}\leq 150$ GeV2 and $10^{-6}\leq x\leq 10^{-2}$. Similar results are obtained with the dipole of Albacete et al. (AAMQS) Albacete et al. (2015). The AAMQS dipole includes further theoretical improvements to the approximate solution to the BK equations including the running coupling constant.

Each of the approaches gives a prediction for the small-$x$ electromagnetic structure function for $x$ below the measured regime. Fig. 1 shows four approaches for $Q^{2}=0.25$ GeV2 and $Q^{2}=1.5$ GeV2 , along with HERA data combined from H1 and ZEUS Aaron et al. (2010). The ALLM parameterization of $F_{2}$ has a steeper growth with small-$x$ than the other approaches. The BDHM parameterization, the Soyez dipole and AAMQS dipole do well in comparison with the small-$x$ HERA data and have similar small-$x$ behaviors. Since the Soyez and AAMQS dipoles have similar small-$x$ behavior, we use the Soyez dipole as representative of dipole models in our evaluation of the tau fluxes below.

It may be eventually possible to distinguish between the ALLM parameterizations and other approaches from LHC data for low invariant mass and forward rapidity kinematic regions. These data constrain low-$x$ parton distribution functions, as discussed in, e.g., Ref. Gauld and Rojo (2017). At least for sufficiently large $Q^{2}$, the parton distribution functions can be combined to form the electromagnetic structure functions at smaller $x$ values than accessible with HERA, where the ALLM curve is distinct from the other curves in Fig. 1.

For $\beta_{\ell}^{\rm nuc}$ for rock, with $Z=11$ and $A=22$, nuclear corrections are required. One can define

[TABLE]

the nuclear structure function normalized by $A$ times the proton structure function. For the ALLM and BDHM parameterizations, we use a shadow factor of the form $S(A,x,Q^{2})\simeq S(A,x)$,

[TABLE]

based on Fermilab E665 data Adams et al. (1992a, b, 1995). For the dipole model, one can use the Glauber-Gribov (GG) prescription to incorporate the dipole cross section for protons to the one for nuclei, where

[TABLE]

This requires the nuclear density $\rho_{A}$ and nuclear profile $T_{A}(b)$, normalized to unity, with

[TABLE]

and $a^{2}=2R_{A}^{2}/3$. The quantity $R_{A}$ is the nuclear radius $R_{A}=1.12\,A^{1/3}-0.86/A^{1/3}$ fm. Finally, for the Soyez dipole, nuclear effects can be included with a modification of the saturation scale as discussed by Armesto, Salgado and Wiedemann (ASW) Armesto et al. (2005).

Our focus is on tau energy loss, however, we show the muon $\beta_{\mu}(E)$ to show the differences between approaches. In Figs. 2 and 3, we show the photonuclear contributions to $\beta_{\ell}$ for muons and taus as a function of energy for rock. For comparison, we also show the Bugaev and Shlepin evaluation of $\beta_{\ell}^{\rm nuc}$ Bugaev and Shlepin (2003), an update of Bezrukov and Bugaev (1981) with both soft and hard contributions, and the pair production and bremsstrahlung (muons only) contributions to $\beta_{\ell}$. For muons, the BDHM and Soyez evaluations of $\beta_{\mu}^{\rm nuc}$ fall in a similar range, while the ALLM and Bugaev-Shlepin results for $\beta_{\mu}^{\rm nuc}$ rise more quickly with energy. For taus, again the ALLM paramterization and Bugaev-Shlepin result for $\beta_{\tau}^{\rm nuc}$ rises more steeply with energy than the dipole model and BDHM parameterization results. In our discussion below, we will use the ALLM parameterization as representative of the steeper growth of $\beta_{\tau}^{\rm nuc}$ with tau energy.

As noted above, for the dipole model, we can evaluate the nuclear correction with a shadow factor and Glauber-Gribov (GG) correction. The shadow function or GG correction reduces the unshadowed energy loss photonuclear $\beta_{\mu}^{\rm nuc}$ by a factor of $0.73-0.75$ for $E=10^{6}-10^{12}$ GeV. For the Soyez dipole model, we can also use the ASW prescription, which reduces $\beta_{\mu}^{\rm nuc}$ by a factor of about 10% relative to the uncorrected value, for the same energy range. We note that the Soyez-ASW $\beta_{\mu}^{\rm nuc}$ coincides with the AAMQS result using the GG correction or shadow factor. For taus, a similar nuclear suppression is obtained: a factor of $0.73-0.79$ for $\beta_{\tau}^{\rm nuc}$ with the shadow function or GG correction, and $0.90-0.94$ for the ASW correction for the same energy range for the Soyez dipole cross section.

For taus, energy loss effects can be largely determined by $\beta_{\tau}$ rather than a full Monte Carlo when considering incident tau neutrino fluxes Chirkin and Rhode (2004); Dutta et al. (2005); Koehne et al. (2013), as we discuss below. In Ref. Dutta et al. (2005), we found that a reasonable parameterization of $\beta_{\tau}$ is

[TABLE]

for $E_{0}=10^{10}$ GeV and initial tau $E=10^{8}-10^{12}$ GeV. The bremsstrahlung and pair contributions in this energy range are essentially constant, so $\beta_{0}=\beta_{0}^{\rm brem}+\beta_{0}^{\rm pair}+\beta_{0}^{\rm nuc}$ and $\beta_{1}=\beta_{1}^{\rm nuc}$. Parameterizations for $\beta_{\tau}^{\rm nuc}$ are listed in Table 1. In this energy range, the Soyez and BDHM results are such that the ratio of the numerical evaluation to the fit form gives $0.97<\beta_{\tau}^{\rm nuc}/\beta_{\tau}^{\rm fit,nuc}<1.06$. The fits are best above $E=10^{9}$ GeV, within $\pm 3\%$. We have also included in the Table the parameterization of $\beta_{\tau}^{\rm nuc}$ for the AAMQS dipole calculation of the average energy loss of the tau. The ALLM $\beta_{\tau}^{\rm nuc}$ is not well fit with two parameters. A three parameter fit to $\beta_{\tau}^{\rm nuc}$ is a better match to the calculated $\beta_{\tau}^{\rm nuc}$ for all the models, as discussed in Appendix B. For the ALLM evaluation of $\beta_{\tau}^{\rm nuc}$, the three parameter fit does well for $E=10^{8}-10^{11}$ GeV. We focus our discussion in the text on the two parameter fit. For the flux results below, the inclusion of $\beta_{2}$ in the fit does not impact flux predictions except for the ALLM case. We use the three parameter fit for the ALLM case.

II.2 Survival probability and average range

We start with discussion of the average lepton range to allow for comparisons with other authors. It has been known for a long time that approximating

[TABLE]

does not accurately represent the average muon range because of the effect of energy fluctuations in electromagnetic interactions Lipari and Stanev (1991). Muons are essentially stable at these energies, however, the tau’s lifetime is important here. Following the prescription described in Lipari and Stanev (1991) and others Antonioli et al. (1997); Sokalski et al. (2001); Chirkin and Rhode (2004); Bigas et al. (2008); Koehne et al. (2013), we use a Monte Carlo program to evaluate the tau survival probability. For $y<10^{-3}$, the energy loss is evaluated in the continuous approximation, but for $y>10^{-3}$, it is evaluated stochastically. The average range $R(E_{\tau}^{i})$, the average distance the tau travels, is determined by integrating the survival probability $P_{\rm surv}(E_{\tau}^{i},X)$ as a function of column depth for taus with an initial energy $E_{\tau}^{i}$ that survive to an energy larger than $E_{\tau}^{f}$,

[TABLE]

In the Monte Carlo evaluation of the average range, we use $E_{\tau}^{f}=50$ GeV for the tau, however, the short lifetime of the $\tau$ dominates its range below $E_{\tau}^{i}\sim 10^{8}$ GeV. Our Monte Carlo evaluation of the average tau range is shown in Fig. 4 and in Table 2. Our results, based on our Monte Carlo developed for Ref. Dutta et al. (2001) are in good agreement with Ref. Bigas et al. (2008). They are in qualitative agreement with those in Refs. Chirkin and Rhode (2004) and Koehne et al. (2013) where energies over many orders of magnitude are shown. For incident tau neutrino fluxes that fall with energy, the average range, without detailed final energy information, is less relevant than the survival probability as a function of final energy, as we discuss in Sec. IV.

In the approximations of Eqs. (15) and (16), just accounting for energy loss over a distance $x=X/\rho$, the initial tau energy $E_{\tau}^{i}$, final tau energy $E_{\tau}$ and distance traveled are related by integrating $dE/dX$ to get

[TABLE]

for the 2 parameter fit to $\beta_{\tau}(E)$. For our discussion of tau fluxes, it will be useful to also write $E_{\tau}^{i}$ in terms of the final $E_{\tau}$ after travelling distance $x$:

[TABLE]

Appendix B shows the relation between $E_{\tau}$, $E_{\tau}^{i}$ and $x$ for the three parameter fit.

With the lifetime, one can approximate the survival probability Reya and Rodiger (2005); Bigas et al. (2008),

[TABLE]

An explicit expression for the survival probability for the two parameter fit to $\beta_{\tau}(E)$, in the first order expansion in terms of $\beta_{1}/\beta_{0}$, is shown in Ref. Dutta et al. (2005). In this paper, we present the survival probability for the three parameter fit in Appendix B. Integration of Eq. (20) to get an approximation to the average tau range, as in the muon case, overestimates the average range. Nevertheless, as we discuss below, the flux of tau is well approximated by using the analytic formulas for the survival probability.

We remark that the tau weak interaction energy losses are small compared electromagnetic energy loss except at the very highest energies. Above approximately $E=10^{10}$ GeV, tau charged current interactions are more important than tau decays for eliminating taus from the flux, however, in this regime, electromagnetic energy loss dominates the tau interactions. Weak neutral current contributions to the tau energy loss are quite small Kuzmin et al. (2007). We include weak interaction effects in the Monte Carlo program for tau propagation. Without the tau charged current interaction, the tau range is between 1.3-8% (0.8-3.5%) higher for $E_{\tau}=10^{8}-10^{12}$ GeV for the evaluation of the energy loss using the BDHM (ALLM) parameterization.

II.3 Charged lepton spectra

Before turning to the full process of tau neutrino conversion to taus, followed by tau propagation to exit the Earth, we use a tau flux as a proxy for the neutrino flux, as was done in Ref. Bigas et al. (2008). This is a useful exercise to see the features of the tau propagation without the additional element of neutrino interactions, energy transfer and neutrino flux attenuation. Our results confirm the conclusions of Ref. Dutta et al. (2005); Bigas et al. (2008), namely, that except at high energy, evaluations of the propagation of taus using the analytic survival probability and Eq. (19) reliably reproduce the results of the full Monte Carlo propagation for an incident $(E_{\tau}^{i})^{-2}$ flux propagating over distances of 1–20 km.

In Fig. 5, we show the tau flux with an $(E_{\tau}^{i})^{-2}$ incident spectrum scaled by $E_{\tau}^{2}$ (solid line), and the resulting flux after propagation of 1 km, 2 km, 10 km and 20 km of rock with the BDHM parameterization of the photonuclear contribution (dotted and dashed lines). The Monte Carlo results are shown with dashed lines, which cut off at high energies compared with the analytic approximation shown with the dotted lines. The Monte Carlo generated an input spectrum for $E_{\tau}^{i}=10^{8}-10^{12}$ GeV. Our transport code was developed with look-up tables that extend to $10^{12}$ GeV, so the resulting tau flux from the Monte Carlo propagation shown here comes from an input flux with a sharp cutoff at $E_{\tau}^{i}=10^{12}$ GeV.

The dotted lines in the figure show the analytic formula for the propagated flux at column depth $X$:

[TABLE]

where $E_{\tau}^{i}$ comes from Eq. (19) in the 2 parameter model for $\beta_{\tau}$. Here, we show the analytic result for the $(E_{\tau}^{i})^{-2}$ tau flux without a cutoff. A sharp cutoff for $E_{\tau}^{i}$ translates directly to a sharp cutoff in $E_{\tau}$ for a fixed distance. The colored dots in Fig. 5 show $E_{\tau}^{\rm max}(E_{\tau}^{i\,{\rm max}},x)$ for $x=1,\ 2,\ 10$ and 20 km in rock for a sharp cutoff at $E_{\tau}^{i\,{\rm max}}=10^{12}$ GeV, the same cutoff used in the Monte Carlo. If one would impose the same sharp cutoff in analytic evaluation as was done in our MC, the tau flux would have a steep drop for energies above those labeled with the colored dots. For $E_{\tau}>E_{\tau}^{\rm max}$, $E_{\tau}^{2}dN/dE_{\tau}\to 0$ in the analytic evaluation with a cutoff for $E_{\tau}^{i}$.

Fig. 5 shows the role of the lifetime, which suppresses the lowest energy. With increasing depth, the flux becomes steeper. A priori, the agreement between the analytic and the Monte Carlo results is not obvious for energies between where the lifetime does not dominate and sufficiently lower than shown by the colored dots in Fig. 5. The Monte Carlo results include energy fluctuations, and as noted above, the average range evaluated from the Monte Carlo is shorter than what would be calculated by integrating the survival probability.

The agreement between the analytic and Monte Carlo evaluations comes from competing effects, including the fact that the flux decreases with energy. For a fixed incident tau energy and propagation distance, the Monte Carlo transport yields a slightly higher average final tau energy than the analytic final energy evaluated with Eq. (18). Fig. 5 shows the flux as a function of the final tau energy, so instead, the comparison should be of the initial tau energy. The average corresponding incident tau energy in the Monte Carlo approach is lower than the corresponding analytic value (see Eq. (19)). The survival probability as evaluated by the Monte Carlo approach is also lower than from the analytic formula in Eq. (20). Schematically, Eq. (21) shows that decreases in the incident tau energy and survival probability will tend to compensate. Our numerical comparison shows that this is the case for the $(E_{\tau}^{i})^{-2}$ incident flux. Thus comparison of the tau flux (i.e. final tau energy) in the two approaches effectively corresponds to different values of the $E_{\tau}^{i}$ and somewhat different survival probabilities.

To evaluate the tau flux due to incident neutrinos passing through a depth of rock, we need the neutrino cross section, and the differential cross section as a function of $y$. In the next section we show results for neutrino cross sections as evaluated by different models, then use the Monte Carlo to propagate the neutrinos and produced taus through rock.

III Neutrino cross sections

The flux of taus coming from tau neutrino conversions in the Earth depends on the neutrino cross section, both for neutrino attenuation and for neutrino conversion to taus, and on the differential distribution as a function of inelasticity $y$. As shown in Section II, different models for the electromagnetic structure function $F_{2}$ have different low-$x$ predictions. The neutrino cross sections at high energies probe the low-$x$ behavior of the weak structure functions, albeit at a higher value of $Q^{2}$, namely $Q^{2}\sim M_{W}^{2}$ (see, e.g., Refs. Gandhi et al. (1998); Reno (2005)).

For a fixed column depth, considering both the probability for neutrino interactions to produce the tau lepton, and the tau lepton energy loss, there are competing effects due to the low-$x$ behavior of the structure functions. For the ALLM parameterization, the neutrino cross section is enhanced relative to the BDHM parameterization result, as shown in Fig. 6 with the dotted and solid lines. However, there is also a greater energy loss of the tau at high energies because of the same small-$x$ enhancement. For the results below, we keep the same model for both the neutrino interaction and the tau interaction.

Detailed formulae for the neutrino weak interactions are shown in Appendix C. For the dipole model, the wave functions for $W$ or $Z$ fluctuations to a quark-antiquark pair are functions of the quark masses Fiore and Zoller (2012); Barone et al. (1993); Kutak and Kwiecinski (2003); Arguelles et al. (2015). We include them in Appendix C, along with the relevant couplings. Given, e.g., the Soyez dipole cross section, the ultrahigh energy neutrino cross section is predicted, as shown in Fig. 6 with the (green) dashed curve.

For $E>10^{7}$ GeV, valence contributions, neglected here, are small. Fig. 6 shows that these three approaches yield nearly the same neutrino-nucleon cross sections at $E_{\nu}=10^{6}-10^{7}$ GeV. As noted in Ref. Armesto et al. (2008), the ALLM electromagnetic structure function lies somewhat below the HERA data for $Q^{2}\sim 200-8000$ GeV2 and $x>10^{-2}$. We note that measurements of the neutral current structure function at $Q\sim M_{W}$ are available only for $x>0.1$ Aaron et al. (2010), the scale relevant to the ultrahigh energy neutrino cross section Gandhi et al. (1998). The ALLM parameterization yields neutrino cross sections that are too low below $E_{\nu}\sim 10^{5}$ GeV. Since we focus here on energies even higher than $10^{6}$ GeV, we are probing very small-$x$ values where the ALLM parameterization is useful as a contrast to other approaches with its stronger dependence on $x$.

The explicit inclusion of quark channels in the dipole model through the wave function squared cannot be done with the electromagnetic structure function parameterizations of ALLM and BDHM. Instead, we rescale the electromagnetic structure function by the sum of electric charges squared for five quark flavors, to yield

[TABLE]

We use the Callan-Gross relation $2xF_{1}=F_{2}$ and set $F_{3}=0$ for neutrino interactions when we use the BDHM and ALLM parameterizations of $F_{2}$. The charged current cross sections for these models are also shown in Fig. 6, along with the next-to-leading order perturbative evaluation of the neutrino charged-current cross section from Ref. Jeong and Reno (2010). At $E_{\nu}=10^{6}$ GeV, the cross sections are within $\pm 5\%$ of each other, increasing to $\pm 12\%$ for $E_{\nu}=10^{8}$ GeV. As the the neutrino energy further increases, the spread in cross sections increase to a factor of $\sim 1.7$ for the BDHM and Soyez approach compared to the ALLM extrapolation at $10^{12}$ GeV. For reference, recent evaluations of the neutrino cross sections with a broader range of models Cooper-Sarkar et al. (2011); Fiore and Zoller (2012); Goncalves and Gratieri (2015); Albacete et al. (2015) show a similar spread in cross sections at these energies.

The tau flux also requires a knowledge of the energy transferred from the neutrino to the tau, namely, the distribution as a function of $y\equiv(E_{\nu}-E_{\tau})/E_{\nu}$. The average $y$ as a function of energy is shown in Fig. 7. We note that the average $y$ varies from $-3\%$ to $7\%$ at $10^{6}$ GeV and $-15\%$ to $23\%$ at $10^{11}$ GeV with respect to the values from the Soyez approach. This range of the average $y$ represents theoretical uncertainty in the elasticity due to different approaches that we consider. In the Monte Carlo evaluation, we use look-up tables for the neutrino cross section and distribution in $y$ up to a maximum $E_{\nu}=10^{12}$ GeV, however, we also do a semi-analytic evaluation using $\langle y\rangle$. We show in the next section that $\langle y\rangle$ and a parameterization of the neutrino cross section works well to reproduce the Monte Carlo results. We include here a parameterization of the neutrino charged current cross section and average $y$. For BDHM and Soyez, we find

[TABLE]

while for ALLM and a NLO perturbative evaluation (JR in Ref. Jeong and Reno (2010)),

[TABLE]

for $E$ in GeV. The average $y$ can be written as

[TABLE]

The parameters $\sigma_{i}$ and $y_{i}$ are shown in Tables 3 and 4. The cross section parameterizations reproduce the charged current cross sections to within $\pm 5\%$. The parameterizations reproduce $\langle y\rangle$ to within $\pm 1\%$.

In evaluating the tau neutrino induced tau flux, we also need the neutrino neutral current cross section for the attenuation factor. For the range of $E_{\nu}=10^{8}-10^{12}$ GeV, the ratio of the neutral current to charged current cross section ranges between 0.42-0.44 for the BDHM parameterization, 0.44-0.45 for the Soyez dipole, and 0.41-0.42 for the ALLM parameterization.

IV Results

We begin this section by comparing a Monte Carlo evaluation of neutrino interactions followed by tau propagation through rock of various depths. To quantify the effects of the neutrino cross section to produce taus, attenuation of the neutrino flux and tau energy loss, we plot the ratio of the emerging tau flux to the incoming neutrino flux – the transmission function $F(\tau)/F(\nu)$ – where

[TABLE]

is the number of leptons per unit energy per unit area per sec, assuming an isotropic incident flux. We only show results for incident neutrinos. For $E_{\nu}$ larger than $\sim 10^{7}$ GeV, $\sigma_{\nu N}\simeq\sigma_{\bar{\nu}N}$ to a very good approximation.

In the next section, we consider distances less than 10 km, relevant to using mountains as neutrino converters. In Section IV.B, we evaluate the transmission functions for larger distances where neutrino attenuation effects come into play.

IV.1 Distances less than 10 km

In Fig. 8, we present the results from the analytic evaluation and Monte Carlo evaluation for the depths of 2, 5, and 10 km of the standard rock ($\rho=2.65\ {\rm g/cm^{3}}$) for all models considered here. The incident neutrino flux is assumed to have $E_{\nu}^{-2}$ spectrum. The analytic calculation uses the survival probability for the tau lepton (Eq. (B)) for the three parameter fit to $\beta_{\tau}(E)$ written in Eq. (35). The tau flux, for depth D, is

[TABLE]

where $E_{\nu}=E_{\tau}^{i}/(1-y)$. While evaluating the analytic fluxes, we extended the neutrino cross sections for the energies above $10^{12}$ GeV, while for the Monte Carlo results, the neutrino flux is cut off at $10^{12}$ GeV, leading to tau fluxes that start to steeply decrease at energies that depend on the depth, as shown in Fig. 8.

Monte Carlo results are well matched to the analytic calculations for energies below the maximum tau energy in the analytic approach. There are less than $10\%$ errors for $E\leq 10^{11}$ GeV for the BDHM and the Soyez dipoles, where the impact of the neutrino energy cutoff at $10^{12}$ GeV starts to become evident in the Monte Carlo evaluation. For the ALLM parameterization, the similar difference appears between the Monte Carlo and analytic results at $10^{10}$ GeV. The lower energy of the onset of the cutoff for the ALLM parameterization comes from the higher level of electromagnetic energy loss in this model. Since the Monte Carlo results are so well represented by the analytic approximation below $E\sim 10^{10}-10^{11}$ GeV depending on the depth, we do the rest of our analysis using the analytic tau survival probability integrated with the neutrino differential cross section and the attenuation factor. By comparing the Monte Carlo results (dashed curves) with the analytic curves with no neutrino energy cutoff (solid lines), Fig. 8 clearly shows a feed-down effect of the neutrino energy cutoff at $E_{\nu}=10^{12}$ GeV. At these short distances, the cutoff effect will scale with the cutoff energy. In the remaining figures in Sec. IV.A, we use the analytic approximation without a cutoff for the evaluation of the transmission function.

In the previous section, we presented the parameterizations of $\langle y\rangle$ and the neutrino charged current cross sections. For reference, we compare the resulting tau fluxes using $\langle y(E)\rangle$ and the parameterized cross section with the numerical evaluation of the neutrino differential cross sections. We evaluate $\langle y(E)\rangle$ using $E=E_{\tau}$. The tau fluxes with the parameterizations deviate from the fully evaluated fluxes by at most $4\%$ using Eq. (26) for $\langle y\rangle$ and by at most $\pm 5\%$ using both the parameterized $\langle y(E)\rangle$ and $\sigma_{\nu N}(E)$ between $E=10^{7}-10^{12}\ {\rm GeV}$.

One of the features of the analysis presented here is that the neutrino cross section evaluations use the same form for the small-$x$ behavior as the evaluation of the electromagnetic energy loss. At fixed depth, one can see evidence of this effect in Fig. 9. The solid lines in the figure show the ratio of the outgoing tau flux to the incoming tau neutrino flux when the parameterization (BDHM for upper curves and ALLM for lower curves) is used for both effects. The dashed curves have the BDHM or ALLM parameterizations used for the electromagnetic energy loss, but where the neutrino cross section is evaluated using NLO perturbative QCD Jeong and Reno (2010). The impact of the perturbative QCD choice is significant in the flux ratio at high energies when combined with the BDHM parameterization, since the small-$x$ behaviors are so different. The ALLM electromagnetic energy loss with the perturbative neutrino cross section shows less variation with a complete ALLM evaluation because the ALLM structure function and parton distribution function small-$x$ extrapolations are similar. We discuss the impact of using the perturbative QCD neutrino nucleon cross section below, however, for the remainder of this section, we keep consistent the approach to the small-$x$ behavior for both tau energy loss and neutrino nucleon scattering.

Fig. 10 compares the flux ratio from the different approaches for the various depths. The HAWC Vargas et al. (2016) and Ashra Asaoka and Sasaki (2013) configurations are sensitive to Earth skimming tau neutrinos at lower energies than the Pierre Auger Observatory. For example, the tau neutrino energy range for the proposed Ashra experiment is PeV–EeV. Therefore, we extend our results to lower energies, i.e., $10^{6}$ GeV, for those experiments. In the low energy range, below $10^{8}$ GeV, the tau decay process dominates, and the fluxes are hardly affected by the energy loss. Therefore, the difference of the fluxes in this energy range reflects mostly the ratio of the neutrino cross sections from the different models.

As presented in Fig. 10, the BDHM and the Soyez dipoles give the comparable results all over the energy range. There is at most a 25% difference, which is its maximum at $E=10^{8}$ GeV. The wider error band from all the different models depends on the behavior of the fluxes from the ALLM model compared to the BDHM/Soyez approaches. For the Soyez dipole, we only show the results with nuclear shadowing included by a multiplicative factor, the same factor also used with the BDHM and ALLM models. It is almost the same as the predicted flux ratio with GG nuclear effects. The ASW nuclear corrections for the Soyez flux evaluation yield results at most 10 % higher.

As shown in Fig. 3, the energy loss from the ALLM parameterization increases more quickly with energy than the other models. This reduces the tau flux at high energies with a higher rate, and its effect appears more significantly at larger depth. Thus, the uncertainty range is the largest at $10^{12}$ GeV for the $10$ km depth, where the flux from the Soyez dipole is the largest. It is $83\%$ higher than that from the ALLM parameterization. For shorter depths, tau energy loss is less important, so the higher neutrino cross section from the ALLM parameterization shown in Fig. 6 results in the higher flux ratio for the ALLM evaluation.

In Fig. 11, we present the tau fluxes produced from a representative incident cosmogenic neutrino flux. The figure is in the same format as Fig. 10. We take the cosmogenic neutrino flux from the pure proton cosmic ray model with the AGN evolution scenario shown in Ref. Aloisio et al. (2015). Other assumptions about the cosmic ray composition and evolution scenarios primarily change the overall normalization of the neutrino flux, so $F(\tau)/F(\nu)$ is representative for all of the cosmogenic neutrino fluxes in Ref. Aloisio et al. (2015). The tau fluxes from the BDHM and Soyez models have the similar relations to the case of the $E^{-2}$ flux for incident neutrinos. However, the ALLM transmission function is relatively lower, and it gives rise to the larger uncertainty for the depths of 5 km and 10 km, with the largest error is still about a factor of 2. In this case, the Soyez dipole approach yields a factor of 2.1 larger flux than the ALLM approach at the highest energies.

Figs. 8-11 show the neutrino propagation through short distances relative to the tau range. The tau range goes up to $\sim 50$ km for tau propagation in rock for $E_{\tau}^{i}<10^{12}$ GeV, as shown in Fig. 4. Requiring the emerging tau to have at least 10% of the incident tau energy reduces the tau range to $\sim 10$ km at $E_{\tau}^{i}=10^{12}$ GeV. The tau lifetime sets another distance scale, namely, $\gamma c\tau=4.9\ {\rm km}\cdot(E/{10^{8}\ {\rm GeV}})$.

For distances less than 10 km in rock, the increase in distance in rock translates to an increase in number of targets for neutrino interactions. At a distance of 10 km of rock, the transmission function for an incident $E_{\nu}^{-2}$ flux is a factor of between 2.6-3.2 larger than for 2 km of rock for the BDHM parameters. Tau energy loss and decay, depending on the energy, account for the discrepancy between the transmission function and a factor of 5 coming from the ratio of the two depths, 2 km and 10 km. Neutrino attenuation is not significant, even at the highest energies, since the interaction length in rock is $\sim 1,500-150$ km for $E_{\nu}=10^{8}-10^{12}$ GeV.

Fig. 12 shows the overall uncertainties on the ratio of the emerging tau flux to the incident neutrino flux from the different approaches for the neutrino cross section, the tau energy loss and the neutrino inelasticity, where a single approach is not used consistently for both the weak and electromagnetic structure functions. The solid curves are calculated in consistent approaches, hence the yellow and the orange bands present the uncertainties by the approaches as in Fig. 10 and Fig. 11. We also investigated how big is the difference on the results due to the neutrino cross section and the inelasticity from the different approaches using their different combinations. In the figures, we present the results that give the largest differences, with the dashed curves for the respective models. The largest contribution to the uncertainties due to the evaluation in the mixed approaches is from the neutrino cross section. The range of neutrino cross sections results in 30% (50%) contribution towards the transmission function uncertainty at $E=10^{8}$ GeV ($E=10^{11}$ GeV), while the mean neutrino interaction inelasticity gives about 5% and 10% difference at the same energies. The overall uncertainties are about 30% (60%) at $10^{8}$ GeV and a factor of 3.3 (3.8) at $10^{11}$ GeV for the diffuse $E^{-2}$ flux (cosmogenic flux).

IV.2 Larger distances

For depths larger than 10 km, one begins to see a saturation in the transmission function, then attenuation effect for even larger depths. In this section, we use the $E_{\nu}^{-2}$ flux to demonstrate these effects.

In the upper plot of Fig. 13, we show with the dashed lines the Monte Carlo result for 30 km (upper), 50 km (middle) and 100 km (lower) for $E_{\nu}^{-2}$ with $E_{\nu}^{\rm max}=10^{12}$ GeV using the BDHM approach. For an energy $\sim 10^{8}$ GeV, the transmission function is nearly equal for 30, 50 and 100 km of rock. For this energy, the transmission function is dominated by the decay length $\lambda_{\tau}(E)$. For much larger distances, the produced tau will decay inside the rock. The relevant target column depth for the incident neutrino is $\sim\rho\lambda_{\tau}(E)$ which grows linearly with energy. As the tau energy and decay lengths increase, tau energy loss becomes more important. Then the relevant neutrino target column depth is $\sim R_{\tau}(E)\rho$. The transmission function for $30-50$ km at $E\sim 10^{9}-10^{10}$ GeV is nearly the same, growing with energy because the neutrino cross section increases with energy faster than the energy loss effects degrade the tau energy. For 100 km of rock, neutrino attenuation comes in as well, so that the transmission function is $\sim 10-20\%$ reduced relative to 30-50 km. At higher energies, the curves for transmission function begin to further separate as more neutrinos interact before the final column depth of $\sim R_{\tau}(E)\rho$.

The dotted lines in the upper plot of Fig. 13 show the evaluation of the transmission function where the Monte Carlo program is used only for the neutrino interaction and the maximum neutrino energy is $E_{\nu}=10^{12}$ GeV. The tau energy $E_{\tau}^{i}$ is generated stochastically, smearing the initial tau energy and therefore the maximum $E_{\tau}$ for the analytic survival probability calculation. For these three distances, the dotted lines match the Monte Carlo results to within 10%. This hybrid method of Monte Carlo simulation for the neutrino interaction and analytic result for the tau propagation speeds up the numerical comparisons for longer distances, where most of the taus do not emerge after longer distances.

The lower plot in Fig. 13 shows the transmission function for three different approaches to the energy loss: BDHM, Soyez and ALLM. Increased tau energy loss and some attenuation come into play in the ALLM evaluation compared to the BDHM and Soyez evaluations.

The Monte Carlo results are in good agreement with the results obtained analytically, within about 5% and 10% at $E=10^{10}$ GeV for the BDHM and Soyez approaches, respectively. Above $10^{10}$ GeV, the Monte Carlo results start decreasing due to the cutoff of $E_{\nu}^{max}=10^{12}$ GeV: the difference between Monte Carlo and analytic evaluations with no cutoff, at $E=10^{11}$ GeV, is about 20% for BDHM and about 35% for Soyez.

Fig. 14 shows the transmission function calculated analytically for the incident cosmogenic neutrino flux. For the 100 km, the neutrino attenuation results in the decrease of the transmission function at the high energies relative to 30 km.

Fig. 15 shows the overall uncertainties on the ratio of the emerging tau flux to an incident neutrino flux as in Fig. 12, but for the larger distance of 100 km. The band shows the impact of including the range of neutrino cross sections and neutrino inelasticity obtained in all three approaches for the tau energy loss obtained with the Soyez approach. The overall uncertainties in this figure, for the energy range of $10^{8}$ GeV$\leq E\leq 10^{11}$ GeV, are within 35% for both $E_{\nu}^{-2}$ and cosmogenic neutrino fluxes. The largest difference is at energy of about $10^{10}$ GeV, about 32% and 34% for $E_{\nu}^{-2}$ and cosmogenic neutrino fluxes. At $E=10^{11}$ GeV, the uncertainty reduces to 24% and 30%, respectively. From Fig. 13 we note that for $E_{\nu}^{-2}$ neutrino flux and distance of 100 km rock, the different approaches to tau energy loss give larger uncertainty than the uncertainty due to neutrino cross section or inelasticity. For $E=10^{8}$ GeV, the transmission functions differ by a factor of 1.3, while for $E=10^{10}$ GeV, the BDHM gives a factor of 2.0 larger result that the one evaluated in ALLM approach.

A distance of 100 km of rock corresponds to an angle of approximately $0.5^{\circ}$ below the horizon. Fig. 16 shows schematically the chord length $D$ as a function of nadir angle $\theta$, where

[TABLE]

For a nadir angle of $85^{\circ}$, the chord length is more than 1,000 km. For neutrinos with energies larger than a few times $10^{8}$ GeV, the neutrino interaction length is smaller than the chord length for nadir angles much less than $85^{\circ}$, so the transmission function will be appreciable only for nadir angles close to $90^{\circ}$. We define the angle $\alpha$ relative to the horizon:

[TABLE]

For reference, we show the chord lengths and maximum depth of the neutrino/tau trajectory below the Earth’s surface (sagitta) in Table 5.

In Fig. 17, we show a plot of $F(\tau)/F(\nu)$ for different angles $\alpha$ relative to the horizontal using $R_{\oplus}=6.37\times 10^{3}$ km. The surface rock-mantle interface is several hundred kilometers below the Earth’s surface Palomares-Ruiz et al. (2006), so for $\alpha\leq 6^{\circ}$, $\rho=2.65$ g/cm3 is a good approximation for a detector on land. Each curve represents the transmission function for a given approach and the shaded band shows the spread in predictions for $F(\tau)/F(\nu)$ for tau neutrinos skimming with $\alpha=1^{\circ},\ 2^{\circ},\ 4^{\circ}$ and $6^{\circ}$ below the horizon. Similar to the shorter trajectories, the span of predictions at $E=10^{8}$ GeV is less than at $E=10^{10}$ GeV and higher. At $E=10^{8}$ GeV, $F(\tau)/F(\nu)$ from the BDHM approach is a factor of $1.2-1.4$ higher than for the ALLM calculation for all of the angles shown. By $E=10^{10}$ GeV, the BDHM transmission function is a factor of $2.2\ (1^{\circ})$ to $13\ (6^{\circ})$ higher.

How the uncertainty in the transmission function translates to a particular measurement depends on the relative importance of different energy regions and angles. If one considers air showers that must develop before an altitude $h=20$ km and taking the time dilated decay length as the distance the tau travels after emerging from the Earth at an angle $\alpha=1^{\circ}-6^{\circ}$, then $E_{\tau}<8.3\times 10^{9}$ GeV. For neutrino fluxes that decrease with energy, there is a further suppression of the high energy regime.

Before we turn to the specific example of the Pierre Auger Observatory, we discuss the impact of different spectral indices on the range of predictions for the transmission function. In Fig. 18, we present the transmission function for the incident neutrino fluxes scaling as $E_{\nu}^{-2.5}$ and $E_{\nu}^{-2}$ taking the Soyez approach for the reference. The transmission function for the steeper neutrino flux is smaller. For example, the $E_{\nu}^{-2.5}$ incident neutrino flux gives transmission function about 23% (33%) lower than in case of incident neutrino flux being $E_{\nu}^{-2}$ at $E=10^{8}$ GeV for the 10 km (100 km) distance. This reduction effect becomes larger as energy increases, at $E=10^{11}$ GeV it is about 30% (40%). We find effect to be the same for the BDHM approach, while for the ALLM approach it is 27% (45%) at the corresponding energies for the 10 km distance.

V Application: Pierre Auger Observatory

Our results are directly applicable to an approximate evaluation of the effective aperture of the Pierre Auger Observatory. Following Refs. Aramo et al. (2005); Miele et al. (2006); Abraham et al. (2009), we use the analytic survival probability for the kernel $K(E_{\nu},D,E_{\tau})$, involving an integration of Eq. (27) to get

[TABLE]

To simplify the expression, we use $\langle y\rangle=y_{0}=0.15$ for all energies. The $\delta$-functions relating initial and final tau energies and the neutrino energy to the initial tau energy simplify the expression of the kernel. The kernel for the 2-parameter fit of $\beta$ appears in Appendix B.

The aperture depends on the tau energy and lifetime. Ref. Abraham et al. (2009) suggests that the relevant altitude for detection is about $x_{0}=10$ km along the shower axis from the point of the tau decay. The kernel is translated to the effective aperture (approximately) by first evaluating the differential probability to have a tau decay a distance $x-x_{0}$ from exiting from the Earth, as a function of altitude $h$ above the Earth’s surface at $x$. The configuration is shown schematically in Fig. 19.

The distance $x$ and altitude $h$ are related to the nadir angle $\theta$, as is the chord length $D$ that first the tau neutrino, then the tau, travels through the Earth. We take a constant Earth density $\rho=2.65$ g/cm3 for the angles of relevance, namely $\theta=\pi/2-\alpha\to\pi/2$ for $\alpha=0.1$ rad. The distances $x$ and $h$, and nadir angle $\theta$, are related by

[TABLE]

The probability density is

[TABLE]

where $\lambda_{dec}=E_{\tau}c\tau_{\tau}/m_{\tau}$ is the usual decay length.

We evaluate the effective aperture by integrating over the fraction of the solid angle that is important, over the tau energies and over altitudes weighted by an approximate efficiency for detection, a product of the identification and trigger efficiencies. The maximum efficiency for identification and trigger is each at 0.84 since the muonic decay of the tau does not produce a shower. In Ref. Abraham et al. (2009), efficiencies for selected energies are shown and described as being primarily functions of altitude and tau energy. We approximate the effective efficiency by an energy dependent step function for $h$,

[TABLE]

This gives the effective aperture in the analytic approximation,

[TABLE]

The evaluation of the effective apertures for the Pierre Auger Observatory with $A=3,000$ km2 are shown with the solid lines for Soyez, BDHM and ALLM approaches from top to bottom in Fig. 20. Our results here, with the crude approximation for the efficiency, are within a factor of $\sim 4$ of those of Ref. Miele et al. (2006), where a careful accounting of the local topography near the Observatory was included. Our approximations serve to test the impact of the different approaches to energy loss and neutrino cross sections. In the range of $E_{\nu}=10^{9}-10^{12}$ GeV, the Soyez and BDHM evaluations differ by between $2-7\%$. At $E=10^{9}$ GeV, the ALLM curve is about $6\%$ lower, however ALLM aperture drops with higher energies, so that the Soyez aperture is a factor of 1.7 times larger than the ALLM aperture.

The dashed lines in the figure show the impact of including a parameterized NLO perturbative neutrino nucleon cross section from Ref. Jeong and Reno (2010). For $E_{\nu}=10^{9}-10^{12}$ GeV, the ratios of Soyez to ALLM apertures are from 1.04 to 2.0. For the event rates for an $E^{-2}$ incident tau neutrino flux, the peak of the differential event rate is at $E_{\tau}$ a few times $10^{8}$ GeV, where differences in the approaches are small. Our approximate calculations have the total event rates differing by only a few percent.

VI Conclusions

We have calculated the emerging tau flux generated from tau neutrino interactions followed by the tau propagation through rock. We have evaluated neutrino charged current cross section and the tau energy loss in various QCD-motivated approaches, such as dipole approach as well as the more phenomenological approach that determines neutrino cross section using the same $F_{2}$ parameterization used in evaluating the tau energy loss parameter $\beta_{\tau}^{\rm nuc}$. Our comparisons of results obtained analytically with those obtained with the Monte Carlo treatment of the tau energy loss show good agreement for fluxes that scale like $E_{\nu}^{-2}$. While not a fundamental result, the analytic treatment of tau propagation in materials is a useful tool to survey different approaches to calculate tau electromagnetic energy loss. We find that the tau electromagnetic energy loss at PeV range can be well described by $\beta_{\rm fit}(E)=\beta_{0}+\beta_{1}\ln(E/E_{0})+\beta_{2}\ln^{2}(E/E_{0})$ with the parameters in Tables 1 and 6. For the BDHM and dipole model treatment of electromagnetic energy loss, $\beta_{2}\simeq 0$ is a good approximation.

We have also parameterized our results for the neutrino charged current cross section for each approach we considered, reproducing the numerical result to within $\pm 5\%$. Finally comparisons of the $\nu_{\tau}$ fluxes after propagation over various distances in rock, between Monte Carlo simulations and analytic calculations for each theoretical approach, are in good agreement in the energy regime where the Monte Carlo upper limit on the neutrino energy has no impact. The parameterization of the survival probability and energy relations are useful tools to determine the implications of new models or approaches without a full Monte Carlo simulation.

The different approaches to the tau energy loss calculation have larger impacts at high energies and longer trajectories in the Earth. For example, different approaches to electromagnetic energy loss can give transmission functions that differ by a factor of $\sim 10$ for $D\sim 1,000$ km. Nevertheless, our results for the Pierre Auger Observatory aperture, with a simplified efficiency, show that event rates with an $E_{\nu}^{-2}$ incident neutrino flux are minimally uncertain due to the energy loss and the neutrino cross section, because at $E_{\tau}\simeq 10^{9}$ GeV, the approaches considered here give very similar results.

In addition to applying our results to evaluating the effecting aperture for the Pierre Auger Observatory, our results can be applied to the High Altitude Water Cherenkov (HAWC) observatory and the All-sky Survey High Resolution Air-shower detector (ASHRA). For HAWC, the propagation distance is 0.35 to 6 km (width of the volcano at summit and at the observatory elevation) Vargas et al. (2016). We have compared results using the Monte Carlo simulation to analytic calculation of the tau flux exiting the Earth performed using the approaches presented in this paper for the above distances. Additionally, one can also consider the effective aperture for HAWC, given the approximate distance of 4 km from the edge of the volcano to the HAWC array and the array surface area of $22,000$ m2. Since HAWC is sensitive to showers that pass through detector elements coming from neutrino conversions in the volcano, the effective solid angle is much smaller than for the Pierre Auger Observatory. Similar conditions can be found in the Ashra configuration proposed to be 30 km from Mauna Kea, with a propagation distance of approximately 20 to 25 km. One can adopt our Auger aperture calculation with $\beta_{\tau}$ parameterizations along with the analytic computations for the tau neutrino propagation to predict the induced events at the detector arrays. The calculational technique discussed here could also be applied to a space-based telescope pointing towards the Earth as discussed in, e.g., Ref. Neronov et al. (2017).

By using a consistent low-$x$ treatment of both energy loss and the cross section, we have shown that the different parameterizations of $F_{2}$ and dipole approach yield similar emerging tau fluxes for short distances and for tau energies $\sim 10^{8}-10^{9}$ GeV. By using a separate NLO parton model calculation of the neutrino cross section, the range of emerging tau flux predictions is broader for the approaches considered here. Overall, for large distances, the differences in modeling the tau energy loss account for most of the differences in the calculated transmission functions. As forward, low invariant mass processes are studied in more detail, small-$x$ parton distribution functions will help constrain $F_{2}$ at $Q^{2}$ values where the parton model applies. The ALLM parameterization gives a larger tau energy loss at high energies than the dipole model or BDHM parameterization. If the ALLM parameterization or the BDHM parameterization were ruled out experimentally, the theoretical uncertainty associated with ultrahigh energy skimming tau neutrinos would be reduced.

Acknowledgements.

We thank J. Alvarez-Muniz and G. Parente for discussions. This research was supported in part by the US Department of Energy contracts DE-SC-0010113, DE-FG02-13ER-41976 and DE-SC0009913.

Appendix A Dipole model

In the dipole model, the electromagnetic structure function $F_{2}(x,Q^{2})$ is written in terms of the photon wave functions and the dipole cross section Nikolaev and Zakharov (1991)

[TABLE]

The photon wave functions, in terms of the transverse separation $r$ and the fractional momentum $z$, are given by

[TABLE]

using the modified Bessel functions $K_{0}$ and $K_{1}$, the quark electric charge $e_{f}$ and mass $m_{f}$. Here, $\bar{Q}_{f}^{2}=z(1-z)Q^{2}+m_{f}^{2}$, and we use $m_{u}=m_{d}=m_{s}=0.14$ GeV, $m_{c}=1.4$ GeV and $m_{b}=4.5$ GeV. The electromagnetic fine-structure constant is labeled $\alpha_{em}$ and $N_{c}=3$ is the number of colors in Eqs. (33-34). For the dipole cross section, $\sigma_{\rm dip}(r,x)$, we use the parameterization in Ref. Soyez (2007) by Soyez.

Appendix B Three parameter fit to $\beta_{\tau}(E)$

In order to extend our investigation to lower energies, a better parameterizations for the energy loss $\beta_{\tau}$ than Eq. (15) is useful. Adding a $\ln^{2}E$ term, we use

[TABLE]

with $E_{0}=10^{10}$ GeV. Results for the parameters are in Table 6. The results for the BDHM and Soyez approaches are matched well with the output from the direct calculation. The differences are mostly less than 1% between $E=10^{6}-10^{12}$ GeV. The ALLM fit is not as good as the other models, especially at low energies, but it works well, to within 5 % between $E=2\times 10^{7}-5\times 10^{8}$ GeV. Between $E=5\times 10^{8}-10^{11}$ GeV, it works better, to within 3 %.

Starting from the energy dependence of $\beta_{\tau}$, the relation between the initial tau energy $E_{\tau}^{i}$ and distance traveled $z^{\prime}$ to the final tau energy $E_{\tau}$ is

[TABLE]

with ${\cal B}=\sqrt{\beta_{0}\beta_{2}-\beta_{1}^{2}}$. The survival probability for the three parameter fit to $\beta_{\tau}$ is

[TABLE]

where the exponential integral ${\rm Ei}(x)=-\int^{\infty}_{x}e^{-t}/t\ dt$.

Setting $\beta_{2}=0$ and $d\sigma_{CC}(E_{\nu})/dE_{\tau}=\sigma_{CC}(E_{\nu})\delta(E_{\tau}-(1-y_{0})E_{\nu})$, the expression for the kernel for the effective aperture is

[TABLE]

where

[TABLE]

Appendix C Weak interactions

It is convenient to write the differential neutrino-nucleon cross section as

[TABLE]

where $F_{T}^{V}=(F_{L}^{V}+F_{R}^{V})/2$, and $2xF_{3}^{V}=F_{L}^{V}-F_{R}^{V}$. The relative coupling of the $W$ and $Z$ are absorbed into the definition of the structure functions for the $Z$. For weak interactions, $\lambda=L,R,S$ with

[TABLE]

For the dipole model Fiore and Zoller (2012); Barone et al. (1993); Kutak and Kwiecinski (2003); Arguelles et al. (2015), where $\bar{z}=1-z$ and $\overline{Q}=z\bar{z}Q^{2}+m_{1}^{2}z+m_{2}^{2}\bar{z}$,

[TABLE]