text:
# 1 Unfolding correction factor

February 5, 2001

A Compilation of High Energy Atmospheric
Muon Data at Sea Level

Thomas Hebbeker

Institut für Physik der Humboldt-Universität zu Berlin,

Invalidenstr. 110, D-10115 Berlin, Germany

Charles Timmermans

Katholieke Universiteit Nijmegen,

Toernooiveld 1, N-6525 ED Nijmegen, Netherlands

ABSTRACT

We collect and combine all published data on the vertical atmospheric muon flux and the muon charge ratio for muon momenta above 10 GeV. At sea level the world average of the momentum spectra agrees with the flux calculated by E.V. Bugaev et al. within 15 %. The observed shape of the differential flux versus momentum is slightly flatter than predicted in this calculation. The experimental accuracy varies from 7% at 10 GeV to 17% at 1 TeV. The ratio of fluxes of positive to negative muons is found to be constant, at a value of $`1.268`$, with relative uncertainties increasing from approximately $`1\%`$ at low momenta to about $`6\%`$ at 300 GeV.

1. Introduction

We collect measured atmospheric muon flux data and charge ratios as a function of momentum and compute world averages. Only measurements at sea level or low altitudes and for (near) vertical incidence are taken into account, since several data sets are available for these experimental conditions. Here we consider only data with muon momenta above 10 GeV. At lower momenta geomagnetic effects and solar influences play a significant role and make the interpretation of the data more difficult. A recent compilation of charge ratio data at low muon momenta can be found in reference .

A precise knowledge of the muon spectrum and charge ratio allows to constrain the primary flux and the models of atmospheric showers so that also the atmospheric neutrino fluxes can be calculated with a good precision. This is a very important issue, since the Superkamiokande experiment and others have seen indications for a disappearance of atmospheric muon neutrinos. So far this interpretation is based on the angular distribution and on the ratio of muon neutrino to electron neutrino fluxes. It is very important to compare also directly the measured and calculated absolute muon neutrino fluxes; until now this was prevented by the large model uncertainties.

2. Effects relevant for spectrum and charge ratio

The following effects might influence the measurements of the muon flux and the charge ratio. It is possible that the published data need to be corrected accordingly in order to arrive at a meaningful comparison between the various measurements.

* Geomagnetic effects
For near vertical incidence the geomagnetic cutoff for primary protons is below $`10\mathrm{GeV}`$ for all latitudes at which the cosmic ray measurements were made (exceptions are discussed below). Geomagnetic effects can therefore be neglected.
* Solar modulation
Using the parameterization given in reference  we estimate that the primary proton flux at 50 GeV (100 GeV) decreases by $`3\%`$ (1.6 %) at maximum solar activity compared to the minimum. The mean primary proton momentum resulting in 10 GeV muons at sea level exceeds $`100\mathrm{GeV}`$. Using the air shower program CORSIKA , we found that about $`80\%`$ of those protons have a momentum larger than $`50\mathrm{GeV}`$. This results in an uncertainty of $`\pm 1\%`$ for the muon flux at a momentum of 10 GeV. Similarly, one can estimate a flux uncertainty of $`\pm 0.5\%`$ at 20 GeV and less at higher momenta. At 10 GeV the charge ratio is expected to change by about $`\pm 0.2\%`$. At higher momenta the effect is even smaller. We do not correct the data for time dependent solar effects.
* Altitude dependence
Not all experiments measure at sea level. In order to investigate the dependence of flux and charge ratio on the altitude we used the air shower simulation program CORSIKA and also apply the empirical formula found by De Pascale et al..

For muon momenta above 10 GeV and altitudes less than about 1000 m the vertical muon flux can be parameterized by

$`{\displaystyle \frac{\mathrm{\Phi }(h)}{\mathrm{\Phi }(h=0)}}=e^{h/L}\pm 0.003`$ (1)
where $`h`$ = altitude, $`L=4900\mathrm{m}+750\mathrm{m}\frac{p}{\mathrm{GeV}}`$ and $`p`$ = muon momentum.

The form of the parameterization is similar to the one used in. The uncertainty of $`\pm 0.003`$ reflects the quality of the parameterization and the comparison to the measurements. Example: For $`h=1000\mathrm{m}`$ and $`p=10\mathrm{GeV}`$ we obtain the flux $`\mathrm{\Phi }(h)=1.08\mathrm{\Phi }(h=0)`$. Note: Caprice data disagree with both  and CORSIKA for higher momenta; here they have not been taken into account.

The charge ratio is not affected, it changes by less than $`0.005`$ for $`h<1000\mathrm{m}`$ and $`p>10\mathrm{GeV}`$.

We do correct all published fluxes using formula (1).
* Zenith angle dependence
The muon data are normally collected within a certain cone around the vertical direction, including zenith angles up to $`\theta ^{max}`$. With help of CORSIKA we find that the zenith angle dependence can be parameterized in the form

$`{\displaystyle \frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}\mathrm{cos}\theta }}1+a(p)(1\mathrm{cos}\theta )`$ (2)
with a momentum dependent coefficient $`a(p)`$. Accordingly we estimate the following flux reduction factors

$`g(\theta )={\displaystyle \frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}\mathrm{cos}\theta }}(\theta )/{\displaystyle \frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}\mathrm{cos}\theta }}(0)`$ (3)
| $`p/\mathrm{GeV}`$ | $`a`$ | $`g(5^0)`$ | $`g(10^0)`$ | $`g(20^0)`$ |
| --- | --- | --- | --- | --- |
| 10 | -1.50 | 0.994 | 0.978 | 0.910 |
| 30 | -1.28 | 0.995 | 0.981 | 0.925 |
| 100 | -0.94 | 0.996 | 0.986 | 0.944 |
| 300 | -0.61 | 0.998 | 0.991 | 0.963 |
| 1000 | -0.22 | 0.999 | 0.997 | 0.987 |

Note: the entries are differential values, they have not been integrated over $`\theta `$.

Since not all experiments quote the range of accepted zenith angles, we cannot correct for this effect. We have to keep in mind that this might cause a bias, especially at low momenta.
* Atmospheric pressure/temperature profile
Previous calculations and measurements indicate that the relative muon flux variation $`\mathrm{\Delta }\mathrm{\Phi }`$ at ground level is related to the temperature-distribution in the atmosphere via

$`{\displaystyle \frac{\mathrm{\Delta }\mathrm{\Phi }}{\mathrm{\Phi }}}=\alpha {\displaystyle \frac{\mathrm{\Delta }T_{eff}}{T_{eff}}}.`$ (4)
$`\mathrm{\Phi }`$ is the integral flux above a certain muon momentum threshold $`p_{th}`$. $`T_{eff}`$ is the absolute effective temperature of the higher atmosphere. $`\alpha `$ is the temperature coefficient, which is a function of zenith angle and muon energy. For zenith angles $`\theta 0`$:

$`\alpha =\left[1+{\displaystyle \frac{70\mathrm{GeV}}{p_{th}}}\right]^1`$ (5)
Simulations using CORSIKA arrive at similar conclusions. Example: At a threshold of $`70\mathrm{GeV}`$ the formula yields $`\alpha =0.5`$. Since the atmospheric temperature, with a typical value of $`220\mathrm{K}`$, varies over the year by up to $`\pm 5\mathrm{K}`$, this implies a flux change of $`\pm 1\%`$.

For muon momenta above 10 GeV the pressure at ground level is not expected to show a significant correlation with the flux.

Unfortunately, most experiments do not report the atmospheric temperature, nor do they correct for this effect. It is even not clear, how to define the reference value. Therefore, we cannot correct for atmospheric effects.
* Unfolding of the momentum spectrum
The measured muon spectrum agrees with the true spectrum only if the momentum resolution is small compared to the momenta being investigated. Otherwise, the steepness of the spectrum, which falls off approximately according to

$`{\displaystyle \frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}p}}p^3`$ (6)
leads to an asymmetric distortion, an enhancement of the measured flux at high momenta. Thus, the measured spectrum needs to be unfolded for experimental resolution effects. In the most simple approach \- assuming the spectrum is roughly known - this can be achieved by a simple correction factor, which has been calculated in . The authors assume the spectrum (6) and a Gaussian error distribution in the variable $`1/p`$ with width $`\sigma _{1/p}`$. Often the experimental resolution is given in terms of the ‘Maximum Detectable Momentum’ $`p_{MDM}`$, defined as the momentum value for which the integral over the Gauss distribution becomes $`1/2`$:

$`E({\displaystyle \frac{1/p_{MDM}}{\sqrt{2}\sigma _{1/p}}})={\displaystyle \frac{1}{2}}\mathrm{with}E(x){\displaystyle \frac{2}{\sqrt{\pi }}}{\displaystyle _0^x}e^{t^2}𝑑t`$ (7)
Thus,

$`{\displaystyle \frac{1}{p_{MDM}}}=0.6745\sigma _{1/p}`$ (8)
The ratio of the measured and true spectra is then given by

$`R({\displaystyle \frac{p_{MDM}}{p}})=E(0.4769{\displaystyle \frac{p_{MDM}}{p}})+1.1829{\displaystyle \frac{p}{p_{MDM}}}\mathrm{exp}(0.2275{\displaystyle \frac{p_{MDM}^2}{p^2}}))`$ (9)
The measured flux must be multiplied by $`1/R`$ to correct for the experimental resolution. Figure 1 shows the dependence of R on $`1/p`$. For $`p<0.3p_{MDM}`$ the correction amounts to less than $`1\%`$ and can be neglected. For higher momenta the correction rises strongly and must be taken into account.

We have assumed that the experimenters have corrected their data for momentum resolution effects or that they can be neglected. However, several papers are not very clear on this point. Therefore, some published spectra might be biased towards too high flux values at large momenta.

3. Experimental Data

Only published results are taken into account. In appendix A we summarize the characteristics of all relevant experiments/publications, in chronological order. The spectrum and charge ratio data used in this compilation are listed explicitly in appendices B and C.

4. Absolute muon flux

There are two aspects to the measurement of the absolute muon flux, namely the shape of the spectrum as a function of energy and the absolute normalization. Some experiments only measure the relative muon flux as a function of momentum. Therefore we will analyze the data in two steps. First we check the spectral shape, leaving the normalization as a free parameter. Secondly we determine the absolute normalization of the spectrum.

4.1 The shape of the muon spectrum

A whole range of experiments are performed to measure the muon flux, the measurements used are listed in appendix B. We have corrected the datasets for altitude, which is a small correction in most cases. In order to be able to compare the datasets, we fit each set to a reference shape, using the data with momenta above $`10\mathrm{GeV}`$. In this fit, and in the following, we assume that the measurement performed in each momentum bin is independent of the other momentum bins. The reference shape is taken from the theoretical calculation by Bugaev et al , leaving the normalization as a free parameter. In general this shape provides a good description of the datasets, as can be seen below. The results of our fit are listed in table 1.

In this table, we separated the data from Rastin into two sets, as different normalizations are used in their paper. Next to fitting the normalization, we also calculated the normalization by comparing the integrated flux above $`10\mathrm{GeV}/\mathrm{c}`$ (or the lower cutoff of the experiment whichever is higher) to an integrated flux calculation using the reference shape. The data published in references ,  and  do not allow for this normalization method. As these papers are normalized to Rossi  we recalculate this normalization point with the reference shape. The data of  are normalized to the differential flux at 10 GeV/c. The results of this calculation are shown in the last column of table 1. In general, both normalizations are in good agreement.

The high energy part of the Rastin data is shown in figure 2. Here and in the following we present all spectra weighted with $`p^3`$, a common practice to compensate for the steep fall-off with momentum. Figure 2 nicely shows that the reference shape fits the data rather well, which justifies the use of the Bugaev curve as a reference. However, the data has the tendency to be slightly higher than the normalized curve at the higher momentum values.
The $`\chi ^2`$ of the fit as listed in this table made us re-check five datasets; the first four are shown in figure 3. The data of Holmes clearly show

that a simple re-normalization will not work. The data points do not follow the reference shape, especially at higher momenta. In their paper Holmes et al. apply additional corrections to the highest two data points, indicating that these are close to the MDM of the detector. Unfortunately, the value of this maximal momentum is not mentioned.
The Appleton data are scattered a lot around the curve. With the value of $`\chi ^2/NDF`$ being only slightly less than 2, this plot suggests that the errors could be underestimated.
The Allkofer data have a completely different shape. The data rise faster than the reference shape and plateau at a lower value. This plateau also seems to be wider than suggested by the reference distribution.
The Nandi data rise to a significantly higher value than predicted by the reference shape. This and the low value of the third data point create the large $`\chi ^2`$.

The Ayre data (left side of fig. 4) start off below the curve and continue to rise longer than expected. Therefore, the peak is at a higher value, but the drop-off rate seems to be similar as predicted on a log(p) scale. The curve suggests that the momentum could be over-estimated. A best fit of the momentum scale leads to a scaling of the momenta by a factor of 0.825, see right side of fig. 4. The $`\chi ^2/NDF`$ improves from 348/44 to 136/44. It naturally changes the normalization. The muon spectrum closely follows a $`p^3`$-dependence, thus the normalization is changed to about 56 % of the original. This is in fact what we observe. Even using the modified momentum, for which we cannot find a justification, the Ayre data do not fit the curve very well.

If we ignore these five datasets for the moment, we can compare the remaining data to the reference curve. We do this by applying the normalization calculated as outlined before, and listed in table 1 as ‘normalization from integration’. The result is shown in figure 5.

The top part of this plot shows a direct comparison of all the datasets to the theory, while the bottom part shows the relative difference between these measurements and the description. The larger differences at higher momenta are mainly due to the data by Rastin et al. . If we include the remaining five datasets the difference between the shape of the data and the reference shape increases, especially at higher momenta. This can be seen in figure 6.
The relative differences are plotted in 8 bins per decade, equidistant in log(p). This choice represents a good compromise taking into account experimental uncertainties and the rate of change of the $`p^3`$-weighted spectrum with momentum. The $`\chi ^2/NDF`$ values and the average flux values in the units $`\mathrm{m}^2\mathrm{sr}^1\mathrm{s}^1\mathrm{GeV}^1`$ are shown for each bin in table 2.

The five datasets discussed above have a large impact on the $`\chi ^2`$ of the relative difference in these bins. Therefore, we will exclude them when adjusting the shape according to the measurements. We fit a third degree polynomial to the logarithm of the flux as a function of the logarithm of momentum. We parameterize this function as follows:

$$\begin{array}{cc}\hfill H(y)& =H_1(y^3/25y^2/2+3y)\hfill \\ & +H_2(2y^3/3+3y^210y/3+1)\hfill \\ & +H_3(y^3/6y^2/2+y/3)\hfill \\ & +S_2(y^3/32y^2+11y/32)\hfill \\ \hfill y& ={}_{}{}^{10}\mathrm{log}(p/\mathrm{GeV})\hfill \\ \hfill F(p)& =10^{H(y)}\mathrm{m}^2\mathrm{sr}^1\mathrm{s}^1\mathrm{GeV}^1\hfill \end{array}$$
(10)
This parameterization is similar to the one used by , however the fit variables are chosen such that they have a simple interpretation: $`H_1`$, $`H_2`$, and $`H_3`$ represent the logarithm of the differential flux at 10, 100 and 1000 GeV, $`S_2`$ represents the exponent of the differential flux at 100 GeV. The $`\chi ^2/NDF`$ of this fit is 8/12, the correlation matrix is shown in appendix D. The fitted parameters are listed below, as well as the equivalent values from our reference shape.

|  | $`H_1`$ | $`H_2`$ | $`H_3`$ | $`S_2`$ |
| --- | --- | --- | --- | --- |
| Reference | 0.127 | -2.539 | -5.86 | -2.00 |
| Fit | $`0.135\pm 0.002`$ | $`2.529\pm 0.004`$ | $`5.76\pm 0.03`$ | $`2.10\pm 0.03`$ |

Figure 7 contains a $`1\sigma `$ error band, which can be approximated by

$`\delta =0.003+0.00015{\displaystyle \frac{p}{\mathrm{GeV}}}`$ (11)
It represents the size of the combined relative experimental uncertainties as a function of momentum bin.

The shape uncertainty at reference momenta of 10, 100 and 1000 GeV are 0.5%, 1.8% and 15% respectively. Above 200 GeV the uncertainty rapidly increases, indicating that more measurements are needed at these momenta.

4.2 The absolute normalization of the muon spectrum

In section 4.1 we fitted a functional shape to all renormalized datasets. The renormalization was such that the integral flux above 10 GeV corresponds to the calculation by Bugaev. We will now fit the functional shape obtained in 4.1 to the datasets of those experiments providing an absolute flux measurement, while leaving the normalization as a free parameter. Therefore we fit the function

$`F(p)`$ $`=C10^{H(y)}`$ (12)
The result is shown in table 3.

We again ignore the datasets with a very high $`\chi ^2/NDF`$ (Allkofer, Ayre). The three remaining data sets with the largest normalization factors (Bateman, Green and Tsuji) are measurements performed with solid iron magnet spectrometers, whereas the other three (Kremer 1994 and 1997 and De Pascale) use the same superconducting magnet.
We will first average the normalizations performed by the same collaboration (Bateman and Green, Kremer), and afterwards calculate the normalization measurements performed by the solid iron magnet spectrometers and the superconducting magnet spectrometers, which gives the following results:

| Solid iron magnets | $`0.937\pm 0.012`$ |
| --- | --- |
| Superconducting magnet | $`0.811\pm 0.007`$ |

These values are clearly not in agreement. We have no explanation for this observation. We will simply take the average of these two values to be our normalization and half the difference to be the uncertainty. We arrive to our final value of $`0.874\pm 0.063`$, thus a normalization with a relative uncertainty of 7 %.

4.3 The muon spectrum

In the preceding sections we have parameterized the muon spectrum at sea-level. We summarize the parameters obtained:

| $`C`$ | $`H_1`$ | $`H_2`$ | $`H_3`$ | $`S_2`$ |
| --- | --- | --- | --- | --- |
| $`0.874\pm 0.063`$ | $`0.135\pm 0.002`$ | $`2.529\pm 0.004`$ | $`5.76\pm 0.03`$ | $`2.10\pm 0.03`$ |

Our description, as well as the calculation from Bugaev, and the measurements used for the normalization are shown in figure 8. The error is given be the estimated normalization uncertainty of $`7\%`$ and the shape error in (11), added in quadrature.

If we compare our description of the differential flux to the theoretical description of Bugaev  we get the following :

| Momentum: | 10 GeV | 100 GeV | 1000 GeV |  |
| --- | --- | --- | --- | --- |
| Our Result | $`1.19\pm 0.08`$ | $`(2.59\pm 0.19)10^3`$ | $`(1.52\pm 0.26)10^6`$ | $`\mathrm{m}^2\mathrm{sr}^1\mathrm{s}^1\mathrm{GeV}^1`$ |
| Bugaev et al. | $`1.34`$ | $`2.8910^3`$ | $`1.3910^6`$ | $`\mathrm{m}^2\mathrm{sr}^1\mathrm{s}^1\mathrm{GeV}^1`$ |

At 10 GeV the measured flux is 89% of the calculation of Bugaev. However, the measured shape is slightly less steep, and at 1 TeV we arrive to a value which is close to the predicted one.

5. Charge Ratio

The charge ratio $`R_\mu `$ is defined as the ratio of vertical fluxes for positive and negative muons at sea level.

The measured charge ratios together with the published uncertainties are listed in appendix C. Figure 9 shows all values as a function of momentum.

In order to study the momentum dependence we have grouped all 15 data sets into momentum bins chosen to be equidistant in $`\mathrm{log}p`$. We have combined the different measurements by assuming that they are uncorrelated. The bin size is relatively large, since a strong momentum dependence is not expected. The result is shown in figure 10.

The two data points around $`500\mathrm{GeV}`$ and the single measurement above $`1\mathrm{TeV}`$ have huge uncertainties ($`20\%`$) and are therefore not included in the figure.

For all eight momentum bins the $`\chi ^2`$ values are good or at least acceptable; this implies the various experimental data agree among each other.
The third momentum yields the highest charge ratio. However, this result can not be attributed to a single ‘outlier’.

To see if there is a momentum dependence we have performed the following fits to the charge ratio values shown in figure 10 as a function of $`\mathrm{log}p`$:

a) $`f(\mathrm{log}p)=R_\mu ^0=\mathrm{const}`$
This gives a good fit with $`\chi ^2/NDF=158/142`$. The resulting charge ratio of

$`R_\mu ^0=1.270\pm 0.003`$ (13)
is displayed in figure 10 as horizontal line.

b) $`f(\mathrm{log}p)=R_\mu ^1+S_\mu ^1\mathrm{log}(p/\mathrm{GeV})`$
Naturally this fit is satisfactory, too. The slope comes out as

$`S_\mu ^1=0.006\pm 0.011`$ (14)
which is compatible with zero.

Therefore, the measured charge ratios are consistent with the hypothesis of being momentum independent in the range $`10\mathrm{GeV}p300\mathrm{GeV}`$.

We have looked at the data in more detail and tried to answer the following questions:
i) Do the different experiments agree with each other ?
ii) Is the ‘peak’ at about $`30\mathrm{GeV}`$ significant ?

The previous statistical analyses and figure 9 seem to imply the answer ‘yes’ to the first question. However, when separately plotting the two (by far) most precise data sets (Baxendale 1975 and Rastin 1984), one finds the discrepancy displayed in figure 11.

Averaged over all momenta the mean values

$`R_\mu ^0(Baxendale)`$ $`=`$ $`1.2799\pm 0.0042`$ (15)
$`R_\mu ^0(Rastin)`$ $`=`$ $`1.2511\pm 0.0058`$ (16)
disagree on the $`4\sigma `$ level. To reduce the discrepancy to about $`1\sigma `$ we assume - in the spirit of the Particle Data Group \- that for all experiments an additional systematic error of $`\pm 0.015`$ must be added, in form of a scale uncertainty common to all measurements of one experiment, independent of momentum. Clearly, this is a crude model!

Including these errors results in the charge ratios displayed in figure 12. Note that the values are quite close to those in figure 10, while the error bars are enlarged. The corresponding numbers are listed in table 4. The momenta are calculated from the central values of the logarithmic bins.

There is no final answer to question ii). While the measurements by Baxendale (and also others, with larger errors) indicate an increase of the charge ratio at momenta around $`30\mathrm{GeV}`$, the data by Rastin do not support this hypothesis. For the moment the measurements are consistent with the simple hypothesis of a momentum independent charge ratio.

We try to summarize the charge ratio measurements and their uncertainty ($`68\%`$ CL) with the following formula:

$`R_\mu =1.268\pm \left[0.008+0.0002{\displaystyle \frac{p}{\mathrm{GeV}}}\right]`$ (17)
in the momentum range $`10300\mathrm{GeV}`$. Figure 12 shows the corresponding mean value and the error band. The central value is the mean of all measurements, taking into account the additional systematic error of $`0.015`$. The momentum dependent error is estimated such that it is roughly of the same size as the uncertainties of the corresponding data points in figure 12 and table 4.

Theoretical models of atmospheric showers must be able to reproduce these data, the calculated charge ratios should fall into the band given in equation (17) and figure 12. Clearly, at high momenta precise data is still lacking. Momenta above a few hundred GeV are of particular interest, since a growing influence of kaons and a resulting increase of the charge ratio is predicted.

6. Summary and conclusions

We have combined the published data on the vertical muon spectrum and charge ratio at sea-level. In this comparison we have found that the differential spectrum can be described using a simple formula. The shape of the momentum spectrum is well measured at momenta below 100 GeV. Above 200 GeV only a few data points exist, therefore the uncertainty increases to 17 % at 1 TeV.

Several experiments measure the absolute normalization of the spectrum. Our combined result is compared to the calculation by Bugaev et al. At 10 GeV the measured flux is 11 % below the calculated one.

The charge ratio is reported by many experiments. The combined result favors a momentum independent value of the charge ratio of 1.268. The error on the charge ratio increases rapidly above 200 GeV due to a lack of precise experimental data in that region.

7. Acknowledgments

We thank Dieter Heck for helping us to run CORSIKA and Albert van Mil for collecting charge ratio data from the literature.

References

APPENDIX A Measurements.

The following three tables list the experiments/publications we considered, in chronological order, together with the most important parameters.

Some entries are missing since the corresponding figures are not published. In particular the zenith angle regime is characterized frequently only by verbal expressions like ‘near vertical’. We distinguish three magnet types: solid iron and conventional coil, air gap magnet with conventional coil and air gap magnet with superconducting coil. The period of data taking is indicated by year and month, e.g. 59/11 stands for November 1959.

The following remarks refer to the experiment numbers in column 1 of the tables. The figures mentioned are those of the corresponding publication.

1) Not used, since data set is a subset of 3).

2) Spectrum data read off from figure 1; charge ratio taken from 3).

6) Spectrum normalized to Rossi; data read off from figures 8 and 15.

7) Spectrum normalized to Rossi.

8) Spectrum normalized to Rossi; data read off from figures 1 and 3.

9) The two values at $`p=240\mathrm{GeV}`$ are not statistically independent. We have calculated the arithmetic mean of the two figures and took the smaller of the two statistical errors as the uncertainty.

10) We do not use these data, which were obtained in the equator region, where the geomagnetic cutoff is large (14 GeV) and may influence the muon flux.

11) The spectrum data form a subset of those published in 13). Nevertheless we consider it separately, since the normalization procedures are slightly different. When calculating world averages we exclude these data. We do not use the charge ratio data, since they are included in the superset published in 13).

13) Spectrum: no absolute flux determination, only normalization to previous measurements by other experiments. Charge ratio: A few measurements are given with slightly asymmetric errors; they have been ‘symmetrized’ by shifting the central value to the center of the error interval.

15) Apparatus under concrete shelter of $`868\mathrm{g}/\mathrm{cm}^2`$.

16) We use only the charge ratio value obtained at the town of Kiel; for the other measurements, made at muon momenta close to and below 10 GeV in the equator region, the geomagnetic cutoff is large (14 GeV) and may influence the muon flux.

20) Resolution correction is based upon MDM = 100 GeV; if MDM of 350 GeV is used, spectrum is flatter and better consistent with 18) . Points up to 100 GeV are considered reliable.

21) Spectrum data are normalized to an integral intensity at 5 GeV. The spectrum provided is the theoretical spectrum that fits the data best.

22) Above 50 GeV authors question results due to resolution.

24, 25) Same magnet as 22). Charge ratio and its error read off from figure 2.

| No | author(s), | name | location | coordinates | alt. | magnet | detector(s) | period | MDM | zenith | spec. | ratio |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
|  | reference |  |  |  | /m |  |  |  | /GeV |  |  |  |
| 1 | B.G. Owen and |  | Manchester | 53<sup>0</sup>N, 2<sup>0</sup>W | 50 | air | Geiger, |  | 30 | $`<8^0`$ | no | yes |
|  | J.G. Wilson, 1949 |  | (Great Britain) |  |  |  | flash tubes |  |  |  |  |  |
| 2 | D.E. Caro |  | Melbourne | 38<sup>0</sup>S, 145<sup>0</sup>E | 50 | air | counters |  | $``$ 50 |  | yes | yes |
|  | et al, 1950 |  | (Australia) |  |  |  |  |  |  |  |  |  |
| 3 | B.G. Owen and |  | Manchester | 53<sup>0</sup>N, 2<sup>0</sup>W | 50 | air | Geiger, |  | 30 | $`<8^0`$ | no | yes |
|  | J.G. Wilson, 1951 |  | (Great Britain) |  |  |  | flash tubes |  |  |  |  |  |
| 4 | I. Filosofo |  | Agordo | 46<sup>0</sup>N, 12<sup>0</sup>E | 600 | iron | counters |  |  |  | no | yes |
|  | et. al, 1954 |  | (Italy) |  |  |  |  |  |  |  |  |  |
| 5 | B.G. Owen and |  | Manchester | 53<sup>0</sup>N, 2<sup>0</sup>W | 50 | air | counters |  |  | $`<10^0`$ | yes | no |
|  | J.G. Wilson, 1955 |  | (Great Britain) |  |  |  |  |  |  |  |  |  |
| 6 | J. Pine et al, |  | Cornell | 42<sup>0</sup>N, 76<sup>0</sup>W | 500 | air | Geiger, |  | 175 |  | yes | yes |
|  | 1959 |  | (USA) |  |  |  | cloud ch. |  |  |  |  |  |
| 7 | J.E.R. Holmes |  | Manchester | 53<sup>0</sup>N, 2<sup>0</sup>W | 50 | air | cloud ch., | 53-55 |  | $`<10^0`$ | yes | yes |
|  | et al, 1961 |  | (Great Britain) |  |  |  | Geiger |  |  |  |  |  |
| 8 | W. Pak et al, |  | Cornell | 42<sup>0</sup>N, 76<sup>0</sup>W | 500 | air | Geiger, |  | 175 |  | yes | yes |
|  | 1961 |  | (USA) |  |  |  | hodosc. |  |  |  |  |  |
| 9 | P.J. Hayman and A.W. |  | Durham | 54<sup>0</sup>N, 1<sup>0</sup>W | 70 | air | Geiger, | 59/11-60/03 | 657 |  | yes | yes |
|  | Wolfendale, 1962 |  | (Great Britain) |  |  |  | flash tubes |  |  |  |  |  |
| 10 | O.C. Allkofer |  | near | 29<sup>0</sup>W, 0<sup>0</sup>N | 0 | air | spark ch. |  |  |  | yes | yes |
|  | et al, 1968 |  | equator | 22<sup>0</sup>W, 1<sup>0</sup>S |  |  |  |  |  |  |  |  |

| No | author(s), | name | location | coordinates | alt. | magnet | detector(s) | period | MDM | zenith | spec. | ratio |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
|  | reference |  |  |  | /m |  |  |  | /GeV |  |  |  |
| 11 | S.R. Baber |  | Nottingham | 53<sup>0</sup>N, 1<sup>0</sup>W | 52 | iron | Geiger | 64/06-65/05 | 360 |  | yes | yes |
|  | et al, 1968 |  | (Great Britain) |  |  |  |  |  |  |  |  |  |
| 12 | A.M. Aurela and A.W. |  | Durham | 54<sup>0</sup>N, 1<sup>0</sup>W | 70 | air, | Geiger, | 64/06-65/01 |  |  | yes | no |
|  | Wolfendale, 1967 |  | (Great Britain) |  |  | iron | flash |  |  |  |  |  |
| 13 | I.C. Appleton |  | Nottingham | 53<sup>0</sup>N, 1<sup>0</sup>W | 52 | iron | flash | 64-68 | 360 |  | yes | yes |
|  | et al, 1971 |  | (Great Britain) |  |  |  | tubes |  |  |  |  |  |
| 14 | B.J. Bateman | AMH | College Station | 31<sup>0</sup>N, 96<sup>0</sup>W | 80 | iron | spark ch., |  |  |  | yes | no |
|  | et al, 1971 |  | (USA) |  |  |  | scint. |  |  |  |  |  |
| 15 | O.C. Allkofer |  | Kiel | 54<sup>0</sup>N, 11<sup>0</sup>E | 10 | iron | spark ch., |  | $`>1000`$ |  | yes | no |
|  | et al, 1971 |  | (Germany) |  |  |  | scint. |  |  |  |  |  |
| 16 | O.C. Allkofer and |  | Kiel | 54<sup>0</sup>N, 11<sup>0</sup>E | 10 | iron | spark ch., |  |  |  | no | yes |
|  | W.D. Dau, 1972 |  | (and equator) |  |  |  | scint. |  |  |  |  |  |
| 17 | B.C. Nandi |  | Durgapur | 24<sup>0</sup>N, 87<sup>0</sup>E | 70 | iron | flash tubes, | 69/02-70/02 | $`985`$ |  | yes | yes |
|  | et al, 1972 |  | (India) |  |  |  | Geiger |  |  |  |  |  |
| 18 | C.A. Ayre | MARS | Durham | 54<sup>0</sup>N, 1<sup>0</sup>W | 70 | iron | scint., | 72/05-73/01 | 670 |  | yes | no |
|  | et al, 1975 |  | (Great Britain) |  |  |  | flash tubes |  |  |  |  |  |
| 19 | J.M. Baxendale |  | Durham | 54<sup>0</sup>N, 1<sup>0</sup>W | 70 | iron | scint., | 72/02-72/12 |  |  | yes | yes |
|  | et al, 1975 |  | (Great Britain) |  |  |  | flash tubes |  |  |  |  |  |
| 20 | P.J. Green | AMH | Houston | 30<sup>0</sup>N, 95<sup>0</sup>W | 10 | iron | spark ch., |  | 345 | $`<9^0`$ | yes | no |
|  | et al, 1979 |  | (USA) |  |  |  | scint. |  |  |  |  |  |

| No | author(s), | name | location | coordinates | alt. | magnet | detector(s) | period | MDM | zenith | spec. | ratio |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
|  | reference |  |  |  | /m |  |  |  | /GeV |  |  |  |
| 21 | B.C. Rastin |  | Nottingham | 53<sup>0</sup>N, 1<sup>0</sup>W | 52 | iron | flash tubes, | 74/09-78/05 | 3400 |  | yes | yes |
|  | 1984 |  | (Great Britain) |  |  |  | scint. |  |  |  |  |  |
| 22 | M.P. De Pascale | MASS | Prince Albert | 53<sup>0</sup>N, 106<sup>0</sup>W | 600 | air | MWPC, | 89/08 | 118 |  | yes | yes |
|  | et al, 1993 |  | (Canada) |  |  | (superc.) | scint., TOF |  |  |  |  |  |
| 23 | S. Tsuji |  | Okayama | 34<sup>0</sup>N, 134<sup>0</sup>E | 5 | iron | drift, | 92/09-97/12 | 270 | 0<sup>0</sup>-1<sup>0</sup> | yes | no |
|  | et al, 1998 |  | (Japan) |  |  |  | scint. |  |  |  |  |  |
| 24 | J. Kremer | CAPRICE | Lynn Lake | 57<sup>0</sup>N, 101<sup>0</sup>W | 360 | air | prop., drift, | 94/07 | 175 | 0<sup>0</sup>-20<sup>0</sup> | yes | yes |
|  | et al, 1999 |  | (Canada) |  |  | (superc.) | TOF, RICH |  |  |  |  |  |
| 25 | J. Kremer | CAPRICE | Fort Sumner | 34<sup>0</sup>N, 104<sup>0</sup>W | 1270 | air | drift, | 97/04-97/05 | 175 | 0<sup>0</sup>-20<sup>0</sup> | yes | yes |
|  | et al, 1999 |  | (USA) |  |  | (superc.) | TOF, RICH |  |  |  |  |  |

APPENDIX B Muon flux data.

The following lists contain all spectrum data for momenta above $`10\mathrm{GeV}`$. The flux is given in ($`\mathrm{m}^2\mathrm{sr}^1\mathrm{s}^1\mathrm{GeV}`$). Each line contains the momentum in $`\mathrm{GeV}`$ together with the published value and uncertainty of the flux. Altitude corrections have been applied later and are not included in the figures listed here.

* D.E. Caro et al, 1950

| 12.2 | (6.46 $`\pm `$ 0.75)$`10^1`$ |
| --- | --- |
| 16.9 | (2.59 $`\pm `$ 0.26)$`10^1`$ |
| 17.8 | (2.38 $`\pm `$ 0.25)$`10^1`$ |
| 26.5 | (7.4 $`\pm `$ 0.9)$`10^2`$ |
| 51 | (1.3 $`\pm `$ 0.2)$`10^2`$ |
* B.G. Owen and J.G. Wilson, 1955

| 10.0 | 1.09 $`\pm `$ 0.03 |
| --- | --- |
| 15.0 | (4.36 $`\pm `$ 0.11)$`10^1`$ |
| 20.0 | (2.20 $`\pm `$ 0.06)$`10^1`$ |
* J.Pine et al, 1959

| 11.6 | (8.26 $`\pm `$ 0.99)$`10^1`$ |
| --- | --- |
| 13.7 | (5.24 $`\pm `$ 0.84)$`10^1`$ |
| 16.3 | (3.05 $`\pm `$ 0.49)$`10^1`$ |
| 18.1 | (2.79 $`\pm `$ 0.33)$`10^1`$ |
| 19.8 | (1.91 $`\pm `$ 0.31)$`10^1`$ |
| 25.5 | (1.16 $`\pm `$ 0.19)$`10^1`$ |
| 29.8 | (7.68 $`\pm `$ 0.92)$`10^2`$ |
| 35.4 | (4.95 $`\pm `$ 0.79)$`10^2`$ |
| 46.0 | (1.99 $`\pm `$ 0.32)$`10^2`$ |
| 59.1 | (9.29 $`\pm `$ 1.49)$`10^3`$ |
| 83.2 | (4.60 $`\pm `$ 0.87)$`10^3`$ |
| 125. | (1.34 $`\pm `$ 0.33)$`10^3`$ |
* W.Pak et al, 1961

| 12.9 | (6.14 $`\pm `$ 0.68)$`10^1`$ |
| --- | --- |
| 15.8 | (3.35 $`\pm `$ 0.37)$`10^1`$ |
| 19.8 | (1.89 $`\pm `$ 0.21)$`10^1`$ |
| 24.1 | (1.35 $`\pm `$ 0.19)$`10^1`$ |
| 26.4 | (1.03 $`\pm `$ 0.11)$`10^1`$ |
| 37.4 | (3.77 $`\pm `$ 0.41)$`10^2`$ |
| 58.9 | (1.29 $`\pm `$ 0.14)$`10^2`$ |
* J.E.R. Holmes et al, 1961

| 11 | (7.88 $`\pm `$ 0.32)$`10^1`$ |
| --- | --- |
| 13 | (5.90 $`\pm `$ 0.30)$`10^1`$ |
| 16 | (3.68 $`\pm `$ 0.18)$`10^1`$ |
| 19 | (2.33 $`\pm `$ 0.12)$`10^1`$ |
| 23 | (1.67 $`\pm `$ 0.12)$`10^1`$ |
| 28 | (1.04 $`\pm `$ 0.07)$`10^1`$ |
| 36 | (5.75 $`\pm `$ 0.40)$`10^2`$ |
| 49 | (2.02 $`\pm `$ 0.20)$`10^2`$ |
| 67 | (8.65 $`\pm `$ 1.30)$`10^3`$ |
| 89 | (4.38 $`\pm `$ 0.70)$`10^3`$ |
| 134 | (1.42 $`\pm `$ 0.26)$`10^3`$ |
| 271 | (2.8 $`\pm `$ 0.6)$`10^4`$ |
| 1160 | (4.8 $`\pm `$ 2.3)$`10^6`$ |
* P.J. Hayman and A.W. Wolfendale 1962

| 10.8 | (8.51 $`\pm `$ 0.26)$`10^1`$ |
| --- | --- |
| 12.4 | (6.14 $`\pm `$ 0.17)$`10^1`$ |
| 14.6 | (4.35 $`\pm `$ 0.11)$`10^1`$ |
| 17.8 | (2.52 $`\pm `$ 0.07)$`10^1`$ |
| 22.6 | (1.39 $`\pm `$ 0.04)$`10^1`$ |
| 31.3 | (5.85 $`\pm `$ 0.15)$`10^2`$ |
| 42.3 | (2.88 $`\pm `$ 0.11)$`10^2`$ |
| 56.1 | (1.22 $`\pm `$ 0.05)$`10^2`$ |
| 72.5 | (5.75 $`\pm `$ 0.46)$`10^3`$ |
| 88.1 | (3.27 $`\pm `$ 0.29)$`10^3`$ |
| 112 | (1.36 $`\pm `$ 0.14)$`10^3`$ |
| 153 | (5.18 $`\pm `$ 0.93)$`10^4`$ |
| 244 | (1.14 $`\pm `$ 0.22)$`10^4`$ |
| 413 | (1.98 $`\pm `$ 0.99)$`10^5`$ |
| 894 | (1.84 $`\pm `$ 1.01)$`10^6`$ |
* A.M. Aurela and A.W. Wolfendale 1967

| 15.1 | (4.25 $`\pm `$ 0.16)$`10^1`$ |
| --- | --- |
| 41.5 | (3.40 $`\pm `$ 0.44)$`10^2`$ |
| 82.1 | (4.10 $`\pm `$ 0.35)$`10^3`$ |
* S.R. Baber et al, 1968

| 11.60 | (7.77 $`\pm `$ 0.26)$`10^1`$ |
| --- | --- |
| 15.22 | (4.22 $`\pm `$ 0.21)$`10^1`$ |
| 19.20 | (2.42 $`\pm `$ 0.12)$`10^1`$ |
| 24.00 | (1.39 $`\pm `$ 0.07)$`10^1`$ |
| 33.5 | (5.78 $`\pm `$ 0.35)$`10^2`$ |
| 50.0 | (1.90 $`\pm `$ 0.16)$`10^2`$ |
| 81.0 | (4.59 $`\pm `$ 0.60)$`10^3`$ |
| 127.0 | (1.14 $`\pm `$ 0.18)$`10^3`$ |
| 266.0 | (1.00 $`\pm `$ 0.24)$`10^4`$ |
| 810.0 | (2.11 $`\pm `$ 0.55)$`10^6`$ |
* I.C. Appleton et al, 1971

| 12.84 | (2.92 $`\pm `$ 0.04)$`10^1`$ |
| --- | --- |
| 17.2 | (1.46 $`\pm `$ 0.02)$`10^1`$ |
| 24.3 | (5.78 $`\pm `$ 0.13)$`10^2`$ |
| 33 | (2.54 $`\pm `$ 0.06)$`10^2`$ |
| 43.4 | (1.09 $`\pm `$ 0.10)$`10^2`$ |
| 45.6 | (1.07 $`\pm `$ 0.09)$`10^2`$ |
| 48.2 | (9.00 $`\pm `$ 0.81)$`10^3`$ |
| 51.0 | (6.95 $`\pm `$ 0.67)$`10^3`$ |
| 54.2 | (6.84 $`\pm `$ 0.63)$`10^3`$ |
| 57.9 | (5.20 $`\pm `$ 0.52)$`10^3`$ |
| 62.0 | (4.79 $`\pm `$ 0.47)$`10^3`$ |
| 66.9 | (2.80 $`\pm `$ 0.33)$`10^3`$ |
| 72.6 | (2.76 $`\pm `$ 0.31)$`10^3`$ |
| 79.3 | (2.02 $`\pm `$ 0.23)$`10^3`$ |
| 87.5 | (1.57 $`\pm `$ 0.17)$`10^3`$ |
| 97.6 | (1.17 $`\pm `$ 0.13)$`10^3`$ |
| 110.3 | (6.97 $`\pm `$ 0.94)$`10^4`$ |
| 129.9 | (3.00 $`\pm `$ 0.54)$`10^4`$ |
| 149.4 | (2.93 $`\pm `$ 0.29)$`10^4`$ |
| 181.5 | (1.34 $`\pm `$ 0.24)$`10^4`$ |
| 230.9 | (6.61 $`\pm `$ 1.41)$`10^5`$ |
| 316.1 | (2.72 $`\pm `$ 0.59)$`10^5`$ |
| 491.5 | (7.67 $`\pm `$ 1.69)$`10^6`$ |
| 1000.0 | (5.21 $`\pm `$ 1.34)$`10^7`$ |
* B.J. Bateman et al, 1971

| 10.0 | 1.12 $`\pm `$ 0.03 |
| --- | --- |
| 13.0 | (6.63 $`\pm `$ 0.13)$`10^1`$ |
| 15.0 | (4.56 $`\pm `$ 0.09)$`10^1`$ |
| 20.0 | (2.31 $`\pm `$ 0.05)$`10^1`$ |
| 27.0 | (1.08 $`\pm `$ 0.03)$`10^1`$ |
| 35.0 | (5.45 $`\pm `$ 0.16)$`10^2`$ |
| 40.0 | (3.78 $`\pm `$ 0.11)$`10^2`$ |
| 46.0 | (2.55 $`\pm `$ 0.13)$`10^2`$ |
| 53.0 | (1.70 $`\pm `$ 0.08)$`10^2`$ |
* O.C. Allkofer et al, 1971

| 11.4 | 1.13 $`\pm `$ 0.01 |
| --- | --- |
| 14.8 | (6.04 $`\pm `$ 0.08)$`10^1`$ |
| 20.5 | (2.51 $`\pm `$ 0.03)$`10^1`$ |
| 31.4 | (8.01 $`\pm `$ 0.13)$`10^2`$ |
| 52.3 | (1.89 $`\pm `$ 0.05)$`10^2`$ |
| 93.0 | (3.38 $`\pm `$ 0.14)$`10^3`$ |
| 175.0 | (5.19 $`\pm `$ 0.37)$`10^4`$ |
| 329.0 | (7.84 $`\pm `$ 1.12)$`10^5`$ |
| 642.0 | (6.40 $`\pm `$ 1.92)$`10^6`$ |
* B.C. Nandi and M.S. Sinha, 1972

| 11.8 | (9.43 $`\pm `$ 0.15)$`10^1`$ |
| --- | --- |
| 14.0 | (6.38 $`\pm `$ 0.14)$`10^1`$ |
| 16.4 | (4.21 $`\pm `$ 0.09)$`10^1`$ |
| 19.7 | (2.72 $`\pm `$ 0.06)$`10^1`$ |
| 24.2 | (1.41 $`\pm `$ 0.04)$`10^1`$ |
| 29.6 | (1.01 $`\pm `$ 0.03)$`10^1`$ |
| 37.1 | (5.40 $`\pm `$ 0.17)$`10^2`$ |
| 46.9 | (2.99 $`\pm `$ 0.13)$`10^2`$ |
| 60.0 | (1.45 $`\pm `$ 0.07)$`10^2`$ |
| 84.0 | (5.58 $`\pm `$ 0.28)$`10^3`$ |
| 118 | (2.04 $`\pm `$ 0.17)$`10^3`$ |
| 167 | (6.09 $`\pm `$ 0.61)$`10^4`$ |
| 260 | (1.96 $`\pm `$ 0.25)$`10^4`$ |
| 467 | (2.69 $`\pm `$ 0.67)$`10^5`$ |
| 1109 | (1.03 $`\pm `$ 0.36)$`10^6`$ |
* C.A. Ayre et al, 1975

| 21.3 | (2.096 $`\pm `$ 0.029)$`10^1`$ |
| --- | --- |
| 22.1 | (1.909 $`\pm `$ 0.027)$`10^1`$ |
| 23.1 | (1.708 $`\pm `$ 0.024)$`10^1`$ |
| 24.1 | (1.574 $`\pm `$ 0.022)$`10^1`$ |
| 25.1 | (1.432 $`\pm `$ 0.020)$`10^1`$ |
| 26.3 | (1.224 $`\pm `$ 0.017)$`10^1`$ |
| 27.7 | (1.067 $`\pm `$ 0.015)$`10^1`$ |
| 29.3 | (9.130 $`\pm `$ 0.128)$`10^2`$ |
| 31.0 | (7.968 $`\pm `$ 0.112)$`10^2`$ |
| 33.1 | (6.947 $`\pm `$ 0.097)$`10^2`$ |
| 35.3 | (5.704 $`\pm `$ 0.080)$`10^2`$ |
| 38.3 | (4.547 $`\pm `$ 0.068)$`10^2`$ |
| 40.8 | (4.208 $`\pm `$ 0.046)$`10^2`$ |
| 41.7 | (3.663 $`\pm `$ 0.055)$`10^2`$ |
| 42.8 | (3.420 $`\pm `$ 0.041)$`10^2`$ |
| 44.8 | (2.962 $`\pm `$ 0.036)$`10^2`$ |
| 45.8 | (2.797 $`\pm `$ 0.042)$`10^2`$ |
| 47.1 | (2.628 $`\pm `$ 0.032)$`10^2`$ |
| 49.3 | (2.217 $`\pm `$ 0.027)$`10^2`$ |
| 50.7 | (2.086 $`\pm `$ 0.033)$`10^2`$ |
| 52.1 | (2.014 $`\pm `$ 0.024)$`10^2`$ |
| 55.2 | (1.646 $`\pm `$ 0.021)$`10^2`$ |
| 57.0 | (1.525 $`\pm `$ 0.024)$`10^2`$ |
| 58.9 | (1.434 $`\pm `$ 0.019)$`10^2`$ |
| 63.0 | (1.123 $`\pm `$ 0.015)$`10^2`$ |
| 65.3 | (1.023 $`\pm `$ 0.017)$`10^2`$ |
| 67.9 | (9.216 $`\pm `$ 0.129)$`10^3`$ |
| 73.7 | (7.084 $`\pm `$ 0.099)$`10^3`$ |
| 76.6 | (6.585 $`\pm `$ 0.118)$`10^3`$ |
| 80.0 | (5.753 $`\pm `$ 0.081)$`10^3`$ |
| 88.3 | (4.149 $`\pm `$ 0.062)$`10^3`$ |
| 93.0 | (3.616 $`\pm `$ 0.072)$`10^3`$ |
| 98.3 | (3.252 $`\pm `$ 0.052)$`10^3`$ |
| 112.0 | (2.037 $`\pm `$ 0.035)$`10^3`$ |
| 118.0 | (1.842 $`\pm `$ 0.042)$`10^3`$ |
| 128.0 | (1.454 $`\pm `$ 0.026)$`10^3`$ |
| 145.0 | (9.603 $`\pm `$ 0.192)$`10^4`$ |
| 160.0 | (7.459 $`\pm `$ 0.201)$`10^4`$ |
| 177.0 | (5.352 $`\pm `$ 0.123)$`10^4`$ |
| 214.0 | (2.893 $`\pm `$ 0.078)$`10^4`$ |
| 236.0 | (2.171 $`\pm `$ 0.083)$`10^4`$ |
| 274.0 | (1.260 $`\pm `$ 0.044)$`10^4`$ |
| 358.0 | (5.328 $`\pm `$ 0.250)$`10^5`$ |
| 367.0 | (4.843 $`\pm `$ 0.349)$`10^5`$ |
| 442.0 | (2.764 $`\pm `$ 0.221)$`10^5`$ |
* P.J. Green et al, 1979

| 12.18 | (8.33 $`\pm `$ 0.30)$`10^1`$ |
| --- | --- |
| 19.20 | (2.96 $`\pm `$ 0.12)$`10^1`$ |
| 31.40 | (8.14 $`\pm `$ 0.50)$`10^2`$ |
| 52.40 | (1.77 $`\pm `$ 0.16)$`10^2`$ |
| 87.10 | (4.79 $`\pm `$ 0.78)$`10^3`$ |
| ( 249.90 | (3.95 $`\pm `$ 0.53)$`10^4`$ ) |
* B.C. Rastin, 1984

| 10.69 | 1.156 $`\pm `$ 0.008 |
| --- | --- |
| 11.94 | (9.05 $`\pm `$ 0.06)$`10^1`$ |
| 13.58 | (6.72 $`\pm `$ 0.04)$`10^1`$ |
| 15.81 | (4.70 $`\pm `$ 0.03)$`10^1`$ |
| 19.05 | (2.97 $`\pm `$ 0.02)$`10^1`$ |
| 24.14 | (1.63 $`\pm `$ 0.01)$`10^1`$ |
| 28.35 | (1.03 $`\pm `$ 0.02)$`10^1`$ |
| 29.30 | (9.3 $`\pm `$ 0.2)$`10^2`$ |
| 30.32 | (8.5 $`\pm `$ 0.2)$`10^2`$ |
| 31.42 | (7.8 $`\pm `$ 0.2)$`10^2`$ |
| 32.60 | (7.1 $`\pm `$ 0.2)$`10^2`$ |
| 33.88 | (6.3 $`\pm `$ 0.1)$`10^2`$ |
| 35.27 | (5.8 $`\pm `$ 0.1)$`10^2`$ |
| 36.79 | (5.1 $`\pm `$ 0.1)$`10^2`$ |
| 38.44 | (4.4 $`\pm `$ 0.1)$`10^2`$ |
| 40.25 | (3.91 $`\pm `$ 0.09)$`10^2`$ |
| 42.25 | (3.44 $`\pm `$ 0.08)$`10^2`$ |
| 44.47 | (3.03 $`\pm `$ 0.07)$`10^2`$ |
| 46.94 | (2.62 $`\pm `$ 0.06)$`10^2`$ |
| 49.71 | (2.23 $`\pm `$ 0.06)$`10^2`$ |
| 52.84 | (1.87 $`\pm `$ 0.05)$`10^2`$ |
| 56.40 | (1.56 $`\pm `$ 0.04)$`10^2`$ |
| 60.49 | (1.25 $`\pm `$ 0.03)$`10^2`$ |
| 65.23 | (1.04 $`\pm `$ 0.03)$`10^2`$ |
| 70.80 | (7.6 $`\pm `$ 0.2)$`10^3`$ |
| 77.42 | (6.3 $`\pm `$ 0.2)$`10^3`$ |
| 85.43 | (4.4 $`\pm `$ 0.1)$`10^3`$ |
| 95.34 | (3.2 $`\pm `$ 0.1)$`10^3`$ |
| 107.88 | (2.12 $`\pm `$ 0.08)$`10^3`$ |
| 124.27 | (1.42 $`\pm `$ 0.05)$`10^3`$ |
| 146.62 | (8.8 $`\pm `$ 0.4)$`10^4`$ |
| 178.85 | (4.8 $`\pm `$ 0.2)$`10^4`$ |
| 229.36 | (2.2 $`\pm `$ 0.1)$`10^4`$ |
| 319.72 | (7.5 $`\pm `$ 0.5)$`10^5`$ |
| 525.82 | (1.4 $`\pm `$ 0.1)$`10^5`$ |
| 1288.74 | (5.9 $`\pm `$ 0.8)$`10^7`$ |
* M.P. De Pascale et al, 1993

Positive Muons:

| 10.19 | (5.983 $`\pm `$ 0.233)$`10^1`$ |
| --- | --- |
| 14.42 | (2.523 $`\pm `$ 0.144)$`10^1`$ |
| 20.36 | (1.246 $`\pm `$ 0.071)$`10^1`$ |
| 28.80 | (4.709 $`\pm `$ 0.414)$`10^2`$ |
| 40.64 | (1.430 $`\pm `$ 0.162)$`10^2`$ |
| 70.16 | (5.176 $`\pm `$ 0.554)$`10^3`$ |

Negative Muons:

| 10.19 | (5.025 $`\pm `$ 0.216)$`10^1`$ |
| --- | --- |
| 14.42 | (2.124 $`\pm `$ 0.132)$`10^1`$ |
| 20.36 | (8.653 $`\pm `$ 0.580)$`10^2`$ |
| 28.80 | (3.788 $`\pm `$ 0.375)$`10^2`$ |
| 40.64 | (1.389 $`\pm `$ 0.158)$`10^2`$ |
| 70.16 | (2.423 $`\pm `$ 0.383)$`10^3`$ |
* S. Tsuji et al, 1998

| 12.1 | (8.37 $`\pm `$ 0.17)$`10^1`$ |
| --- | --- |
| 17.2 | (3.75 $`\pm `$ 0.12)$`10^1`$ |
| 22.3 | (2.04 $`\pm `$ 0.09)$`10^1`$ |
| 27.3 | (1.17 $`\pm `$ 0.07)$`10^1`$ |
| 34.3 | (6.12 $`\pm `$ 0.37)$`10^2`$ |
| 44.5 | (3.21 $`\pm `$ 0.28)$`10^2`$ |
| 54.6 | (1.39 $`\pm `$ 0.20)$`10^2`$ |
| 64.6 | (8.68 $`\pm `$ 1.64)$`10^3`$ |
| 74.7 | (8.07 $`\pm `$ 1.68)$`10^3`$ |
| 84.7 | (3.73 $`\pm `$ 1.18)$`10^3`$ |
| 94.7 | (1.69 $`\pm `$ 0.85)$`10^3`$ |
| 119.8 | (1.00 $`\pm `$ 0.32)$`10^3`$ |
| 171.2 | (2.84 $`\pm `$ 2.01)$`10^4`$ |
| 222.0 | (3.59 $`\pm `$ 2.54)$`10^4`$ |
* J. Kremer et al, 1999 1994 data

Positive Muons:

| 12.42 | (3.89 $`\pm `$ 0.08)$`10^1`$ |
| --- | --- |
| 18.85 | (1.38 $`\pm `$ 0.04)$`10^1`$ |
| 26.68 | (6.3 $`\pm `$ 0.3) $`10^2`$ |
| 36.69 | (2.8 $`\pm `$ 0.1) $`10^2`$ |
| 51.47 | (9.9 $`\pm `$ 0.7) $`10^3`$ |
| 72.08 | (3.6 $`\pm `$ 0.3) $`10^3`$ |
| 100.96 | (1.4 $`\pm `$ 0.2) $`10^3`$ |

Negative Muons:

| 12.42 | (3.09 $`\pm `$ 0.07)$`10^1`$ |
| --- | --- |
| 18.85 | (1.08 $`\pm `$ 0.03)$`10^1`$ |
| 26.68 | (4.6 $`\pm `$ 0.2) $`10^2`$ |
| 36.69 | (1.9 $`\pm `$ 0.1) $`10^2`$ |
| 51.47 | (7.1 $`\pm `$ 0.6) $`10^3`$ |
| 72.08 | (3.0 $`\pm `$ 0.3) $`10^3`$ |
| 100.96 | (1.2 $`\pm `$ 0.2) $`10^3`$ |
* J. Kremer et al, 1999 1997 data

Positive Muons:

| 12.42 | (4.14 $`\pm `$ 0.09)$`10^1`$ |
| --- | --- |
| 18.85 | (1.54 $`\pm `$ 0.04)$`10^1`$ |
| 26.68 | (6.4 $`\pm `$ 0.2) $`10^2`$ |
| 36.69 | (2.8 $`\pm `$ 0.1) $`10^2`$ |
| 51.47 | (10.2 $`\pm `$ 0.5) $`10^3`$ |
| 72.08 | (4.2 $`\pm `$ 0.3) $`10^3`$ |
| 100.96 | (1.5 $`\pm `$ 0.1) $`10^3`$ |

Negative Muons:

| 12.42 | (3.20 $`\pm `$ 0.07)$`10^1`$ |
| --- | --- |
| 18.85 | (1.16 $`\pm `$ 0.03)$`10^1`$ |
| 26.68 | (4.5 $`\pm `$ 0.2) $`10^2`$ |
| 36.69 | (2.03 $`\pm `$ 0.08)$`10^2`$ |
| 51.47 | (7.7 $`\pm `$ 0.4) $`10^3`$ |
| 72.08 | (3.2 $`\pm `$ 0.2) $`10^3`$ |
| 100.96 | (1.1 $`\pm `$ 0.1) $`10^3`$ |

APPENDIX C Charge ratio data.

The following lists contain all charge ratio data for momenta above $`10\mathrm{GeV}`$. Each line contains the momentum in $`\mathrm{GeV}`$ together with the published value and uncertainty of the charge ratio.

* D.E. Caro et al, 1950

| 35.0 | 1.6 $`\pm `$ 0.2 |
| --- | --- |
* B.G. Owen et al, 1951

| 11.5 | 1.229 $`\pm `$ 0.036 |
| --- | --- |
* I. Filosofo et al, 1954

| 21.0 | 1.232 $`\pm `$ 0.016 |
| --- | --- |
* J. Pine et al, 1959

| 19.6 | 1.303 $`\pm `$ 0.031 |
| --- | --- |
| 22.8 | 1.29 $`\pm `$ 0.10 |
| 34.8 | 1.222 $`\pm `$ 0.052 |
| 48.6 | 1.15 $`\pm `$ 0.12 |
* J.E.R. Holmes et al, 1961

| 6.7 | 1.39 $`\pm `$ 0.08 |
| --- | --- |
| 11.0 | 1.35 $`\pm `$ 0.08 |
| 18.0 | 1.29 $`\pm `$ 0.08 |
| 36.0 | 1.29 $`\pm `$ 0.14 |
| 98.0 | 1.02 $`\pm `$ 0.14 |
* W. Pak et al, 1961

| 13.1 | 1.252 $`\pm `$ 0.029 |
| --- | --- |
| 18.1 | 1.237 $`\pm `$ 0.064 |
| 25.3 | 1.262 $`\pm `$ 0.050 |
| 49.3 | 1.137 $`\pm `$ 0.093 |
* P.J. Hayman and A.W. Wolfendale, 1962

| 10.4 | 1.223 $`\pm `$ 0.038 |
| --- | --- |
| 17.5 | 1.233 $`\pm `$ 0.037 |
| 35.0 | 1.268 $`\pm `$ 0.051 |
| 77.0 | 1.37 $`\pm `$ 0.16 |
| 120.0 | 1.45 $`\pm `$ 0.23 |
| 240.0 | 1.51 $`\pm `$ 0.38 |
* I.C. Appleton et al, 1971

| 12.7 | 1.312 $`\pm `$ 0.039 |
| --- | --- |
| 17.2 | 1.263 $`\pm `$ 0.038 |
| 28.3 | 1.306 $`\pm `$ 0.044 |
| 50.0 | 1.285 $`\pm `$ 0.085 |
| 81.0 | 1.165 $`\pm `$ 0.14 |
| 127.0 | 1.266 $`\pm `$ 0.20 |
| 288.0 | 1.105 $`\pm `$ 0.25 |
* O.C. Allkofer and W.D. Dau, 1972

| 11.4 | 1.22 $`\pm `$ 0.10 |
| --- | --- |
* B.C. Nandi et al, 1972

| 10.8 | 1.263 $`\pm `$ 0.030 |
| --- | --- |
| 15.2 | 1.268 $`\pm `$ 0.040 |
| 19.7 | 1.293 $`\pm `$ 0.057 |
| 26.6 | 1.290 $`\pm `$ 0.050 |
| 37.1 | 1.209 $`\pm `$ 0.079 |
| 46.9 | 1.257 $`\pm `$ 0.104 |
| 60.0 | 1.235 $`\pm `$ 0.114 |
| 84.0 | 1.430 $`\pm `$ 0.172 |
| 142.0 | 1.363 $`\pm `$ 0.188 |
| 260.0 | 1.364 $`\pm `$ 0.383 |
| 566.0 | 1.259 $`\pm `$ 0.460 |
* J.M. Baxendale et al, 1975

| 11.5 | 1.259 $`\pm `$ 0.029 |
| --- | --- |
| 11.9 | 1.269 $`\pm `$ 0.029 |
| 12.3 | 1.229 $`\pm `$ 0.028 |
| 12.8 | 1.371 $`\pm `$ 0.031 |
| 13.3 | 1.219 $`\pm `$ 0.027 |
| 13.9 | 1.261 $`\pm `$ 0.027 |
| 14.5 | 1.292 $`\pm `$ 0.028 |
| 15.3 | 1.291 $`\pm `$ 0.027 |
| 16.1 | 1.299 $`\pm `$ 0.027 |
| 17.0 | 1.256 $`\pm `$ 0.026 |
| 18.2 | 1.319 $`\pm `$ 0.028 |
| 19.5 | 1.279 $`\pm `$ 0.027 |
| 21.2 | 1.283 $`\pm `$ 0.021 |
| 22.7 | 1.300 $`\pm `$ 0.018 |
| 24.1 | 1.293 $`\pm `$ 0.035 |
| 25.2 | 1.301 $`\pm `$ 0.022 |
| 26.3 | 1.289 $`\pm `$ 0.036 |
| 27.8 | 1.285 $`\pm `$ 0.028 |
| 29.3 | 1.286 $`\pm `$ 0.036 |
| 31.5 | 1.325 $`\pm `$ 0.023 |
| 33.1 | 1.324 $`\pm `$ 0.037 |
| 35.3 | 1.300 $`\pm `$ 0.038 |
| 36.8 | 1.262 $`\pm `$ 0.028 |
| 39.6 | 1.304 $`\pm `$ 0.019 |
| 43.2 | 1.283 $`\pm `$ 0.019 |
| 45.2 | 1.258 $`\pm `$ 0.024 |
| 47.1 | 1.222 $`\pm `$ 0.029 |
| 49.3 | 1.265 $`\pm `$ 0.031 |
| 50.7 | 1.358 $`\pm `$ 0.043 |
| 52.1 | 1.361 $`\pm `$ 0.034 |
| 54.9 | 1.250 $`\pm `$ 0.022 |
| 58.2 | 1.271 $`\pm `$ 0.025 |
| 63.0 | 1.223 $`\pm `$ 0.033 |
| 65.3 | 1.312 $`\pm `$ 0.046 |
| 68.1 | 1.272 $`\pm `$ 0.025 |
| 73.7 | 1.286 $`\pm `$ 0.037 |
| 76.6 | 1.323 $`\pm `$ 0.049 |
| 80.4 | 1.235 $`\pm `$ 0.028 |
| 88.3 | 1.215 $`\pm `$ 0.038 |
| 93.0 | 1.238 $`\pm `$ 0.050 |
| 98.3 | 1.245 $`\pm `$ 0.039 |
| 112.0 | 1.243 $`\pm `$ 0.044 |
| 118.0 | 1.270 $`\pm `$ 0.058 |
| 128.0 | 1.287 $`\pm `$ 0.047 |
| 145.0 | 1.268 $`\pm `$ 0.052 |

| 160.0 | 1.250 $`\pm `$ 0.068 |
| --- | --- |
| 177.0 | 1.221 $`\pm `$ 0.056 |
| 214.0 | 1.298 $`\pm `$ 0.071 |
| 236.0 | 1.327 $`\pm `$ 0.101 |
| 274.0 | 1.316 $`\pm `$ 0.093 |
| 358.0 | 1.586 $`\pm `$ 0.154 |
* B.C. Rastin, 1984

| 10.69 | 1.239 $`\pm `$ 0.016 |
| --- | --- |
| 11.94 | 1.247 $`\pm `$ 0.016 |
| 13.58 | 1.251 $`\pm `$ 0.016 |
| 15.81 | 1.285 $`\pm `$ 0.016 |
| 19.05 | 1.263 $`\pm `$ 0.016 |
| 24.14 | 1.233 $`\pm `$ 0.016 |
| 28.35 | 1.267 $`\pm `$ 0.053 |
| 29.30 | 1.166 $`\pm `$ 0.050 |
| 30.32 | 1.182 $`\pm `$ 0.051 |
| 31.42 | 1.250 $`\pm `$ 0.054 |
| 32.60 | 1.325 $`\pm `$ 0.058 |
| 33.88 | 1.253 $`\pm `$ 0.055 |
| 35.27 | 1.277 $`\pm `$ 0.057 |
| 36.79 | 1.238 $`\pm `$ 0.056 |
| 38.44 | 1.252 $`\pm `$ 0.059 |
| 40.25 | 1.206 $`\pm `$ 0.057 |
| 42.25 | 1.277 $`\pm `$ 0.060 |
| 44.47 | 1.291 $`\pm `$ 0.062 |
| 46.94 | 1.388 $`\pm `$ 0.069 |
| 49.71 | 1.289 $`\pm `$ 0.064 |
| 52.84 | 1.160 $`\pm `$ 0.059 |
| 56.40 | 1.273 $`\pm `$ 0.067 |
| 60.49 | 1.240 $`\pm `$ 0.067 |
| 65.23 | 1.310 $`\pm `$ 0.073 |
| 70.80 | 1.157 $`\pm `$ 0.066 |
| 77.42 | 1.203 $`\pm `$ 0.071 |
| 85.43 | 1.256 $`\pm `$ 0.080 |
| 95.34 | 1.207 $`\pm `$ 0.082 |
| 107.88 | 1.320 $`\pm `$ 0.096 |
| 124.27 | 1.194 $`\pm `$ 0.091 |
| 146.62 | 1.235 $`\pm `$ 0.101 |
| 178.85 | 1.161 $`\pm `$ 0.104 |
| 229.36 | 1.234 $`\pm `$ 0.127 |
| 319.72 | 1.30 $`\pm `$ 0.16 |
| 525.82 | 1.32 $`\pm `$ 0.22 |
| 1288.74 | 1.14 $`\pm `$ 0.29 |
* M.P. De Pascale et al, 1993

| 19.89 | 1.292 $`\pm `$ 0.075 |
| --- | --- |
| 55.87 | 1.409 $`\pm `$ 0.173 |
* J. Kremer et al, 1999 1994 data

| 12.42 | 1.257 $`\pm `$ 0.038 |
| --- | --- |
| 18.85 | 1.269 $`\pm `$ 0.055 |
| 26.68 | 1.372 $`\pm `$ 0.086 |
| 36.69 | 1.466 $`\pm `$ 0.110 |
| 51.47 | 1.384 $`\pm `$ 0.144 |
| 72.08 | 1.212 $`\pm `$ 0.169 |
| 100.96 | 1.235 $`\pm `$ 0.233 |
* J. Kremer et al, 1999 1997 data

| 12.42 | 1.291 $`\pm `$ 0.027 |
| --- | --- |
| 18.85 | 1.335 $`\pm `$ 0.038 |
| 26.68 | 1.427 $`\pm `$ 0.058 |
| 36.69 | 1.383 $`\pm `$ 0.068 |
| 51.47 | 1.320 $`\pm `$ 0.090 |
| 72.08 | 1.337 $`\pm `$ 0.120 |
| 100.96 | 1.337 $`\pm `$ 0.169 |

APPENDIX D Correlation matrix for the fit of the spectrum shape:

|  | $`H_3`$ | $`H_2`$ | $`H_1`$ | $`S_2`$ |
| --- | --- | --- | --- | --- |
| $`H_3`$ | 1.000 | 0.319 | -0.298 | -0.771 |
| $`H_2`$ | 0.319 | 1.000 | -0.105 | -0.558 |
| $`H_1`$ | -0.298 | -0.105 | 1.000 | 0.672 |
| $`S_2`$ | -0.771 | -0.558 | 0.672 | 1.000 |


format: text

text:
# A nonlinear detection algorithm for periodic signals in gravitational wave detectors

## Introduction

Kilometer-scale gravitational wave interferometric detectors will become operational by the end of 2001, allowing observations with unprecedented sensitivity at frequencies ranging from a few tens of hertz to approximately one kilohertz . Periodic sources form an interesting subclass of potentially detectable astrophysical objects, due to the possibility to improve the “visibility” of a weak periodic signal in noisy data by increasing the observation time. Good candidates for such sources are spinning galactic neutron stars with a deformation misaligned with their rotation axis; optimistic estimates suggest that LIGO I could detect such objects with an observation time of the order of one year .

Frequency modulations due to the Doppler shift induced by the detector motions, and intrinsic frequency variations due to the loss of angular momentum to gravitational waves, are two easily parameterizable examples of properties of a realistic signal that will dramatically increase the difficulty of the data analysis. In fact, for all but the simplest types of searches, the data analysis will be sub-optimal because of computational limitations. More complications are likely to arise if the source is in a binary, due to the additional frequency modulation  and possibly larger or random intrinsic frequency variations expected if the source is accreting.

These considerations have driven serious efforts to define data analysis schemes that would maximize the sensitivity of a search (to be defined in section I B), given a maximal available computational power. Matched filtering analysis, together with an in-depth discussion of signal properties, is presented in a series of papers ; for all-sky searches, this technique would be computationally prohibitive. The so-called coherent search, discussed in , uses a resampling technique to demodulate the signal and to correct for intrinsic frequency variations: the detection is performed by taking the Fourier transform with respect to the resampling time, which amounts to a change of frame of reference from the detector to the solar system barycenter and to a “stretching” of time for frequency variations. This technique is generally not practical for all-sky searches over long time periods (longer than a few days), but a refinement of it called the stack-slide technique  is considered as one of the best of all known techniques: the initial data are segmented into shorter intervals, and their power spectra are built using the coherent method. For these shorter intervals, a smaller number of points in parameter space have to be used. The spectra are added (incoherently) after having been shifted according to a finer gridding of parameter space. A related technique uses Hough transforms  to track the frequency evolution of peaks in the spectra from the shorter intervals, which is known in advance for some choice of the source parameters. In general, it is possible to improve the sensitivity of these techniques by integrating them into more general hierarchical searches : a first search in parameter space is done at reduced sensitivity but low confidence level (so that spurious events are likely); the vicinity of candidate events is then scrutinized at full sensitivity to get higher confidence detections.

We have explored yet another approach that exploits a particular symmetry of the problem in an attempt to reduce its complexity. Firstly, we note that the Doppler shift in the frequency of a given source due to the rotation of the earth, independently of the source position, is exactly equal and opposite at any two instants of time separated by half a sidereal day (we will only consider linear terms in the Doppler shift, an excellent approximation). Secondly, in the approximation that the earth orbital motion around the sun is circular, the Doppler shift in frequency has the same property as above, but for a time separation of half a sidereal year. We shall use this approximation for illustrative purposes; corrections for the small eccentricity of the earth orbit and for the influence of the moon and of other planets would have to be included in practice. Signal phases at different times can be added by properly multiplying the (analytic extensions of the) signal time series, and therefore a signal free of frequency modulations may in principle be built. We will give a more precise formulation of this idea in section I.

The multiplications of different stretches of the data will strongly affect the noise characteristics, principally by mixing noise at different frequencies; a useful implementation will have to bandpass filter the data to minimize this effect, and therefore a search over source parameters (frequency, frequency derivatives, source position) will be required. An example of a possible implementation will be presented in section I A. The achievable sensitivity (to be defined in section I B) will be studied numerically as a function of computational power requirements in section II, and comparisons with the stack-slide method mentioned above will be presented.

## I Algorithm

We will use the following model for the signal $`s(t)`$ at the detector:

$`s(t)=h\mathrm{cos}[\varphi (t)],`$ (1)
$`\varphi (t)=2\pi f_0[t+M_1(t)+M_2(t)],`$ (2)
$`M_1(t)={\displaystyle \frac{A_1\mathrm{d}}{2\pi }}\mathrm{cos}\left({\displaystyle \frac{2\pi t}{\mathrm{d}}}+\alpha \right),`$ (3)
$`M_2(t)={\displaystyle \frac{A_2\mathrm{yr}}{2\pi }}\mathrm{cos}\left({\displaystyle \frac{2\pi t}{\mathrm{yr}}}+\psi \right),`$ (4)
$`A_1={\displaystyle \frac{2\pi R_{}}{\mathrm{d}c}}\mathrm{cos}(\lambda )\mathrm{cos}(\delta ),`$ (5)
$`A_2={\displaystyle \frac{2\pi \mathrm{A}.\mathrm{U}.}{\mathrm{yr}c}}\mathrm{cos}(\delta +ϵ),`$ (6)
where $`ϵ`$ is the earth oblicity, $`\lambda `$ is the detector latitude, $`c`$ is the speed of light, $`R_{}`$ is the earth radius, A.U. is an astronomical unit, d is a sidereal day length, yr is a sidereal year length, $`(\alpha ,\delta )`$ is the source angular position, $`f_0`$ is the signal frequency in the source rest frame (it may be a function of time), and $`\psi `$ is a fixed number that defines the origin of time. For simplicity, since our focus is on the phase modulation of the signal, the explicit time dependence of $`h`$ on the detector response pattern will not be considered, although its averaged effect will be included in our definition of sensitivity (section I B). If $`f_0`$ is independent of time, a good approximation to the signal half-bandwidth $`\mathrm{\Delta }f`$ over a time interval of length $`t`$ (i.e., half the intrinsinc frequency change over $`t`$), is given by:

$$\mathrm{\Delta }f(t)=f_0[A_1\mathrm{min}(1,t/\mathrm{d})+A_2\mathrm{min}(1,t/\mathrm{yr})].$$
(7)
The numerical values of the terms multiplying the cosine terms in eq. 5 and eq. 6 are approximately $`1.5510^6`$ and $`9.9410^5`$, respectively. Hence, the half-bandwidth of a 1 kHz source would be roughly 2 mHz over one day, and 0.1 Hz over a year.

The phase function $`\varphi (t)`$ is sufficiently well-behaved that we can use the Hilbert transform of $`s(t)`$ to build its quadrature , $`\widehat{s}(t)=h\mathrm{sin}[\varphi (t)]`$, and therefore construct the analytic signal $`S(t)`$:

$$S(t)=he^{i\varphi (t)}.$$
(8)
Assuming that $`f_0`$ is independent of time, we have that

$`S_4(t)=S(t)S(t+\mathrm{d}/2)S(t+\mathrm{yr}/2)S(t+\mathrm{d}/2+\mathrm{yr}/2)`$ (9)
$`=h^4e^{8\pi if_0(t+\mathrm{d}/4+\mathrm{yr}/4)},`$ (10)
that is, $`S_4`$ is monochromatic with frequency $`4f_0`$. If $`f_0`$ is a function of time, however, $`S_4`$ as defined above will not be demodulated because of the coupling of $`f_0`$ with the modulations; its phase will be

$`{\displaystyle \frac{\varphi _4(t)}{2\pi }}=f_0(t)t+f_0(t_1)t_1+f_0(t_2)t_2+f_0(t_3)t_3+`$ (11)
$`M_1(t)[f_0(t)f_0(t_1)]+M_1(t_2)[f_0(t_2)f_0(t_3)]+`$ (12)
$`M_2(t)[f_0(t)f_0(t_2)]+M_2(t_1)[f_0(t_1)f_0(t_3)],`$ (13)
where $`t_1=t+\mathrm{d}/2`$, $`t_2=t+\mathrm{yr}/2`$ and $`t_3=t+\mathrm{d}/2+\mathrm{yr}/2`$. The signal so obtained still presents some residual phase modulation from the detector motions, but it results in a much smaller contribution $`\mathrm{\Delta }f_{\mathrm{mod}}`$ to the signal bandwidth; the terms multiplying the modulation functions $`M_1`$ and $`M_2`$ are now of order of the intrinsic frequency change over half a day and half a year, respectively, rather than the source frequency itself. We will show below that our algorithm best performances are obtained by dividing the total data from an observation of length $`t_{\mathrm{obs}}`$ into smaller segments of length $`T`$ shorter than one day. Consequently, we impose the condition that the signal bandwidth produced by the residual modulations in eq. 13 be smaller than the frequency resolution for a segment of length $`T`$, $`\mathrm{\Delta }f_{\mathrm{mod}}<1/T`$, assuming a linear frequency model: $`f_0(t)=f_0[1t/\tau ]`$, where $`\tau `$ is the source spin-down characteristic time. We find that this condition is:

$$\tau >f_0T^2(A_1+A_2),$$
(14)
i.e.,

$$\tau >\left(\frac{f_0}{1\mathrm{k}\mathrm{H}\mathrm{z}}\right)\left(\frac{T}{1\mathrm{d}}\right)^225\mathrm{y}\mathrm{r}.$$
(15)
This estimate shows clearly that for most realistic sources, the demodulation by multiplication is performed correctly, so that we can assume that spin down enters the phase function only through the combination $`f_0(t)t+f_0(t_1)t_1+f_0(t_2)t_2+f_0(t_3)t_3`$.

Noise will be seriously affected by the multiplications of data segments at different times, through nonlinear mixing of signal with noise and of noise with noise at different frequencies. The latter of these effects is so severe that we will have to minimize it as much as possible by bandpass filtering individually the data segments around the signal, or the technique will be of no utility. An efficient filtering requires a knowledge of the frequency content of the signal as a function of time, and will therefore involve a search over source frequencies, sky positions and spin-down parameters; this will greatly increase the computational costs. This search will be implemented by meshing this parameter space with a certain resolution, i.e. with a certain number of points $`N_p`$. Obviously, if a given signal has parameters that are not exactly coincident with any of these points, the sensitivity of our search to this signal will be reduced. We will explore below how the density of the mesh couples to the achievable mean sensitivity, and therefore give an estimate of the usefulness of our technique.

### A Implementation

From the general guidelines established in the previous section, we are able to discuss an implementation of our algorithm. This implementation has parameters that allow an optimization of the method, subject to computational power constraints.

The first step will be to prepare for the analysis data obtained from an observation of length $`t_{\mathrm{obs}}`$. The dataset is compressed to contain only information for frequencies of interest, i.e., frequencies between $`f_{\mathrm{min}}`$ and $`f_{\mathrm{max}}`$. The analytic signal is then constructed from this reduced dataset, using Hilbert transforms to build the quadrature, as described in section I. The resulting data is divided into equal segments of length $`T`$. The computational cost of these operations is negligible, especially because they are only performed once at the beginning of the search.

Next, a bank of bandpass filtered segments is built. We choose a filter half-bandwidth $`B>\mathrm{\Delta }f`$, where $`\mathrm{\Delta }f`$ is the maximum of the signal half-bandwidth over time $`T`$, and for every segment, we construct $`N_f`$ smaller filtered segments by using a heterodyne technique . These filtered segments are arranged to cover completely the frequency interval from $`f_{\mathrm{min}}`$ to $`f_{\mathrm{max}}`$. If we impose the requirement that no signal power should be lost due to the filtering, the number $`N_f`$ is a function of $`B\mathrm{\Delta }f`$: the whole frequency interval is divided into subintervals of length $`2B`$, displaced from each other by $`2(B\mathrm{\Delta }f)`$, so that any signal with half-bandwidth $`\mathrm{\Delta }f`$ is guaranteed to be fully contained within one subinterval, and $`N_f=(f_{\mathrm{max}}f_{\mathrm{min}})/2(B\mathrm{\Delta }f)`$. In the heterodyne technique, the signal is multiplied by a complex exponential $`\mathrm{exp}2\pi i(Bf_0)t`$, low-pass filtered (for reasons of antialiasing), and resampled at frequency $`8B`$. The resulting segments have $`8BT`$ (complex) data points, and their Nyquist frequency is $`4B`$. The cost of the filtering, if implemented by an infinite impulse response digital filter, is of order $`2(2n1)`$ floating point operation per input point, where $`n`$ is the number of poles in the filter, which is a function of the filter performance; we shall denote this cost by $`C_\mathrm{f}`$. Resampling will have negligible computational cost compared to the multiplication by the complex exponential; assuming that the latter is not computed at run time, the cost for the construction of the filtered segments bank is

$$(C_\mathrm{f}+6)t_{\mathrm{obs}}\frac{(f_{\mathrm{max}}f_{\mathrm{min}})^2}{2(B\mathrm{\Delta }f)}\mathrm{flop}.$$
(16)
This cost can be adjusted by varying $`B`$ and $`T`$, which enters the formula through the dependence of $`\mathrm{\Delta }f`$ on it.

Once the filtered segments bank is constructed, for each of the $`N_p`$ points in parameter space, and for every segment at time $`t<t_{\mathrm{obs}}T\mathrm{d}/2\mathrm{yr}/2`$, we select the appropriate filtered segments from the bank, i.e. filtered segments where the signal is expected for times $`t`$, $`t_1=t+\mathrm{d}/2`$, $`t_2=t+\mathrm{yr}/2`$ and $`t_3=t+\mathrm{d}/2+\mathrm{yr}/2`$, and we multiply them together. At this point, the modulations induced from the detector motion have been removed; what remain to be applied are corrections for the spin-down of the source. A refined gridding of the spin-down portion of parameter space, for the already quite narrow region used for the choice of filtered segments from the bank, is used to resample the product of the four segments according to the new time

$$t^{}=\frac{tf_0(t)+t_1f_0(t_1)+t_2f_0(t_2)+t_3f_0(t_3)}{f_0(0)},$$
(17)
where the origin of $`t`$ was arbitrarily set to zero. The resampling will only affect a few points in the time series, so its cost is negligible in itself. The power spectrum is constructed by taking the norm of the Fourier transform. Finally, all the spectra corresponding to a given point in parameter space (from the partitioning of time in segments of length $`T`$) are added together to form the final spectrum, which is searched for significant peaks up to the Nyquist frequency $`4B`$. All these operations are done on small segments of length $`8BT`$, and the total cost is the sum of the cost for building the products for all $`N_p`$ points and of the cost for the spin-down corrections, power spectra construction and sum for the $`N_p^{}`$ points corresponding to the finer grid:

$$Bt_{\mathrm{obs}}N_p[72+20N_p^{}(1+\mathrm{log}_28BT)]\mathrm{flop}.$$
(18)
Again, this can be adjusted by varying $`B`$ and $`T`$, since $`N_p`$ and $`N_p^{}`$ depend on both $`B`$ and $`T`$. The total cost of the search is the sum of equations 16 and 18.

The number $`N_p`$ will be the number of non-intersecting (hyper-) volume elements in parameter space required to cover completely the region of interest, such that all points within a given volume element correspond to the same choice of filtered segments from the bank, at each of the four times of interest in this problem. Manifestly, $`N_p`$ will depend on $`B`$ through the way it partitions the frequency interval $`[f_{\mathrm{min}},f_{\mathrm{max}}]`$, and on $`T`$, which influences the size of a volume element over which the source frequency varies by some fixed amount. We solve the problem of computing $`N_p`$ numerically, by building a mean volume element: we select at random a point in parameter space, and compute its associated volume element. Repeating this procedure, we construct the average of a volume element, and divide the parameter space volume by it; the result is an estimate of $`N_p`$.

The refined gridding should only be necessary for searches with large spin-down rates (see below), and its effect will be to force to repeat $`N_p^{}`$ times the operations following and including the resampling above, where $`N_p^{}`$ is the number of points on the refined grid. The computation of $`N_p^{}`$ is not of the same nature as the one for $`N_p`$: the problem now only involves the spin-down parameters, and the goal is to minimize the losses in signal power from residual frequency drifts, or misalignments of the spectra that are summed together. In that sense, it is similar to the problem of computing the number of points in parameter space for the stack-slide technique , and should be solved using the same geometrical approach. However, we shall rather take $`N_p^{}=1`$, and justify this approximation by the relatively large spin-down time ($`\tau =5000\mathrm{y}\mathrm{r}`$) to be used in the numerical example of section II. It should be kept in mind that this may not be sufficient for faster spin-down, in the sense that under this approximation $`N_p^{}`$ may be under-estimated. With the ability we now have to compute the required computational cost as a function of $`B`$ and $`T`$, the only thing that remains to be described is how we characterize the sensitivity of the search, as a function of the same choice of parameters.

### B Sensitivity

From a statistical point of view, the results in  are presented in a somewhat non-standard way, so to ease comparison of our results with theirs, we first briefly reformulate their definition of sensitivity. By setting the average noise power to be equal to one, we can express the signal amplitude $`h`$ in units of $`\sqrt{S_n(f)/t_{\mathrm{obs}}}`$, where $`S_n(f)`$ is the noise spectral density at frequency $`f`$. In order to be consistent with , we use the definition that the averaged signal power (also the square of the signal to noise ratio in amplitude) is $`P_s=F_+^2h^2`$, where $`F_+^2=1/5`$ is the average of the detector response over all angles and polarizations . This definition is useful to approximately account for variations in $`h`$ over time that are produced by the motion of the source in the detector response pattern. For a true signal power $`P_s`$, computing the power spectrum of a time series as a sum of $`n`$ of its subseries spectra (this describes the statistics of both the coherent and stack-slide techniques) will give a distribution $`p_n(P|P_s)`$ of the power $`P`$ . Given a false alarm probability $`1\alpha `$, one defines a threshold power $`P_\alpha `$ such that the integral of $`p_n(P|0)`$ from 0 to $`P_\alpha `$ is $`\alpha `$. Observations with power exceeding $`P_\alpha `$ will be called “detections with confidence level $`\alpha `$***The authors of  define the false alarm probability as $`1\alpha /N_p`$, where $`N_p`$ is the number of points in parameter space.. One further defines a false dismissal probability $`\beta `$ as a function of the signal power, as the integral of $`p_n(P|P_s)`$ with respect to $`P`$, from 0 to $`P_\alpha `$. It has been customary in the field  to use the relation $`P=P_s+n`$ for the average observed power $`P`$ to interpret $`P_\alpha `$ as a lower bound on the detectable signal amplitude, by setting $`P_\alpha =P`$. For instance, performances of sub-optimal techniques are frequently compared to the optimal case defined by $`n=1`$; for $`\alpha =0.99`$, the minimal detectable amplitude as defined above (the optimal amplitude) is $`h=4.2`$, and the corresponding false dismissal rate is $`\beta 0.57`$. In fact, it is easily checked that as $`n`$ gets larger, this definition gives a value of $`\beta `$ that approaches $`1/2`$ from above. Therefore, for the ease of comparison with other works, we will define our minimal detectable amplitude as the smallest value of $`h`$ that corresponds to a false alarm probability $`1\alpha `$ and a false dismissal rate $`\beta =1/2`$. The sensitivity will be defined as $`\mathrm{\Theta }=4.2/h`$, and this definition will be approximately compatible with those of .

We use numerical analysis to determine $`h`$ as a function of $`B`$ and $`T`$. We first set the signal amplitude to zero, and generate a large number of filtered segments of pure white noise, that are multiplied together and then added by groups of $`(t_{\mathrm{obs}}T\mathrm{d}/2\mathrm{yr}/2)/T`$ (i.e., the number of segments in the portion of the observation period that can be used for multiplications). By constructing the cumulative probability function for the resulting segments, we then evaluate the threshold $`P_\alpha `$. We repeat the procedure for non-zero signal amplitudes: picking at random in a uniform distribution the parameters of the source, we construct the final spectrum, and evaluate the power in the expected frequency channel. For each choice of signal amplitude, we repeat this a large number of times for different source parameters, and build a distribution of the power in the expected channels; the values of $`\beta `$ as a function of signal power are then deduced from these distributions and the previously calculated value of $`P_\alpha `$. Therefore, the sensitivity $`\mathrm{\Theta }`$ is effectively averaged over all possible source parameters. To reduce computational burden, spin-down effects are not directly included, in the sense that we only consider them in determining the bandwidth of the signal for a given $`T`$, and do not include them in the initial simulated segments that are multiplied together. By doing so, we might reduce the non-linear mixing of signal with noise, and therefore deduce a sensitivity that is too high; however, the magnitude of this mixing is such that the error so induced is expected to be small.

## II Numerical Simulations

We opted for a relatively “easy” broad-band all-sky search for our explorations. We choose $`f_{\mathrm{min}}=40\mathrm{H}\mathrm{z}`$, $`f_{\mathrm{max}}=1\mathrm{k}\mathrm{H}\mathrm{z}`$, $`\tau =5000\mathrm{y}\mathrm{r}`$, $`t_{\mathrm{obs}}=1\mathrm{y}\mathrm{r}`$, and confidence level $`\alpha =99\%`$. Note that $`t_{\mathrm{obs}}=1\mathrm{y}\mathrm{r}`$ is the minimal possible observation time, and that the data at times $`t>\mathrm{d}/2+\mathrm{yr}/2`$ can not be demodulated. This obviously limits the sensitivity of the search, although it does so in a diminushing proportion as $`t_{\mathrm{obs}}`$ grows above one year. A different narrow-band all-sky search was also considered, with parameters as for the broad-band search, except for $`f_{\mathrm{min}}=450\mathrm{H}\mathrm{z}`$ and $`f_{\mathrm{max}}=500\mathrm{H}\mathrm{z}`$, perhaps corresponding to the case of an interferometer made narrow-band using signal recycling. For definitiveness, we choose the latitude of the detector to correspond to the LIGO Hanford detector. As described above, we computed numerically the value of $`N_p`$ for different values of $`T`$ and $`B`$ (from now on, $`B`$ will be expressed in units of $`\mathrm{\Delta }f(T)`$ as defined in eq. 7, with $`f_0=f_{\mathrm{max}}`$, and $`A_1`$, $`A_2`$ averaged over all angles $`\delta `$); see figure 1. Using the method described in section I B, we computed the dependence of the relative sensitivity, $`\mathrm{\Theta }`$, on the same variables; the results are shown in figure 2.

To get a meaningful notion of the computational power $`P`$ from the cost of the search (the sum of equations 16 and 18), we imposed the requirement that the analysis should be done in a time $`t_{\mathrm{obs}}`$. Our map of $`N_p`$ vs. $`B`$ and $`T`$ was then reexpressed as a map of $`P`$ versus $`B`$ and $`T`$. We observed that for plausible filter implementations, the contribution to $`P`$ from eq. 18 was at least two orders of magnitude larger than the contribution from eq. 16 in the region of interest of the $`T,B`$–plane; our results are highly insensitive to the real cost of constructing the filtered segments bank. Finally, we optimized our algorithm for a given computational power by choosing the values of $`B`$ and $`T`$ that maximize $`\mathrm{\Theta }`$ along a line of constant $`P`$. The results are shown in figure 3, together with the corresponding sensitivity for the stack-slide technique.

### A Discussion

Fig. 3 shows that the sensitivity of the algorithm described in this paper is roughly 30% of the sensitivity of the stack-slide technique in a broad-band search, for computational power ranging from 0.4 Tflops to 8 Tflops. For the same range of computational power, the minimum amplitude detectable with 50% efficiency ranges approximately from 25 to 60 times the optimal amplitude.

Also shown in fig. 3 is the sensitivity of an easier narrow-band all-sky search. When the computational power reaches 4 Tflops, the sensitivity eventually gets better than that of the stack-slide technique, as computed for the harder, broad-band search. However, it is not clear how the performances of the stack-slide technique scale with the bandwidth of the search, although they do scale strongly with the upper frequency. It is not possible, for instance, to apply the stack-slide technique on a heterodyned dataset without having to search over the additional source frequency parameter, because of the coupling of the source frequency to both the detector rest frame time and to the phases introduced by the detector motion ($`M_1`$ and $`M_2`$ in equation 2). The comparison of our algorithm to the stack-slide technique is therefore not trivial, and in particular there might be some narrow-band searches where the stack-slide performances are approached or outperformed.

Our results suggest that the distortions in the noise produced by the nonlinearities inherent to our algorithm can not be alleviated by narrow-banding the signal such that the algorithm performs in a manner that is satisfactory for a realistic broad-band all-sky search. The need for narrow-banding adds one dimension (the source frequency) to the parameter space, and this greatly increases the computational burden. Moreover, the dependence of the sensitivity on the bandwidth of the bandpass filter, which is the most important parameter in determining the number of points in parameter space, is strong enough that efficient detection is computationally intensive. Consequently, our principal conclusion is that, although the symmetry described in section I simplifies the problem from a formal point of view, it does so by strongly increasing the level of noise, principally by mixing components of noise of different frequencies. We have presented an algorithm that minimizes this effect, but have found that for a realistic broad-band search, it did not perform better than the stack-slide technique at the same computational power.

Most detection algorithms for periodic sources, and in particular the stack-slide technique, probably don’t have a strong scaling of their performances with the bandwidth of the search, essentialy because making use of the reduced bandwidth adds complexity by forcing one to consider the frequency as an explicit parameter of the search. It is therefore possible that our algorithm compares more advantageously to these other algorithms in narrow-band searches, which can be considered as more natural problems for it, because it necessarily involves a search over the source frequency.

###### Acknowledgements.

The author wishes to thank Rainer Weiss for helpful discussions and comments on the work presented here. This work was supported by the National Science Foundation under cooperative agreement PHY-9210038, and by the Fonds pour la Formation de Chercheurs et l’Aide à la Recherche of the Province of Québec, Canada.