Gamma-ray glow preceding downward terrestrial gamma-ray flash

Two types of high-energy events have been detected from thunderstorms. One is “terrestrial gamma-ray flashes” (TGFs), sub-millisecond emissions coinciding with lightning discharges. The other is minute-lasting “gamma-ray glows”. Although both phenomena are thought to originate from relativistic runaway electron avalanches in strong electric fields, the connection between them is not well understood. Here we report unequivocal simultaneous detection of a gamma-ray glow termination and a downward TGF, observed from the ground. During a winter thunderstorm in Japan on 9 January 2018, our detectors caught a gamma-ray glow, which moved for ~100 s with ambient wind, and then abruptly ceased with a lightning discharge. Simultaneously, the detectors observed photonuclear reactions triggered by a downward TGF, whose radio pulse was located within ~1 km from where the glow ceased. It is suggested that the highly-electrified region producing the glow was related to the initiation of the downward TGF. Thunderstorms are thought to produce two types of high-energy emissions, terrestrial gamma-ray flashes and gamma-ray glows however due to the difficulty in their observation the exact relation between the two is still not well-understood. Here, the authors report the simultaneous detection of a gamma-ray glow and a downward terrestrial gamma-ray flash suggesting the origin of the two phenomena are related.

S ince McCarthy and Parks 1 found radiation-dose enhancements inside thunderclouds with an airborne detector in 1980s, high-energy phenomena associated with thunderstorms have been detected inside the Earth's atmosphere and from space. Terrestrial gamma-ray flashes (TGFs) are burst-like emission with their photon energy extending up to 20 MeV that last for several hundred microseconds, coincident with lightning discharges. They were first detected from space by Compton Gamma-Ray Observatory 2 , and since then have been reported by many other satellites [3][4][5][6][7][8] . Similar phenomena but going downward have been found in recent years at ground level [9][10][11][12][13][14][15][16][17] . They, now called "downward TGFs", share several features with TGFs observed from space, such as coincidence with lightning, submillisecond durations, and energy spectra extending to >10 MeV. Downward TGFs that contains enough photons above 10 MeV have been experimentally shown to trigger atmospheric photonuclear reactions, namely producing neutrons and positronemitting radioactive nuclei 13,14 . These photoneutrons can be observed as a short-duration gamma-ray burst lasting for several hundreds of milliseconds, as they are absorbed by atmospheric nuclei via neutron-capture processes 14,18 .
Although TGFs and gamma-ray glows are distinguished clearly by duration, brightness, and timing with regard to lightning discharges, both of them are thought to originate from a common fundamental mechanism, called relativistic runaway electron avalanches (RREAs 39,40 ). According to Wilson's hypothesis 41 , seed electrons (provided by, e.g., cosmic rays) can be accelerated up to an energy of tens of MeVs in strong electric fields, producing secondary electrons. The number of multiplied and accelerated electrons exponentially increases, and the accelerated electrons finally emit bremsstrahlung gamma rays as they interact with ambient atmospheric nuclei. Dwyer 42 proposed additional electron-seeding processes by positrons and backscattered gamma rays into the RREA mechanism, called "relativistic feedback model". This model can achieve a higher multiplication factor than that of a RREA alone, and thus are thought to explain extraordinarily high brightness of TGFs.
Despite an increasing number of respective observation samples of TGFs and gamma-ray glows, connections between them remain poorly understood. This is primarily because there has been no report of simultaneous detection of both, except for a very recent short report on a marginal detection 17 . In this paper, we report the first unequivocal simultaneous detection of them at sea level and discuss its implications.

Results
Observation of high-energy phenomena in winter thunderstorms. The Gamma-ray Observation of Winter Thunderclouds (GROWTH) collaboration 31,32,35,43 has been engaged with a multi-point observation campaign of atmospheric high-energy phenomena in coastal areas of Japan Sea 14,44 . Winter thunderstorms in Japan are ideal targets to observe this type of phenomena due to their unique characteristics; most notably typical altitude of clouds is significantly lower than ordinary 38,45,46 , which makes sea-level observations of gamma-ray glows viable.
We have developed portable radiation detectors dedicated to the multi-point observation. They have a 25 cm × 8 cm × 2.5 cm Bi 4 Ge 3 O 12 (BGO) scintillation crystal coupled with two photomultiplier tubes (PMTs; HAMAMATSU R1924A). Outputs from the PMTs are amplified, and then read out by a 50 MHz digitiser onboard a data acquisition system. The data acquisition system stores 20-μs waveforms of the amplified analogue outputs once a pulse is detected, and extracts the maximum and minimum value as well as the timing of the pulse (see also Detector calibration). The maximum value corresponds to the energies of the pulse, and the minimum the analogue baseline voltage. The data acquisition system also records counts of discarded photon events due to buffer overflow, which are used for dead-time correction. Three detectors were deployed at three observation sites in Kanazawa City, the capital of Ishikawa Prefecture, by the Japan Sea coast (Fig. 1) and have been operated since October 2016.
Lightning discharges were monitored by a broadband lowfrequency (LF: 0.8-500 kHz) lightning mapping network (hereafter LF network), for which detectors are installed along Toyama Bay and in Noto peninsula. Another receiver in the extreme-lowfrequency band (ELF: 1-100 Hz) is installed at Kuju, as summarised in the section Radio observations. We also utilise lightning location data of Japanese Lightning Detection Network (JLDN) operated by Franklin Japan Co., Ltd.
Detection of gamma-ray glow and downward TGF. On 9 January 2018, two of our detectors shown in Fig. 1 recorded gamma-ray glows. Figure 2a,  An snapshot image of the X-band radar network at 17:55 shows a heavy precipitation area, corresponding to a thundercloud, located between detectors A and B (Fig. 1a). The radar data suggest that the thundercloud passed over the two detectors towards east-northeast with a speed of 19.3 ± 1.4 m s −1 (see Wind estimation with X-band radar). Since the temporal separation between the glow detection by the two detectors is consistent with the time for the thundercloud to travel the distance between the two detectors, we consider that the gamma-ray glows recorded by the two detectors are from the same cloud and hence of the same origin.
At the same time as the glow termination and the lightning discharge, both detectors A and B recorded a short-duration radiation burst lasting for~200 ms simultaneously. The countrate profiles of the 200-ms-lasting short burst shown in Fig. 2c, d exhibit a steep rise and decay with time constants of 52.0 ± 4.9 and 59.2 ± 1.7 ms for detectors A and B, respectively. Combining the timing analysis with spectral analysis (see Gamma-ray emission originating from neutrons), the short burst is found to originate from neutron captures by atmospheric nitrogen nuclei, which Rutjes et al. 18 predicted as "TGF afterglow", and Enoto et al. 14 observationally demonstrated. In addition, detector B recorded a faint annihilation emission at 511 keV for 10 s after the short burst (see Positron production by beta-plus decay). These features imply that atmospheric photonuclear reactions such as 14 N + γ→ 13 N + n and 16 O + γ→ 15 O + n took place coincident with the lightning discharge, as discussed in Bowers et al. 13 and Enoto et al. 14 Fig. 1 Maps of observation sites in Kanazawa. a A precipitation map at 17:55 on 9 January 2018 in coordinated universal time, obtained by the extended radar information network. Orange circles show radiation monitors, the white arrow the wind speed and its direction, and the dashed-line rectangle the region of panel b. b A panoramic photograph around detectors A and B marked by filled orange circles. The lightning locations determined by Japanese Lightning Detection Network and our low-frequency network are marked by red triangles and white stars, respectively. The numbers beside the white stars correspond to the order of the low-frequency pulses. A typical locating error circle (~300 m in radius) is overlaid for the fifth pulse. The magenta dashed line shows a track of the centre of the gamma-ray glow, which ceased at the magenta cross marker. The shaded blue circle shows the region where 30-sbinned count rates in the 3.0-20.0 MeV band reached twice or more higher than the background rates when the glow disappeared, estimated with the simulation (see Simulation of gamma-ray glow). This corresponds to a detection significance of 8.1σ or higher. The radius of the region is 0.88 km. The background photographs are provided by the Geospatial Information Authority of Japan Figure 3a, b shows the maximum and minimum waveform values of photon events during the short burst recorded by detectors A and B, respectively. At the very beginning of the short burst, both detectors A and B recorded saturated pulses (the maximum values exceeding >4 V), and then significant negative values of the baseline (the minimum values) called "undershoot" for~10 ms. Although detector B failed to acquire the main part of the undershoot due to buffer overflow in the data acquisition system, it recorded the saturated pulses and the last part of the undershoot. As demonstrated in Methods: Initial flash of Enoto et al. 14 , this feature manifests the existence of an extremely large energy deposit (much more than hundreds of MeVs) in the scintillation crystal within a few milliseconds, which is a clear sign of a downward TGF. In the following analysis we employ an elapsed time t from the onset of the downward TGF at 17:54:50.308892 UTC, recorded by detector B.
The LF network recorded a consecutive series of waveforms of the lightning discharge lasting for~400 ms (Fig. 3c). The downward TGF coincided with a large-amplitude pulse at the initial phase of the lightning discharge within 10 μs (Fig. 3d). We detected four or so precursory pulses shortly before the largeamplitude pulse. No pulses had been detected before the precursory pulses by the LF network. The ELF measurement also confirmed that the associated ELF pulse was coming from the LF source. In addition, JLDN also reported a negative intracloud/intercloud (IC) discharge of −197 kA at t = −13 μs, which is temporally associated with the large-amplitude pulse. Figure 1b shows the source positions of the large-amplitude and precursory pulses determined by the LF network. At the beginning, the small precursory pulses took place in a southwest region less than 3 km away from detector B. Then, the main large-amplitude pulse (the fifth one in Figs. 1b and 3c) occurred 0.6 km southwest of detector B at t = −5.5 μs. JLDN also located the large-amplitude pulse within 0.9 km from detector B. These temporal and spatial correlations lead us to conclude that the large-amplitude LF pulse is associated with the downward TGF.
Production mechanism of gamma-ray glow. The multi-point observation enables us to investigate characteristics of the gamma-ray glow preceding the lightning initiation and the downward TGF. First, we perform spectral analysis. Figure 4 shows the background-subtracted gamma-ray energy spectra, extracted from −69 s < t < −39 s and −30 s < t < −10 s for detectors A and B, respectively. The detector response function is calculated with the GEANT4 Monte Carlo simulation framework 47 , and is convolved with a model spectrum in spectral fitting using the XSPEC package 48 . The observed spectra, of which instrumental responses are corrected, are found to be well explained by an empirical power-law function with an exponential cutoff, ε −Γ exp[(−ε/ε cut ) α ], where ε, Γ, ε cut , and α are the photon energy (MeV), power-law photon index, cutoff energy (MeV), and cutoff index, respectively. The best-fitting parameters are Γ ¼ 0:90 þ0:06 À0:08 and 1:02 þ0:04 À0:05 , ε cut ¼ 6:4 þ1:0 À1:1 and 8:5 þ0:8 À0:9 MeV,  À0:5 10 À5 and 2:4 þ0:7 À0:6 10 À5 ergs cm −2 s −1 on average over −69 s < t < −39 s and −30 s < t < −10 s integration periods for detectors A and B, respectively. Here and after, all the errors are statistical at 1σ confidence level, unless otherwise mentioned.
We then perform another set of Monte Carlo simulations, using GEANT4, and compare the obtained energy spectra and count-rate histories with the simulated ones to investigate atmospheric interactions and propagation of electrons and gamma rays (see Simulation of gamma-ray glow). We find a model of spatial and energy spectral distribution for avalanche electrons in the RREA region which can reproduce both the obtained gamma-ray spectra and count-rate histories, and summarise the results in Figs. 1 and 4. The best-fit value of the RREA terminus altitude h base is 400 m, which means the electron avalanche took place in the lower part of the winter thundercloud, and the offsets from the centre of the RREA region are 540 and 80 m for detectors A and B, respectively. The electron flux distribution is consistent with being proportional to a function of a distance from the RREA centre r, exp(−r/150 m), providing the circularly symmetric distribution. Figure 1b shows the centre position of the RREA region at the moment of the termination. Normalising the simulation result, we estimate the total production rate of 1-50 MeV avalanche electrons to be 3.66 × 10 12 electrons s −1 . The electron flux F(r, ε) at the terminus of RREA is also estimated to be a function of r and ε This model reproduces the observed count-rate histories and spectra, except the increase in the count rate of detector B during −5 s < t < 0 s. This period is discussed in the section Abrupt increase in count rates of gamma-ray glow before downward TGF.
Let us consider the electron multiplication factor M = F RREA / F seed , where F RREA and F seed are the average electron flux at the RREA terminus and seed electron flux, respectively. Integrating Eq. (1) yields the 0.3-50 MeV average flux within r = 150 m of F RREA = 7.5 × 10 2 electrons cm −2 s −1 . Assuming that the seed electrons are mainly produced by cosmic rays, the 0.3-50 MeV seed electron flux is a function of a vertical acceleration length L and h base given by (see Seed electrons). The multiplication factor M is thus a function of L, with the fixed h base (400 m).
In the RREA region, electron flux is known to increase 39 exponentially as a function of L, F RREA = F seed exp(L/λ), assuming that change of the vertical atmospheric pressure is negligible for the RREA processes at the low altitude. The avalanche length λ is empirically determined (see ref. 49  As Dwyer 50 pointed out, the multiplication factor would not exceed~10 5 in the RREA-only case because thunderclouds cannot maintain an acceleration length required for it. Given that L can reach twice as high as the typical diameter of the RREA region 50 , L < 600 m is required in this case, where the typical radius r = 150 m is employed. The 0.3 MV m −1 case is not plausible because the required acceleration length L = 3240 m cannot be maintained inside the thundercloud. In the other cases, it is necessary to take into account the relativistic feedback processes to explain the estimated avalanche multiplication factor. The relativistic feedback processes are parameterised with a feedback factor γ, the fraction of the seed electrons provided by the steady-state relativistic feedback processes 50 . The flux of runaway electrons is then modified as F RREA = F seed (L) exp(L/λ)/ (1 − γ). Figure 5 shows this relation between L and γ to explain the observed flux at the RREA terminus. To satisfy the condition L < 600 m, γ should be larger than 0.998 and 0.846 for 0.35 and 0.4 MV m −1 , respectively. This suggests that the number of feedback-origin seed electrons is higher than that of cosmic-ray seed electrons by a factor of >5.5 for our event.
Abrupt increase in count rates of gamma-ray glow before TGF. The count-rate history of detector B exhibited an additional increase during −5 < t < −0 s (Fig. 6a). Figure 6b shows the ratio of the simulated model to the observed history. Although the observed history is well reproduced by the simulation up to t = −5 s, the observed count rate is twice as high as the simulation in −5 s < t < −0 s. Figure 6c shows the three energy spectra extracted from the time regions of −10 s < t < −5 s, −5 s < t < −2 s, and −2 s < t < 0 s. All the spectra show a power-law function with an exponential cutoff, indicating that bremsstrahlung is still the main process of gamma-ray production. Since our simulations fail to reproduce this increase in count-rate, we speculate that the increase was caused by a fluctuation of the intrinsic electron fluxes, rather than by the movement of the RREA region with the ambient wind flow.
Based on the working hypothesis of the speculated increase of the accelerated electron flux, at least one of the following is required to have taken place: (1) stronger electric fields of the RREA region, (2) longer acceleration length, and/or (3) increase in the feedback factor γ. However, since lightning did not occur during this period (−5 s < t < 0 s), atmospheric mechanism could not drastically change the meteorological conditions, such as electric fields and acceleration length, within 5 s. We thus The RREA and relativistic feedback processes remained stable until t = −5 s; this state corresponds to the "steady state" of relativistic feedback as defined by Dwyer 50 , namely γ < 1. In general, when γ exceeds 1, an electron flux would spontaneously increase, and an RREA region should collapse. The timescale of the flux increase depends on the types of the relativistic feedback processes. The feedback process by positrons can discharge RREA regions within microseconds 50 . This timescale is close to that of TGFs, and is much shorter than that of the observed abrupt increase (i.e. 5 s). Alternatively, the feedback by backscattered X-rays may trigger a second-order discharge in RREA regions 50 . At present, even though the 5-s abrupt flux rise seems to be of great importance, its origin is yet to be understood.

Discussion
To conclude the relation between the gamma-ray glow and the downward TGF, verifying their temporal and positional coincidence will give a strong clue. Our observation cannot clarify whether the glow termination or the downward TGF took place first because these phenomena seemed to be slightly overlapped. On the other hand, the positional coincidence of the gamma-ray glow and the downward TGF in the present case is precisely determined owing to the multiple gamma-ray detectors and the LF network. The discussion in the section Production mechanism of gamma-ray glow suggests that the gamma-ray glow ceased when the source cloud was moving 130 m southwest of detector B (Fig. 1b). Also, the TGF-associated LF pulse was located within 0.5 km from detector B. Therefore, it is clear that the two phenomena are physically related to one another.
Our interpretation of the observed gamma-ray glow suggests that the electron acceleration site should have electric fields of 0.35 MV m −1 or higher in order to achieve the high electron multiplication factor of >10 5 with a plausible acceleration length. In such highly electrified regions, TGFs are thought to initiate more easily than in other less-electrified regions as Smith et al. 17 suggested. From another point of view, we speculate that the avalanche electrons of the gamma-ray glow can behave as seed electrons of the downward TGF. At the point where the TGFassociated LF pulse was located (point 5 in Fig. 1b), the 0.3-50 MeV electron flux at 400 m altitude is estimated to be 1.7 × 10 2 electrons cm −2 s −1 . By comparing this flux with that of the cosmic-ray-induced seed electrons (the canonical seed electron source), it is suggested that the highly-electrified region responsible for the gamma-ray glow can be the dominant source of seed electrons for the TGF which occurs in the close proximity of the gamma-ray glow. In addition, the abrupt count-rate increase monitored by detector B before the TGF (see section Abrupt increase in count rates of gamma-ray glow before downward TGF) suggests additional production of avalanche electrons for the gamma-ray glow, and might have predicted drastic changes in the electrified region such as the lightning discharge and the TGF.
In the present high-energy event, the discussion above suggests a possibility that the high electron current in the gamma-ray glow assisted the initiation of the downward TGF. However, it still remains observationally unclear how gamma-ray glows and TGFs are related with each other in general. Among an increasing sample of glow terminations, TGF-associated events are still quite rare, i.e. only Smith et al. 17 and the present event. For example, a termination event during a winter thunderstorm in 2017 (ref. 38 ) was associated with an intracloud/intercloud discharge but not related with any signals for TGF-like emissions. As another example, a TGF-like intensive emission associated with photonuclear reactions was reported 14 , where no gamma-ray glows were recorded before the event. In these cases, we lack sufficient evidences due to our present sparse observation sites on the ground to conclude that glow terminations are not always associated with TGFs. Our future gamma-ray monitoring network combined with radio-frequency lightning mapping systems will give a clue to reveal the relation between TGFs and gammaray glows.
In summary, we detected a gamma-ray glow, terminated with a downward TGF which triggered atmospheric photonuclear reactions. The gamma-ray glow was so bright that the relativistic feedback processes are required. Although we cannot determine whether the glow termination or the downward TGF occurred  Fig. 6 Increase in count rates of detector B. a Comparison between the count-rate history of detector B and the exponential model (red dashed line). Background was subtracted. b Ratio of the Gaussian model to the observed count-rate history. c Background-subtracted spectra of detector B extracted from (black) −10 s < t < −5 s, (red) −5 s < t < −2 s, and (blue) −2 s < t < 0 s. The error of each data point is at 1σ confidence level. The spectra in −10 s < t < −5 s and −2 s < t < 0 s are shifted by a factor of 0.1 and 10, respectively, to avoid overlapping. Dead-time and detector responses are corrected first, the two high-energy phenomena in the atmosphere took place in an identical electrified region of a winter thundercloud, and hence are clearly related to each other in the present case.

Methods
Detector calibration. Energy calibration of the detectors was performed to convert the maximum value of a pulse into photon energy. We measured the centre of environmental background lines of 40 K (1.46 MeV) and 208 Tl (2.61 MeV), and built a linear calibration function which is utilised to assign the energy of each photon. All the detectors record 0.4-20.0 MeV gamma rays. See also Instrumental calibration in Enoto et al. 14

for details.
Absolute timing is conditioned by pulse-per-second signals of the Global Positioning System (GPS). The timing-assignment logic employed from 2017 to 2018 winter provides absolute timing accuracy of each photon better than 1 μs. However, detector A failed to receive the GPS signals during the experiment. Instead, we performed the calibration of detector A, using the internal clock time with~1 s accuracy, and then corrected the absolute timing so that the detection time of the downward TGF matches that with detector B.
Wind estimation with X-band radar. We utilised data of eXtended RAdar Information Network (XRAIN). XRAIN is a polarimetric weather radar network in the X band and has a spatial resolution of 280 m (east-west) × 230 m (north-south) mesh. It records two-dimensional precipitation maps with a 1-min interval. XRAIN also obtains three-dimensional maps of radar echoes and particle types with a 5min interval by the constant-altitude plan position indicator technique. However, the three-dimensional data are not utilised in the present paper because the XRAIN observations have a moderate spatial resolution of altitude (≥1 km), which is insufficient to discuss charge structures in the thundercloud.
Wind velocity and direction are estimated by overlaying and shifting precipitation maps at different time. First, 11 maps from 17:50 to 18:00 were extracted in the range of 36.4°N-36.7°N, 136.4°E-136.8°E. We then took a pair of maps with a 5-min interval (six pairs in total), and calculated the sum of precipitation residual at each mesh, given by Σ i;j ðP 1 ij À P 2 ij Þ 2 , where P 1 ij and P 2 ij are precipitation at each mesh on each map, and i and j are mesh indexes. With trial shifting of one map with several steps of the spatial resolution for four directions, we searched for the position which takes the minimum residual sum. The distance and direction for which the cloud moved in 5 min can be estimated from the amount of the map shift at the point of the minimum residual sum. Consequently, the wind direction and velocity at the moment of the glow detection were determined to be west-northwestwards and 19.3 ± 0.9 (systematic) ± 1.1 (statistic) m s −1 , respectively. Here, the quoted statistical error was calculated from the the standard deviation (1σ) of six pairs. The systematic error was determined by the mesh size and temporal interval of the map pair. The wind velocity with the overall error is then calculated to be 19.3 ± 1.4 m s −1 , where the standard error propagation in quadrature between the systematic and statistical errors is assumed to hold. Since the statistic error is smaller than 10% and is comparable with the systematic error, it is reasonable to assume that the wind parameters did not change considerably during the glow observation.
Gamma-ray emission originating from neutrons. Photonuclear reactions such as 14 N + γ→ 13 N + n and 16 O + γ→ 15 O + n expel~10 MeV neutrons from atmospheric nitrogen and oxygen nuclei [51][52][53] . The photoneutrons gradually lose their kinetic energy via elastic scatterings, and are eventually captured by atmospheric nuclei such as 14 N. In the dominant reaction 14 N + n → 15 N + γ, 15 N nuclei in excited states emit various de-excitation gamma-ray lines up to 10.8 MeV. In addition, de-excitation gamma rays from other nuclei such as Si and Al should be also emitted when photoneutrons were captured by ambient nuclei in soil, buildings, and components of the detectors. These de-excitation gamma rays originating from neutron captures are thought to compose the short burst 14,18 .
The timescale of the short burst is determined by neutron thermalisation 13,14,18 . A numerical calculation predicts the neutron-capturing rate of exp(−t/τ) for 5 ms < t < 120 ms, where t is the elapsed time from the onset of the TGF and τ ≈ 56 ms is the decay constant 14 . The count-rate histories of the observed burst have decay constants of 52.0 ± 4.9 and 59.2 ± 1.7 ms for detectors A and B, respectively. These results are consistent with the calculation. Supplementary Fig. 1 shows the energy spectra of the burst with detectors A and B. Enoto et al. 14 simulated the de-excitation emission, considering atmospheric scattering of the gamma rays and moderate energy resolution of BGO crystals. The emission model from 15 N and ambient nuclei, such as Al and Si, well reproduces the results of both detectors A and B. From the spectral and temporal analyses, we confirm that the observed short burst is caused by neutrons produced via atmospheric photonuclear reactions.
Positron production by beta-plus decay. After neutrons are expelled from 14 N and 16 O, unstable nuclei 13 N and 15 O start emitting positrons via β+ decay with half-lives of 10 and 2 min, respectively. Positrons immediately annihilate and emit 511 keV annihilation gamma rays. Supplementary Fig. 2a-d shows count-rate histories in the 0.4-0.65 and 0.65-30.0 MeV bands. Whereas detector A recorded no enhancements after the short burst, detector B recorded an afterglow in the 0.4-0.65 MeV band for the period 0 s < t < 10 s. The count rates decreased with a decay constant of 6.0 ± 2.1 s. The background-subtracted photon count in the 0.4-0.65 MeV band for 1 s < t < 10 s is (2.0 ± 0.4) × 10 2 photons. The backgroundsubtracted energy spectrum is shown in Supplementary Fig. 2e. The centre energy of the line emission is 528 ± 14 keV, which is consistent with 511 keV of the annihilation line within error.
These results lead us to conclude that a positron-emitting region filled with 13 N and 15 O were produced in the atmosphere by the photonuclear reactions, and then passed over detector B flown by the ambient wind flow 14 . Considering that the count-rate history shows a monotonic decrease, the positron source might be generated somewhere above detector B or downwind.
Radio observations. The LF network has five stations ( Supplementary Fig. 3a). Each station has a flat plate antenna sensitive to 0.8-500 kHz. Analogue outputs from the antenna are sampled by a 4 MHz digitiser, whose absolute timing is calibrated with the GPS signals. The LF network can locate radio pulses with the time-of-arrival technique. Supplementary Fig. 3b, c shows the entire LF waveforms of the observed lightning discharge.
The ELF receiver is installed in Kuju (33.059°N, 131.233°E) as a station of the Global ELF Observation Network operated by Hokkaido University. The station has two horizontal search coil magnetometers sensitive to 1-100 Hz magnetic-field perturbations in the east-west and north-south directions. The analogue output is sampled by a 400 Hz digitiser. The direction-of-arrival of the ELF pulses can be confirmed with the magnetic-detection-finder technique. Supplementary Fig. 3d shows the observed waveform in the ELF band.
The JLDN reported two other discharges besides the TGF-associated radiofrequency pulse: an IC of −14 kA at t = 18.7 ms and a CG of −13 kA at t = 228.6 ms. Supplementary Fig. 3b, c shows the corresponding LF pulses. Since these pulses occurred long after the observed TGF, we consider that they were not associated with the high-energy phenomena.
Simulation of gamma-ray glow. We performed Monte Carlo simulations of electron propagation in the atmosphere to reproduce the count-rate histories and energy spectra, using GEANT4 (ref. 47 ). We assume that electron avalanches towards the ground developed in thundercloud, and that the electron spectrum of the RREA at the end of the region has the shape of exp(−ε/7.3 MeV) 49 , where ε is the electron energy. We also assume that the distribution of the electron flux in the avalanche region is circularly symmetric and has no intrinsic time fluctuation. These assumption should be reasonable, given that the count-rate history of detector A is symmetric about the peak, and that the wind velocity was approximately constant (see Wind estimation with X-band radar).
The energy spectra of bremsstrahlung gamma rays from the avalanche electrons approximately follow ε −Γ exp(−ε/7.3 MeV) 54 . The photon index Γ is determined from the source altitude h and offset from the source centre. Count-rate histories depend on the size of the RREA region, wind velocity, and h. The distribution of gamma rays is more diffuse at a higher source altitude due to atmospheric scattering, hence resulting in a longer and fainter gamma-ray glow.
First, we tested a disk-like region with a uniform electron flux, varying h and disk radius in our simulations. Supplementary Fig. 4a shows some examples of the simulation results at various altitudes. Comparing the simulation results with the observation, h = 1500 m is required to reproduce the observed count-rate histories, whereas Γ of the energy spectra indicates h = 900 m. Since any other conditions cannot satisfy both the spectra and count-rate histories, this uniformed-disk model is thus rejected in this analysis.
Then, we considered two disk-like models in which the spatial distribution of the electron flux follows either of the two functions of a distance from the RREA centre l: a Gaussian model, exp(−l 2 /2σ 2 ) and an exponential model, exp(−l/L). The parameters σ and L are free parameters, which denote the spatial extent of the surface brightness of the emission.
We found that both models can reproduce the obtained count-rate histories and spectra; The estimated parameters are h = 600 m and σ = 200 m for the Gaussian model, and h = 400 m and L = 150 m for the exponential model. Comparing these two best models, we found that the exponential model explains the observation better, particularly for the count-rate histories of detector B ( Supplementary  Fig. 4b). Therefore, we employ the exponential model as a working hypothesis to interpret the observation.
Seed electrons. We assume that the seed electrons of the RREA processes are mainly produced by cosmic rays. To calculate the electron fluxes of secondary cosmic rays, we employed Excel-based Program for calculating Atmospheric Cosmic-ray Spectrum (EXPACS) 55,56 , which calculates the flux and spectrum of cosmic-ray particles as a function of an altitude, latitude, longitude, and solar modulation. We extracted electron spectra at an altitude h of 300-2000 m, and then integrated the spectra to obtain the electron fluxes F seed in the energy range of 0.3-50.0 MeV. The electron flux was found to increase exponentially as a positive function of altitude, given by F seed = 2.56 × 10 −3 × exp(h/1890 m) electrons cm −2 s −1 .
Carlson et al. 57 have considered 1 MeV seed electrons produced by cosmic rays. Kelley et al. 22 employed it, and derived the seed flux to be 0.25 cm −2 s −1 at 14.1 km. Our calculation with EXPACS gives the electron flux at 14.1 km of 0.86 cm −2 s −1 . Given that Carlson et al. took a more thorough approach than ours by simulating the effective seeding efficiency for various particles, energies, and geometries, our method might have overestimated the seed electron flux. Regardless of the potential errors in our method, our conclusion that the gammaray glow requires relativistic feedback is unaffected, because overestimation of the seed flux, even if it was the case, would result in an underestimation of the multiplication factor.

Data Availability
The data sets generated and analysed during the current study are available from the corresponding author on request.