Monitoring and forecasting of great radiation hazards for spacecraft and aircrafts by online cosmic ray data

Abstract. We show that an exact forecast of great radiation hazard in space, in the magnetosphere, in the atmosphere and on the ground can be made by using high-energy particles (few GeV/nucleon and higher) whose transportation from the Sun is characterized by a much bigger diffusion coefficient than for small and middle energy particles. Therefore, high energy particles come from the Sun much earlier (8-20 min after acceleration and escaping into solar wind) than the main part of smaller energy particles (more than 30-60 min later), causing radiation hazard for electronics and personal health, as well as spacecraft and aircrafts. We describe here principles of an automatic set of programs that begin with "FEP-Search", used to determine the beginning of a large FEP event. After a positive signal from "FEP-Search", the following programs start working: "FEP-Research/Spectrum", and then "FEP-Research/Time of Ejection", "FEP-Research /Source" and "FEP-Research/Diffusion", which online determine properties of FEP generation and propagation. On the basis of the obtained information, the next set of programs immediately start to work: "FEP-Forecasting/Spacecrafts", "FEP-Forecasting/Aircrafts", "FEP-Forecasting/Ground", which determine the expected differential and integral fluxes and total fluency for spacecraft on different orbits, aircrafts on different airlines, and on the ground, depending on altitude and cutoff rigidity. If the level of radiation hazard is expected to be dangerous for high level technology or/and personal health, the following programs will be used "FEP-Alert/Spacecrafts", "FEP-Alert/ Aircrafts", "FEP-Alert/Ground".


Introduction
It is well known that Flare Energetic Particle (FEP) events in the beginning stage are very anisotropic, especially during great events as in February 1956, July 1959, August 1972, September-October 1989, July 2000, January 2005, and many others (e.g. Dorman, 1957Dorman, , 1963aDorman, ,b, 1978Carmichael, 1962;Dorman and Miroshnichenko, 1968;Duggal, 1979;Shea and Smart, 1990;Stoker, 1994;Miroshnichenko, 2001;Bieber et al., 2002Bieber et al., , 2005. On the basis of experimental data, determination of the properties of the FEP source and parameters of propagation, i.e. to solve the inverse problem, is very difficult, and requires data from many CR stations. By using the "FEP-Search" for each CR station the start of a FEP event can be automatically determined and then the program "FEP-Research/Spectrum" determines for different moments of time, the energy spectrum of FEP out of the atmosphere over the individual CR station. As a result, we may obtain the planetary distribution of FEP intensity over the atmosphere and then, by taking into account the influence of the geomagnetic field on particles trajectories -the FEP angle distribution out of the Earth's magnetosphere. In this way, by using on a the planetary network of CR stations with online registration on a real-time scale the continued online monitoring of great ground observed FEP events (see Dorman et al., 2004;Mavromichalaki et al., 2004). This information, as well as what is obtained for many historical ground level FEP events can take place (mentioned above), will be especially useful in the near future during the development of a method of great radiation hazard forecasting on the basis of ground CR observations using the planetary neutron monitor network. Unfortunately, this procedure cannot be realized now, but we hope it, can be applied in the near future, when on a real-time scale it will provide one-minute data for the main part of the worldwide CR stations' network.
In this paper we describe the first step, practically based on two well established facts: 1) the time of particle acceleration on the Sun and injection into solar wind is very short in comparison with the time of propagation, so the source function can be considered as a delta function from time; 2) the very anisotropic distribution of FEP with the development of the event over time became near isotropic after a few scattering of energetic particles (well-known examples of February 1956, September 1989 andmany others). In the first step of forecasting of a great FEP we will use only one online detector on the ground for high energy particles and one online detector on the satellite for small energies. Therefore, we are based here on the simplest model of generation (delta function in time and in space) and on the simplest model of propagation (isotropic diffusion). In this case we have only a few parameters (time of FEP ejection, spectrum in source, and diffusion coefficient in dependence from energy and distance from the Sun) in determining all of the behavior of FEP event development in time. We obtained algorithms for determining all these parameters on the basis of CR observations on the neutron monitor (NM) and on the satellite. Using these parameters at each moment, FEP development and the expected variation of the FEP intensity in the NM counting rate and on the satellite was forcased, and compared with observations. For the first 10-15 min after the beginning of FEP, we obtained no agreement between the forecasting and the observations. But because overtime we continued to determine parameters and forecast additional new experimental data, and because the FEP distribution became more and more isotropic, the forecasting with time became better and better, and after about 30-35 min the forecast practically coincided with observations, and resulted in good forecasting for more than 50 h (by using data for only half an hour). Of course, it would have been better to obtain a good forecast earlier, maybe by using only the first 10-15 min after the FEP start, but for this we needed to use a network of NM and a much more complicated model of propagation. Now we are working hard at this problem, on the one hand, through the development of the experiment (see Mavromichalaki et al., 2004) and on the other hand, to obtain algorithms for the inverse problem of FEP propagation described by anisotropic diffusion and by the kinetic equation (Dorman and Miroshnichenko, 1968;Dorman and Kats, 1977;Bieber et al., 2002Bieber et al., , 2005Dorman et al., 2003). It is a difficult problem, but we hope that it can be solved step-by-step in near future.
The first of the programs under "FEP-Research" is the program "FEP-Research/Spectrum". We consider two variants: 1) geomagnetic quiet period (no change in cutoff rigidity), 2) disturbed period (characterized with possible changing in cutoff rigidity). We describe the method of determining the spectrum of FEP out of the atmosphere by using the so-called Dorman functions (the method was developed in Dorman, 1957Dorman, , 1969; the last review in Dorman, 2004) in the quiet period (for this we need data for at least two components with different Dorman functions) and in the disturbed period (it needs data for at least three different components).
We show that after determining the FEP spectrum at different moments of time, the time of ejection, diffusion coefficient in the interplanetary space and energy spectrum in source of FEP. By using data for four different moments of time can be determined all three unknown parameters (time of ejection, diffusion coefficient in the interplanetary space and energy spectrum in source of FEP) can be determined. We describe in detail the algorithms of the programs "FEP-Research/Time of Ejection", "FEP-Research /Source" and "FEP-Research/Diffusion". We show how these programs worked on the example of the historical great FEP event in September 1989. On the basis of these programs the time of ejection, diffusion coefficient in the interplanetary space and energy spectrum in source of FEP can be determined online. To extend the obtained information to the region of very small energies, we use simultaneously with the NM data the available satellite one-minute data. We show how, on the basis of these results, forecasting of expected radiation hazard for computers, electronics, solar batteries, technology in space on different distances from the Sun and on different helio-latitudes can be achieved. We show that the same can be done for satellites on different orbits in the magnetosphere by taking into account the change in cutoff rigidities along the orbits (for personal health, solar batteries, computers, electronics, technology). Using the method of the Dorman functions for different altitudes in the atmosphere, we describe the principles of radiation hazard forecasting online for airplanes, depending on altitudes and the cut-off rigidities, and the value of shielding. On-line will be made Forecasting of radiation hazard on the ground for personal health and technology, depending on the cut-off rigidity and atmospheric pressure, will be made online. If for some cases the calculated radiation hazard is expected to be dangerous, special alerts will be sent online.

Data from the past and classification of space weather radiation hazard (NOAA classification and its modernization)
NOAA Space Weather Scale establishes 5 gradations of FEP events, what are called Solar Radiation Storms: from S5 (the highest level of radiation, corresponding to the flux of solar protons with energy >10 MeV about 10 5 protons cm −2 s −1 ) up to S1 (the lowest level, the flux about 10 protons cm −2 s −1 for protons with energy >10 MeV). According to Dorman (2004), in the first, for satellite damage and influence on personal health and technology, on communications by HF radio waves, is the total fluency of FEP during the event is more important than the protons flux, what is used now in NOAA Space Weather Scale; in the second, the level S5 is not maximal. As it was shown by McCracken et al. (2001), the dependence of event probability on fluency can be prolonged at least up to F =2×10 10 protons cm −2 for protons with E k ≥30 MeV, what was observed in the FEP of September 1869, according to nitrate contents data in polar ice. This type of great dangerous event is very rare (about one in a few hundred years). According to McCracken et al. (2001), it is not excluded that in principle very great FEP events can occur with a fluency even 10 times bigger (one in a few thousand years). So, we intend to correct the very important classification, developed by NOAA, in two directions: to use fluency F of FEP during the event instead of flux I, and to extend the levels of radiation hazard. The existing radiation level S3 in this case will correspond at E k ≥30 MeV to F ≈10 7 protons cm −2 (about 1 event per year), S4 -to F ≈10 8 protons cm −2 (1 event per 3-4 years), and S5 -to F ≈10 9 protons cm −2 (1 event per 20-50 years). The supposed additional radiation level S6 at E k ≥30 MeV will correspond to F ≈10 10 protons cm −2 (1 event in several hundred years), and S7 -to F ≈10 11 protons cm −2 (1 event in several thousand years). Let us note that the Radiation Storms S5 -S7 are especially dangerous not only for electronics, navigation, and communication systems, for spacecraft, aircrafts and on the ground, but also for personal health (up to a lethal dose). The forecasting of these rare but very dangerous FEP events is especially important to avoid planetary catastrophic damage. Let us note that the above mentioned expected probabilities of dangerous FEP events are averaged over the solar cycle. In reality, these probabilities are much higher in periods of high solar activity than in periods of low solar activity Pustil'nik, 1995, 1999).

The method of automatically searching for the start of a great FEP event
Let us consider the problem of automatically searching for the start of a great FEP event. The determination of increasing flux is made by comparison with intensity averaged from 120 to 61 min before the present Z-th one-minute data. For each Z minute data, start the program "FEP-Search-1 min" which, for both independent channels A and B, and for each Z-th minute, determines the values where I Ak and I Bk are one-minute total intensities in the sections of neutron super-monitor A and B. If simultaneously the program "FEP-Search-1 min" repeats the calculation for the next Z+1-th minute and if Eq. (3) is satisfied again, the onset of a great FEP is established and the program "FEP-Research/Spectrum" starts.

The probability of false alarms
Given the probability function (2.5) =0.9876, the probability of an accidental increase with amplitude more than 2.5σ in one channel will be (1− (2.5)) 2=0.0062 min −1 , which means one in 161.3 min (in one day we expect 8.93 accidental increases in one channel). The probability of accidental increases simultaneously in both channels will be (1− (2.5)) 2 2 =3.845×10 −5 min −1 which means one in 26007 min ≈18 days. The probability that the increases of 2.5σ will be accidental in both channels in two successive minutes is equal to (1− (2.5)) 2 4 =1.478×10 −9 min −1 , which means one in 6.76×10 8 min., ≈1286 years. If this false alarm (one in about 1300 years) is sent, it is not dangerous, because the first alarm is preliminary and can be cancelled if in the third successive minute there is no increase in both channels larger than 2.5σ (it is not excluded that in the third minute there will also be an accidental increase, but the probability of this false alarm is negligible: (1− (2.5)) 2 6 =5.685×10 −14 min −1 , which means one in 3.34×10 7 years). Let us note that the false alarm can also be sent in the case of a solar neutron event (which really is not dangerous for electronics in spacecraft or for astronauts' health), but this event usually is very short (only a few minutes) and this alarm will be automatically canceled in the successive minute after the end of a solar neutron event.

The probability of missed triggers
The probability of missed triggers depends very strongly on the amplitude of the increase. Let us suppose, for example, that we have a real increase of 7σ (that for ESO corresponds to an increase of about 9.8 %). The trigger will be missed in either of both channels and in either of both successive minutes if, as a result of statistical fluctuations, the increase in intensity is less than 2.5σ . For this the statistical fluctuation must be negative with an amplitude of more than 4.5σ . The probability of this negative fluctuation in one channel in one minute is equal to (1− (4.5)) 2=3.39×10 −6 min −1 , and the probability of a missed trigger for two successive minutes of observation simultaneously in two channels is 4 times larger: 1.36×10 −5 . It means that a missed trigger is expected only to be one per about 70 000 events.

Analytical approximation for coupling functions
Based on the latitude survey data of Alexanyan et al. (1985), Moraal et al. (1989), Clem and Dorman (2000), , the polar normalized coupling functions connected the observed spectrum inside the atmosphere with the spectrum outside of the atmosphere (introduced in Dorman, 1957) for a total counting rate, and different multiplicities m can be approximated by the so-called Dorman function (Dorman, 1969): where m=tot, 1, 2, 3, 4, 5, 6, 7, ≥8. Coupling functions for muon telescopes with different zenith angles θ can be approximated by the same type of functions determined only by two parameters, a m (θ ) and k m (θ ). Let us note that the Dorman polar functions are normalized: at any value of a m and k m . The Dorman function for a point with cutoff rigidity R c , will be The first appoximations of the FEP energy spectrum In the first approximation the spectrum of the primary variation of the FEP event (out of the atmosphere) can be described by the function where is the differential spectrum of galactic cosmic rays before the FEP event and D (R, t) is the spectrum at a later time t. Approximation (6) can be used for describing a limited interval of energies in the sensitivity range detected by the various components.

Online determining of the FEP spectrum from data of a single observatory in the case of a magnetically quiet period
In this case the observed variation δI m (R c , t) ≡ I m (R c , t) I mo (R c ) of some component m can be described in the first approximation according to Dorman (1957Dorman ( , 1975 by function F m (R c , γ ): where m= tot, 1, 2, 3, 4, 5, 6, 7, ≥8 for neutron monitor data (but can also be denote the data obtained by muon telescopes at different zenith angles and data from satellites), and is a known function. Let us compare data for two components, m and n. According to Eq. (7) we obtain δI m (R c , t) δI n (R c , t) = mn (R c , γ ) , where mn (R c , γ ) = F m (R c , γ ) F n (R c , γ ) are calculated by using Eq. (8). Comparison of experimental results with function mn (R c , γ ), according to Eq. (9), gives the value of γ (t), and then from Eq. (7) the value of the parameter b (t).

Determining of the FEP spectrum for a magnetically distrubed period
For magnetically disturbed periods the observed CR variation instead of Eq. (6) will be described by where R c (t) is the change in cutoff rigidity due to the change of the Earth's magnetic field, and W k (R c , R c ) is determined according to Eq. (5) at R=R c : Now we have unknown variables γ (t), b (t), R c (t), and for their determination we need data from at least 3 different components, k=l, m, n in Eq. (11). In accordance with the spectrographic method (Dorman, 1975), let us introduce the function Then from equation lmn (R c , γ ) = W l (R c , R c ) δI m (R c , t) − W m (R c , R c ) δI l (R c , t) W m (R c , R c ) δI n (R c , t) − W n (R c , R c ) δI m (R c , t) .
the value of γ (t) can be determined. Using the determined value of γ (t), for each time t, we determine R c (t) = F l (R c , γ (t)) δI m (R c , t) − F m (R c , γ (t)) δI l (R c , t) F m (R c , γ (t)) δI n (R c , t) − F n (R c , γ (t)) δI m (R c , t) , So, in magnetically disturbed periods, the observed FEP increase for different components also allows for the determination of parameters γ (t) and b (t) ,for the FEP spectrum beyond the Earth's magnetosphere, and R c (t), giving information on the magnetospheric ring currents.

Online simultaneously determination of time of ejection, diffusion coefficient and FEP spectrum in source
As we described in Sect. 1, we suppose that the time variation of FEP flux and energy spectrum can be described, in the first approximation, by the solution of isotropic diffusion from the pointing instantaneous source described by function Q (R, r, t) =N o (R) δ (r) δ (t) . In this case the expected FEP rigidity spectrum on the distance r from the Sun at the time t after ejection will be where N o (R) is the rigidity spectrum of the total number of FEP in the source, t is the time relative to the time of ejection and K (R) is the diffusion coefficient in the interplanetary space in the period of FEP event. Let us suppose that the time of ejection T e , diffusion coefficient K (R), and the source spectrum N o (R) are unknown. In this case for determining online simultaneously the time of ejection T e , diffusion coefficient K (R) and FEP spectrum in source N o (R), we need information on the observed FEP spectrum at least in three moments of time, T 1 , T 2 and T 3 (all times T are in UT scale). In this case we will have for the times after FEP ejection from the Sun into solar wind: where T 2 −T 1 and T 3 −T 1 are known values and x is an unknown value which we need to determine first. From the three equations of type Eq. (17), but for times t 1 , t 2 and t 3 , we obtain After dividing Eq. (19) by Eq. (20), we obtain Equation (21) can be solved by the iteration method: as the first approximation, we can use x 1 =T 1 −T e ≈500 s which is a minimum time of relativistic particle propagation from the Sun to the Earth's orbit. Then by Eq. (22), we determine (x 1 ) and by Eq. (21) determine the second approximation x 2 , and so on. After solving Eq. (21) and determining the time of ejection, we compute very easyly the diffusion coefficient from Eq. (19) or Eq. (20): After determining the time of ejection and diffusion coefficient, it is very easy to determine the FEP spectrum in source:

The model testing by controlling the diffusion coefficient
In order to check above model of FEP propagation in the interplanetary space, we determine the first values of K (R). These calculations we remade according to the procedure described above with the assumption that K (R) does not depend on the distance to the Sun. The results are shown in Fig. 1. It can be seen that in the beginning of the event, the obtained results are not stable due to relatively big statistical errors. A few minutes after the FEP beginning, the amplitude of the CR intensity increase became many times larger than σ , and we see a systematical increase in the diffusion coefficient with time: in reality, it reflects the increasing of K (R) with the distance to the Sun.

The inverse problem for the case when the diffusion coefficient depends on the distance to the Sun
Let us suppose, according to Parker (1963), that the diffusion coefficient K (R, r) =K 1 (R) × r r 1 β . In this case If we know n 1 , n 2 , n 3 at moments of time t 1 , t 2 , t 3 , the final solutions for β, K 1 (R) , and N o (R) will be   × ln n 1 n 2 − In the last Eq. (28) index k=1, 2 or 3.

Simulation of FEP forecasting by using only neutron monitor data
By using the first few minutes of the FEP event NM data we can determine by Eqs. (26)-(28) the effective parameters β, K 1 (R), and N o (R), corresponding to a rigidity of about 7-10 GV, and then by Eq. (25) we determine the forecasting curve of the expected FEP flux behavior for total neutron intensity. We compare this curve with the time variation of the observed total neutron intensity. In reality, we use data for more than three moments of time by fitting the obtained results in comparison with experimental data to reach the minimal residual (see Fig. 2, which contains 8 figures for time moments t=110 min up to t=220 min after 10.00 UT on 29 September 1989). From Fig. 2 it can be seen that using only the first few minutes of NM data (t=110 min) is not enough: the obtained curve forecasts are too low in intensity. For t=115 min the forecast shows a little bigger intensity, but also not enough. Only for t=120 min (15 min after beginning) and later, up to t=140 min (35 min after beginning) can we obtain a stable forecast with good agreement with the observed CR intensity (with accuracy about ± 10 %). In Fig. 3 values of parameter K 1 (R) are shown. From Fig. 3 it can be seen that at the very beginning of the event (the first point), the result is unstable: in this period the amplitude of increase is relatively small, so that the relative accuracy is too low. Let us note that for the very beginning part of the event, the diffusion model can be hardly applied (more natural is the application of the kinetic model of FEP propagation). After the first point we have a stable result with an accuracy ±20%, which is about the same we obtained for parameter β(the first point is unusually big, but after the computation β became stable with an average value β∼0.6). Therefore, we came to the conclusion that the model, described by Eqs. (25)-(28), reflects well the FEP propagation in the interplanetary space.
14 Simulation of FEP forecasting by online using both neutron monitor and satellite data All the described above results, based on NM online data, reflect the situation in FEP behavior in the high energy (more than few GeV) region. For extrapolation of these results to the low energy interval (dangerous for space-probes and satellites), we use satellite online data available through the Internet. The problem is how to extrapolate the FEP energy spectrum from high NM energies to very low energies detected in the GOES satellite. The main idea of this extrapolation is as follows: the source function, time of ejection and diffusion coefficient in both energy ranges are the same. The source function relative to time and space is the δ−function, and relative to energy is the power function, with an energy-dependent index γ (γ =γ o + ln E k E ko with maximum at E k max = E ko exp (−γ o ): In Fig. 4 the results are shown based on the NM and satellite data of forecasting of the expected FEP fluxes also in small energy intervals and the comparison with observation satellite data.
From Fig. 4 it can be seen that by using online data from ground NM in the high energy range and from satellites in the low energy range during the first 30-40 min after the start of the FEP event, it is possible to predict the expected FEP integral fluxes for different energies up to a few days ahead.

Forecasting of expected total FEP fluency
In Fig. 5 we show the results of calculations for the expected total (event-integrated) FEP fluency for E k ≥E o =0.1 GeV; the stable predictions are reached at about 40-50 min after the start.

Alerts in cases where the fluxes and fluency are expected to be dangerous
If the predicted fluxes and fluency are expected to be dangerous, preliminary "FEP-Alert-1/Space", "FEP-Alert-1/Magnetosphere", "FEP-Alert-1/Atmosphere" will be sent in the few minutes after the beginning of the event. As more data become available, better predictions of the expected fluxes will be made and more definitive Alert-2, Alert-3 and so on will automatically be issued. Alerts will give information on the expected time and level of dangerous situation for different objects in space, in the magnetosphere, in the atmosphere, on different altitudes and at different cut-off rigidities.

Conclusion
We show that by using online data from ground NM in the high energy range and from satellites in the low energy range during the first 30-40 min after the start of the FEP event, it is possible to predict the expected FEP integral fluxes for different energies up to a few days ahead. The total (eventintegrated) fluency of the event, and the expected radiation hazards can also be estimated for 30-40 min and the corresponding Alerts to experts operat different objects in space, in the magnetosphere, and in the atmosphere on different altitudes and at different cutoff rigidities, can be sent automatically, and they must decide what to do operationally (for example, for space-probes in space and satellites in the magnetosphere whether to switch-off the electric power for a few   hours to save the memory of the computers and high level electronics; for jets to decrease their altitudes from 10 km to 4-5 km, to protect crew and passengers from great radiation hazard, and so on).