Observation of gamma ray bursts at ground level under the thunderclouds

We observed three $\gamma$-ray bursts related to thunderclouds in winter using the prototype of anti-neutrino detector PANDA made of 360-kg plastic scintillator deployed at Ohi Power Station at the coastal area of the Japan Sea. The maximum rate of the events which deposited the energy higher than $3\,$MeV was $(5.5 \pm 0.1) \times 10^2 {\rm /s}$. Monte Carlo simulation showed that electrons with approximately monochromatic energy falling downwards from altitudes of order $100\,$m roughly produced the observed total energy spectra of the bursts. It is supposed that secondary cosmic-ray electrons, which act as seed, were accelerated in electric field of thunderclouds and multiplied by relativistic runaway electron avalanche. We actually found that the $\gamma$-rays of the bursts entered into the detector from the direction close to the zenith. The direction stayed constant during the burst within the detector resolution. In addition, taking advantage of the delayed coincidence detection of the detector, we found neutron events in one of the bursts at the maximum rate of $\sim 14\pm5\,{\rm /s}$.


Introduction
In the early 1920's, C.T.R. Wilson suggested that strong electric fields in thunderclouds might accelerate free electrons present in the atmosphere to high energies [1]. Since then, radiation associated with thunderstorms attracted the interest as natural particle-acceleration process and many experiments have been attempted to detect these radiations in various environments.
For instance, bursts of γ-rays were observed on orbiting satellites with energy up to tens of MeV and with duration of less than 1 ms. They are called Terrestrial Gamma-ray Flashes(TGF's) [2].
Recently, Dwyer et al. [3] reported unexpected observation of positron bursts, lasting about 0.2 s, by an airborne detector when the aeroplane flew into a thundercloud.
On the other hand, γ-ray flux enhancements of longer duration of order 100 s were reported in limited environments like high mountains [4,5,6,7,8,9,10,11,12,13] and sea level locations in the coastal area of the Japan Sea. They are also called Thunderstorm Ground Enhancements(TGE's) [13]. Japanese groups found that radiation monitoring posts or dedicated scintillation counters in and near nuclear power plants signaled an increase of γ-ray dose which seemed to originate from low altitude winter thunderclouds [14,15,16,17,18,19]. Especially, Torii et al. [20] found that area of γ-ray flux enhancements was moving as the associated thundercloud passed across the observation site.
Gurevich et al. [21,22] developed the runaway electron model to explain the electron acceleration in the electric field of the thunderclouds. The stopping power of air for electrons decreases with increasing electron energy and goes up again by relativistic effects. Therefore, electric field in the thundercloud may accelerate electrons if the the electric force is larger than the minimum stopping power and the electron energy is in the region where the electric force exceeds the stopping power. Such electrons are called runaway electrons. By generating knock-on electrons successively, the runaway electrons can cause an avalanche multiplication process called relativistic runaway electron avalanche (RREA).
Numerical simulations [23,24,25] with models of thundercloud electric field suggested that the avalanche can be produced continuously if energetic seed electrons are provided, for example, by cosmic ray secondaries. A significant flux of relativistic runaway electrons in the lower parts of thunderclouds is capable of producing intensive bremsstrahlung which can reach the Earth's surface or the mountain top to account for the observed flux enhancement.
Recently, the neutron bursts associated with thunderstorms were also observed in various experiments [26,12,27,28,29]. The generation of neutrons is most probably by photoproduction by γ-rays with air nuclei as the detected γ ray spectrum extends above the photonuclear reaction threshold for nitrogen (∼ 10.5MeV) [30]. It may have a significant effect on 14 C dating [31,32] through the neutron capture reaction 14 N(n, p) 14 C.
Our research group have developed prototypes of a reactor neutrino detector "PANDA", which stands for Plastic Anti-Neutrino Detector Array [33,34]. We have originally targeted PANDA at presenting the feasibility of reactor monitoring using neutrinos with a tonne-size detector. γ-rays and neutrons can also be detected by PANDA by Compton scattering and the delayed coincidence of proton recoil and neutron captures. We installed the PANDA detector outside of the reactor building of Ohi power station, which stands near the Japan Sea in Fukui, and tried to watch the reactor operating status via detecting and analyzing the anti-electron neutrinos produced in the reactor core.
We accidentally found that there were intensive increases of γ-ray flux correlated with the winter thunder-storm activity during the measurement. In this paper, we report the investigated properties of these burst events taking advantages of the unprecedented features of the detector including high statistics, good energy response, direction sensitivity and neutron identification.
The module was made of a plastic scintillator bar(10cm × 10cm × 100cm) with effective mass of about 10kg wrapped with aluminized Mylar films and gadolinium (Gd) coated Mylar films (4.9 mg of Gd per cm 2 ). Each bar was connected to acrylic light guides and photomultipliers on both ends (Figure 1).
The light intensity ratio seen by each PMT pair allows one to estimate the position of the hit along the module [33]. Using the position of the hit Each PMT signal was divided into two: about 15% of the original charge was sent to CAEN V792 multi-event Charge-to-Digital-Converters (QDCs) and the other 85% was passed to CAEN V895 leading edge discriminators.
The discriminator outputs were sent to CAEN V1495 general purpose VME board, which has customizable FPGA unit (Altera Cyclone EP1C20). The logic counted the number of pairs of fired PMTs seeing the same scintillator. Whenever the number of the pairs was greater than or equal to two, the logic generated the gate pulses of 400 ns duration for the QDCs.
The timing of the gate pulses and busy signals from the QDCs were recorded by the same FPGA. We used these time stamps to select neutrino events by delayed coincidence method offline.
The PANDA36 detector was loaded on and transported by a 2-tonne dry van. The detector was deployed beside the Unit 2 of Ohi Power Station (35 • 32 ′ 32 ′′ N, 135 • 39 ′ 14 ′′ E and about 10 m above the sea level) of Kansai Electric Power Co., Inc on November 18th, 2011. We continued the measurement for 62 days.
Energy calibrations were carried out before the deployment using the Compton edge of 60 Co γ-rays. Time drifting of gains of each PMT and QDC was corrected using the peak of through-going cosmic muons in the spectrum of the events.

Event-rate increase
In the data acquired by the PANDA36 detector through the neutrino detection experiment, we found unexpected increases of event rate. The trigger rate got twice or higher for a few minutes for the events with total energy deposit larger than 3 MeV independently of the reactor operation.
Temporal variation of the event rate are shown in Figure 2. Burst duration is defined as an interval whose event rates are 5σ greater than average event rate. The date and time, duration and peak event rate are summarized in Table 1. We hereafter call each burst by the name defined in the table.
We searched for the association between thunderclouds and the radiation bursts via "Kaminari Nowcast" data provided by Japan Meteorological Agency. Kaminari Nowcast provides the thunder activity information analyzed from lightning discharges detected by Lightning Detection Network System (LIDEN) and radar observations for every ∼ 1 km grid [35].
Kaminari Nowcast data show that there were more than one grid which have level 2 thunder activity of five levels around the experimental site for all three bursts at almost the same time that the radiation bursts have been observed.
Conversly, we found 22 time-consecutive data sets of level 2 or higher in 20 × 20 grids around the detector. Nevertheless, radiation bursts are observed only three times. It is not strange since the γ-ray emitting region is conceivably much smaller than the above area.

The energy and the height of the source electrons
We performed Monte Carlo simulations to estimate the energy and the height of the source electrons using Geant4 toolkit [36], assuming the γ-rays are caused by high energy electrons emitted downwards from thunderclouds.
At first, we simulated the numbers and spectra of bremsstrahlung photons at the ground level initiated by mono-energetic electrons projected vertically downward from the sky. The height and energy of the source electrons were chosen to be combinations of 100 m, 500 m, 900 m, 1300 m, and 17 MeV, 23 MeV, 29 MeV, 35 MeV. As a generic spectral shape of the source electrons, we estimated it as a sum of these four components. Lower-energy electron components were ignored for the efficiency of analysis. It does not necessarily exclude the theoretically expected exponential electron energy distribution [37]. Lower-energy components of electrons, if exist, are less efficient in producing bremsstrahlung γ-rays on the ground, and therefore flux of those lower-energy electrons could only be determined with poor accuracy.
Then we calculated the detector response to those γ-rays. Next, we calculated the weighted sum of the spectra of the detector responses to the simulated bremsstrahlung γ-rays from four energies of electrons projected from each height.
The summed spectra were compared with the observed spectra of the bursts after background subtraction. The background data were taken in two ten-minute intervals starting thirteen minutes before and three minutes after the burst periods. It is seen that the observed spectra extend up to 15 MeV or more. We fitted the weights to minimize the χ 2 for each projection height. The analyses were done for all the three bursts and the fit results with the heights and the weights which minimize χ 2 are shown in Figure 3. The source electron spectra with the weights at the projection heights are shown in Figure 4.
They indicate that 17 MeV electrons, as a somewhat arbitrary choice, produce a spectrum qualitatively similar to the data. Peak flux of the source electrons of each burst is (1.9 ± 0.1) × 10 5 m −2 s −1 , (2.0 ± 0.1) × 10 5 m −2 s −1 and (8.8 ± 0.5) × 10 4 m −2 s −1 , for burst-20111225, burst-20120102 and burst-20120105, respectively. It should be noted that the estimated source electron fluxes are the lower limits since we ignored the electron components lower than 17MeV. Estimated peak flux of the source electrons is also plotted as a function of the assumed height in Figure 5 to see the dependence of the flux on height.

Arrival direction of γ-rays
The arrival direction of a γ-ray can be investigated, if the γ-ray is Comptonscattered by an electron in a plastic scintillator module and then the scattered  The arrival direction lies along the cone called the Compton cone with its axis being the line connecting two interaction points and its half opening angle α which satisfies the relation, where E γ is the energy of the incident photon and E ′ γ is the energy of the scattered photon.
For the analysis, we chose the events with the total energy deposit in the range 5 MeV ≤ E total ≤ 12 MeV to remove the ordinary environmental γ-rays and to take those γ-rays which deposited enough energy on each module to ensure sufficient position resolution. We assumed E total to be the incident γ-ray energy, E γ . We also assumed that E γ − E ′ γ and E ′ γ correspond to E 1st and E 2nd , the highest and the second highest energy deposit of all the modules, or vice versa. We simply chose the combination which corresponds to the gamma ray coming from the upper hemisphere.
The first assumption that E total = E γ is often not true. To reduce the effect of error in E total , we introduced the selection criterion, We can cut the energy region near the Compton edge by this selection , where the effect of E total error to cos α is large. The selection can effectively cut the event with two Compton scatterings in one module, too. In addition, the positions of E 1st and E 2nd were required to be interleaved with more than two modules along the stacking direction. This selection reduces the effect of the periodical structure of the detector which may discretize the result of the scattering angle calculation.
After above preparations, the Compton cone and the corresponding circle on the unit sphere centered at the detector center was calculated for each event which satisfied the selection criteria. Then each event is counted in the predefined grids in (cos θ, φ) space with a weight of the fraction of the circumference which overlaps the grid. The polar angle θ and azimuth φ are defined with respect to the axis along the length of the modules so that the zenith is (0,π/2). Then we normalized the weighted number of events in each grid by the live time of the burst so that it represents the arrival rate of the selected event from the direction.
Monte Carlo simulations were made for γ-rays isotropically incident on the detector and then the same arrival-direction analysis was made as for the data. We made the following calculation to get rid of a possible bias due to the detector response, Here i represents the grid number, M i,burst , M i,bg are the arrival rate of each grid in the burst period and the background period, respectively. M i,sim is the arrival rate calculated by the Monte Carlo simulation, which is normalized so as to get the average value of all the grids to be 1. Maps of M i , M i,burst −M i,bg and M i,sim for burst-20120105 are shown in Figure 6 as typical examples. We found the arrival direction was close to the zenith as shown in Figure 6 and stayed constant during the periods of all three bursts within the rate of d cos θ dt ≤ (0.2 ± 0.5)/30s and dφ dt ≤ (0.4 ± 0.5)rad/30s.

Neutron event-rate increase
PANDA36 detector is capable of detecting fast neutrons by delayed coincidence method. A neutron entering the detector interacts with a proton in the plastic scintillator. A neutron transfers a part of its energy to the recoil-proton by an elastic scattering. It is referred to as the prompt event. Then, the scattered neutron loses its energy by subsequent multiple scatterings, and after O(10)µs it is captured by a gadolinium (Gd) nucleus in the Gd coated Mylar films wrapped around the scintillator. The neutron capture on Gd results in a gamma cascade emission with total energy of 7.9 MeV for 157 Gd and 8.5 MeV for 155 Gd. It is referred to as the delayed event.
We selected two kinds of delayed coincidence events, correlated events with the delay time window of 8 − 150µs and accidental events with the time window of 1008 − 1150µs. The total energy is required to be 1.5 MeV ≤ E total ≤ 10.0 MeV for both the selections. We calculated the number of neutron events by subtracting the accidental event rate from the correlated event rate.
As a result, we found event-rate increase which is synchronized with the γ-ray burst with maximum rate of 14 ± 5 /s in burst-20120105 ( Figure 7).

Discussion
Our observational results are mostly consistent with a model of the long duration γ-ray bursts from thunderclouds [25]. In the model, seed electrons are provided to the thunderclouds continuously mainly by secondary cosmic rays, and they are multiplied by RREA process with a help of electric field between the negative charge of the lower part of a thundercloud and the positive pocket charge region just below the negative charge [38]. Then, those high energy electrons make elctromagnetic shower in the atmosphere resulting in numerous γ-rays on the ground. The downward vertical directions of γ-ray bursts also reflects the direction of the electric field.
The estimated almost monochromatic energies of the source electrons roughly reproduce the shape of the observed energy spectra. However, the present analysis has little power to examine the theoretically favored model [37] in which the spectrum becomes exponential shape independent of the electric field or the air density.
The duration of the bursts may depend on the movement or the development stage of the clouds. For example, the duration of the burst-20120105, ∼ 180 s can be explained by a thundercloud which have typical velocity of 50 km/h and diameter of 2.5 km, which is somewhat smaller than typical echo size of thunderclouds (which is 4-6 km [38]).
In addition, due to the low temperature, the altitude of the thunderclouds at midwinter is low in the coastal area of the Japan Sea [38], which explains the result that the observed energy spectra implies the low altitude electron source. Similar winter thunderclouds are relatively rare, but are also observed in the west coast of Norway and toward the east coast from the Great Lakes of North America.
RREA and resulting electromagnetic shower may exist also in more common high altitude thunderclouds in summer. However, γ-rays and electrons of shower may totally absorbed by the thick air under the clouds. It is also consistent with the observations that long duration bursts are reported in summer on high mountains, where the relative altitude of the clouds are very low.

Conclusion
We observed three γ-ray bursts in winter with the PANDA detector made of 360-kg plastic scintillator at Ohi Power Station which stands on the coastal area of the Japan Sea. Table. 1 summarizes the features of the observed bursts. The maximum rate of the events with E total ≥ 3 MeV was (5.5 ± 0.1) × 10 2 /s of burst-20120105.
We found that for all the bursts periods, there were active thunderclouds near the detector.
In addition, we found that γ-rays of the bursts entered into the detector from the direction close to the zenith. The arrival direction stayed constant during the burst.
These results indicated that the bursts originated in thunderclouds. Monte Carlo simulation showed that the observed E total spectra of the bursts are reproduced by the bremsstrahlung γ-rays by electrons with more or less monochromatic energy shown in Figure 4 from low altitudes.
The arrival direction of the γ-rays and the estimated energy of the source electrons of over 10 MeV can be described by relativistic runaway electron avalanche (RREA). Namely, the secondary cosmic-ray electrons, which act as seed, were accelerated and amplified in electric field of thunderclouds by avalanche multiplication.
Additionally, neutrons were detected in one of the bursts at the maximum rate of ∼ 14 ± 5 /s with high confidence by the delayed coincidence method. There is probability the event rate increase includes neutrons emitted by the photonuclear reactions on the nitrogen atoms in the air. It could be due to the large flux of bremsstrahlung γ-ray with energy greater than the photonuclear reaction threshold for nitrogen. The observation of fast neutrons on the ground implies that more neutrons are produced in the air between a thundercloud and the ground and even in the cloud itself. However, only small fraction of those neutrons reach the ground because of the short absorption length in the air [39].