3.1 GLE N 68 detected by Aragats Multichannel Muon Monitor

3.2 nalysis of the Residuals (Checking the Gaussian Model)

3.3 Calculation of the Chance Probability

3.4 Effect of the Multiple Attempts in Searches of “Biggest” Deviation From H₀

3.1 GLE N 68 detected by Aragats

Multichannel Muon Monitor

At middle - low latitudes the enhancements of count rates due to additional particles coming from the Sun (Ground Level Enhancements, GLE) usually not exceed 1-2%. Nonetheless reliable detection of these enhancements is extremely important, because only at middle low latitudes it is possible to detect highest energy solar cosmic rays and determine the maximal energy of solar accelerators.

The problem of proving existence of not very large short peaks in time series is extremely complicated for the reason of fluctuation of mean values of particle detector count rates due to numerous random and systematic effects. Among them are atmospheric effects, earth natural radiation effects, disturbances of the magnetosphere and interplanetary magnetic field (IMF) and, of course, instrumental errors.

We will demonstrate the procedures of estimating significance of the observed peaks in time series by the detection of GLE N 68 with particle detectors of Aragats Space Environmental Center (ASEC, Chilingarian et al., 2005). On January 20, 2005 NOAA reported an X7 importance flare with helio-coordinates (14N, 67 W), which started at 6:36 UT with maximal X-ray flux at 7:01 UT. The first results from the space-born spectrometers for the proton energies up to 800 MeV pointed to very hard energy spectra of the Solar Energetic Particle (SEP) event.
Middle and low-latitude neutron monitors cannot be used for the reconstruction of the primary energy spectra well above 10 GeV due to very weak fluxes, rather small sizes of the detectors and overwhelming flux of low energy solar protons. The Aragats Multidirectional Muon Monitor (AMMM) is located at (40.5°N, 44.17°E), altitude 3200 m. a.s.l.; cutoff rigidity 7.6 GV; relative accuracy of measuring 3-minute time-series ~0.17%, much better than of neutron monitor 18NM64, located at the same altitude. The AMMM in 2005 was comprised of 45 plastic scintillators with detecting surface of 1 m2 and thickness of 5 cm each. The Detector is under 14 meters of soil and concrete, plus 12 cm. of iron bars. Only muons with energies greater than 5 GeV can reach this underground detector. 5 GeV muons correspond to an ensemble of primary protons with mode energy ~50 GeV if we assume the power law with spectral index of γ=-2.7 (Galactic Cosmic Rays); and - most probable energy of primary protons is between 23-30 GeV if we assume power index between γ= -4 - -7. The 1- minute time series of the AMMM are presented in Figure 1. Enhancement of the count rate is seen from 7:02 till 7:04 UT with maximum at 7:03 UT. Three out of the 45 one m2 scintillators of the AMMM were not operational at the time, therefore only 42 m2 of muon detectors were in use to measure the high energy muon flux. The estimated mean count rate of the Galactic Cosmic Rays (GCR) as measured by the 42 m2 of the AMMM detector during the 6:30 – 7:35 UT time span, excluding the enhanced interval from 7:02 to 7:04 UT, is 123818 particles per minute. The additional signal at 7:03 UT equals to 863 particles (0.70%). Taking into account that the standard deviation of 1 minute data is 352 (0.29%) we come to the significance of 2.5 σ for the 1 minute peak at 7:03 UT.

Figure 1. 1- minute count rate of secondary muons >5 GeV by AMMM.

To emphasize the peak in the AMMM time series we group the 1 minute date in 3 minute time-intervals (see Figure 2). The 3-minute time series demonstrates a more pronounced peak of 3.93 σ. The mean count rate of GCR again was estimated in the time span of 6:30 – 7:35 UT with the excluded interval of 7:01 – 7:05 UT and it equals 371494 particles per 3 minutes for the 42 m2. The additional signal at 7:02-7:04 UT equals 2354 (0.644%). If we adopt the Poisson standard deviation for the 3 minute time series 0.164% (see detailed discussion on the method of the estimating of the standard deviation and peak significance in the attachment) we come to the significance of 3.93 σ for the 3 minute peak at 7:02 – 7:04 UT. The excess count rate registered at AMMM during the interval 7:02-7:04 UT corresponds to the flux 3.1( +/- 0.8) 10^-5 muons/cm²/sec.
Due to the very short enhancement time span no corrections on the atmospheric pressure and temperature are necessary.

Figure 2.3- minute count rate of secondary muons >5 GeV of the AMMM, expressed in terms of standard deviations.

Due to upmost importance of the detection of solar particles with highest energies we discuss statistical methods used to reveal peak in the time series of AMMM.

3.2 Analysis of the Residuals

(Checking the Gaussian Model)

The difficulty of testing hypothesis of Gaussian nature lies in the slow drift of the mean count rate of time series due to systematic changes of several geophysical and interplanetary parameters. Disturbances of the Interplanetary Magnetic Field (IMF) in the end of January 2005 (triggered by passage of several Interplanetary CMEs at 16-20 January) modulate cosmic ray flux, introducing trend in the secondary cosmic ray fluxes.
To account for the changing mean of the greater than 5 GeV muon flux we calculate the hourly mean count rates and corresponding residuals (fitting errors, differences between observed hourly means and values of 3-minute count rates in this hour; 20 numbers for each hour):

where X_i,j are normalized residuals, C_i,j are 3 minute count rates of the AMMM at j^th hour, are hourly means of the 3-minute time series and

are proxies of root mean square errors and N_h. is number of hours. Statistical distribution (1) represents, so called, multinomial process. Multinomial process consists of sum of j Gaussian random processes; in our case – time series of count rates corresponding to Gaussian process with same variance and different means. In our probabilistic treatment of the problem we normalize time series by the “moving” means and variances s_j², estimated each hour. In this way we plan to obtain a proxy of the standard Gaussian distribution N(0,1) to use later on as a test statistics.

To check our assumptions on Gaussian nature of the distribution (1) we perform calculation of residuals for 20 January 2005 and for whole January 2005. As we describe in (Bostanjyan et al., 2007) we prepare 3-minute time series from the 1 minute ones. Joining 1 minute time series in 3, 5, 10 or 60 minute time series is ordinary operation used by the all groups running the neutron and muon monitors. To account for the arbitrary choice of the start minute we integrate other all 3 possibilities of different starts of the 3-minute time series, therefore number of events in histograms is 3 times more than number of 3-minute count rates.
The resulting histograms of the normalized residuals are shown in the Figure 2 and 3. We see rather good agreement with standard normal distribution N(0,1); values of the χ² test are ~ 1 for degree of freedom. The maximal values of 3.77¹ (see the right tail of histogram in Figure 2 ) corresponds to a peak at 7:02-7:04 UT. The same maximal value remains maximal also for the one-month histogram (Figure 3). The second maximal value for a month histogram is 3.64.
Proceeding from good agreement of histogram with Gaussian curve and from rather large value of the biggest residual, we can accept the hypothesis that there is additional signal superimposed on the galactic cosmic ray background. Of course, within validity of the Gaussian hypothesis this and larger values can encounter by chance, therefore we’ll need additional statistical tests proving that detected peak is caused by the highest energy solar protons

Figure 3. Normalized residuals calculated by 3-minute AMMM time series at 20 January 2005.

In the picture legend are posted the histogram mean and RMS and also fitted curve mean and variance, as well as number of degrees of freedom in the χ² test.

Figure 4. Normalized residuals calculated by 3-minute AMMM time series at January 2005.

In the picture legend are posted the histogram mean and RMS and also fitted curve mean and variance, as well as number of degrees of freedom in the χ² test.

3.3 Calculation of the Chance Probability

As usually in statistical hypothesis testing, the hypothesis we want to check (named H₀ ) consists in the opposition to the hypotheses we are interested, i.e. we will check the hypothesis that there is no additional muons in 3-minute time series (“no-signal” hypothesis) and, therefore, that detected peak is random fluctuation only. To prove the existence of signal, we have to reject H0 with the maximal possible confidence. Detecting large deviations from H₀ , i.e. very low probability of H₀ being true, do not imply that the opposite hypothesis is automatically valid. As was mentioned by P.Astone and G.D’Agostini (Astone et.al., 1999) behind logic of standard hypothesis testing is hidden a revised version of the classical proof by contradiction. “In standard dialectics, one assumes a hypothesis to be true, then looks for a logical consequence which is manifestly false, in order to reject the hypothesis. The ‘slight’ difference introduced in ‘classical’ statistical tests is that the false consequence is replaced by an improbable one”.
If the experimental data will not differ significantly from test distribution obtained under assumption of “no-signal” hypothesis there will be no reason to reject H₀ and therefore we can’t claim that AMMM detected high energy muons of “solar origin”. And if we will be able to reject H0, we can accept with definite level of confidence that there are high energy protons coming from the sun. Usually confidence level is enumerated as “chance probability”, the probability of H₀ hypothesis to be true.
The statistical test for accepting or rejecting hypothesis is based on the maximal deviation from most probable value (3.77 in our case) observed in time series. The probability to obtain this or another maximal deviation depends on the number of events considered, i.e. on the time series length. Therefore, the most appropriate test provides the extreme statistics distribution (Chapmen et al., 2002, Chilingarian et al., 2006):

Where g (x) is standard Gaussian probability density N (0,1) .

is, so called, standard Gaussian distribution’s p–value: the probability to obtain the value of test statistics greater than x ; M is number of attempts we made to find the biggest deviation from H0 (number of elements of considered time-series multiplied by number of attempts we made to find greatest deviation).
To obtain the probability to observe extremely deviation equal to x among M identically distributed random variables (p-value of the distribution c_M (x) we have to integrate c_M (x) in the interval [x, +∞):

C_x^M(x), p-value of the distribution (3), equals to probability that observed test statistics x maximally deviates from the most probable value under assumption that H₀ is valid. And if this probability is low enough we can reject H0 and accept alternative hypothesis that observed deviation is not fluctuation, but a contamination of the distribution of different statistical nature, i.e., a signal.
The probability to observe in one from 480 (i.e. during the day) of 3-minute time-series count rate enhancement of 3.77 equals according to equations (3-5) to:

It means that in absence of any signal when examining daily variations of the 3-minute count rates in one case from 10 it is expected to detect the deviation of the mean value equal to 3.77 . Equivalent statement: approximately once in 10 days only we will detect 3.77 enhancements in the 3-minute time series of AMMM.
However, we have to correlate the expected signal from protons, accelerated at Sun with time of X-ray flare and CME launch. Of course, we cannot expect the signal from solar protons before X-ray flare and an hour after the X-ray flare or/and CME launch occurs. The chance probability to detect a deviation equals to 3.77 σ in one hour equals to
, i.e. only once in 200 cases we can expect such enhancement.
As we can see in the Figure 3 the second maximal monthly deviation equals to 3.64. If we accept hypothesis that 3.77 value was due to solar protons, we have to check if 3.64 is typical monthly maximal deviation. Calculated according (3-5) value of is rather large and we have no reasons to reject H₀ ; i.e. at January 2005 the residual distribution (Figure 3) was Gaussian with only one outlier attributed to high energy solar protons.

3.4 Effect of the Multiple Attempts in Searches

of “Biggest” Deviation From H₀

To check assumption that when calculated significance of signal we should take into account 3 possible starts of time series we perform simulations with simple model of time series.
The model can be described as following:

1. generate 1440 numbers from the standard normal distribution N(0,1);
2. prepare 3 time series summing 3 consequent numbers of the raw, starting from the first, second and the third elements;
3. perform according to equation (1) normalization procedure (subtract the mean and divide to root of variance);
4. determine and store the maximal element of each of 3 normalized time series;
5. determine and store the maximal element among 3 time-series maximums;
6. repeat i-vi 1000 times and consider 4 histograms of extreme values;
7. calculate the frequencies of obtaining values equal or greater than 4 (for simplicity we take 4 instead of 3.77).

Intuitively, when having 3 possibilities physicist will choose one that emphasis the presence of signal (the situation (v) ). But as we can see from the Figure 5d, the probability to obtain the fake signal is dramatically enhanced (approximately by 3 times). From the same picture we can see that obtained in d) chance probability 0.041 is in good agreement with value calculated according to equations (3-5):

Figure 5. Histograms of the extreme statistics. a) - c) – selecting extreme statistics for 3 independent time series (iv);and d) - selecting maximal value among 3 extreme statistics – (v. Black area in the histograms denotes the summation region and number the integral (sum) value from 4 till infinity

References

• P. Astone, G. D’Agostini, Inferring the intensity of Poison processes at the limit of the detector sensitivity, CERN-EP/99-126, 1999.
• N.Kh. Bostanjyan , A.A. Chilingarian, V.S. Eganov, G.G. Karapetyan. ,On the production of highest energy solar protons on 20 January 2005, Advances in Space Research 39, 1454–1457, 2007
• S.C. Chapman, G. Rowlands, N.W.Watkins, Extremum statistics – a Framework for Data Analysis, Nonlinear Processes in Geophysics, 9, 409-418, 2002.
• A.Chilingarian, K.Avakyan, V. Babayan, et. al, Aragats Space-Environmental Center: Status and SEP Forecasting Possibilities. Journal of Physics G:Nucl.Part.Phys., 29, 939-952, 2003.
• A. Chilingarian, K.Arakelyan, K.Avakyan, et al., Correlated measurements of secondary cosmic ray fluxes by the Aragats Space Enviromental Center monitors", NIM, A543, 483-492, 2005.
• A.Chilingarian, G.Gharagyozyan, G.Karapetyan, G.Hovsepyan, Statistical Methods for Signal Estimation of Point Sources of Cosmic Rays, Astroparticle Physics, 25, 269-276, 2006.
• A. Chilingarian, A.Reymers, Particle detectors in Solar Physics and Space Weather Research, Astroparticle Physics (2007), Astropart. Phys., 27, 465-472
  
• S.N. Karpov, Z.M. Karpova, Yu.V. Balabin and E.V. Vashenyuk, Proc. 29th ICRC, Pune, India, vol. 1, 193-196, 2005.
  
• H. Miyasaka, E. Takahashi, S. Shimoda, et al., Proc. of of 29th ICRC, Pune, India 1, 245-2248, 2005.
  
• F.R. Zhu, Y.Q. Tang F.R., Y. Zhang et, et al., Proc of 29th ICRC, Pune, India, 1, 185-188, 2005.