Full Simulation of Space-Based Extensive Air Showers Detectors with ESAF

Full Simulation of Space-Based Extensive Air Showers Detectors with ESAF


Future detection of Extensive Air Showers (EAS) produced by Ultra High Energy Cosmic Particles by means of space based fluorescence telescopes will open a new window on the universe and allow cosmic ray and neutrino astronomy at a level that is virtually impossible for ground based detectors. In the context of the Extreme Universe Space Observatory (EUSO) project, an end-to-end simulation of EAS observation with a spatial detector has been designed (EUSO Simulation and Analysis Framework, ESAF). This paper describes the detailed Monte-Carlo developed to simulate all the physical processes involved in the fluorescence detection technique, from the EAS development to the instrument response. Particular emphasis is given to modeling the light propagation in the atmosphere and the effect of clouds. The simulation is used to assess the performances of EAS spatial detection. Main results on energy threshold and resolution, direction resolution and determination are reported. Results are based on EUSO telescope design, but are also extended to larger and more sensitive detectors.

cosmic rays; space detectors; extensive air showers; atmosphere models;
PACS 95.55Vj 95.55Fw

[Grenoble]C. Berat, [FirenzeINFN]S. Bottai, [Genova,GenovaINFN]D. De Marco, [Grenoble]S. Moreggia, [GenovaINFN,FirenzeINFN,Dubna]D. Naumov, [Genova,GenovaINFN]M. Pallavicini, [Genova,GenovaINFN]R. Pesce, [Genova,GenovaINFN]A. Petrolini, [Grenoble]A. Stutz, [Firenze,FirenzeINFN]E. Taddei, [Genova,GenovaINFN]A. Thea

Corresponding author. E-mail: Roberto.Pesce@ge.infn.itCorresponding author. E-mail: Alessandro.Thea@cern.ch

[AleAt]Now at ETH, Zürich [DanielAt]Now at Bartol Research Institute, University of Delaware, Newark - DE 19716, USA

1 Introduction

Ultra-High Energy Cosmic Particles (UHECP) with energies in excess of eV hit the Earth with a very low flux which falls down to less of one particle km sr century for particles with energy E 10 eV [1].

The observation of UHECP and the interpretation of the related phenomenology is one of the most challenging topics of modern High-Energy Astro-Particle Physics. Direct detection is impossible at these energies, due to the exceedingly low flux, but UHECP can be detected by observing the Extensive Air Showers (EAS) produced by the interaction of the primary particle with the Earth atmosphere (see [2, 3, 4], and references therein, for recent reviews on these topics).

The ground-based Pierre Auger Observatory (PAO) [5, 6, 7, 8] is currently taking data: its south site (surface of 3000 km) in Argentina and its forthcoming north site in the US (20000 km) will provide, in the next few years, a clear understanding of many important topics [1, 9, 10, 11, 12]. However, a next generation of experiments with sensitivity to UHECP flux well above the PAO one at the highest energies is essential but it seems to be difficult to realize by extending the ground-based technique at ever larger arrays.

Ground-based experiments detect EAS either by sampling charged particles at ground by means of suitable counters and/or by detecting the fluorescence light emitted by the atmospheric nitrogen when excited by the charged particles in the shower. While the former technique is intrinsically ground based, the latter is not, being the fluorescence light emission isotropic. However, fluorescence light detection is limited to moonless nights and is affected by weather conditions.

In 1979 J. Linsley [13] proposed for the first time the idea of detecting EAS fluorescence light with a telescope in orbit around the Earth. While this approach is technically, programmatically and economically challenging, it is one of the most promising that might allow a big step forward in this field, opening a new window on the universe. A space based experiment might observe EAS on an effective surface of a few times 10 km, increasing the event statistics by more than an order of magnitude compared to existing and running ground based experiments, even when taking into account that the duty cycle for observing EAS with fluorescence is limited to clear and moonless nights. Moreover, the huge target mass of the observed air volume (> 10 kg), three orders of magnitude more than the current or foreseen projects in the South Pole ice [14] or mediterranean sea [15] ) may open the era of high energy neutrino astronomy, at least above energies of the order of 10 eV or so. However, we do not discuss high energy neutrino detection in this paper.

The relevance of this innovative approach is widely recognized by the scientific community. The Extreme Universe Space Observatory (EUSO) project [16], originally proposed in 2000 as response to F2/F3 European Space Agency (ESA) call for mission successfully completed the phase A conceptual study. Nevertheless it did not continue only for programmatic reasons related to International Space Station and Space Shuttle availability. The inclusion of UHECP physics in the ESA Cosmic Vision 2015-2025 program [17] provides a suitable framework for the study of future space missions. An introduction to the basic issues concerning the detection of UHECP with a space detector can be found in [18].

In this paper we therefore report the results obtained in the context of the EUSO collaboration by means of a full simulation and reconstruction software code named ESAF, standing for EUSO Simulation and Analysis Framework.

Particularly, we report about results concerning energy threshold as function of UHECP direction and detector size and geometry, direction and energy resolution, the dependence on atmospheric conditions and clouds.

This paper is structured as follows: section 2 presents the ESAF design; section 3 describes the simulation part, giving details on each modules, from the EAS development to the detector response; section 4 contains a description of the reconstruction tasks; section 5 presents several results obtained with the ESAF code; the most relevant features are summarized in the conclusions.

2 The ESAF design

ESAF handles the whole simulation and reconstruction chain from primary particle interaction in the Earth atmosphere until the final reconstruction of the event. The code includes:

  • extensive air shower simulation both by means of internally developed algorithms (e.g. Greizen-Ilina-Linsley (GIL) parametrization, [19]) and with interfaces to existing widely used codes (e.g. CORSIKA, [20]);

  • a complete description of the atmosphere, including aerosols111Aerosols are airborne liquid droplets or solid particles, which may exist at low altitude. and clouds

  • fluorescence and Cherenkov light production

  • a complete simulation of photon propagation, from production point up to the telescope, including detailed interactions with ground and atmosphere, and a Monte-Carlo code dealing with multiple scattering

  • telescope optics simulation

  • telescope geometry

  • telescope photodetector simulation

  • electronics and trigger simulation

  • background simulation

  • pattern recognition and shower signal identification above background

  • reconstruction of direction, energy, and the so called slant depth of the shower maximum (), i.e. the distance (expressed in  of traversed material) of the shower maximum from the top of the atmosphere.

Although the simulation was done having in mind the specific EUSO design, the flexibility of the ESAF code allowed us to obtain results that are quite general and might be useful for any future project.

Figure 1: The upper level structure of ESAF. The simulation and reconstruction modules are detailed in Fig. 2 and Fig. 3 respectively.
Figure 2: ESAF simulation structure.
Figure 3: ESAF reconstruction structure.

2.1 General structure of the ESAF code

The ESAF code is written mostly in C++ with a few external Fortran libraries. A complete object-oriented design is used to allow flexibility, modularity and maintenability of the code. The code is divided in two independent sections, simulation and reconstruction, that share a common infrastructure and a set of common libraries. The code is not public, but is available upon request to interested users222http://www.ge.infn.it/pesce/esaf.html. The main structure of the program and the interactions between the two main sections are shown in Fig. 1.

The simulation is structured in 6 main sections (see Fig. 2): shower simulation starting from the expected UHECP flux, fluorescence and Cherenkov photons production in atmosphere, light propagation from the production point up to the orbiting telescope with simulation of occurring photon interactions during the propagation, simulation of telescope optics (lenses and focal surfaces), simulation of front end electronics and trigger. In this paper we do not cover implementation details that can be found in [21, 22, 23].

The output of the simulation module is a ROOT file [24] with a structure that contains the same information that is expected from real data plus a set of "Monte Carlo truth" data. The level of detail is user-configurable, to optimize the output size to the goal of the simulation (i.e. physics, optics, trigger, etc.). To provide the user with an easy access to the data, a set of detector event viewers is provided with ESAF. They allow a coherent display of the run informations, histograms and several 2D and 3D views of the event development in the detector and on the focal surface (some example plots are in section 3.4).

The reconstruction module reads back the ROOT file representing the fully simulated shower with background and reconstructs the events as if they were real data events from the telescope. Several modules allow pattern recognition and identification of EAS signal above incoherent star light and moonlight background, 3D direction reconstruction, energy and reconstruction. The reconstruction structure is shown in Fig. 3.

The reconstruction framework is the main structure that acquires the data (real or simulated) from an input module and builds the chain of modules necessary to reconstruct the events. The event reconstruction is divided into different tasks (e.g. background subtraction, direction reconstruction, ) and for each task different modules may be available. The use of such a structure allows to easily replace or exclude a given module. In this way different algorithms or algorithm combinations can be tested and compared on the same events. As in the simulation part, the reconstruction modules are configurable by the user. The event container stores all the relevant informations about the reconstructed event. Every module can access the event container, reads the existing data stored in it and writes its own results.

2.2 Atmosphere description

In the simulation process, a description of the atmosphere is required for proper simulation of shower development, light production and photons transfer toward the telescope entrance pupil (see Fig. 2). The atmosphere parameters are also required during the reconstruction process.

In ESAF, the atmosphere and Earth surface are spherically shaped, and several models have been introduced, to form database of atmospheric profiles: the US-Standard model and its complements [25] and the NRLMSIS-00 (MSISE) empirical model [26]. The formers provide temperature, pressure, refractive index and density profiles, up to 100 km, as well as ozone and water vapour profiles. Available data are mean values over time of atmospheric data at different latitudes. With MSISE model, atmospheric profiles are provided for given latitude, longitude, day of the year and hour.

Other crucial data concerning atmosphere description have been implemented in ESAF. The presence of clouds and aerosols has an impact on the UHECP detection from space. The cloud and aerosol effects have to be evaluated with the simulation. The cloud mean fraction in the field of view of an orbiting telescope is around 70%. To describe cloud main characteristics, informations collected with the satellite TOVS (TIROS Operational Vertical Sonder) Path-B [27] are used. To describe aerosols ESAF is also supplied with several atmospheric models and parameters which exist in the LOWTRAN7333LOWTRAN7 is a standard program used in the atmosphere community, to compute transmission in atmosphere. software [28].

A set of dedicated configuration files allow to choose one particular atmosphere model and status among all the possibilities available in ESAF, before runs. The atmospheric parameters are then used at each step where they are required during the simulation or the reconstruction process.

3 The simulation code

3.1 Shower simulation

An accurate and physically reliable simulation of EAS is fundamental for the correct understanding of the UHECP phenomenology and the possible experimental outcomes. During the last 20 years several authors have dealt with the very demanding task of preparing detailed UHECP simulation programs. Particularly, CORSIKA [20] and AIRES [29] , are at present widely used and represent the ’state of the art’ in atmospheric shower simulation. Another shower simulator, UNISIM [30], was developed by some ESAF developers, and has the unique feature of simulating neutrino interactions in atmosphere.

The basic ESAF approach, therefore, does not uniquely consist in the development of brand new shower simulators, but also in having the capability to handle external shower simulation programs. This attitude is realized by implementing in ESAF a generic interface by means of a ’Shower Track’ object which contains the structure for an extremely detailed description of a single shower event, a description that is even more accurate than any of the output of the present existing simulation programs. Using some specific modules, the user can easily write an interface which plugs into the ’Shower Track’ the shower details event by event.

The ’Shower Track’ can manage all typical variables (number of charged particles, number of electrons, energy loss..) representing the shower during its longitudinal development. The ’Shower Track’ has a block of general information and a list of ’Shower Step’ objects along the EAS longitudinal development.

Each ’Shower Step’ object contains an information about the shower longitudinal development corresponding to the track length , which may also vary from event to event. For each step there is also the possibility to handle histograms representing the local distributions of particles in the shower according to their energy, distance from the axis, direction in space, arrival time, and taking care of preserving some of the most important correlation between couples of such physical quantities.

At present it exists inside ESAF an interface to UNISIM, a preliminary version of interface to CORSIKA, and an interface to CONEX [31, 32].

In addition to external program interfaces, ESAF can also fill event by event the longitudinal description inside the ’shower object’ using several parametrizations (Greizen [33], Greizen-Ilina-Linsley (GIL) [19], Gaisser-Hillas [34]) tuned on data and Monte - Carlo outcomes. All these parametrizations need the depth of the first interaction point (). In ESAF it can be fixed or sampled according to Nuclei-Air cross section data [35] or randomly. Such fast algorithms are very useful during the phase of detector optimization where the exact representation of the shower to shower fluctuations is less important.

When some important information is not present in a given external simulation program, ESAF can decide either to skip or to replace using some of the predefined parametrizations implemented in the code. Several different parametrizations are available in ESAF for energy distribution of particles inside the shower [36, 37, 38], particles angular distribution [39], and particles lateral distribution [33, 40] as a function of the shower age.

Most of the results presented in this paper are achieved using the GIL parametrization for the longitudinal development of the number of charged particles , as a function of the shower age :


where is the energy of the primary particle, the atomic number of the primary, the depth considered, is the depth of the first interaction, (air radiation length), MeV (critical energy), , and GeV. These values are chosen to match CORSIKA-QGSJET results [41].

3.2 Cherenkov and fluorescence light simulation

Relativistic charged particles (mainly electrons and positrons) in extended air showers are at the origin of both fluorescence and Cherenkov light emissions. Generation of fluorescence and Cherenkov photons has been included in ESAF. The shower simulation module provides for each ’Shower Step’ the number of electrons , the electron energy distribution and the lateral distributions. The light generation module uses them to compute bunches of photons, as explained in section 3.2.3. This approach is dictated by the huge number of photons produced which makes impossible to follow the fate of each individual photon.

3.2.1 Fluorescence light simulation

The number of photons created along by the EAS electrons is:


where is the normalized electron energy spectrum, is the fluorescence yield (in photons/m) at wavelength for an electron of energy . and are the local pressure and temperature considered at the middle of the track segment.

The fluorescence yield dependence on the local temperature and density conditions is considered according to the empirical relation :

where is the atmosphere density, in kgm, and is in K. The experimental values of the two coefficients and are provided with yield measurements from Nagano measurements [43] or Kakimoto ones [44]. In the Kakimoto case, the fluorescence spectrum is completed with data from Bunner PhD thesis [45]. In ESAF, the choice to carry out simulations with one of these two fluorescence yield measurement data sets is provided by the dedicated configuration file.

The yield is assumed to be proportional to the energy loss ( is the slant depth, in gcm) in the atmosphere, thus is given by:


In ESAF, the Berger-Seltzer [46] formula is used to compute the energy loss directly from .
Fluorescence yield at a fixed energy is determined at  K and  hPa either from [43] or from [44, 45], according to the chosen configuration, with respectively  MeV and  MeV. Examples of computed fluorescence yield spectra are presented on Fig. 4.

Figure 4: Fluorescence emission spectra at  hPa and C computed using Nagano measurements [43] (black line histogram) and Kakimoto ones [44] completed with Bunner rays [45] (red filled histogram). The main fluorescence rays are between 330 and 400 nm, and the three most intense rays are at 337, 357 and 391 nm.

In ESAF the normalized electron energy spectrum is provided by the shower module (see section 3.1). Using equation (4), labeling the integrated fluorescence yield over the wavelength spectrum, and computing the averaged electron energy loss,

we can derive from (2) the following formula implemented in ESAF:

Figure 5: Fluorescence yield integrated over the emission spectrum as a function of the altitude (dashed blue line: from [43] photons/m, continuous red line: from [44], photons/m). and profiles are obtained from the US Standard atmosphere model (see section 2.2). The mean altitude of the maximum of shower development is indicated for 3 shower zenith angles. The weak dependence of the yield with altitude is observed.

The resulting total fluorescence yield (Fig. 5) varies weakly along the EAS development (the mean yield is between 4.1 and 4.7 photons/m). The photons longitudinal distribution along the shower axis is therefore similar to the electrons longitudinal profile. The total number of fluorescence photons isotropically emitted along the shower track is in the order of 10 for a primary energy of 10 eV. This number depends on the shower zenith angle, since inclined EAS starting at higher altitude have a longer path in the atmosphere and thus produce more photons compared to vertical EAS.

3.2.2 Cherenkov light simulation

The algorithm to simulate Cherenkov photon generation is similar to the fluorescence one. The number of Cherenkov photons originated by electrons along a shower track segment is:


where is the Cherenkov production energy threshold. The term describes the normalized energy spectrum of the electrons at the emission point. is the total Cherenkov yield integrated over . The atmosphere refractive index verifies ; it varies very weakly with in the UV range of interest, thus only variation with atmosphere density is taken into account. Therefore Cherenkov yield can be computed in ESAF using the formula:


where is the fine structure constant, and . Typically  nm and  nm; these values can be configured. It has to be noticed that Cherenkov yield depends on local atmosphere parameters through the factor present in (3.2.2) and in the energy threshold expression. The normalized energy spectrum of the electrons at the emission point is provided by the shower generation module.

Then, an effective yield corresponding to the integral of over the electron energy spectrum is computed at each step of length and used to compute the number of Cherenkov photons produced along the step :

Figure 6: Cherenkov effective yield as a function of the altitude, for 3 values of the shower age .

Contrary to fluorescence, Cherenkov yield decreases rapidly with altitude (Fig. 6). Nevertheless, the number of Cherenkov photons produced all along the shower track ( at 10 eV) is of the same order as the fluorescence one. It varies also with shower zenith angle (Fig. 7). Cherenkov intensity reaches a maximum for shower zenith angle around . In case of larger zenith angle, development of showers occurs at higher altitudes, where the Cherenkov yield is weaker. Quasi-vertical showers may reach Earth ground before their development ends, and Cherenkov photon production is not complete.

Figure 7: Number of Cherenkov photons produced by an EAS (proton,  eV) as a function of the shower zenith angle .

The Cherenkov emission angle , as the yield and the energy threshold, is related to the atmospheric refractive index and therefore correlated with altitude. In the atmosphere and for shower relativistic electrons, . In order to obtain the angular distribution of the emitted Cherenkov photons, one has to convolute shower electron angular distribution with . In ESAF, the formula established by Elbert et al. [47] is used:

where is a normalization constant such has . Photons emitted at high altitude are beamed close to the shower axis: at 12 km, , at 4 km, .

3.2.3 Bunches of photons

The number of photons produced by a typical shower is huge, of the order of 10. This fact forbids to follow a direct approach, and simulate the fate of each individual photon. For this reason, the light simulation is done by introducing the concept of bunch, i.e. an object dedicated to EAS photon description.

Bunches of photons are created along the EAS axis: the shower longitudinal distribution is split in steps of length , set by the user in the steering configuration file; typical values are between 1 and 10 gcm. At each step, two bunches are formed, one for the fluorescence emission, the other for the Cherenkov one, according to equations (5) and (9) respectively. Each bunch is characterized by a mean position value, a mean direction according to the computed angular distribution, and a creation time. Specific distributions are associated to each bunch of photons, according to the emission type: wavelength spectrum, longitudinal distribution inside the step, lateral and angular distributions (w.r.t the shower axis). So, when a photon has to be randomly sampled from a bunch to be propagated, its characteristics can be retrieved from each distribution associated to the bunch.

3.3 Photon propagation from EAS to the space-based detector

As mentioned above, more than 10 UV photons are produced by a 10 eV EAS. Before reaching the telescope, these photons may interact with atmosphere constituents or Earth ground: the physics processes to be considered are absorption and scattering. The tracking of such an amount of photons through the atmosphere is obviously impossible with “full” Monte-Carlo algorithms. In ESAF, two different approaches have been designed, each of them being adapted to specific simulation needs. In the following, the different types of interactions which can interfere with photon transfer towards the space-based detector are outlined, as well as their implementation in ESAF. Then the two propagation algorithms are described.

3.3.1 Scattering and absorption processes

Photons emitted in the direction of the telescope, and that are neither absorbed nor scattered during their path arrive at the entrance pupil. To evaluate this fraction of light straightly transmitted, the transmission is evaluated from the interaction cross section (or from associated quantities such as the optical depth and the interaction length) related to the different extinction processes, and taking into account each photon characteristics (wavelength, position, …). EAS photons that are not emitted in the direction of the detector may be scattered later in this direction. To quantify this effect, one needs to know as well the angular distributions of the scattering processes, represented by their phase function. To simulate the radiative transfer, both optical depths and phase functions should be computed in ESAF.

When the scattering particle size is much smaller than the photon wavelength, Rayleigh theory is used, as for scattering on air molecules. When the scattering particle size is of the same order or larger than the photon wavelength, the scattering is described by Mie theory [48]. In the case of a space-based detector, this latter type of description is required to treat photon propagation in clouds and in aerosols.

Rayleigh scattering

In ESAF, the optical depth at a given wavelength (in nm) is computed with the formula used in LOWTRAN7:

where is the atmospheric slant depth along the photon path. The dependence, at the first order, of Rayleigh cross-section is dominant here. The term in is partly due to the dependence of the cross-section on the air index. Studies made with ESAF show that vertical transmission between ground and top of atmosphere is around 50% at , 30% at 300 nm, and 80% at 450 nm. Simulations have been performed with several different density profiles: results show that air density profile variations have very small impact on associated vertical transmission value ().

The normalized phase function of Rayleigh scattering is given by:

Rayleigh scattering is independent of the azimuth angle and almost isotropic: the phase function varies by only a factor of two between 0 and 180, that is very weak compared to phase functions associated to other scattering processes (see below, and Fig. 8).

Mie scattering

Below an altitude of around 8 km, clouds are made of tiny water droplets. Due to the spherical shape of droplets, the original Mie theory can be applied, provided knowledge of the distribution of the scattering droplets radius . At a given , determines the phase function. In UV range, clouds do not absorb photons, but scatter them. Scattering can be considered as independent of wavelength in the range 300-450 nm [45]. The same behaviour is observed for high altitude cirrus clouds (cirrostratus) composed of ice crystals. In ESAF a cloud is described as a homogeneous layer, with a finite vertical height (TOVS cloud). Cloud databases provide integrated vertical depth . Laying from 1.2 to 3.3, the mean value of is around 2. For , the transmission of a vertical UV light beam crossing the cloud layer is %.

Phase function values for low altitude clouds, obtained from one particular cumulus model [49], are stored in tables for  nm, as well as those for cirrostratus, for  nm [50]. Contrary to Rayleigh scattering, Mie scattering in clouds is highly anisotropic, and very forward peaked. Phase functions are compared in Fig. 8.

Figure 8: Cloud phase functions (cumulus and cirrostratus) compared to Rayleigh one. is equal to .

Concerning aerosols, the three main models (rural, urban, maritime) provided by LOWTRAN7 have been introduced in ESAF. These models provide the dependence of the extinction coefficient with altitude and wavelength, specifying the respective part of scattering and absorption processes. These aerosol layers are located between 0 and 2 km altitude. More details can be found in [23] (in french) or [51]. The aerosol layer opacity is specified via the visibility. The smaller is the visibility, the greater is the extinction coefficient. Default visibility values in LOWTRAN are 5 km and 23 km, corresponding to a vertical transmission of 75% and 20% respectively. To describe scattering angular distribution, ESAF is supplied with phase functions from LOWTRAN7 for the three mentioned models.

Ozone absorption

In the UV range of interest, the molecular absorption in atmosphere is dominated by ozone (O) absorption. At a fixed altitude , the ozone absorption length is expressed by

Ozone absorption cross section () data provided by LOWTRAN have been introduced in ESAF. The dependence of cross-section on local temperature is taken into account, entailing a dependence on altitude. is determined from available ozone profiles, and air density. Ozone exists mainly between 10 and 30 km, with a maximum absorption around 20 km. It absorbs the UV light below 330 nm, therefore it has a weak impact on fluorescence transmission, since the three main fluorescence lines are above the wavelengths affected by ozone; it affects mainly Cherenkov light.

Ground reflection

The reflection on the Earth surface is a crucial point for Cherenkov signal detection. The Earth albedo, in UV range, has been measured by ULTRA [52] experiment and is less than 10% [53]. Moreover, reflected photon angular distribution depends on the reflection type, which can be Lambertian (reflected light is scattered isotropically by a Lambertian surface) or specular (highly anisotropic, as on the ocean surface). Both types have been implemented in ESAF. To describe specular reflection, a simple model, reproducing qualitatively air-born measurements of the ocean BRDF (Bidirectional Reflectance Distribution Function[54] with dependence on wind speed, has been introduced. For a wind speed of 9 m/s, specular peak amplitude is 5 times lower than for a wind speed of 2 m/s. If specular reflection is dominant, reflected light intensity is not isotropic in , and can be zero on a large zenith angle range (see Fig. 9).

Figure 9: Simple model for light reflection on ocean. The zenith angle of the incident light beam is 20. The reflectance is computed in the plane containing the incident direction, as a function of the reflected light zenith angle , in case of both lambertian or specular reflection (for 2 wind speed values ).

3.3.2 Propagation algorithms

Below are described two original propagation algorithms designed and implemented in ESAF.

“Bunch” algorithm

The algorithm has been designed to handle mainly three components of the light signal received by the telescope. The first component, called ’Direct Fluorescence’, is the fluorescence light reaching directly the lens without undergoing any interaction in the atmosphere. The second component, called ’Reflected Cherenkov’, is the ground reflected Cherenkov light, having interacted only once along its path from emission point to the telescope. The third component, called ’Air Scattered Cherenkov’, represents Cherenkov photons which have been scattered only once, on air molecules.

The extreme low value of the detector solid angle ( for EUSO telescope shipped on ISS) is of major importance: when considering photons emitted or scattered within the detector solid angle, ’single’ photons can be sampled from the related bunches, and be transmitted one by one to the telescope.

In the case of fluorescence photons directly emitted towards the detector, single photons are sampled from bunches created at fluorescence emission point. Then, they are distributed in time, position and wavelength, according to distributions associated to the bunch they originate. For each photon, the probability that it is transmitted up to the telescope is computed taking into account all possible interactions in atmosphere.

Concerning the two other components, they are both simulated by the following algorithm. Since Cherenkov radiation is beamed in the forward direction, it is possible to propagate each Cherenkov photon bunch step by step, until it reaches ground or exits atmosphere. At each step, two operations are performed:

  1. The number of Cherenkov photons scattered once on air molecules, and directed toward the telescope after the interaction, is computed. A corresponding number of single photons are therefore created. Their final transmission to the telescope from their scattering position is then evaluated : if they undergo another interaction along their path to the detector, they are lost.

  2. The wavelength spectrum of the propagated bunch is convoluted by transmission values computed along the step. Then, the possibility to pursue the bunch propagation is evaluated.

When a Cherenkov bunch reaches the ground, the number of photons reflected toward the detector is computed. A corresponding number of single photons are created and their transmission to the telescope is evaluated : if they undergo another interaction, they are lost.

In presence of aerosols or clouds, the bunch algorithm takes into account their effects on the transmission values, considering the photons scattered on these media as lost.

The bunch algorithm is useful to perform fast simulations (see Fig. 10). But such an algorithm can not handle properly single scattering of fluorescence photons, and multiple scattering of both fluorescence and Cherenkov photons. A mandatory point in EAS analysis is to evaluate the detailed contributions of scattered photons to the detected signal, as this component may introduce systematic uncertainties if not taken into account during reconstruction procedure. To achieve this, a more sophisticated algorithm has been designed to deal with single and multiple scattering of the light produced by EAS.

Figure 10: Examples of time photon distributions, for a shower with  eV, , and a footprint at the telescope nadir - The different contributions are superimposed: fluorescence photons transmitted directly to the telescope, Cherenkov photons scattered once by the atmosphere molecules, and Cherenkov photons reflected by the ground in case of clear sky (top panel) or scattered by the cloud in case of cloudy sky (bottom panel).
“Reduced” Monte-Carlo algorithm

First of all, the simulation of ’Direct Fluorescence’ component does not raise any problem, and it is treated in the same way as presented in the previous section. Then, the aim is to simulate fluorescence and Cherenkov light propagation by performing a Monte-Carlo ray-tracing algorithm, propagating photons one by one. Let us consider photons that reach the telescope after exactly scattering processes. To simulate the transfer of each of them, an ideal algorithm has to follow the steps described below, as illustrated on Fig. 11.

  1. Photon position, direction and wavelength have to be randomly generated, taking into account related distributions defined during light production process (see section 3.2).

  2. Interaction type and position are randomly determined according to the different physical processes able to alter photon path. If the photon is absorbed or if is going out of the atmospheric medium before interacting, it is considered as lost.

  3. Previous step is repeated on the scattered photon, until its loss or until -th scattering.

  4. At the th scattering, a random sampling is performed considering the probability :


    in order to determine whether the photon is scattered toward the telescope or not. is the phase function of the -th scattering, and the solid angle associated to telescope pupil entrance. If the photon is not scattered toward the telescope, it is rejected.

  5. Then, transmission from the -th scattering position to the detector is computed. If the photon undergoes another interaction it is lost; if not, it reaches the telescope lens.

“full” MC: photons“reduced” MC: photonsone photon(…)scattered ?norejected photonyesnoscattered times?yesin the detectorsolid angle ?norejected photonprobability (“full” MC)probability (“reduced” MC)yestransfered up to detector ?norejected photonyesPhoton on theentrance pupil
Figure 11: Scheme of the Monte-Carlo (MC) algorithm principle to simulate photons scattered exactly times.

As previously stated, with the ideal algorithm described above, it will be unaffordable CPU time consuming to propagate each of the UV photons produced by an EAS initiated by a primary of 10eV (see section 3.2). About 10 photons reach the telescope. The high decrease between these two numbers is due to the solid angle. If this reducing factor could be applied before photon tracking, a ray-tracing simulation should be possible. The key point comes from the particular situation of a space-based detector : it stands outside the scattering medium.

Above an altitude of 30 km, atmospheric density is very weak, and photon scattering can be neglected. As a consequence, one can define a minimum distance below which photons cannot be scattered. This distance and the associated solid angle are shown on Fig. 12. In case of a detector shipped on ISS, its altitude is , thus sr.

Figure 12: Minimum distance between the top of atmosphere (considered here at an altitude of 30 km) and the telescope and the associated largest solid angle .

The probability defined by the equation (10) can be written as the product of two independent factors defined as:


is the maximum over all phase functions involved in simulation. By normalising the phase function of the th scattering (which type is a priori not known) by , the condition is insured. is small enough to insure . Random draw of probability is split into two independent draws of and . Since has been set up to be independent of all other stochastic processes involved in photon propagation, random computation of can be performed at the very beginning of simulation chain. Moreover, considering that it can be identically applied to all photons, only one evaluation is performed, according to a binomial law with (, ) parameters. is the total number of photons produced by the shower, . Then the reduced number of photons to propagate when simulating the th scattering, is given by:


where the Poisson distribution is the binomial distribution limit, in case of a large sample size and low probability. The described algorithm is then applied, taking probability instead of in the 4th step. All interaction types described in section 3.3.1 - molecular absorption, scattering on air molecules, aerosols, clouds, and ground - are treated by this algorithm.

The treatment presented above is valid for a fixed scattering order. To obtain the contributions of different scattering orders and , the described algorithm has to be performed for one set of generated photons , until the th scattering, then for another set , until the th scattering, in a full independent way. Before computation, a maximum scattering order has to be fixed. The contribution of the scattering processes to the signal is the sum of the photons scattered exactly 1 time, 2 times, 3 times, …, times. “Zero” order corresponds to fluorescence photons produced in the detector solid angle, and transmitted without interacting up to the telescope. Considering the first scattering order in clear sky conditions, the different contributions come from:

zero order
  • fluorescence photons straightly transmitted

first order
  • Cherenkov photons reflected on ground

  • Cherenkov photons undergoing molecular scattering

  • fluorescence photons reflected on ground

  • fluorescence photons undergoing molecular scattering.

The main advantage of this algorithm is to be able to simulate radiative transfer of photons emitted by any type of light source, taking into account the multi-scattered photons at any order; but its main limitation is the computing time, which can be high, despite the applied reducing factor. The computing time is strongly dependent on maximum scattering order , as well as on number of photons to propagate . The latter is related to the maximum solid angle, fixed by detector setup, to the number of emitted photons, related to shower parameters () and to the maximum value of the phase function. can vary from 0.3 in case of clear sky condition to 1000 in case of cloudy condition. Simulation runs show that computing time is not affordable in the latter case.

To bypass this limitation, the cloud phase function involved in the last scattering process is approximated. Let us highlight that the cloud phase function involved in previous scattering processes is not approximated ; it concerns only the last scattering process, if it occurs within a cloud. Two kinds of approximation have been studied. The first one uses an analytical form, the Double Henyey-Greenstein function [55]. The second one uses an averaged phase function peak over a limited range of 10.

In the current state of ESAF development, aerosol phase functions are treated in a simpler way, using the analytical approximation for every scattering process occurring on aerosol layer.

Figure 13: LIDAR shot at 1064 nm from ground to 1 km distant cumulus cloud, 300 m deep and with an extinction coefficient of 17.25 km (no atmosphere contribution, cloud only). Contribution of the multiple scattering component for a receiver field of view of 10 mrad (top) and 1 mrad (bottom). ESAF results are obtained in simulating scattering process until order 8. In both cases the ESAF results are superimposed to MUSCLE ones, showing that they are in very good agreement.

3.3.3 Validation tests

To proof that both propagation algorithms are efficient to properly simulate radiative transfer, validation tests were carried out. First, the propagation in clear sky conditions of photons emitted by a given shower (fixed energy and zenith angles) has been simulated by the two designed algorithms. Results concerning Cherenkov photons propagation (molecular single scattering, ground reflection) were compared. The time distribution of photons at telescope entrance pupil and their last scattering position in FoV are similar, clearly proving the consistency of the two algorithms.

In a second step, studies have been undertaken to test the validity of the most efficient algorithm, namely the “reduced” Monte-Carlo algorithm. The aim is to compare it to simulation tools developed within the atmosphere physicists community. One study, handled by MacKee et al. [56], concerns cloud albedo change with optical depth, predicted by a Monte-Carlo simulation and checked experimentally. Results got with ESAF are in agreement with those obtained by MacKee et al. The second set of tests were done to reproduce studies of the MUSCLE (MUltiple SCattering Lidar Experiments) working group [57]. The aim is to predict LIDAR shot response after multiple scattering in clouds. The results obtained with ESAF are compared to the MUSCLE ones: a good agreement has been found (see Fig. 13 and reference  [23] for more details) that actually proves the validity of the “reduced” Monte-Carlo propagation algorithm implemented in ESAF.

Let us specify that the two kinds of cloud phase function approximations have been compared during these tests. It has been shown the use of an averaged peak is the best solution, and will be used further in this paper.

With this propagation package in ESAF, generated fluorescence and Cherenkov lights are transfered up to the space-based detector. Thanks to the designed algorithms, careful studies of the signal at detector pupil entrance are possible dealing with different atmospheric situations: clear sky, aerosol layer, clouds.

3.4 Detector response simulation

The ESAF detector simulation is the part of the framework which is more tightly connected to the EUSO design. Most of the implementation is based on a monocular cylindrical telescope with a photomultiplier (PMT) based photodetector. The main scheme of the detector is shown in Fig. 14. Nevertheless, the ESAF detector simulation structure retains large degree of flexibility and can be adapted to a different detector concept.

The telescope simulation part is divided in 2 main modules:

  1. The Monte-Carlo photon tracking engine through the optical system.

  2. The electronics and trigger response.

3.4.1 Optics and focal surface

The detector simulation starts from the photons in front of the entrance pupil. One by one the photons are propagated through the optical elements of the detector; in Fig. 14 the baffle, the optical system and the focal surface are visible. Three simulated optical systems are provided within ESAF:

  1. The EUSO Phase A baseline optics, a double Fresnel system. The tracking through this optics is simulated by tracking each individual photon entering the pupil up to the focal surface. A detailed description of the fine structure of the Fresnel lenses allows to reproduce and study the optics distortions. As shown in Fig. 15, a large fraction of the light in the Point Spread Function (PSF) of the double-Fresnel optics is spread all over the focal surface: the so-called stray-light. A commercial program (CodeV, [42]) was used to simulate outside ESAF the exact ray tracing properties of the lenses and the result of this CodeV simulation was used to create an accurate parametrization of the PSF, of the amount of stray light and of the aberrations. Within ESAF, this parametrization is used to trace the photons. We underline that the light that contributes to the signal and, particularly, to the trigger, is the one that is contained in a circle of the size of the focal surface pixel. Fig. 16 shows that in a typical pixel of about 5 mm diameter (EUSO design), about 60% of the photons contribute to the trigger, while the other 40% is lost as stray light. Fig. 16 shows also the (weak) dependence of this number as a function of the angle that describes the position of the shower in the field of view (see Fig. 23).

  2. A second multi-purpose optical system, whose characteristics are implemented through a parametrization of the optical response (PSF, stray light, aberrations) and are defined by a set of configuration files. Highly flexible, it allows an effective fast simulation without tracking individual photons. It is particularly useful for studies where the performances of the optics have to be tuned according to the overall requirements, or as a reference to compare the performances of the full simulation, or to speed up the simulation for large statistics studies. The parametrization and the configuration files were obtained by running the full simulation. Most of the work in this paper was obtained with this faster implementation but we have checked that the results obtained with the full simulation are always consistent.

  3. A Geant4 based implementation of the optics and of the ray tracing (see next paragraph). This Geant4 implementation is new and was not used in this work.

Ideally the efficiency of the optics scales with the angle in the FoV, , as ; for the EUSO optics the mean efficiency is of about 45%. To this number, another 60% light loss must be added because of the optical aberrations and stray light, as determined from the simulations.

Figure 14: Main scheme of the detector.
Figure 15: An example of the irregular PSF of the Fresnel optics generated using an incident flux at an angle in the FoV. (a) Distribution of photons on the whole focal surface. (b) Zoom on the spot of the previous picture. The red circle corresponds to a bucket of of radius.
Figure 16: Fraction of energy (fraction of photons) contained in a circle as a function of the circle diameter.
(a) (b)
Figure 17: Focal surface layouts in ESAF. (a) Cartesian layout. (b) Polar layout.
Figure 18: Image of a EAS track () on the focal surface. No background is simulated. (a) 2D view. The small squares are the pixels, the big ones the PMT. (b) 3D view including also the time counted in Gate Time Units (GTU) of . The and axis are the coordinates of the pixels while the axis is the time from the first signal photon counted (in GTU). On the -plane is plotted the profile of the pixels with signal.

The focal surface module implements the geometry of the EUSO photodetector. The two available layouts of the focal surface are illustrated in Fig. 17. The focal surface is based on the modular structure of small autonomous units called elementary cells [58], each one being a matrix of multi-anode photomultipliers. An optional optical adapter can be placed in front of the PMTs (efficiency of about 95%). The focal surface elementary cell layout geometry is loaded at runtime from the layout configuration file. The simulations show that the focal surface, in the polar layout, has a “filing factor” efficiency of . The tracking engine propagates the photons up to the PMT photocathodes and converts them into an electronic signal, stored for the electronics simulation.

Any PMT (or optical adapter) that fulfills the elementary cell specifications and constraints can be easily added to the focal surface simulation. Currently two models444In the PMT models ESAF takes into account the quantum efficiency as a function of the wavelength, the collection efficiency, the gain, the cross-talk between pixels, the PMT geometry, the time width of the signal and the dark noise rate. The typical PMT efficiency is 25% (quantum) plus 70% (collection). are available, the 36-pixel Hamamatsu R8900 and the 64-pixel Hamamatsu R7600.

In the Fig. 18a and 18b there is an example of how an EAS appears on the telescope focal surface.

3.4.2 GEANT4 simulation of optics in ESAF

A detector simulation chain with the possibility to trace particles through the detector using the GEANT4 package [59] has been implemented in ESAF. We have kept the original way of ray-tracing for checks and compatibility. This upgrade makes ESAF much flexible for a simulation of any detector and not necessarily an EUSO-like detector. As two examples we implemented an EUSO like detector with two Fresnel lenses and a diffractive surface between them and a TUS-like detector [60] with a Fresnel mirror. Both examples are shown in Fig. 19.

Figure 19: EUSO Fresnel lens (top) an TUS Fresnel mirror (bottom) detectors built with Geant4. Fluxes of photons focused on their respective focal planes are also displayed.

Almost all segments of EUSO Fresnel lenses have parabolic surfaces. Every segment is done with boolean operations on the basic shapes (G4Tubs and G4Paraboloid). More complex surface shapes were approximated by cones. Approximation precision can be chosen on demand. For the current simulation it is equal 0.1 m. A piece of a EUSO Fresnel lens implemented with GEANT4 is shown in Fig. 20.

Figure 20: A piece of a EUSO Fresnel lens simulated with help of GEANT4

The following standard optical processes are defined in the GEANT4 implementation: Fresnel refraction, Fresnel reflection, total internal reflection and absorption. In addition to the standard optical processes, ESAF uses a ”special“ process parametrisation: the diffraction on the second Fresnel lens. In fact, the second lens of the EUSO detector has two surfaces, one with Fresnel structure and the second one with a normal diffractive surface. The role of the diffractive surface is to reduce the chromatic aberrations. For the simulation of the diffraction process on the second lens we use a parametrization. The momentum direction correction depends on the material refraction index and on transverse component of position on the lens surface. To implement it in the Geant4 simulation we have created a special class that inherits from G4VDiscreteProcess. We have added a thin vacuum layer just after the flat lens surface, and have assigned proper material parameters. Thus, rays going out of the lens are affected by ordinary boundary processes, and, if not reflected back, they go through that thin vacuum layer and are therefore handled by the new class.

The GEANT4 implementation of the simulation was not used in the work described in this paper. For this reason, we do not give a detailed account of GEANT4 - ESAF capabilities. However this capability will be very important for future detector simulations and it will be documented in a future work.

3.4.3 Electronics and trigger

The electronics simulation converts all the PMT signals into front-end electronics (FEE) counts, simulates the front-end chip response and applies the trigger patterns.

In the front-end the PMT signals are sorted by time, divided in time units called Gate Time Units (GTUs) and counted. In the baseline configuration of EUSO the duration of a GTU is . An amplitude threshold and the double hit resolution of the front-end chip are taken into account as well, according to the specifications of the EUSO Phase A design report [16]. The electronics has a typical efficiency of about 90%.

The atmospheric background contribution is added directly after the front-end simulation as additional FEE counts. The motivations for this solution are described in detail in section 3.5.

The trigger in ESAF is implemented as a set of ’Trigger Engine’ objects, one for every trigger algorithm. Several engines can be applied to the event, giving the possibility to compare their efficiency on an event-by-event basis. The specific trigger algorithm used for the studies presented in this paper is described in the section 3.4.4.

3.4.4 The level-1 trigger algorithm: Tracking

The image of a developing EAS on the focal surface is a bright spot moving on a straight line whose brightness rises up to the maximum and gradually fades away. At a given time, the spot is mainly contained in a single pixel, as required by the optics design. Integrated over the development, the shower appears as a track in space and time whose length and brightness depend on the primary energy and incident angle but also on the EAS position in the FoV.

In a space-based experiment, only those EAS that produce a track on the focal surface, standing out from the surrounding nightglow background, can be detected and provide any useful information. Without the track structure visible, the EAS signal is not distinguishable from the background fluctuations.

Therefore two conditions can be established for an EAS to be detectable:


The signal in a given pixel and GTU must be well above the average background counts.


The bright, or active, pixels must be contiguous in space and time.

At pixel level, to set a “brightness” requirement, a background dependent threshold can be established:


where is the mean number of background counts per pixel per GTU, and a coefficient depending on the requested background rejection efficiency. Since the random background is a poissonian distribution with mean , can be regarded as a number of standard deviations. This parametrized threshold is intrinsically adaptable: for fixed, can be continuously measured and the threshold can be automatically updated.

The contiguity criteria requests that the signal must form a track in space and time. Starting from the first active pixel, the next active one in the next GTU must be the same pixel or one of the nearest neighbors. The minimum number of consecutive active pixels, the track length , can be assumed as the second requirement a signal must satisfy to be recognized as a possible EAS candidate.

We have studied the trigger efficiency using the contiguity tracking trigger (CTT) developed during the EUSO Phase A. It works at the Elementary cell (144 pixels) front-end level, that is an ensemble of 22 closely packed multianode PMT (36 channels each). The CTT is a simple algorithm that checks the signal contiguity in space and time, as explained before. Thought to be integrated in the front-end chip as first level trigger, it can be implemented using purely combinatorial circuits with low cost in term of power and mass budget. Due to the large background non-uniformities, cannot have the same value for every channel on the focal surface. Instead, we used as one of the variables of the trigger. In the following we shall denote a particular CTT configuration with the notation .

The rate of fake triggers (over the whole detector) due to the random background is easy to compute [21, 22]:


where is the number of elementary cells (1386 in the EUSO Phase A configuration), is the duration of a gate time unit and

is the poissonian probability to have a random track of length with a pixel threshold .

The maximum acceptable overall fake trigger rate is constrained mainly by the available telemetry resources. Thus, the maximum fake CTT rate is then limited by the time the second level trigger further processes the data. The final number depends of course on the details on the second level trigger, but a realistic estimation, based on the current technology, sets this limit to .

3.5 Simulation of the random background

The main source of background to the EAS signal is the atmospheric nightglow (see [61] and references therein). Other sources like moon light, lightenings and man-made lights are not constant in time. Currently, ESAF is able to simulate only the effect of an uniform background. This is not a major issue since the main interest is to study the detector capabilities with various levels of average background. The mean value of the background radiance used in this work is , as measured by different experiments (e.g. BABY [62], TATIANA [63]).

Assumed uniform and isotropic on the EAS development timescale, it floods the entrance pupil with a constant flux of photons. The resulting flux on the focal surface is not uniform, because of the optics effects. The average background counts per pixel per GTU, , varies from the focal surface center to the border with a profile that depends on the optics used. sums up the contributions of photons from all incoming directions and therefore the stray-light (see Fig. 15) cannot be neglected.


where is the pixel position, the number of background photons per surface unit, the photo-detection efficiency and the pixel side. The distribution is obtained by simulating several million photons with an uniform incident direction distribution, according to the background radiance. For such a purpose ESAF has a special running mode where only the optics is loaded, optimizing the memory and CPU consumption.

Once this has been done, and accordingly with the user choice, ESAF can simulate the random background over the whole focal surface only in those elementary cells that have at least one pixel with signal. In each pixel the simulation adds a random number of background counts for every GTU, following the Poisson distribution with mean , where is the GTU time length.

4 Shower reconstruction

After the observation of one event, the following questions should be addressed: what are the direction, the energy and the nature (nucleus, gamma, neutrino, ) of the cosmic ray at the origin of the detected EAS ?

ESAF has a reconstruction package which consists of special modules aiming to recover these informations from an event. Let us consider below the underlying principles of the reconstruction software which consists essentially of the following modules:

Pattern recognition:

the goal of this module is to find the “shower track”, i.e. a set of spatially correlated pixels on the focal surface that describe a moving point as a function of time (the projection of the EAS development on the focal surface). The main problem here is to disentangle the pixels that contain some signal to those that are pure background.

3D Direction reconstruction:

the goal of this module is to find the direction of the primary particle in space. The 3-dimensional reconstruction is possible by using both the direction of the projected track on the focal surface and the arrival time of the photons.

and Energy reconstruction:

the goal of this module is to find the altitude in the atmosphere where the EAS reaches its maximum of development and the corresponding atmospheric slant depth, and to reconstruct the energy of the primary particle. The fluorescence light emitted by an EAS is proportional to the deposited energy. However, the collected light is not because there are important corrections that depend on direction, position on the focal surface, and detector efficiency. The and its variation with the energy are very important for primary particle identification. Since the energy and are correlated, it is not possible to infer a precise determination of each variable in an independent way. Therefore, this module reconstructs the best estimate of the energy and of the .

4.1 Pattern recognition

Once a shower track has been triggered by the apparatus, the pattern recognition module selects, within every considered time interval (GTU), the pixels containing signal from the pixels containing background only. Since the background light is expected to be randomly distributed both in direction and time, the space-time correlation in the focal plane is the key tool to select clusters of pixels which belong to the shower track on the focal surface.

The pixels selection is implemented in ESAF by a simple multi-step analysis which, within the focal plane, searches for clusters of pixels in space and time separately. A following module also analyzes the X versus time and Y versus time correlations inside the focal plane by means of a linear fit, followed by the exclusion of the points far away from the fitted line. A more sophisticated method using similar concepts, but based on the Hough transform technique [64] is also successfully implemented in ESAF.

4.2 Direction reconstruction

Let us introduce a standard reference system (right-handed) with the origin in the center of the detector focal surface and the -axis pointing towards the nadir of the instrument (see Fig. 23). 555We are assuming that there is no tilt between the nadir and the orientation of the instrument optical axis.

The EAS axis is described in this reference system by an unit vector in terms of two polar angles and :


In an analogous way we define the unit vector of the direction that is seen by the -th front-end pixel in the FoV:


Therefore, for each selected pixel on the focal plane there is an associated direction in space where the fluorescence light activating such pixel is expected to come from.666For any detector configuration, a map that associates each pixel with a unique direction in the FoV is built simulating an uniform photon flux in the FoV through the detector.

Such directions in space are expected to belong, within instrument resolution and neglecting the shower lateral extension, to the plane containing the shower axis and the detector itself. Using the information of the light arrival time in each pixel, it is then possible to look for the unique track in space associated to the measured direction-time pattern. Such pattern depends strongly on the direction of the shower axis in the detector-shower plane, but very weakly on the distance of the shower maximum to the detector because of non-proximity effect. This is one of the advantages of space-based observation: the altitude of the shower maximum may vary up to 20 km, depending on the shower inclination, but it is always a small distance compared to the satellite typical altitude of 400 km, which means that the collected light changes at most by 10%.

With a monocular detector, it is not possible to determine by means of a fitting procedure the detector-shower distance, unless the resolution is so good to allow a measurement of the spatial width of the focal surface track. The shower detector distance can be simply fixed to the expected average value for each zenith angle of the incoming shower. However, such an information can be inferred by the arrival time of the Cherenkov light reflected on the ground, if it is available, since it is strongly correlated to the distance of the shower from Earth and detector. With the detector-shower distance inferred by such procedure, the fit of the direction-time pattern yields a better determination of the shower-axis direction inside the shower-detector plane. In case the Cherenkov light is not detected, the detector-shower distance can be inferred by the spatial width of the fluorescence light profile (see section 4.3). In both cases, the determination of the detector-shower distances and of the zenith angle is made more precise using an iterative reconstruction process.

The first step that leads to the determination of the EAS direction is the reconstruction of the Shower-Detector Plane (SDP), i.e. the ideal plane that contains both the EAS track and the detector. The SDP can be reconstructed by ESAF fitting the “moving” EAS image on the planes -time and -time. In this way one obtains a projection of the shower track on the plane tangent to the Earth surface at the telescope nadir point. The SDP is then the plane that contains this projection and the detector itself.

The second step is the fit of the track on the SDP in order to infer the direction of the shower axis. In ESAF two different classes of track fitting are implemented [65].

Numerical Exact (NE)

For each signal FEE count we compute its expected arrival time, , as a function of the EAS axis direction (unknown), the pixel direction in the FoV (see above) and other free parameters777These free parameters depend on the particular track fitting algorithm. They can be, for example, the position in FoV, the time and the distance from the detector of the EAS maximum., . These times can be exactly calculated using only geometrical and physical relations. Since we know the experimental arrival times, , we can minimize the function

where the are the uncertainties on the measured arrival times, via numerical minimization methods. From the minimization results we know the EAS direction.

Analytic Approximate (AA)

In this case we use approximate formulas, e.g. assuming that the EAS track has a constant angular speed on the focal surface, that can be analytically solved using linear fit procedures.

4.3 and Energy reconstruction

The principle of the method described in this section is to assume a functional form for the longitudinal development of the shower, to calculate the expected number of photoelectrons at the telescope as a function of time and to vary the parameters of the shower function until an agreement with the measured signal is obtained. From these parameters, one can deduce the energy and of the EAS.

To calculate the number of photons emitted along the track and the fraction of them that reach the telescope, one has to know the shower track position in space and the local atmospheric parameters. Moreover to obtain the expected signal, one has to estimate the detector efficiency.

In what follows and as functions of time are needed to fix the direction in space where the fluorescence light is expected to come from as a function of time. These functions are obtained by fitting and of each active pixel with presumably a signal by a linear function of time. An example of such fit is shown in Fig. 21.

Figure 21: An example: vs time (upper plot) and vs time (bottom plot) (each GTU is 2.5 s).

4.3.1 Background and signal estimations

A cluster of pixels belonging to the shower track is taken from the pattern recognition module to assign for each time stamp a list of pixels presumably enriched by the shower signal. For the background study, a similar second assigning map is created which contains the pixels around the signal cluster within a fixed margin, corresponding, in an EUSO like detector, to the size of two adjacent pixels. For each pixel of the track, the background can be estimated in two complementary ways: either by looking at the light collected by the adjacent pixels at the same time, or by looking at the light collected by the very same track pixels before and after the shower signal. One obtains the observed time profile of the signal by adding at each time stamp the signal from the shower pixels after subtraction of the background.

4.3.2 Optics efficiency

The focal surface may have a significant amount of holes which distort the detected profile as it is e.g. the case for the polar design (see Fig. 17b) of an EUSO like detector. Fig. 22 demonstrates an example of how the shower signal looks like on this type of focal surface; the corresponding distribution of numbers of photoelectrons as a function of time is shown in Fig. 24. One can notice at least three relatively large and rapid losses of detected photoelectrons (around times 13, 24 and 40 GTU) which correspond to holes between macrocells. The dead spaces between PMTs in one macrocell are smaller in size but also produce distortions in the profile.

Figure 22: An example of distribution of photoelectrons on the focal surface for a shower energy of 5. 10 eV and zenith angle =60. Upper plot: large scale; bottom plot: zoomed picture, in a region where the dead space among PMTs modules is not negligible.

These distortions are taken into account in the shower reconstruction process by building for each arrival time an effective optics efficiency . The previously reconstructed and of the track are used to compute:


where is the function returning 1 if the considered point with coordinates is inside a pixel, and 0 otherwise. is the optics point spread efficacy function on the focal surface at position for photons incident direction determined by . The PSF and the optics transmission are taken into account to build the optics efficacy function. The PMT efficiencies are carried separately.

4.3.3 reconstruction

As already mentioned in section 4.2, the determination of the detector-shower distance can be estimated in two ways: either by using the Cherenkov echo and the information encoded in the fluorescence signal profile, or, when the former is not available, by using fluorescence signal profile only. The latter case is, unfortunately, the most common (see section 5). In both methods, all contributions to the signal other than direct fluorescence and reflected Cherenkov lights are neglected. This allows to obtain a first estimation of the altitude of the fluorescence signal maximum . In what follows the measured signal profile is corrected at each arrival time by for optical loss and focal surface inefficiencies.

reconstruction with the Cherenkov peak

This method needs to infer from the time profile of the signal the difference () between the arrival time of the Cherenkov peak and that of the maximum of the fluorescence light . is used to estimate .

From Fig. 23, one can find that:




with the known altitude of the reflected Cherenkov surface, the telescope altitude, and the previously reconstructed parameters: the shower direction (with zenith angle ) and the unit vector pointing to the telescope and defined by and .

Figure 23: Geometry of the Extensive Air Shower seen by the detector. The Cherenkov point, when available, gives the position on the impact point of the shower on the reflecting surface (ground or cloud).
reconstruction without the Cherenkov peak

The altitude of the EAS maximum can be reconstructed relying on the fluorescence light information only. The idea of the method [66][67] is based on the simple observation that the time width of the EAS development strongly depends on air density (and thus on the altitude) in the atmosphere at which EAS develops.

A longitudinal profile of the number of shower electrons at given atmosphere depth can be conveniently described by an analytic function which depends on a dimensionless parameter as

where is the depth of the atmosphere at which EAS starts to develop and is the air radiation length. An example of is given in section 3.1 for the GIL parametrization, while the consideration below is more general. The EAS profile can be expressed as a function of detection time . Around the EAS maximum one can write

In this way we have obtained a gaussian distribution of times around with standard deviation . In order to calculate , we first obtain taking into account the following:

Thus, is expressed by:


This general equation can be simplified further for the case of GIL parametrization for which and the time width becomes


with only logarithmically depending on UHECP energy which can be easily taken into account in an iterative procedure or fixed to a mean expected value.

Because the fluorescence yield is almost constant between 0 and 20 km (it varies between 4.1 and 4.7 photons/m, see 3.2.1) and the atmospheric attenuation is a slowly varying function along the shower track, the shape of the shower electrons time profile is at first order similar to that of the total collected signal. This way one could measure and thus of the shower maximum using observable informations and reconstructed .

In principle, when the altitude of the Cherenkov scattering surface is not known, this second method may allow to recover it just as:

4.3.4 Profile reconstruction

The aim of this part of the reconstruction is to estimate the atmospheric depth profile of the shower electrons from the light flux observed by the telescope.

A functional form for this profile is assumed, a possible choice is the GIL function (1), which has three free parameters: E the energy of the shower, the atmospheric depth where the number of shower electrons reach its maximum and the depth of the first interaction. gives the electrons of the shower as a function of the slant depth along the shower track. In a first estimation of , and can be set to their mean expected value.

For each arrival time the expected signal at telescope is computed as a function of the number of shower electrons (unknown) to minimize the following function via numerical minimization methods.

where is the measured signal (background subtracted) at time and is the uncertainty on it.

Direct fluorescence signal observed at time is produced at a slant depth in the atmosphere. Once a reference in altitude and are reconstructed, the EAS track in the atmosphere can be found unambiguously and computed as an integral along the EAS track from , the shower position at time , to top of atmosphere:

this requires to know the precise density profile of the atmosphere for each event.

For each time stamp the number of fluorescence photons produced in a slant depth interval by the electrons of the shower is computed using the following formula:


where is the atmospheric density at position , is the local fluorescence yield (number of photons per unit of deposited energy) for atmospheric conditions at point , is the average energy deposit per unit depth per electron at shower age (a parametrization of can be found in [36] or set to its mean value of 2.2 MeV/() ), and is the length of the time unit.

The fraction of photons produced at that reach the telescope due to atmospheric attenuation and the distance between the shower point and the telescope are also computed. Thus the expected direct fluorescence signal is expressed by:

where is the previously computed detector efficiency at time defined in (18) (see section 4.3.2) and A the telescope aperture.

Although this has not yet been implemented inside ESAF a more precise reconstruction process would require to compute the expected reflected Cherenkov signal, the Cherenkov scattered signal and the fluorescence scattered signal, but at first order the expected total signal can only be expressed by .

Since the reconstructed shower track depends on the reference in altitude of the shower maximum, the fit is done iteratively until a convergence is obtained. Typically 2-3 iterations are enough to get a stable fit result.

An example of the fit is shown in Fig. 24.

Figure 24: An example: Fit of temporal distribution. Points with error bars are observed number of p.e., blue curve - estimated background, red curve is the sum of fluorescence and background. The Cherenkov echo seen in the last bins of the temporal distribution was identified by the code and excluded from the fit. Energy and direction are the same of Fig. 22

5 Simulation results

All the simulations discussed in the present section are based on a standard configuration which derives from the EUSO Phase A design, i.e. a 30 FoV space-born telescope installed on board of the International Space Station at a typical orbital height ranging from 360 to 450 km.

However, depending on particular cases, the diameter of the entrance pupil is varied. In order to study the light flux at the pupil, independently of the detector, the size is fixed at 5 m.

When the performances of an EUSO-like detector are involved, the pupil diameter is fixed at the standard EUSO diameter of 2.5 m. In the trigger efficiency studies different pupil diameters have been tested.

5.1 Light flux at the pupil in clear sky conditions

Results presented in this section concern light flux arriving on a 5 m pupil diameter.

Simulation results were obtained from showers parameterized with the GIL formula and a fluorescence yield based on the Kakimoto model completed with Bunner rays. The reference detector altitude is assumed to be 430 km and the atmosphere model is a US-Standard clear sky. The Earth ground is assumed to be a lambertian surface, with an albedo of 8%. Propagation of light is performed with the “Reduced" Monte-Carlo algorithm.

Figure 25: Number of photons reaching the pupil as a function of time. The signal is integrated within the whole FoV. Top histograms refer to the last photon interaction and bottom histograms refer to the number of scattering interactions before reaching the pupil. The simulation has been carried out with 4 orders of scattering in a US-Standard clear sky atmosphere. The energy of the shower is 10 eV and the zenith angle is 60.
Figure 26: Distribution of collected photons positions in the FoV. The simulation has been performed with 4 orders of scattering in a US-Standard clear sky atmosphere. The energy of the shower is 10 eV and the zenith angle is 60.

Fig. 25 and 26 show an example of the simulated light at the pupil with the different light components: direct fluorescence, ground reflected lights (fluorescence and Cherenkov), Rayleigh scattered fluorescence and Cherenkov lights. The contribution of the first 4 orders of scattering are also shown. The signal as a function of time has been integrated within the whole FoV. As one can see from Fig. 25 the direct fluorescence light contributes only to 37% of the total signal. This result shows the importance to include the scattered light contamination in the simulation program. The contribution of scattered light decreases by a factor of 2 at each scattering order from 35% at order 1 to only 4% at order 4. Although scattered photons are numerous, their positions in the FoV constitute a large halo around the main track as shown in Fig. 26. Besides this, they reach the pupil later than direct signal and thus do not affect the first part of the profile. As a consequence, during the reconstruction process, appropriate pixels selection around the main track in space - and time will reduce the scattered light contribution as it will be explained later in section 5.1.3.

5.1.1 Direct fluorescence light

Fluorescence light is related to the shower development in atmosphere, therefore direct fluorescence light is the signal from which the energy and the are derived. We report here its main characteristics.

As underlined in section 3.2, inclined showers produce more fluorescence photons compared to vertical ones. Traveling in a less dense atmosphere these photons are moreover less scattered or absorbed, therefore direct fluorescence light flux at pupil depends strongly on the shower zenith angle. Fig. 27 shows this dependence for showers impacting at nadir. The amount of direct fluorescence light increases by more than a factor 6 between vertical showers and horizontal ones (zenith angle above 75 degrees).

Figure 27: Number of direct fluorescence photons at pupil as a function of shower zenith angle for showers at the telescope nadir. Numbers of photons are normalized to the number of photons for a vertical shower, which is around 8000 for a primary energy of 10 eV.

At fixed zenith angle, direct fluorescence light reaching the telescope depends only slightly on the shower position in the FoV. This weak dependence of signal intensity within the FoV is a specific feature of a space-based detector, due to the large distance between showers and telescope. Specific studies carried out with ESAF prove that the signal from showers at the edge of the FoV is reduced by 40% compared to the nadir case.

As already mentioned, at higher altitude, where the atmosphere is less dense, the shower is longer than at lower altitude. On the contrary, because the fluorescence yield is almost constant between 0 and 20 km, the number of photons at the maximum depends weakly on the altitude. Therefore at fixed zenith angle the ratio between the total and the maximum number of fluorescence photons is an indicator of the shower altitude as illustrated in Fig. 28 and could be used when the reflected Cherenkov signal cannot be observed.

Figure 28: Ratio of total over maximum number of fluorescence photons as a function of the shower maximum altitude at different zenith angles. As is expressed in GTU, the ratio of over is expressed in GTU. The spread of the curves is due to the random first interaction point of UHECP used in the simulation. No correction is applied for the shower distance.

5.1.2 Ground reflected Cherenkov

Cherenkov light that is beamed around shower axis, travels with shower front down to Earth and then back to the detector for the part that is reflected. It provides an extrapolation of shower axis to Earth, therefore the time difference between Cherenkov peak and fluorescence light maximum is an indicator of the shower altitude.

Fig. 29 shows results obtained at different zenith angles. Nevertheless, this signal can only be used if it is sufficiently intense.

Figure 29: Time difference between ground reflected Cherenkov and direct fluorescence peaks as a function of the altitude of the shower maximum at different zenith angles. The spread of the curves is due to the random first interaction point of UHECP used in the simulation. No correction is applied for the shower distance.

In case of inclined showers, the Cherenkov beam has a longer path to go down to the ground and is therefore more attenuated than in the case of vertical showers. Moreover Cherenkov photons are more dispersed around the shower impact due to the photon angular distribution. As a consequence the Cherenkov peak amplitude decreases when the zenith angle increases (Fig. 30). Being quite constant up to 40 the signal is reduced by a factor of 30 at 75.

Figure 30: Amplitude of the ground reflected Cherenkov peak as a function of shower zenith angle. Numbers of photons are normalized to the number of photons for a vertical shower, which is around 2900 for a primary energy of 10 eV.

In most cases Earth surface cannot be modeled by a Lambertian surface and specular reflection has to be used.

Taking into account showers with an impact at nadir, in the case of a specular reflection on oceans with a wind speed of 2 m/s (respectively 9 m/s), one can deduce from Fig. 9 that events with a zenith angle greater than 20 (respectively 30) will give a weaker Cherenkov peak than in Lambertian case. On the contrary, for smaller zenith angles, the Cherenkov peak should be greatly enhanced.

The anisotropy of the specular reflection has consequences for the dependence of the reflected Cherenkov light on the shower position in the FoV.

Figure 31: Amplitude of the ground reflected Cherenkov peak as a function of the shower impact position in the FoV for showers with a zenith angle of 20 in case of specular reflection over ocean with a wind speed of 2 ms. The total albedo is 8%.

As an example, Fig. 31 shows the amplitude of the ground reflected Cherenkov peak as a function of the position of the shower impact in the FoV for showers with 20 zenith angle. For a wind speed of 2 m/s, 76% of the showers that impact in the FoV have a reflected Cherenkov light lower than the one expected in the Lambertian case. Furthermore, for showers with zenith angle greater than 30, the specular peak does not point toward the telescope, regardless of the shower impact position in the FoV. Over the oceans the hypothesis of a Earth Lambertian surface seems to be too optimistic and the anisotropy of the specular reflection will be responsible of a dependence of the reflected Cherenkov light with the position in the FoV.

We conclude, therefore, that in case of specular reflection on Earth, the Cherenkov light of most of the showers will not be detected. Simulation results have been produced with a constant albedo of a few percents. In case of specular reflection, the Cherenkov peak could be enhanced only if, by chance, the shower angle and its impact position in the FoV allows the light to be reflected towards the telescope.

5.1.3 Scattered light

In this section, the term “scattered light” represents all photons that have been scattered by Earth surface (assumed lambertian) or by atmosphere in clear sky conditions. The signal component previously named “ground reflected Cherenkov” is not included in this scattered light component, but is considered as a component of the basic signal. As already said, scattered photons are slightly time delayed compared to the direct signal and spread around the main track in the FoV. Therefore they can partly mix with the basic signal as a noise.

Figure 32: Number of photons reaching the pupil as a function of time at 4 zenith angles. The signal is integrated in the whole FoV. In gray are the basic signals: direct fluorescence and ground reflected Cherenkov. The single scattered signal is superimposed in light red and the multi-scattered signal in dark blue.

Fig. 32 shows time profiles of photons reaching the pupil at different zenith angles. The plotted signal is integrated in the whole FoV. Although the scattered light increases with shower zenith angle (by a factor of 5 for single scattered photons between 0 and 80 and by a factor of 4 for multi-scattered photons), its contribution is of the same order of magnitude as the direct fluorescence light. However part of the scattered signal is not correlated with the track in the FoV because most of the photons come from a last scattering position (in the atmosphere) far from the shower axis. Thus most of the scattered photons will be rejected in the same way as background photons.

Figure 33: Integrated single scattered photon number as a function of the distance to the shower axis for a vertical shower. Fluorescence and Cherenkov photons are compared. Distributions are normalized to 1.
Figure 34: Contributions of selected scattered photons from different origins related to direct fluorescence light as a function of shower zenith angle in clear sky conditions. Selected scattered photons are those correlated in space and time with the direct fluorescence signal.

As an example Fig. 33 shows the lateral distribution of the last scattering position in the atmosphere for vertical showers. As one can see, due to the fact that Cherenkov light is beamed at emission, 90% of single scattered Cherenkov photons have their last scattering position in the atmosphere located less than 0.4 km from the shower axis, in contrast to only 10% of single scattered fluorescence photons. Simulations performed to study multiple scattering effects prove that only 5% of the photons scatter more than once before reaching the detector.

If not taken into account during the reconstruction process, the scattered light component would entail systematic errors on the shower parameters. The amount of light is increased and the detected shower profile is distorted.

To estimate precisely the contribution of scattered light to the detected signal, we have selected the scattered photons, at the photo-electron level, using a simple detector simulation. It uses EUSO configuration (pupil diameter of 2.5 m), applying a global efficiency coefficient (5%) to convert photons into photo-electrons. The FoV being pixelized and time profile being discretized in GTU, one selects at a given GTU the scattered signal standing in pixels containing the direct fluorescence signal.

Figure 35: Contributions of scattered photons to the integral of the signal , the integral of the first half part of the signal and the maximum of the signal .

It appears that the contribution of scattered photons represents less than 30% of the direct fluorescence light. This is shown in Fig. 34 where contributions of scattered photons from different origins are represented as a function of the shower zenith angle. The main contribution is that of single scattered Cherenkov photons that represents slightly more than 20% of the direct fluorescence photons, in contrast to only less than 4% for single scattered fluorescence photons and 1% for multi-scattered photons.

It has to be noticed that presence of scattered photons has consequences not only on the total number of photons reaching the pupil but also on the shape of the reconstructed profile. This can be illustrated by comparison of scattered photon contributions to the integral of the signal, to the integral of the first half part of the signal and to the maximum of the signal as shown in Fig. 35.

5.2 Triggering efficiency and energy threshold

All the simulations discussed in the present section are based on a standard configuration which derives from the EUSO Phase A design at a typical orbital height .

From this basic configuration, four different detectors were implemented, identical in all parts but for the diameter of the entrance pupil: (EUSO Phase A configuration), , and . In this way ESAF is able to evaluate the performances of EUSO and of possible larger detectors.

Figure 36: Trigger efficiency curves for various trigger configurations ( detector): (a) FTR above ; (b) FTR below (c) Efficiency with fixed and .

A sample of simulated EAS with energy in the range (flat energy distribution in a logarithmic scale and angular distribution uniform in solid angle) was used for all the configurations.

Figure 37: Efficiency of the [3;4] and [4;4] trigger patterns on three ranges in . The relative weight of the three samples are (a) , (b) , (c) . The red line at is drawn as a reference.

The showers were generated by GIL (see section 3.1) in clear sky and US-Standard Atmosphere (no clouds). A parameterized optics with a 2D gaussian PSF of RMS, a GTU of and an overall photo-detector efficiency of 10% were used according to the EUSO Phase A design.

Figure 38: Trigger efficiency curves for various sizes of the detector using the pattern [3;4].

Because of the limits imposed by the available telemetry and computing resources in the EUSO project, only the contiguity tracking triggers that satisfy were considered useful. Future detectors may have less stringent requirements.

As shown in Fig. 36, if the condition is not satisfied, the efficiency curve flattens at low energies to a value that depends on (see section 3.4.4 for the definition of ), else flattens to zero at low energies. Furthermore, when the value of is raised to a fixed , the efficiency curve is shifted toward higher energies (without changes of shape), i.e. the full efficiency is obtained at a larger value of .

The dependency of the trigger efficiency on the zenith angle is shown in Fig. 37. The efficiency curve is shifted toward higher energies for and toward lower energies for .

To compare the detectors on the same grounds, the continuity trigger pattern must satisfy for all detectors. was the selected trigger for the further studies being the less restrictive pattern satisfying this additional condition.

The results are shown in Fig. 38. As expected the curves are shifted toward lower energies when increasing the entrance pupil diameter. The energy threshold of the experiment scales as the inverse of the diameter or, equivalently, as the inverse of the square root of the photo-detection efficiency.

Simulations with ESAF show that an EUSO like detector has full detection efficiency around 2. 10 eV and 50% efficiency at 8. 10 eV. Better efficiency can be obtained only by enlarging the detector size or by using better photo-detection devices than commercially available MAPMTs. We recall, however, that these results are obtained in ideal conditions, not taking into account variable weather and dead pixels in the focal surface.

This result is really very weakly dependent on the specific triggering algorithm. The lack of efficiency is related to lack of signal at low energy.

It is worth recalling that these conclusions apply to photo-detector similar to MAPMTs, that have a global photo-detection efficiency (quantum efficiency PMT practical factor geometric factors) of about 10%. To our knowledge, better devices do not exist yet, but new photocathodes with much improved quantum efficiency are already available and might yield more efficient MAPMTs in the near future. This of course may have a significant impact on efficiency. However, higher photo-detection efficiency also means higher nightglow background, so the gain would not be linear. To be more precise, the gain in photon collection efficiency scales as the pupil area in Fig. 38. For instance, referring to a 2.5 m diameter pupil (right-most curve in Fig. 38), an increase of efficiency of a factor 2.56 (significant improvement!) corresponds to 4 m diameter (second curve from right).

5.3 Angular resolution

Fig. 39 shows the angular resolution concerning the identification of the shower axis direction in space in presence of typical random background of 500 ph mns sr (see section 3.5)888We recall here that the typical night glow background in EUSO is 500 ph mns sr, reaching 1000 with 1/4 of moon. With more than 1/4 moon, the observation of EAS is not possible. ; in the figure, is the angle in space between the simulated and reconstructed directions which defines the cone containing 68% of the events, i.e. one standard deviation in angle; is the zenith angle of the incoming shower calculated at the ground impact point. Standard EUSO parameters are used for this simulation. The showers are simulated with UNISIM in clear sky conditions neglecting the multiple scattering effects.

Figure 39: Angular resolution. is the angle in space between the simulated and reconstructed directions which defines the cone containing 68 % of the events.

The curves labeled refers to exact formulas and numerical methods, while refers to approximated formulas and analytical methods (see section 4.2). As expected the direction reconstruction is far better for very inclined showers, since for such inclinations the track length is greater and consequently the total number of detected photons and the number of hit pixels are higher (see section 5.1.1).

At 10 eV the angular resolution is close to 10 for zenith angles smaller than 40, it goes down to a few degrees in the angular range 50-70 and finally approaches about 1 for very inclined showers with zenith angle greater than 80. The angular resolution reaches 1 already at 70 zenith angle for 2 10 eV and at 60 zenith angle for 5 10 eV.

In a conservative way, the detection of Cherenkov light is not taken into account here and the detector shower distance is set in the reconstruction at the expected average value. Shower to shower fluctuations yields approximately 0.5 km rms error on the altitude.

The relation between the uncertainty in the estimation of the shower maximum height (and so of the detector-shower distance) and the uncertainty in the reconstructed direction has been studied by Monte-Carlo simulation and analytical calculation. In average we find an uncertainty of 0.2 in the reconstructed zenith angle for every km of error in the altitude of the EAS maximum.

5.4 Energy and reconstruction performances

As previously seen, the event reconstruction starts with the determination of the shower track. Then a functional form for the shower longitudinal development is assumed and the expected collected light is calculated. The energy and of the EAS are deduced by varying the parameters of the longitudinal shower function until an agreement with the measured signal is obtained.

Several sources of uncertainties can be identified at each step of this reconstruction process. ESAF is a useful tool to evaluate most of them. Different samples of proton showers have been generated. In each sample, either one shower parameter (the energy, the zenith angle, the shower maximum altitude, the ground impact position in the FoV) or the atmosphere profile, or one telescope characteristics is varied. Possible uncertainty sources can be studied using each of these simulated samples.

5.4.1 Uncertainties on the shower track altitude

Beside the angular resolution described in section 5.3, the precision on the shower track geometry is related to the knowledge of the shower track altitude. As explained, two different methods may be applied to determine , either using the reflected Cherenkov signal, or when the information can not be retrieved, using the shower track length only. In both cases, the shower track zenith angle must be known to infer .

At the first order, according to geometrical considerations only, the impact on of an uncertainty on is estimated. Whereas an uncertainty induces an underestimation of when the Cherenkov information is used (relatively increasing with but less than 2% for shower with zenith angle ), it induces an overestimation of when the Cherenkov information is not used (relatively decreasing with but less than 2% for shower with zenith angle ).

The results are obtained simulating about 12600 events (of which 11200 reconstructed) in clear sky (US-Standard atmosphere) with LOWTRAN and UNISIM and a telescope aperture of 2.5 m (parametric optics).

ESAF results show that in clear sky and US-Standard Atmosphere, can be reconstructed with an error of 100-300 m for events with both fluorescence and Cherenkov echo components for , and with an error of 400-800 m when the Cherenkov echo is either missing or just ignored for . The absolute error grows with for both methods, while the relative error is rising for the method based on the Cherenkov echo (from 5% to 8%) and decreases for the method based on the time width of the p.e. distribution (from 22% to 5%). Therefore these methods can be considered as complementary.

5.4.2 Effects of geometric uncertainties

Uncertainties on the shower geometry affect the factors used in the expected light profile calculation. At each shower point , the reconstruction process involves the fluorescence yield , the light transmission , and a solid angle factor . ESAF allows to evaluate uncertainties induced on them related to shower geometry. A 1 km shift on altitude affects the solid angle computation by less than 1%. Despite the fact that the light transmission is the more affected factor by an uncertain altitude, an uncertainty on the reconstructed altitude of 1 km leads to less than 7% uncertainty on the direct transmission in clear sky. The fluorescence yield dependence with altitude is almost canceled by the dependence of the ionization energy deposit per meter of track length, thus an uncertainty of 1 km on altitude affects the fluorescence yield per meter by less than 1%.

Furthermore, the uncertainty on the shower geometry affects the slant depth calculation. The integral of the reconstructed profile is thus affected by an uncertainty of 10% for every km of error in altitude, and by around  g/cm. In the same way an uncertainty on shower zenith angle leads to an uncertainty on of around  g/cm for but to only 0.2% for on the integral of the shower profile.

5.4.3 Effects of atmospheric uncertainties

ESAF allows also to evaluate uncertainties on , and slant depth calculation related to atmosphere knowledge. During simulation runs, atmosphere parameters are varied according to all available MSISE profiles, and in the reconstruction process, only US Standard atmosphere model is considered. These studies show that the uncertainty on transmission and on the fluorescence yield per meter stays below 2% (respectively 2.5%). They also demonstrate that uncertainties induced by slant depth calculation lead to a shift on from to , and to an uncertainty on the reconstructed profile integral around 3% which could reach 4% for shower detected close to the equator.

5.5 Energy and resolution

Apart from the uncertainties on the atmospheric profile, different contributions described above are taken into account to perform preliminary studies of the relative resolution of energy for an EUSO detector with 2.5 m entrance pupil diameter. Showers are simulated with UNISIM, standard EUSO parameters are used for the detector and simulations have been performed in clear sky conditions with a US-Standard atmosphere, neglecting the multiple scattering effects. Results are displayed in Fig. 40 as a function of UHECP energy in various bins of zenith angle for proton initiated EAS.

Figure 40: Relative energy resolution as a function of energy in various bins of zenith angle for proton initiated EAS.

An energy resolution between % is expected at eV. It worsens at lower energy due to a low signal statistics and improves at higher energies. There is a strong dependence on the zenith angle of developing EAS, also due to the increasing number of collected photons as increases.

As already said, an uncertainty on the shower altitude induces an important uncertainty on the reconstructed , that is why results on resolution are presented depending on the method used for the first altitude estimation. For events with both fluorescence and Cherenkov signals can be reconstructed with resolution better than for . At larger angles the Cherenkov signal becomes suppressed, and the method based on the fluorescence profile must be used. In this case the resolution on is degraded, but it is still in the range between and , depending on the shower inclination.

Figure 41: The plot shows the error (in g/cm) with which is reconstructed when the Cherenkov pulse is available (lower curve, blue points) and when only fluorescence profile is used (upper curve, red points).

5.5.1 Effect of scattering on E and resolution

Simulation results show that photon scattering distorts the fluorescence signal. Studies described in section 5.1.3 are used to evaluate systematics related to the fluorescence signal contamination. In case of clear sky conditions, the total number of photons at the entrance pupil is increased by 10% to 25%, depending on the shower zenith angle (see Fig. 34). The contamination is mainly due to the scattered Cherenkov light. If this contribution is not subtracted to the signal, it leads to an uncertainty on in the same range, and to a shift of the maximum of the profile (in time); is then overestimated by around . Scattering of fluorescence photons as well as multi-scattering (more than once) of the Cherenkov photons have weaker impacts on the energy (less than 5%) and () resolutions. The subtraction of the diffused Cherenkov light is possible, as it was already done in HiRes experiment (see [68]).

5.6 Effect of clouds on the detected signal

5.6.1 Reflected peak

Because clouds and aerosols are responsible for Mie scattering, an adequate simulation including multi-scattering must be carried out. As examples Fig. 42 and 44 show simulation results at the pupil entrance of the same shower as in section 5.1 but with presence of rural aerosols or clouds. As one can see, contributions from scattering orders greater than 4 cannot be neglected as in the clear sky case. Therefore we carried out the simulations taking into account up to 8 scattering orders.

Figure 42: Number of photons reaching the pupil as a function of time. Signal integrated in the whole FoV. Top histograms refer to the last photon interaction and bottom histograms refer to the number of scattering interactions before reaching the pupil. The simulation has been carried out with 8 orders of scattering in a US-Standard atmosphere including a rural aerosol model of visibility 5 km. The energy of the shower is 10 eV and the zenith angle is 60.

In presence of aerosols, fluorescence photons being emitted mostly above aerosols layers, direct fluorescence light is only weakly affected. On the contrary, ground reflected light is strongly attenuated since it has to travel twice through the dense layers of atmosphere before being detected. Aerosols being mostly a scattering medium, a multi-scattering peak dominated by the contribution of Mie scattered Cherenkov photons appears in the total photon time distribution (Fig. 42). Its maximum is between time corresponding to ground impact and that corresponding to aerosol top layer impact. However most of these photons are not contiguous with the shower track. To estimate the shape of the multi-scattering peak likely to be detected, an appropriate selection of the scattered signal component at the photo-electron level is applied, using a simple detector simulation (already mentioned in section 5.1.3). Only the scattered signal which is contiguous in space and time with the shower track is selected, keeping in mind that selecting too many pixels reduces the signal to noise ratio.

This selection has been applied to the considered event and the resulting signal is shown on Fig. 43. As one can see the amplitude of the multi-scattering peak is weak compared to statistical fluctuations. Let us notice moreover that the peak maximum after pixel selection in space and time still does not correspond to the impact on the aerosol top layer. In a general way, its position in time should depend on the aerosols properties. Therefore the use of the multi-scattering peak (if detected) instead of ground reflected Cherenkov appears in the reconstruction procedure to be more difficult than expected.

Figure 43: Mean photoelectrons number as a function of time after selection of the shower track contiguous pixels. The energy of the shower is 10 eV and the zenith angle is 60. The atmosphere model is US-Standard with rural aerosols of visibility 5 km. The vertical bars indicate respectively the impact on aerosols top layer and the impact on ground.
Figure 44: Number of photons reaching the pupil as a function of time. Top histograms refer to the last photon interaction and bottom histograms refer to the number of scattering interactions before reaching the pupil. The simulation has been carried out with 8 orders of scattering in a US-Standard atmosphere including an homogeneous cumulus cloud model, characterized by a top altitude of 3 km, a vertical extension of 1 km and an optical depth of 2. The shower energy is 10 eV and the zenith angle is 60.
Figure 45: Mean photoelectron number as a function of time after selection of shower track contiguous pixels. The shower energy is 10 eV and the zenith angle is 60 degrees. The atmosphere model is US-Standard with an homogeneous cumulus cloud model, characterized by a top altitude of 3 km, a vertical extension of 1 km and an optical depth of 2. The vertical bars indicate respectively the impact on cloud layer and the impact on ground.

Contrary to aerosols, clouds can be found at higher altitudes, and thus distort strongly direct fluorescence light depending on their altitude, optical depth and the shower zenith angle. However, in some cases of low altitude clouds, when major part of fluorescence light (including the maximum of development) is emitted above clouds, direct fluorescence light is only weakly distorted. Nevertheless reflected Cherenkov is strongly attenuated and scattered photons are added to signal. As in case of aerosols, a multi-scattering peak dominated by Mie scattered Cherenkov photons appears in total photon time distribution (Fig. 44). Mie scattered photons are also spread in FoV therefore the amplitude of the multi-scattering peak is greatly reduced and low compared to statistical fluctuations when pixel selection is applied (Fig. 45). We have studied several configurations as a function of shower zenith angle, cloud top altitude, cloud thickness and cloud optical depth. In some particular situations, for dense clouds located at low altitude and showers weakly inclined (), the multi-scattering peak may be clearly visible. Two peaks may even appear in very specific situations, the first related to the impact on the cloud layer, the second related to the impact on ground. In a general way, the peak shape and its position in time depend strongly on clouds properties. Therefore, in the reconstruction procedure, the use of the multi-scattering peak (if detected) instead of ground reflected Cherenkov appears to be more difficult than expected.

An on board atmospheric sounding device associated with the telescope is thus required in order to interpret the signal distortions due to aerosols and cloud layers.

5.6.2 Effect of clouds on acceptance

Presence of clouds and aerosols in the FoV makes signal more complex to exploit and modifies also detector acceptance. Despite the fact that multi-scattering until order 8 is needed in cloud simulations and because the “Reduced” Monte-Carlo algorithm for radiative transfer is time-consuming, a compromise using the “Bunch” algorithm has been done for preliminary studies concerning acceptance with clouds. The main effect of clouds on the acceptance comes from the attenuation of the detected fluorescence light when a cloud layer stands between a shower and the detector. For this study, the scattering component has not been calculated, only the attenuation effect has been considered, using a simplified version of the “Bunch” algorithm.

Figure 46: Trigger efficiency curves in a modeled atmosphere with and without clouds.
Figure 47: Acceptance for triggering showers that have their maximum within the FoV as a function of the zenith angle in clear sky conditions (top) and cloudy conditions (bottom). The plot also shows the contribution of events with (in blue) and without (in red) Cherenkov peak. The shower energy spectrum is flat (in a logaritmic scale) between 10 eV and 10 eV. To guide the eyes, the curves giving an estimate of the geometric acceptance (showers having their maximum in the FoV) are also superimposed. Percentage are related to the simulated geometric acceptance whose integrated value is 633 000 km sr (the estimated one gives a 4% lower value than the simulated one).

Two samples of 210 proton EAS with flat energy distribution (in a logarithmic scale between 10 eV and 10 eV), isotropic angular distribution at top of atmosphere and a maximum of shower development in FoV have been used. Showers were generated with GIL parametrization in an atmosphere distribution obtained from MSISE model combined with the International Space Station Trajectory data. TOVS distribution and description of clouds (probability of presence, opacity and altitude) have been included in one of the two samples. The Earth surface is considered as a lambertian one, with an albedo of 5%. The fluorescence yield is based on the Kakimoto model completed with Bunner rays. The detector is derived from the EUSO Phase A design and the trigger efficiency is computed using the pattern [4;4].

In Fig. 46 the effect of clouds on the trigger efficiency is illustrated: it decreases when presence of clouds is taken into account. At the highest energies, between 10eV and 10eV, the mean efficiency is only 80% when clouds are taken into account instead of 96% in clear sky conditions.

In Fig. 47 we also show the trigger acceptance as a function of the shower zenith angle in clear sky conditions and in the presence of an atmosphere with clouds. In clear sky conditions the data nicely agree with the expectation that is at first order proportional to sin(2). The contribution to the total number of events of near horizontal proton initiated showers is negligibly small; however, the possibility to reconstruct and identify such showers might possibly open a new window for the detection of high energy neutrinos. We have not considered in this paper the problems related to neutrino detection.

The main effect of clouds, as already said, is that the global trigger acceptance is decreased from 96% to 80%, but also that events having a detectable Cherenkov peak is reduced from 79% down to 36%. To obtain this last result we have defined a “detectable Cherenkov peak” as a peak with maximum amplitude greater than 4 standard deviations above background.

In case of clouds presenc