Table Of ContentA GPU-based Calculation Method for Near Field Effects of Cherenkov Radiation
Induced by Ultra High Energy Cosmic Neutrinos
Chia-Yu Hu,1,2,∗ Chih-Ching Chen,1,2,† and Pisin Chen1,2,3,4,‡
1Graduate Institute of Astrophysics, National Taiwan University, Taipei, Taiwan 10617
2Leung Center for Cosmology and Particle Astrophysics (LeCosPA),
National Taiwan University, Taipei, Taiwan, 10617
3Department of Physics, National Taiwan University, Taipei, Taiwan 10617
4Kavli Institute for Particle Astrophysics and Cosmology,
SLAC National Accelerator Laboratory, Menlo Park, CA 94025, U.S.A.
The radio approach for detecting the ultra-high energy cosmic neutrinos has become a mature
field. The Cherenkov signals in radio detection are originated from the charge excess of particle
showers due to Askaryan effect. The conventional way of calculating the Cherenkov pulses by
0
makingFraunhoferapproximationfailswhenthesizesoftheelongatedshowersbecomecomparable
1
with the detection distances. We present a calculation method of Cherenkov pulses based on the
0
finite-difference time-domain (FDTD) method, and attain a satisfying effeciency via the GPU-
2
acceleration. Ourmethodprovidesastraightforwardwayofthenearfieldcalculation,whichwould
n beimportantforultrahighenergyparticleshowers,especaillytheelectromagneticshowersinduced
a by the high energy leptons produced in the neutrino charge current interactions.
J
9 PACSnumbers:
2
] I. INTRODUCTION rent or neutral current. The former also produces a lep-
E
tonwithcorrespondingflavor. Boththehighenergylep-
H tonsandthehadronicdebrisinduceparticleshowers. As
Cosmic neutrinos, as a probe of the universe to the
. proposed by Askaryan in the 1960’s [5], the high energy
h highest energy regime, are wonderful in many respects.
particle shower develops in a dense medium would have
p Due to their extremely small interaction cross section,
- net negative charges. This charge imbalance appears as
they can penetrate through galactic infrared (IR) and
o a result of the knocked-off electrons being part of the
r cosmic microwave background (CMB) photons, while
st photons of energy above 10 TeV would be attenuated. shower, as well as the positrons in the shower annihilat-
ingwiththeelectronsofthemedium. Thenetchargesof
a Furthermore, being uncharged, they propagate along
[ theshowers,typically20%oftotalshowerparticles,serve
straightlinesandthereforeareabletopointdirectlyback
asasourceemittingtheCherenkovradiationswhenthey
1 to their sources, while protons or other charged particles
travel in the medium. The sizes of the showers are quite
v would be deflected by the magnetic field in the universe.
localized (tens of cm in radial and few meters in longi-
1
4 Ultra-high energy cosmic rays (UHECRs) have been tudinal development) compared to those develope in the
3 observedupto≈1019.6eV.Thesourceofsuchamazingly air(kmscale),andthereforeresultincoherentradiations
5 energeticeventshaveremainedamystery. Abovethisen- for wavelengths longer than the shower sizes. The cor-
. ergyscale,UHECRsinteractwithCMBphotonsthrough responding coherent wavelength turns out to be in the
1
0 the Greisen-Zatsepin-Kuzmin(GZK) processes [1]. The radio band, from hundreds of MHz to few GHz. This
0 GZK cut-off of the cosmic ray energy spectrum has been Askaryan effect has been confirmed in a series of exper-
1 first observed by the High Resolution Fly’s Eye Exper- iments at Stanford Linear Accelerator Center (SLAC),
: iment [2] and later confirmed by the Pierre Auger Ob- using different dense media such as silica sand, rock salt
v
i servatory [3], so the corresponding GZK neutrinos are and ice [6–9].
X almost guaranteed to exist. Nevertheless, none of these
r have been observed so far. Detecting the GZK neutrinos
a
provides critical informations for unraveling the mystery
II. COHERENT CHERENKOV PULSES
oftheoriginandevolutionofthecosmicaccelerators,and
will be one of the utmost tasks in the coming decade [4].
In the radio detection experiment, the signals come
One promising way of detecting UHE neutrinos is the
from the Cherenkov radiations of the net charges in the
radio approach. When an ultra-high energy cosmic neu-
shower. The key concept which makes this technique
trino interacts with ordinary matters on the Earth, it
possible is the coherent emission. In fact, the Cherenkov
would lead to a hadronic debris, either by charged cur-
radiation is a broad band emission and the intensity in-
creases as frequency. For a single charged particle, the
Cherenkov signal in the radio band should be the weak-
est in the spectrum. It is the compact size of the shower
∗Electronicaddress: [email protected]
†Electronicaddress: [email protected] thatmakesradiosignalsospecial. Thecoherentemission
‡Electronicaddress: [email protected] greatly enhences the signal strength in the radio band.
2
TheelectricfieldofCherenkovradiationscanbecalcu- III. NUMERICAL METHOD
lated by solving the inhomogeneous Maxwell equations,
as it has been demonstrated in the paper of Zas, Halzen,
Numerical algorithms for calculating electromagnetic
and Stanev [10]. The vector potential can be obtained
fields have been existed for decades. However, it was
by the Green’s function method:
notuntiltherecentrapidgrowthofcomputationalpower
A(cid:126) ((cid:126)x)= 1 (cid:90)d3x(cid:48) exp(ik|(cid:126)x−(cid:126)x(cid:48)|)(cid:90)dt(cid:48) exp(iωt(cid:48))J(cid:126)(t(cid:48),(cid:126)x(cid:48)) that this approach became wildly adopted. Among all
ω 4π(cid:15)0c2 |(cid:126)x−(cid:126)x(cid:48)| theexistingalgorithms,thefinite-differencetime-domain
(1) (FDTD hereafter) method has several advantages:
where k is the wavenumber, (cid:126)x the position of the de-
tector, (cid:126)x(cid:48) the position of the shower particles, and J(cid:126)the • It is exceptionally simple to be implemented by
current sources. Adopting the Fraunhoffer approxima- computer programs.
tion
exp(ik|(cid:126)x−(cid:126)x(cid:48)|) exp(ik|(cid:126)x|−ikxˆ·(cid:126)x(cid:48)) • It is a time-domain approach that is well suitable
≈ , (2)
|(cid:126)x−(cid:126)x(cid:48)| R for an impulse signal. A broad band impulse can
where R is the absolute value of (cid:126)x, the integration in be calculated in one single run.
Eq.( 1) can be greatly simplified and therefore enhances
the computational effeciency in the Monte Carlo simula- • The algorithm itself is inherently parallel and the
tion[10–12,17]. Thevalidityofthisapproximationrelies effeciencycanbelargelyimprovedviaparellelcom-
on several length scales: the detection distance (R), the puting.
spatial size of the shower (l) and the wavelengths of in-
terest (λ). The Fraunhoffer approximation works well The idea of FDTD was first proposed by Yee in the
under the condition 1960’s [21], and has been in use for many years for
l2 the electromagnetic impulse modeling. Like most of the
sin2θ(cid:28)1 (3)
λR numerical finite difference methods, the space are dis-
cretizedintosmallgridsandfieldsarecalculatedoneach
where θ is the observational angle between the shower
grid by solving Maxwell equations. Adopting a special
axis and the observational direction.
latticearrangment(knownastheYeelattice),theE-field
However, for the ultra-high energy showers the lon-
andH-fieldarestaggeredinbothspaceandtimeandcan
gitudinal development is longer, especially for the
be calculated in a leapfrog time-marching way.
electromagnetic showers that suffer from the Landau-
Pomeranchuk-Migdal (LPM) suppression [13, 14]. Elec-
tromagnetic showers can be produced by the charge cur-
A. Algorithm
rent generated leptons. For a electromagnetic shower
of EeV-scale energy, The impact of LPM effect on the
The Maxwell curl equations in differential forms are
shower development has been investigated by Monte
Carlo simulations [15–17]. Electromagnetic showers of ∂H
primary energies 1020 eV can be extended to about 200- ∇×E=−µ , (4)
∂t
m long with great fluctuations. In such cases, the far
∂E
field condition cannot be satisfied for distance up to sev- ∇×H=(cid:15) +σE . (5)
∂t
eral kilometers, while the typical detection distance for
ground array detectors is about 1 km due to the attenu-
The FDTD method approximates derivatives by finite
ation length of radio signals in ice. Under these circum- differences. The central difference is adopted to achieve
stances,theFraunhofferapproximationisclearlyinvalid, 2nd order accuracy in both spatial and temporal deriva-
and one has to deal with the complicated integration in tives. We assume cylindrical symmetry along the shower
Eq.( 1). axis (defined as z-axis), and therefore all the derivatives
In the paper of Buniy and Ralston [18], the correction with respect to φ vanish. In addition, due to the polar-
has been made by the saddle point approximation, while izationpropertyofCherenkovradiations,theHr,Hz,Eφ
components also vanish. This can save large amounts
it still cannot cope with the extreme cases for R∼l, the
of computer memories as well as calculation time. Fig-
nearfieldregime. Wehandlethisproblembyanumerical
ure(1)showstheconfigurationofthelatticeunderthese
method based on first principle so that the near field
assumptions. Maxwell equtions in a cylindrical coordi-
radiations can effeciently obtained. Although far field
nate are reduced as:
radiations would be more time consuming, it is not our
∂H ∂E
focus in this paper. − φ =(cid:15) r +σE (6)
∂z ∂t r
Hadronic showers, on the other hand, are less affected
1∂(rH ) ∂E
by LPM effect [19, 20], since the sources of electromag- φ =(cid:15) z +σE (7)
r ∂r ∂t z
netic components in hadronic showers are the decay of
∂E ∂E ∂H
neutralpions, which tendto be interactwith mattersin- r − z =−µ φ (8)
∂z ∂r ∂t
stead of decay at energy above 6.7 PeV. The far field
condition is well fulfilled for hadronic showers in most We can use Eq.( (6)) and Eq.( (7)) to update the E-
parctical cases. field, and then use Eq.( (8)) to update the H-field. The
3
spatial and temporal grid sizes have been chosen such Cherenkov radiations. We define a spherical coordinate
thatthenumericalstabilityissatisfiedandthenumerical whoseoriginliesontheintersectionofthezaxisandthe
dispersion is controlled at an acceptable level [22, 23]. shower maximum. Choosing the desired detector posi-
tions by varying detection angle (θ) and detection dis-
tance (R), and record the electric field values as time
evolves. Afteronesinglerun,wecanhaveallthesimula-
tionresultswewant. Hereonecanseethemeritofusing
a time-domain calculation method.
FIG. 1: The lattice configuration we adopted.
B. Simulation Setup
InordernottolosethefocusontheRFcalculation,we
FIG. 2: Cartoon for our simulation set up.
simplyassumeashowermodelforthelongitudinaldevel-
opment rather than invoking the Monte Carlo packages.
For electromagnetic showers, the well known Nishimura-
Kamata-Greisen(NKG)parametrizationformula[24]de-
scribing the number of shower particles is IV. PERFORMANCE IMPROVEMENT VIA
THE GPU PARALLELIZATION
(cid:20)(cid:18) (cid:19) (cid:21)
0.31 3 X
N (X)= exp 1− lns (9)
e (cid:112)lnE0/Ec 2 XR Recently, graphical processing units (GPU) have be-
came more and more important in the field of high per-
where E is the energy of the primary particle and E
0 c formancecomputation. Undertheneedsinthecomputer
the critical energy, X the slant depth, X depth as
max
showerreachingitsmaximumnumber,sthe(dimension- game markets, the GPUs are currently under booming
less) shower age defined as s = 3X/(X +X ). The developments. The GPUs have become programmable
max
NKG formula can be fit by a Gaussian distribution as: via the Compute Unified Device Architecture (CUDA)
providedbyNVIDIA,andhavebeenimplementedinvar-
z2 ious scientific areas such as molecular dynamics, gravi-
N (z)=N exp(− ) (10)
e max 2l2 tational N-body simulations and lattice QCD [25–29].
where N is the particle number at shower maximum CUDA is an extension of C language, one of the most
max
and l is the longitudinal shower length. The NKG for- popular high level languages in the world. Programmer
mula for E = 10 TeV has l ∼ 2 m in ice. For charge familiar with C language can utilize GPU computations
0
distribution in a snapshot, we assume the Gaussian dis- by simply calling the functions from CUDA. The gen-
tribution in both radial (r) and longitudinal (z) direc- eral parallelization strategies can be found in the CUDA
tions: Programming Guide (available on their webpage).
N (cid:18) (z−X)2(cid:19) (cid:18) r2 (cid:19) Because of the multicore architecture, GPUs are
n(z,r)= e exp − exp − (11) ideal for implementing parallel algorithms. The FDTD
(2π)1.5σ σ2 2σ2 2σ2
z r z r
method is therefore the ideal candidate to benefit from
wherez andr arethecylindricalcoordinatesintheunit GPUcomputing. Thetexturememoriesprovidethepos-
ofg/cm2,andσz andσr arethestandarddeviationofthe sibility to realize the memory cache speedup. Since the
distribution in z and r direction respectively. We choose texturememoriesinCUDAareread-only,webindthe1D
σr = 5 cm in our simulations according to the Moliere cuda array to E-field and H-field alternatively, and write
radius in ice. The σz is generally even smaller than the backtoglobalmemoriesfortheunboundE/H-field. The
Moliere radius and thus has no significant effect on the basic steps are like the following:
radiation pattern.
Figire (2) is a cartoon that demonstrates the setup. 1. Allocate global memories to store E-field and H-
The shower travels in the +z direction emitting field.
4
2. Bind texture to H-field. whiletheydeviatefromeachotheratlowfrequencies. In
Fig.( 4), the spectra at R = 25 m and R = 50 m are
3. Calculate and store the updated E-field by reading
shown with the latter scaled by a factor of two. In this
the cached H-field.
case,thetwospectramatchwellatlowfrequencies,while
they deviate at high frequencies. According to these two
4. Unbind H-field.
examples, one can see that the behavior in high frequen-
√
5. Bind texture to E-field. ciessuggeststhatEω ∝1/ R,namelyacylindricalwave
behavior. On the other hand, the behavior in low fre-
6. CalculateandstoretheupdatedH-fieldbyreading quencies suggests E ∝ 1/R, a spherical wave behavior.
ω
the cached E-field.
7. Unbind E-field.
The different behaviors in different frequency regimes
8. Repeat steps (2) - (7). can be interpreted in a physical way. It is a known fact
thatwavesofhigherfrequenciesarelesslikelytodiffract.
We have investigated the performances of our codes Therefore, the high frequency waves in Cherenkov radi-
withonesingleNVIDIAGTX285graphiccard,whichhas ations are more confined in the θ-direction, and their
240 cores and 933 floating point operations per second energies can only spread into the Cherenkov cone. Their
(GFlOPS) of theoretical peak performance. The follow- energiesgolike1/Rduetogeometricalreason,andhence
√
ingtableshowstheperformancesofourcodeswithdiffer-
1/ Rforthefields. Forthelowfrequencywaves,diffrac-
enttotalgridnumbers(N),comparedtoIntelQuadCore
tion allows another direction (the θ-direction) for their
i7 920 at 2.66GHz (sequential code without CPU paral-
energies to disperse, and therefore decrease faster.
lelization). The GPU system is provided by the Center
for Quantum Science and Engineering of National Tai-
wan University (CQSE).
The fact that the higher frequency regime decreases
slower implies that there is a shift of the peak frequecy
N t t speed up GPU performance
grid CPU GPU in different R. The peak will migrate to the higher
(sec) (sec) (tCPU/tGPU) (GFLOPS) frequency regime as distance moves further away. The
2562 5.15 0.040 128.75 57.0 shapeofthespectrumisR-dependent,andasimplescal-
5122 34.94 0.186 187.85 98.1 ing relation of R is no longer valid here.
10242 294.41 1.167 252.28 125.1
20482 2122.80 8.372 253.56 139.5
40962 15751.17 67.617 232.95 138.2
81922 120617.48 681.082 177.11 109.8
TABLE I: Performance of CPU and GPU codes.
The algorithm of FDTD method is inherently paral-
lel since each grid can be updated independently, so the
advantage of GPU systems can be maximally utilized in
this algorithm. The performance test shows a tremen-
dous speed up using GPU parallelization, which allows
the calculation to be done in very reasonable time, with
a cost-efficient computer resource.
V. RESULTS
A. At the Cherenkov Angle
First we compare the E-fields fixed at Cherenkov an-
gle (θ ) with various distances (R). The magnitude of
c
the E decreases as R increases. However, the decreas-
ω
ing speeds in different frequencies are not the same. For FIG. 3: Eω spectra at R =√25 m and R = 50 m with the
example,Fig.(3)showstheE spectraatR=25mand latter scaled by a factor of 2; these two spectra match
ω
R = 50 m, but with the latter multiplied by a factor of well at hig√h frequencies, implying a cylindrical behavior
two. These two spectra match well at high frequencies (Eω ∝ 1/ R).
5
FIG. 4: E spectra at R = 25 m and R = 50 m with the FIG. 5: The E - R relation of 100MHz. The blue curve
ω ω
latterscaledbyafactorof2;thesetwospectramatchwell is the E-field. The green and red curve show the cylin-
√
at low frequencies, implying a spherical behavior (E ∝ drical (E ∝ 1/ R) and spherical (E ∝ 1/R) behavoir
ω ω √ω
1/R). respectively. TheE-fieldgoeslike1/ R inthenearfield
regime, while approaches to (1/R)-behavoir in the far
field regime. Note the smooth transition from near field
to far field.
B. Angular Distribution
The diffraction effect can also be seen in the angular
distribution. Figures( 6) and ( 7) show the angular dis-
tributions of E at R = 4 m and R = 64 m respectively.
ω
For the case at R = 4 m, the distance is too close for
the waves to diffract. In fact, the radiation pattern in
the near field (before diffraction happen) is just a fuzzy
image of the radiating source. We can see the waves of
higherfrequencieshavestrongermagnitudes,whichisthe
In principle, all waves will eventually diffract and be- inherentcharacterofCherenkovradiations. Notethedis-
have like spherical waves no matter how high their fre- tribution is not symmetric on two sides, since the θ <θ
c
quencies are, if we set the distance R infinitely large. part is closer to the shower axis than the θ >θ part is.
c
Therefore, the terms ”high” and ”low” frequencies are This asymmetry can be understood as we are using the
only relative concepts. From the far field condition in spherical coordinate to describe a sysytem that is actu-
Eq.( 3) we can see that all three length scales couple allycylindricalsymmetric. ForthecaseatR=64m,the
with each other. The waves start to diffract after they detection distance is long enough for waves diffract into
propagate to the distance large enough such that the far the θ-direction. The lower the frequency is, the wilder
fieldconditionissatisfied. Thischaractercanbedemon- ofitsangulardistributionis,whichisthestandardprop-
strated more clearly if we plot the E - R relation with erty in diffraction. In time domain, it can be seen in
ω
onesinglefrequency. Figure(5)showshowE decreases Figure( 8) that the pulse at detection angle more away
ω
withR. Atthedistanceveryclosetotheshower,E goes from θ has wilder width. In principle, if the frequency
√ ω c
like1/ R. Asthedistanceincreases,forlargeenoughR, goes to infinity, the angular distribution will be a delta
the radiation can be viewed as a point source and thus function. However the destructively interference of the
have E ∝ 1/R. We can see a smooth transition from lateral distribution suppress these high frequecy compo-
ω
cylindrical behavoir to spherical behavoir. At R large nents. Thedistributionnowlooksmuchmoresymmetric,
enough such that all the waves of different frequencies in since the differences of the distances to the shower axis
theCherenkovpulsereachthefarfieldregime,theshape between θ <θ and θ >θ parts are negligible compares
c c
of the spectrum is fixed and independent of R, and R to R. Namely, as R goes further, the system, originally
becomes just a normalization factor. cylindrical symmetric, becomes more and more spherical
6
symmetric.
FIG. 8: Time domain E-fields at different angles θ at R
= 32 m. The pulse is wilder at angles more away from
θ , implying the fact that lower frequency components
c
diffract more.
FIG.6: E angulardistributionsforR=4m(nearfield).
ω
The diffraction has not yet happened and the radiation
pattern follows the shower image.
C. Comparison
We compare our results with the conventional far field
formula. The one dimensional approach in [11] should
be a reasonable approximation except at the Cherenkov
angle. Substituting our shower model in Eq.( 10), the
E-field can be obtained:
e eikR (cid:90)
E(cid:126)(ω,(cid:126)x)= iω sinθ nˆ dz(cid:48) N (z(cid:48)) eipz(cid:48)
4π(cid:15) c2 R ⊥ e
√ 0
= 4π2(cid:15)πce2 iω sinθ eiRkRNmaxl e−l22p2 nˆ⊥ (12)
0
where the parameter p(θ,ω) = (1−ncosθ) ω/c. Figure
( 9) shows the comparison between the spectra of the
far field formula and our simulation results at R = 64 m
anddifferentθ. Atlowerfrequencyparttheyareingood
agreements, while at high frequency part there are sig-
nificant differences between them. It can be understood
FIG.7: E angulardistributionsforR=64m(farfield). asthedisagreementparthasnotyetreachedthefarfield
ω
The higher the frequency is, the narrower of its angular regime and thus decreases slower. If we set larger R, the
distribution is. Diffraction phenomenon is quite obvious disagreement part will enter the far field regime and will
here. match with the formula.
7
rower than the ordinary ones, and the detection solid
angle is considered to be small. However the far field as-
sumption neglects the shower size and treat it as a point
source which is not fair. A shower of hundred meters
long would in fact generates signals spaneed also hun-
dred meters, and is surely as possible to be detected as
those with compact size.
The idea of using staggered grid configuration to solve
the two coupled first-ordered PDEs is not limited to
the electromagnetic problems. Recently there are also
applications of FDTD method in the acoustic simula-
tions [30, 31] solving the coupled pressure fields and ve-
locity fields. It is possible to simulate signals in the neu-
trino detection experiments using acoustic approaches,
which is another potential field in UHE neutrino detec-
tion [32, 33].
Inthenextgenerationneutrinodetectorsapplyingthe
ground array layout, it is possible to simultaneously de-
tectthehadronicshowerandtheelectromagneticshower
FIG. 9: Comparison between the spectra of the far field
that are induced by one single charged current neutrino
formula and our simulation results at R = 64 m. From
toptobottomarethespaetraatθ +5◦,θ +10◦,θ +15◦, interaction. If both the hadronic shower and the electro-
c c c
θ +20◦ respectively. Thesolidcurvesareourresultsand magneticshowercanbecorrectlyreconstructed,itopens
c
an opportunity to distinguish the electron neutrino from
thedashedcurvesareobtainedfromthefarfieldformula.
others [34], since the two shower vertexes are nearly lo-
cated at the same places. However, the features of near
field are very different from far field, and will face some
VI. SUMMARY AND CONCLUSIONS
detection difficulties. For example, the normal way to
reconstruct the direction of incoming Cherenkov pulses
We have developed a numerical code to calculate the
bythearrivaltimedifferencesbetweenantennasisbased
radiation patterns of Cherenkov signals from near field
on the assumption that the shower is a point source, i.e.
to far field based on the FDTD method. By utilizing
the far field assumption. For a extended shower this as-
GPUparallelcomputation,theeffeciencyofthecodecan
sumption fails and therefore requires a new reconstruc-
be greatly improved to a satisfying level on a comercial
tionmethod. Anygroundbasedneutrinodetectorhasto
graphic card NVIDIA GTX285. This will be useful in
take this near field effect into account in order to recon-
studying the signals originate from an elongated shower
struct signals from extended showers.
of its size comparable to the detection distance, where
the traditional Fraunhofer approximation does not ap-
ply. Signals from the ultra-high energy electromagnetic
showers induced by the electrons produced in neutrino VII. ACKNOWLEDGEMENTS
charged-current interaction are the typical examples, for
they suffer from severe impact of the LPM effect. Our WethankMelinHuangforusefuldiscussionsandTing-
result shows a smooth transition between near field and Wai Chiu for the help in GPU-calculation. This re-
far field pattern. In fact, the FDTD method is more search is supported by Taiwan National Science Coun-
suitableforthecalculationofnearfieldpatternsincethe cil(NSC)underProjectNo. NSC98-2811-M-002-501,No.
far field pattern may challange the computer resources. NSC98-2119-M-002-001,theCenterforQuantumScience
The spectrum and the angular distribution of near field and Engineering of National Taiwan University(NTU-
patternhavequitecomplicatedRdependencesinsteadof CQSE) under Nos. 98R0066-65, 98R0066-69, and US
simple scalings in the case of far field. Department of Energy under Contract No. DE-AC03-
Inthecasesoffarfield,theangulardistributionofsig- 76SF00515. We would also like to thank Leung Center
nals induced by LPM-elongated showers are much nar- forCosmologyandParticleAstrophysicsforthesupport.
[1] K. Greisen, Phys. Rev. Lett. 16, 748 (1966). G. T. Zat- [4] P. Chen and K. D. Hoffman, Astronomy Decadal
sepin and V. A. Kuzmin, JETP. 4, 114 (1966). Survey (2010-2020) Science White Paper, (2009),
[2] R.U.Abbasietal.,Phys.Rev.Lett.100,101101(2008). arXiv:0902.3288.
[3] Yamamoto, T. 2008, International Cosmic Ray Confer- [5] Askaryan,G.A.,Zh.Eksp.Teor.Fiz.41,616(1961)[So-
ence, 4, 335 viet Physics JETP 14, 441 (1962)].
8
[6] D. Saltzberg, P. Gorham, D. Walz, et al., Phys. Rev. Journal, March 1967, pp. 215-234
Lett. 86, 2802 (2001). [23] A. Taflove and S. C. Hangess: Computational Electro-
[7] P. W. Gorham, D. Saltzberg, R. C. Field, et al., Phys. dynamics,TheFinite-DifferenceTime-DomainMethod,
Rev. D 72, 023002 (2005). Third Edition, (Artech House, London, 2005).
[8] P. Miocinovic, et al., Phys. Rev. D 74, 043002 (2006). [24] K. Greisen, in: Prog. of Cosmic Ray Phys., Ed. J.G.
[9] P. W. Gorham, et al., [ANITA Colla.], Phys. Rev. Lett. Wilson,Vol.III,(North-HollandPubl.Co.,Amsterdam,
99, 171101 (2007). 1956)p.1;K.Kamata,J.Nishimura,Prog.Tbeor.Phys.
[10] Zas, E., Halzen, F., and Stanev, T., Phys. Rev. D 45, (Kyoto) Suppl. 6 (1958) 93.
362 (1992). [25] Anderson,J.A.,Lorenz,C.D.,Travesset,A.,J.Comput.
[11] J.Alvarez-Mun˜iz,R.A.Va´zquez,andE.Zas,Phys.Rev. Phys. 227 (10), 5342 (2008).
D D 62, 063001 (2000). [26] E.Gaburov,S.Harfst,andS.PortegiesZwart.,NewAs-
[12] J. Alvarez-Mun˜iz, E. Marqu´es, R. A. Va´zquez, and E. tron. 14, 630 (2009).
Zas, Phys. Rev. D 74, 023007 (2006). [27] Belleman, R.G., Bedorf, J., Portegies Zwart, S.F., New
[13] L.D. Landau, I. Ya. Pomeranchuk, Dokl. Akad. Nauk Astron. 13, 103 (2008).
SSSR,92,535,735(1953);A.B.Migdal,Phys.Rev.103, [28] Hsi-Yu Schive, Yu-Chih Tsai, and Tzihong Chiueh,
1811 (1956). arXiv:0907.3390
[14] S. Klein , Rev. Mod. Phys., 71, 1501 (1999). [29] Ting-Wai Chiu et al., arXiv:0911.5029
[15] V. Niess and V. Bertin Astropart. Phys. 26, 243 (2006) [30] D. Botteldooren, Acoustical finite-difference time-
[16] J. Bolmont et al., Proceedings of the 30th ICRC, July, domain simulation in a quasi-Cartesian grid, J. Acoust.
2007. Soc. Amer., vol. 95, pp. 2313V2319, May 1994.
[17] J.Alvarez-Muniz, C.W.James, R.J.Protheroe, andE. [31] J. G. Maloney and K. E. Cummings, Adaptation of
Zas, Astropart. Phys. 4, 100 (2009). FDTD techniques to acoustic modeling, in 11th Annual
[18] R.V.Buniy,J.P.Ralston,Phys.Rev.D65016003(2002). Review of Progress in Applied Computational Electro-
[19] J.Alvarez-Mun˜iz,E.Zas,Phys.Lett.B434,396(1998) magnetics,vol.2,Monterey,CA,Mar.1995,pp.724V731
[20] S. Bevan et al., Astropart. Phys. 28, 366 (2007). [32] J.Vandenbroucke,G.Gratta,andN.Lehtinen,ApJ.621,
[21] K. S . Yee, Numerical solution of initial boundary value 301 (2005).
problems involving Maxwells equations in isotropic me- [33] V. Aynutdinov et al., Proceedings of the 31st ICRC,
dia, IEEE Trans. Antennas Propagat., vol. AP-14, pp. Lodz, Poland, July 2009, arXiv:0910.0678.
302-307, May 1966. [34] J.Alvarez-Mun˜iz,R.A.Va´zquez,andE.Zas,Phys.Rev.
[22] R. Courant, K. Friedrichs and H. Lewy, ”On the par- D 61, 023001 (1999)
tial difference equations of mathematical physics”, IBM
