Thundercloud Project: Exploring high-energy phenomena in thundercloud and lightning

We designed, developed, and deployed a distributed sensor network aiming at observing high-energy ionizing radiation, primarily gamma rays, from winter thunderclouds and lightning in coastal areas of Japan. Starting in 2015, we have installed, in total, more than 15 units of ground-based detector system in Ishikawa Prefecture and Niigata Prefecture, and accumulated 551 days of observation time in four winter seasons from late 2015 to early 2019. In this period, our system recorded 51 gamma-ray radiation events from thundercloud and lightning. Highlights of science results obtained from this unprecedented amount of data include the discovery of photonuclear reaction in lightning which produces neutrons and positrons along with gamma rays, and deeper insights into the life cycle of a particle-acceleration and gamma-ray-emitting region in a thundercloud. The present paper reviews objective, methodology, and results of our experiment, with a stress on its instrumentation.


Introduction
Lightning discharges and thunderclouds have been known as electrical phenomena in the atmosphere since the discovery by Benjamin Franklin in 1752. Thanks to recent observational and theoretical studies, they have been also found to be closely associated with high-energy phenomena comprising high-energy photons, electrons, neutrons, etc., and a new academic field "high-energy atmospheric physics" has been established. In 1925, Wilson [1] proposed the first idea that strong electric fields in thunderclouds can accelerate β-particles or electrons of cosmic-ray origin to MeV energies, even in the dense atmosphere. Electrons accelerated in electric fields emit bremsstrahlung photons by colliding with atmospheric nuclei. This Wilson's runaway electron scheme was developed with multiplication processes into relativistic runaway electron avalanches (RREA) by Gurevich et al. [2]; secondary electrons produced by accelerated electrons via ionization loss processes also become seed electrons and are accelerated in electron fields.
The first reports of high-energy atmospheric phenomena were made by Parks et al. [3] and McCarthy and Parks [4]. They utilized an X-ray counter onboard an F-108 aircraft and detected enhancements of count rates lasting for tens of seconds while flying in thunderclouds. This phenomenon is now called "gamma-ray glow". Gamma-ray glows originate from electron acceleration and multiplication in thunderclouds. Their duration ranges from seconds to tens of minutes; their life cycle is thought to be connected to the stability of electric fields inside thunderclouds. So far, gamma-ray glows have been observed by aircrafts [5][6][7], balloons [8], and mountain-top experiments [9][10][11][12][13][14]. When they are detected by ground-based facilities, they are also referred to as thunderstorm ground enhancements (TGEs) [15]. In particular, the observatory at Mount Aragats in Armenia has observed the largest number of TGEs by cosmic-ray monitors [15][16][17]. Gamma-ray glows are sometimes quenched by lightning discharges [4-6, 10, 18, 19]. This is evidence that electric fields responsible for gamma-ray glows can be destroyed by lightning currents.
Besides airborne and mountain-top observations of gamma-ray glows, experiments during winter thunderstorms in Japan are of great importance. In coastal areas facing the Sea of Japan, northern seasonal winds blow and provide heavy snow with lightning discharges. These winter thunderstorms in Japan are distinctive comparing to typical thunderstorms, in particular cloud bases. While typical summer thunderstorms develop above an altitude of 3-km or higher, winter thunderclouds in Japan have a cloud base of lower than 1 km [20]. Gamma-ray photons are absorbed in the atmosphere typically within 1 km. Therefore, we need in-situ measurements by airborne detectors or getting closer to thunderclouds by putting detectors on mountain tops, to observe gamma rays from summer thunderstorms. On the other hand, winter thunderstorms allow us to observe high-energy atmospheric phenomena at sea level. Torii et al. [21] reported gamma-ray glows lasting for ∼1 minute during winter thunderstorms for the first time, recorded by dosimeters installed at a nuclear power facility in a coastal area of the Sea of Japan. Another measurement with multiple dosimeters succeeded in tracking a gamma-ray glow moving with a thundercloud and ambient wind flow [22].
Motivated by the initial findings in 1990s and early 2000s, we launched the Gamma-Ray Observation of Winter Thunderclouds (GROWTH) experiment in 2006. The GROWTH experiment is a ground-based measurement of gamma rays and high-energy particles aiming at detecting and exploring high-energy atmospheric phenomena during winter thunderstorms in Japan. The experiment started with a suite of gamma-ray and electron detectors installed at Kashiwazaki-Kariwa Nuclear Power Station of Tokyo Electric Power Holdings in Niigata Prefecture, Japan. The power station faces the Sea of Japan, and frequently encounters lightning discharges in winter seasons. Tsuchiya et al. [43] reported the first detection of a gamma-ray glow lasting for ∼40 sec in the Kashiwazaki-Kariwa site. Its energy spectrum, a continuum extending up to 10 MeV originating from bremsstrahlung of electrons, suggested that a thundercloud continuously accelerated electrons to 10 MeV or higher energy. Combining Monte-Carlo simulations, glows in the site were found to originate at an altitude of<1 km [44]. Tsuchiya et al. [45] reported a glow abruptly terminated with a lightning discharges. The energy spectrum of the glow gradually became hard, i.e. the ratio of >10 MeV to 3-10 MeV photons was increasing as the lightning discharge was drawing near. Umemoto et al. [46] reported an enigmatic enhancement of electron-positron annihilation gamma rays after a lightning discharge.
During the first decade of the GROWTH experiment in 2006-2015, one or two observation points were maintained at the power station. However, the sparse distribution of detectors is not sufficient to delve deeper into the nature of gamma-ray glows such as on-ground distribution of particles and the life cycle of their acceleration site, i.e. how particle acceleration is initiated, develops, and comes to an end. Therefore, we launched a new campaign of the GROWTH experiment with multiple gamma-ray detectors and observation sites, called "Thundercloud Project" in 2015. The initial scientific results of the campaign have been already reported as journal publications [47][48][49][50][51]. In this paper, we describe the design and the performance of our gamma-ray detector system followed by details of completed observation campaigns. Highlights of scientific achievements obtained based on data from these observation campaigns are also summarized.
Throughout the paper, gamma-ray glow (gamma-ray emission from the thunder cloud, typically lasting for a few minutes) and short-duration gamma-ray burst caused by a downward TGF (lasting for a fraction of second) are collectively called thundercloud radiation 3/31 bursts (TRBs). When we put a stress on the time scale of gamma-ray burst events from the observational point of view, we also call minute-lasting gamma-ray glow as "long-duration gamma-ray burst".

Experiment setup: gamma-ray detector system
At high level, our detector system is a conventional photon-counting gamma-ray spectrometer based on scintillation crystal and photo-multiplier tube (PMT). For realizing a distributed observation network of TRBs, we set miniaturization of an entire system as a primary design goal to allow easy handling and deployment in rooftop/outdoor environments. Keeping the cost of the system as low as possible is also essential because otherwise the number of detector systems manufactured would not be large due to the tight research budget and the scale of the observation network would be limited. In addition, we deploy the detector system to multiple locations (distances between detectors varying from a few to hundreds of km), frequent on-site maintenance is not an option, and remote-monitoring and remote-control capabilities are indispensable. Figure 1 shows a high-level block diagram of the detector system, which consists of 1) scintillation crystal viewed with photo-multiplier tube (hereafter sensor assembly), 2) Detector-control and data-acquisition electronics subsystem (hereafter DAQ subsystem), 3) telecommunication subsystem, and 4) mechanical support structure and waterproof enclosure. The detector is supplied with AC 100 V from the commercial power line, and a switching regulator generates DC voltages (12 V and 5 V) required by the electronics subsystem and the telecommunication subsystem. The telecommunication subsystem provides internet connectivity via a cellular network, and is used for telemetry transmission from the detector system, and remote-login to a computer in the DAQ subsystem via secure shell (ssh). Typical power consumption of the entire system is about 7 W.
In the following subsections, detailed specifications of individual subsystems are described.

Sensor assembly
For detailed temporal and spectral analyses, it is critically important to detect gamma rays from thundercloud and lightning at as high photon counts as possible. The only way to achieve this is to make the effective area of a sensor larger and to select sensor materials with high stopping power against gamma rays with energies of MeV to a few tens of MeV. Bismuth germanite (Bi 4 Ge 3 O 12 ; hereafter BGO) scintillation crystal is one of optimal crystals in the thundercloud gamma-ray observation due to its high stopping power and environmental durability (no deliquescence). In the pilot observation campaign in 2015, we have employed cylindrical BGO crystals each with a diameter of 7.62 cm and a height of 7.62 cm. The standard BGO crystals that we used in the regular observation campaign since 2016 have dimensions of 25 cm × 8 cm × 2.5 cm. One crystal is viewed with two HAMAMATSU R1924A PMTs, and outputs from the two PMTs are combined in the analog stage, and then amplified and digitized as a single signal. Each set of a crystal and two PMTs are enclosed in a 2-mm-thick aluminum case. We have used 15 of this BGO-based sensor assemblies since 2016.
During our detector development, low-cost Thallium-doped Cesium Iodide crystals, or CsI (Tl) for short, that were extracted from a terminated accelerator experiment project became available, and we have purchased a dozen of 30   area of the CsI-based sensor assembly is slightly smaller than that of BGO-based ones, but helped to expand our observation network at a moderate increment of the manufacturing cost. Figure 2 illustrates the effective area of each scintillation crystal over gamma-ray energies 0.2-20 MeV, which is a typical energy range our detectors observe. Table 1 summarizes the energy resolution of the BGO and CsI crystals, measured via the laboratory calibration (0.662 MeV from 137 Cs isotope) and using the environmental background signal (1.46 MeV and 2.61 MeV from 40 K and 208 Tl, respectively).

Detector-control and data acquisition (DAQ) subsystem
We developed a data-acquisition and detector-control system based on 1) an analog front-end board, 2) digital signal-processing (DSP) board, and 3) commercial-off-the-shelf single-board computer Raspberry Pi. The analog front-end board is a custom board designed by our group specifically for the present experiment. The DSP board is a general purpose Field Programmable Gate Array (FPGA) board with 4-ch waveform-sampling Analog-to-Digital Converters (ADCs). We developed the FPGA/ADC board in collaboration with Shimafuji Electric, primarily for our experiment but also aiming for broader applications in other projects. Uniformly-distributed gamma rays arriving from the direction of the normal of the detection surface (25 cm × 8 cm face of BGO and 30 cm × 5 cm face of CsI) are assumed in the simulation. Photo absorption, Compton scattering, and electron-positron pair creation are the physical processes involved in the simulation, and interactions that deposited energies larger than 40 keV in the crystal were considered detectable.
As shown in Fig. 3, these boards are vertically stacked using 2.54-mm-pitch board-to-board connectors, forming a standalone data acquisition system within a cube of 10×10×10 cm 3 , excluding protruding high-voltage power supply connectors. This design was chosen to save footprint of the system, and also to reduce required cabling during fabrication and integration at each observation site. Though the entire DAQ system is compact, it fully implements analog and digital signal processing required to function as a gamma-ray spectrometer and autonomously collect data for several months. Since we consider this miniaturized DAQ system as one of key enablers of our multi-point observation campaign, the design of the system is detailed in the following paragraphs. The high-level technical specification of the system is also summarized in Table 2. 2.2.1. Analog front-end board. The analog front-end board carriers high-voltage power supply (HVPS) modules, amplifier chains, a Global Positioning System (GPS) receiver, and an organic light-emitting-diode (OLED) display. The board also implements a combined temperature, pressure, and humidity sensor BME-280 for providing house-keeping information.
We selected the OPTON-1.5PA HVPS module from Matsusada Precision as our system, because of its small footprint and volume (44×30×16 mm 3 ). The board can carry up to two HVPS modules and high-voltage outputs from the modules are routed to two Safe High Voltage (SHV) connectors. The reference voltage signals of the HVPS modules are connected to a 2-ch 12-bit digital-to-analog converter MCP4822-E/MS on the digital signal processing board, so that output voltages can be flexibly controlled from software on Raspberry Pi via Serial Peripheral Interface (SPI). The amplifier chain consists of a simple charge-integration amplifier that converts a charge output of the PMTs to a voltage signal, followed by a differentiator-integrator band-pass filter and a linear amplifier. Figure 4 shows a circuit diagram of the chain. Four copies of the same amplifier chains are implemented. When a pulse of charge with a decay time of ∼300 ns is fed from the sensor assembly (BGO and PMT) to the first-stage charge-integration amplifier, an output pulse from the band-pass-filter amplifier should look like a uni-polar pulse with a ∼1 µs rise and ∼4 µs fall timescales as shown in the right panel of the figure. These timescales are sufficiently slow compared with the sampling frequency of the waveformsampling ADC on the digital signal processing board (see the next section), and therefore, the peak pulse height, which is proportional to the energy deposit in the scintillation crystal, can be accurately measured.
The OLED display is connected to the Inter-integrated Circuit (I2C) bus of Raspberry Pi via board-to-board connectors, and is controlled by a simple Python program running on Raspberry Pi that prints a status and parameters of the system, such as observation mode, high-voltage output values, an Internet Protocol (IP) address, and so on. Although the size of the display is small (∼1 inch diagonal) and the resolution is very limited (128×64 pixels), the display turned out to be very helpful in understanding the state of the DAQ system during in particular outdoor deployment works thanks to its high visibility.    in a self-trigger mode, was developed in Hardware Description Language (HDL) and programmed to FPGA. Figure 5 illustrates the high-level block diagram of the FPGA logic. Once the input voltage exceeds the trigger threshold, a predefined number of ADC values are recorded as a "waveform" in Waveform Buffer (typically covering 10 µs since the trigger), and various properties of the waveform is then computed (maximum pulse height as well as supplementary data such as ADC values of the first/last/minimum pulse heights, sample index of the maximum pulse height, and maximum derivative of waveform values). The maximum pulse height is converted to energy deposit in the crystal in the post processing, and the supplementary data can be used to verify the normal operation of the electronics (PMT, amplifier, ADC) when necessary (see e.g. [47] for use of the supplementary data in addition to the pulse height data). The derived properties are then packed into a certain data packet structure, and stored in Event Packet, and then read by the data acquisition program running on the single-board computer via USB. The source code is publicly available on our project's online repository 1 .
The analog front-end board and Raspberry Pi are connected to the DSP board via two 20×2-pin 2.54-mm pitch connectors placed near two edges of the board. The PMT output signal amplified by the analog front-end board and the 1 PPS/NMEA output from the GPS module are routed to ADC and FPGA, respectively, via the connector. The I2C signal, HVPS reference voltage from the slow DAC, and output enable signals are also passed through the connector from the DSP board to the analog front-end board.
The I2C from the analog front-end board and the SPI communication signals of the slow ADC/DAC, and the HVPS output status control signals are connected to Raspberry Pi via the other 20×2 pin header connector.
The DSP board was manufactured by Shimafuji Electric, and is available for customers as a general-purpose waveform-sampling ADC/FPGA board.  Fig. 5 High-level block diagram of the FPGA logic.

Single-board computer.
We selected the Raspberry Pi single-board computer as our platform to run software programs that control the HVPS output mode and output voltages, collect gamma-ray event data from the DSP board, read housekeeping data from the house-keeping sensors, and transmit the house-keeping data and status information to the internet. Primary reasons of the selection include its small size (85×56×17 mm 3 ), low price (<US$100 including a power adapter and an SD card), and sufficiently high performance with a quad-core ARM Cortex-A53 processor running at 1.4 GHz and 1 gigabytes (GB) main memory. The data collection program is the only performance-critical program as the processing speed of it limits the number of gamma-ray events that the entire detector can record, and is written in C++. When the detector system is powered on, the program configures the FPGA logic on the DSP board; for example, it sets enabled ADC channels, trigger threshold values, and the number of waveform samples to be recorded per trigger. When a data collection is started, the program continuously reads the data stored on the Event Packet Buffer of the DSP board, and save the read data to a file as an event list; supported output file formats are CERN/ROOT [52] and FITS [53]. Ruby and Python are used for other non-performance critical programs to expedite the development by leveraging existing software libraries provided in the ecosystem of these scripting languages such as an OLED display controller library, a digital general-purpose input/output library, and so on. 10/31 The process monitoring framework God 2 was used to run these programs as resident processes (so-called daemons) after power on. Configurations of the programs, such as HVPS output voltages, trigger threshold, enabled ADC channels, and data collection mode (i.e. whether to start data collection and HVPS output automatically after boot) are stored in a file on the non-volatile memory (micro SD card) along with the programs and the Linux operating system.
Data collected by the programs, for example, gamma-ray event-list data, house-keeping data, are also stored on the micro SD card. An external flash memory disk connected via Universal Serial Bus (USB) was also used as a back-up storage, and the data are regularly copied from the micro SD card to the external disk.

Telecommunication subsystem
Raspberry Pi in the DAQ subsystem was connected via Ethernet to a mobile WiFi router (Aterm MR04LN from NEC) that is connected to the internet over a cellular network. Due to the stringent monthly data limitation (1 GB per 1 months) of the cellular plan that was allowed by the research grant expenditure regulation, it was infeasible to transfer all the gamma-ray event list data that amount ∼5-10 GB per month to a remote data-storage server, and therefore the connectivity was primarily used to transmit the low-data-rate telemetry sent every 300 s and a digest report of gamma-ray data such as binned count-rate histories and time-integrated energy spectra.
The telemetry data were sent to a cloud-based data base, and this allowed centralized monitoring of the status of the distributed detector systems using a web-browser-based data visualization tool. Figure 6 shows an example screen shot of the temperature telemetry. During observation campaigns, occasional stoppages of Raspberry Pi, which is thought to arise from instantaneous AC power failure, were noticed as absence of telemetry data, enabling prompt actions such as power cycling (reboot) by local support personnel.
The digest report of gamma-ray data were also useful in rapidly identifying gamma-ray enhancement events originating from thundercloud and/or lightning; when an enhancement event candidate is noticed, we remote-logged in to Raspberry Pi via ssh and manually transferred a limited number of data files for in-depth analyses.
Having "bi-directional" connectivity to individual detector systems thus helped a day-today operation during observation campaigns and also contributed to reduce latency between observation and data analysis, and to expedite publication of the data. If we were unable to retrieve data remotely, the time/human resource/financial costs of frequent data retrieval, for example, once per month, would have been impractically expensive. Therefore, we consider that the cost of installing a mobile WiFi router (∼US$100) and purchasing cellular data plan for each detector system (∼US$10 per month) have been well paid off.

Mechanical structure
For operating detectors in outdoor environments where snow and sea wind are the norm, we selected the Takachi Electronics Enclosure's water-proof and dust-tight plastic enclosure family BCAR as containers of our detector systems. The dimensions of a standard enclosure we used are 35×45×20 cm 3 . A water-proof power connector was attached to one side of an enclosure to pass through AC 100 V power. When integrating an entire detector system, a sensor assembly (or multiple of them, depending on configuration), a DAQ subsystem, and a telecommunication subsystem were screw-mounted on an aluminum base plate for ruggedization, and then the base plate was screw-mounted to the base of a water-proof enclosure, as shown in Fig. 7.

Calibration and offline data analysis
After each observation campaign in winter, data stored on the detector are retrieved from each detector system, and the energy and the timing calibrations are applied as detailed in §3.1 and §3.2. Based on the energy-and time-calibrated data, gamma-ray enhancement events, both long-and short-duration ones, are searched using a count-history-based algorithm that is described in §3.3.

Energy scale calibration
During outdoor observations, the energy scale changes over time as ambient temperature and temperature of the scintillation crystal vary. Instead of actively compensating this change by for example dynamically adjusting the PMT or the analog amplifier gain, we let the detector 12/31 operate at the pre-determined fixed gain, and corrected energy scale in the offline analysis. The correction was made by fitting the prominent gamma-ray lines seen in the environmental background radiation spectrum such as lines at 1.46 MeV ( 40 K), and 2.61 MeV ( 208 Tl) in the ADC channel space, and by constructing the best-fit linear function which returns energy in MeV for a given ADC channel.
During the outdoor observation campaign, the temperature of the detector system (measured on the DSP board) varied between 25-60 • C (the high temperature occurred under the clear skies, due to heating by the direct sun light). Even with this temperature variation, energy scale did not change significantly; typical shifts of 1.46 MeV ( 40 K) and 2.61 MeV ( 208 Tl) peak centers in the ADC channel (i.e. raw voltage value before energy scale correction) were less than 3% and 4%, respectively. The 0.609 MeV line from 214 Bi, which is clearly visible when there is precipitation, is used to validate the derived energy scale, and it was confirmed that the accuracy of the linear function is better than 2% at 0.609 MeV for the BGO scintillation crystals. An example count history is shown in Fig. 8, and spectra of the environmental background radiation during fair weather and precipitation are plotted in Fig. 9. and stored in an output data file. The information is then used by the offline data processing pipeline to assign absolute time to each gamma-ray pulse which is recorded with a (freerunning) time-counter value at trigger (i.e. when its pulse height exceeded the threshold value). Figure 10 shows an example of minute-scale variation of the local clock reconstructed based on the recorded GPS-based absolute time and the free-running time counter. When the receiver tracks sufficient number of GPS satellites, the accuracy of the 1 PPS signal generated by the module is reported to be ∼10 ns based on its data sheet. The sampling of ADC (50 MHz) governs the time resolution of trigger time of each gamma-ray pulse signal to be 20 ns. The time scale of scintillation photon emission (de-excitation) in the scintillation crystals (∼a few hundred ns to 1 µs depending on crystals) and that of the band-pass-filtered pulse (∼2 µs) are longer than the 1 PPS timing accuracy and the ADC sampling interval, and jitter of these components could potentially worsen the overall time accuracy. However, based on time correlation study between our gamma-ray measurement and radio-frequency observations (for example, [49]) confirmed that an absolute time accuracy better than 1 µs is achieved in this GPS-supported time assignment mode.

Time assignment
Occasionally, the GPS receiver did not generate navigation solution (thus no time information) due to low number of satellites in the field of view. In such a case, the pulse trigger time was converted to absolute time based on the system time of Raspberry Pi which was 14 9 Energy spectra of the environmental background extracted from the time periods without precipitation (blue) and with precipitation (red), as indicated in Fig. 8. Statistical errors are plotted in the figure, but can be hardly visible due to high counting statistics.
synchronized to the public NTP server via a cellular network. An absolute of time assignment in this mode is thought to be on the order of 10-100 ms, depending mostly on the round trip time of the cellular network.

Search of TRB events
Gamma-ray enhancement events are usually found as an excess from the environmental background gamma-ray radiation, while the background itself is also variable. As shown in Figures 8 and 9, background gamma-ray count rate (∼6.5 counts s −1 ) varies significantly below 3 MeV depending on presence of precipitation, and this variation can lower sensitivity of the search. In contrast, the > 3 MeV energy range is dominated by cosmic-ray induced signals, of which count rate is almost stable, and relatively lower than of the < 3 MeV range. Therefore, to lower the contamination from low-energy (< 3 MeV) time-variable background signal, and to increase the signal-to-noise ratio, we implemented the following processes in the search algorithm; 1) count-rate history of photons with energies above 3 MeV is generated for each 30-min data chunk, 2) the 30-min count-rate history is further binned to a histogram, and a standard deviation is computed, 3) the maximum count rate in the 30-min data chunk is divided by the standard deviation to derive "significance" value, after the mean count rate is subtract, and then, 4) a potential TRB event is reported when the "significance" exceeds a threshold value. In the nominal batch analyses, we used a time bin width of 10 s and a significance threshold of 5 standard deviation; i.e. when a count-history bin contains gamma-ray counts which is more than 5 σ apart (higher count rate) from the mean of the histogram, the bin is flagged for further examination by humans.
To illustrate this event search process, Fig. 11 presents two example 10-second-binned 30min count histories, one with no significant count increase, and the other with a gamma-ray glow being detected.

Observation campaign
In 2015, we developed 4 prototype detectors, and started a multi-point observation campaign in Kanazawa City, Ishikawa Prefecture, with 3 detectors deployed in the city. The detectors were installed on the rooftop of a building of the observation sites, as exemplified in Fig. 12. In later years, we increased the number of detectors, and deployed in more observation sites in Ishikawa and Niigata Prefectures. Figure 13 presents the locations of each observation site. Table 3

Number of detected TRB events
We have applied the event search algorithm described in Section 3.3 to the data collected through the observation campaigns in the past 4 winter seasons (late 2015 to early 2019), and detected 46 long-duration bursts and 5 short-duration bursts. The two short-duration bursts 16 Table 3 The number of detectors deployed in each observation campaign, in each observation area, and the duration of each observation campaign in days.

Prefecture
Area

Multi-point detection of TRB
The primary objective of the present experiment is to measure TRB events (both shortand long-duration gamma-ray bursts) with multiple detectors located in different sites, to study the physical extent of the gamma-ray emitting region in the cloud and its potential temporal/spatial variability as well as the movement of the cloud in detail. In fact, 14 events of all the detected TRB events were simultaneously observed by multiple detectors. For example, Fig. 16 shows a gamma-ray glow event detected by two detectors in Komatsu City, Ishikawa Prefecture at ∼17:54:00−18:00:00 of December 7th, 2016 (UTC). In this event, the detector in Komatsu High School first detected enhanced gamma-ray counts starting at ∼17:54:00 and ending at ∼17:58:00. A minute later, at around 17:55:00, a similar count-rate increase was recorded by the detector in Science Hills Komatsu (the art science museum), and lasted till 18:00:00. These 3-15 MeV count-rate time profiles are well described by a gaussian function plus a constant (corresponding to background signal), where t is time, a, b, c, and d are a normalization factor, a peak-center time, a width of the gaussian, and the environmental background count rate, respectively. Table 4 lists the bestfit parameters. The normalization factors and the widths yield the total counts of gamma rays from the gamma-ray glow, in 3-15 MeV, are estimated to be ∼ 755 ± 36 counts and ∼ 3310 ± 58 counts, in Komatsu High School and Science Hills Komatsu, respectively. The separation of the two gaussian centers is 114 ± 3 s. Errors are 1 standard deviation. The location of the two detectors deployed at Komatsu High School and Science Hills Komatsu are plotted in Fig. 17, with radar echo images taken during this time period being 18  overlaid. The straight-line distance of the two sites is 1.36 km. By tracking the movement of the precipitation feature in the radar image, we estimated a wind speed of 10.9 ± 1.2 m s −1 and wind direction as shown in Fig. 18. The wind direction is consistent with a hypothesis that a gamma-ray emitting region in the thundercloud was moving from west northwest to east southeast, first traveling over Komatsu High School, and arriving Science Hills Komatsu after that. Based on the estimated wind speed (10.9 ± 1.2 m s −1 ) and the distance measured  along the wind (1.20 km), a hypothetical travel time of the gamma-ray emission region can be estimated to be 110 ± 12 s. This value is consistent within errors with the peak-time difference based on the gaussian fitting (114 ± 3 s), and therefore we consider that the wind speed and direction estimated based on the radar images are sufficiently accurate to be used in interpreting the temporal and the geometrical aspects of this particular gamma-ray glow event. 20 As mentioned above, the total gamma-ray count of Science Hills Komatsu is larger than that of Komatsu High School by a factor of 4.4. Based on this combined with the wind direction, we infer that Science Hills Komatsu was (laterally) closer to the electron acceleration region in the thundercloud, and observed less attenuated gamma rays than the other.
The high counting statistics of the Science Hills Komatsu data allowed us to extract an energy spectrum of the gamma-ray glow event, as shown in Fig. 19. The spectrum of the glow event was extracted from a time range 17:56:30-18:00:00 (UTC). The environmental background signals were extracted using two 60-s chunks of data before and after the glow event, and subtracted from that of the glow event. Based on previous spectral studies [44], we tried to characterize the spectral shape by fitting it with a power law with an exponential cutoff: where E is gamma-ray energy in MeV, and N , Γ, and E c are a normalization factor in photons cm −2 s −1 MeV −1 , a power-law photon index, and a scaling factor for the exponential cutoff, respectively. An energy response function of the detector was generated based on a Monte-Carlo simulation using the particle transport framework Geant4 [54][55][56], and was convolved with the model function during the fitting which happened in the detector countrate dimension. The χ 2 value, which was computed as a square sum of difference between the model and the data divided by the statistical error, was minimized using the Levenberg-Marquardt algorithm in the SciPy software package. With the best-fit model parameters listed in Table 5, the model reproduced the data reasonably well with no particular structure in the fit residual (middle panel of Fig. 19), with a null hypothesis probability of 7.5%. When the same spectrum was fitted with a simple power law, a significant "convex"-shaped systematic residual was seen with a large (unacceptable) χ 2 value of 450 for 47 degrees of freedom, 21/31 supporting the presence of a spectral cutoff feature. Because the electron bremssstrahlung is thought to be the primary emission process in the gamma-ray glow, the cutoff energy should have a close relation with the maximum energy of accelerated electrons (e.g. [37] for a detailed Monte-Carlo simulation study of the acceleration and the emission processes). In addition, we anticipate that statistical analyses of the spectral shape and their temporal evolution based on multi-point observation data will allow us to better constrain the properties of the electron acceleration (electric field strength, lateral extent of the acceleration region), and plan to publish a consolidated result elsewhere. 18 Fig. 19 Top panel: Gamma-ray energy spectrum of the gamma-ray glow event recorded at Science Hills Komatsu (red filled circles). Solid red curve is the best-fit power-law with exponential cutoff (Eq. 2). The model is convolved with the detector's energy response function so that the fitting is performed in the detector count-rate dimension. Middle panel: Fit residual computed as (data − model)/error. Bottom panel: The same best-fit model function as that of the top panel, but without being convolved with the detector energy response function. Note that ordinate is in units of photons cm −2 s −1 MeV −1 , which represents the photon flux arriving at the detector. 24/31

Science highlights
In this section, we review new findings and advancements of our understanding on highenergy radiation from lightning and thundercloud based on our publications which utilized data collected with our detector system.

Photonuclear reaction triggered by a downward TGF
Enoto et al. [47] reported a sub-millisecond intense gamma-ray flash (downward TGF) and a subsequent short-duration gamma-ray burst lasting for ∼ 200 ms, recorded on February 6th, 2017, at our observation site in Kashiwazaki-Kariwa nuclear power station in Niigata Prefecture.
As shown in Fig. 20, the energy spectrum of the short-duration burst consisted of an extremely "flat" or "hard" continuum (a photon indices Γ ∼ 0.5 when fitted with a powerlaw function of N × (E/1 MeV) −Γ where N and E are a normalization factor and gamma-ray energy), associated with an abrupt cutoff at ∼10 MeV. These features made the spectrum look very different from typical energy spectra of bremssstrahlung emission seen e.g. in typical gamma-ray glows (e.g. Fig. 19).
The short-duration burst was followed by a minute-lasting gamma-ray burst. The energy spectrum of this distinctive emission, in turn, predominantly consisted of electron-positron annihilation gamma-ray line at 511 keV and its Compton scattered continuum signals.
After extensive spectral, temporal, and simulation studies, we showed unequivocally that a lightning discharge emitted huge amount of energetic (> 10 MeV) gamma rays, and neutrons were produced via atmospheric photonuclear reactions (such as γ+ 14 N → 13 N+n). The shortduration burst was interpreted well as a superposition of nuclear gamma-ray lines emitted from nuclei that underwent neutron capture, and the peculiar minute-lasting annihilation gamma-ray line emission was explained as a result of β+ decay of unstable nuclei (again, produced via photonuclear reaction).
Production of neutrons via the photonuclear reaction has been suggested based on observational results [57][58][59], and theoretical studies [60][61][62], there have been multiple reports on potential detection of neutron signals from thundercloud-and lightning-related high energy radiation (for complete reference list, see [47]). Our observation provided multi-point timeresolved data that confirm a) intense gamma-ray flash that caused neutron via photonuclear reaction, b) presence of unbound neutrons (via gamma-ray lines from neutron capture), and c) 511 keV annihilation lines from β + -decay radioisotopes generated by photonuclear reaction. These formed the first comprehensive observational evidence of such an exotic photonuclear reaction happening in the Earth's dense atmosphere.

Physical properties of downward TGF
On November 24th, 2017, three of our detectors deployed in Niigata, Japan, detected four bunches of intense short-duration ( 1 ms) gamma-ray flashes (TGFs), followed by exponentially-decaying ∼ 200-ms signals, which is, again, considered to be a result of photonuclear reaction (gamma-ray signals from de-excitation of isotopes generated via neutron capture).
We analysed time-resolved gamma-ray signals from our detectors, integrated radiation dose measured by argon ionization chambers, and low-frequency radio (LF) observations. 25 Our scintillation-crystal based detectors were heavily saturated by the intense gamma-ray signals from the four pulses of the downward TGFs, and therefore could not provide the total number of gamma-rays that entered the detector nor spectral information of the TGF. Though, we were able to derive arrival times of the four TGF events with an accuracy of ∼ 200 µs. Comparison of these TGF event times against LF time-series data showed clear correlation between TGFs and positive unipolar pulses (first and second gamma-ray flashes) or bipolar pulses (third and forth ones). Compared to a scintillation-crystal based photon counter, an ionization chamber is much tolerant to high-flux radiation when effective area is approximately the same because the latter measures the amount of integrated ionization at the cost of fine time and energy resolutions. In the TGF event in question, the ionization chambers successfully provided accurate dose information at 5 locations (400-1900 m horizontally from the estimated location of TGF), as anticipated. These dose data, combined with a Monte-Carlo simulation of gamma-ray emission and propagation in the atmosphere, were used to estimate an altitude of electron acceleration to be 2.5±0.5 km from the sea level. Based on the altitude and the measured radiation dose, the total number of avalanche electrons (> 1 MeV) was computed to be 8 +8 −4 × 10 18 , which is approximately in the same range as those of accelerated electrons estimated from space-based observations of upward TGFs (4 × 10 16 -3 × 10 19 by [29]), while many of TGFs observed in space are thought to originate at altitudes higher than 8 km [29,63].

The end of gamma-ray glow from thundercloud
One of key questions the GROWTH project is set to answer is how stable electronaccelerating region starts to form, evolves over time, and disappears in thundercloud, or in other words, the life cycle of the source of gamma-ray glow. When the close phase of 26/31 the life cycle is concerned, multiple previous measurements reported abrupt termination of gamma-ray glow that coincided with lightning discharge ( [18,45,64] and references therein).
For revealing precise relationship between lightning discharge and cessation of gamma-ray glow, Wada et al. [50] analysed an abrupt-termination event that was observed in Ishikawa Prefecture, on February 11th, 2017, by combining gamma-ray data collected by our detector and one from the GODOT project [59] as well as LF data collected by multiple receivers located ∼ 50 km from the gamma-ray observation site. Although there have been previous reports of abruptly-terminated gamma-ray glows coinciding with radio frequency observations of lightning discharges that triggered the termination [18], the nature of single-site measurements of radio signal did not allow a detailed position and time correlation study between gamma-ray glows and lightning discharges.
However, in our study [48], a multi-site LF observation provided, for the first time, a fine time-and position-resolved structure of an intercloud/intracloud lightning discharge that coincided with the gamma-ray glow termination and extended over ∼ 60 km lateral area with a 300 ms duration. Time association with the LF data and the gamma-ray data revealed that the termination happened when a stepped leader of the lightning discharge passed over the gamma-ray observation site with a horizontal distance of 0.7 km. Since the discharge started prior to the abrupt termination of gamma ray about 15 km away from the gamma-ray observation site, causality in this event is obvious; lightning discharge affected the electric field structure and effectively disabled acceleration. Still, due to the long distance between the event and the LF observation sites (∼ 50 km), we were unable to resolve vertical structure of the discharge. Continued simultaneous observation in gamma ray and radio frequency is anticipated to shed light on the charge structure in the cloud in such events in the future.

Conclusion
• Aiming at multi-point observation of particle acceleration and high-energy gamma-ray emission of thundercloud and lightning, we launched a new experiment campaign called "Thundercloud Project" in 2015, and developed a new, compact gamma-ray detector system (35×45×20 cm 3 in size) each carrying BGO or CsI scintillation crystal. • We have deployed 15 detectors to four cities in Ishikawa Prefecture and Niigata Prefecture in Japan in four winter seasons in 2015-2019, and accumulated 46 long-duration and 5 short-duration gamma-ray burst events, respectively. • Some of these events, for example the short-duration burst on February 6th, 2017 in Niigata, allowed us to record the whole process of downward TGF followed by photonuclear reaction and a traveling positron-emitting isotope cloud. • On long-duration burst, we have revealed that the long-duration gamma-ray burst can be abruptly terminated by a passage of a developing lightning leader (separated by 700 m horizontally) based on February 11th, 2017 data collected in Ishikawa Prefecture [48]. This is another stepping stone for understanding the life cycle of particle acceleration region in a thundercloud. • With accurate timing information with GPS, we have been able to correlate our gammaray data with radio frequency observations, enabling multi-messenger studies of highenergy activities of thundercloud and lightning. • We will continue observation campaigns in coming winter seasons. 27/31