TY - JOUR AU - Ikkatai,, Yuko AB - Abstract We have designed, developed, and deployed a distributed sensor network to observe high-energy ionizing radiation, primarily gamma rays, from winter thunderclouds and lightning in coastal areas of Japan. Starting in 2015, we have installed a total of more than 15 ground-based detector systems 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 the 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 the objective, methodology, and results of our experiment, with an emphasis on its instrumentation. 1. 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 also been found to be closely associated with high-energy phenomena comprising high-energy photons, electrons, neutrons, etc., and a new academic field of “high-energy atmospheric physics” has been established. In 1925, Wilson [1] proposed the first idea that strong electric fields in thunderclouds can accelerate |$\beta$|-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. 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 aircraft [5–7], balloons [8], and mountain-top experiments [9–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–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 in comparison with 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 to get 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 |$\sim$|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]. Another important class of high-energy atmospheric phenomena is “terrestrial gamma-ray flash” (TGF). TGFs are transient emissions coinciding with lightning discharges. Their energy spectrum extends up to |${>} 20$| MeV [23,24], and their duration is typically several hundreds of microseconds [25]. Since their discovery by the Compton Gamma-Ray Observatory [26], they have been routinely detected by in-orbit experiments such as RHESSI [23], AGILE [24,27], Fermi [28,29], and ASIM [30,31]. TGFs are thought to be produced by |$10^{16}$|–|$10^{19}$| energetic electrons above 1 MeV [29,32]. While several models have been proposed [33–37], the mechanism that produces such an enormous number of energetic electrons is still in debate. TGFs detected from space are upward-going, namely emitted from thunderclouds into space. More recently, downward-going ones called “downward TGFs” have been detected by ground-based experiments [38–42]. Motivated by the initial findings in the 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 aimed 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 the 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. Tsuchiya et al. [43] reported the first detection of a gamma-ray glow lasting for |$\sim$|40 s at the Kashiwazaki-Kariwa site. Its energy spectrum, a continuum extending up to 10 MeV originating from electron bremsstrahlung, suggested that a thundercloud continuously accelerates electrons to 10 MeV or higher energy. Combining with 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 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 the 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 the “Thundercloud Project,” in 2015. The initial scientific results of the campaign have already been reported [47–51]. In this paper we describe the design and the performance of our gamma-ray detector system and give details of the completed observation campaigns. Highlights of the scientific achievements based on data from these observation campaigns are also summarized. Throughout the paper, gamma-ray glow (gamma-ray emission from the thundercloud, typically lasting for a few minutes) and short-duration gamma-ray bursts caused by a downward TGF (lasting for a fraction of second) are collectively called thundercloud radiation bursts (TRBs). When we put a stress on the timescale of gamma-ray burst events from the observational point of view, we also call minute-lasting gamma-ray glow a “long-duration gamma-ray burst.” 2. Experiment setup: gamma-ray detector system At high level, our detector system is a conventional photon-counting gamma-ray spectrometer based on a scintillation crystal and photo-multiplier tube (PMT). To realize a distributed observation network of TRBs, we set miniaturization of the 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 kilometers), 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) the scintillation crystal viewed with a photo-multiplier tube (hereafter the sensor assembly); (2) the detector-control and data-acquisition electronics subsystem (hereafter the DAQ subsystem); (3) the telecommunication subsystem; and (4) the mechanical support structure and waterproof enclosure. The detector is supplied with 100 V AC from the commercial power line, and a switching regulator generates the DC voltages (12 V and 5 V) required by the electronics and telecommunication subsystems. The telecommunication subsystem provides internet connectivity via a cellular network, and is used for both telemetry transmission from the detector system and remote login to a computer in the DAQ subsystem via secure shell (ssh). The typical power consumption of the entire system is about 7 W. Fig. 1. Open in new tabDownload slide Exploded view of the CAD drawing (left) and system block diagram (right) of the gamma-ray detector system. The size of the detector is 35 cm (depth) |$\times$| 45 cm (width) |$\times$| 20 cm (height). Fig. 1. Open in new tabDownload slide Exploded view of the CAD drawing (left) and system block diagram (right) of the gamma-ray detector system. The size of the detector is 35 cm (depth) |$\times$| 45 cm (width) |$\times$| 20 cm (height). The following subsections describe the detailed specifications of the individual subsystems. 2.1. Sensor assembly For detailed temporal and spectral analyses, it is critically important to detect gamma rays from thundercloud and lightning at as high a photon count as possible. The only way to achieve this is to make the effective sensor area 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) is one of the optimal scintillation crystals for thundercloud gamma-ray observation due to its high stopping power and environmental durability (no deliquescence). In the pilot observation campaign in 2015, we 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 |$\times$| 8 cm |$\times$| 2.5 cm. One crystal is viewed with two HAMAMATSU R1924A PMTs; the 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 is enclosed in a 2 mm-thick aluminum case. We have used 15 of these 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 purchased a dozen 30 cm |$\times$| 5 cm |$\times$| 5 cm crystals. The effective area of the CsI-based sensor assembly is slightly smaller than the BGO-based ones, but they helped to expand our observation network at a moderate increment to the manufacturing cost. Figure 2 illustrates the effective area of each scintillation crystal over gamma-ray energies of 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 during laboratory calibration (0.662 MeV from a |$^{137}$|Cs isotope) and using the environmental background signal (1.46 MeV and 2.61 MeV from |$^{40}$|K and |$^{208}$|Tl, respectively). Fig. 2. Open in new tabDownload slide Effective area of each scintillation crystal calculated by a Monte Carlo simulation. Uniformly distributed gamma rays arriving from the direction of the normal to the detection surface (25 cm |$\times$| 8 cm face of BGO and 30 cm |$\times$| 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. Fig. 2. Open in new tabDownload slide Effective area of each scintillation crystal calculated by a Monte Carlo simulation. Uniformly distributed gamma rays arriving from the direction of the normal to the detection surface (25 cm |$\times$| 8 cm face of BGO and 30 cm |$\times$| 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. Table 1. Typical energy resolution of the detector. Crystal . Resolution . . 0.662 MeV . 1.46 MeV . 2.61 MeV . BGO 19% 12% 9% CsI(Tl) 12% 9% 7% Crystal . Resolution . . 0.662 MeV . 1.46 MeV . 2.61 MeV . BGO 19% 12% 9% CsI(Tl) 12% 9% 7% Open in new tab Table 1. Typical energy resolution of the detector. Crystal . Resolution . . 0.662 MeV . 1.46 MeV . 2.61 MeV . BGO 19% 12% 9% CsI(Tl) 12% 9% 7% Crystal . Resolution . . 0.662 MeV . 1.46 MeV . 2.61 MeV . BGO 19% 12% 9% CsI(Tl) 12% 9% 7% Open in new tab 2.2. Detector-control and data-acquisition subsystem We developed a data-acquisition and detector-control system based on (1) an analog front-end board, (2) a digital signal-processing (DSP) board, and (3) a commercial-off-the-shelf single-board computer (a 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 four-channel waveform-sampling analog-to-digital converters (ADCs). We developed the FPGA/ADC board in collaboration with Shimafuji Electric, primarily for our experiment but also aimed at broader applications in other projects. 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 \times 10 \times 10$| cm|$^3$|⁠, excluding the protruding high-voltage power supply connectors. This design was chosen to reduce the footprint of the system, and also to reduce the cabling required during fabrication and integration at each observation site. Though the entire DAQ system is compact, it fully implements the 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. Fig. 3. Open in new tabDownload slide The manufactured DAQ subsystem. From top to bottom, the front-end analog signal-processing board, the FPGA/ADC board, and the Raspberry Pi computer can be seen. The size of the stacked subsystem is about 10 cm (depth) |$\times$| 10 cm (width) |$\times$| 10 cm (height). Fig. 3. Open in new tabDownload slide The manufactured DAQ subsystem. From top to bottom, the front-end analog signal-processing board, the FPGA/ADC board, and the Raspberry Pi computer can be seen. The size of the stacked subsystem is about 10 cm (depth) |$\times$| 10 cm (width) |$\times$| 10 cm (height). Table 2. Specification of the DAQ system. Function . Specification . Analog front-end board Amplifier Custom design; four channels; pass band |$\sim$|2 |$\mu$|s HVPS 2 |$\times$| OPTON-1.5PA (Matsusada); up to 1.5 kV GPS FGPMMOPA6H; on-chip patch antenna or external antenna via SMA connector OLED display |$128 \times 64$| pixels; 0.9 in; I2C Env. sensor Temperature, humidity, pressure with BME280 (Bosch Sensortec); I2C Stack connector |$20 \times 2$|-pin header Digital signal processing board FPGA Artix-7 XC7A35T-1FTG256C (Xilinx); system clock 50 MHz GPIO HVPS output enable; GPS 1PPS and NMEA data reception; six GPIO pins to Raspberry Pi ADC 2 |$\times$| AD9231BCPZ-65 (Analog Devices); four-channel, 12-bit, 50 Msps sampling; input range |$\pm$|5 V Slow ADC MCP3208-BI/SL (Microchip); four-channel, 12-bit sampling; SPI; four channels left for user Slow DAC MCP4822-E/MS (Microchip); two-channel, 12-bit sampling; SPI; used for HVPS reference voltage USB interface FT2232HL (FTDI Chip); USB Micro-B connector Temp. sensor LM60BIM3 (Texas Instruments); provides FPGA and DC/DC converter temperatures Current sensor LT6106HS5 (Linear Technology); I2C; provides 12 V, 5 V, 3.3 V current consumption Stack connector |$20 \times 2$|-pin socket for analog front-end board; |$20 \times 2$|-pin header for Raspberry Pi Power 12 V via 2.1 mm jack; |$\sim$|7 W power consumption in the nominal observation mode Dimensions |$9.5 \times 9.5 \times 2.9$| cm|$^3$| Raspberry Pi 3 CPU Quad-core 1.2 GHz ARM Cortex-A53 RAM 1 GB Storage 32 GB, Class 10 SD card USB Four ports Ethernet 100BASE Ethernet Function . Specification . Analog front-end board Amplifier Custom design; four channels; pass band |$\sim$|2 |$\mu$|s HVPS 2 |$\times$| OPTON-1.5PA (Matsusada); up to 1.5 kV GPS FGPMMOPA6H; on-chip patch antenna or external antenna via SMA connector OLED display |$128 \times 64$| pixels; 0.9 in; I2C Env. sensor Temperature, humidity, pressure with BME280 (Bosch Sensortec); I2C Stack connector |$20 \times 2$|-pin header Digital signal processing board FPGA Artix-7 XC7A35T-1FTG256C (Xilinx); system clock 50 MHz GPIO HVPS output enable; GPS 1PPS and NMEA data reception; six GPIO pins to Raspberry Pi ADC 2 |$\times$| AD9231BCPZ-65 (Analog Devices); four-channel, 12-bit, 50 Msps sampling; input range |$\pm$|5 V Slow ADC MCP3208-BI/SL (Microchip); four-channel, 12-bit sampling; SPI; four channels left for user Slow DAC MCP4822-E/MS (Microchip); two-channel, 12-bit sampling; SPI; used for HVPS reference voltage USB interface FT2232HL (FTDI Chip); USB Micro-B connector Temp. sensor LM60BIM3 (Texas Instruments); provides FPGA and DC/DC converter temperatures Current sensor LT6106HS5 (Linear Technology); I2C; provides 12 V, 5 V, 3.3 V current consumption Stack connector |$20 \times 2$|-pin socket for analog front-end board; |$20 \times 2$|-pin header for Raspberry Pi Power 12 V via 2.1 mm jack; |$\sim$|7 W power consumption in the nominal observation mode Dimensions |$9.5 \times 9.5 \times 2.9$| cm|$^3$| Raspberry Pi 3 CPU Quad-core 1.2 GHz ARM Cortex-A53 RAM 1 GB Storage 32 GB, Class 10 SD card USB Four ports Ethernet 100BASE Ethernet Open in new tab Table 2. Specification of the DAQ system. Function . Specification . Analog front-end board Amplifier Custom design; four channels; pass band |$\sim$|2 |$\mu$|s HVPS 2 |$\times$| OPTON-1.5PA (Matsusada); up to 1.5 kV GPS FGPMMOPA6H; on-chip patch antenna or external antenna via SMA connector OLED display |$128 \times 64$| pixels; 0.9 in; I2C Env. sensor Temperature, humidity, pressure with BME280 (Bosch Sensortec); I2C Stack connector |$20 \times 2$|-pin header Digital signal processing board FPGA Artix-7 XC7A35T-1FTG256C (Xilinx); system clock 50 MHz GPIO HVPS output enable; GPS 1PPS and NMEA data reception; six GPIO pins to Raspberry Pi ADC 2 |$\times$| AD9231BCPZ-65 (Analog Devices); four-channel, 12-bit, 50 Msps sampling; input range |$\pm$|5 V Slow ADC MCP3208-BI/SL (Microchip); four-channel, 12-bit sampling; SPI; four channels left for user Slow DAC MCP4822-E/MS (Microchip); two-channel, 12-bit sampling; SPI; used for HVPS reference voltage USB interface FT2232HL (FTDI Chip); USB Micro-B connector Temp. sensor LM60BIM3 (Texas Instruments); provides FPGA and DC/DC converter temperatures Current sensor LT6106HS5 (Linear Technology); I2C; provides 12 V, 5 V, 3.3 V current consumption Stack connector |$20 \times 2$|-pin socket for analog front-end board; |$20 \times 2$|-pin header for Raspberry Pi Power 12 V via 2.1 mm jack; |$\sim$|7 W power consumption in the nominal observation mode Dimensions |$9.5 \times 9.5 \times 2.9$| cm|$^3$| Raspberry Pi 3 CPU Quad-core 1.2 GHz ARM Cortex-A53 RAM 1 GB Storage 32 GB, Class 10 SD card USB Four ports Ethernet 100BASE Ethernet Function . Specification . Analog front-end board Amplifier Custom design; four channels; pass band |$\sim$|2 |$\mu$|s HVPS 2 |$\times$| OPTON-1.5PA (Matsusada); up to 1.5 kV GPS FGPMMOPA6H; on-chip patch antenna or external antenna via SMA connector OLED display |$128 \times 64$| pixels; 0.9 in; I2C Env. sensor Temperature, humidity, pressure with BME280 (Bosch Sensortec); I2C Stack connector |$20 \times 2$|-pin header Digital signal processing board FPGA Artix-7 XC7A35T-1FTG256C (Xilinx); system clock 50 MHz GPIO HVPS output enable; GPS 1PPS and NMEA data reception; six GPIO pins to Raspberry Pi ADC 2 |$\times$| AD9231BCPZ-65 (Analog Devices); four-channel, 12-bit, 50 Msps sampling; input range |$\pm$|5 V Slow ADC MCP3208-BI/SL (Microchip); four-channel, 12-bit sampling; SPI; four channels left for user Slow DAC MCP4822-E/MS (Microchip); two-channel, 12-bit sampling; SPI; used for HVPS reference voltage USB interface FT2232HL (FTDI Chip); USB Micro-B connector Temp. sensor LM60BIM3 (Texas Instruments); provides FPGA and DC/DC converter temperatures Current sensor LT6106HS5 (Linear Technology); I2C; provides 12 V, 5 V, 3.3 V current consumption Stack connector |$20 \times 2$|-pin socket for analog front-end board; |$20 \times 2$|-pin header for Raspberry Pi Power 12 V via 2.1 mm jack; |$\sim$|7 W power consumption in the nominal observation mode Dimensions |$9.5 \times 9.5 \times 2.9$| cm|$^3$| Raspberry Pi 3 CPU Quad-core 1.2 GHz ARM Cortex-A53 RAM 1 GB Storage 32 GB, Class 10 SD card USB Four ports Ethernet 100BASE Ethernet Open in new tab 2.2.1. Analog front-end board The analog front-end board carries 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) to provide housekeeping information. We selected the OPTON-1.5PA HVPS module from Matsusada Precision as our system, because of its small footprint and volume (⁠|$44 \times 30 \times 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 two-channel 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 the Raspberry Pi via the Serial Peripheral Interface (SPI). The amplifier chain consists of a simple charge-integration amplifier that converts the 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 |$\sim$|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 unipolar pulse with |$\sim$|1 |$\mu$|s rise and |$\sim$|4 |$\mu$|s fall timescales, as shown in the right panel of the figure. These timescales are sufficiently slow compared with the sampling frequency of the waveform-sampling 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. Fig. 4. Open in new tabDownload slide Left: A circuit diagram of the analog amplifier chain. Bypass capacitors used for operational amplifiers are not shown in this diagram for simplicity. Right: SPICE-simulated pulse shapes. The blue and pink voltage waveforms are measured at Test Points 0 and 1 of the circuit diagram, respectively, against a typical charge supplied by a BGO+PMT assembly. Fig. 4. Open in new tabDownload slide Left: A circuit diagram of the analog amplifier chain. Bypass capacitors used for operational amplifiers are not shown in this diagram for simplicity. Right: SPICE-simulated pulse shapes. The blue and pink voltage waveforms are measured at Test Points 0 and 1 of the circuit diagram, respectively, against a typical charge supplied by a BGO+PMT assembly. The OLED display is connected to the Inter-integrated Circuit (I2C) bus of the Raspberry Pi via board-to-board connectors, and is controlled by a simple Python program running on the Raspberry Pi that prints the status and parameters of the system, such as observation mode, high-voltage output values, its Internet Protocol (IP) address, and so on. Although the size of the display is small (⁠|$\sim$|1 in diagonal) and the resolution is very limited (⁠|$128 \times 64$| pixels), the display turned out to be very helpful in understanding the state of the DAQ system, in particular during outdoor deployment works, thanks to its high visibility. 2.2.2. Digital signal-processing board The DSP board is a custom-made digitizer consisting of a Xilinx Artix-7 FPGA (XC7A35T-1FTG256C) and two dual 12-bit ADCs (Analog Devices AD9231BCPZ-65) that operate at 50 million samples per second (Msps), temperature and current sensors, a USB interface (FTDI Chip FT2232H), a slow ADC/DAC, and DC/DC converters. Custom hardware logic that collects the timing and pulse height of gamma-ray signals in a self-trigger mode was developed in Hardware Description Language (HDL) and programmed to the FPGA. Figure 5 shows 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 the Waveform Buffer (typically covering 10 |$\mu$|s since the trigger), and various properties of the waveform are then computed (maximum pulse height, as well as supplementary data such as the ADC values of the first/last/minimum pulse heights, the sample index of the maximum pulse height, and the maximum derivative of the 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., Ref. [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, stored in an event packet, and then read by the data acquisition program running on the single-board computer via USB. The source code is publicly available from our project’s online repository.1 Fig. 5. Open in new tabDownload slide High-level block diagram of the FPGA logic. Fig. 5. Open in new tabDownload slide High-level block diagram of the FPGA logic. The analog front-end board and Raspberry Pi are connected to the DSP board via two |$20 \times 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 the 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, the SPI communication signals of the slow ADC/DAC, and the HVPS output status control signals are connected to the Raspberry Pi via the other |$20 \times 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. 2.2.3. Single-board computer We selected the Raspberry Pi single-board computer as our platform to run programs that control the HVPS output mode and output voltages, collect gamma-ray event data from the DSP board, read housekeeping data from the housekeeping sensors, and transmit the housekeeping data and status information to the internet. The primary reasons for the selection include its small size (⁠|$85 \times 56 \times 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 GB main memory. The data collection program is the only performance-critical program, as its processing speed limits the number of gamma-ray events that the entire detector can record; it 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 data collection is started, the program continuously reads the data stored in the event packet buffer of the DSP board, and saves the data to a file as an event list; the 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. The process monitoring framework God2 was used to run these programs as resident processes (so-called daemons) after power-on. The configurations of the programs, such as the 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 in the non-volatile memory (microSD card) along with the programs and the Linux operating system. Data collected by the programs, for example the gamma-ray event-list data or housekeeping data, are also stored on the microSD card. An external flash memory disk connected via USB was also used as back-up storage, and the data are regularly copied from the microSD card to the external disk. 2.3. Telecommunication subsystem The Raspberry Pi in the DAQ subsystem is 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 month) 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 (which amounts to |$\sim$|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 database, 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 the Raspberry Pi, thought to arise from instantaneous AC power failure, were noticed as an absence of telemetry data, enabling prompt action such as power cycling (reboot) by local support personnel. Fig. 6. Open in new tabDownload slide Web-browser-based telemetry visualization tool. This example screen shot shows the temperature of the DC/DC module on the FPGA/ADC board over a four-day period in February 2018. The top and the bottom panels are for the detector systems deployed in Niigata Prefecture and Ishikawa Prefecture, respectively. Fig. 6. Open in new tabDownload slide Web-browser-based telemetry visualization tool. This example screen shot shows the temperature of the DC/DC module on the FPGA/ADC board over a four-day period in February 2018. The top and the bottom panels are for the detector systems deployed in Niigata Prefecture and Ishikawa Prefecture, respectively. The digest reports 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 was noticed, we remotely logged in to the Raspberry Pi via ssh and manually transferred a limited number of data files for in-depth analyses. Having “bidirectional” connectivity to individual detector systems thus helped day-to-day operation during observation campaigns, and also contributed to reduce the latency between observation and data analysis, and thus to expedite publication of the data. If we were unable to retrieve data remotely, the time, human resource, and financial costs of frequent data retrieval, for example once per month, would have been impractically expensive. Therefore, we consider the cost of installing a mobile WiFi router (⁠|$\sim$|US|${\$}$|100) and purchasing a cellular data plan for each detector system (⁠|$\sim$|US|${\$}$|10 per month) were good investments. 2.4. Mechanical structure For operating detectors in outdoor environments where snow and sea wind are the norm, we selected the Takachi Electronics Enclosure waterproof and dust-tight plastic enclosure family BCAR as containers for our detector systems. The dimensions of the standard enclosure we used are |$35 \times 45 \times 20$| cm|$^3$|⁠. A waterproof power connector was attached to one side of the enclosure to pass through 100 V AC power. When integrating an entire detector system, a sensor assembly (or more than one, depending on the 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 the waterproof enclosure, as shown in Fig. 7. Fig. 7. Open in new tabDownload slide Interior of the waterproof enclosure. The BGO-based sensor assembly, wrapped with pink bubble wrap, and the DAQ subsystem are located in the center and in the left bottom corner, respectively. The small black module adjacent to the sensor assembly is the mobile WiFi router. A power strip and AC adapters for the DAQ subsystem and the mobile WiFi router are placed in the top right corner. Fig. 7. Open in new tabDownload slide Interior of the waterproof enclosure. The BGO-based sensor assembly, wrapped with pink bubble wrap, and the DAQ subsystem are located in the center and in the left bottom corner, respectively. The small black module adjacent to the sensor assembly is the mobile WiFi router. A power strip and AC adapters for the DAQ subsystem and the mobile WiFi router are placed in the top right corner. 3. 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 Sects. 3.1 and 3.2. Based on the energy- and time-calibrated data, gamma-ray enhancement events, both long- and short-duration ones, are sought using a count-history-based algorithm that is described in Sect. 3.3. 3.1. Energy scale calibration During outdoor observations, the energy scale changes over time as the ambient temperature and the temperature of the scintillation crystal vary. Instead of actively compensating for this change, for example by dynamically adjusting the PMT or the analog amplifier gain, we let the detector operate at the predetermined fixed gain, and corrected the 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 the 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|$^{\circ}$|C and 60|$^{\circ}$|C (the high temperature occurred under clear skies, due to heating by direct sunlight). Even with this temperature variation, the 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. the 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, was 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. Fig. 8. Open in new tabDownload slide Example count history of the environmental gamma-ray background, recorded by one of the BGO-based detectors in Ishikawa in February 2019, in the |${>} 400$| keV (top panel) and |${>} 3$| MeV (bottom panel) energy bands. Three days’ worth of data are shown. The count rate of the lower energy band (top panel) varies significantly after |$10^5$| s due to gamma rays from radioisotope washout due to precipitation, while that of the |${>} 3$| MeV band stays almost constant. The blue and red rectangles indicate the time periods of fair weather and intermittent rain for which the energy spectra are shown in Fig. 9. Fig. 8. Open in new tabDownload slide Example count history of the environmental gamma-ray background, recorded by one of the BGO-based detectors in Ishikawa in February 2019, in the |${>} 400$| keV (top panel) and |${>} 3$| MeV (bottom panel) energy bands. Three days’ worth of data are shown. The count rate of the lower energy band (top panel) varies significantly after |$10^5$| s due to gamma rays from radioisotope washout due to precipitation, while that of the |${>} 3$| MeV band stays almost constant. The blue and red rectangles indicate the time periods of fair weather and intermittent rain for which the energy spectra are shown in Fig. 9. Fig. 9. Open in new tabDownload slide Energy spectra of the environmental background extracted from the time periods without precipitation (blue) and with precipitation (red), as indicated in Fig. 8. The statistical errors are plotted in the figure, but are hardly visible due to the high count statistics. Fig. 9. Open in new tabDownload slide Energy spectra of the environmental background extracted from the time periods without precipitation (blue) and with precipitation (red), as indicated in Fig. 8. The statistical errors are plotted in the figure, but are hardly visible due to the high count statistics. 3.2. Time assignment The analog daughter board carries a single-frequency GPS receiver with an external patch antenna. The navigation message output and the 1 pulse-per-second (PPS) signal of the module are routed to the FPGA on the main board via the stacking connector. On the rising edge of the 1 PPS signal, the FPGA logic registers the absolute time information in the navigation message, along with the value of the free-running 48-bit time counter which is incremented at 100 MHz. The registered information is read by the DAQ software every 30 s, and stored in an output data file. The information is then used by the offline data processing pipeline to assign an absolute time to each gamma-ray pulse which is recorded with a (free-running) time-counter value at trigger (i.e. when its pulse height exceeded the threshold value). Figure 10 shows an example of the 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 a sufficient number of GPS satellites, the accuracy of the 1 PPS signal generated by the module is reported to be |$\sim$|10 ns based on its data sheet. The ADC sampling (50 MHz) governs the time resolution of the trigger time of each gamma-ray pulse signal, making it 20 ns. The timescale of scintillation photon emission (de-excitation) in the scintillation crystals (⁠|$\sim$|a few hundred nanoseconds to 1 |$\mu$|s, depending on the crystals) and that of the band-pass-filtered pulse (⁠|$\sim$|2 |$\mu$|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 a time correlation study between our gamma-ray measurement and radio-frequency observations (for example, Ref. [49]) an absolute time accuracy of better than 1 |$\mu$|s is achieved in this GPS-supported time assignment mode. Fig. 10. Open in new tabDownload slide Top: Typical relation of a counter incremented by the free-running 100 MHz clock fed to the FPGA versus the GPS time over 30 min (one observation interval). Middle: The same counter value, but with the best-fit linear model being subtracted to visualize the time-variable drift of the local oscillator. A counter value of 1000 corresponds to 10 |$\mu$|s. Bottom: The temperature measured on the FPGA board. Fig. 10. Open in new tabDownload slide Top: Typical relation of a counter incremented by the free-running 100 MHz clock fed to the FPGA versus the GPS time over 30 min (one observation interval). Middle: The same counter value, but with the best-fit linear model being subtracted to visualize the time-variable drift of the local oscillator. A counter value of 1000 corresponds to 10 |$\mu$|s. Bottom: The temperature measured on the FPGA board. Occasionally, the GPS receiver did not generate a navigation solution (thus no time information) due to a low number of satellites in the field of view. In such a case, the pulse trigger time was converted to an absolute time based on the system time of the Raspberry Pi, which was synchronized to a public NTP server via the cellular network. The absolute 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. 3.3. Search for 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 Figs. 8 and 9, the background gamma-ray count rate (⁠|$\sim$|6.5 counts s|$^{-1}$|⁠) varies significantly below 3 MeV depending on the presence of precipitation, and this variation can lower the sensitivity of the search. In contrast, the |${>} 3$| MeV energy range is dominated by cosmic-ray-induced signals, the count rate for which is almost stable, and relatively lower than for the |${<} 3$| MeV range. Therefore, to lower the contamination from the 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) The 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 the standard deviation is computed. (3) The maximum count rate in the 30 min data chunk is divided by the standard deviation to derive the “significance” value, after the mean count rate is subtracted. (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 deviations; i.e. when a count-history bin contains gamma-ray counts more than 5 |$\sigma$| different (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 s-binned 30 min count histories, one with no significant count increase, and the other with a gamma-ray glow being detected. Fig. 11. Open in new tabDownload slide Example of 24 hr count histories of photons with energies above 3 MeV, without (top panels) and with (bottom panels) the count rate exceeding the event detection threshold. The histograms in the right panels show the distribution of the count rate, with the event-detection threshold being shown as red dashed lines. The top and bottom rows show data from December 1 and December 6 2016, respectively, obtained by a detector deployed in Komatsu City, Ishikawa Prefecture. Fig. 11. Open in new tabDownload slide Example of 24 hr count histories of photons with energies above 3 MeV, without (top panels) and with (bottom panels) the count rate exceeding the event detection threshold. The histograms in the right panels show the distribution of the count rate, with the event-detection threshold being shown as red dashed lines. The top and bottom rows show data from December 1 and December 6 2016, respectively, obtained by a detector deployed in Komatsu City, Ishikawa Prefecture. 4. Results 4.1. Observation campaign In 2015 we developed four prototype detectors, and started a multi-point observation campaign in Kanazawa City, Ishikawa Prefecture, with three detectors deployed in the city. The detectors were installed on the rooftop of a building at the observation sites, as shown in Fig. 12. In later years we increased the number of detectors, and deployed at more observation sites in Ishikawa and Niigata Prefectures. Figure 13 presents the locations of each observation site. Table 3 and Fig. 14 summarize the number of detectors that were deployed during annual observation campaigns since 2015. An annual observation campaign typically extends over five to six months from October or November to March the following year; for example, the 2016 observation campaign started in November 2016 and ended in March 2017. Fig. 12. Open in new tabDownload slide Photograph of a detector system deployed on the roof of one of the observation sites. Fig. 12. Open in new tabDownload slide Photograph of a detector system deployed on the roof of one of the observation sites. Fig. 13. Open in new tabDownload slide Locations of the observation sites of the experiment. Fig. 13. Open in new tabDownload slide Locations of the observation sites of the experiment. Fig. 14. Open in new tabDownload slide Number of detectors deployed to each observation area. Fig. 14. Open in new tabDownload slide Number of detectors deployed to each observation area. 4.2. Number of detected TRB events We applied the event search algorithm described in Sect. 3.3 to the data collected through the observation campaigns in the past four winter seasons (late 2015 to early 2019), and detected 46 long-duration bursts and 5 short-duration bursts. The two short-duration bursts detected in the 2017 campaign happened during simultaneously observed long-duration bursts. Figure 15 presents a yearly histogram of the detected events. Fig. 15. Open in new tabDownload slide Top panel: Number of TRB events detected in each observation campaign. The blue and red rectangles correspond to long- and short-duration TRBs, respectively. Bottom panels: The same as the top panel, but events detected in each observation area are shown separately. Fig. 15. Open in new tabDownload slide Top panel: Number of TRB events detected in each observation campaign. The blue and red rectangles correspond to long- and short-duration TRBs, respectively. Bottom panels: The same as the top panel, but events detected in each observation area are shown separately. 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 . Year . . . 2015 . 2016 . 2017 . 2018 . Ishikawa Kanazawa 3 3 6 9 Komatsu 0 2 2 2 Suzu 0 1 0 0 Niigata Kashiwazaki 0 4 4 4 Duration (days) 94 189 127 141 Total (days) 551 Prefecture . Area . Year . . . 2015 . 2016 . 2017 . 2018 . Ishikawa Kanazawa 3 3 6 9 Komatsu 0 2 2 2 Suzu 0 1 0 0 Niigata Kashiwazaki 0 4 4 4 Duration (days) 94 189 127 141 Total (days) 551 Open in new tab 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 . Year . . . 2015 . 2016 . 2017 . 2018 . Ishikawa Kanazawa 3 3 6 9 Komatsu 0 2 2 2 Suzu 0 1 0 0 Niigata Kashiwazaki 0 4 4 4 Duration (days) 94 189 127 141 Total (days) 551 Prefecture . Area . Year . . . 2015 . 2016 . 2017 . 2018 . Ishikawa Kanazawa 3 3 6 9 Komatsu 0 2 2 2 Suzu 0 1 0 0 Niigata Kashiwazaki 0 4 4 4 Duration (days) 94 189 127 141 Total (days) 551 Open in new tab 4.3. Multi-point detection of TRB The primary objective of the present experiment is to measure TRB events (both short- and long-duration gamma-ray bursts) with multiple detectors located at different sites in order 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 out 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 |$\sim$|17:54:00–18:00:00 on December 7 2016 (UTC). In this event, the detector at Komatsu High School first detected enhanced gamma-ray counts starting at |$\sim$|17:54:00 and ending at |$\sim$|17:58:00. A minute later, at around 17:55:00, a similar count-rate increase was recorded by the detector at Science Hills Komatsu (the art science museum), and lasted until 18:00:00. These 3–15 MeV count-rate time profiles are well described by a Gaussian function plus a constant (corresponding to the background signal), $$\begin{equation} f(t) = a\times\exp\left(-\frac{(t-b)^2}{c^2}\right) + d [\mathrm{counts \, s}^{-1}] , \label{eq:count_history} \end{equation}$$(1) where |$t$| is time, and |$a$|⁠, |$b$|⁠, |$c$|⁠, and |$d$| are a normalization factor, the peak-center time, the width of the Gaussian, and the environmental background count rate, respectively. Table 4 lists the best-fit parameters. The normalization factors and the widths yield estimated total counts of gamma rays from the gamma-ray glow, in 3–15 MeV, of |${\sim} 755\pm36$| counts and |${\sim} 3310\pm58$| counts at Komatsu High School and Science Hills Komatsu, respectively. The separation of the two Gaussian centers is |$114\pm3$| s. The errors are one standard deviation. Fig. 16. Open in new tabDownload slide Gamma-ray count-rate time histories recorded by our detectors at Komatsu High School (black filled circles) and Science Hills Komatsu (red filled circles) in the 3–15 MeV energy band, with a bin size of 10 s. The error bars show the statistical errors. The solid black and red curves are the best-fit “Gaussian + constant” model functions [see Eq. (1)]. Fig. 16. Open in new tabDownload slide Gamma-ray count-rate time histories recorded by our detectors at Komatsu High School (black filled circles) and Science Hills Komatsu (red filled circles) in the 3–15 MeV energy band, with a bin size of 10 s. The error bars show the statistical errors. The solid black and red curves are the best-fit “Gaussian + constant” model functions [see Eq. (1)]. Table 4. Best-fit parameters of the count-rate time-profile model [Eq. (1)]. The errors are one standard deviation, and the count rates are in the 3–15 MeV energy band. Location . |$a$| . |$b$| . |$c$| . |$d$| . |$\chi^2$| (n.d.o.f.)|$^{*}$| . . counts s|$^{-1}$| . UTC . s . counts s|$^{-1}$| . . Komatsu High School |$6.7\pm0.4$| 17:56:26|$\pm4$| |$66\pm4$| |$2.0\pm0.1$| 143.2 (133) Science Hills Komatsu |$35.6\pm0.8$| 17:58:19|$\pm1$| |$52\pm1$| |$2.7\pm0.1$| 148.8 (133) Location . |$a$| . |$b$| . |$c$| . |$d$| . |$\chi^2$| (n.d.o.f.)|$^{*}$| . . counts s|$^{-1}$| . UTC . s . counts s|$^{-1}$| . . Komatsu High School |$6.7\pm0.4$| 17:56:26|$\pm4$| |$66\pm4$| |$2.0\pm0.1$| 143.2 (133) Science Hills Komatsu |$35.6\pm0.8$| 17:58:19|$\pm1$| |$52\pm1$| |$2.7\pm0.1$| 148.8 (133) |$^{*}$| Chi-square value of the fit, with the number of degrees of freedom in parentheses. Open in new tab Table 4. Best-fit parameters of the count-rate time-profile model [Eq. (1)]. The errors are one standard deviation, and the count rates are in the 3–15 MeV energy band. Location . |$a$| . |$b$| . |$c$| . |$d$| . |$\chi^2$| (n.d.o.f.)|$^{*}$| . . counts s|$^{-1}$| . UTC . s . counts s|$^{-1}$| . . Komatsu High School |$6.7\pm0.4$| 17:56:26|$\pm4$| |$66\pm4$| |$2.0\pm0.1$| 143.2 (133) Science Hills Komatsu |$35.6\pm0.8$| 17:58:19|$\pm1$| |$52\pm1$| |$2.7\pm0.1$| 148.8 (133) Location . |$a$| . |$b$| . |$c$| . |$d$| . |$\chi^2$| (n.d.o.f.)|$^{*}$| . . counts s|$^{-1}$| . UTC . s . counts s|$^{-1}$| . . Komatsu High School |$6.7\pm0.4$| 17:56:26|$\pm4$| |$66\pm4$| |$2.0\pm0.1$| 143.2 (133) Science Hills Komatsu |$35.6\pm0.8$| 17:58:19|$\pm1$| |$52\pm1$| |$2.7\pm0.1$| 148.8 (133) |$^{*}$| Chi-square value of the fit, with the number of degrees of freedom in parentheses. Open in new tab The locations of the two detectors deployed at Komatsu High School and Science Hills Komatsu are plotted in Fig. 17, overlaid with radar echo images taken during this time period. The straight-line distance between 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\pm1.2$| m s|$^{-1}$|⁠, with the wind direction as shown in Fig. 18. The wind direction is consistent with the 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 then arriving at Science Hills Komatsu. Based on the estimated wind speed (⁠|$10.9\pm1.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\pm12$| s. This value is consistent within errors with the peak-time difference based on the Gaussian fit (⁠|$114\pm3$| 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. Fig. 17. Open in new tabDownload slide The location of the two observation sites in Komatsu, Ishikawa Prefecture (filled circles). Precipitation intensity maps obtained by XRAIN are shown for four 5 min intervals on December 7 2016 (UTC). Fig. 17. Open in new tabDownload slide The location of the two observation sites in Komatsu, Ishikawa Prefecture (filled circles). Precipitation intensity maps obtained by XRAIN are shown for four 5 min intervals on December 7 2016 (UTC). Fig. 18. Open in new tabDownload slide Close-up view of the aerial photograph of Komatsu. The filled circles indicate the observation sites. The magenta arrow shows the wind direction estimated from the radar image analysis. The white double-headed arrow is the hypothesized shortest distance as traveled by the brightest part of the gamma-ray emission region which yielded the highest peaks in the two count-rate histories. Fig. 18. Open in new tabDownload slide Close-up view of the aerial photograph of Komatsu. The filled circles indicate the observation sites. The magenta arrow shows the wind direction estimated from the radar image analysis. The white double-headed arrow is the hypothesized shortest distance as traveled by the brightest part of the gamma-ray emission region which yielded the highest peaks in the two count-rate histories. As mentioned above, the total gamma-ray count at Science Hills Komatsu was 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 count for the gamma-ray glow event, as shown in Fig. 19. The spectrum of the glow event was extracted from the 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: $$\begin{equation} f(E) = N \times E^{-\Gamma} \exp(-E/E_\mathrm{c})\label{eq:cutoffpl}~\mathrm{photons}\,\mathrm{cm}^{-2}\,\mathrm{s}^{-1}\,\mathrm{MeV}^{-1} , \end{equation}$$(2) where |$E$| is the gamma-ray energy in MeV, and |$N$|⁠, |$\Gamma$|⁠, and |$E_\mathrm{c}$| are the normalization factor in photons cm|$^{-2}$| s|$^{-1}$| MeV|$^{-1}$|⁠, the power-law photon index, and the scaling factor for the exponential cutoff, respectively. An energy response function for the detector was generated based on a Monte Carlo simulation using the particle transport framework Geant4 [54–56], and was convolved with the model function during the fitting, which happened in the detector count-rate dimension. The |$\chi^2$| value, which was computed as the 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) |$\chi^2$| value of 450 for 47 degrees of freedom, supporting the presence of a spectral cutoff feature. Because 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 (see, e.g., Ref. [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 we plan to publish a consolidated result elsewhere. Fig. 19. Open in new tabDownload slide Top panel: Gamma-ray energy spectrum of the gamma-ray glow event recorded at Science Hills Komatsu (red filled circles). The 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 fit 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 the top panel, but without being convolved with the detector energy response function. Note that the ordinate is in units of photons cm|$^{-2}$| s|$^{-1}$| MeV|$^{-1}$|⁠, which represents the photon flux arriving at the detector. Fig. 19. Open in new tabDownload slide Top panel: Gamma-ray energy spectrum of the gamma-ray glow event recorded at Science Hills Komatsu (red filled circles). The 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 fit 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 the top panel, but without being convolved with the detector energy response function. Note that the ordinate is in units of photons cm|$^{-2}$| s|$^{-1}$| MeV|$^{-1}$|⁠, which represents the photon flux arriving at the detector. Table 5. Result of energy-spectral model fitting with a power law with an exponential cutoff [Eq. (2)] to the Science Hills Komatsu data. The errors are at the 90% confidence level. N . |$\Gamma$| . |$E_\mathrm{c}$| . Energy flux|$^{*}$| . |$\chi^2$| (n.d.o.f.)|$^{\dagger}$| . ph cm|$^{-2}$| s|$^{-1}$| MeV|$^{-1}$| . . MeV . MeV cm|$^{-2}$| s|$^{-1}$| . . |$0.158^{+0.015}_{-0.016}$| |$0.26^{+0.14}_{-0.15}$| |$4.10^{+0.51}_{-0.33}$| 1.18 60.5 (46) N . |$\Gamma$| . |$E_\mathrm{c}$| . Energy flux|$^{*}$| . |$\chi^2$| (n.d.o.f.)|$^{\dagger}$| . ph cm|$^{-2}$| s|$^{-1}$| MeV|$^{-1}$| . . MeV . MeV cm|$^{-2}$| s|$^{-1}$| . . |$0.158^{+0.015}_{-0.016}$| |$0.26^{+0.14}_{-0.15}$| |$4.10^{+0.51}_{-0.33}$| 1.18 60.5 (46) |$^{*}$| Energy flux in the 3–15 MeV energy band. |$^{\dagger}$| Chi-square value of the fit, with the number of degrees of freedom in parentheses. Open in new tab Table 5. Result of energy-spectral model fitting with a power law with an exponential cutoff [Eq. (2)] to the Science Hills Komatsu data. The errors are at the 90% confidence level. N . |$\Gamma$| . |$E_\mathrm{c}$| . Energy flux|$^{*}$| . |$\chi^2$| (n.d.o.f.)|$^{\dagger}$| . ph cm|$^{-2}$| s|$^{-1}$| MeV|$^{-1}$| . . MeV . MeV cm|$^{-2}$| s|$^{-1}$| . . |$0.158^{+0.015}_{-0.016}$| |$0.26^{+0.14}_{-0.15}$| |$4.10^{+0.51}_{-0.33}$| 1.18 60.5 (46) N . |$\Gamma$| . |$E_\mathrm{c}$| . Energy flux|$^{*}$| . |$\chi^2$| (n.d.o.f.)|$^{\dagger}$| . ph cm|$^{-2}$| s|$^{-1}$| MeV|$^{-1}$| . . MeV . MeV cm|$^{-2}$| s|$^{-1}$| . . |$0.158^{+0.015}_{-0.016}$| |$0.26^{+0.14}_{-0.15}$| |$4.10^{+0.51}_{-0.33}$| 1.18 60.5 (46) |$^{*}$| Energy flux in the 3–15 MeV energy band. |$^{\dagger}$| Chi-square value of the fit, with the number of degrees of freedom in parentheses. Open in new tab 5. Science highlights In this section we review the new findings and advancement in our understanding of high-energy radiation from lightning and thundercloud based on our publications utilizing data collected with our detector system. 5.1. Photonuclear reaction triggered by a downward TGF Enoto et al. [47] reported a submillisecond intense gamma-ray flash (downward TGF) and a subsequent short-duration gamma-ray burst lasting for |${\sim} 200$| ms, recorded on February 6 2017, at our observation site at the 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 (photon index |$\Gamma\sim0.5$| when fitted with a power-law function of |$N \times (E/1\,\mathrm{MeV})^{-\Gamma}$|⁠, where |$N$| and |$E$| are the normalization factor and gamma-ray energy), associated with an abrupt cutoff at |$\sim$|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). Fig. 20. Open in new tabDownload slide Observed gamma-ray spectra of the short-duration gamma-ray burst (black filled circle) and model spectrum constructed based on the Monte Carlo simulation of neutron-induced nuclear gamma rays (the green dash-dotted, blue dashed, and purple dotted curves). The red solid curve shows the sum of the individual model components. Fig. 20. Open in new tabDownload slide Observed gamma-ray spectra of the short-duration gamma-ray burst (black filled circle) and model spectrum constructed based on the Monte Carlo simulation of neutron-induced nuclear gamma rays (the green dash-dotted, blue dashed, and purple dotted curves). The red solid curve shows the sum of the individual model components. The short-duration burst was followed by a minute-long gamma-ray burst. The energy spectrum of this distinctive emission, in turn, predominantly consisted of the 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 a huge amount of energetic (⁠|$>10$| MeV) gamma rays, and neutrons were produced via atmospheric photonuclear reactions (such as |$\mathrm{\gamma + ^{14}N \to ^{13}N + n}$|⁠). The short-duration burst was interpreted well as a superposition of nuclear gamma-ray lines emitted from nuclei that underwent neutron capture, and the peculiar minute-long annihilation gamma-ray line emission was explained as a result of |$\beta+$| decay of unstable nuclei (again, produced via a photonuclear reaction). Production of neutrons via the photonuclear reaction has been suggested based on observational results [57–59] and theoretical studies [60–62]; there have been multiple reports on potential detection of neutron signals from thundercloud- and lightning-related high-energy radiation (for a complete reference list, see Ref. [47]). Our observation provided multi-point time-resolved data that confirm (a) an intense gamma-ray flash that caused neutrons via the photonuclear reaction, (b) the presence of unbound neutrons (via gamma-ray lines from neutron capture), and (c) 511 keV annihilation lines from |$\beta^{+}$|-decay radioisotopes generated by the photonuclear reaction. These formed the first comprehensive observational evidence of such an exotic photonuclear reaction happening in the Earth’s dense atmosphere. 5.2. Physical properties of downward TGF On November 24 2017, three of our detectors deployed in Niigata, Japan, detected four bunches of intense short-duration (⁠|$\ll1$| ms) gamma-ray flashes (TGFs), followed by exponentially decaying |${\sim} 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 analyzed time-resolved gamma-ray signals from our detectors, the integrated radiation dose measured by argon ionization chambers, and low-frequency radio (LF) observations. Our scintillation-crystal-based detectors were heavily saturated by the intense gamma-ray signals from the four downward TGF pulses, and therefore could not provide the total number of gamma-rays that entered the detector nor any spectral information on the TGF, though we were able to derive arrival times of the four TGF events with an accuracy of |${\sim}200$||$\mu$|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 fourth ones). Compared to a scintillation-crystal-based photon counter, an ionization chamber is much more tolerant to high-flux radiation when the 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 five locations (400–1900 m horizontally from the estimated location of the 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 the altitude of electron acceleration to be 2.5|$\pm$|0.5 km from 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}\times10^{18}$|⁠, which is approximately in the same range as the accelerated electrons estimated from space-based observations of upward TGFs (⁠|$4\times10^{16}$|–|$3\times10^{19}$| in Ref. [29]), while many TGFs observed in space are thought to originate at altitudes higher than 8 km [29,63]. 5.3. The end of gamma-ray glow from thundercloud One of key questions the GROWTH project set out to answer is how a stable electron-accelerating 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 closing phase of the life cycle is concerned, multiple previous measurements reported abrupt termination of gamma-ray glow that coincided with lightning discharge (see Refs. [18,45,64] and references therein). To reveal the precise relationship between lightning discharge and the cessation of gamma-ray glow, Wada et al. [50] analyzed an abrupt termination event that was observed in Ishikawa Prefecture on February 11 2017, by combining gamma-ray data collected by our detector and from the GODOT project [59], as well as LF data collected by multiple receivers located |${\sim} 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 radio signal measurements 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 gamma-ray glow termination and extended over a |${\sim} 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 rays about 15 km away from the gamma-ray observation site, causality in this event is obvious: the 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 (⁠|${\sim} 50$| km), we were unable to resolve the vertical structure of the discharge. Continued simultaneous observation in gamma-ray and radio frequencies is anticipated to shed light on the charge structure in the cloud in such events in the future. 6. Conclusion ◦ Aimed at multi-point observation of particle acceleration and high-energy gamma-ray emission of thundercloud and lightning, we launched a new experiment campaign called the “Thundercloud Project” in 2015, and developed a new, compact gamma-ray detector system (⁠|$35 \times 45 \times 20$| cm|$^3$| in size), each carrying a BGO or CsI scintillation crystal. ◦ We have deployed 15 detectors to four cities in Ishikawa Prefecture and Niigata Prefecture in Japan over 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 6 2017 in Niigata, allowed us to record the whole process of downward TGF followed by photonuclear reaction and a traveling positron-emitting isotope cloud. ◦ We have revealed that long-duration gamma-ray bursts can be abruptly terminated by the passage of a developing lightning leader (separated by 700 m horizontally) based on February 11 2017 data collected in Ishikawa Prefecture [48]. This is another stepping stone to understanding the life cycle of particle acceleration regions in a thundercloud. ◦ With accurate timing information from GPS, we have been able to correlate our gamma-ray data with radio frequency observations, enabling multi-messenger studies of high-energy activities of thundercloud and lightning. ◦ We will continue observation campaigns in coming winter seasons. Acknowledgements We deeply thank D. Yonetoku and T. Sawano (Kanazawa University), K. Watarai (Kanazawa University High School), K. Yoneguchi and Kanazawa Izumigaoka High School, K. Kimura (Komatsu High School), K. Kitano (Science Hills Komatsu), K. Kono (Ishikawa Plating Industry Co., Ltd), Kanazawa Nishi High School, Industrial Research Institute of Ishiwaka, Sakaida Fruits, Kanazawa Institute of Technology, Television Kanazawa Corporation, Ishikawa Prefectural University, Ishikawa Prefectural Institute of Public Health and Environmental Science, Kanazawa University Noto School, Norikura Observatory of Institute of Cosmic-Ray Research, The University of Tokyo, Non-Profit Organization Mount Fuji Research Station, Niigata Prefectural Radiation Monitoring Center, and the Radiation Safety Group of the Kashiwazaki-Kariwa nuclear power station, Tokyo Electric Power Company Holdings, Inc. for support of the detector deployment. We are grateful to H. Sakurai, M. Niikura, and the Sakurai group members at the Graduate School of Science, The University of Tokyo for providing the BGO scintillation crystals. T. Nakano, T. Tamagawa (RIKEN), A. Bamba, H. Odaka (The University of Tokyo), D. Umemoto (Kobe University), M. Sato and Y. Sato (Hokkaido University), H. Nanto, G. Okada (Kanazawa Institute of Technology), M. Kamogawa (University of Shizuoka), G. S. Bowers (Los Alamos National Laboratory), D. M. Smith (University of California Santa Cruz), T. Morimoto (Kindai University), Y. Nakamura (Kobe City College of Technology), A. Matsuki and M. Kubo (Kanazawa University), T. Ushio (Osaka University), and H. Sakai (University of Toyama) contributed to the project through detector development, data provisioning, and discussions on results. S. Otsuka (The University of Tokyo) and H. Kato (RIKEN) also supported detector developments. The Monte Carlo simulations were performed on the HOKUSAI GreatWave and BigWaterfall supercomputing systems operated by the RIKEN Advanced Center for Computing and Communication. The maps and aerial photos used in the figures are taken from the Geospatial Information Authority of Japan. Data of XRAIN are provided by the Ministry of Land, Infrastructure, Transport and Tourism via the Data Integration and Analysis System (DIAS), The University of Tokyo. This project is supported by RIKEN Special Postdoctoral Researchers Program, Japan Society for the Promotion of Science (JSPS)/MEXT KAKENHI grants (15K05115, 15H03653, 16H04055, 16H06006, 16K05555, 17K12966, 18J13355, 19H00683), the Hakubi project and SPIRITS 2017 of Kyoto University, and the joint research program of the Institute for Cosmic Ray Research (ICRR), The University of Tokyo. The bootstrapping phase of this project was supported by a crowdfunding campaign named Thundercloud Project on the academic crowdfunding platform “academist.” We thank Y. Shikano, Y. Araki, M. T. Hayashi, N. Matsumoto, T. Enoto, K. Hayashi, S. Koga, T. Hamaji, Y. Torisawa, S. Sawamura, J. Purser, S. Suehiro, S. Nakane, M. Konishi, H. Takami, T. Sawara, and all other supporters of the crowdfunding campaign, and Adachi design laboratory for their support of the creation of digital content and visuals. Footnotes 1https://github.com/growth-team/ 2http://godrb.com/ References [1] Wilson C. T. R. , Math. Proc. Camb. Phil. Soc. 22 , 534 ( 1925 ). ( http://dx.doi.org/10.1017/S0305004100003236 ) Crossref Search ADS [2] Gurevich A. V. , Milikh G. M., and Roussel-Dupre R., Phys. Lett. A 165 , 463 ( 1992 ). ( http://dx.doi.org/10.1016/0375-9601(92)90348-P ) Crossref Search ADS [3] Parks G. K. , Mauk B. H., Spiger R., and Chin J., Geophys. Res. Lett. 8 , 1176 ( 1981 ). ( http://dx.doi.org/10.1029/GL008i011p01176 ) Crossref Search ADS [4] McCarthy M. and Parks G. K., Geophys. Res. Lett. 12 , 393 ( 1985 ). ( http://dx.doi.org/10.1029/GL012i006p00393 ) Crossref Search ADS [5] Kelley N. A. , Smith D. M., Dwyer J. R., Splitt M., Lazarus S., Martinez-McKinney F., Hazelton B., Grefenstette B., Lowell A., and Rassoul H. K., Nature Commun. 6 , 7845 ( 2015 ). ( http://dx.doi.org/10.1038/ncomms8845 ) Crossref Search ADS [6] Kochkin P. , van Deursen A. P. J., Marisaldi M., Ursi A., de Boer A. I., Bardet M., Allasia C., Boissin J.-F., Flourens F., and Østgaard N., J. Geophys. Res. Atmos. 122 , 12801 ( 2017 ). ( https://doi.org/10.1002/2017JD027405 ) Crossref Search ADS PubMed [7] Østgaard N. et al. , J. Geophys. Res. Atmos. 124 , 7236 ( 2019 ). ( https://doi.org/10.1029/2019JD030312 ) Crossref Search ADS PubMed [8] Eack K. B. , Beasley W. H., Rust W. D., Marshall T. C., and Stolzenburg M., J. Geophys. Res. Atmos. 101 , 29637 ( 1996 ). ( http://dx.doi.org/10.1029/96JD01705 ) Crossref Search ADS [9] Chubenko A. P. , Antonova V. P., Kryukov S. Yu., Piskal V. V., Ptitsyn M. O., Shepetov A. L., Vildanova L. I., Zybin K. P., and Gurevich A. V., Phys. Lett. A 275 , 90 ( 2000 ). ( http://dx.doi.org/10.1016/S0375-9601(00)00502-8 ) Crossref Search ADS [10] Alexeenko V. V. , Khaerdinov N. S., Lidvansky A. S., and Petkov V. B., Phys. Lett. A 301 , 299 ( 2002 ). ( http://dx.doi.org/10.1016/S0375-9601(02)00981-7 ) Crossref Search ADS [11] Tsuchiya H. et al. , Phys. Rev. Lett. 102 , 255003 ( 2009 ). ( http://dx.doi.org/10.1103/PhysRevLett.102.255003 ) Crossref Search ADS PubMed [12] Tsuchiya H. et al. , Phys. Rev. D 85 , 092006 ( 2012 ). ( https://doi.org/10.1103/PhysRevD.85.092006 ) Crossref Search ADS [13] Torii T. , Sugita T., Tanabe S., Kimura Y., Kamogawa M., Yajima K., and Yasuda H., Geophys. Res. Lett. 36 , L13804 ( 2009 ). ( http://dx.doi.org/10.1029/2008GL037105 ) Crossref Search ADS [14] Bowers G. S. et al. , Phys. Rev. D 100 , 043021 ( 2019 ). ( http://dx.doi.org/10.1103/PhysRevD.100.043021 ) Crossref Search ADS [15] Chilingarian A. , Daryan A., Arakelyan K., Hovhannisyan A., Mailyan B., Melkumyan L., Hovsepyan G., Chilingaryan S., Reymers A., and Vanyan L., Phys. Rev. D 82 , 043009 ( 2010 ). ( http://dx.doi.org/10.1103/PhysRevD.82.043009 ) Crossref Search ADS [16] Chilingarian A. and Mkrtchyan H., Phys. Rev. D 86 , 072003 ( 2012 ). ( http://dx.doi.org/10.1103/PhysRevD.86.072003 ) Crossref Search ADS [17] Chilingarian A. , Mkrtchyan H., Karapetyan G., Chilingaryan S., Sargsyan B., and Arestakesyan A., Sci. Rep. 9 , 6253 ( 2019 ). ( http://dx.doi.org/10.1038/s41598-019-42786-7 ) Crossref Search ADS PubMed [18] Chilingarian A. , Khanikyants Y., Mareev E., Pokhsraryan D., Rakov V. A., and Soghomonyan S., J. Geophys. Res. Atmos. 122 , 7582 ( 2017 ). ( http://dx.doi.org/10.1002/2017JD026744 ) Crossref Search ADS [19] Chilingarian A. , Soghomonyan S., Khanikyanc Y., and Pokhsraryan D., Astropart. Phys. 105 , 54 ( 2019 ). ( http://dx.doi.org/10.1016/j.astropartphys.2018.10.004 ) Crossref Search ADS [20] Kitagawa N. and Michimoto K., J. Geophys. Res. 99 , 10713 ( 1994 ). ( http://dx.doi.org/10.1029/94JD00288 ) Crossref Search ADS [21] Torii T. , Takeishi M., and Hosono T., J. Geophys. Res. Atmos . 107 , ACL 2-1 ( 2002 ). ( http://dx.doi.org/10.1029/2001JD000938 ) Crossref Search ADS [22] Torii T. , Sugita T., Kamogawa M., Watanabe Y., and Kusunoki K., Geophys. Res. Lett. 38 , L24801 ( 2011 ). ( http://dx.doi.org/10.1029/2011GL049731 ) Crossref Search ADS [23] Smith D. M. , Lopez L. I., Lin R. P., and Barrington-Leigh C. P., Science 307 , 1085 ( 2005 ). ( http://dx.doi.org/10.1126/science.1107466 ) Crossref Search ADS PubMed [24] Tavani M. et al. [AGILE Team], Phys. Rev. Lett. 106 , 018501 ( 2011 ). ( http://dx.doi.org/10.1103/PhysRevLett.106.018501 ) Crossref Search ADS PubMed [25] Foley S. et al. , J. Geophys. Res. Space Phys. 119 , 5931 ( 2014 ). ( http://dx.doi.org/10.1002/2014JA019805 ) Crossref Search ADS [26] Fishman G. J. et al. , Science 264 , 1313 ( 1994 ). ( http://dx.doi.org/10.1126/science.264.5163.1313 ) Crossref Search ADS PubMed [27] Marisaldi M. et al. , J. Geophys. Res. Space Phys. 115 , A00E13 ( 2010 ). ( https://doi.org/10.1029/2009JA014502 ) Crossref Search ADS [28] Briggs M. S. et al. , J. Geophys. Res. Space Phys. 115 , A07323 ( 2010 ). ( https://doi.org/10.1029/2009JA015242 ) Crossref Search ADS [29] Mailyan B. G. , Briggs M. S., Cramer E. S., Fitzpatrick G., Roberts O. J., Stanbro M., Connaughton V., McBreen S., Bhat P. N., and Dwyer J. R., J. Geophys. Res. Space Phys. 121 , 346 ( 2016 ). ( http://dx.doi.org/10.1002/2016JA022702 ) Crossref Search ADS [30] Neubert T. et al. , Science 367 , 183 ( 2019 ). ( http://dx.doi.org/10.1126/science.aax3872 ) PubMed [31] Østgaard N. et al. , J. Geophys. Res. Atmos., 124 , 14024 ( 2019 ). ( https://doi.org/10.1029/2019JD031214 ) Crossref Search ADS [32] Dwyer J. R. and Smith D. M., Geophys. Res. Lett. 32 , L22804 ( 2005 ). ( https://doi.org/10.1029/2005GL023848 ) [33] Celestin S. and Pasko V. P., J. Geophys. Res. Space Phys. 116 , A03315 ( 2011 ). ( http://dx.doi.org/10.1029/2010JA016260 ) Crossref Search ADS [34] Celestin S. , Xu W., and Pasko V. P., J. Geophys. Res. Space Phys. 120 , 10712 ( 2015 ). ( http://dx.doi.org/10.1002/2015JA021410 ) Crossref Search ADS [35] Moss G. D. , Pasko V. P., Liu N., and Veronis G., J. Geophys. Res. Space Phys. 111 , A02307 ( 2006 ). ( http://dx.doi.org/10.1029/2005JA011350 ) Crossref Search ADS [36] Dwyer J. R. , Phys. Plasmas 14 , 042901 ( 2007 ). ( http://dx.doi.org/10.1063/1.2709652 ) Crossref Search ADS [37] Dwyer J. R. , J. Geophys. Res. Space Phys. 117 , A02308 ( 2012 ). ( http://dx.doi.org/10.1029/2011JA017160 ) [38] Dwyer J. R. et al. , Geophys. Res. Lett. 31 , L05119 ( 2004 ). ( https://doi.org/10.1029/2003GL018771 ) [39] Hare B. M. et al. , J. Geophys. Res. Atmos. 121 , 6511 ( 2016 ). ( http://dx.doi.org/10.1002/2015JD024426 ) Crossref Search ADS [40] Tran M. D. , Rakov V. A., Mallick S., Dwyer J. R., Nag A., and Heckman S., J. Atmos. Sol.-Terr, Phys. 136 , 86 ( 2015 ). ( http://dx.doi.org/10.1016/j.jastp.2015.10.010 ) [41] Abbasi R. U. et al. , J. Geophys. Res. Atmos. 123 , 6864 ( 2018 ). ( http://dx.doi.org/10.1029/2017JD027931 ) Crossref Search ADS [42] Ringuette R. , Case G. L., Cherry M. L., Granger D., Guzik T. G., Stewart M., and Wefel J. P., J. Geophys. Res. Space Phys. 118 , 7841 ( 2013 ). ( http://dx.doi.org/10.1002/2013JA019112 ) Crossref Search ADS PubMed [43] Tsuchiya H. et al. , Phys. Rev. Lett. 99 , 165002 ( 2007 ). ( https://doi.org/10.1103/PhysRevLett.99.165002 ) Crossref Search ADS PubMed [44] Tsuchiya H. et al. , J. Geophys. Res. 116 , D09113 ( 2011 ). ( https://doi.org/10.1029/2010JD015161 ) [45] Tsuchiya H. et al. , Phys. Rev. Lett. 111 , 015001 ( 2013 ). ( http://dx.doi.org/10.1103/PhysRevLett.111.015001 ) Crossref Search ADS PubMed [46] Umemoto D. et al. , Phys. Rev. E 93 , 021201(R) ( 2016 ). ( https://doi.org/10.1103/PhysRevE.93.021201 ) Crossref Search ADS [47] Enoto T. et al. , Nature 551 , 481 ( 2017 ). ( http://dx.doi.org/10.1038/nature24630 ) Crossref Search ADS PubMed [48] Wada Y. et al. , Geophys. Res. Lett. 45 , 5700 ( 2018 ). ( http://dx.doi.org/10.1029/2018GL077784 ) Crossref Search ADS [49] Wada Y. et al. , Commun. Phys. 2 , 67 ( 2019 ). ( http://dx.doi.org/10.1038/s42005-019-0168-y ) Crossref Search ADS [50] Wada Y. , Enoto T., Nakazawa K., Furuta Y., Yuasa T., Nakamura Y., Morimoto T., Matsumoto T., Makishima K., and Tsuchiya H., Phys. Rev. Lett. 123 , 061103 ( 2019 ). ( http://dx.doi.org/10.1103/PhysRevLett.123.061103 ) Crossref Search ADS PubMed [51] Wada Y. et al. , J. Geophys. Res. Atmos. 125 , e2019JD031730 ( 2020 ). ( https://doi.org/10.1029/2019JD031730 ) Crossref Search ADS [52] Brun R. and Rademakers F., Nucl. Instr. Meth. Phys. Res. A 389 , 81 ( 1997 ). ( http://dx.doi.org/10.1016/S0168-9002(97)00048-X ) Crossref Search ADS [53] Wells D. C. , Greisen E. W., and Harten R. H., A&AS 44 , 363 ( 1981 ). ( https://ui.adsabs.harvard.edu/abs/1981A%26AS...44..363W/abstract ) [54] Agostinelli S. et al. , Nucl. Instr. Meth. Phys. Res, A 506 , 250 ( 2003 ). ( https://doi.org/10.1016/S0168-9002(03)01368-8 ) [55] Allison J. et al. , IEEE Trans. Nucl. Sci. 53 , 270 ( 2006 ). ( http://dx.doi.org/10.1109/TNS.2006.869826 ) Crossref Search ADS [56] Allison J. et al. , Nucl. Instr. Meth. Phys. Res. A 835 , 186 ( 2016 ). ( http://dx.doi.org/10.1016/j.nima.2016.06.125 ) Crossref Search ADS [57] Gurevich A. V. et al. , Phys. Rev. Lett. 108 , 125001 ( 2012 ). ( http://dx.doi.org/10.1103/PhysRevLett.108.125001 ) Crossref Search ADS PubMed [58] Chilingarian A. , Bostanjyan N., and Vanyan L., Phys. Rev. D 85 , 085017 ( 2012 ). ( http://dx.doi.org/10.1103/PhysRevD.85.085017 ) Crossref Search ADS [59] Bowers G. S. , Smith D. M., Martinez-McKinney G. F., Kamogawa M., Cummer S. A., Dwyer J. R., Wang D., Stock M., and Kawasaki Z., Geophys. Res. Lett. 44 , 10063 ( 2017 ). ( http://dx.doi.org/10.1002/2017GL075071 ) Crossref Search ADS [60] Babich L. P. , JETP Lett. 84 , 285 ( 2006 ). ( http://dx.doi.org/10.1134/S0021364006180020 ) Crossref Search ADS [61] Babich L. P. , Geomag. Aeron. 47 , 664 ( 2007 ). ( http://dx.doi.org/10.1134/S0016793207050155 ) Crossref Search ADS [62] Carlson B. E. , Lehtinen N. G., and Inan U. S., J. Geophys. Res. Space Phys. 115 , A10324 ( 2010 ). ( https://doi.org/10.1029/2010JA015647 ) [63] Cummer S. A. , Lyu F., Briggs M. S., Fitzpatrick G., Roberts O. J., and Dwyer J. R., Geophys. Res. Lett. 42 , 7792 ( 2015 ). ( http://dx.doi.org/10.1002/2015GL065228 ) Crossref Search ADS [64] Chilingarian A. , Khanikyants Y., Rakov V. A., and Soghomonyan S., Atmos. Res. 233 , 104713 ( 2020 ). ( http://dx.doi.org/10.1016/j.atmosres.2019.104713 ) Crossref Search ADS Author notes Currently, an individual researcher. © The Author(s) 2020. Published by Oxford University Press on behalf of the Physical Society of Japan. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. TI - Thundercloud Project: Exploring high-energy phenomena in thundercloud and lightning JF - Progress of Theoretical and Experimental Physics DO - 10.1093/ptep/ptaa115 DA - 2020-10-09 UR - https://www.deepdyve.com/lp/oxford-university-press/thundercloud-project-exploring-high-energy-phenomena-in-thundercloud-wGiW37Gzz8 VL - 2020 IS - 10 DP - DeepDyve ER -