W D Apel

J C Arteaga

K Bekk

M Bertaina

J Blümer

H Bozdog

I M Brancus

P Buchholz

C Büttner

E Cantoni

A Chiavassa

F Cossavella

K Daumiller

V De Souza

F Di Pierro

P Doll

R Engel

J Engler

M Finger

D Fuhrmann

P L Ghia

H J Gils

R Glasstetter

C Grupen

A Haungs

D Heck

J R Hörandel

T Huege

P G Isar

K.-H Kampert

D Kang

D Kickelbick

H O Klages

K Link

M Ludwig

P Łuczak

H J Mathes

H J Mayer

M Melissas

J Milke

B Mitrica

C Morello

G Navarra

S Nehls

R Obenland

J Oehlschläger

S Ostapchenko

S Over

N Palmieri

M Petcu

T Pierog

H Rebel

M Roth

G Schatz

H Schieler

F Schröder

O Sima

G Toma

G C Trinchero

H Ulrich

A Weindl

J Wochele

M Wommer

J Zabierowski

P Luczak

Muon Production Height Studies with the Air Shower Experiment KASCADE-Grande

Keywords: Cosmic rays, Air shower, Muons, Muon production height

A large area (128 m 2 ) muon tracking detector, located within the KASCADE experiment, has been built with the aim to identify muons (E µ > 0.8 GeV) and their angular correlation in extensive air showers by track measurements under 18 r.l. shielding. Orientation of the muon track with respect to the shower axis is expressed in terms of the radial and tangential angles, which are the basic tools for all muon investigations with the tracking detector. By means of triangulation the muon production height is determined. Distributions of measured production heights are compared to CORSIKA shower simulations. Analysis of these heights reveals a transition from light to heavy cosmic ray primary particles with increasing shower energy in the energy region of the 'Knee' of the cosmic ray spectrum.

Introduction

Ultra High Energy (UHE) astrophysics includes, as one of its topics, the investigation of UHE Cosmic Ray (CR) particles interacting with the nuclei of the atmosphere and creating Extensive Air Showers (EAS) which are studied by means of ground based experiments of large detection area. To understand the nature of UHE particles it is necessary to measure as many components of the EAS cascade as possible. Particularly important is the measurement of the muon component, because muons, being a result of decays of mesons -the most numerous products of hadronic interactions in an air shower -carry information about those interactions to the observation level. In turn, knowledge of UHE hadronic interactions is a necessary condition for answering many questions formulated in CR research. The angular correlation information on EAS muons measured with the tracking detector is especially closely correlated to the parameters of the interactions in question. A recent compilation of muon study results is given in Ref. [START_REF] Grieder | Extensive Air Showers[END_REF].

In the past, various tracking detectors for measuring CR particles have been put into operation either on the surface of Earth [START_REF] Gress | [END_REF]- [START_REF] Petrukin | NEVOD-DECOR experiment[END_REF] or deep under-ground [START_REF] Ahlen | [END_REF]. Currently active installations for muon and neutrino studies are IceCube [START_REF] Gaisser | IceCube Coll., Proc. 30 th ICRC[END_REF], ANTARES [START_REF] Ageron | [END_REF] and MINOS [12]. Tracking detectors provide the possibility to study the muon production height which is the subject of this work.

Muon information has usually been integrated over a large sample of showers and over the whole longitudinal profile. However, muons have some advantage compared with optical photons [13]- [15]: they reflect the development of the nuclear cascade with no mediation from electromagnetic cascades and, similar to the radio emission [16], they can be measured 24 hours a daynot only on clear moonless nights. Their evident disadvantage is that muons are less numerous than photons and are, therefore, subject to large fluctuations. In addition, they are charged particles and are subject to deflection in the geomagnetic field and also suffer from multiple Coulomb scattering in the atmosphere and detector shielding.

Muon tracking allows to measure the composition sensitive profile of a shower in the 'Knee' region, where it is not possible by the fluorescence technique used, e.g. by Fly's Eye [17] and Auger [18] experiments. It enables as well the study of hadron interactions [START_REF] Milke | Proc. 29 th ICRC 2005[END_REF][START_REF] Antoni | [END_REF]. Muons have never been used up to now, with sufficient accuracy and large statistics, for the reconstruction of the longitudinal development of the EAS hadron component of individual showers due to the difficulty of building large area, ground-based muon telescopes. The reconstruction of the longitudinal development of the muon component by means of triangulation [21]- [24] provides a powerful tool for primary mass estimation and for study of high-energy hadron interactions with atmospheric nuclei.

The KASCADE-Grande air-shower experiment [25], set-up in Karlsruhe, Germany provides information on individual EAS up to primary energy 10 18 eV. Complex analyses of KASCADE data [26] resulted in CR flux spectra for five groups of primary masses over the 'Knee' region. However, these spectra exhibit a strong dependence on the hadronic interaction models used in the data analysis. Therefore, an independent investigation of additional shower observables like muon angular correlations with respect to the shower direction is of a great importance. To address this task a Muon Tracking Detector (MTD) [27] was constructed providing additional information on muons for the measured EAS.

In this work we present the experimental investigation of the hadronic cascade in EAS using tracks of muons measured at the KASCADE-Grande experiment. Precise angular information on those tracks is used to determine muon production heights. The sensitivity of this quantity to the mass and energy of CR primary particles is shown. The validity of hadronic interaction models used in Monte Carlo simulations is discussed. Due to the energy dependence of muon production height a transformation to muon production depth is used to investigate, in a model independent way, the CR composition. An estimate of the composition covering two decades of primary energy around the 'Knee' is given.

Muon Tracking Detector (MTD) in KASCADE-Grande

The original KASCADE experiment [28], with its array of scintillation The MTD utilizes streamer tube (ST) gas detectors grouped in modules.

Four modules, three positioned on horizontal planes (top, middle, bottom) and one arranged vertically (wall), form a muon telescope, called a detector tower. The whole MTD comprises 16 towers arranged in two rows of a total detection area of 128 m 2 . A multi-layer filter of six 3 cm thick iron plates, separated each by 5 cm of sand, covered with soil and resting on 30 cm of steel-concrete roof, absorbs a large fraction of accompanying lowenergy particles, thus enhancing the identification of muons with an energy exceeding 0.8 GeV. The details of the design and performance of the MTD can be found elsewhere [27,29].

The MTD shows a good stability of operation. Hereby, the stability of the gas composition on the sub-percent level is of great importance. However, the open gas system is influenced by the atmospheric pressure and temperature (internal in the detector T int and atmospheric T atm ), which is observed in the measured free monitor muon count rate. After applying corrections, the free monitor rate (R) has a remaining variation of about 1%. The resulting temperature and pressure coefficients of the rate, together with the gas composition stability parameters, are given in Table 1. They provide the parameters for the MTD track efficiency determination and correction for temperature and pressure, i.e. for varying atmospheric conditions at KASCADE-Grande. The analysis is based on 3-hit tracks, which are derived from x and y (wire and strip) hits in the 3 horizontal modules for each tower. Using 3-hit tracks only, together with the track quality parameter discussed in Ref. [27],

makes the amount of tracks uncorrelated with a shower being below 1 %.

The modules show, after correction for geometry, pressure and temperature variations, a hit efficiency close to 100%, varying on average by less than 1%.

The track efficiency ( ) stability parameter is given in Table 1. The cluster sizes are limited to 10 cm at the wires and to 16 cm at the strips to control punch-through effects and to achieve tracks of a good quality. More details on track quality and its influence on the direction determination and comparison with shower and detector simulations can be found in Refs. [27], [START_REF] Büttner | Proc. 28 th ICRC 2003[END_REF] and [START_REF] Obenland | [END_REF].

Previously, the anode wire readout in the MTD was operated by combining two neighboring wire cells. To improve the angular resolution a modification was introduced [START_REF] Obenland | KASCADE Coll[END_REF], enabling to identify even and odd numbered wires, using the same common readout electronic channel. The individual wire cells are identified by means of the pulse shape, as illustrated in Fig. 1.

This procedure improves the angular resolution by 30%.

Tangential and Radial Angles

The KASCADE-Grande data, both from the KASCADE and 'Grande' detector arrays, allow us, after detailed analysis steps [START_REF] Antoni | [END_REF], [START_REF] Chiavassa | KASCADE-Grande Coll[END_REF], [START_REF] Apel | [END_REF] to determine the direction (θ and φ), the shower core position, the total number of electrons -electron size lg(N e ), and the total number of muons -muon size lg(N µ ), for each shower. In KASCADE the so-called 'truncated' muon Determination of primary cosmic ray energy from the parameters of registered air showers is of primary interest for all EAS experiments. In KASCADE-Grande total energy E 0 of the primary CR particle, in GeV, can be expressed in an almost mass independent manner through lg(N e ) and lg(N µ ). Based on extensive simulation studies [START_REF] Weber | KASCADE Coll., Proc. 25 th ICRC 1997[END_REF] with CORSIKA [START_REF] Heck | FZKA-Report 6019[END_REF],

using QGSJet01 [START_REF] Kalmykov | [END_REF] the following formula is found

lg(E 0 ) = a• lg(N e ) + b• lg(N µ ) + c , (1)
where the parameters have the following values. For the KASCADE array data analysis a = 0.19, b = 0.79 and c = 2.33, and in addition, lg(N tr µ ) is used instead of lg(N µ ). For the Grande array analysis, following Ref. [START_REF] Glasstetter | Proc. 29 th ICRC 2005[END_REF], we have a = 0.31, b = 0.67 and c = 1.85. Note, that because of the large coefficient for lg(N µ ) this parameter is sometimes used to demonstrate the behavior of various shower quantities as a function of primary energy, where an increase in lg(N µ ) means increase in E 0 . Formula (1) gives ≈ 10-20 % (depending on the primary energy range) uncertainty in energy determination which is sufficient for this study. Investigation of an improved determination of E 0 in KASCADE-Grande is in progress and first results can be found in [START_REF] Haungs | Proc. 31 st ICRC[END_REF] and references therein. For the determination of the muon production height h µ in EAS, orientation of the muon track with respect to the shower axis has to be investigated. Due to their transverse momentum parent mesons travel away from the shower axis. To investigate angular correlation of the muon track with respect to the shower axis, the angle in space between both directions is decomposed into two components: the radial (ρ) and the tangential (τ ) angles [27], [START_REF] Bernlöhr | [END_REF], shown in Fig. 2. The ρ angle is defined in the radial plane, subtended by the shower axis, the shower core position and the point where the muon hits the MTD detector plane. The τ angle is defined in the tangential plane which is perpendicular to the radial plane and parallel to the shower axis, and goes also through the point where the muon hits the detector plane.

The value of τ reflects predominantly the amount of muon scattering in the atmosphere and in any relevant absorber/detector material, and to some extent also the fact that muons can be produced off the shower axis. A contribution from the deflection of high-energy muons in the geomagnetic field (locally < 0. narrower with increasing muon momentum [43], which is expected from the momentum dependence of multiple scattering and from geomagnetic bending.

The ρ angle value, being defined in the shower coordinate system, is dominantly correlated with the transverse momentum of the parent meson.

It is also substantially larger than the corresponding τ values. The multiple scattering of the muon in the atmosphere and muon filter above the MTD contributes also to the value of ρ, but has a relatively minor effect on it.

An example of ρ distributions are shown in Fig. 3 [START_REF] Obenland | [END_REF]. They exhibit an asymmetric shape, being shifted with increasing lg(N tr µ ) towards larger ρ values. This shift is due to the deeper penetration of higher energy showers, but it is not so strongly pronounced because of the wide range of shower core distances selected for Fig. 3. Negative ρ values are due to a finite resolution of the ρ angle (≈ 0.2 • ) and due to the possible muon scattering towards the shower axis.

Muon Production Height

In previous theoretical studies [44] the longitudinal development of the muon component in EAS was investigated by an integral transformation from lateral particle densities to a more shower development related observable like the muon production heights. Due to the very different shower development for light and heavy primary CR particles, their corresponding mean production heights are expected to be different, and thus to be indicative for the primary mass (see Fig. 4).

Conditions of the experiment and the response of the MTD were extensively studied with experimental data and simulations ( [27], [START_REF] Büttner | Proc. 28 th ICRC 2003[END_REF], [START_REF] Obenland | [END_REF])

to ensure the best tracking resolution and to clean the sample from non-

shower-muon tracks. The effects of the non-central location of the MTD within KASCADE were also taken into account.

In order to decrease the effect of multiple muon scattering on their path to the detector |τ | angles smaller than 0.7 • are considered. This cut enriches the sample with high energy muons, above a few GeV [START_REF] Büttner | Proc. 28 th ICRC 2003[END_REF], having smaller multiple scattering angles (σ τ ≈ 0.2 • ) and smaller bending in the geomagnetic field, which improves considerably the determination of the correct production height.

Elimination of tracks not belonging to shower muons or originating in the absorber above the MTD (in addition to the cuts ensuring good track quality, mentioned in Section 2) is achieved by limiting the radial angle value to ρ < 8 • . For muon production height analysis in this work only positive radial angles are considered.

To ensure a nearly full azimuthal symmetry of the measured muons, events with 40 m -80 m core distance to the MTD are selected (the location of the MTD is shifted by 54.65 m to the north from the KASCADE center because this detector was not considered from the beginning). In addition, this selection lowers the influence of the geomagnetic field on the mean muon directions. Below 40 m the punch-through effects forbid a valuable muon tracking and above 80 m the highly asymmetric azimuthal event distribution (see Fig. 6.3 in Ref. [START_REF] Obenland | [END_REF]) is difficult to correct. For Grande reconstructed high-energy events the distance range 250 m -360 m was chosen to avoid extended punch-through and trigger efficiency problems for large distances. For KASCADE (Grande) selections a combined Array-MTD analysis is only reliable above shower size lg(N e ) ≥ 4.8 (6.0) and lg(N µ ) ≥ 4.1 (5.5).

In this work the experimental data on showers with zenith angles θ < 18

•
were analyzed. The mean height values < h µ > were calculated in the altitude range up to 12 km along the shower axis, where the majority of muons with energy exceeding 0.8 GeV at ground are produced. For example, as CORSIKA simulations show, in a 10 16 eV proton induced shower about 90% of such muons are produced in this altitude range. Despite of the fact that the first encounters with nuclei of the atmosphere occur at even higher altitudes, they are not considered in this work, both in data and simulations, because of their large uncertainty. The value of 12 km was chosen because at this altitude, and for the selected distance range of investigated muons to the shower core, the uncertainties in the ρ angle values become comparable to the values themselves. Employing the ρ and τ angles to define the orientation of muon track with respect to the shower axis together with the relative distance R µ from the muon hit to the shower core position, the height h µ can be calculated. This can be done along the muon track, which is the appropriate way to consider absorption or interaction of muons in showers [45,46].

Another approach, used in the present work, is the determination of the height h µ along the shower axis. One way is to use the relation

h µ = R µ / tan ρ
in the shower coordinate system. We could also employ, for angles much smaller than one radian, instead of ρ, the angle ζ = (ρ 2 + τ 2 ) [43], which is the angle in space between muon track and the shower direction.

As shown in Fig. 2 both reconstruction methods lead to muon production heights different from the true one by ∆h ρ µ in case of the calculation with ρ, and by ∆h ζ µ when using the angle ζ. In the former case the muons are assumed to be produced on the shower axis (which is not always the case) while in the latter, the muon production point is reconstructed on a cylindrical surface around the shower axis. The actual value of this difference in both reconstruction methods (being one of the sources of systematic errors)

varies slightly with the type of primary and its energy, leading in most cases to an underestimation of the mean production height by a few percent.

In Table 2 a comparison of these differences is given for two fixed primary energies and two primary species. The results for muons with E µ >0.8 GeV are obtained with CORSIKA simulated vertical showers (500 showers for 10 15 eV and 158 for 10 16 eV primary energy) using the QGSJet01 [START_REF] Kalmykov | [END_REF] and GHEISHA [START_REF] Fesefeldt | [END_REF] interaction models. As seen in columns 3 and 4 the uncertainties range from ≈ 2% to ≈ 8%, showing slightly better results for the 'ρ-method'.

One can expect smaller errors if a correction for multiple scattering in the atmosphere was made in the calculations. The actual amount of scattering is, of course, not known. However, for the selection used in the analysis 2. The investigation of h µ in this work was done using the following relation

h µ = R µ tan (ρ -|τ |) . ( 2
)
The mean muon production height reconstructed with the relation [START_REF] Gress | [END_REF] differs at most by ± 4 % from the true value.

In Fig. 5 an example of muon production height distributions h µ (from Ref. [START_REF] Obenland | [END_REF]), obtained with KASCADE data, shows the behavior of this quantity With increasing shower zenith angle θ, < h µ > is shifting somewhat towards larger values, but the change amounts to less than hundred meters in the whole range of θ. Similar shifts are observed for an analysis using full simulation, only. This effect is small compared to the difference of < h µ > for proton and iron induced showers, which amounts to about 1000 m (see Fig. 4), when using the cuts employed in our analysis.

Muon Production Heights in Showers Induced by Light and

Heavy Primary Masses

The electron size lg(N e ) and muon size lg(N tr µ ) provide an opportunity to separate light from heavy primary CR initiated showers. The size parameters lg(N e ) and lg(N tr µ ), as determined from the KASCADE array, are for the following studies corrected for their shower angle dependent attenuation, see

Ref. [START_REF] Antoni | [END_REF], employing the formula lg(

N 0 i ) = lg(N i (θ))+0.43•(1022/Λ i )•(sec θ- 1)
, where Λ e = 175 g • cm -2 , Λ µ = 823 g • cm -2 , and i stands for e or µ, respectively. The ratio of the corrected parameters, namely lg(N tr,0 µ )/ lg(N 0 e ), turned out to be sensitive to the mass composition of primary CR [START_REF] Antoni | [END_REF][START_REF] Weber | KASCADE Coll., Proc. 25 th ICRC 1997[END_REF].

A ratio lg(N tr,0 µ )/ lg(N 0 e ) of 0.74, which corresponds roughly to Nitrogen, is used to divide the experimental shower data into samples enriched with light and heavy primaries [START_REF] Weber | KASCADE Coll., Proc. 25 th ICRC 1997[END_REF][START_REF] Obenland | [END_REF]. Grouping the showers in the twoparameter space lg(N 0 e ) vs. lg(N tr,0 µ ), according to lg(N tr,0 µ ) > 0.74 lg(N 0 e )

or lg(N tr,0 µ ) ≤ 0.74 lg(N 0 e ), heavy or light primary CR particles are selected, respectively. Using this coarse selection of the showers their muon tracks can 13 be used to calculate the muon production height h µ along the shower axis. Besides the dynamical shift of the radial angle distributions with increasing muon number, the figures show an increasing sensitivity to the mass of the CR particle. The drop to zero yield at small h µ relates to the cuts on ρ and τ angles and the µ-hit shower core distance range employed in the analysis.

At high energies large muon production altitudes are predominantly created by heavy primaries. Narrower distributions, differently for light and heavy CR primaries, indicate a decrease of fluctuations in the < h µ > parameter with increasing primary energy. Fig. 7 compares the muon production height distributions of Fig. 6 to simulation results for proton (triangles) and iron (squares) primaries. The CORSIKA simulations were performed with ver. 6.307 of the code using the QGSJetII and, for interaction energies below 200 GeV, FLUKA2002.4 models [50].

Detector simulations were based on GEANT [START_REF] Brun | GEANT Detector Description and Simulation Tool[END_REF]. Identical cuts for data and simulations were used and the slope of the simulated energy spectrum was weighted to -2.7 below and to -3.1 above the 'Knee'. In Fig. 7 (as well as in Fig. 6) the plots of h µ are normalized to integral yield equal one in the full range up to 12 km, but they are shown to 9 km only in order to expand the low production height region. However, the following discussion of the features seen at high values of h µ is also true for the production heights above 9 km. In this expanded region we see up to ≈ 3.5 km that data are embraced by the simulation results. Muons observed there, up to lg(N µ ) ≈ 5.0, stem from hadronic interactions of energies below 200 GeV, modeled in simulations with FLUKA2002.4 code. At higher lg(N µ ) -the lowest panel in Fig. 7 -contribution from higher interaction energies in this region increases.

This suggests that low-energy interaction model FLUKA describes the data well. Muons produced above ≈ 3.5 km have parent mesons predominantly created at interaction energies larger than 200 GeV, which are modeled in simulations by the QGSJetII code. We observe here that the high-energy model has problems in describing the data.

The comparison reveals more muons at high production heights in the simulations (the distributions are shifted to the right). This excess of muons in the simulations at high altitudes with respect to the data may indicate that muons produced higher up have too high an energy and do not decay, surviving to the observation level, what is not observed in the measurements.

The same effect will occur when the simulated mesons have too small an energy in the region of the first or second interaction, and would decay earlier than in reality. Shifting the maximum of shower development in the models deeper in the atmosphere will act in the direction of reducing the observed discrepancy with the data. suggests that production heights corrected for an appropriate elongation rate will exhibit a clear remaining dependence on the CR particle mass. Mean muon production heights for light and heavy mass enriched showers exhibit a height difference of about 12%, which is 3-4 times the systematic error quoted in Table 2.

Shaded bands show how < h µ > changes when the boundary between 'heavy' and 'light' is moved by ± 0.01 in the lg(N tr,0 µ )/ lg(N 0 e ) ratio, which is 10-20 % of the peak position difference in the distributions of this ratio for proton and iron primary induced showers [START_REF] Obenland | [END_REF]. One observes that the light sample is more sensitive to such a change giving a broader band, what can be explained by larger fluctuations in the shower development than in case of showers initiated by heavy CR primaries.

Mean Muon Production Depth vs. Electron and Muon Shower

Sizes

In the previous section (Figs. 6, 7 and 8) the muon production height was shown to depend on the mass of the CR primary and its energy. In the following we will study this mass dependence of the muon production height in more detail. For this purpose we will use production height values referenced to the top of the atmosphere, called muon production depth, denoted as H µ , and expressed in units of atmospheric depth (g•cm -2 ) using the US-Standard atmosphere [START_REF] Heck | FZKA-Report 6019[END_REF] in data -corrected for temperature and pressureand simulations. Seasonal differences between local atmospheric conditions near the KASCADE site and the US Standard are smaller than ± 2 % (see Fig. 2.1 in Ref. [START_REF] Heck | FZKA-Report 6019[END_REF]) having little influence on our final results.

Usually, shower development is described by the evolution of its electromagnetic component, where X max is understood as the atmospheric depth at which the number of electrons and photons of the air shower reaches its maximum. X max is considered to be a primary mass sensitive parameter.

Concerning muons which stem dominantly from π ± decays, the corresponding production height at which most muons are created may also provide a mass sensitive observable.

To reveal the primary mass sensitivity of the H µ parameter one has to subtract the energy dependence expressed in terms of the elongation rate, taking into account the shower size observables.

Subtracting from the production depth H µ for each track in a shower the energy dependent penetration depth we use the following relation for the

mass dependent H A µ H A µ = H µ -D µ (lg(N tr,0 µ ) -3.6) + D e (lg(N e ) -4.8) , ( 3
)
where the parameter offsets originate from the analysis thresholds.

For the following analysis the elongation rate components (the coefficients in equation ( 3)), similar to the energy dependence of X max observed in EAS simulations, were assigned the following values: Therefore, we may argue that in the corrections of H µ for the elongation rate fluctuations will cancel to some extent and thus, the resulting mass dependent muon production depth H A µ represents a stable observable. In Fig. 9 muon production depth distributions are compared for three N µ size bins before (H µ ) and after (H A µ ) correction for the elongation rate, according to equation ( 3). The borders between different regions are marked for some cases with solid or dashed lines with a certain slope in the lg(N 0 e ) vs. lg(N µ ) plane. In the middle ridge the solid line has the previously employed slope for selecting light or heavy primary particles, being now equal to 0.83 instead of 0.74, due to the transformation to lg(N µ ). For regions away from the middle different slopes of the ridges may be observed, as in the case of the dashed line which marks the ridge of one of the heavier groups of primaries.

D µ =70 g•cm

Using equation (3) the presentation shown in

Data boundaries are curved because of the limited statistics of showers and CR primary types, as shown in Fig. 22 in Ref. [26].

An integral number of muons for a nucleus A induced shower is

N A µ ∼ A(E 0 /A) β , (4)
where β ≈ 0.85 [53]. Taking into account that X max ∼ lg(E 0 /A) + const and assuming that X A max,µ exhibits a similar lg(E 0 /A) dependence as X max , using equation (4), we obtain: e ) vs. lg(N µ ) matrix with effective mean muon production depth < H A µ > along the z-axis. Borders between different regions are marked with solid or dashed lines (see text). Here, lg(N µ ) ≈ lg(N tr,0 µ )+0.5 (see Section 3). Above lg(N µ )=5.5 data from the Grande array in a distance range 250 -360 m are considered and normalized to KASCADE data.

X A max,µ ∼ 1/β[lg(N A µ ) -(1 -β) lg(A)] . (5)
Fig. 9 shows that in the H A µ range of our analysis (> 200 g•cm -2 ), X A max,µ peak position) is larger than the < H A µ > value, due to the tails in H A µ distributions towards the small values. However, this X A max,µ in Fig. 9 is close to the largest < H A µ > value in Fig. 10.

A scale for < H A µ > regions corresponding to 'iron' and 'proton' primaries can be estimated based on an analogy to the electromagnetic cascade, ignoring the difference between X A max,µ and < H A µ >. For a 10 16 eV proton shower we obtain the number of hadronic cascade generations to be n ≈ 5.5 from X p max,µ ≈ 660 g•cm -2 in Fig. 10, having in mind that X p max,µ = n•λ π-air ≈ 5.5

• 120 g•cm -2 ≈ 660 g•cm -2 , and that a fraction of pions decay immediately into muons. For heavy CR particles we read from Fig. 10 a corresponding value X F e max,µ ≈ 540 g•cm -2 ≈ 4.5 • 120 g•cm -2 , leading to 4.5 generations and providing the scale of equation ( 5). For proton primaries the normalization of equation ( 5) is achieved by setting 1/β ≈ 138 g•cm -2 . The other extreme value for equation (5) for iron would lead to X F e max,µ ≈ X p max,µ -29.5

g•cm -2 • ln(56) ≈ 540 g•cm -2 .
The slope of the lines in Fig. 10 can be also interpreted using the muon number vs. energy relation in equations ( 4) and (1). The exponent β can be connected to the amount of inelasticity involved in the processes of A-air collisions [53]. A comparatively steeper slope observed in the dashed line may correspond to an increased inelasticity [53] with respect to the light CR primaries, which may be a subject of further studies. In the present analysis -due to the detection threshold of the MTD -a small fraction of tracks (≈ 15 %) may be missing at small lg(N µ ) (Fig. 11), leading to a bias in the light particle 'mass' interpretation. Unfortunately, our result cannot contribute to the discussion on the behavior of the total CR flux spectrum at higher energies, because the MTD detection area is too small and the amount of collected data limited. The required corrections from Fig. 11 would become prohibitively large. However, as indicated by the dashed line in Fig. 12, the total flux spectrum is in good agreement with the KASCADE spectrum obtained by an unfolding technique [26].

Table 3 provides a collection of the fluxes for different muon production depth H A µ windows showing results for every second bin in Fig. 12. The errors of the fluxes are dominated by systematic errors. detected by the MTD are selected for specific < H A µ > regions (see Ref. [27] or Ref. [START_REF] Bertaina | Proc. 29 th ICRC 2005[END_REF] for comparable information). Data show variations due to the fluctuations in the lg(N 0 e ) parameter. These multiplicities rise with lg(N µ ) but with different slopes for different < H A µ > regions. Therefore, the mean muon multiplicity rise with primary energy can be used for the primary CR mass estimation.

Summary and Outlook

The KASCADE-Grande experimental setup, where the electron and muon components of EAS are measured with high accuracy by arrays of scintillator detectors simultaneously with the precise measurement of muon directions in the tracking detector, allows us to approach the estimation of primary CR composition by means of investigation of the muon production heights in EAS. (Fig. 12). These spectra have features similar to the energy spectra of primary mass groups from Ref. [26]. One can state that the behavior of larger production depths H µ describes the 'lighter' masses and the smaller production depths correspond to 'heavier' ones. A coarse mass scale may be provided by X A max,µ ∼ ln(A).

Muon multiplicity distributions provide an extra parameter in determining the CR primary mass and future tests for high energy interaction models employed in EAS simulations.

Muons are direct messengers of hadronic interactions, and it is a common understanding that they are very well suited for testing the interaction models used in simulations and interpretations of the EAS data. The results obtained in this work are to some extent based on the QGSJet01 model (energy conversion formula (1)) and confirm the effect seen in the spectra in Ref. [26], namely, that the lightest masses (here the largest atmospheric depths) have a break in the spectrum at lower energies compared to the heavier (smaller depths) ones. However, the obtained relative abundances of various production depth ranges (Fig. 12) are model independent.

In addition to the investigation of EAS development using triangulation of muons the mean pseudorapidity of muons in EAS [START_REF] Zabierowski | KASCADE Coll., Proc. 31 st ICRC[END_REF][START_REF] Zabierowski | [END_REF][START_REF] Zabierowski | Proc. 29 th ICRC 2005[END_REF] is a very promising tool for tests of the interaction models. This parameter is closely related to the rapidity of pions -the most abundant products of high energy interactions -and it is highly sensitive to the longitudinal development of the shower in the atmosphere. Simulations show [START_REF] Zabierowski | KASCADE Coll., Proc. 31 st ICRC[END_REF] that muons produced at a certain height H µ carry to the observation level a certain mean pseudorapidity < η >, fairly independent from primary mass and energy. Therefore, H µ provides an average measure of < η > and an excess of certain H µ values can be related to the excess of certain pseudorapidity values in the production of mesons at given heights. KASCADE-Grande with its MTD offers the possibility of such investigation.

The discrepancy in the muon production height distributions between QGSJet simulations and measurements (Fig. 7) points to the necessity of further investigation of high energy interaction models.

Concerning the muon production height, shower zenith angle ranges larger than 18 • will be investigated in the near future. Also very inclined muons studied using the wall modules of the MTD are subject of ongoing investigations. Lines connect the data points to guide the eyes.

Figure captions for the article

  counters registering the electron and muon components in EAS has been enlarged into KASCADE-Grande in 2001 -2003 by adding the extended 'Grande' array of scintillation counters. The MTD was put into operation and started regular data taking in the beginning of 2003. It provides additional information on muons for the two parts of the experiment: KASCADE and Grande. Both arrays trigger the MTD.

Figure 1 :

 1 Figure 1: Anode wire amplitude spectra for individual odd/even numbered wire cells.

Figure 2 :

 2 Figure 2: Illustration of the angular correlations in EAS and the definition of the radial (ρ) and tangential (τ ) angles. The differences between the true and the reconstructed h µ values for the two reconstruction approaches -see text -are also shown.

Figure 3 :

 3 Figure 3: An example of radial angle distributions from Ref. [31] shows the behavior of this quantity in a wide range of lg(N tr µ ), i.e. of primary energy. The selection used comprises a radial distance range 40-120 m, shower zenith angles θ < 30 • and |τ | <0.7 • . Relative yields reflect the CR flux in different lg(N tr µ ) intervals. Lines connect the data points to guide the eyes.

Figure 4 :

 4 Figure4: Mean muon production heights observed at different distances from the shower core in CORSIKA simulated vertical showers for proton and iron primaries and two fixed primary energy values[START_REF] Zabierowski | KASCADE Coll., Proc. 31 st ICRC[END_REF].

( 0

 0 < ρ <8 • and |τ | <0.7 • ) the simulations predict an average scattering angle at the level of 0.2 • . It was checked with Monte Carlo simulations that up to 12 km altitude and a distance range of 40-80 m between a muon and the shower core, a subtraction of the absolute value of τ from ρ or ζ when calculating h µ provides an improvement of the accuracy. Reconstructed negative h µ values are below 5-7 %. The difference between the reconstructed mean value of the positive h µ and the true mean muon production height is illustrated in columns 5 and 6 of Table

Further systematic effects are

  to be considered. In the region of small lg(N tr µ ) in Fig.5, the measured h µ values can have a bias. For a successful reconstruction of the shower parameters on the KASCADE array level, the showers must have a certain electron particle number N e at ground. At small energies only showers which fluctuate to lower production heights have enough particles and are fully reconstructed. Those showers are very likely proton induced showers, and hence a bias may be present in the analysis for small lg(N tr µ ) due to the smaller trigger efficiency for iron, being of about 80% at the lowest lg(N tr µ ) bin.

Figure 5 :

 5 Figure 5: An example of muon production height distributions from Ref. [31] shows the behavior of this quantity in a wide range of lg(N tr µ ), i.e. primary energy. The selection used comprises a radial distance range of 40-120 m and shower zenith angles θ <30 • . Radial and tangential angle cuts are the same as used in this article (ρ=0-8 • and |τ | <0.7 • ). Relative yields reflect the CR flux in different lg(N tr µ ) intervals. Lines connect the data points to guide the eyes.

Figure 6 :

 6 Figure 6: Radial angle distributions (left panel) and muon production height distributions h µ along the shower axis (right panel) for light and heavy CR primary mass enriched showers and for different lg(N tr,0 µ ) intervals. Shower core distances range is 40-80 m, and shower angles θ <18 • . Lines connect the data points to guide the eyes.

Fig. 6

 6 Fig.6shows the distributions of the radial angle and muon production height (left and right panel, respectively) for the selected distance range 40-80 m for different primary energies (expressed in terms of the muon number).

Figure 7 :

 7 Figure 7: Muon production height distributions along the shower axis h µ for light (open circles) and heavy (full circles) enriched shower samples and a shower core distance range 40-80 m, and shower zenith angles θ <18 • compared to CORSIKA simulations employing QGSjetII +FLUKA models for proton (H) and iron (Fe) primaries. Here lg(N µ ) ≈ lg(N tr,0 µ ) + 0.5 (see Section 3). Lines connect the data points to guide the eyes.

Figure 8 :

 8 Figure 8: Experimental values of the mean muon production height along the shower axis < h µ > vs. lg(N tr,0 µ ) for light and heavy primary mass enriched showers and a shower core distance range 40-80 m, and shower angles θ <18 • . The bands bracketing the data points represent variations for the lg(N tr,0 µ )/ lg(N 0 e ) ratio from 0.73 to 0.75.

  -2 and D e =20 g•cm -2 per decade. These values have been varied iteratively by up to ± 20 g•cm -2 to obtain the slope of the middle H A µ ridge in the parameter lg(N 0 e ) vs. lg(N µ ) representation reproducing the lg(N µ )/ lg(N 0 e ) ratio of heavy-light separation (Section 5). We assume that the elongation rates D µand D e are independent of the primary mass.This two-component correction on H µ is similar to the one for lg(E 0 ) in equation (1). There, contributions from lg(N e ) and lg(N µ ), reflecting the complementary information on shower development, are also employed. The shower development leads also to different fluctuations in those shower parameters. It is known from earlier studies that the lg(N e ) parameter exhibits fluctuations to large values, in agreement with the simulations, while the lg(N tr µ ) parameter shows little fluctuations. On the contrary, the H µ parameter (independent of shower angle) is fluctuating towards smaller values.

Fig. 10 hasFigure 9 :

 109 Figure 9: Muon production depth distributions are compared for three lg(N µ )-size bins before (open symbols) H µ and after (full symbols) H A µ correction for the elongation rate, according to the equation (3). The depth cut corresponds to 12 km.

Figure 10 :

 10 Figure 10: lg(N 0 e ) vs. lg(N µ ) matrix with effective mean muon production depth < H A µ > along the z-axis. Borders between different regions are marked with solid or dashed lines (see text). Here, lg(N µ ) ≈ lg(N tr,0 µ )+0.5 (see Section 3). Above lg(N µ )=5.5 data from the Grande array in a distance range 250 -360 m are considered and normalized to KASCADE data.

Figure 11 :

 11 Figure 11: Flux correction factors due to the acceptance and the geometry differences between the array and the MTD (see text for details). The curve represents a polynomial fit.

Figure 12 :

 12 Figure 12: Energy spectra for primaries which produce muons at different effective muon production depth H A µ ; above lg(E 0 /GeV ) = 7.7 the Grande array data analysed in a separate distance range are considered. The dashed line reproduces the CR spectrum as measured by KASCADE -Ref. [26].

Figure 13 :

 13 Figure 13: Differential muon multiplicity spectra in the MTD for different effective muon production depth H A µ . Above lg(N µ ) =5.5 Grande multiplicities are normalized to KAS-CADE multiplicities in the overlapping region.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Anode wire amplitude spectra for individual odd/even numbered wire cells.

Fig. 3 .

 3 Fig.3. An example of radial angle distributions from Ref.[START_REF] Obenland | [END_REF] shows the behavior of this quantity in a wide range of lg(N tr µ ), i.e. of primary energy. The selection used comprises a radial distance range 40-120 m, shower zenith angles θ < 30 • and |τ | <0.7 • . Relative yields reflect the CR flux in different lg(N tr µ ) intervals. Lines connect the data points to guide the eyes.

Fig. 4 .

 4 Fig.4. Mean muon production heights observed at different distances from the shower core in CORSIKA simulated vertical showers for proton and iron primaries and two fixed primary energy values[START_REF] Zabierowski | KASCADE Coll., Proc. 31 st ICRC[END_REF].

Fig. 5 .Fig. 6 .

 56 Fig.5. An example of muon production height distributions from Ref.[START_REF] Obenland | [END_REF] shows the behavior of this quantity in a wide range of lg(N tr µ ), i.e. primary energy. The selection used comprises a radial distance range of 40-120 m and shower zenith angles θ <30 • . Radial and tangential angle cuts are the same as used in this article (ρ=0-8 • and |τ | <0.7 • ). Relative yields reflect the CR flux in different lg(N tr µ ) intervals. Lines connect the data points to guide the eyes.

Fig. 7 .Fig. 8 .Fig. 9 .

 789 Fig. 7. Muon production height distributions along the shower axis h µ for light (open circles) and heavy (full circles) enriched shower samples and a shower core distance range 40-80 m, and shower zenith angles θ <18 • compared to CORSIKA simulations employing QGSjetII +FLUKA models for proton (H) and iron (Fe) primaries. Here lg(N µ ) = lg(N tr,0 µ )+0.5 (see Section

Fig. 10 . 5 (

 105 Fig. 10. lg(N 0 e ) vs. lg(N µ ) matrix with effective mean muon production depth < H A µ > along the z-axis. Borders between different regions are marked with solid or dashed lines (see text). Here, lg(N µ ) ≈ lg(N tr,0 µ ) + 0.5 (see Section 3). Above lg(N µ )=5.5 data from the Grande array in a distance range 250 -360 m are considered and normalized to KASCADE data.

Fig. 11 .

 11 Fig. 11. Flux correction factors due to the acceptance and the geometry differences between the array and the MTD (see text for details). The curve represents a polynomial fit.

Fig. 12 .

 12 Fig. 12. Energy spectra for primaries which produce muons at different effective muon production depth H A µ ; above lg(E 0 ) = 7.7 the Grande array data analysed in a separate distance range are considered. The dashed line reproduces the CR spectrum as measured by KASCADE -Ref. [26].

Fig. 13 .

 13 Fig. 13.Differential muon multiplicity spectra in the MTD for different effective muon production depth H A µ . Above lg(N µ ) =5.5 Grande multiplicities are normalized to KASCADE multiplicities in the overlapping region.

Table 1

 1

	: Variation of track efficiency ( ), Monitor Rate (R) and Gas Composition Stability
	Parameters of the MTD
	d /dT int	+(0.16±0.02) %/ • C
	d /dp int	-(0.12±0.04) %/mb
	dR/dT atm	-(0.31±0.02) %/ • C
	dR/dp atm	-(0.12±0.04) %/mb
	Argon	13.8±0.1%
	CO 2	54.0±0.1%
	Isobutane	30.2±0.1%
	Ethanol	2.0±0.2%
	number lg(N tr µ ), is evaluated in the distance range 40-200 m from the shower
	core. In 'Grande' the muon size is determined in a different way [35]. The
	'truncated' muon number is approximately related to the total muon number
	in the following way: lg(N µ ) ≈ lg(N tr µ ) + 0.5 [36]. This relation, obtained by
	simulations, results for the same showers reconstructed with KASCADE and
	Grande separately in similar muon shower size lg(N µ ) values within 10-15 %.

Table 2 :

 2 Dependence of the systematic mean muon production height errors (in % of the true value) on the reconstruction method

	1	2	3	4	5	6
	Primary E 0 [eV]	∆h ρ µ	∆h ζ µ	∆h ρ µ (cor) ∆h ζ µ (cor)
	H	10 15	-3.1% -4.0%	+2.1%	+2.4%
	Fe	10 15	-7.3% -8.4%	-3.5%	-3.1%
	H	10 16	-1.8% -2.7%	+4.2%	+4.3%
	Fe	10 16	-5.0% -5.7%	-0.5%	+0.02%
	in a wide range of primary energy (expressed by lg(N tr

µ )). Event selection differs from the one used for the analysis in this work in order to show some general features of the distributions. These distributions indicate that the production heights may extend to very large values. For large lg(N tr µ ) values (large primary energies) the deeper shower development is mostly demonstrated by relatively fewer muons from large h µ .

Table 3 :

 3 Fluxes Φ and systematic errors ∆ for different "mass" groups (in (m 2 •sr•s•GeV -1.5 ) -1 ).

	lg(E 0 )[GeV] Φ 620 ∆ 620	Φ 600	∆ 600	Φ 560	∆ 560 Φ 540 ∆ 540
	6.1	551.6 194.1 1641.4 130.5 510.4 324.6	0.0	0.0
	6.3	576.9 133.8 1316.4 109.3 447.5 243.1	3.4	0.0
	6.5	573.6 87.8 1124.7 83.3 375.6 164.9 22.2 11.4
	6.7	312.3 67.1	980.4 106.1 465.7 169.7	6.8	3.6
	6.9	211.7 70.8	669.7	34.7 458.4 105.5 63.5	0.0
	7.1	86.8	22.7	457.7 143.9 471.3	86.4 119.9 80.1
	7.3	72.8	78.6 489.4	10.4 308.5	39.4 157.2 107.6
	7.5	47.7	72.3 131.8	34.2 139.6	72.2 502.8 110.3
	7.7	< 8.7	5.0 < 45.8 31.8 < 47.1 19.1 717.6 194.8
	7.9	< 7.3 16.7 < 48.9	2.2 < 30.9	8.4 496.3 246.0
	8.1	< 4.8 14.6 < 13.2	6.9 < 14.0 14.5 424.1 305.0
	It has been shown that certain muon production depth H A µ regions create
	bands in the two parameter lg(N

e ) vs. lg(N µ ) space, which can be transformed into CR energy spectra in the range from ≈ 10 15 eV to ≈ 10 17 eV

Acknowledgements

The authors would like to thank the members of the engineering and technical staff of the KASCADE-Grande collaboration, who contributed to the success of the experiment. Our special thanks go to all who over past years have contributed to the design and build-up of the MTD. The KASCADE-Grande experiment is supported by the BMBF of Germany, the MIUR and INAF of Italy, the Polish Ministry of Science and Higher Education (this work by grant for the years 2009-2011), PPP-DAAD project for 2009-2010, and by the Romanian Authority for Scientific Research CNCSIS-UEFISCSU (grant PNII-IDEI no.461/2009, code 1442/2008 and project PN 09 37 01 05).


M Ludwig

T Huege

REAS3: Monte Carlo simulations of radio emission from cosmic ray air showers using an "end-point" formalism

Keywords: cosmic rays, extensive air showers, electromagnetic radiation from moving charges, computer modeling and simulation, radio emission, endpoint formalism

In recent years, the freely available Monte Carlo code REAS for modelling radio emission from cosmic ray air showers has evolved to include the full complexity of air shower physics. However, it turned out that in REAS2 and all other time-domain models which calculate the radio emission by superposing the radiation of the single air shower electrons and positrons, the calculation of the emission contributions was not fully consistent. In this article, we present a revised implementation in REAS3, which incorporates the missing radio emission due to the variation of the number of charged particles during the air shower evolution using an "end-point formalism". With the inclusion of these emission contributions, the structure of the simulated radio pulses changes from unipolar to bipolar, and the azimuthal emission pattern becomes nearly symmetric. Remaining asymmetries can be explained by radio emission due to the variation of the net charge excess in air showers, which is automatically taken into account in the new implementation. REAS3 constitutes the first self-consistent time-domain implementation based on single particle emission taking the full complexity of air shower physics into account, and is freely available for all interested users.

Introduction

High energy cosmic rays initiate extensive air showers when entering the Earth's atmosphere. The electromagnetic component of the air shower produces radio emission which can be measured to study the characteristics of the cosmic rays. Radio detector arrays like LOPES [1], [2] and CODALEMA [START_REF] Ardouin | CODALEMA Collaboration[END_REF], [4] have verified that radio emission in the air is dominated by a geomagnetic effect, and they have studied correlations of the radio signal with shower parameters in great detail. With AERA [START_REF] Van Den Berg | Proceedings of the 31 st International Cosmic Ray Conference[END_REF] within the Pierre Auger Observatory [START_REF] Abraham | [END_REF] and LOFAR [7], radio detection will be applied on larger scales and new "super-hybrid" techniques for measuring cosmic ray air showers will be employed. To understand the measurements and to learn more about the physics of cosmic rays using radio signals from air showers, a solid theoretical understanding of the radio emission process is needed. Presently, two major approaches exist, both of which are mainly based on geomagnetic effects [START_REF] Huege | Proc. of the ARENA 2008 conference[END_REF]. On the one hand, there is a model based on the calculation of emission from the deflection of single particles in the geomagnetic field which was mainly investigated by Huege et al. [9], [10], [11], [12], in particular with the implementation in the REAS Monte Carlo code. The treatment of the radio emission with Monte Carlo techniques makes it straight-forward to couple it with detailed Monte Carlo air shower simulations. In case of REAS, this is done with CORSIKA [START_REF] Heck | FZKA Report 6019[END_REF] in which the needed particle distributions are histogrammed by COAST [12]. It is interesting to note that the original frequency-domain calculation of this implementation predicted spectra decaying to zero at small frequencies [9], the time-domain implementation in the simulation code REAS, however, hitherto predicted spectra levelling off at a finite field strength for small frequencies, leading to essentially unipolar pulses [11], [12]. On the other hand, the macroscopic description of geomagnetic radiation (MGMR) developed by Scholten et al. [START_REF] Scholten | [END_REF] constitutes a modern implementation of the approach for radio emission modelling by Kahn and Lerche in 1966 [15]. Due to the Lorentz force, moving electrons and positrons are separated in the Earth's magnetic field. In MGMR, this is described by a time-variable net electric current in the electron-positron plasma, which is moving through the atmosphere with approximately the speed of light. The emission of an electromagnetic pulse is caused by these time-dependent transverse currents. Due to the dependence on the evolution of the charged particles in the shower with time, the MGMR model predicts a bipolar structure of the pulse shape [START_REF] Scholten | [END_REF]. The differences in the results of both models were studied and led to the conclusion that at least one model was not complete. The main reason for the different results of both approaches was that in REAS2 (and all other timedomain approaches based on single particle emission, e.g. [START_REF] Duvernois | Proc. of the 29th ICRC[END_REF], [START_REF] Suprun | [END_REF]), emission due to the variation of the number of charged particles within the shower was not considered. A detailed discussion of this issue can be found in [START_REF] Huege | Proc. of the ARENA conference 2010[END_REF] to which we kindly refer the reader for further details. To correct the implementation in REAS, these missing radiation contributions had to be taken into account. It should be stressed that the resulting implementation, based on the "endpoint" formalism, consistently describes the radiation of the complete underlying particle motion, not just "synchrotron"-like emission from the particle acceleration (see [START_REF] James | [END_REF] for a discussion of the universal nature of the end-point formalism).

Simulation algorithm of REAS

From REAS2 to REAS3, the simulation algorithm is unchanged and the air shower information is provided by CORSIKA and COAST in the same way as before. Therefore, we here only give a short overview of the technical implementation and for details, we kindly refer the interested reader to [12]. First, the shower is simulated with CORSIKA using the air shower parameters of interest (such as primary energy, magnetic field, mass of primary, incoming direction, etc.). Using COAST, the information of the electrons and positrons is saved in histograms. These histograms contain information about the atmospheric depth of the particle, the particle arrival time, the lateral distance of the particle from the shower axis, the particle energy and the particle momentum direction. In the next step, REAS is generating individual electrons and positrons randomly according to the histogrammed distributions. These particles are then tracked analytically in the geomagnetic field. Note that the trajectories of the REAS simulation do not represent real trajectories of particles, i.e. one long particle trajectory is represented by an ensemble of several shorter, unrelated trajectories in the code. The length of the trajectories is determined by a parameter λ which is explained in chapter 5. In REAS3, the trajectory is distributed symmetrically around the position at which the particle is generated.

Contributions due to charge variation

In REAS2 [12], the radiation of single particles of an air shower is calculated as

E( x, t) = e ⎡ ⎣ n -β γ 2 (1 -β • n) 3 R 2 ⎤ ⎦ ret + e c ⎡ ⎣ n × [( n -β) × ˙ β] (1 -β • n) 3 R ⎤ ⎦ ret , ( 1
)
where e indicates the particle charge, β = v(t)/c is given by the particle velocity, R(t) = | R(t)| describes the vector between particle and observer position,

n(t) = R(t)/R(t)
is the line-of-sight direction between particle and observer, and γ is the Lorentz factor of the particle. The index "ret" means that the equation needs to be evaluated in retarded time.

The electrons and positrons emit radiation continuously along their track. To get a consistent description of radio emission in air showers, however, not only contributions due to the deflection of the particles in the magnetic field have to be taken into account, but also contributions due to the variation of the number of charged particles. REAS2 treats radiation processes only along the trajectories, but not at the end or the beginning of trajectories. Strictly speaking, this is equivalent to the situation that the particle arrives with the velocity v ≈ c given by CORSIKA, enters the Earth's magnetic field where it is deflected and describes a short curved track and finally flies out of the influence of the magnetic field still with a velocity v ≈ c. Figure 1 shows a sketch for such a particle trajectory. It is obvious that this is not describing the real situation in an air shower. Consequently, to complement the description in the Monte Carlo code, radiation at the beginning and at the end of each particle trajectory has to be calculated. If at a given atmospheric depth more particle trajectories start than end, i.e. the number of charged particles grows, this results in a net contribution. The same is true if the number of charged particles declines, i.e. more particle tracks end than start. Note that a net contribution occurs as well due to the change of the geometrical distribution caused by the spatial separation of the charged particles. In REAS3, a particle is treated as if it was created at rest and became "instantaneously" accelerated to v ≈ c, flew on a short curved track through the Earth's magnetic field and got decelerated to rest again (see figure 2). The acceleration process at the injection as well as the removal of an electron or positron takes place on time scales short in comparison with the frequencies of interest (ν observed ≤ 100 -1000 MHz), i.e. δt 1 ν observed . Hence, only the time-averaged process is of interest, which will give a discrete contribution in contrast to the continuous emission along the curved particle trajectory. To calculate contributions for the start-point, we consider the far field term of the general radiation equation and integrate over the injection time δt. The "static" term of the radiation formula (1), the velocity field, can be neglected, because the "radiation" term completely dominates the signal for distances R relevant in practical applications. The relation for the retarded time used for transforming the integral from dt to dt is derived from t = t -R(t )/c. With dt = (1β • n)dt the integral is solved as shown in the following calculation:

E( x, t)dt = e c t 1 t 0 n × [( n -β) × ˙ β] (1 -β • n) 3 R ret dt = e cR t 1 t 0 n × [( n -β) × ˙ β] (1 -β • n) 2 dt = e cR t 1 t 0 d dt ⎛ ⎝ n × ( n × β) (1 -β • n) ⎞ ⎠ dt = e cR ⎡ ⎣ 1 • n × ( n × β) (1 -β • n) ⎤ ⎦ t 1 t 0 - e cR t 1 t 0 0dt = e cR ⎛ ⎝ n × ( n × β) (1 -β • n) ⎞ ⎠ (2)
Likewise, one gains the electric field for the end-point of the trajectory:

E( x, t)dt = - e cR ⎛ ⎝ n × ( n × β) (1 -β • n) ⎞ ⎠ (3)
It is important to first transform the integration time into retarded time before calculating the integral. The electric fields of the end-points1 and of the radiation along the track then have to be added to get a self-consistent implementation for modelling radio emission in an air showers.

Continuous vs. discrete calculation and incorporation of endpoint contributions

Adding the discrete endpoint contributions to the continuous contributions along the tracks may produce problems. The radiation associated with the end-points not only contains the emission due to the tangential acceleration. Due to the change in the direction of particle movement between the beginning and the end of the trajectory, radiation associated with the perpendicular acceleration, which was so far treated in the continuous description, is also contained. Combining the two descriptions therefore exhibits a risk of double-counting. In order to avoid such problems, it is preferable to change the calculation for the continuous contributions along the curved particle tracks to a discrete representation. Hence, the chosen representation in REAS3 is completly discrete to ensure that the calculations are self-consistent. To describe Fig. 3. Sketch of the trajectory with a discrete description and end-points.

radio emission contributions along the particle trajectories in a discrete picture, the trajectories of the particles are split in straight track fragments joined by "kinks" in which the velocity of the particles is changing instantaneously. A sketch of this description is given in figure 3. The instantaneous change of velocity at the kinks leads to radiation. With particle velocity β 1 • c before and β 2 • c after the kink, radiation for one kink of the trajectory is:

E( x, t)dt = t 2 t 1 e c n × [( n -β) × ˙ β] (1 -β • n) 3 R ret dt = F (t 2 ) -F (t 1 ) = e cR ⎛ ⎝ n × ( n × β 2 ) (1 -β 2 n) ⎞ ⎠ - e cR ⎛ ⎝ n × ( n × β 1 ) (1 -β 1 n) ⎞ ⎠ (4)
To verify that the continuous and the discrete calculations of emission contributions along the trajectories are equivalent and that both descriptions produce the same results, the REAS2 code was changed to calculate emission using the discrete approach of straight track segments connected with kinks (without the additional contributions of end-points). Analytically, the equivalence of the two approaches for the frequency domain has been shown by Konstantinov et al. [20]. Several tests with the REAS code have confirmed that the implementation of the discrete description is equivalent to the continuous one as can be seen in figure 4. The figure illustrates the equivalence for a vertical proton-induced air shower with primary energy of E p = 10 17 eV. The observer position is 100 m north of the shower core. The advantage of the discrete calculation is the consistency of the description of emission contributions along the tracks and emission at the endpoints which is making the incorporation of radiation at the endpoints canonical. To complement the former implementation in the Monte Carlo code with the emission due to the variation of the number of charged particles, it is therefore convenient to use the discrete description. In the discrete picture, contributions at the beginning and the end of a track are just kinks with β 1 = 0 and β 2 = 0, respectively. This self-consistent emission model has been incorporated in REAS3, taking the radiation at the beginning and the end of a trajectory as well as along the curved trajectory into account. The obtained results for the radio signal is discussed in section 6, but first the numerical stability will be demonstrated.

Numerical stability

As already mentioned in chapter 2, particle trajectories are represented by ensembles of several shorter, unrelated trajectories in the code. This ensures that the phase-space distribution of the analytically propagated particles stays consistent with the underlying particle distribution without having to treat energy losses during the propagation explicitly (cf. section 4 of [12]). The length of the short segments is controlled by a parameter λ. This parameter λ determines the length of the short tracks into which the real trajectories of the particles are divided (it does not denote the sampling density of kinks on a track). If λ is chosen inadequately large, the discrepancies between the particle distributions recreated in REAS and the distributions histogrammed in CORSIKA are getting too large. Thus, to avoid these discrepancies λ has to be chosen small enough. In REAS2, it was recommended to set this parameter to 1 g cm 2 , even though it is just a technical parameter, i.e. the result does not depend on the exact value of λ as long as λ is set small enough. Also, in REAS3 the tracklength of the single short trajectories should have no influence on the physics results as long as it is chosen small enough. This cross check was made for a vertical proton-induced air shower with E p = 10 17 eV. As shown in figure 5 for an observer 100 m north of the shower core, the result converges with decreasing values of λ. The same result is obtained for an observer 400 m north of the shower core which is displayed in the right column of figure 5. The comparison of both observer positions demonstrates that far away from the shower core, the result converges faster than close to the shower core. For λ = 0.1 g cm 2 , a stable result is obtained in both cases. For observer positions close to the shower core as for 100 m distance, it is also visible that for larger values of λ, e.g. for λ = 1.0 g cm 2 , there is already a finite radio contribution for negative times. These unphysical contributions appear if partciles are created with symmetric trajectories around the place of particle generation and λ is chosen inadequately large. For small enough values of λ, however, the symmetric trajectories ensure a faster convergence of the result than the asymmetric trajectories used in REAS2. With smaller parameters of λ, however, the high frequency noise increases as well as the computing time. This first effect is larger for observers far away as can be seen in figure 5. To optimize the calculation in REAS3, the path depth λ of the electron and positron trajectories is therefore chosen as a dynamical parameter depending on the lateral distance of the observer. Therefore, it is recommended to set λ to 0.1 g cm 2 at the core and increase this value linearly every 100 m as it is done by default in REAS3. The advantage of this implementation is on the one hand to gain stable results for all observers and on the other hand to avoid high frequency noise for observers at larger lateral distances. The difference in the recommended value of λ in REAS2 and REAS3 is due to the fact that the emission model in REAS3 requires a precise description of the particle momenta during the propagation, which requires a more fine-grained treatment.

Results

In chapter 5, we have verified that REAS3 is producing stable results and that the endpoints were implemented correctly. In chapter 6.1, we focus on the results obtained with REAS3 and in section 6.2 we discuss the influence of the charge excess in air showers on the radio emission. In both chapters, the radio emission was calculated for a proton-induced vertical air shower with primary energy of E p = 10 17 eV. The shower itself was generated with CORSIKA 6.7 and COAST. The positions of the observers were chosen at sea level with different lateral distances and relative observer orientations to the shower core. For the simulations with REAS2 and REAS3, identical showers were taken. This is easily possible because the histogramming approach allows an easy seperation between the air shower modelling and the radio emission calculation.

REAS2 vs. REAS3

To study the changes introduced by the implementation of emission due to the variation of the number of charged particles, a simple, vertical shower geometry was chosen as specified above. The magnetic field was taken as horizontal with a field strength of 0.23 Gauss to get a geomagnetic angle of 90 • . Figure 6 shows a comparison between unfiltered (i.e. unlimited bandwidth) pulses and frequency spectra of REAS2 and REAS3 for observers with different positions with respect to the shower core. The spectra show the total field strength for two azimuthal observer directions: the thick line displays the total field strength for an observer east of the shower core and the thin line for an observer north of the shower core. It is clear that the strength of the pulses as well as the time structure of the pulses has changed, while the changes in the pulse amplitudes are dependent on the observer azimuth angle. The time structure of the pulses has changed from unipolar to bipolar. This can also be seen in the frequency spectra (lowest row of figure 6) where in case of REAS3, the spectral field strengths drop to zero for frequency zero, as it was the case in the analytical implementation [9]. It can be argued from basic physical arguments that the spectral field strength has to drop to zero for small frequencies because the source of the coherent emission exists only for a finite time in a finite region of space and thus the the zero-frequency component of the emission, which corresponds to an infinite time-scale, can contain no power (cf. [START_REF] Scholten | Proc. of the ARENA 2008 conference[END_REF]). While the emission pattern of REAS2 was azimuthally asymmetric, REAS3 simulations are much more symmetric, as is expected for the given shower geometry. The increased symmetry can also be seen in the frequency spectra where the characteristics of the spectral field strengths with frequency are getting very similar for the two different azimuthal observer positions in case of REAS3. The remaining asymmetry will be discussed in section 6.2. To get a more general impression of the changes in the amplitudes of the signal with REAS3, the lateral dependence of the emission is studied. Figure 7 shows this dependence for the unfiltered, full bandwith amplitudes. In REAS2, there was a much stronger signal for observers in the north or south than for observers in the east or west. In REAS3, the signal pattern is much more symmetric, but obervers in the eastern direction measure a higher absolute field strength and observers in the west see a lower field strength than observers with other azimuthal positions. In general, the amplitudes of the field strength got lower from REAS2 to REAS3, while the changes for observers north and south of the shower core are much larger than for observers in the east or the west. To compare simulations with experimental data (e.g. of LOPES), the REAS simulations have to be filtered to a finite observing bandwidth. This is done by the helper application REASPlot which is included in the REAS3 package. In this paper, REASPlot was used with an idealised rectangular 43 MHz-76 MHz bandpass filter which can lead to acausal contributions at negative times due to the idealisation of the filter. In figure 8, the filtered pulses for observers with lateral distance of 100 m east and north of the shower core are shown.

The quick oscillations are determined by the selected filter bandwidth. The differences in the pulse strengths are much smaller for the filtered than for the raw pulses because the strongest changes between REAS2 and REAS3 occur at low frequencies (Please note that figure 6 is plotted on a log-log scale). In general, e.g. for different geometries, the differences in the amplitudes of the filtered pulses can be larger than in this example (cf. [START_REF] Ludwig | Huege -Proc. of the ARENA conference 2010[END_REF]). Again, the increased azimuthal symmetry of REAS3 is observable. To get an overall impression of the change from REAS2 to REAS3 and the influence of the net charge excess on the east-west asymmetry, we show the contour plots of the 60 MHz field strength in figure 9. The ellipticity which is seen in REAS2 is replaced by a nearly azimuthally symmetric pattern in REAS3. The "clover" pattern seen earlier in the north-south polarisation component is no longer visible. In the vertical polarisation component (not shown here), there is no significant flux for either of the two simulations, as expected for a vertical shower. In the contour plots of REAS3, the east-west asymmetry is visible as well. The resulting REAS3 emission pattern can be interpreted as a superposition of a circularly symmetric contribution with a v × B (in this case thus pure east-west) polarisation and a radially polarised emission contribution caused by the time-varying net charge excess. In summary, the incorporation of radiation due to the variation of the num- ber of charged particles in the form of endpoint contributions results in a clear change from REAS2 to REAS3. The revised, self-consistent implemented model in the Monte Carlo code REAS3 predicts bipolar pulses with a mostly symmetrical emission footprint. In addition to these changes, a further emission contribution arises from the variation of the net charge excess, which explains the remaining asymmetries. This contribution will be discussed in the following section.

Discussion of charge excess emission

The observed east-west asymmetry mentioned in section 6.1 arises from the fact that more electrons than positrons exist in an extensive air shower [START_REF] Bergmann | [END_REF]. This net charge excess of order 10-20% leads to a contribution in the radio signal even in the absence of any magnetic field. Hence, in this section an air shower was simulated with CORSIKA for the geometry and primary characteristics as already mentioned at the beginning of this chapter, but in contrast to section 6.1, the magnetic field strength was set to 0 Gauss. Although the showers calculated with B = 0.23 Gauss and B = 0 Gauss do not represent the exact same particle distributions, we can interpret the predictions for the shower with B = 0 Gauss as the contribution due to the net charge excess of the shower used in section 6.1. This allows us to test whether indeed the eastwest asymmetry can be associated with the pure charge excess. In REAS2, radio emission is produced due to deflection of charged particles in the magnetic field. Hence, there is no radiation without magnetic field. In contrast, REAS3 takes also emission due to the variation of the net charge excess into account. A radially polarised component for emission due to charge excess is expected, as seen also in the macroscopic approach [24]. The radiation pattern for a shower with B = 0 Gauss is indeed radially polarised as illustrated in the contour plots of the 60 MHz field strength in figure 10. For the vertical component (not shown here) there is again no significant flux, as expected. Again, the closest observer position to the shower core is 50 m. Please note that the contour levels for the simulation without magnetic field are smaller than for the simulation with magnetic field, i.e., the relative field strength of the charge excess emission at 60 MHz is small in the distance range up to 200 m. To study the influence of the net charge excess emission on the overall radio signal, it is helpful to look at the polarisation vectors in the plane perpendicular to the shower axis. For the pure geomagnetic emission, the polarisation vectors at all observer positions point in the same direction as illustrated in the left sketch of figure 11. The right sketch illustrates the polarisation vector for the Fig. 11. Sketch of the polarisation vector. Left: uniform pattern as it is the case for pure geomagnetic emission. Right: radial pattern as it is the case for the net charge excess emission.

emission due to the variation of the net charge excess. The direction in which the vector points is changing with the observer position, following a radial pattern. Hence, for an observer in the east the total signal S E is given by S E = S gm + S ce [START_REF] Van Den Berg | Proceedings of the 31 st International Cosmic Ray Conference[END_REF] where S gm is the pure geomagnetic contribution and S ce the net charge excess contribution to the signal. For an observer in the west, the total signal S W is composed of

S W = S gm -S ce (6)
With the signals measured east and west from the shower axis, the signal for the pure geomagnetic emission and for the net charge excess can be calculated by

S gm = 1 2 (S E + S W ) and S ce = 1 2 (S E -S W ). (7)
To verify the assumption made in equations ( 5) and ( 6), we calculated the signal of the charge excess as it is resulting from above for the shower with B = 0.23 Gauss and compared this with the emission of the shower with B = 0 Gauss. Figure 12 illustrates that both pulses match. Therefore, the eastwest asymmetry in the azimuthal emission pattern seen in the figures of section Fig. 12. Comparison of emission due to net charge excess without magnetic field (solid line) and the calculated signal for the net charge excess (dashed line) from a shower with magnetic field. Displayed is in each case the east-west polarisation. Left: 100 m distance from shower core. Right: 400 m from shower core.

6.1 is completely reducible to the emission of the net charge excess in an air shower.

Finally, it is interesting to quantify the relative strength of the charge excess emission with respect to the pure geomagnetic radio emission. Analyses studying the dependence of the radio signal on pure geomagnetic emission have already been done and have shown that there might be discrepancies between a pure geomagnetic model and the measured data [START_REF] Isar | Proc. of the 31st ICRC[END_REF]. The ratio of the net charge excess signal and the pure geomagnetic radiation can be calculated from equation (7) to quantify the relative influence of the net charge excess. Figure 13 illustrates the ratio for the unfiltered full bandwidth amplitudes and the 43 to 76 MHz filtered bandwidth amplitudes for the east-west polarisation component. For the filtered case the ratio in the plot is shown only for observers up to 400 m lateral distance. The reason is that the frequency spectra for a vertical shower drop fast with increasing lateral distance as was also seen in figure 6. Consequently, in the used frequency range the signal is not anymore distinguishable from numerical noise at very large distances. For the amplitudes of the unfiltered pulses the charge excess has more and more influence with larger distances, from a few % close to the core to around 90% of the pure geomagnetic emission at 1200 m. However, for the filtered pulses the ratio is almost constant over the whole range, at a level of ∼10%. It is important to clarify that the emission of the net charge excess The lines may not represent the correct interpolation between the points due to the logarithmic scale of the x-axis. The ratio of the filtered data is only shown for distances up to 400 m, because then the signal is not anymore distinguishable from noise in this frequency range.

occurs due to the variation of the number of charged particles and not due to Cherenkov-like emission. Both processes have been described in the pioneering articles of Askaryan [START_REF] Askaryan | [END_REF], [27], but the term "Askaryan radiation" is today generally interpreted as Cherenkov emission in dense media. For the inclusion of Cherenkov-like emission, a refractive index is needed which is not "unity". REAS3 approximates the index of refraction to be unity so far.

Conclusion

Using an "end-point" formalism, we have implemented a fully consistent Monte Carlo model for radio emission from air showers, which in particular includes also the emission contributions due to the variation of the number of mov-ing charges. For the implementation into the simulation code REAS, the description was changed from a continuous treatment to a discrete description, because in the latter representation endpoint contributions can be added canonically. REAS3 produces stable results. The pulse shape has changed from unipolar pulses in REAS2 to bipolar pulses in REAS3. This is also apparent in the frequency spectra of REAS3 simulations which drop to zero for small frequencies. Several results illustrate the increased symmetry in the azimuthal emission pattern predicted by REAS3. Due to charge excess in extensive air showers it is evident that a small asymmetry has to remain. This effect was shown by calculating the pure geomagnetic and charge excess contributions from the asymmetries present in the full simulations and simulating radio emission in the complete absence of a magnetic field which agree well within numerical uncertainties. In the absence of a magnetic field, REAS3 predicts a radially polarized emission pattern. Due to the presence of this charge excess emission it is also obvious that the radio emission is not purely of geomagnetic origin and thus cannot be described by a pure v × B dependence, but that the signal polarisation depends on the exact observer position relative to the core. For a vertical shower, the relative strength of the emission due to the variation of the net charge excess with respect to the pure geomagnetic emission increases from a few % close to the core up to ∼90% at lateral distances of 1200 m for the unlimited bandwidth pulse amplitudes. This net charge excess affects the east-west symmetry as well as the polarisation of the radio emission. For the filtered 43 MHz-76 MHz bandwidth amplitudes the relative strength of the net charge excess with respect to the pure geomagnetic emission is constant at a level of approximately 10% for different lateral distances. Although the changes introduced with REAS3 are significant, many qualitative results obtained with earlier simulations are still valid (e.g., approximately exponential lateral distributions, a dependence of the lateral slope on X max , the coherent scaling of the pulse amplitudes, and the fall-off of the frequency spectra to higher frequencies and larger distances). Certainly, the absolute amplitudes, pulse shapes and values of scaling parameters have changed with REAS3. With the implementation of emission due to the variation of the number of charged particles in REAS3, for the first time a self-consistent timedomain model exists which takes the full complexity of air shower physics as provided by CORSIKA simulations into account. The source code of REAS3 will be freely available on request. REAS3 therefore constitutes a widely usable radio simulation tool which can be employed to study radio emission from air showers initiated by arbitrary primary particles, using any of the available hadronic interaction models, realistic atmospheric profiles and arbitrary shower geometries (including near-horizontal showers, since the atmosphere in REAS3 is treated with a curved geometry).

Fig. 1 .

 1 Fig. 1. Sketch of the the trajectory how it was implemented in REAS2.

Fig. 2 .

 2 Fig. 2. Sketch of the the trajectory with start-and end-point as needed for a consistent description of radio emission from EAS.

Fig. 4 .

 4 Fig. 4. Comparison of discrete and continuous description of radio emission (without endpoint contributions) for an observer 100 m north of shower core of a vertical 10 17 eV air shower. The results are identical.

Fig. 5 .

 5 Fig. 5. Influence of different values of the technical parameter λ. Left: observer 100 m north of shower core. Right: observer 400 m north of shower core.

Fig. 6 .

 6 Fig.6. Upper row: raw pulses for observers 100 m (left) and 400 m (right) north of shower core. Middle row: raw pulses for observers 100 m (left) and 400 m (right) east of shower core. Lowest row: frequency spectra for observers 100 m (left) and 400 m (right) east and north of shower core. For the raw pulses the east-west polarisation of the electric field is shown, whereas for the frequency spectra the total spectral field is shown.

Fig. 7 .

 7 Fig. 7. Lateral dependence with full bandwidth amplitudes. Left: REAS2. Right: REAS3.

Fig. 8 .

 8 Fig.8. Filtered pulses (a simple rectangular 43-76 MHz bandpass filter is used) for observers with 100 m distance to the shower core. Left: the observer is north of the shower core. Right: the observer is east of the shower core.

Fig. 9 .

 9 Fig. 9. Contour plots of the 60 MHz field strength for emission from a 10 17 eV vertical air shower. From left to right: total field strength, north-south and east-west polarisation component. Contour levels are 0.1 µVm -1 MHz -1 apart. The closest position of the simulated observers to the shower core is 50 m. Upper row: REAS2. Lower row: REAS3

Fig. 10 .

 10 Fig. 10. Contour plots of the 60 MHz field strength for emission from a 10 17 eV vertical air shower without any magnetic field. Contour levels are 0.03 µVm -1 MHz -1 apart. The closest observer position to the shower core is 50 m.From left to right: total field strength, north-south and east-west polarisation component. The "spike" in the lower-left part of the contours is associated with noise in the simulation.

Fig. 13 .

 13 Fig.[START_REF] Heck | FZKA Report 6019[END_REF]. Comparsion of the charge excess and the pure geomagnetic contribution on the radio signal of an vertical air shower for the east-west polarisation component. The lines may not represent the correct interpolation between the points due to the logarithmic scale of the x-axis. The ratio of the filtered data is only shown for distances up to 400 m, because then the signal is not anymore distinguishable from noise in this frequency range.

In the following we use the term "end-points" as a general term for both, startpoints and end-points, as both are treated in the exact same way.

Acknowledgements

The authors would like to thank S. Buitink, R. Engel, H. Falcke, O. Scholten and K. D. de Vries for very helpful discussions and A. Haungs for very useful comments to the manuscript. This research has been supported by grant number VH-NG-413 of the Helmholtz Association.