# An analytical approach to the multiply scattered light in the optical images of the extensive air showers of ultra-high energies

###### Abstract

One of the methods for studying the highest energy cosmic rays is to measure the fluorescence light emitted by the extensive air showers induced by them. To reconstruct a shower cascade curve from measurements of the number of photons arriving from the subsequent shower track elements it is necessary to take into account the multiple scatterings that photons undergo on their way from the shower to the detector. In contrast to the earlier Monte-Carlo work, we present here an analytical method to treat the Rayleigh and Mie scatterings in the atmosphere. The method consists in considering separately the consecutive ’generations’ of the scattered light. Starting with a point light source in a uniform medium, we then examine a source in a real atmosphere and finally - a moving source (shower) in it. We calculate the angular distributions of the scattered light superimposed on the not scattered light registered from a shower at a given time. The analytical solutions (although approximate) show how the exact numerical results should be parametrised what we do for the first two generations (the contribution of the higher ones being small). Not allowing for the considered effect may lead to an overestimation of shower primary energy by and to an underestimation of the primary particle mass.

###### keywords:

ultra high energy extensive air showers, cosmic rays, fluorescence light, shower reconstruction^{†}

^{†}journal: Astroparticle Physics

## 1 Introduction

One of the methods for studying extensive air showers of high energies ( eV) is to register their images in the optical (mainly fluorescence) light. This can be done by observing showers from the side in order to avoid the more intense Cherenkov light emitted rougly in the shower direction. The observations are made by a number of optical telescopes, each containing a large mirror and a camera with a matrix of photomultipliers (PMTs) placed at the focus of the optical system (HiRes hires (), The Pierre Auger Observatory pao (), The Telescope Array tarr ()), so that photons arriving from a given direction on the sky are focused on a particular PMT (pixel). Photon arrival time can also be measured if the time structure of the PMT signals is recorded.
A cosmic ray induced shower produces at a given time a light spot on the camera which moves across it as the shower develops in the atmosphere so that succeeding PMTs are being hit.

It would have been ideal if the light producing a shower image had contained the fluorescence photons only. This is because there exists experimental evidence that the number of fluorescence photons induced by a charged electron in the atmosphere is proportional to the energy lost by it for ionisation dedx (). As practically all primary particle energy is eventually used for ionisation, this energy can be determined by measuring the fluorescence light emitted along the shower track in the atmosphere.

However, there are several problems in deriving the flux of the fluorescence light emitted by a shower from that arriving at the detector.
Firstly, the arriving light contains not only the fluorescence but also Cherenkov photons. If the viewing angle (the angle between the line of sight and the shower direction) is large (say, ) then it is mainly the Cherenkov light scattered in the atmosphere region just passed by the shower and observed by the detector (this light, before being scattered, travels roughly along the directions of the shower particles). Typically its fraction at the detector is about of the fluorescence flux. For smaller viewing angles it is the Cherenkov light produced at the observed part of the shower that may dominate even the fluorescence signal. The contribution of the Cherenkov light, which has to be subtracted from the total signal, has been extensively studied espec (); simil (); nerl (); latch ().

The subject of this paper is another phenomenon, affecting shower images, most commonly called the multiple scattering (MS) of light. Photons produced at the observed shower element, whatever their origin (fluorescence or Cherenkov), may undergo scattering in the atmosphere on their way from the shower to the detector, causing an attenuation of the light flux arriving at the detector and a smearing of the image. This scattering may take place on the air molecules (Rayleigh scattering) or on larger transparent particles, aerosols (Mie scattering). Most of the scattered photons change their directions, so that they no longer arrive at the pixel registering the not scattered (direct) light. Moreover, they arrive later having longer path lengths to pass. On the other hand, photons emitted by the shower at earlier times and scattered somewhere, may fall in the field of view of the pixels just registering the direct photons emitted at a later time. The net effect is that the scattered light forms its own instantaneous image superimposed on that in the direct light.

Our aim is to calculate the shower images in the multiply scattered light, so that this effect could be allowed for (subtracted) when determining the shower primary energy from the PMT signals. This problem was already studied by Roberts roberts (), by us in several, short conference contributions giller () and more recently by Pȩkala et al pekala (). The approach of the other authors was based on Monte Carlo simulations of
photons emitted by a shower. Photons were followed up to 5-6 scatterings and their arrival directions and time were registered by the detector. Many shower simulations were needed to obtain the MS images for various distances, heights, viewing angles of the observed shower parts. Finally, a phenomenological parametrisation of the number of MS photons was made as a function of the parameters found as relevant.

In contrast to Roberts and Pȩkala et al this approach is based on an analytical treatment. The main idea is to consider the arriving MS light as a sum of the photons scattered only once (the first generation), of those scattered two times (the second generation) and so on and calculate separately the angular and temporal distributions for each generation.

We start (Section 2) with a consideration of the simplest situation when a point source of isotropic light flashes for a very short time in an uniform medium. We derive analytical expressions for the angular and temporal distributions of the first and next generations of light arriving at a particular distance from the source.

As our aim is to apply our results to cosmic ray showers we need to consider a non-uniform medium like the atmosphere. Assuming an exponential distribution of the gas density and similarly for aerosols we show that an effective scattering length between any two points in the atmosphere can be easily calculated analytically. Signals of the first two generations arriving at a particular detector within a given angle to the direction to the source are found as a function of time (Section 3).

Using these it is straightforward to derive the corresponding distributions if the source moves across the atmosphere, integrating the point source distributions over changing distance and time of light emission.
In Section 4 we consider a moving light source, modelling a distant shower. We calculate angular distributions of MS light arriving at a detector at the same time as the direct (not scattered) photons emitted by the shower. This particular approach is quite natural because the data from optical detectors consist of the recorded signals by the camera PMTs within short time intervals so that one needs to know how much of the MS light has to be subtracted from the main, direct signal. The method for calculating images in the MS light simultaneous with images in the direct light is relatively simple (for the first and the second generations) as it is based on the geometry of the scattered photons in the particular generation. It does not require time-consuming Monte-Carlo simulations that were done for various shower-observer geometries.
Making some approximations we derived analytical formula for a shower image produced by the first generation (Section 4.2). Our analytical derivations allowed us to choose easily the variables on which and how to parametrise the MS signals. They also made us to realise that the dependence of it on the viewing angle was different for Rayleigh and Mie scatterings. Thus, we introduced a new, simple parametrisation of the fraction of the scattered photons arriving at the telescope within a given viewing cone, depending on the viewing angle of the shower (what has not been done before), separately for the two scatterings. For the same reason we also parametrised separately the second generation (Section 4.3).

A discussion of the results and the implication of the MS effect on the derivation of shower parameters is given in Section 5. The last Section (6) contains a summary and conclusions.

## 2 Point source flashing isotropically in uniform medium

At any fixed time a distant shower can be regarded as a point source emitting isotropically fluorescence light (about Cherenkov light see later). As explained above, we treat the light scattered in the medium as a sum of consecutive generations consisting of photons scattered only once, twice, and so on, on their way from the source to the observation point.

### 2.1 First generation

Let us consider the first generation, consisting of photons scattered once only. We shall calculate the flux of these photons, , at a distance from the source, such that is the number of photons scattered only once, arriving at time after the flash, within a solid angle at the surface (perpendicular to the arrival direction) located at a distance . To do this we shall calculate first the number of photons crossing the sphere of radius (from inside) at an angle () with respect to the normal, at time (); see Fig. 1. The average number of photons, per one photon emitted, interacting at a distance () from the source and scattered at an angle within , equals

(1) |

where is the mean scattering path length, is the probability that, once the scattering has occured, the scattering angle is within . To each pair of variables there corresponds another pair related to the former by

(2) |

and

(3) |

where and is the speed of light.

The Jacobian of the transformation gives

(4) |

Thus, we obtain

(5) |

Finally, the number of photons arriving at a unit surface at an angle (all azimuths) at time equals

(6) |

and

(7) |

#### 2.1.1 Rayleigh scattering

For the Rayleigh scattering we have

(8) |

Expressing as a function of and (Eq. 2) we obtain

(9) |

where . Thus, the flux of the first generation, defined above, equals

(10) |

where is the mean free path length for the Rayleigh scattering. However, in the exponent depends on all the scattering processes active. In general it is determined by

(11) |

for processes. Thus, if both molecular and aerosol scatterings are active but one wants to calculate the flux of photons scattered by the Rayleigh process only,
in the exponent equals
but in the denominator one has .

From Eq. 10 one can find the number of photons arriving per unit time at a unit surface within a given angle , as a function of time.

We have

(12) |

The integral can be found analytically, giving the result

(13) | |||||

where , , , , , ,
.

We have also calculated analytically a similar distribution if the scattering was isotropic. i.e. if (Appendix A).

One can also find analytically the angular distribution of the arriving light (integrated over time), but for small angles only (Appendix B). The result is

(14) | |||||

where , , if there is no Mie scattering and is the Euler constant.

The ratio of all photons arriving within a small angle to those not scattered , equals

where terms have been neglected.

#### 2.1.2 Mie scattering

In the next paragraph we shall consider the scattering of light emitted by showers developing in the real atmosphere, i.e. with the density depending on height. In addition to the Rayleigh process one has to take into account the Mie scattering occurring on particles (aerosols) larger that the light wavelength.
The Mie angular distribution is concentrated at rather small angles, in contrast to the Rayleigh case. Moreover,
in the deeper parts of the atmosphere the mean free path length for the Mie scattering may be comparable to that for Rayleigh, so that it is necessary to calculate the distribution of the light scattered by the Mie process only.

As before, we start with a simpler case - a uniform medium. The angular distribution of light scattered on particles with sizes larger
than the light wavelength depends on the distribution of the sizes and is not a well known function. Roberts roberts () adopts a function of the form

(16) |

Here, however, we prefer an expression allowing us to perform some integrations analytically. Most crucial is to have the number of numerical integrations for the second generation as few as possible. We shall see that to find for the Rayleigh scattering (Section 2.3) there is only one integration (over ) to be done numerically since the form of enables one to integrate analytically over and (Eq. 25 and 26). Thus, we adopt the following form for the Mie angular distribution:

(17) |

where , , .

This function is normalised as follows

(18) |

It describes quite reasonably the distribution used by Roberts.

From Eq. 6 and 7 we have

(19) |

where if

and if

Since

(20) |

we obtain

(21) |

where , and is determined as before. In principle, it is possible to find analytically the number of photons arriving within an angle after time per unit time. However, each of the nine integrals

(22) |

contains many terms itself, so that an analytical dependence on and/or would be practically lost.
Thus, we have found by integrating the flux (Eq. 21) numerically (similarly to Eq. 12).

Finally, the total flux of the first generation is the sum of the two fluxes arising from the two active mechanisms
of the scattering

(23) |

### 2.2 The second generation

These are the photons scattered exactly two times. We shall consider first a general case when there are more than one scattering processes (as Rayleigh and Mie). Let us call the process of the first scattering as and this of the second one as . Both and can be either Rayleigh or Mie. We denote their mean scattering path lengths by and , and the angular distribution functions of the scattering by and , correspondingly. The light source flashes isotropically at the centre of a sphere with radius at time (Fig. 2). As before we want to calculate the number of photons crossing the surface of the sphere from inside at a given angle , at time , per unit time.

Let us consider the photons scattered for the second time at a distance from the source. The number of photons incident on a small surface at an angle (within ) at time and scattered within a spherical shell of thickness by an angle
(within ) equals

(24) |

where . As now both processes are active the meaning of in the factor in the expression for is the effective mean path length for both processes.

The direction of the scattered photons is at an angle to the radius of the sphere and , where is the azimuth of the photons scattered for the second time. For any given direction before and after the second scattering, the scattering angle fulfils the relation
.
The only function depending on the azimuth angle of the incident photons is . Denoting

(25) |

and integrating (25) over we obtain for the number of photons incident on and scattered within into the solid angle the following expression

(26) | |||||

where the function is defined by the integral in the l.h.s. of (26).

The pair of fixed variables and defines another pair and , where is the photon path length after the second scattering. The Jacobian of the transformation gives the relation

(27) |

Putting , the contribution of photons scattered for the second time at a distance to arrive at an angle at the sphere with radius equals

(28) | |||||

The factor multiplied by in the expression for gives , independent of . Integration over gives the total number of the above photons

where .

The maximum value of results from fixing time . We have that

(30) |

Thus, the flux of the second generation equals

and, integrated over for , gives . With the Rayleigh and Mie processes active we must take into account all four cases or and or . Finally, the flux of the combined second generation photons is a sum of all specific fluxes

(32) |

Some of the integrals defined in this Section can be found as analytical functions (Appendix C). It is of some importance when calculating higher generations (see the next Section).

### 2.3 The next generations

Any next generation of the scattered photons can be calculated in the same way as the second one has been found from the previous one (the first). To calculate the flux of the generation, given we proceed as before when calculating the second generation from the first one (Eq. 24). The number of photons, incident on at an angle within at time and scattered along into equals:

(33) |

If there are two scattering processes, , then

(34) |

is the scattering probability by an angle by any process per unit distance per unit solid angle. The rest of the derivation of is the same as in the previous Section. However, for each next generation the number of numerical integrations increases, unless values of are stored as a 3-dimension matrix. Thus, it is convenient to find analytical solutions of the integrals and/or , if possible.

### 2.4 Results of calculations

Fig. 3a shows the number of photons arriving at the detector within an angle per unit area per unit as a function of .
The upper curves refer to the first generation, the lower - to the second one. We also compare here the time distributions obtained for the Rayleigh with those for the isotropic scattering. The distance detector- source equals to one scattering length ().

First of all we notice that for short times () the first generation dominates over the second one, and (as we can guess) over the higher ones. It can be seen from the formulae for the isotropic scattering (Appendix A) that the ratio for any given time should be proportional to , so that the importance of the second (and the higher) generation will be bigger for larger .

We can also see that the number of photons arriving within an opening angle reaches the dependence only at later times. This reflects the fact that the initial angular distribution of light is steep and becomes almost flat at times or so.
When comparing the Rayleigh curves with the isotropic ones one can see that the latter are slightly flatter for shorter times, as might be expected but become parallel to the former for later times.

Next, we consider a situation when there are two scattering processes, Rayleigh and Mie with quite different angular distributions (as discussed before). We adopt

(35) |

and for simplicity.

The result for the first and the second generation depending on time is shown in Fig. 3b. There is now more light at earlier times than in the previous case (Fig. 3a) due to the strong Mie scattering in the forward directions. However, the flux of the first generation decreases about 3 times quicker over the considered time region. Although the ratio of the second to the first generation is practically the same at in both cases . the importance of the second one is reached at later times when the Mie scattering is present.

## 3 Point light source in the atmosphere

Now we shall study the situation when a point light source flashes in a non-uniform medium, such as the atmosphere. We assume that the atmosphere is composed of two sorts of matter, molecules and aerosols, each having its density decreasing with height exponentially with a different scale heights, - for molecules and for aerosols, and having the corresponding mean path lengths for scattering at the ground and . It is not difficult to derive that the effective mean free path for a scattering for light travelling between two arbitrary points and equals

where and are the heights of points and above the level, where the Rayleigh and Mie scattering path lengths are correspondingly and . If the source (point ) is at a distance from the detector (point on the ground) then the ratio , of to the mean free path length along equals

(37) | |||||

It can be seen that increasing the distance to infinity (keeping constant) the ratio reaches its maximum finite value

(38) |

This situation is illustrated in Fig. 4. Here a vertical cross-section of the atmosphere is shown. Detector is at and lines represent constant values of corresponding to the straight path from the detector to the point on the line. We have adopted the following values: , and . These values describe approximately the atmospheric conditions at the Pierre Auger Observatory pao ().

From Fig. 4 one can also deduce that relevant values of , if light sources are at distances (as extensive air showers seen by Auger) are . Thus, this will be the region of our interest.

### 3.1 The first generation

As the medium is non-uniform, we cannot use the idea of a sphere to be crossed by the scattered photons, as in Section 2. Now their flux will depend on the zenith angle of the source and on the azimuth angle around the direction towards it. Fig. 5 shows a trajectory of a first generation photon scattered at . We want to calculate the angular distribution of the first generation as a function of time, for a fixed and , , crossing a unit area perpendicular to the direction towards the source.

We notice that for a fixed arrival direction of photons and time , the scattering point is uniquely determined. To arrive at the detector at angles within after time , photons have to cross the surface (shaded in the figure) and be scattered along a path length by the angle determined by Eq.2. The number of such photons equals

where , , is the angle between the normal to the surface and the direction of the incident photons , is the solid angle determined by the unit area at and the scattering point (). There is no need to calculate because , so that it cancels out. It can be shown that

(40) |

Inserting this into (40) we obtain

One can see that this formula is practically the same as (7) for the uniform medium, the only difference being in the scattering path lengths depending not only on distances but also on the geometry.

The height of the scattering point necessary to calculate and equals:

We calculate numerically the time distributions of light , arriving at the detector at angles smaller than for different zenith angles of the source. The results, in the form of the ratio:

(43) |

are presented in Figs. 7, 8 and 9, where are the distributions obtained in the previous section for a uniform medium.

To compare the light flux obtained for the real atmosphere with that for a uniform medium we adopt the same value of and the same distance from the source to detector for both cases. To understand the effect of a purely exponential atmosphere we start with considering only the Rayleigh scattering (Fig. 7), neglecting Mie (). Understanding the behaviour of the curves in this figure is easier with the help of Fig. 6. Each ellipse is a cross-section of a rotational ellipsoid with the symmetry axis determined by point - the detector and point - the light source. These are the focal points of all the ellipses. Photons arriving at the detector at after some fixed time (in units of the distance PD) must have been scattered on the surface of an ellipsoid with the eccentricity equal

(44) |

The ellipses refer to and keeping the right proportions. The detector field of view cuts only a part of the ellipsoid surface where the photons registered after time must have been scattered. We notice that all ratios in Fig. 7 are smaller than 1 and decrease with time (although those for are practically constant) and ratios for are larger than those for . All this becomes clear when inspecting Fig. 6 and the corresponding scattering sites. For example - the constancy of )
for times reflects the fact that the scattering takes place very close to point during all this time and starts to move away from it (higher in the atmosphere) only for larger times i.e . Since the scattering path length in the uniform medium is chosen equal to the effective path length in the atmosphere (Eq. 36), we have that and the scattering probability at is smaller in the exponential atmosphere.

Let us consider now a more realistic atmosphere with both processes, Rayleigh and Mie at work.
Fig. 8 shows time dependence of the ratio for and and zenith angles . Let us take a closer look at the case and (upper solid curve). It may seem strange that the ratio increases since in the case of the Rayleigh scattering only it decreases, although very slowly. We have checked that a similar slow decrease takes place if the scattering is only Mie. The behaviour of the curves in this figure would be hard to understand without the presentation of the scattering sites in Fig. 6. It can be seen that for smaller than a few the scattering sites are close to the source at . Thus, we may approximate the ratio as follows

(45) |

where are some effective angular distributions of photons scattered by Rayleigh (Mie) on the cut surface . As time increases, none of the changes much. However, the typical scattering angles increase, what affects much more than , so that decreases. Since

(46) |

the numerator decreases by a smaller factor than the denominator so that the ratio increases. It can be seen from Fig. 6 that for the scattering angles of the registered photons do not change much (the denominator stays constant) but now and have to be substituted by and , where is an effective scattering point with growing height. As , the ratio decreases and so does .

A different behaviour of can be explained also with the help of Fig. 6. At first decreases (the smaller - the stronger decrease) because the detector field of view cuts out a growing part of the deep atmosphere where the scattering is strong in the real atmosphere. At starts to increase for the same reason as just described in the case . It must finally decrease since the scattering takes places further and further behind the source, where and are growing in the real atmosphere.

In Fig. 9 we show for for several intermediate values of , and for two values and . Note that changing must result in changing - the distance to the source. This is the main reason why the curves for are lower than those for . It can be seen from Fig. 5 that for (the scales on both axes are the same so in the figure the angles are correctly represented) the distance is much shorter for than for , implying that the corresponding heights of the source differ considerably. Inspecting Fig. 6 we can estimate that for and the curves for should be down with respect to those for by a factor

(47) |

where and are positions of the source referring to and 1 respectively, with the meaning of all as defined before. In our example , , , so that the above factor equals whereas the exact ratio from Fig. 9 .

### 3.2 The second generation

We proceed like in the case of the first generation (Fig. 5), but now point refers to the second scattering. Photons scattered only once arrive at the surface from all directions according to , where . Thus, the number of photons incident on and scattered for the second time towards the detector (to arrive there within after time equals

(48) | |||||

where the integration has to be done over full solid angle , is the path length for the second scattering to occur () and is the angle of the second scattering. Now the pair of variables, and , does not determine uniquely the positions of the second scattering, since the times elapsed from photon emission to their arrival at (or strictly speaking, at ) have some distribution. However, can not be smaller than , thus . Expressing as a function of , and only we obtain

(49) |

and the distribution of the second generation equals

Finally, the number of photons within an angle equals

(51) |

The above integrals have been calculated numerically.

The ratio as a function of in the real atmosphere, where the source is at at the distance corresponding to , within two opening angles of the detector and , is presented in Fig. 10. For comparison we have also drawn there the same
ratio for the uniform medium, for with Rayleigh scattering active only (following from Fig. 3a for ). First we see that the contribution of the second generation increases with time as it should be expected. Next, it is considerably larger for the uniform medium than for the real atmosphere, which may not be so obvious at first sight. In fact, at the curve for the uniform medium is times higher than that for the real atmosphere. This factor equals ( for simplicity we neglect the operators )

The numerator as explained in the previous paragraph. The denominator refers to photons scattered exactly two times. For small the scatterings must take place close to the source or along the field of view ( in our example). Thus we should have that

(53) |

Our calculations show that for few so that , in agreement with the exact calculations.

## 4 A moving point source - a shower

So far we have been interested in the light signals from a flash of a stationary point source in a detector at a given distanc