A.1 Simulation of the Cascade Development
Already in a previous analyses, the simulation of the energy deposit and the corre-sponding Cherenkov light emission was adapted to reflect the longitudinal devel-opment of cascades [98]. Because the software used at that time was not applicable in this work, a new cascade simulation software was developed. For energies below O(1 PeV) it follows the existing algorithm, where a single cascade is represented by a sequence of sub-cascades (Figure 3.8 on Page 34). Each sub-cascade is dis-placed in space and time and attributed with an energy according to the energy loss profile of low energetic cascades given by equation (3.16). The shift in space
∆x follows the direction of the parent particle and is set to a few X0, typically three radiation length. The corresponding time shift is ∆t = ∆x/c. The time integrated total light flux Φ emitted at shower depth x is given by the sum of fluxesφi from each sub-cascades at depth xi
Φ(x, t) =
n
X
i=0
φi(xi, t)E(xi)H(t−ti), (A.1)
wherenis the number of sub-cascades,E(xi) is the integrated energy loss between two sub-cascades spaced by ∆x: E(xi) =Rxxi
i−∆x dE(x)
dx dx0i andH(t−ti) the Heavi-side function shifted byti=t0+i∆t. The Heaviside function accounts for the fact that only those sub-cascades contribute which have already emerged, hence, not only the spacial distribution of light changes due to the longitudinal development, but also the time distribution is affected.
A.2 Simulation of the longitudinal Energy Deposit at extremely-high Energies
For the simulation of the energy deposit profile at extremely-high energies, differ-ential cross sections for bremsstrahlung and pair creation processes, including the LPM suppression effects derived by Midgal are used. These are shown in Figure 3.7 and 3.6 for different energies. The lower plots show the radiation length and mean free path, respectively, which have been obtained by numerical integration of the cross sections at distinct energies. For the simulation a cubic spline function is used to interpolate between different energies.
Cross section formulas
The cross section formulas used for the simulation have been taken from [91], where they are given per nucleus with charge Z. For a compound material like water the number of nuclei per unit volume is calculated taking into account the average charge per nucleus hZ/Ai which is∼0.56 for water:
N = NA A· hZ/Ai ρ,
where NA is Avogadro’s number, A the atomic weight, and ρ the density of the material. In the cross section formulas the atomic charge appears with exponents.
For compound materials the averaged charge, has to be used where components are weighted with the exponent.
Simulation program flow
The simulation is performed in one dimension, all particles are propagated and secondaries are created until they fall below a threshold energy ofO(10 TeV) where the energy deposit of the particle can be well described by the parametrization of the energy loss profile given by equation (3.16). The total energy loss is then calculated by summation of all individual energy loss profiles appropriately shifted in space. The program flow of the simulation is given as follows:
1. determine the interaction point given by an exponential distribution with a scale parameter given by the radiation length (mean free path);
2. create bremsstrahlung photon (electron-positron pair) and sample energy transfer from the differential cross section;
3. a) remaining particles which are below the energy threshold are deposited and the energy loss profile obtained from the parametrization is added to the total energy loss at the current depth;
b) remaining particles above the threshold are propagated further, start-ing at 1.
A difficulty is the sampling of the energy transfer according to the differential cross section. Since the parameter space is large it is not feasible to use the inverted histogram representation of the cross section at different energies to draw random variates from it. Here, a Markov chain algorithm is used to generate random variates following the differential cross sections.
Metropolis-Hastings algorithm
The Metropolis-Hastings [114, 146] algorithm allows to sample from an arbitrary probability density function (PDF), requiring only that a function P(x) propor-tional to the density can be calculated at x. A Markov chain is constructed in which each state xt+1 depends only on the previous state xt. It uses a proposal density function Q(x0|xt), which depends on the current state xt, to generate a new proposal sample x0. The transition from xt to x0 is accepted according to a likelihood given by the ratio
L= P(x0)Q(xt|x0)
P(xt)Q(x0|xt). (A.2)
With this ratio, one determines whether the new proposalx0 is more likely (L ≥1) than the current state xt, and if so, it is accepted as the new state xt+1. If it is not, it is only accepted with the probability ofL.
The proposal statex0 can be drawn from an arbitrary distributionQ(x0), how-ever, if the proposal state is drawn from a distribution which is similar to the target PDF, less proposals are rejected and the algorithm is very efficient. Here, the beta distribution [119, 136, 141] is chosen for Q. Using appropriate shape parameters, the distribution matches the cross section functions quite well, n.b. it does not depend on the previous state1. The shape parameters have been obtained by fitting the beta function to the cross sections. The beta function is related to the beta distribution and has identical shape parameters. Although the agreement is not sufficient in terms of a χ2 test, it is good enough to draw proposal states from it.
Results of the simulation are depicted in Figure 3.9 on 35.
1which is then called an independence chain Metropolis-Hastings algorithm