Simulation of Ultra High Energy Neutrino Interactions in Ice and Water
October 31, 2018
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
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.
(Submitted to Astroparticle Physics)
http://arxiv.org/abs/0704.1025v1
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 Mi-
crowave 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 asymme-
try 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 fol-
lowed down to very low kinetic energies (∼ 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 en-
ergy is needed. Once the electromagnetic energy in the shower reaches the MeV level (electron
range ∼ 1 cm) the energy can be simply added to the total deposited energy and the simula-
tion 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 105 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 1012 GeV are discussed. The distributions are generated using the air
shower program CORSIKA [10] modified to work in a sea water medium. The salt compo-
nent of the sea water has a negligible effect1 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 pa-
rameterised 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 modifica-
tion of the stopping power formula to allow for the density effect in water 2. 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 com-
ponent 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 3. 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
1The shower maximum was observed to peak at a depth 2.4 ± 1.1% less in sea water than in fresh water with
the same peak energy deposited, for protons of energy 105 GeV.
2The stopping power was computed using the Bethe-Bloch formula [20] and the density effect from the formu-
lae of Sternheimer et al [21].
3The level was set at 1 TeV compared to the characteristic energy for water ELPM = 270 TeV [20].
RMS peak energy deposit to the mean peak energy deposit was observed to be 14% at 105 GeV
reducing to 4% at 1011 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 105 GeV. Hence
the comparison is restricted to energies below this.
Figure 2 shows the longitudinal distributions of proton showers at energies of 104 and 105
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 ∼ 5% at 104 GeV rising to ∼ 10% at energy 105
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 ∼ 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 ∼ 0). It is relatively insensitive to the shape of the longitudinal distribution.
3.2 Comparison with the simulation of Alvarez-Muñiz and Zas
The CORSIKA simulation was also compared with the longitudinal shower profile for protons
computed in the simulation by Alvarez-Muñ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 ∼ 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 simula-
tions for energy less than 105 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 stud-
ied 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 ∼ 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 (Z0), 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, µ and τ 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 µ and τ 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 τ 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,
σ, [27–29]. In the quark parton model of the nucleon for the single vector boson exchange pro-
cess, the differential cross section for CC interactions can be expressed in terms of the measured
structure functions of the target nucleon F2 and xF3 as
dQ2dy
Q2 +M2W
(F2(x,Q
2)(1− y + y2/2)± y(1− y/2)xF3(x,Q
2)) (1)
where GF is the Fermi weak coupling, MW is the mass of the weak vector boson, Q
2 is the
square of the four momentum transferred to the target nucleon, y = ν/E where ν is the energy
transferred to the nucleon (ν = E−E ′ with E and E ′ the energies of the incident and scattered
leptons) and x = Q2/2Mν 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 104 to 1013 GeV.
The structure functions F2 and xF3 are the sum of the quark distribution functions which have
been parameterised by fitting data [30, 31]. It can be shown that Q2 = 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 . M2W/s i.e. at values
well outside the region of the fits to the parton distribution functions (PDFs) which have been
performed for x & 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 ∼ 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/σdσ/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. Fig-
ure 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 con-
structed 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 gener-
ated, 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 105 and 1012 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 105 to 1012 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 · 1011
GeV. These were divided into the shower energy intervals 0.5−2 ·1010, 2−4 ·1010, 4−7.5 ·1010,
0.75− 1.3 · 1011 and 1.3− 2 · 1011. 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 · 1010 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 EW ) 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 COR-
SIKA (averaged over 100 showers depositing the same total energy) is described. Other avail-
able 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
= L(z, EL) · R(r, z, EL) (2)
where the function L(z, EL) represents the longitudinal distribution of deposited energy and
R(r, z, EL) the radial distribution. Here EL is log10E with E the total shower energy.
The function L(z, EL) = dE/dz is a modified
4 version of the Gaisser-Hillas function [37].
This function represents the longitudinal distribution of the energy deposited.
L(z, EL) = P1L
z − P2L
P3L − P2L
(P3L−P2L)
P4L+P5Lz+P6Lz
P3L − z
P4L + P5Lz + P6Lz2
Here the parameters PnL were fitted to quadratic functions of EL = log10E with values
= 2.760 · 10−3 − 1.974 · 10−4EL + 7.450 · 10
−6E2L (4)
P2L = −210.9− 6.968 · 10
−3EL + 0.1551E
L (5)
P3L = −41.50 + 113.9EL − 4.103E
L (6)
P4L = 8.012 + 11.44EL − 0.5434E
L (7)
P5L = 0.7999 · 10
−5 − 0.004843EL + 0.0002552E
L (8)
P6L = 4.563 · 10
−5 − 3.504 · 10−6EL + 1, 315 · 10
−7E2L. (9)
The parameter P1L represents the peak energy deposited and P3L the depth in the z coordinate
at this peak while P2L, P4L, P5L and P6L are related to the shower width and shape in z.
The radial distribution was represented by the NKG function [37]
R(r, z, EL) =
)(P2R−1)
)(P2R−4.5)
where the integral
)(P2R−1)
)(P2R−4.5)
dr = P1R
Γ(4.5− 2P2R)Γ(P2R)
Γ(4.5− P2R)
4 The modification is to replace the shape parameter λ in equation 3.5 of reference [37] by the quadratic
expression in z in equation 3.
The parameter P1R was found to vary strongly with depth while P2R was only a weak function
of depth. The parameters PnR (with n = 1,2) were each represented by the quadratic form
PnR = A +Bz + Cz
2 (11)
and the quantities A,B,C parameterised as quadratic functions of EL. This gave for P1R
A = 0.01287E2L − 0.2573EL + 0.9636 (12)
B = −0.4697 · 10−4E2L + 0.0008072EL + 0.0005404 (13)
C = 0.7344 · 10−7E2L − 1.375 · 10
−6EL + 4.488 · 10
−6 (14)
and for the parameter P2R
A = −0.8905 · 10−3E2L + 0.007727EL + 1.969 (15)
B = 0.1173 · 10−4E2L − 0.0001782EL − 5.093 · 10
−6 (16)
C = −0.1058 · 10−7E2L + 0.1524 · 10
−6EL − 0.1069 · 10
−8. (17)
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
Fi −Di
Fi +Di
where Fi and Di refer to the fitted value and the value observed in the ith bin from the COR-
SIKA showers, respectively. In order to improve the fit at small radii the contributions to χ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 frac-
tional deviations was 3.4% for radii less than 51.25 g cm−2 and for energies greater than 106.5
GeV. The fit becomes poorer at lower energies and greater radii than these. Integrating the pa-
rameterisation shows that the fraction of the total energy computed from the fit within the fit
range was 91% averaged over the deposited energy range 107 to 1012 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
= Ek(
)t exp (t− z/λ) 2πrρ(r) (19)
where zmax = 0.9X0 ln(E/Ec) is the maximum shower depth, X0 = 36.1 g cm
−2 is the radia-
tion length and Ec = 0.0838GeV. The constants t = zmax/λwhere λ = 130−5 log10(E/10
4GeV)
g cm−2 and k = tt−1/ exp (t)λΓ(t). The radial density is given by
ρ(r) =
as−2(1 + a)s−4.5
Γ(4.5− s)
2πΓ(s)Γ(4.5− 2s)
where a = r/rM with rM = 9.04 g cm
−2, the Moliè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 approxi-
mately 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 P1R (equivalent to RM in equation
20) were strongly depth dependent and much lower than the Molière radius in water, assumed
by the SAUND collaboration. In addition, the value of P2R (equivalent to s in equation 20)
while relatively constant tended to be at a higher value (∼ 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 105
GeV and electromagnetic showers to higher energies by Niess and Bertin [15,16]. The hadronic
showers were parameterised as follows.
= rf(z)g(r, z) (21)
f(z) =
(bz′)a−1 exp−bz′
where E is the energy of the hadron shower, X0 is the radiation length in water, z
′ = z/X0,
b = 0.56 as determined from the fit and a is chosen to satisfy z′max = (a − 1)/b. Here z
max is
the depth in radiation lengths at which the shower maximum occurs. This is parameterised as
z′max = 0.65 log(
) + 3.93 (23)
with Ec = 0.05427 GeV. The radial distribution function is parameterised as
g(r, z) = g0
where ri = 3.5 cm, n = n1 = 1.66 − 0.29(z/zmax) for r < ri and n = n2 = 2.7 for r > ri.
The constant g0 is chosen to be (2− n1)(n2 − 2)/((n2 − n1)r
i ) so that the integral of the radial
distribution is unity.
Figure 10 shows the radial distributions from this parameterisation compared with the pre-
dictions 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 Fig-
ure 2). Furthermore, both this and the SAUND parameterisation (Section 5.2) assume a linear
variation of the shower peak depth with logE 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 ǫ = (1/E)(1/2πr)d2E/drdz at a point distant d from the
volume, dV , follows the form [13]
P (d, t) =
δ(t− d/c)
dV (25)
where the integral is over the total volume of the shower. Here β = 2.0 · 10−4 is the thermal
expansion coefficient of the medium at 14◦C, Cp = 3.8 · 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 2n bins (in this case n = 10 is chosen) centred on the mean flight time and with a suitable
bin width, τ (chosen here to be 1µs). The counts in each bin of the histogram are divided by τ
yielding the function Exyz(t). The Fourier transform of the pressure wave is then
P (ω) =
Exyz(t)e
−iωtdt =
Exyz(t)e
−iωtdt =
iωExyz(ω)
using the standard Fourier transform theorem, that taking the derivative in the time domain is the
same as multiplying by iω in the frequency domain. The Fourier transform Exyz(ω) at angular
frequency ω is evaluated numerically by a fast Fourier Transform (FFT) from the histogram
Exyz(t). A correction is applied for attenuation in the water by a factor A(ω) = e
−α(ω)r where
α(ω) 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 τ i.e. 1 MHz in this case). This gives
P (t) =
n=511
n=−512
P (ωn)A(ωn)e
inΩ (27)
where Ω = 2π/1024 radians and ωn/2π = nΩ/2π MHz is the nth frequency. The attenuation
coefficient α(ω) 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 |Pmax| − |Pmin|/|Pmax| + |Pmin|.
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 (105 to 1012 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
(< 105 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 logE 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 Van-
denbroucke for valuable discussions.
References
[1] Proceedings of the Workshop on Acoustic and Radio EeV Neutrino Detection Activities
(ARENA), DESY, Zeuthen (May 2005), Editors R. Nahnhauer and S. Böser
[2] K. Griesen, Phys. Rev. Lett.16 (1966) 748,
G.T. Zaptsepin, V.A. Kuzmin, JETP Lett. 4 (1966) 78.
[3] E. Waxman and J. Bahcall, Phys. Rev. D59 (1999) 023002, (hep-ph/9807282).
[4] R.D. Engel, D. Seckel and T.Stanev, Phys. Rev. D64 (2001) 093010 (astro-ph/0101216).
[5] See for example M. Ackermann et al., (astro-ph/0412347) Phys. Rev. D71 (2005) 077102.
[6] See for example Nucl. Instrum. Meth. A567 (2006) 438
[7] J.A. Aguilar et al., astro-ph/0606229.
[8] See for example G. Aggouras et al., Nucl. Instrum. and Meth. A552 (2005) 420.
[9] See for example Nucl. Phys. Proc. Suppl. 143 (2005) 373.
[10] “CORSIKA: A Monte Carlo Code to Simulate Extensive Air Showers”, D. Heck et al.,
Karlsruhe Report FZKA 6019. (http://www-ik.fzk.de/corsika).
[11] G. Askar’yan, Soviet Physics JETP 14 (1962) 441 and 21 (1965) 658.
[12] J. Alvarez-Muñiz, E. Marqués, R.A. Vázquez and E. Zas Phys. Rev. D68 (2003) 043001
(astro-ph/0206043).
[13] J.G. Learned Phys. Rev. D19 (1979) 3293.
[14] J. Alvarez-Muñiz and E.Zas, Phys. Lett. B434 (1998) 396 (astro-ph/9806098).
[15] V. Niess and V. Bertin astro-ph/0511617 and V. Niess, PhD Thesis, CPPM, Marseille.
[16] V. Niess, PhD Thesis, CPPM, Marseille, see equations 1-55 and 1-56.
http://arxiv.org/abs/hep-ph/9807282
http://arxiv.org/abs/astro-ph/0101216
http://arxiv.org/abs/astro-ph/0412347
http://arxiv.org/abs/astro-ph/0606229
http://www-ik.fzk.de/corsika
http://arxiv.org/abs/astro-ph/0206043
http://arxiv.org/abs/astro-ph/9806098
http://arxiv.org/abs/astro-ph/0511617
[17] Geant4, J. Allison et al., Nucl. Inst. and Meths. in Phys. Research A506 (2003) 250 and
IEEE Transactions on Nucl. Science 53 (2006) 270.
[18] “The EGS4 Code System” W.R. Nelson, H.Hirayama and D.W.O. Rogers, report number
SLAC-265.
[19] L.D. Landau and I.J. Pomeranchuk, Dokl. Akad. Nauk. SSSR 92 (1953) 535 and 92 (1953)
735. These papers are available in English in L. Landau, “The Collected Papers of L.D.
Landau”, Pergamon Press 1965.
A.B. Migdal, Phys. Rev. 103 (1956) 1811.
[20] Particle data table, Phys. Lett. 592 (2004) 1.
[21] R. M. Sternheimer, S.M. Seltzer and M.J.Berger, Atomic Data and Nuclear Data Tables
30 (1984) 261.
[22] D. Heck et al., Forschungszentrum Karlsruhe GmbH, Karlsruhe, Report number FZKA
6019 (1998).
[23] N.N. Kalmykov and S. Ostapchenko, Phys. Atom. Nucl. 56 (1993) 346, N.N. Kalmykov
et al., Nucl. Phys. Proc. Suppl. 52B (1997) 17.
[24] J. Alvarez-Muniz and E. Zas, Phys. Lett. B434 (1998) 396 (astro-ph/9806098).
[25] Influence of Hadronic Interaction Model on the Development of EAS in Monte Carlo
Simulations, D. Heck, J.Knapp and G. Schatz, Nucl. Phys. B (Proc. Suppl.) 52B (1997)
139-141.
[26] “An Introduction to the Physics of Quarks and Leptons” by P. Renton (published by Cam-
bridge University Press, 1990)
[27] J. Kwiecinski, A.D. Martin and A.M. Stasto Acta Phys. Polon. B31 (2000) 1273
(hep-ph/0004109).
[28] A.Z. Gazizov and S.I. Yanush Phys. Rev. D65 (2002) 093003 (hep-ph/0105368)
[29] R. Ghandi, C. Quigg, M.H. Reno, I. Sarcevic, Astroparticle Physics 5 (1996) 81.
[30] A.D. Martin, R.G. Roberts, W.J.Stirling and R.S. Thorne, Eur. Phys. J. C14 (2000) 133
(hep-ph/9907231).
[31] http://www.phys.psu.edu/∼cteq/
[32] J. Kwiecinski, A.D. Martin and A.M. Stasto, Phys. Rev. D59 (1999) 093002.
[33] A.D. Martin, M.G.Ryskin and A.M. Stasto, Acta Phys. Polon. B34 (2003)3273.
[34] R.S. Thorne, private communication.
[35] O. Pisanti, private communication, see also M. Ambrosio et al., astro-ph/0302062.
[36] HERWIG, G. Corcella et al., hep-ph/0011363.
http://arxiv.org/abs/astro-ph/9806098
http://arxiv.org/abs/hep-ph/0004109
http://arxiv.org/abs/hep-ph/0105368
http://arxiv.org/abs/hep-ph/9907231
http://www.phys.psu.edu/~cteq/
http://arxiv.org/abs/astro-ph/0302062
http://arxiv.org/abs/hep-ph/0011363
[37] “Introduction to Ultra High Energy Cosmic Rays” by P. Sokolsky (published by Addison-
Wesley, 1989).
[38] F. James and M. Roos, “Minuit, A System for Function Minimization and Analysis of the
Parameter Errors and Correlations”, Comput. Phys. Commun. 10 (1975) 343.
[39] J. Vandenbroucke, G. Gratta, N. Lehtenin, Astrophys. J. 621 (2005) 301.
(astro-ph/0406105).
[40] J. Vandenbroucke, private communication.
[41] N.G. Lehtinen, S. Adam, G.Gratta, T.K. Berger and M.J. Buckingham (astro-ph/0104033).
[42] M.A. Ainslie and J.G. McColm, J.Acoust. Soc. Am 103 (1998) 1671.
http://arxiv.org/abs/astro-ph/0406105
http://arxiv.org/abs/astro-ph/0104033
Interaction length of photons in water (CORSIKA)
Pair production length
with LPM Effect
(dash dotted curve)
Total Interaction length
including photonuclear
interaction (solid curve)
9/7 X0 Water (No LPM effect)
Log10 Photon Energy/GeV
Corsika - LPM Effect On
Corsika - LPM Effect Off
4 5 6 7 8 9 10 11 12
Figure 1: The interaction length for high energy gamma rays versus the photon energy measured
in CORSIKA (data points with statistical errors). The dash dotted curve shows the pair produc-
tion length computed from the LPM effect using the formulae of Migdal [19]. The solid curve
shows the computed total interaction length, including both pair production and photonuclear
interactions with the cross section from CORSIKA. The dashed line labelled 9/7X0 shows the
expected pair production length without the LPM effect. Here X0 is the radiation length of the
material.
z (cm)
0 200 400 600 800 1000 1200 1400 1600 1800 2000
Average de/dz
Geant4  seawater
Corsika seawater
Average dE/dz of 100 proton showers
z (cm)
0 200 400 600 800 1000 1200 1400 1600 1800 2000
Average de/dz
Geant4  seawater
Corsika seawater
Average dE/dz of 100 proton showers
Figure 2: Averaged longitudinal energy deposited per unit path length of 100 proton showers
at energy 104 GeV (upper plot) and 105 GeV (lower plot) generated in Geant4 and CORSIKA
versus depth in the water.
Geant4  seawater
Corsika seawater
Depth = 250cm
210 Depth = 450cm
Depth = 650cm
210 Depth = 850cm
0 5 10 15 20 25 30 35 40
210 Depth = 1050cm
Geant4  seawater
Corsika seawater
Depth = 250cm
10 Depth = 450cm
) Depth = 650cm
10 Depth = 850cm
0 5 10 15 20 25 30 35 40
10 Depth = 1050cm
Figure 3: Averaged radial energy deposited per 20 g cm−2 vertical slice per unit radial distance
for 100 proton showers at energy 104 GeV (left hand plots) and 105 GeV (right hand plots)
generated in Geant4 and CORSIKA versus distance from the axis in the water for different
depths of the shower.
Average y for different structure function extrapolations
Standard extrapolation (solid)
Scaled extrapolations (dashed)
log Eν/GeV
4 6 8 10 12
Figure 4: The mean value of y as a function of energy for νµ interactions computed according
to the standard model with the PDFs of MRS99 [30], extrapolating x and Q2 out of the fit
range from x = 10−4 linearly on a log-log scale. The upper dashed curve shows the result of
multiplying the PDFs by 1.32log(10
−4/x) for PDFs with x < 10−4 and the lower dashed curve by
dividing by this factor. The deviations of the dashed curve from the solid one is an indication
of the precision of the standard model.
Differential cross section for νµ Scattering
E=104 GeV
E=1013 GeV
E=106 GeV
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Figure 5: The fraction of events per unit y interval for different νµ energies computed by inte-
grating the expressions for the CC and NC cross sections.
x 10 5
10000
x 10 5
x 10 6
x 10 6
Mean EW = 10
10 GeV
Proton (dashed) νµ (solid) showers of same energy
   EW = 3 10
10 GeV
    EW = 5.75 10
10 GeV
    EW = 10
11 GeV
    EW = 1.65 10
11 GeV
Depth gm cm-2
x 10 6
200 400 600 800 1000 1200 1400 1600 1800
Figure 6: The longitudinal distribution of the deposited energy for neutrino showers (solid)
generated by the Herwig-CORSIKA package and proton showers (dashed) scaled to the same
values of shower energy EW . The scaling factors applied to the average of the protons showers
with energy 1010 GeV were 1.0 and 3.0 for EW = 10
10 GeV and EW = 3 · 10
10 GeV, respec-
tively. Those applied to proton showers with energy 1011 GeV were 0.575, 1.0 and 1.65 for
EW = 5.75 · 10
10 EW = 10
11 and EW = 1.65 · 10
11 GeV, respectively.
Depth 250 g/cm2
Proton (dashed) νµ (solid) EW = 3 10
10 GeV
Depth 450 g/cm2
Depth 650 g/cm2
Depth 850 g/cm2
Depth 1050 g/cm2
Radius(cm)
5 10 15 20 25 30 35 40
Depth 250 g/cm2
Proton (dashed) νµ (solid) EW = 1.65 10
11 GeV
Depth 450 g/cm2
Depth 650 g/cm2
Depth 850 g/cm2
Depth 1050 g/cm2
Radius(cm)
5 10 15 20 25 30 35 40
Figure 7: The solid curves show the averaged radial energy deposited per 20 g cm−2 vertical
slice per unit radial distance for 70 neutrino showers with energy transfer EW = 3 · 10
10 GeV
(left hand plots) and EW = 1.65 · 10
11 GeV (right hand plots). The incident neutrino energy
was 2 ·1011 GeV. For comparison the dashed curves show the distributions from proton showers
scaled to these energies. In the left (right) hand plots protons of energy 1010 GeV (1011 GeV)
were scaled by a factor of 3 (1.65).
−100 −80 −60 −40 −20 0 20 40 60 80 100
time (µs)
 Pulse @1km, z=8m (109GeV Primary)
0−1 cm
0−2 cm
Total
Figure 8: The acoustic signal at a distance of 1 km from the shower axis in the median plane
computed from the average of 100 CORSIKA showers each depositing a total energy of 109
GeV in the water. The dotted, dashed and solid curves shows the signals computed from the
deposited energies within cores of radius 1.025 g cm−2, 2.05 g cm−2 and the whole shower
(solid curve), respectively. It can be seen that most of the amplitude of the signal comes from
the energy within a core of radius 2.05 g cm−2.
Solid CORSIKA, dash SAUND parameterisation
Depth 250 g/cm2
106 GeV protons 100 shower average.
Depth 450 g/cm2
Depth 650 g/cm2
Depth 850 g/cm2
Depth 1050 g/cm2
Radius(cm)
5 10 15 20 25 30 35 40
Solid CORSIKA, dash SAUND parameterisation
Depth 250 g/cm2
1011 GeV protons 100 shower average.
Depth 450 g/cm2
Depth 650 g/cm2
Depth 850 g/cm2
Depth 1050 g/cm2
Radius(cm)
5 10 15 20 25 30 35 40
Figure 9: The radial distributions of the deposited energy at different depths from CORSIKA
compared to the parameterisation used by the SAUND collaboration [39] for 106 GeV and 1011
GeV proton induced showers.
Solid CORSIKA, dash NB parameterisation
Depth 250 g/cm2
106 GeV protons 100 shower average.
Depth 450 g/cm2
Depth 650 g/cm2
Depth 850 g/cm2
Depth 1050 g/cm2
Radius(cm)
5 10 15 20 25 30 35 40
Solid CORSIKA, dash NB parameterisation
Depth 250 g/cm2
1011 GeV protons 100 shower average.
Depth 450 g/cm2
Depth 650 g/cm2
Depth 850 g/cm2
Depth 1050 g/cm2
Radius(cm)
5 10 15 20 25 30 35 40
Figure 10: The radial distributions of the deposited energy at different depths from CORSIKA
compared to the parameterisation used by the Niess and Bertin [15, 16] (labelled NB parame-
terisation) for 106 GeV and 1011 GeV proton induced showers.
CORSIKA
SAUND
Log10E/GeV
6 7 8 9 10 11 12
Figure 11: The depth of the shower peak as a function of log10 E from CORSIKA (black points)
for the showers starting at z = 0. The solid curve shows the parameterisation according to
equation 6. The dashed (dash dotted) lines show the values assumed by the SAUND (Niess-
Bertin) Collaborations.
−100 −50 0 50 100
−0.05
time (µs)
The Acoustic Pulse from a 1011GeV Shower at 1km
ACoRNE
SAUND
0 5 10 15 20 25 30 35 40
Frequency (kHz)
Relative Power Spectrum of a Pulse from a 1011GeV Shower at 1km
ACoRNE
SAUND
Figure 12: The left hand plot shows acoustic pulses generated from the parameterisation de-
scribed in section 5.1 labelled ACoRNE, the parameterisation from reference [15, 16] labelled
NB and that from reference [39,40] labelled SAUND. These pulses were evaluated for a hadron
shower from a neutrino interaction depositing hadronic energy of 1011 GeV 1 km distant from
an acoustic detector in a plane perpendicular to the shower axis at the shower maximum (the
median plane). The right hand plot shows the frequency decomposition of the pulses in the left
hand plot.
−6 −4 −2 0 2 4 6
Angle (degrees)
Peak Pressure vs Angle from a 1011GeV Shower at 1km
ACoRNE
SAUND
0 0.5 1 1.5 2 2.5 3
−0.15
−0.05
Angle (degrees)
The Asymmetry of a Pulse from a 1011GeV Shower
ACoRNE
SAUND
Figure 13: The left hand plot shows the variation of the peak pressure in the pulse with angle
from the median plane at 1 km from the shower. The right hand plot shows the variation of
the pulse asymmetry with this angle at the same distance. The curves were computed from the
parameterisations labelled.
Distance (meters)
Peak Pressure vs Distance from a 1011GeV Shower
ACoRNE
SAUND
Distance (Meters)
The Asymmetry of a Pulse from a 1011GeV Shower
ACoRNE
SAUND
Figure 14: The left hand plot shows the decrease of the pulse peak pressure and the right hand
plot the pulse asymmetry, both in the median plane, as a function of distance from the shower
computed from the parameterisations.
0 5 10 15 20 25 30 35 40
Frequency (kHz)
Relative Power Spectrum of a Pulse from a 1011GeV Shower at 1km
0 5 10 15 20 25 30 35 40
Frequency (kHz)
Relative Power Spectrum of a Pulse from a 1011GeV Shower at 1km
Figure 15: The left hand plot shows the frequency decomposition of the acoustic signal, com-
puted from the parameterisation of the CORSIKA showers, at different angles to the median
plane at a distance of 1km from the shower and the right hand plot shows the cumulative fre-
quency spectrum i.e. the integral of the left hand plot.
0 5 10 15 20 25 30 35 40
Frequency (kHz)
Relative Power Spectrum of a Pulse from a 1011GeV Shower
1000m
2000m
0 5 10 15 20 25 30 35 40
Frequency (kHz)
Relative Power Spectrum of a Pulse from a 1011GeV Shower
1000m
2000m
Figure 16: The left hand plot shows the frequency decomposition of the acoustic signal, com-
puted from the parameterisation of the CORSIKA showers, at different distances from the hy-
drophone in the median plane and the right hand plot shows the cumulative frequency spectrum
i.e. the integral of the left hand plot.
April 8, 2007
Simulation of Ultra High Energy Neutrino Interactions in Ice
and Water
(the ACoRNE Collaboration)a
S. Bevan1, S. Danaher2, J. Perkin3, S. Ralph3y, 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
y 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
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.
(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 Mi-
crowave 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 asymme-
try 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 fol-
lowed down to very low kinetic energies (� 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 en-
ergy is needed. Once the electromagnetic energy in the shower reaches the MeV level (electron
range � 1 cm) the energy can be simply added to the total deposited energy and the simula-
tion 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 105 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 1012 GeV are discussed. The distributions are generated using the air
shower program CORSIKA [10] modified to work in a sea water medium. The salt compo-
nent of the sea water has a negligible effect1 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 pa-
rameterised 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 modifica-
tion of the stopping power formula to allow for the density effect in water 2. 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 com-
ponent 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 3. 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
1The shower maximum was observed to peak at a depth 2:4� 1:1% less in sea water than in fresh water with
the same peak energy deposited, for protons of energy 105 GeV.
2The stopping power was computed using the Bethe-Bloch formula [20] and the density effect from the formu-
lae of Sternheimer et al [21].
3The level was set at 1 TeV compared to the characteristic energy for water E
= 270 TeV [20].
RMS peak energy deposit to the mean peak energy deposit was observed to be 14% at 105 GeV
reducing to 4% at 1011 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 105 GeV. Hence
the comparison is restricted to energies below this.
Figure 2 shows the longitudinal distributions of proton showers at energies of 104 and 105
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 � 5% at 104 GeV rising to � 10% at energy 105
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 � 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 � 0). It is relatively insensitive to the shape of the longitudinal distribution.
3.2 Comparison with the simulation of Alvarez-Muñiz and Zas
The CORSIKA simulation was also compared with the longitudinal shower profile for protons
computed in the simulation by Alvarez-Muñ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 � 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 simula-
tions for energy less than 105 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 stud-
ied 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 � 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 (Z0), 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, � and � 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 � and � 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 � 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,
�, [27–29]. In the quark parton model of the nucleon for the single vector boson exchange pro-
cess, the differential cross section for CC interactions can be expressed in terms of the measured
structure functions of the target nucleon F
and xF
)(1� y + y
=2)� y(1� y=2)xF
)) (1)
where G
is the Fermi weak coupling, M
is the mass of the weak vector boson, Q2 is the
square of the four momentum transferred to the target nucleon, y = �=E where � is the energy
transferred to the nucleon (� = E�E 0 with E and E 0 the energies of the incident and scattered
leptons) and x = Q2=2M� 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 104 to 1013 GeV.
The structure functions F
and xF
are the sum of the quark distribution functions which have
been parameterised by fitting data [30, 31]. It can be shown that Q2 = 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 . M2
=s i.e. at values
well outside the region of the fits to the parton distribution functions (PDFs) which have been
performed for x & 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 � 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=�d�=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. Fig-
ure 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 con-
structed 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 gener-
ated, 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 105 and 1012 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 105 to 1012 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 � 1011
GeV. These were divided into the shower energy intervals 0:5�2�1010, 2�4�1010, 4�7:5�1010,
0:75� 1:3 � 10
11 and 1:3� 2 � 1011. 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 � 1010 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
) 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 COR-
SIKA (averaged over 100 showers depositing the same total energy) is described. Other avail-
able 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
= L(z; E
) �R(r; z; E
) (2)
where the function L(z; E
) represents the longitudinal distribution of deposited energy and
R(r; z; E
) the radial distribution. Here E
is log
E with E the total shower energy.
The function L(z; E
) = dE=dz is a modified4 version of the Gaisser-Hillas function [37].
This function represents the longitudinal distribution of the energy deposited.
L(z; E
) = P
z � P
z + P
Here the parameters P
were fitted to quadratic functions of E
= log
E with values
= 2:760 � 10
� 1:974 � 10
+ 7:450 � 10
= �210:9� 6:968 � 10
+ 0:1551E
= �41:50 + 113:9E
� 4:103E
= 8:012 + 11:44E
� 0:5434E
= 0:7999 � 10
� 0:004843E
+ 0:0002552E
= 4:563 � 10
� 3:504 � 10
+ 1; 315 � 10
: (9)
The parameter P
represents the peak energy deposited and P
the depth in the z coordinate
at this peak while P
and P
are related to the shower width and shape in z.
The radial distribution was represented by the NKG function [37]
R(r; z; E
�4:5)
where the integral
�4:5)
dr = P
�(4:5� 2P
�(4:5� P
4The modification is to replace the shape parameter � in equation 3.5 of reference [37] by the quadratic expres-
sion in z in equation 3.
The parameter P
was found to vary strongly with depth while P
was only a weak function
of depth. The parameters P
(with n = 1,2) were each represented by the quadratic form
= A +Bz + Cz
2 (11)
and the quantities A;B;C parameterised as quadratic functions of E
. This gave for P
A = 0:01287E
� 0:2573E
+ 0:9636 (12)
B = �0:4697 � 10
+ 0:0008072E
+ 0:0005404 (13)
C = 0:7344 � 10
� 1:375 � 10
+ 4:488 � 10
�6 (14)
and for the parameter P
A = �0:8905 � 10
+ 0:007727E
+ 1:969 (15)
B = 0:1173 � 10
� 0:0001782E
� 5:093 � 10
�6 (16)
C = �0:1058 � 10
+ 0:1524 � 10
� 0:1069 � 10
: (17)
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
where F
and D
refer to the fitted value and the value observed in the ith bin from the COR-
SIKA showers, respectively. In order to improve the fit at small radii the contributions to �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 frac-
tional deviations was 3:4% for radii less than 51.25 g cm�2 and for energies greater than 106:5
GeV. The fit becomes poorer at lower energies and greater radii than these. Integrating the pa-
rameterisation shows that the fraction of the total energy computed from the fit within the fit
range was 91% averaged over the deposited energy range 107 to 1012 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
= Ek(
exp (t� z=�) 2�r�(r) (19)
where z
= 0:9X
ln(E=E
) is the maximum shower depth, X
= 36:1 g cm�2 is the radia-
tion length andE
= 0:0838GeV. The constants t = z
=�where � = 130�5 log
(E=10
g cm�2 and k = tt�1= exp (t)��(t). The radial density is given by
�(r) =
(1 + a)
s�4:5
�(4:5� s)
2��(s)�(4:5� 2s)
where a = r=r
with r
= 9:04 g cm�2, the Moliè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 approxi-
mately 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
(equivalent to R
in equation
20) were strongly depth dependent and much lower than the Molière radius in water, assumed
by the SAUND collaboration. In addition, the value of P
(equivalent to s in equation 20)
while relatively constant tended to be at a higher value (� 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 105
GeV and electromagnetic showers to higher energies by Niess and Bertin [15,16]. The hadronic
showers were parameterised as follows.
= rf(z)g(r; z) (21)
f(z) =
exp�bz
where E is the energy of the hadron shower, X
is the radiation length in water, z0 = z=X
b = 0:56 as determined from the fit and a is chosen to satisfy z0
= (a � 1)=b. Here z0
the depth in radiation lengths at which the shower maximum occurs. This is parameterised as
= 0:65 log(
) + 3:93 (23)
with E
= 0:05427 GeV. The radial distribution function is parameterised as
g(r; z) = g
where r
= 3:5 cm, n = n
= 1:66 � 0:29(z=z
) for r < r
and n = n
= 2:7 for r > r
The constant g
is chosen to be (2� n
� 2)=((n
) so that the integral of the radial
distribution is unity.
Figure 10 shows the radial distributions from this parameterisation compared with the pre-
dictions 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 Fig-
ure 2). Furthermore, both this and the SAUND parameterisation (Section 5.2) assume a linear
variation of the shower peak depth with logE 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 � = (1=E)(1=2�r)d2E=drdz at a point distant d from the
volume, dV , follows the form [13]
P (d; t) =
Æ(t� d=
)
dV (25)
where the integral is over the total volume of the shower. Here � = 2:0 � 10�4 is the thermal
expansion coefficient of the medium at 14ÆC, C
= 3:8 � 10
3 J kg�1 K�1 is the specific heat
capacity and 
 = 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 2n bins (in this case n = 10 is chosen) centred on the mean flight time and with a suitable
bin width, � (chosen here to be 1�s). The counts in each bin of the histogram are divided by �
yielding the function E
(t). The Fourier transform of the pressure wave is then
P (!) =
using the standard Fourier transform theorem, that taking the derivative in the time domain is the
same as multiplying by i! in the frequency domain. The Fourier transform E
(!) at angular
frequency ! is evaluated numerically by a fast Fourier Transform (FFT) from the histogram
(t). A correction is applied for attenuation in the water by a factor A(!) = e��(!)r where
�(!) 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 � i.e. 1 MHz in this case). This gives
P (t) =
n=511
n=�512
 (27)
where 
 = 2�=1024 radians and !
=2� = n
=2� MHz is the nth frequency. The attenuation
coefficient �(!) 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 jP
j � jP
j + jP
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 (105 to 1012 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
(< 105 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 logE 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 Van-
denbroucke for valuable discussions.
References
[1] Proceedings of the Workshop on Acoustic and Radio EeV Neutrino Detection Activities
(ARENA), DESY, Zeuthen (May 2005), Editors R. Nahnhauer and S. Böser
[2] K. Griesen, Phys. Rev. Lett.16 (1966) 748,
G.T. Zaptsepin, V.A. Kuzmin, JETP Lett. 4 (1966) 78.
[3] E. Waxman and J. Bahcall, Phys. Rev. D59 (1999) 023002, (hep-ph/9807282).
[4] R.D. Engel, D. Seckel and T.Stanev, Phys. Rev. D64 (2001) 093010 (astro-ph/0101216).
[5] See for example M. Ackermann et al., (astro-ph/0412347) Phys. Rev. D71 (2005) 077102.
[6] See for example Nucl. Instrum. Meth. A567 (2006) 438
[7] J.A. Aguilar et al., astro-ph/0606229.
[8] See for example G. Aggouras et al., Nucl. Instrum. and Meth. A552 (2005) 420.
[9] See for example Nucl. Phys. Proc. Suppl. 143 (2005) 373.
[10] “CORSIKA: A Monte Carlo Code to Simulate Extensive Air Showers”, D. Heck et al.,
Karlsruhe Report FZKA 6019. (http://www-ik.fzk.de/corsika).
[11] G. Askar’yan, Soviet Physics JETP 14 (1962) 441 and 21 (1965) 658.
[12] J. Alvarez-Muñiz, E. Marqués, R.A. Vázquez and E. Zas Phys. Rev. D68 (2003) 043001
(astro-ph/0206043).
[13] J.G. Learned Phys. Rev. D19 (1979) 3293.
[14] J. Alvarez-Muñiz and E.Zas, Phys. Lett. B434 (1998) 396 (astro-ph/9806098).
[15] V. Niess and V. Bertin astro-ph/0511617 and V. Niess, PhD Thesis, CPPM, Marseille.
[16] V. Niess, PhD Thesis, CPPM, Marseille, see equations 1-55 and 1-56.
[17] Geant4, J. Allison et al., Nucl. Inst. and Meths. in Phys. Research A506 (2003) 250 and
IEEE Transactions on Nucl. Science 53 (2006) 270.
[18] “The EGS4 Code System” W.R. Nelson, H.Hirayama and D.W.O. Rogers, report number
SLAC-265.
[19] L.D. Landau and I.J. Pomeranchuk, Dokl. Akad. Nauk. SSSR 92 (1953) 535 and 92 (1953)
735. These papers are available in English in L. Landau, “The Collected Papers of L.D.
Landau”, Pergamon Press 1965.
A.B. Migdal, Phys. Rev. 103 (1956) 1811.
[20] Particle data table, Phys. Lett. 592 (2004) 1.
[21] R. M. Sternheimer, S.M. Seltzer and M.J.Berger, Atomic Data and Nuclear Data Tables
30 (1984) 261.
[22] D. Heck et al., Forschungszentrum Karlsruhe GmbH, Karlsruhe, Report number FZKA
6019 (1998).
[23] N.N. Kalmykov and S. Ostapchenko, Phys. Atom. Nucl. 56 (1993) 346, N.N. Kalmykov
et al., Nucl. Phys. Proc. Suppl. 52B (1997) 17.
[24] J. Alvarez-Muniz and E. Zas, Phys. Lett. B434 (1998) 396 (astro-ph/9806098).
[25] Influence of Hadronic Interaction Model on the Development of EAS in Monte Carlo
Simulations, D. Heck, J.Knapp and G. Schatz, Nucl. Phys. B (Proc. Suppl.) 52B (1997)
139-141.
[26] “An Introduction to the Physics of Quarks and Leptons” by P. Renton (published by Cam-
bridge University Press, 1990)
[27] J. Kwiecinski, A.D. Martin and A.M. Stasto Acta Phys. Polon. B31 (2000) 1273 (hep-
ph/0004109).
[28] A.Z. Gazizov and S.I. Yanush Phys. Rev. D65 (2002) 093003 (hep-ph/0105368)
[29] R. Ghandi, C. Quigg, M.H. Reno, I. Sarcevic, Astroparticle Physics 5 (1996) 81.
[30] A.D. Martin, R.G. Roberts, W.J.Stirling and R.S. Thorne, Eur. Phys. J. C14 (2000) 133
(hep-ph/9907231).
[31] http://www.phys.psu.edu/ cteq/
[32] J. Kwiecinski, A.D. Martin and A.M. Stasto, Phys. Rev. D59 (1999) 093002.
[33] A.D. Martin, M.G.Ryskin and A.M. Stasto, Acta Phys. Polon. B34 (2003)3273.
[34] R.S. Thorne, private communication.
[35] O. Pisanti, private communication, see also M. Ambrosio et al., astro-ph/0302062.
[36] HERWIG, G. Corcella et al., hep-ph/0011363.
[37] “Introduction to Ultra High Energy Cosmic Rays” by P. Sokolsky (published by Addison-
Wesley, 1989).
[38] F. James and M. Roos, “Minuit, A System for Function Minimization and Analysis of the
Parameter Errors and Correlations”, Comput. Phys. Commun. 10 (1975) 343.
[39] J. Vandenbroucke, G. Gratta, N. Lehtenin, Astrophys. J. 621 (2005) 301. (astro-
ph/0406105).
[40] J. Vandenbroucke, private communication.
[41] N.G. Lehtinen, S. Adam, G.Gratta, T.K. Berger and M.J. Buckingham (astro-ph/0104033).
[42] M.A. Ainslie and J.G. McColm, J.Acoust. Soc. Am 103 (1998) 1671.
Interaction length of photons in water (CORSIKA)
Pair production length
with LPM Effect
(dash dotted curve)
Total Interaction length
including photonuclear
interaction (solid curve)
9/7 X0 Water (No LPM effect)
Log10 Photon Energy/GeV
Corsika - LPM Effect On
Corsika - LPM Effect Off
4 5 6 7 8 9 10 11 12
Figure 1: The interaction length for high energy gamma rays versus the photon energy measured
in CORSIKA (data points with statistical errors). The dash dotted curve shows the pair produc-
tion length computed from the LPM effect using the formulae of Migdal [19]. The solid curve
shows the computed total interaction length, including both pair production and photonuclear
interactions with the cross section from CORSIKA. The dashed line labelled 9=7X
shows the
expected pair production length without the LPM effect. Here X
is the radiation length of the
material.
z (cm)
0 200 400 600 800 1000 1200 1400 1600 1800 2000
Average de/dz
Geant4  seawater
Corsika seawater
Average dE/dz of 100 proton showers
z (cm)
0 200 400 600 800 1000 1200 1400 1600 1800 2000
Average de/dz
Geant4  seawater
Corsika seawater
Average dE/dz of 100 proton showers
Figure 2: Averaged longitudinal energy deposited per unit path length of 100 proton showers
at energy 104 GeV (upper plot) and 105 GeV (lower plot) generated in Geant4 and CORSIKA
versus depth in the water.
Geant4  seawater
Corsika seawater
Depth = 250cm
210 Depth = 450cm
Depth = 650cm
210 Depth = 850cm
0 5 10 15 20 25 30 35 40
210 Depth = 1050cm
Geant4  seawater
Corsika seawater
Depth = 250cm
10 Depth = 450cm
) Depth = 650cm
10 Depth = 850cm
0 5 10 15 20 25 30 35 40
10 Depth = 1050cm
Figure 3: Averaged radial energy deposited per 20 g cm�2 vertical slice per unit radial distance
for 100 proton showers at energy 104 GeV (left hand plots) and 105 GeV (right hand plots)
generated in Geant4 and CORSIKA versus distance from the axis in the water for different
depths of the shower.
Average y for different structure function extrapolations
Standard extrapolation (solid)
Scaled extrapolations (dashed)
log Eν/GeV
4 6 8 10 12
Figure 4: The mean value of y as a function of energy for �
interactions computed according
to the standard model with the PDFs of MRS99 [30], extrapolating x and Q2 out of the fit
range from x = 10�4 linearly on a log-log scale. The upper dashed curve shows the result of
multiplying the PDFs by 1:32log(10
=x) for PDFs with x < 10�4 and the lower dashed curve by
dividing by this factor. The deviations of the dashed curve from the solid one is an indication
of the precision of the standard model.
Differential cross section for νµ Scattering
E=104 GeV
E=1013 GeV
E=106 GeV
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Figure 5: The fraction of events per unit y interval for different �
energies computed by inte-
grating the expressions for the CC and NC cross sections.
x 10 5
10000
x 10 5
x 10 6
x 10 6
Mean EW = 10
10 GeV
Proton (dashed) νµ (solid) showers of same energy
   EW = 3 10
10 GeV
    EW = 5.75 10
10 GeV
    EW = 10
11 GeV
    EW = 1.65 10
11 GeV
Depth gm cm-2
x 10 6
200 400 600 800 1000 1200 1400 1600 1800
Figure 6: The longitudinal distribution of the deposited energy for neutrino showers (solid)
generated by the Herwig-CORSIKA package and proton showers (dashed) scaled to the same
values of shower energy E
. The scaling factors applied to the average of the protons showers
with energy 1010 GeV were 1.0 and 3.0 for E
10 GeV and E
= 3 � 10
10 GeV, respec-
tively. Those applied to proton showers with energy 1011 GeV were 0.575, 1.0 and 1.65 for
= 5:75 � 10
11 and E
= 1:65 � 10
11 GeV, respectively.
Depth 250 g/cm2
Proton (dashed) νµ (solid) EW = 3 10
10 GeV
Depth 450 g/cm2
Depth 650 g/cm2
Depth 850 g/cm2
Depth 1050 g/cm2
Radius(cm)
5 10 15 20 25 30 35 40
Depth 250 g/cm2
Proton (dashed) νµ (solid) EW = 1.65 10
11 GeV
Depth 450 g/cm2
Depth 650 g/cm2
Depth 850 g/cm2
Depth 1050 g/cm2
Radius(cm)
5 10 15 20 25 30 35 40
Figure 7: The solid curves show the averaged radial energy deposited per 20 g cm�2 vertical
slice per unit radial distance for 70 neutrino showers with energy transfer E
= 3 � 10
10 GeV
(left hand plots) and E
= 1:65 � 10
11 GeV (right hand plots). The incident neutrino energy
was 2 �1011 GeV. For comparison the dashed curves show the distributions from proton showers
scaled to these energies. In the left (right) hand plots protons of energy 1010 GeV (1011 GeV)
were scaled by a factor of 3 (1.65).
−100 −80 −60 −40 −20 0 20 40 60 80 100
time (µs)
 Pulse @1km, z=8m (109GeV Primary)
0−1 cm
0−2 cm
Total
Figure 8: The acoustic signal at a distance of 1 km from the shower axis in the median plane
computed from the average of 100 CORSIKA showers each depositing a total energy of 109
GeV in the water. The dotted, dashed and solid curves shows the signals computed from the
deposited energies within cores of radius 1.025 g cm�2, 2.05 g cm�2 and the whole shower
(solid curve), respectively. It can be seen that most of the amplitude of the signal comes from
the energy within a core of radius 2.05 g cm�2.
Solid CORSIKA, dash SAUND parameterisation
Depth 250 g/cm2
106 GeV protons 100 shower average.
Depth 450 g/cm2
Depth 650 g/cm2
Depth 850 g/cm2
Depth 1050 g/cm2
Radius(cm)
5 10 15 20 25 30 35 40
Solid CORSIKA, dash SAUND parameterisation
Depth 250 g/cm2
1011 GeV protons 100 shower average.
Depth 450 g/cm2
Depth 650 g/cm2
Depth 850 g/cm2
Depth 1050 g/cm2
Radius(cm)
5 10 15 20 25 30 35 40
Figure 9: The radial distributions of the deposited energy at different depths from CORSIKA
compared to the parameterisation used by the SAUND collaboration [39] for 106 GeV and 1011
GeV proton induced showers.
Solid CORSIKA, dash NB parameterisation
Depth 250 g/cm2
106 GeV protons 100 shower average.
Depth 450 g/cm2
Depth 650 g/cm2
Depth 850 g/cm2
Depth 1050 g/cm2
Radius(cm)
5 10 15 20 25 30 35 40
Solid CORSIKA, dash NB parameterisation
Depth 250 g/cm2
1011 GeV protons 100 shower average.
Depth 450 g/cm2
Depth 650 g/cm2
Depth 850 g/cm2
Depth 1050 g/cm2
Radius(cm)
5 10 15 20 25 30 35 40
Figure 10: The radial distributions of the deposited energy at different depths from CORSIKA
compared to the parameterisation used by the Niess and Bertin [15, 16] (labelled NB parame-
terisation) for 106 GeV and 1011 GeV proton induced showers.
CORSIKA
SAUND
Log10E/GeV
6 7 8 9 10 11 12
Figure 11: The depth of the shower peak as a function of log
E from CORSIKA (black points)
for the showers starting at z = 0. The solid curve shows the parameterisation according to
equation 6. The dashed (dash dotted) lines show the values assumed by the SAUND (Niess-
Bertin) Collaborations.
−100 −50 0 50 100
−0.05
time (µs)
The Acoustic Pulse from a 1011GeV Shower at 1km
ACoRNE
SAUND
0 5 10 15 20 25 30 35 40
Frequency (kHz)
Relative Power Spectrum of a Pulse from a 1011GeV Shower at 1km
ACoRNE
SAUND
Figure 12: The left hand plot shows acoustic pulses generated from the parameterisation de-
scribed in section 5.1 labelled ACoRNE, the parameterisation from reference [15, 16] labelled
NB and that from reference [39,40] labelled SAUND. These pulses were evaluated for a hadron
shower from a neutrino interaction depositing hadronic energy of 1011 GeV 1 km distant from
an acoustic detector in a plane perpendicular to the shower axis at the shower maximum (the
median plane). The right hand plot shows the frequency decomposition of the pulses in the left
hand plot.
−6 −4 −2 0 2 4 6
Angle (degrees)
Peak Pressure vs Angle from a 1011GeV Shower at 1km
ACoRNE
SAUND
0 0.5 1 1.5 2 2.5 3
−0.15
−0.05
Angle (degrees)
The Asymmetry of a Pulse from a 1011GeV Shower
ACoRNE
SAUND
Figure 13: The left hand plot shows the variation of the peak pressure in the pulse with angle
from the median plane at 1 km from the shower. The right hand plot shows the variation of
the pulse asymmetry with this angle at the same distance. The curves were computed from the
parameterisations labelled.
Distance (meters)
Peak Pressure vs Distance from a 1011GeV Shower
ACoRNE
SAUND
Distance (Meters)
The Asymmetry of a Pulse from a 1011GeV Shower
ACoRNE
SAUND
Figure 14: The left hand plot shows the decrease of the pulse peak pressure and the right hand
plot the pulse asymmetry, both in the median plane, as a function of distance from the shower
computed from the parameterisations.
0 5 10 15 20 25 30 35 40
Frequency (kHz)
Relative Power Spectrum of a Pulse from a 1011GeV Shower at 1km
0 5 10 15 20 25 30 35 40
Frequency (kHz)
Relative Power Spectrum of a Pulse from a 1011GeV Shower at 1km
Figure 15: The left hand plot shows the frequency decomposition of the acoustic signal, com-
puted from the parameterisation of the CORSIKA showers, at different angles to the median
plane at a distance of 1km from the shower and the right hand plot shows the cumulative fre-
quency spectrum i.e. the integral of the left hand plot.
0 5 10 15 20 25 30 35 40
Frequency (kHz)
Relative Power Spectrum of a Pulse from a 1011GeV Shower
1000m
2000m
0 5 10 15 20 25 30 35 40
Frequency (kHz)
Relative Power Spectrum of a Pulse from a 1011GeV Shower
1000m
2000m
Figure 16: The left hand plot shows the frequency decomposition of the acoustic signal, com-
puted from the parameterisation of the CORSIKA showers, at different distances from the hy-
drophone in the median plane and the right hand plot shows the cumulative frequency spectrum
i.e. the integral of the left hand plot.
	Introduction
	Adaptation of the CORSIKA program to a water medium
	Comparison with other simulations
	Comparison with Geant4
	Comparison with the simulation of Alvarez-Muñiz and Zas
	Simulation of neutrino induced showers
	Neutrino-nucleon interaction cross sections.
	A simple generator for neutrino interactions.
	The HERWIG neutrino generator.
	Parameterisation of showers
	Parameterisation of the CORSIKA Showers
	The parameterisation used by the SAUND Collaboration
	The parameterisation used by Niess and Bertin
	The acoustic signals from the showers.
	Conclusions
	Acknowledgments
