Simulation of Ultra High Energy Neutrino Interactions in Ice and Water

TL;DR

This paper adapts the CORSIKA simulation program to water and ice mediums to model ultra high energy neutrino interactions, enabling improved prediction of acoustic signals for neutrino detection.

Contribution

It introduces a novel adaptation of CORSIKA for water and ice, allowing detailed simulation of neutrino-induced showers and acoustic signals in these media.

Findings

Simulated hadronic showers in water and ice.

Detailed spatial energy distribution data.

Enhanced acoustic signal prediction accuracy.

Abstract

The CORSIKA program, usually used to simulate extensive cosmic ray air showers, has been adapted to work in a water or ice medium. The adapted CORSIKA code was used to simulate hadronic showers produced by neutrino interactions. The simulated showers have been used to study the spatial distribution of the deposited energy in the showers. This allows a more precise determination of the acoustic signals produced by ultra high energy neutrinos than has been possible previously. The properties of the acoustic signals generated by such showers are described.

Figures

Equations

$$\frac{d^{2}\sigma}{dQ^{2}dy}=\frac{G_{F}^{2}}{2\pi y}\bigg{(}\frac{M_{W}^{2}}{Q^{2}+M_{W}^{2}}\bigg{)}^{2}(F_{2}(x,Q^{2})(1-y+y^{2}/2)\pm y(1-y/2)xF_{3}(x,Q^{2}))$$

$$\frac{d^{2}E}{drdz}=L(z,E_L)\cdot R(r,z,E_L)$$

$$L(z,E_{L})=P_{1L}\bigg{(}\frac{z-P_{2L}}{P_{3L}-P_{2L}}\bigg{)}^{\frac{(P_{3L}-P_{2L})}{P_{4L}+P_{5L}z+P_{6L}z^{2}}}\exp\bigg{(}{\frac{P_{3L}-z}{P_{4L}+P_{5L}z+P_{6L}z^{2}}}\bigg{)}$$

$$\frac{P_{1L}}{E}=2.760\cdot 10^{-3}-1.974\cdot 10^{-4}{E}_L+7.450\cdot 10^{-6}{E}_L^2$$

$$P_{2L}=-210.9-6.968\cdot 10^{-3}{E}_L+0.1551E_L^2$$

$$P_{3L}=-41.50+113.9E_L-4.103E_L^2$$

$$P_{4L}=8.012+11.44E_L-0.5434E_L^2$$

$$P_{5L}=0.7999\cdot 10^{-5}-0.004843E_L+0.0002552E_L^2$$

$$P_{6L}=4.563\cdot 10^{-5}-3.504\cdot 10^{-6}{E}_L+1,315\cdot 10^{-7}{E}_L^2.$$

$$R(r,z,E_{L})=\frac{1}{I}\bigg{(}\big{(}\frac{r}{P_{1R}}\big{)}^{(P_{2R}-1)}\big{(}1+\frac{r}{P_{1R}}\big{)}^{(P_{2R}-4.5)}\bigg{)}$$

$$I=\int_{0}^{\infty}\bigg{(}\big{(}\frac{r}{P_{1R}}\big{)}^{(P_{2R}-1)}\big{(}1+\frac{r}{P_{1R}}\big{)}^{(P_{2R}-4.5)}\bigg{)}~{}dr=P_{1R}\frac{\Gamma(4.5-2P_{2R})\Gamma(P_{2R})}{\Gamma(4.5-P_{2R})}.$$

$$P_{nR}=A+Bz+C{z}^2$$

$$A=0.01287E_L^2-0.2573E_L+0.9636$$

$$B=-0.4697\cdot 10^{-4}{E}_L^2+0.0008072E_L+0.0005404$$

$$C=0.7344\cdot 10^{-7}{E}_L^2-1.375\cdot 10^{-6}{E}_L+4.488\cdot 10^{-6}$$

$$A=-0.8905\cdot 10^{-3}{E}_L^2+0.007727E_L+1.969$$

$$B=0.1173\cdot 10^{-4}{E}_L^2-0.0001782E_L-5.093\cdot 10^{-6}$$

$$C=-0.1058\cdot 10^{-7}{E}_L^2+0.1524\cdot 10^{-6}{E}_L-0.1069\cdot 10^{-8}.$$

$$\chi^{2}=\sum_{i}\bigg{(}\frac{F_{i}-D_{i}}{F_{i}+D_{i}}\bigg{)}^{2}$$

$$\frac{d^{2}E}{drdz}=Ek(\frac{z}{z_{max}}{)}^t\exp (t-z/\lambda )2\pi r\rho (r)$$

$$\rho (r)=\frac{1}{r_M^2}{a}^{s-2}(1+a{)}^{s-4.5}\frac{\Gamma (4.5-s)}{2\pi \Gamma (s)\Gamma (4.5-2s)}$$

$$\frac{d^{2}E}{drdz}=rf(z)g(r,z)$$

$$f(z)=\frac{E}{X_0}b\frac{(b{z}^{\prime }{)}^{a-1}\exp -b{z}^{\prime }}{\Gamma (a)}$$

$$z_{max}^{\prime }=0.65\log (\frac{E}{E_c})+3.93$$

$$g(r,z)=g_{0}\bigg{(}\frac{r_{i}}{r}\bigg{)}^{n}$$

$$P(d,t)=\frac{E\beta}{4\pi C_{p}}\int\frac{\epsilon}{d}\frac{d}{dt}\bigg{(}\delta(t-d/c)\bigg{)}dV$$

$$P(\omega )=\frac{1}{r}{\int }_{-\infty }^{\infty }\frac{E\beta }{4\pi C_p}\frac{d}{dt}{E}_{xyz}(t)e^{-i\omega t}dt=\frac{1}{r}\frac{E\beta }{4\pi C_p}i\omega {\int }_{-\infty }^{\infty }{E}_{xyz}(t)e^{-i\omega t}dt=\frac{1}{r}\frac{E\beta }{4\pi C_p}i\omega E_{xyz}(\omega )$$

$$P(t)=\frac{1}{1024}\sum _{n=-512}^{n=511}P({\omega }_n)A({\omega }_n)e^{in\Omega }$$

Full text

**Simulation of Ultra High Energy Neutrino Interactions in Ice and Water

(the ACoRNE Collaboration)a

**

S. Bevan1, S. Danaher2, J. Perkin3, S. Ralph3†, C. Rhodes4, L. Thompson3, T. Sloan5b and D. Waters1.

1* Physics and Astronomy Dept, University College London, UK.

2* School of Computing Engineering and Information Sciences, University of Northumbria, Newcastle-upon-Tyne, UK.

3* Dept of Physics and Astronomy, University of Sheffield, UK.

4* Institute for Mathematical Sciences, Imperial College London, UK.

5* Department of Physics, University of Lancaster, Lancaster, UK

†* Deceased

a* Work supported by the UK Particle Physics and Astronomy Research Council and by the Ministry of Defence Joint Grants Scheme

b Author for correspondence, email t.sloan@lancaster.ac.uk

The CORSIKA program, usually used to simulate extensive cosmic ray air showers, has been adapted to work in a water or ice medium. The adapted CORSIKA code was used to simulate hadronic showers produced by neutrino interactions. The simulated showers have been used to study the spatial distribution of the deposited energy in the showers. This allows a more precise determination of the acoustic signals produced by ultra high energy neutrinos than has been possible previously. The properties of the acoustic signals generated by such showers are described.

(Submitted to Astroparticle Physics)

1 Introduction

In recent years interest has grown in the detection of very high energy cosmic ray neutrinos [1]. Such particles could be produced in the cosmic particle accelerators which make the charged primaries or they could be produced by the interactions of the primaries with the Cosmic Microwave Background, the so called GZK effect [2]. The flux of neutrinos expected from these two sources has been calculated [3, 4]. It is found to be very low so that large targets are needed for a measurable detection rate. It is interesting to measure this neutrino flux to see if it is compatible with the values expected from these sources, incompatibility implying new physics.

Searches for cosmic ray neutrinos are ongoing in AMANDA [5], IceCube [6], ANTARES [7] and NESTOR [8], detecting upward going muons from the Cherenkov light in either ice or water. In general, these experiments are sensitive to lower energies than discussed here since the Earth becomes opaque to neutrinos at very high energies. The experiments could detect almost horizontal higher energy neutrinos but have limited target volume due to the attenuation of the light signal in the ice. The Pierre Auger Observatory, an extended air shower array detector, will also search for upward and almost horizontal showers from neutrino interactions [9]. In addition to these detectors there are ongoing experiments to detect the neutrino interactions by either radio or acoustic emissions from the resulting particle showers [1]. These latter techniques, with much longer attenuation lengths, allow very large target volumes utilising either large ice fields or dry salt domes for radio or ice fields, salt domes and the oceans for the acoustic technique.

In order to assess the feasibility of each technique the production of the particle shower from neutrino interactions needs to be simulated. Since experimental data on the interactions of such high energy particles do not exist it is necessary to use theoretical models to simulate them. The most extensive ultra high energy simulation program which has so far been developed is CORSIKA [10]. However, this program has been used previously only for the simulation of cosmic ray air showers. The program is readily available [10].

Different simulations are necessary for the radio and acoustic techniques. Radio emission occurs due to coherent Cherenkov radiation from the particles in the shower, the Askaryan Effect [11]. The emitted energy is sensitive to the distribution of the electron-positron asymmetry which develops in the shower and which grows for lower energy electromagnetic particles. Hence, to simulate radio emission, the electromagnetic component of the shower must be followed down to very low kinetic energies ($\sim 100$ keV)[12]. In contrast, an acoustic signal is generated by the sudden local heating of the surrounding medium induced by the particle shower [13]. Thus to simulate the acoustic signal the spatial distribution of the deposited energy is needed. Once the electromagnetic energy in the shower reaches the MeV level (electron range $\sim 1$ cm) the energy can be simply added to the total deposited energy and the simulation of such particles discontinued. Extensive simulations have been carried out for the radio technique [14]. However, the simulations for the acoustic technique are less advanced. Some work has been done [15, 16] using the Geant4 package [17]. However, this work is restricted to energies less than $10^{5}$ GeV for hadron showers since the range of validity of the physics models in this package does not extend to higher energy hadrons.

In this paper the energy distributions of showers produced by neutrino interactions in sea water at energies up to $10^{12}$ GeV are discussed. The distributions are generated using the air shower program CORSIKA [10] modified to work in a sea water medium. The salt component of the sea water has a negligible effect111The shower maximum was observed to peak at a depth $2.4\pm 1.1\%$ less in sea water than in fresh water with the same peak energy deposited, for protons of energy $10^{5}$ GeV. and the results are presented in distance units of g cm*-2*, hence they should be applicable to ice also. The computed distributions have been parameterised and this parameterisation is used to develop a simple program to simulate neutrino interactions and the resulting particle showers. The properties of the acoustic signals from the generated showers are also presented.

2 Adaptation of the CORSIKA program to a water medium

The air shower program, CORSIKA (version 6204) [10], has been adapted to run in sea water i.e. a medium of constant density of 1.025 g cm*-3* rather than the variable density needed for an air atmosphere. Sea water was assumed to consist of a medium in which $66.2\%$ of the atoms are hydrogen, $33.1\%$ of the atoms are oxygen and $0.7\%$ of the atoms are made of common salt, NaCl. The salt was assumed to be a material with atomic weight and atomic number A=29.2 and Z=14, the mean of sodium and chlorine. The purpose of this is to maintain the structure of the program as closely as possible to the air shower version which had two principal atmospheric components (oxygen and nitrogen) with a trace of argon. The presence of the salt component had an almost undetectable effect on the behaviour of the showers.

Other changes made to the program to accommodate the water medium include a modification of the stopping power formula to allow for the density effect in water 222The stopping power was computed using the Bethe-Bloch formula [20] and the density effect from the formulae of Sternheimer et al [21].. This only affects the energy loss for hadrons since the stopping powers for electrons are part of the EGS [18] package which is used by CORSIKA to simulate the propagation of the electromagnetic component of the shower. Smaller radial binning of the shower was also required since shower radii in water are much smaller than those in air. In addition the initial state energy for electrons and photons above which the LPM effect [19] was simulated in the program was reduced to the much lower value necessary for water 333The level was set at 1 TeV compared to the characteristic energy for water $E_{LPM}=270$ TeV [20].. The LPM effect suppresses pair production from high energy photons and bremsstrahlung from high energy electrons. Similarly, the interactions of neutral pions had to be simulated at lower energy than in air because of the higher density water medium. In all about 100 detailed changes needed to be made to the CORSIKA program to accommodate the water medium.

To test the implementation of the LPM effect [19] in the program 100 showers from incident gamma ray photons at several different energies were generated and the mean depth of the first interaction (the mean free path) calculated. The observed mean free path was found to be in agreement with the expected behaviour when both the suppression of pair production and photonuclear interactions were taken into account (see Figure 1). This showed that the LPM effect had been properly implemented in CORSIKA.

Considerable fluctuations between showers occurred. These are expressed in terms of the ratio of the root mean square (RMS) deviation of a given parameter to its mean value: the RMS peak energy deposit to the mean peak energy deposit was observed to be $14\%$ at $10^{5}$ GeV reducing to $4\%$ at $10^{11}$ GeV, that for the depth of the peak position varied from $19\%$ to $7.4\%$ and for the full width at half maximum of the shower from $63\%$ to $18\%$. To smooth out such fluctuations averages of 100 generated showers will be taken in the following. The statistical error on the averages is then given by these RMS values divided by 10. The hadronic energy contributes only about $10\%$ to the shower energy at the shower peak, the remainder being carried by the electromagnetic part of the shower.

The simulations were all carried out in a vertical column of sea water 20 m long. The deposited energy generated by CORSIKA was binned into 20 g cm*-2* slices longitudinally and 1.025 g cm*-2* annular cylinders radially for $0<r<10.25$ g cm*-2* and 10.25 g cm*-2* for $10.25<r<112.75$ g cm*-2* where $r$ is the distance from the vertical axis. To reduce computing times, the thinning option was used i.e. below a certain fraction of the primary energy (in this case $10^{-4}$) only one of the particles emerging from the interaction is followed and an appropriated weight is given to it [22]. The simulation of particles continued down to cut-off energies of 3 MeV for electromagnetic particles and 0.3 GeV for hadrons. When a particle reached this cut-off, the energy was added to the slice where this occurred. The QGSJET [23] model was used to simulate the hadronic interactions.

3 Comparison with other simulations

3.1 Comparison with Geant4

Proton showers were generated in sea water using the program Geant4 (version 8.0) [17] and compared with those generated in CORSIKA. Unfortunately, the range of validity of Geant4 physics models for hadronic interactions does not extend beyond an energy of $10^{5}$ GeV. Hence the comparison is restricted to energies below this.

Figure 2 shows the longitudinal distributions of proton showers at energies of $10^{4}$ and $10^{5}$ GeV (averaged over 100 showers) as determined from Geant4 and CORSIKA. The showers from CORSIKA tend to be slightly broader and with a smaller peak energy than those generated by Geant4. The difference in the peak height is $\sim 5\%$ at $10^{4}$ GeV rising to $\sim 10\%$ at energy $10^{5}$ GeV. Figure 3 shows the radial distributions. The differences in the longitudinal distributions are reflected in the radial distributions. However, the shapes of the radial distributions are very similar between Geant4 and CORSIKA, with CORSIKA producing $\sim 10\%$ more energy near the shower axis at depths between 450 and 850 g cm*-2* where most of the energy is deposited. The acoustic signal from a shower is most sensitive to the radial distribution, particularly near the axis ($r\sim 0$). It is relatively insensitive to the shape of the longitudinal distribution.

3.2 Comparison with the simulation of Alvarez-Muñiz

and Zas

The CORSIKA simulation was also compared with the longitudinal shower profile for protons computed in the simulation by Alvarez-Muñiz and Zas (AZ) [24]. There was a reasonable agreement between the longitudinal shower shapes from CORSIKA and those shown in Figure 2 of ref. [24]. However, the numbers of electrons and positrons at the peak of the CORSIKA showers was $\sim 20\%$ lower than those from ref. [24]. This number is sensitive to the energy below which these particles are counted and this is not specified in [24]. Hence the agreement between CORSIKA and their simulation is probably satisfactory within this uncertainty.

In conclusion, the modifications made to CORSIKA to simulate high energy showers in a water medium give results which are compatible with the predictions from the Geant4 simulations for energy less than $10^{5}$ GeV and the simulation of AZ within $20\%$. This is taken to be the accuracy of the simulation program assuming that there are no unexpected and unknown interactions between the centre of mass energy explored at current accelerators and those studied in these simulations. Studies of the sensitivity of the CORSIKA simulation to the different models of the hadronic interactions have been reported in reference [25]. They find that the peak number of electrons plus positrons varies by $\sim 20\%$ for proton showers in air depending on the choice of the hadron interaction model used. These differences are similar in magnitude to the differences between the AZ, Geant4 and CORSIKA simulations reported here. Hence the observed differences between the Geant4, AZ and CORSIKA simulations in water could be within the uncertainties of the hadronic interaction models.

4 Simulation of neutrino induced showers

Neutrinos interact with the nuclei of the detection medium by either the exchange of a charged vector boson ($W^{+}$), i.e. charged current (CC) interactions or the exchange of the neutral vector boson ($Z^{0}$), i.e. neutral current (NC) deep inelastic scattering interactions (see for example [26]). The ratio of the CC to NC interaction cross sections is approximately 2:1. The CC interactions produce charged secondary scattered leptons while the NC interactions produce neutrinos. The hadron shower carries a fraction $y$ of the energy of the incident neutrino and the scattered lepton the remaining fraction $1-y$. We assume that the neutrino flavours are homogeneously mixed when they arrive at the Earth by neutrino oscillations. Hence in the CC interactions electrons, $\mu$ and $\tau$ leptons will be produced as the scattered leptons in equal proportions. At the energies we shall consider, these particles behave in a manner similar to minimum ionising particles for $\mu$ and $\tau$ leptons. This is almost true also for electrons for which the bremsstrahlung process will be suppressed by the LPM effect. Hence the charged scattered leptons contribute little to the energy producing an acoustic signal. In the case of NC interactions there is no contribution to this energy from the scattered lepton. For these reasons the contribution of the scattered lepton to the shower profile is ignored beyond $z=20$ m in what follows.

It is interesting to note that a $\tau$ lepton can decay to hadrons or a very high energy electron or muon can produce bremsstrahlung photons at large distances from the interaction point. These can initiate further distant showers, the so called “double bang” effect. The stochastic nature of such electron showers is studied in [15, 16]. These effects are not considered in this study.

4.1 Neutrino-nucleon interaction cross sections.

A number of groups have computed the high energy neutrino-nucleon interaction cross sections, $\sigma$, [29, 27, 28]. In the quark parton model of the nucleon for the single vector boson exchange process, the differential cross section for CC interactions can be expressed in terms of the measured structure functions of the target nucleon $F_{2}$ and $xF_{3}$ as

[TABLE]

where $G_{F}$ is the Fermi weak coupling, $M_{W}$ is the mass of the weak vector boson, $Q^{2}$ is the square of the four momentum transferred to the target nucleon, $y=\nu/E$ where $\nu$ is the energy transferred to the nucleon ($\nu=E-E^{\prime}$ with $E$ and $E^{\prime}$ the energies of the incident and scattered leptons) and $x=Q^{2}/2M\nu$ is the fraction of the momentum of the target nucleon carried by the struck quark (here $x$ and $y$ are defined for a stationary target nucleon). The plus (minus) sign is for neutrino (anti-neutrino) interactions. It can be seen that $y$ is the fraction of the neutrino’s energy which is converted into the energy of the hadron shower. A similar expression can be written down for the NC interaction (see for example [26]) which has a ratio to the CC cross section varying from 0.33 to 0.41 as the neutrino energy increases from $10^{4}$ to $10^{13}$ GeV. The structure functions $F_{2}$ and $xF_{3}$ are the sum of the quark distribution functions which have been parameterised by fitting data [30, 31]. It can be shown that $Q^{2}=sxy$ where $s=2ME$ is the square of the centre of mass energy ($M$ is the target nucleon mass). To compute the cross sections the structure functions must be calculated at values of $x\lesssim M_{W}^{2}/s$ i.e. at values well outside the region of the fits to the parton distribution functions (PDFs) which have been performed for $x\gtrsim 10^{-5}$, the range of current measurements. The extrapolation outside the measurement range is discussed in [27], [29] and [32, 33]. Here we adopt the procedure of extrapolating linearly on a log-log scale from the parameterised parton distribution functions of [30] computed at $x=10^{-4}$ and $x=10^{-5}$. By considering various theoretical evolution procedures it is estimated in [29] that the procedure has an accuracy of $\sim 32\%$ per decade and we use this as an estimate of the accuracy of the calculation. However, this could be an underestimate [34].

The expression in equation 1 for charged current interactions and the one for neutral current interactions were integrated to obtain the total neutrino-nucleon interaction cross section, the value of the fraction of events per interval of $y$, $1/\sigma d\sigma/dy$, and the mean value of $y$. The total cross section was found to be in good agreement with the values in [27, 29] and in reasonable agreement with [28] which is based on a model different from the quark parton model. Figure 4 shows the mean value of $y$ obtained from this procedure (solid curve) and the effect of multiplying or dividing the PDFs by a factor 1.32 per decade (dashed curves) as an indication of the possible range of uncertainties in the extrapolation of the PDFs. Figure 5 shows the $y$ dependence of the cross section for different neutrino energies.

4.2 A simple generator for neutrino interactions.

A simple generator for neutrino interactions in a column of water of thickness 20 m was constructed as follows. The neutrino interacts at the top of the water column (z=0, with the z axis along the axis of the column). The energy fraction transferred, $y$, for the interaction was generated, distributed according to the curve for the energy of the neutrino shown in Figure 5. This allows the energy of the hadron shower to be calculated for the event. The assumption was made that these hadron showers will have approximately the same distributions as those of a proton interaction at z=0 (see Section 4.3 for a test of this assumption). A series of files of 100 such proton interactions were generated at energies in steps of half an order of magnitude between $10^{5}$ and $10^{12}$ GeV. The hadron shower for each neutrino interaction was selected at random from the 100 showers in the file at the proton energy closest to the energy of the hadron shower. The deposited energy in each bin was then multiplied by the ratio of the energy of the hadron shower to that of the proton shower. This is made possible because the shower shapes vary slowly with shower energy. For example, the ratio of the peak energy deposit per 20 g cm*-2* slice to the shower energy varies from 0.037 to 0.030 as the proton shower energy varies from $10^{5}$ to $10^{12}$ GeV.

4.3 The HERWIG neutrino generator.

The CORSIKA program has an option to simulate the interactions of neutrinos at a fixed point [35]. The first interaction is generated by the HERWIG package [36]. This option was adapted to our version of CORSIKA in sea water. Some problems were encountered with the $y$ dependence of the resulting interactions due to the extrapolation of the PDFs to very small $x$ at high energies. This only affects the rate of the production of the showers at different $y$ and the distribution of the hadrons produced in the interaction at a given $y$ should be unaffected.

A total of 700 neutrino interactions were generated at an incident neutrino energy of $2\cdot 10^{11}$ GeV. These were divided into the shower energy intervals $0.5-2\cdot 10^{10}$, $2-4\cdot 10^{10}$, $4-7.5\cdot 10^{10}$, $0.75-1.3\cdot 10^{11}$ and $1.3-2\cdot 10^{11}$. The showers in which the scattered lepton energy disagreed with the shower energy by more than $20\%$ were eliminated leading to a loss of $17\%$ of the events with shower energy greater than $0.5\cdot 10^{10}$ GeV. This is due to radiative effects and misidentification of the scattered lepton. Approximately 70 events remained in each energy interval. The energy depositions from these were averaged and compared to the averages from proton showers scaled by the ratio of the shower energy to the proton energy. Figure 6 shows the longitudinal distributions of the hadronic shower energy deposited for the different energy intervals (labelled $E_{W}$) compared to the scaled proton distributions. Figure 7 shows a sample of the transverse distributions.

There is a good consistency between the proton and neutrino induced showers. The proton showers peak, on average, 20 g cm*-2* shallower in depth with a peak energy $2\%$ larger than the neutrino induced showers. This is small compared to the overall uncertainty. The slight shift in the longitudinal distribution is reflected as a normalisation shift in the radial distributions. We conclude therefore that to equate a proton induced shower starting at the neutrino interaction point to that from a neutrino is a satisfactory approximation.

5 Parameterisation of showers

In this section a parameterisation of the energy deposited by the showers generated by CORSIKA (averaged over 100 showers depositing the same total energy) is described. Other available parameterisations will then be compared with the showers generated by CORSIKA.

The acoustic signal generated by a hadron shower depends mainly on the energy deposited in the inner core of the shower. This is illustrated in figure 8 which shows the contribution to the acoustic signal from cores of different radii. This figure shows that it is crucial to represent the deposited energy well at radius less than 2.05 g cm*-2*. The calculation of the acoustic signal from the deposited energy is described in section 6.

5.1 Parameterisation of the CORSIKA Showers

The differential energy deposited was parameterised as follows

[TABLE]

where the function $L(z,E_{L})$ represents the longitudinal distribution of deposited energy and $R(r,z,E_{L})$ the radial distribution. Here $E_{L}$ is $\log_{10}E$ with $E$ the total shower energy.

The function $L(z,E_{L})=dE/dz$ is a modified444 The modification is to replace the shape parameter $\lambda$ in equation 3.5 of reference [37] by the quadratic expression in $z$ in equation 3. version of the Gaisser-Hillas function [37]. This function represents the longitudinal distribution of the energy deposited.

[TABLE]

Here the parameters $P_{nL}$ were fitted to quadratic functions of $E_{L}=\log_{10}E$ with values

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The parameter $P_{1L}$ represents the peak energy deposited and $P_{3L}$ the depth in the $z$ coordinate at this peak while $P_{2L}$, $P_{4L}$, $P_{5L}$ and $P_{6L}$ are related to the shower width and shape in $z$.

The radial distribution was represented by the NKG function [37]

[TABLE]

where the integral

[TABLE]

The parameter $P_{1R}$ was found to vary strongly with depth while $P_{2R}$ was only a weak function of depth. The parameters $P_{nR}$ (with $n=1$,$2$) were each represented by the quadratic form

[TABLE]

and the quantities $A,B,C$ parameterised as quadratic functions of $E_{L}$. This gave for $P_{1R}$

[TABLE]

[TABLE]

[TABLE]

and for the parameter $P_{2R}$

[TABLE]

[TABLE]

[TABLE]

The fit was made in a depth range where $dE/dz$ was greater than $10\%$ of the peak value defined by equation 4. The program MINUIT [38] was used to minimise the squared fractional deviations

[TABLE]

where $F_{i}$ and $D_{i}$ refer to the fitted value and the value observed in the $i$th bin from the CORSIKA showers, respectively. In order to improve the fit at small radii the contributions to $\chi^{2}$ were arbitrarily weighted by 10 for $r<2.05$ g cm*-2*, 4 for $2.05<r<3.075$ g cm*-2*, unity for $3.075<r<51.25$ g cm*-2* and 0.25 for $r>51.25$ g cm*-2*. The RMS value of the fractional deviations was $3.4\%$ for radii less than 51.25 g cm*-2* and for energies greater than $10^{6.5}$ GeV. The fit becomes poorer at lower energies and greater radii than these. Integrating the parameterisation shows that the fraction of the total energy computed from the fit within the fit range was $91\%$ averaged over the deposited energy range $10^{7}$ to $10^{12}$ GeV. The corresponding fraction directly from the CORSIKA distributions was $92.5\%$, averaged over the same energy range. When applying this parameterisation at depths with smaller energy deposit than $10\%$ of the peak value, the energy was assumed to be confined to an annular radius of 1.025 g cm*-2*. There was a good agreement (within $5\%$ at the peak) between the acoustic signal computed using this parameterisation and that taken directly from the CORSIKA showers.

5.2 The parameterisation used by the SAUND Collaboration

The SAUND Collaboration [39] uses the following parameterisation [40], based on the NKG formulae (e.g. see reference [37]), for the energy deposited per unit depth, $z$, and per unit annular thickness at radius $r$ from a shower of energy $E$

[TABLE]

where $z_{max}=0.9X_{0}\ln(E/E_{c})$ is the maximum shower depth, $X_{0}=36.1$ g cm*-2* is the radiation length and $E_{c}=0.0838$ GeV. The constants $t=z_{max}/\lambda$ where $\lambda=130-5\log_{10}(E/10^{4}\mathrm{GeV})$ g cm*-2* and $k=t^{t-1}/\exp{(t)}\lambda\Gamma(t)$. The radial density is given by

[TABLE]

where $a=r/r_{M}$ with $r_{M}=9.04$ g cm*-2*, the Molière radius in water, and $s=1.25$. Figure 9 shows the radial distributions from CORSIKA compared with the absolute predictions of this parameterisation.

There is qualitative agreement between the parameterisation and the CORSIKA results. The difference in normalisation is explained by the somewhat different longitudinal profiles of the CORSIKA showers from the SAUND parameterisation. The latter are broader with a lower peak energy deposit and a depth of the maximum which is larger than the CORSIKA showers. CORSIKA predicts more energy at small $r$ than the SAUND parameterisation. Quantitatively, $51\%$ of the shower energy is contained within a cylinder of radius 4 cm for the CORSIKA showers compared to $35\%$ from the SAUND parameterisation. These fractions are approximately independent of energy. Hence, in acoustic detectors a harder frequency spectrum for the acoustic signals is predicted by CORSIKA than by the SAUND parameterisation. Note that in the fit described in Section 5.1 the values of the parameter $P_{1R}$ (equivalent to $R_{M}$ in equation 20) were strongly depth dependent and much lower than the Molière radius in water, assumed by the SAUND collaboration. In addition, the value of $P_{2R}$ (equivalent to $s$ in equation 20) while relatively constant tended to be at a higher value ($\sim 1.9$) than that assumed by SAUND.

5.3 The parameterisation used by Niess and Bertin

Hadron showers, generated by Geant4 (version 4.06 p03), were studied up to energies of $10^{5}$ GeV and electromagnetic showers to higher energies by Niess and Bertin [15, 16]. The hadronic showers were parameterised as follows.

[TABLE]

with

[TABLE]

where $E$ is the energy of the hadron shower, $X_{0}$ is the radiation length in water, $z^{\prime}=z/X_{0}$, $b=0.56$ as determined from the fit and $a$ is chosen to satisfy $z_{max}^{\prime}=(a-1)/b$. Here $z_{max}^{\prime}$ is the depth in radiation lengths at which the shower maximum occurs. This is parameterised as

[TABLE]

with $E_{c}=0.05427$ GeV. The radial distribution function is parameterised as

[TABLE]

where $r_{i}=3.5$ cm, $n=n_{1}=1.66-0.29(z/z_{max})$ for $r<r_{i}$ and $n=n_{2}=2.7$ for $r>r_{i}$. The constant $g_{0}$ is chosen to be $(2-n_{1})(n_{2}-2)/((n_{2}-n_{1})r_{i}^{2})$ so that the integral of the radial distribution is unity.

Figure 10 shows the radial distributions from this parameterisation compared with the predictions of CORSIKA. There is quite good agreement between the two. There is a difference in the normalisation with depth since Geant4, on which this parameterisation is based, produces showers which tend to develop more slowly with depth than those from CORSIKA (see Figure 2). Furthermore, both this and the SAUND parameterisation (Section 5.2) assume a linear variation of the shower peak depth with $\log E$ whereas CORSIKA gives a clear parabolic shape (see equation 6). This is illustrated in Figure 11. The Niess-Bertin parameterisation predicts that $56\%$ of the shower energy is contained within a cylinder of radius 4 cm in reasonable agreement with the value of $51\%$ from CORSIKA (these values are almost independent of energy).

6 The acoustic signals from the showers.

The pressure, $P$, from a hadron shower depositing total energy $E$ at time $t$ resulting from the deposition of relative energy density $\epsilon=(1/E)(1/2\pi r)d^{2}E/drdz$ at a point distant $d$ from the volume, $dV$, follows the form [13]

[TABLE]

where the integral is over the total volume of the shower. Here $\beta=2.0\cdot 10^{-4}$ is the thermal expansion coefficient of the medium at $14^{\circ}$C, $C_{p}=3.8\cdot 10^{3}$ J kg*-1* K*-1* is the specific heat capacity and $c=1500$ ms*-1* is the velocity of sound in the sea water.

Acoustic signals seen by an observer at distance $r$ from the shower centre are computed from equation (25) as follows. Points are produced randomly throughout the volume of the shower with density proportional to the deposited energy density and the time of flight from every produced point to the observer calculated. The flight times to the observer are histogrammed over $2^{n}$ bins (in this case $n=10$ is chosen) centred on the mean flight time and with a suitable bin width, $\tau$ (chosen here to be 1$\mu$s). The counts in each bin of the histogram are divided by $\tau$ yielding the function $E_{xyz}(t)$. The Fourier transform of the pressure wave is then

[TABLE]

using the standard Fourier transform theorem, that taking the derivative in the time domain is the same as multiplying by $i\omega$ in the frequency domain. The Fourier transform $E_{xyz}(\omega)$ at angular frequency $\omega$ is evaluated numerically by a fast Fourier Transform (FFT) from the histogram $E_{xyz}(t)$. A correction is applied for attenuation in the water by a factor $A(\omega)=e^{-\alpha(\omega)r}$ where $\alpha(\omega)$ is the frequency dependent attenuation coefficient. The pressure as a function of time is then evaluated numerically by an inverse FFT using frequency steps from zero to the sampling frequency (the inverse of the bin width $\tau$ i.e. 1 MHz in this case). This gives

[TABLE]

where $\Omega=2\pi/1024$ radians and $\omega_{n}/2\pi=n\Omega/2\pi$ MHz is the $n$th frequency. The attenuation coefficient $\alpha(\omega)$ is computed either according to the formulae in [42] or using the complex attenuation given in [15, 16]. This method of calculation was computationally much faster than the evaluation of the space integral given in equation 18 of reference [13] and gave identical results.

Acoustic pulses, computed with the complex attenuation described in [15, 16], using the parameterisations of the shower profile given above are shown in Figure 12. It can be seen that the parameterisation developed here gives similar results to that described in [15, 16] despite the fact that the latter was an extrapolation from low energy simulations. The parameterisation used by SAUND [39, 40] gives smaller signals concentrated at somewhat lower frequencies.

Further properties of the acoustic signals are shown in Figures 13 to 16. The pulses tend to be somewhat asymmetric with the asymmetry defined by $|P_{max}|-|P_{min}|/|P_{max}|+|P_{min}|$. The complex nature of the attenuation enhances this asymmetry. This is most evident in the far field conditions e.g. at 5km where non complex attenuation would yield a totally symmetric pulse. Figure 13 shows the angular dependence of the peak pressure. Here the angle is that subtended by the acoustic detector relative to the plane, termed the median plane, through the shower maximum at right angles to the axis of the shower. The parameterisation derived here gives a somewhat narrower angular spread than the others. This could be due to the slightly longer showers predicted by CORSIKA than the others. Figure 13 also shows the asymmetry of the pulse as a function of this angle. The pulse initially becomes more symmetric moving out of the median plane and then the asymmetry becomes negative at larger angles. Figure 14 shows the decrease of the pulsed peak pressure with distance from the shower in the median plane and the asymmetry with distance in this plane. Figures 15 and 16 show the frequency composition of the pulses at different angles to the median plane at 1 km from the shower and at different distance in the median plane, respectively.

7 Conclusions

The simulation program for high energy cosmic ray air showers, CORSIKA, has been modified to work in a water or ice medium. This allows both hadron and neutrino showers to be generated in the medium over a wide range of energy ($10^{5}$ to $10^{12}$ GeV). The properties of hadronic showers in water simulated by CORSIKA agree with those from other simulations to within $10-20\%$. A similar uncertainty has been noted previously from the variations in CORSIKA showers in air generated by different models of the hadron interactions. However, none of the other available simulations for water cover the range of energies accessible to CORSIKA. The hadronic showers produced by neutrino interactions are shown to have similar profiles to proton showers which deposit the same amount of energy to that from the neutrino and which start at the interaction point of the neutrino. The properties of the neutrino interactions are described. A parameterisation of the shower profiles generated by CORSIKA is given. There is reasonable agreement with the parameterisation based on the Geant4 simulations at low energy ($<10^{5}$ GeV) developed by Niess and Bertin. However, the agreement with the parameterisation used by the SAUND Collaboration, which is based on the NKG formalism, is less good. The position of the shower maximum, determined from the CORSIKA program, is found to vary quadratically with $\log E$ rather than linearly as assumed in the latter two parameterisations.

The acoustic signals generated by neutrino interactions using CORSIKA and by the two other parameterisations are described and their properties are studied. The acoustic signal is found to be very sensitive to the energy deposited close to the shower axis.

7.1 Acknowledgments

We wish to thank Ralph Engel, Dieter Heck, Johannes Knapp and Tanguy Pierog for their assistance in modifying the CORSIKA program. We also thank Valentin Niess and Justin Vandenbroucke for valuable discussions.