Estimation of the number of muons with muon counters
aa r X i v : . [ a s t r o - ph . H E ] F e b Estimation of the number of muons with muon counters
A. D. Supanitsky a, ∗ a Instituto de Tecnolog´ıas en Detecci´on y Astropart´ıculas (CNEA, CONICET, UNSAM), Centro At´omico Constituyentes, San Mart´ın, 1650 Buenos Aires, Argentina.
Abstract
The origin and nature of the cosmic rays is still uncertain. However, a big progress has been achieved in recent years due tothe good quality data provided by current and recent cosmic-rays observatories. The cosmic ray flux decreases very fast withenergy in such a way that for energies & eV, the study of these very energetic particles is performed by using ground baseddetectors. These detectors are able to detect the atmospheric air showers generated by the cosmic rays as a consequence of theirinteractions with the molecules of the Earth’s atmosphere. One of the most important observables that can help to understand theorigin of the cosmic rays is the composition profile as a function of primary energy. Since the primary particle cannot be observeddirectly, its chemical composition has to be inferred from parameters of the showers that are very sensitive to the primary mass.The two parameters more sensitive to the composition of the primary are the atmospheric depth of the shower maximum and themuon content of the showers. Past and current cosmic-rays observatories have been using muon counters with the main purpose ofmeasuring the muon content of the showers. Motivated by this fact, in this work we study in detail the estimation of the number ofmuons that hit a muon counter, which is limited by the number of segments of the counters and by the pile-up e ff ect. We consideras study cases muon counters with segmentation corresponding to the underground muon detectors of the Pierre Auger Observatorythat are currently taking data, and the one corresponding to the muon counters of the AGASA Observatory, which stopped takingdata in 2004. Keywords:
Cosmic rays, Chemical Composition, Muon counters
1. Introduction
The cosmic ray energy spectrum extends over several ordersof magnitude in energy. The highest primary energies observedat present are of the order of 10 eV. Above ∼ eV thecosmic ray flux is so small that these very energetic particlesare studied by means of ground-based detectors, which are ableto detect the atmospheric air showers generated by the cosmicray interactions with the molecules of the atmosphere. Sincethe primary particle is not observed directly, its energy, arrivaldirection, and chemical composition have to be inferred fromthe shower information obtained by the detectors.The origin of the cosmic rays is still unknown. The threemain observables used to study their nature are: The energyspectrum, the distribution of their arrival directions, and thechemical composition. The composition of the primary parti-cle has to be inferred from di ff erent properties of the showers.The most sensitive parameters to the primary mass are the at-mospheric depth at which the shower reaches its maximum de-velopment and the muon content of the showers [1, 2, 3]. Ingeneral, the muon density at a given distance to the shower axisis used in composition analyses.The composition profile as a function of energy is very im-portant to understand several aspects of the cosmic-ray physics.In particular, at the highest energies the composition infor-mation plays an important role to find the transition between ∗ Corresponding author: [email protected] the galactic and extragalactic components of the cosmic rays[4, 5] and to elucidate the origin of the suppression observedat ∼ . eV [6]. The composition analyses are subject tolarge systematic uncertainties originated by the lack of knowl-edge of the hadronic interactions at the highest energies (seefor instance [7]). The composition is determined by com-paring experimental data with simulations of the atmosphericshowers and the detectors (when it corresponds). The showersare simulated by using high-energy hadronic interaction mod-els that extrapolate low-energy accelerator data to the highestenergies. This practice introduces large systematic uncertain-ties even when models updated using the Large Hadron Col-lider data are considered. Moreover, experimental evidence hasbeen found recently about a muon deficit in shower simulations[8, 9]. Even though this is an important limitation for com-position analyses, it is expected that mass-sensitive parameters,obtained with the next generation of high-energy hadronic inter-action models, present smaller di ff erences allowing for a reduc-tion of the systematic uncertainties introduced by those models.Past and current cosmic-rays experiments have been measur-ing muons by using di ff erent types of detectors [9]. A particularclass of detector is the muon counter. This type of detectors hasbeen used in the past in the Akeno Giant Air Shower Array(AGASA) [10] and at present in the Pierre Auger Observatory[11]. The muon counters are designed to count muons through asegmented detector. The segments of the Auger muon countersare scintillator bars whereas the segments of the AGASA muoncounters were proportional counters. The limitation to measure Preprint submitted to Astroparticle Physics February 9, 2021 given number of incident muons is given by the number ofsegments of the counters.In general, the principle of operation of the muon counters isbased on a binary logic in which each channel of the electronics,associated to a given segment of the detector, is able to di ff er-entiate between a state in which the signal is larger than a giventhreshold level and the one in which it is smaller. The thresholdlevel is chosen in such a way that almost all muons can be iden-tified, i.e. the e ffi ciency of each segment is close to 100 %. Forthe case in which the signal is larger than the threshold level, thesegment is said to be on and otherwise o ff . The time structureof the signal corresponding to one muon limits the time intervalin which it is possible to identify single muons. This leads tothe definition of a time interval, usually called inhibition win-dow, in which it is decided whether a given segment is on or o ff . As a consequence, if one or more muons hit the same seg-ment in a time interval of the order of the one corresponding tothe inhibition window, the segment is tagged as on , losing theinformation about the number of muons that hit that segment.Therefore, when a given number of muons hit a muon counterin a time interval of the order of the one corresponding to theinhibition window, a number k of segments on is obtained. Ifthe number of incident muons is much smaller than the numberof segments, the random variable k is close the the number ofmuons that hit the counter. However, if the number of incidentmuons is close to the number of segments, the variable k be-comes much smaller than the number of incident muons. Thise ff ect is known as pile up [1].In this work, we find an analytic expression for the distribu-tion function of k , the number of segments on , given the num-ber of incident muons, n µ . In this case, k is the random variableand n µ is taken as a parameter of the distribution. The expres-sions for the mean value and the variance of k are inferred byusing the new formula, these expressions are equal to the onesobtained in Ref. [1] by using a di ff erent approach. We alsostudy how to estimate the parameter n µ from measured valuesof k and how to obtain a confidence interval. These studies aredone for 192 segments, which correspond to the total numberof segments of the Auger muon counters and 50 segments thatcorrespond to the number of segments of muon counters usedin AGASA.It is worth mentioning that the main purpose of the muoncounters is the reconstruction of the muon lateral distributionfunction (MLDF), i.e. the muon density as a function of thedistance to the shower axis, which is proportional to the meanvalue of n µ . Even though the estimation of the mean value of n µ ,which is studied in Ref. [12], is closely related to the determi-nation of the MLDF, the estimation of n µ is also important fordi ff erent types of applications. In particular, it is important forthe method used to reconstruct the MLDF developed in [1], inwhich an estimator of the number of muon that hit a given muoncounter is inserted in the Poisson likelihood that approximatesthe exact likelihood in a given range of n µ . Also, the estimationof n µ is necessary to obtain the calibration curve of the integra-tor (a complementary acquisition mode that muon counter canhave), which is given by the mapping of the number of incidentmuons into the integrated signal [11]. The estimation of n µ is also relevant for studies related to the signal fluctuations andsystematic uncertainties performed by using twin muon coun-ters [13].
2. Distribution function of the number of segments on andmuon number estimation The mean value of the number of muons, n µ , that hit a muoncounter is given by, λ = A ρ µ cos θ, (1)where A is the area of the muon counters, ρ µ is the muon densityat a given distance to the shower axis, and θ is the zenith angleof the shower. Since the muon counters sample the MLDF at agiven position in the shower plane, the number of muons thathit a given muon counter is a random variable that follows thePoisson distribution. The distribution of k given λ has beenobtained in Ref. [12] and is given by P ( k | λ ) = n s k ! exp( − λ ) (cid:2) exp ( − λ/ n s ) − (cid:3) k , (2)where n s is the number of segments of the muon counter. Notethat this distribution does not depend on the number of muonsthat hit the muon counters, which means that P ( k | λ ) corre-sponds to the marginalization of the joint distribution P ( k , n µ | λ )with respect to n µ .On the other hand, since all segments are equal, the probabil-ity distribution function corresponding to a given configurationof the number of muons that hit each segment in a time intervalcorresponding to one inhibition window is given by the multi-nomial distribution, P ( n , . . . , n n s ) = n µ ! n ! . . . n n s ! n s ! n µ , (3)where n i is the number of muons that hit the i − th segment. Here { n , . . . , n n s } are such that P n s i = n i = n µ .The total number of segments on can be written in the fol-lowing way, k = n s X i = e Θ ( n i ) , (4)where e Θ ( n ) = n = e Θ ( n ) = n ≥
1. Note that k ∈ N .The mean value and the variance of k can be calculated byusing Eq. (3), as reported in Ref. [1]. The mean value and thevariance of k are given by (see Appendix A for details of thederivation), h k i ( n µ ) = n s " − − n s ! n µ , (5)Var [ k ] ( n µ ) = n s − n s ! n µ " + ( n s − − n s − ! n µ − n s − n s ! n µ . (6)Beyond the fact that the mean value and the variance of k canbe calculated from Eq. (3) in a relatively direct way, an analytic2xpression of the distribution function of k given n µ , P ( k | n µ ), ismore di ffi cult to obtain. However, it can be inferred in a quitestraightforward way by using the distribution P ( k | λ ) given inEq. (2). As mentioned before, P ( k | λ ) is obtained marginalizing P ( k , n µ | λ ) with respect to n µ . Besides, P ( k , n µ | λ ) is given by theproduct of P ( k | n µ ) with the Poisson distribution. Therefore, P ( k | λ ) = ∞ X n µ = P ( k | n µ ) exp( − λ ) λ n µ n µ ! = n s k ! exp( − λ ) (cid:2) exp ( − λ/ n s ) − (cid:3) k . (7)After introducing the following expression in Eq. (7), (cid:2) exp ( − λ/ n s ) − (cid:3) k = k X j = kj ! ( − j exp ( λ ( k − j ) / n s ) , exp ( λ ( k − j ) / n s ) = ∞ X n µ = λ n µ n µ ! k − jn s ! n µ , (8)the form of P ( k | n µ ) is obtained comparing the terms propor-tional to λ n µ / n µ ! in both sides of Eq. (2), P ( k | n µ ) = n s k ! S ( n µ , k ) k ! n n µ s . (9)Here S ( n µ , k ) is the Stirling number of second kind, which isgiven by S ( n µ , k ) = k ! k X j = kj ! ( − j ( k − j ) n µ . (10)The random variable k ranges from 1 to n s . In the case inwhich n µ < n s and k > n µ , the condition P ( k | n µ ) = S ( n µ , k ) = n µ < n s and k > n µ . For that purpose, let us consider thefollowing expression of the Stirling number of second kind, S ( n µ , k ) = ( − k k ! k X i = ki ! ( − i i n µ . (11)From this expression, it is easy to see that, S ( n µ , k ) = ˆ O ( n µ ) f k ( x ) (cid:12)(cid:12)(cid:12) x = , (12)where ˆ O ( n µ ) corresponds to the di ff erential operator ˆ O = x d / dx applied n µ times to the function f k ( x ) = (1 − x ) k . By using thedefinitions given above, the following expression is obtained,ˆ O ( n µ ) f k ( x ) = n µ X i = a i x i d i f k dx i ( x ) , (13)where a i are constant numbers. If k > n µ , all derivatives inEq. (13) are proportional to a non-null and positive power of(1 − x ). Therefore, it follows that ˆ O ( n µ ) f k ( x ) | x = =
0, whichproves that S ( n µ , k ) = k > n µ and then P ( k | n µ ) = k > n µ . Since P ( k | n µ ) is a distribution function, it has to be normal-ized. To see that this is true for the expression in Eq. (9), letus consider the following property of the Stirling numbers ofsecond kind [14], x n µ = n µ X k = S ( n µ , k ) ( x ) k , (14)where ( x ) k = x ( x − ... ( x − k +
1) is the falling factorial. If x = n s , Eq. (14) becomes, n n µ s = n µ X k = S ( n µ , k ) ( n s ) k n µ ≤ n sn s X k = S ( n µ , k ) ( n s ) k + n µ X k = n s + S ( n µ , k ) ( n s ) k n µ > n s . (15)Using that S ( n µ , k ) = k > n µ and that ( n s ) k = k > n s ,Eq. (15) can be written as n n µ s = n s X k = n s k ! S ( n µ , k ) k ! , (16)which implies that n s X k = P ( k | n µ ) = n s X k = n s k ! S ( n µ , k ) k ! n n µ s = . (17)It is worth mentioning that from Eq. (9) it is possible to cal-culate the mean value and the variance of k which have to beequal to the expressions obtained by using the multinomial dis-tribution (see Eqs. (5) and (6)). The mean value of k is given by h k i = n s X k = n s k ! S ( n µ , k ) k k ! n n µ s . (18)From the recurrence satisfied by the Stirling numbers of secondkind, S ( n µ , k ) = k S ( n µ − , k ) + S ( n µ − , k −
1) [14], it followsthat k S ( n µ , k ) = S ( n µ + , k ) + S ( n µ , k − h k i = n s X k = n s k ! S ( n µ + , k ) k ! n n µ s − n s n s − X j = n s − j ! S ( n µ , j ) j ! n n µ s , (19)where the second term of this equation is obtained by doing thechange of variable j = k −
1. From Eq. (16), it is easy to seethat h k i = n n µ + s / n n µ s − n s ( n s − n µ / n n µ s , which after some algebrabecomes h k i = n s (1 − (1 − / n s ) n µ ), i.e. the same expression ob-tained by using the multinomial distribution. In a similar way,but in this case using the recurrence relation satisfied by theStirling number of second kind twice, it is possible to calculatethe variance of k , whose expression obtained in this way is thesame as the one obtained by using the multinomial distribution.The Auger muon counters installed in each position of thearray are composed of three modules of 64 segments each sum-ming a total of 192 segments (see Ref. [11] for details). Figure1 shows P ( k | n µ ) as a function of k for di ff erent values of n µ . As3
20 40 60 80 100 120 140 160 180 k ) µ P ( k | n = 20 µ n = 50 µ n = 100 µ n = 190 µ n = 350 µ n = 1000 µ n Figure 1: P ( k | n µ ) as a function of k for n s = k = expected, the maximum of the distribution is shifted towardslarger values of k for increasing values of n µ . Note that even for n µ = P ( k | n µ ) presents a maximum.As a result of a measurement, a given value of k is ob-tained. From this value of k it is possible to estimate the pa-rameter n µ and to determine a confidence interval at a givenconfidence level. The maximum likelihood estimator of n µ , ˆ n µ ,is obtained by finding the maximum of the likelihood function L ( n µ ) = P ( k | n µ ), i.e.,ˆ n µ = arg max n µ ∈ N P ( k | n µ ) . (20)In this case, ˆ n µ has to be calculated numerically for each partic-ular value of k measured. Figure 2 shows the likelihood func-tion as a function of n µ for k =
100 and k = k =
100 the likelihoodpresents a well-defined maximum as for all other allowed val-ues of k except for k = n µ → ∞ . This means that when the value k = n µ can be found, as shown below. Note that P ( k = n s | n µ ) → n µ → ∞ .Figure 3 shows ˆ n µ as a function of k obtained by maximizingthe likelihood function for n s = n µ (cid:27) k for small values of k . Moreover, ˆ n µ starts todeviate from k in more than 10 % at k (cid:27)
40. For larger valuesof k , ˆ n µ starts to increase faster in such a way that ˆ n µ = k = n µ can befound by using the expression corresponding to the mean valueof k given in Eq. (5). Inverting Eq. (5), the expression e n µ = ln − kn s ! ln − n s ! , (21)
100 120 140 160 180 µ n ) µ P ( k = | n
200 400 600 800 1000 1200 1400 1600 1800 2000 µ n ) µ P ( k = | n Figure 2: Likelihood function L ( n µ ) = P ( k | n µ ) as a function of n µ for n s = k =
100 (top panel) and k =
192 (bottom panel). The solid lines joiningthe discrete points are added to guide the eye. is obtained, where the mean value of k is replaced by the ran-dom variable k . The bottom panel of Fig. 3 shows the di ff erencebetween the maximum likelihood estimator of n µ , ˆ n µ , and theapproximated expression of Eq. (21). It can be seen that e n µ dif-fers in less than one from ˆ n µ in the whole range of the variable k , which shows that e n µ is a very good approximation of ˆ n µ .The confidence belt of the distribution function P ( k | n µ ) isconstructed by finding the values of k , for a given n µ , that satisfy X k P ( k | n µ ) ≥ − α, (22)where 1 − α is the confidence level (CL). In this work, the order-ing proposed by Feldman and Cousins [15] is considered for thecalculation of the confidence belt. The top panel of Fig. 4 showsthe confidence belt of P ( k | n µ ) for n s =
192 and 1 − α = . n µ as a function of k is also shown. From the figure it canbe seen that the width of the confidence belt increases with k , asexpected. It can also be seen that for k = n s =
192 only a lowerlimit can be obtained since the case k = n s , i.e., all segmentsof the counter on is compatible with a semi-infinite set of n µ values at any confidence level.4
20 40 60 80 100 120 140 160 180 k µ n k − − − µ n ~ - µ n Figure 3: Top panel: Maximum likelihood estimator of n µ , ˆ n µ , as a function of k . The solid line joining the discrete points is added to guide the eye. Bottompanel: Di ff erence between ˆ n µ and e n µ (approximated expression of ˆ n µ given byEq. (21)) as a function of k . The number of segments considered is n s = The relative uncertainty of ˆ n µ is defined here as ε = n max µ − n min µ n µ , (23)where n max µ and n min µ are the maximum and the minimum valuesof n µ , respectively, for a given value of k obtained at a givenCL. The bottom panel of Fig. 4 shows ε as a function of ˆ n µ . Asexpected, ε increases with ˆ n µ taking values smaller than ∼
16 %for ˆ n µ .
201 and reaching values of the order of 30 % in theregion where 700 . ˆ n µ .
900 (close to k = λ , given by Eq. (1), is af-fected by the segmentation of the detector in combination withthe pile-up e ff ect and also by the Poisson fluctuations. Figure5 shows the relative uncertainty corresponding to the estimatorof the mean value of the Poisson distribution as a function of n ,a measured valued of a Poisson random variable, which in thiscase coincides with the maximum likelihood estimator of themean value. This relative uncertainty is obtained following thesame procedure used to calculate the relative uncertainty corre-sponding to the ˆ n µ estimator. Comparing the bottom panel ofFig. 4 with Fig. 5 it can be seen that for values of ˆ n µ smaller than k µ n µ nConfidence belt0 200 400 600 800 1000 µ n [ % ] ε Figure 4: Top panel: Confidence belt of P ( k | n µ ) and ˆ n µ as a function of k .Bottom panel: ε as a function of ˆ n µ . The number of segments considered is n s =
192 and the confidence level used in the calculation is 1 − α = . ∼
113 the uncertainty on the determination of the mean value ofthe number of muons is dominated by the Poisson fluctuationsbut for values of ˆ n µ larger than ∼
113 the dominant uncertaintyis the one introduced by the segmentation of the detector incombination with the pile-up e ff ect.The AGASA experiment used muon counters of 50 segmentsto measure the muon content of the showers. Figure 6 showsthe distribution function P ( k | n µ ) as a function of k for di ff er-ent values of n µ and for n s =
50. As expected, for n µ = k =
50, whichmeans that these counters saturate with a much smaller numberof incident muons compared to the ones with 192 segments.The top panel of Fig. 7 shows the maximum likelihood es-timator of n µ calculated numerically by using Eq. (20) for n s =
50. As for the n s =
192 case, ˆ n µ (cid:27) k for small valuesof k in such a way that for k >
20, the departure of ˆ n µ fromˆ n µ (cid:27) k becomes larger than 10 %. For larger values of k , ˆ n µ starts to increase faster in such a way that ˆ n µ =
194 for k = n µ − e n µ as a function of k .Also in this case, the absolute value of this di ff erence is smallerthan one, which indicates that e n µ , given in Eq. (21), is a verygood approximation of ˆ n µ also for the n s =
50 case.The top panel of Fig. 8 shows the confidence belt of P ( k | n µ )5
200 400 600 800 1000 n [ % ] ε Figure 5: Relative uncertainty corresponding to the maximum likelihood esti-mator of the mean value of the Poisson distribution as a function of n , a mea-sured value of a Poisson random variable. The confidence level used in thecalculation is 1 − α = . k ) µ P ( k | n = 5 µ n = 15 µ n = 33 µ n = 70 µ n = 250 µ n Figure 6: P ( k | n µ ) as a function of k for n s =
50. The vertical dashed linecorresponds to k =
50. The solid lines joining the discrete points are added toguide the eye. for n s =
50 and 1 − α = . n µ as a function of k is alsoshown in the plot. As in the case corresponding to n s = k .This behavior becomes more evident from the plot in the bottompanel of the same figure, in which it can be seen that the relativeuncertainty of ˆ n µ is smaller than ∼
16 % for k ≤
67 and takesvalues close to 24 % in the region where 140 . k . n µ smaller than ∼
55 the uncertainty on the determination of the mean value ofthe number of muons is dominated by the Poisson fluctuationsbut for values of ˆ n µ larger than ∼
55 the dominant uncertaintyis the one introduced by the segmentation of the detector incombination with the pile-up e ff ect.Note that the relative uncertainty on the determination of ˆ n µ for the AGASA muon counters cannot be compared straight-forwardly with the one corresponding to Auger, since the areaof the AGASA muon counters is 25 m whereas the one corre- k µ n k − − − µ n ~ - µ n Figure 7: Top panel: Maximum likelihood estimator of n µ , ˆ n µ , as a function of k . The solid line joining the discrete points is added to guide the eye. Bottompanel: Di ff erence between ˆ n µ and e n µ (approximated expression of ˆ n µ given inEq. (21)) as a function of k . The number of segments considered is n s = sponding to the Auger muon counters is 30 m . Therefore, thenumber of muons that hit an AGASA muon detector is, on av-erage, ∼
17 % smaller than the one corresponding to an Augermuon detector, provided that the muon flux that hit the detectorsis the same. In any case, the determination of the muon densityat a given distance to the shower axis done by using the Augermuon detectors should be better than the one corresponding tothe AGASA muon detectors, since the Auger muon detectorshave a larger area and a larger number of segments ( A =
30 m and n s = A = and n s = ff erent e ff ects that have to be taken into account in ordernot to introduce biases in the estimated number of muons. Forinstance, there are three main e ff ects that can introduce biasesin the Auger muon detectors. The first one is the noise producedby the dark rate of the silicon photomultipliers, which can reachthe discriminator threshold due to the inner-cells crosstalk. Thise ff ect can be mostly reduced choosing a proper counting strat-egy (see Ref. [11] for details). The second one is the e ffi ciencyof each segment, which can be smaller than 100 %. In this6
10 20 30 40 50 k µ n µ nConfidence belt0 20 40 60 80 100 120 140 160 180 200 µ n [ % ] ε Figure 8: Top panel: Confidence belt of P ( k | n µ ) and ˆ n µ as a function of k .Bottom panel: ε as a function of ˆ n µ . The number of segments considered is n s =
50 and the confidence level used in the calculation is 1 − α = . case a correction can be obtained from the estimated e ffi ciency,which is measured in the laboratory [11]. Note that the e ffi -ciency of the Auger muon detectors is ∼ . ff ect can be estimated form detailedsimulations of the detector [16]. Note that the sources of biaseson the estimation of the number of incident muons depend onthe specific design of the muon detector under considerationand have to be studied in detail for each particular case.
3. Conclusions
In this work we have studied in detail the estimation of thenumber of muons that hit a muon counter from the number ofsegments on , k , which is the random variable that is measuredin an experiment. For that purpose we have found an analyticexpression for the distribution function of k , given a number ofincident muons. We have considered the number of segmentscorresponding to the muon counters of Auger and also the onecorresponding to the muon counters of AGASA.We have found that for small values of k , compared with thenumber of segments, the estimator of the muon number is closeto k but increases much faster for larger values of k . We have also found that the relative uncertainty in the determination ofthe number of muons is small for small values of k and that itincreases relatively fast with k reaching values close to 24 and30 %, for n s =
50 and n s =
192 respectively, in the region where k is close to the total number of segments.The main motivation of these studies is the measurement ofthe muon content of air showers initiated by cosmic rays, whichis intimately related to the chemical composition of the pri-mary particle, an open problem of the high-energy astrophysics.However, it is worth mentioning that the methods developed inthis work can be relevant in other applications. Appendix A. Calculation of h k i and Var[ k ] from the multi-nomial distribution In this section, the steps that lead to the expressions for h k i and Var[ k ] (see Eqs. (5) and (6)), obtained in Ref. [1] by usingthe multinomial distribution are given.Let us start with the calculation of the mean value of k . FromEq. (4) it can be seen that, h k i = * n s X i = e Θ ( n i ) + = n s De Θ ( n ) E . (A.1)The distribution function of n is given by the binomial distri-bution, i.e. P ( n ) = n µ n ! n s ! n − n s ! n µ − n . (A.2)Therefore, h k i = n s n µ X n = e Θ ( n ) P ( n ) = n s n µ X n = P ( n ) = n s (1 − P (0)) . (A.3)Combining Eqs. (A.2) and (A.3), the expression for the meanvalue of k given by Eq. (5) is obtained.The variance of k is calculated in a similar way. For thatpurpose, let us first calculate the mean value of k , which isgiven by h k i = * n s X i = n s X j = e Θ ( n i ) e Θ ( n j ) + = n s De Θ ( n ) E + n s ( n s − De Θ ( n ) e Θ ( n ) E . (A.4)To calculate the averages in Eq. (A.4) besides P ( n ), P ( n , n )is required. From Eq. (3) it can be seen that P ( n , n ) = n µ ! n ! n !( n µ − n − n )! n s ! n + n − n s ! n µ − n − n , (A.5)where n + n ≤ n µ . In a similar way to the one followed toobtain Eq. (A.3), the next expression for the mean value k isobtained h k i = n s (1 − P (0)) + n s ( n s − × − P (0 , − n µ X n = P ( n , . (A.6)7y using that, n µ X n = P ( n , = − n s ! n µ − − n s ! n µ , (A.7)the following expression is obtained, h k i = n s − n s (2 n s − − n s ! n µ + n s ( n s − − n s ! n µ . (A.8)From Eqs. (A.8) and (5), the expression for the variance of k given by Eq. (6) is obtained. Acknowledgements
A. D. S. is member of the Carrera del Investigador Cient´ıficoof CONICET, Argentina. This work is supported by ANPCyTPICT-2015-2752, Argentina. The author thanks the membersof the Pierre Auger Collaboration, specially C. Dobrigkeit forreviewing the manuscript.