Catalog of Gamma-ray Glows during Four Winter Seasons in Japan

In 2015 the Gamma-Ray Observation of Winter Thunderstorms (GROWTH) collaboration launched a mapping observation campaign for high-energy atmospheric phenomena related to thunderstorms and lightning discharges. This campaign has developed a detection network of gamma rays with up to 10 radiation monitors installed in Kanazawa and Komatsu cities, Ishikawa Prefecture, Japan, where low-charge-center winter thunderstorms frequently occur. During four winter seasons from October 2016 to April 2020, in total 70 gamma-ray glows, minute-lasting bursts of gamma rays originating from thunderclouds, were detected. Their average duration is 58.9 sec. Among the detected events, 77% were observed in nighttime. The gamma-ray glows can be classified into temporally-symmetric, temporally-asymmetric, and lightning-terminated types based on their count-rate histories. An averaged energy spectrum of the gamma-ray glows is well fitted with a power-law function with an exponential cutoff, whose photon index, cutoff energy, and flux are $0.613\pm0.009$, $4.68\pm0.04$ MeV, and $(1.013\pm0.003)\times10^{-5}$ erg cm$^{-2}$ s$^{-1}$ (0.2-20.0MeV), respectively. The present paper provides the first catalog of gamma-ray glows and their statistical analysis detected during winter thunderstorms in the Kanazawa and Komatsu areas.

tal Center in Armenia [4,17,18], Yangbajing in Tibet, China [8], Mt. Norikura in Japan [19], Tien Shan in Kazakhstan [20], as well as the weather station at Mt. Fuji in Japan [13]. In particular the Aragats Space Environmental Center have observed the largest number of TGEs in the world [9].
The other opportunities to detect gamma-ray glows at ground level are winter thunderstorms in Japan. When dried and cold seasonal wind comes from a Siberian cold air mass over the Tsushima warm current, an ac-tive ascending air current, then a thunderstorm can be formed in the north-west coast of the Japan main island called the "Hokuriku" region [26,27]. One of the important characteristics of these winter thunderstorms distinguished from summer ones is low charge center. Precipitation particles that contribute to charging (i.e. ice crystals, graupels) are produced at an altitude whose temperature ranges from −10 • C to 0 • C or higher [28,29]. It typically corresponds to ∼3 km for summer thunderstorms, but less than 1 km for winter ones [30]. There- fore, gamma rays emitted from winter thunderclouds can be detected even at sea-level. High-energy emission form winter thunderstorms in Japan has been discovered unexpectedly and investigated by radiations monitors in nuclear power stations [2,31]. Following these findings, we have performed a ground-based radiation measurement called the Gamma-Ray Observation of Winter Thunderclouds (GROWTH) experiment in the Hokuriku region since 2006 [3,6,22,[32][33][34].
Gamma-ray glows are thought to originate from strong electric fields formed within thunderclouds, not directly from lightning discharges. On the other hand, they are often quenched by nearby discharges since lightning currents could discharge the corresponding electric fields [1,6,7,18,32]. More rare cases of the association between gamma-ray glows and lightning discharges are ter-restrial gamma-ray flashes (TGFs). TGFs are instantaneous gamma-ray emission coincident with lightning discharges [35][36][37][38]. Their duration is typically a few hundreds of microseconds [38,39], and their energy spectrum also extends up to >20 MeV [36,40,41]. TGFs have been detected mainly by space-borne detectors such as CGRO [35], RHESSI [36], AGILE [37], Fermi [42], and ASIM [43]. While the TGFs observed from space are upward-going, downward TGFs have been recently detected by ground-based experiments [44][45][46][47]. In addition, a few cases of glow terminations coincided with downward TGFs [24,48]. The connection between gammaray glows and downward TGFs is one of the intriguing questions in this field.
One of the widely-known theories of electron acceleration in the atmosphere is the relativistic runaway electron avalanches (RREA) [49]. When the strength of electric fields exceeds the RREA threshold (0.284 MV m −1 at 1 atm [50,51]), seed electrons of hundreds of keV or several MeV can be accelerated up to tens of MeV, and produce secondary electrons that are also accelerated. Therefore the number of electrons can be multiplied exponentially in such a high electric field. The accelerated electrons can emit bremsstrahlung photons as they interact with ambient atmospheric atoms. The energetic seed electrons are thought to be provided by cosmic rays, for example. The multiplication processes can be further boosted by the relativistic feedback processes [51,52], and some bright glows are though to require the feedback mechanism in fact [7,24]. Even when the electric fields are less than the RREA threshold, electrons of cosmic-ray origin can gain energies from the field while electron multiplication hardly takes place. The excess in energies from the normal cosmic-ray spectra can also be detected as fainter gamma-ray glows than those produced by RREA. This mechanism is called modification of the energy spectrum (MOS) [17,53,54], and is the other important process to explain gamma-ray glows [12,55].
In 2015, we launched a mapping observation campaign with multiple radiation monitors in the Hokuriku region [22]. After the launch of the campaign, we have reported the discovery of photonuclear reactions triggered by a downward TGF [34], an observation of a glow termination with radio-frequency lightning mapping and electric field measurements [6], a simultaneous detection of a glow termination and a downward TGF [24], and so on. Besides, more and more gamma-ray glows have been detected as the number of detectors increases. The present paper provides the first catalog of gamma-ray glows observed by the GROWTH collaboration during four winter seasons in Japan. Statistical investigations of occurence rates, temporary and spectral features of gamma rays are presented. Throughout this paper, errors are at a confidence level of 1 standard deviation unless otherwise noted.

II. OBSERVATION
Since 2006 we have operated an observation site at the Kashiwazaki-Kariwa nuclear power station of Tokyo Electric Power Company Holdings in Niigata Prefecture. In addition to the Kashiwazaki site, we established Kanazawa and Komatsu observation areas in Ishikawa Prefecture for the mapping observation campaign after 2015. Their locations are presented in Figure 1. In the present paper we focus on the results obtained at the Kanazawa and Komatsu areas because very few gammaray glows have been detected at the Kashiwazaki site   Table I. Four sites are located at high schools, two at universities, two at public facilities, and the rest at a private company. Limited by the number of radiation detectors, up to 10 observation sites were simultaneously in operation (in FY2018 and FY2019). Figure 2 shows the time histories of the number of detectors in operation during four winter seasons. In each winter season, ob-servations typically start in October or November, then end in March or April. Therefore, the average operation duration is 4-5 months in each year.
Our radiation monitors consist of two types of scintillation crystals; a 25.0 × 8.0 × 2.5 cm 3 bismuth germanate (Bi 4 Ge 3 O 12 ) or a 30.0 × 5.0 × 5.0 cm 3 cesium iodide (CsI) crystal. While the BGO crystal is coupled with two Hamamatsu R1924A photomultipliers, the CsI crystal is with a Hamamatsu R6231 photomultiplier. Signals from the photomultipliers are amplified via a preamplifier and a shaping amplifier, then digitally sampled by an analog-to-digital convertor. Information of gamma rays, including energy and detection timing, is stored by Raspberry Pi. Details of the detection system are described in Yuasa et al. [22]. The energy range of radiation monitors varies because configurations (high voltage and lower threshold) can be changed in each year, but a 2019-2020 The absolute timing is conditioned by global positioning system (GPS) signals. When GPS signals are successfully received, the absolute timing accuracy is ±1 µs. When signals are temporarily lost (during up to several hours), we interpolate the absolute timing with the GPS signals received before and after the lost. In this case the precision is better than 1 ms. When GPS signals are completely lost, the absolute timing is assigned by the internal clock of Raspberry Pi, conditioned by the network timing protocol (NTP) service. This case provides the timing precision of 1 second. If neither GPS signal nor internet connection is lost, we still rely on the internal clock of Raspberry Pi, but the absolute timing accuracy is unknown. In FY2016, the absolute accuracy is ±0.5 sec even when GPS signals are successfully received due to a hardware issue. This accuracy can be improved by further individual calibration as reported in Wada et al. [6], but we retain this accuracy in the present paper because this catalog analysis does not require more precision. Throughout this paper, we use Japan Standard Time (JST: UTC+9) to associate gamma-ray glows with daily variations of meteorological conditions in the local time.

A. Event detection
Gamma-ray glows can be identified as count-rate enhancements from background levels. The event scan has made by the method introduced in Yuasa et al. [22]. First, count-rate histories with 10-sec bins are produced. We employ the energy range above 3 MeV because this range is not affected by radon-decay-chain background variations. With the 10-sec-binned count-rate histories, the standard deviation and the mean count rate are computed every 30 minutes. Then the significance of each bin is calculated; the count rate is divided by the standard deviation after the mean count rate is subtracted. If the significance of a bin exceeds the 5 standard deviations (5σ), it is considered as a potential glow event. Finally glow events are added to our event list by checking the count-rate histories manually. In this method, fainter events with a significance of less than 5σ are usually ignored. On the other hand, glows with a significant of less that 5σ can be coincidentally identified when searching over the data around the automatically-detected events. Such fainter glows are also added to the present event list. Most of the events were detected on thunder days in Kanazawa, confirmed by reports of Japan Meteorological Agency. The others were confirmed by a meteorological radar to coincide with a passage of a strong radar echo above the gamma-ray detector, During the four winter seasons, in total 70 gamma-ray glows were detected. The event list is shown in Table II. Among the 70 gamma-ray glows, 66 glows have a significance of more than 5σ, and the remaining 4 glows have less significance. Twenty-four events have been reported in our previous publications [21][22][23][24][25]. There are simultaneous detections with nearby detectors (events 2 and 3 for example). They can be interpreted as an identical gamma-ray glow that passed over multiple radiation detectors with ambient wind flow [22,24]. In the present paper, however, we consider a gamma-ray enhancement recorded by a radiation detector as an event in order to avoid interpretations in the event list. Figure 3 shows the number of detected glow per day. The largest number of gamma-ray glows were detected on 17 February 2020 (8 events). The second largest number is 7 events on 10 January 2018 and 17 December 2018. However, this data does not indicate the occurrence rates of gamma-ray glows because the number of deployed detectors differs as shown in Figure 2. Figure 4 presents the same data as Figure 3 but normalized by the number of radiation detectors in operation shown in Figure 2. On 10 January 2018, 8 detectors were in operation and 7 events were detected. Hence in average 0.88 events were detected by a detector in the day. Nearly one event can be detected by one detector on days under heavy thunderstorms. The first peak of events typically comes in the beginning or middle of December. Also there is sometimes the second peak in the middle of January.
The monthly variation of glow detections is shown in Figure 5. The upper panel presents accumulated data for 4 seasons. The largest number of glows were detected in December, while no events in October and November. This data should be also normalized by the number of detectors in operation as in the lower panel of Figure 5. December has still the largest number of glows, and most of the events were detected during December to February. In average every detector registered nearly one event in December. Though the detectors continue observation, glow events in March are quite rare (only event 70 on 16 March 2020). A detailed discussion is given in the Discussion section.
As well, the yearly variation of event detections is shown in Figure 6. In FY2017 (the 2017-2018 season), the largest number of events (25 glows) were detected. This trend is still true even after normalized by the number of deployed detectors as in the lower panel. FY2017 was the most significant season for glow observations in four years; in average 3-4 events were detected by a detector during the season. The detection rates in the other three seasons were similar; 1-2 events by a detector during a season in average.
The variations in the number of events by observation sites are shown in Figure 8. As mentioned in Wada et al. [21], the largest number of gamma-ray glows were detected at site 5  shows the normalized distribution by operational days. The highest frequency of event detections is at site 10 (3.68 events/120 days). This site was operated in FY2018 and FY2019, but its operation was accidentally suspended in January in FY2019. Therefore the frequency is high for site 10 even though the total number of detection is only 5.

B. Temporal analysis
This subsection presents temporal analyses of individual events based on count-rate histories. Count-rate histories of all the 70 events are shown in Appendix A. First of all, the gamma-ray glows can be apparently classified into three types based on the shape of the histories; temporally-symmetric type with a single peak, temporally-asymmetric one with single or multiple peaks, and lightning-terminated one. The classification is included in Table II. Among the 70 detected glows, 34 events are the temporally-symmetric type, 17 are the temporally-asymmetric one, and the remaining 19 events are terminated with lightning discharges.
The temporally-asymmetric type might be further divided into two types: single-peak and multiple-peak types. Clear examples of the multiple-peak type are Events 22 and 52. We can easily identify two peaks in an event. On the other hand, the boundary between the two types is vague. For example, Event 29 seems to have two peaks, but they are not well separated and hence can be also considered as a single-peaked asymmetric type. The lightning-terminated type has a sudden decrease in count rates. For most of the termination events except Event 38, Japanese Lightning Detection Network (JLDN) reported one or multiple nearby lightning currents within 10 km from the observation sites at the moment of termination, detected by the low-frequency measurement. Furthermore, the termination of Events 25, 58, and 59 coincided with a downward TGF. Events 25 and 58 were reported in Wada et al. [24] and Hisadomi et al. [25], respectively. The other cases coincided with no downward TGFs.
To examine the brightness of the gamma-ray glows, the number of detected photons are calculated. To eliminate the background variations by radon-decay chain, here again signals above 3 MeV are counted. The method to constrain the number of photons is as follows; first a source window that covers the entire event is set, and photons above 3 MeV in the windows are counted up. Next, another window that have the same duration as the source window is set to a background period, and photons above 3 MeV in the windows are also counted up. The background period of each event is described in the public materials (see the Acknowledgment section). Then, the number of photons in the background window is subtracted from that in the source window. The uncer-tainty of this method is the difference in the photons in the background window and the actual background photons in the source window. Therefore, we consider the Poisson error of the photon number in the background window as the error of this estimation. This calculation is applied to the glows with > 5σ detection significance. For the glows terminated with a downward TGF, photons after the TGFs are not included to properly evaluate the gamma-ray glows.
The results are shown in Table II and Figure 9. Due to the differences in detector responses / effective areas of the BGO and CsI crystals, the results are separately displayed. The largest number is 24081 photons and 4771 photons for detectors with BGO and CsI crystals, respectively. Its geometric mean is 1105 photons and 458 photons, and the median is 1024 photons and 383 photons, respectively.
Another key parameter is duration. Here we define T80 duration; the duration window starts when the accumulated photon number exceeds 10% of the total detected photon number calculated above, and ends when the accumulated number exceeds 90% of the total photon number. A similar method T90 is often used for gamma-ray bursts in astrophysics [56]. The calculation was executed with 0.5-sec steps to consider deadtime correction. The results are listed in Table II and also displayed in Figure 9. As well as the photon number, this calculation is applied to the glows with > 5σ detection significance, and the results from the BGO and CsI crystals are separated. The longest duration is 161.5 sec and 110.5 sec for detectors with BGO and CsI crystals, respectively. The average duration is 63.2 sec and 49.7 sec, and the median is 56.0 sec and 48.0 sec, respectively. When the results of BGO and CsI crystals are mixed up, the average duration is 58.9 sec.
Besides the calculation of the photon number and duration, we test a function fitting of the count-rate histories. Since the count-rate histories, in particular the temporally-symmetric type, may be approximated by a Gaussian function, Gaussian fitting to the histories are examined. The fitting function F (t) is where t is time, C bgd a constant value of the background count rate, a the peak count rate, T the peak time, σ i the standard deviation of the Gaussian function, and n the number of Gaussian functions to be added. n should be as small as possible. In most cases of the temporallysymmetric type, one Gaussian function can reproduce the count-rate histories, and hence n = 1. On the other hand, n can be more than 1 for a part of the temporallysymmetric type that cannot be approximated by a Gaussian function, and for the temporally-asymmetric type. This fitting is also applied to a part of the lightningterminated type if the events have a count-rate peak before the termination. In this case, a step function is folded to Equation 1 to imitate the termination. The moment of termination is set to be that of the lightning discharge reported by JLDN. It is noted that there is no physical reason to utilize the Gaussian function, but this is a tentative attempt to evaluate the shape of the count-rate histories.
The fitting results are overlaid in Figure 17 to 23. Most of the symmetric-type and a part of lightning-terminated events can be approximated by a single Gaussian function plus a background constant. A symmetric-type event (Event 26) and most of the asymmetric-type events can be approximated by two Gaussian functions plus a constant. Events 22 and 60 require three Gaussian components. The moment of termination in Event 15 is different from that reported by JLDN. Because the radiation monitor in Site 3 lost both GPS signals and the internet connection at that moment, the absolute timing is unreliable, the discrepancy is not solid and thus we ignore it.

C. Spectral analysis
We then examine spectral analysis of the gamma-ray glows. The analysis is performed with XSPEC, a spectral analysis framework for high-energy astrophysics [57]. XSPEC can unfold detector responses when an appropriate fitting function and the response matrix are given to derive the incident photon fluxes.. The response matrices for the BGO and CsI scintillators were produced by Monte-Carlo simulations with the simulation framework Geant4 [58][59][60], as shown in Yuasa et al. [22]. Among the 70 events, we focus on 28 events with >1000 photons above 3 MeV in order to avoid analyses with poor photon statistics. Also, a background spectrum is subtracted from the source spectrum. Background spectra can change below 3 MeV as this energy band is affected by the radon-decay chain emission. It can lead a large systematic uncertainty of the background-subtracted spectra below 3 MeV if the intrinsic glow flux is low. Therefore, the spectral analysis is applied to the events with high photon statistics. The source spectrum for each event is extracted from the T80 window, namely the window between two bluedashed lines in Figure 17 to 23 of Appendix A. The window for the background spectrum is set to have the same duration as the source window, not to include the gamma-ray glow, but to include a nearby background period. The background-subtracted spectra are fitted by a power-law function with an exponential cutoff S(E) where E, A, Γ, and E cut are energy of gamma-ray photons, normalization factor, power-law index, and cutoff energy, respectively. This is an approximated form of a bremsstrahlung spectrum produced by RREA electrons [61].
by extrapolating the gamma-ray spectrum down to 0.2 MeV with best-fit parameters. The fitting energy range varies by events because the energy range depends on detectors, due to differences in amplifier gain, lower and upper thresholds. The fitting is performed to cover the full energy range of each event, but uses above 0.6 MeV for all the events because the range below 0.6 MeV can be affected by the annihilation line at 0.511 MeV and hence not be fitted only with the powerlaw function. All the response-folded and response-unfolded spectra are shown in Figure 24 to Figure 29. The best-fit functions are also overlaid. The fitting parameters are listed in Table IV. Histograms of the parameters are also shown in Figure 10. The power-law index is 0.50 on average, with a standard deviation of 0.28. As well the average of the cutoff energy is 4.41 MeV, with a standard deviation of 0.41 MeV. The energy flux varies from 1.5 × 10 −6 to 8.1 × 10 −5 erg cm −2 s −1 , with the geometric mean of 8.4 × 10 −6 erg cm −2 s −1 . Scatter plots between two parameters of three are shown in Figure 11. There is no significant correlations between the energy flux and  the power-law index nor between the energy flux and the cutoff energy. This feature is different from observations in Aragats [53], which reports a clear positive correlation between power-law index and intensity at 1 MeV. The power-law index and the cutoff energy have a moderate negative correlation with a correlation coefficient of 0.47.
We then extract an averaged spectrum to derive typical parameters of glow spectra. Here the energy spectra of 21 glows detected by the BGO crystals are accumulated because the spectra of BGO and CsI crystals cannot be summed up due to the difference in detector responses. Also Events 3 and 6 were excluded from this analysis because the energy range of these events are narrower than the others, and hence they would reduce the quality of the analysis. The fitting range is set to 0.5-20.

D. Meteorological conditions
Meteorological conditions, including ambient wind and atmospheric temperature are of great interest for investigation of gamma-ray glows. The wind speed and its direction at the moment of glow detections can be derived from meteorological radar data, by calculating mov- ing vector with two time-separated radar contours. The detail method of the calculation is introduced in Wada et al. [24]. This method works under the assumption that a cloud or an echo region moves with ambient wind. We utilize the data of the eXtended RAdar Information Network (XRAIN), operated by Ministry of  Table II. Also their histograms are shown in Figure 13. Most of the gamma-ray glows occurred during west or westsouthwest winds, and a few during west-northwest winds. Gamma-ray glows never occurred during north, south, nor east winds. Atmospheric temperature during glows and its histogram are also presented in Table II and Figure 14, respectively. The temperature was measured at Komatsu AMeDAS Station (N36.382 • , E136.436 • ; for the events at sites 1 and 2) and Kanazawa Meteorological Observatory (N36.589 • , E136.634 • ; for the other events) of Japan Meteorological Agency. The histories of the temperature measurement can be retrieved from the web site of Japan Meteorological Agency. The average temperature is 4.6 • C with a standard deviation of 2.4 • C. Most cases occurred at a temperature above 0 • C, but two cases were at below 0 • C (Events 7 and 30).

IV. DISCUSSION
Almost all the glows were detected in December to February, except one events in March 2020. In particular, the detection frequency is the highest in December as shown in the lower panel of Figure 5. Winter thunderstorms in the Hokuriku region usually starts at the end of November or the beginning of December, and ends at the beginning of March. According to the records by Japan Meteorological Agency, the averaged thunder days in Kanazawa during the four winter seasons (FY2016 to FY2019) are 2.3, 5.5, 6.0, 3.0, and 2.8 days in November, December, January, February, and March, respectively. Therefore, the monthly variation of glow detections roughly corresponds to that of thunder days, even though the relation between them is not necessarily proportional. The thunder days in December to February are 17, 18, 10, and 10 days in FY2016, FY2017, FY2018, and FY2019, respectively. It is also consistent that FY2017, in which the largest number of glows were detected, had the largest number of thunder days. We note that during 1981-2020, the averaged thunder days in three months from December to February are 21.1 days (according to Japan Meteorological Agency), significantly higher than in the four recent years presented here.
As in the Introduction section, a typical cloud base of winter thunderclouds is several hundreds meters, significantly lower than summer ones, and a typical cloud-top altitude is 5-6 km. [6,21,30] Since the charge center is close to the ground, gamma rays can reach detectors at sea level before they are attenuated in the atmosphere. The average temperature at the moment of glow detections is 4.6 • C with a standard deviation of 2.4 • C, which is significantly lower than one when summer thunderstorms occur. Also the westerly wind is predominant in most cases. In the Hokuriku region, including Kanazawa and Komatsu, the north wind typically causes cold air to flow in, causing the surface temperature to drop below freezing and cause snowfall. At this time, since the altitude at 0 • C where charging occurs reaches the ground, lightning does not occur in many cases, with some exceptions. Maintaining a surface temperature of around 5 • C with the westerly wind is a typical meteorological condition for gamma-ray glow. We note that it is necessary to discuss whether the conditions under which winter thunderstorms occur and ones under which gamma-ray glows occur are the same or different. A further study of win- ter thunderstorms themselves is needed to answer this question.
As described in the Results section, the majority of the gamma-ray glows has a symmetric count-rate history that can be approximated by one Gaussian function. This is suggestive of the idea that the irradiation area has a circle or ellipse shape, and passes above a radiation detector with a constant wind flow. A schematics of this idea is shown in Figure 15. In fact, Events 2 and 3 can be explained by a passage of a thundercloud, presented in Yuasa et al. [22] For the temporally-asymmetric type, there are several interpretations. One is that multiple electron-acceleration regions adjoin each other, as in Figure 15. This idea is supported by the result that the temporally-asymmetric glows are approximated by two or three Gaussian functions. In fact, several temporallyasymmetric glows have two peaks clearly (Events 22, 52, and 60) and they can be explained by this idea. Another one is that the irradiation area does not have a circle nor ellipse shape, but a more complex shape. This can be caused by a complex shape or a tilt of electronacceleration region. It is also possible that an avalanche region rapidly glows and decays [25,62]. It is difficult to distinguish these ideas for the temporally-asymmetric type with a single peak. As mentioned before, using Gaussian functions is tentative and empirical just to evaluate the shape of the gamma-ray glows, without enough physical reasoning. Therefore, some temporallyasymmetric events could be explained and evaluated by other functions.
The spectral fittings allow us to reconfirm that the energy spectra of gamma-ray glows with high photon statistics can be approximated by the power-law function with an exponential cutoff, and hence can be explained by bremsstrahlung emission from electrons. The average cutoff energy, 4.41 MeV with a standard deviation of 0.41 MeV, is close to the average energy of RREA electrons (∼7.3 MeV [61]). As in Fig. 3 of Dwyer et al. [61], the average energy of RREA electrons with an electric field below 0.5 MV m −1 can be lower than 7.3 MeV. Since gamma-ray glows are considered to occur in an electric field close to the RREA threshold, being different from TGFs, the obtained cutoff energy might reflect the electric field strength inside thunderclouds. However, the spectral shape is also affected by atmospheric scattering and absorption, and hence we need further investigations  considering atmospheric propagation of gamma rays with Monte-Carlo simulations. As seen in the lower panel of Figure 11, there is a positive correlation between power-law index and cutoff energy. However, these two parameters are not independent. Figure 16 shows a confidence area at one standard deviation in the space of the power-law index and the cutoff energy, derived from the averaged energy spectrum.
As the confidence area is a narrow ellipse, the parameters are coupled. It means that a function with a larger index and a larger cutoff energy is similar to that with a smaller index and a smaller cutoff energy. Therefore, this positive correlations is caused by the characteristics of the fitting function, and do not necessarily have physical meanings.
Chilingarian et al. [9] published a catalog of gammaray glows/TGEs detected at Mt. Aragats in Armenia in 2017. It is valuable to compare it with the present work. In 2017, 44 TGEs with high-energy particles (above 3 MeV) were detected at the Aragats Space Environment Center (ASEC). One of the significant differences is the seasons. They are all connected to summer thunderstorms from April to November while our cases were during winter thunderstorms. Since ASEC is located at an altitude of 3200 m, a typical temperature at the moment of TGE detections ranges from −3 • C to 3 • C even in summer, lower than the average temperature of our cases.
Another significant difference is the number of detections; 44 TGEs from April to November, and hence 8.8 TGEs per month at a single site in their case. In the present analysis, the detection frequency is as many as 0.93 events per month per observation site. Therefore, TGEs/gamma-ray glows were detected more frequently at ASEC than in the Hokuriku region. This is potentially caused by the difference in meteorological conditions, namely the frequency of lightning discharges at Aragats and Hokuriku. Besides that, there is also a difference in atmospheric temperature. At ASEC, TGEs were detected frequently even below 0 • . In these cases, the cloud base can get closer to the ground, and it is possible that sometimes the detectors are inside thun-derclouds. This indicates that the detectors at ASEC is often closer to the acceleration region than our cases, and they could have more opportunities to detect TGEs due to less atmospheric absorption of gamma rays.

V. CONCLUSION
Thanks to the mapping observation campaign with up to 10 radiation monitors for four winter seasons, we observed 70 gamma-ray glows during winter thunderstorms in Japan, which allows us a statistical analysis. Detections of almost all the glows concentrated on December, January, and February, in which winter thunderstorms are frequent in the Hokuriku region, and 77% were detected in nighttime. A typical T80 duration is minutescale, 58.9 sec on average. A half of the detected glows has a temporally-symmetric (Gaussian-like) count-rate variation, suggestive of a stable electron accelerator passing above the detectors with a constant wind flow. A quarter has a temporally-asymmetric variation, implying multiple acceleration regions or one region with a complex shape. The rest quarter terminated with lightning discharges. Energy spectra are well approximated by the power-law function with an exponential cutoff, whose average power-law index, cutoff energy, and energy flux are 0.50, 4.41 MeV, and 8.4 × 10 −6 erg cm −2 s −1 , respectively. During the glow detections, winds flow typically eastward, and the temperature is 4.6 • C on average. It will be valuable if the enriched data of the present catalog is further investigated with other observational, theoretical and simulation works.  Figure 24 to Figure 29, with responsefolded and response-unfolded forms. All the data points are disclosed on the public data archive for the scientific community and the general public as mentioned in the  Count rate (counts s -1 ) Figure 18: Count-rate histories of Events 11-20 above 3 MeV in the same format as Figure 17. The bin width is 5, 3, 3, 3, 3, 3, 5, 10, 3, and 3 sec for Events 11-20, respectively. Count rate (counts s -1 ) Figure 19: Count-rate histories of Events 21-30 above 3 MeV in the same format as Figure 17. The bin width is 3, 3, 5, 3, 3, 3, 3, 3, 3, and 3 sec for Events 21-30, respectively.  Locat time (JST) 23 Count rate (counts s -1 ) Figure 21: Count-rate histories of Events 41-50 above 3 MeV in the same format as Figure 17. The bin width is 3, 3, 3, 3, 3, 3, 3, 5, 2, and 2 sec for Events 41-50, respectively. Count rate (counts s -1 ) Figure 22: Count-rate histories of Events 51-60 above 3 MeV in the same format as Figure 17. The bin width is 3, 3, 3, 3, 2, 2, 2, 2, 2, and 3 sec for Events 51-60, respectively. Count rate (counts s -1 ) Figure 23: Count-rate histories of Events 61-70 above 3 MeV in the same format as Figure 17. The bin width is 5, 3, 5, 3, 3, 3, 3, 3, 3, and 3 sec for Events 61-70, respectively.