Spectral Analysis of Individual Terrestrial Gamma‐ray Flashes Detected by ASIM

at source. In the thermal runaway mechanism sufficiently many runaway electrons are produced in the high Abstract The Atmosphere-Space Interactions Monitor (ASIM) is the first instrument in space specifically designed to observe terrestrial gamma-ray flashes (TGFs). TGFs are high energy photons associated with lightning flashes and we perform the spectral analysis of 17 TGFs detected by ASIM. The TGF sample is carefully selected by rigorous selection criteria to keep a clean sample suitable for spectral analysis, that is, suitable count statistics, low instrumental effects, and reliable source location. Monte Carlo modeling of individual TGFs has been compared to the observed energy spectra to study the possible source altitudes and beaming geometries. A careful model of the instrumental effects has been developed and validated. Several combinations of source altitudes and beaming geometries are accepted by the statistical tests for all the TGFs in the sample resulting in a large uncertainty in the estimate of the intrinsic source luminosity. The analyzed TGFs show significant variations in observed fluence independent of the distance between source and ASIM. A lower limit on the maximum photon energy produced by TGFs is estimated to be 24 MeV for the analyzed TGFs. The intrinsic limitations of TGF spectral analysis from space are also investigated and it is found that the ability to constrain the source altitude and beaming geometries of TGFs strongly depends on the distance between source and satellite. With the current generation of instruments with effective areas in the range of few hundreds cm 2 , it is very difficult to constrain reliably the source properties without the help of simultaneous measurements in the radio band.


Introduction
Terrestrial gamma-ray flashes (TGFs) originate from Earth and are associated with lightning flashes. They are observed by satellites as sub-millisecond bursts of energetic photons and were first reported by Fishman et al. (1994) using data from the BATSE experiment onboard the Compton Gamma-Ray observatory. TGFs have also been detected by RHESSI , AGILE (Marisaldi et al., 2010), the Fermi Gamma-ray Burst Monitor (GBM) , and BeppoSAX (Ursi et al., 2017). The first space mission specifically designed to detect TGFs is the Atmosphere-Space Interactions Monitor (ASIM) operational since June 2018 .
The TGF source has mostly been characterized by three observables: fluence, source altitude, and beam opening angle (sometimes referred as angular distribution). As most satellites detecting TGFs were not designed to handle the very high count rate of TGFs, a large fraction of the photons are not detected due to dead-time effects. To increase count statistics Dwyer and Smith (2005) combined data from RHESSI TGFs to build a cumulative spectrum that they compared to Monte Carlo models. They estimated the production source altitude of TGFs to be 15-21 km. Hazelton et al. (2009) used RHESSI TGFs and took into account the radial distance using lightning sferics detected by the World Wide Lightning Location Network (WWLLN) to estimate the geolocation of the TGFs. The radial distance is the distance following the surface of the Earth between the TGF foot-point and the sub-satellite point. For the cumulative spectrum of TGFs within 300 km from the subsatellite point, their most likely model had a wide-beam geometry and a source altitude of 15 km. For TGFs produced with a radial distance larger than 300 km they showed that a narrow-beam model at ≥ 20 km altitude was unlikely to reproduce the observed data. Gjesteland et al. (2011), also using RHESSI data, concluded that an isotropic emission cone with a half angle between 30° and 40° was the best fit to the data.
Spectral analysis of individual TGFs is complicated due to the photon transport through the atmosphere and high count rates. Østgaard et al. (2008) estimated a production altitude up to 40 km using BATSE data and Monte Carlo modeling, however Grefenstette et al. (2009) realized that the measured brightness of TGFs, observed by RHESSI and BATSE, were underestimated due to dead-time effects. When Gjesteland et al. (2010) revised the analysis on the BATSE data, taking dead-time into account, the new production altitude estimates were below 26 km. Mailyan et al. (2016) performed spectral analysis on 46 individual TGFs, detected by Fermi-GBM, finding spectral diversity. This implies that cumulative TGF spectra miss important information as spectral diversity gets smeared out. The observed data were compared to Monte Carlo simulations of a large scale RREA model, including the propagation through the atmosphere. They tested production altitudes of 11.6, 13.4, 16, and 20.2 km, with narrow-and wide-beam models. The narrow beam model is modeling the intrinsic angular distribution of bremsstrahlung photons in the RREA region. In the wide-beam model the photons are emitted isotropically in a cone with half angle 45°. For most of the TGFs they could not distinguish between narrow-and wide-beam models. Models produced at 11.6 km gave frequently a better fit than higher altitude models, however 13.6 km could not be rejected in those cases. There are also events where 20.2 km models gave the best fit, but lower altitude models could not be rejected. Of the 46 TGFs 6 had poor fits, but it is speculated that this is due to instrumental effects. Their TGF sample had between 21 and 53 counts per TGF in one BGO detector. A follow-up study by Mailyan et al. (2019) compared also lightning leader models to TGFs detected by Fermi-GBM. They found that for most of the TGFs one mechanism was not favored over the other.
In addition to spectral analysis from space, low frequency (LF) radio measurements can provide an independent estimate of the production altitude. These independent LF measurements estimate TGF production altitudes between 10 and 15 km (Cummer et al., 2014;Pu et al., 2019). Very low frequency radio measurements of sferics produced by lightning can be used to geolocate TGFs. TGFs and sferics are simultaneous within few hundred microseconds (Connaughton et al., , 2013Lindanger et al., 2020).
In this work we present the individual spectral analysis of 17 TGFs detected by ASIM with a number of counts between 82 and 254 in the High Energy Detector (HED). 11 of the TGFs also have data in the Low Energy Detector (LED) down to 60 keV. We will present the ASIM instrument and discuss instrumental effects and its impact on the energy measurements during TGF detection. The forward modeling from the source to the detection of the photons in ASIM is explained in detail, and we will show normalized energy spectra of the different models observed at different radial distances. The results of the spectral analysis are presented in detail for one TGF, while the results of the other TGFs are presented in table format with further documentation in Supporting Information S1. We will then discuss the source properties of TGFs and we will also address the intrinsic limitations of the spectral analysis approach.  ;Østgaard, Balling, et al. (2019) give a thorough description of the ASIM mission and the instruments onboard.

The MXGS and MMIA Instruments Onboard ASIM
The scientific payload onboard ASIM consists of the Modular X-and Gamma-ray Sensor (MXGS) and the Modular Multispectral Imaging Assembly (MMIA). MMIA includes two cameras and three high-speed photometers used for lightning and TLE detection. The main instrument for this analysis is the MXGS, which consists of two detectors: HED and LED. HED consists of 12 Bismuth-Germanium-Oxide (BGO) scintillators, each connected to a photomultiplier tube (PMT). Each BGO is sensitive to photons between 300 keV and >30 MeV. The time resolution of HED is 27.8 ns and each BGO-PMT detector has a dead-time of ∼550 ns.
LED consists of pixelated Cadmium-Zink-Telluride (CZT) detector crystals that are sensitive to photons with energies from 20 to 400 keV. The time resolution of LED is ∼1 µs. A total of 16,384 pixels together with a coded mask in front of the detector makes LED capable of imaging the incident TGFs by using a deconvolution technique. TGFs detected in LED are routinely analyzed in search for source location using the imaging techniques. The pixels are divided into 16 independent readout chains, and each chain is arranged into 8 Application Specific Integrated Circuits (ASIC) with a dead-time of 1.4 µs. If two hits occur in different ASICs in the same chain within 1.4 µs the two hits will be distinguishable, however their energies will be added together. These multi-hit events are flagged. The LED energy calibration was performed for each individual pixel, as described in Østgaard, Balling, et al. (2019), by the pre-launch calibration and the onboard energy calibration system, and has been found to be stable. The energy resolution of LED is <10% at 60 keV. LED and MMIA are only operative during nighttime.
The absolute timing accuracy of MXGS is between ∼0 and ∼30 ms due to a non-optimal timing interface between the ASIM payload and the ISS. The relative timing uncertainty between MMIA and MXGS was ±80 µs prior to April 2019. A software update in April 2019 reduced the relative timing uncertainty to ±5 µs . When MMIA is active and there is lightning activity in its field of view, it is possible to align MMIA photometer data with several ground detected sferics to correct the absolute timing of MMIA down to less than 2 ms. We can then correct the absolute timing of MXGS accordingly.

Measurements and Instrumental Effects in HED
The energy resolution of HED is <20% at 511 keV (Østgaard, Balling, et al., 2019). The energy channels are converted to keV by a quadratic fit to the 511 keV, 1,275 keV, and 31 MeV peaks in the background energy spectrum. 511 and 1,275 keV are produced by the onboard 22 Na calibration source, and 31 MeV is the most probable value of the energy deposition produced by cosmic protons. The energy deposit of the cosmic protons was estimated by the Geant4 ASIM mass model (Agostinelli et al., 2003;Østgaard, Balling, et al., 2019). The dead-time of each individual HED detector is 550 ns and if there is more than one hit within this dead-time, only one count will be recorded and the energy measurement will result from a combination of the two pulses. This is called pulse-pileup. If we have a new count in HED after the 550 ns dead-time, the recorded energy will be the sum of the energy of the new count and the tail of the previous count. This is called tail-pile-up and is well illustrated in Figure 6c in Østgaard, Balling, et al. (2019). The HED data acquisition system is keeping track of tail-pile-up. The counts are labeled as either a normal-, valley-or fast events. The fast events are tail-pile-up events. Normal events are events that are not sitting on the tail of a previous pulse. The valley event is the local minimum between the fast event pulse and the previous pulse and is only associated with fast events. Valley events are not physical counts, but can be used to reconstruct the signal amplitude of fast events. For fast events, a lower energy threshold is implemented. This threshold has been adjusted a few times during the mission.
In the HED PMTs we have a voltage drop over the PMT dynode chain after high energy counts are detected. The time before the high voltage restores itself depends on the energy of the previous count and can be up to ∼30 μs for counts above 30 MeV, and ∼0 μs for counts less than 400 keV when there is no voltage drop. As the voltage is not constant over the dynode chain, the gain varies. Therefore, the channel to keV energy calibrations is not valid during the voltage drop as the pulse-heights of the following counts are measured incorrectly.
LINDANGER ET AL. We account for this by implementing an energy dependent "safety time", which is the time the high voltage in the dynode chains needs for recovering to a level where the energy is within 20% accuracy. Counts that do not meet this criterion are not used for spectral analysis, but are still valid counts for light curves and fluence. The safety time does not bias the energy spectrum if we assume a constant source spectrum during the production of the TGF. The safety time and the voltage drop are further discussed in Appendix A. As HED consists of 12 independent BGO-PMT detectors, a voltage drop introduced by a high energy count in one of the detectors does not affect the other 11 detectors. Dead-time, energy resolution, fast events, pulse-pile-up and safety time are taken into account in the instrument model discussed in Section 3.2.

TGF Sample Selection
Three criteria have to be fulfilled in spectral analysis of individual TGFs: (a) a bright TGF with good count statistics, (b) the geolocation of the TGF is known, (c) minimal instrumental effects. The first criterion is due to that we need good count statistics of the TGF so that the statistical testing is not limited by low count statistics in the detectors. We set this criterion to a minimum of 80 counts in HED.
The second criterion is necessary as we need to know the TGF production location in order to have the correct propagation path from source to MXGS. This is done by correlating lightning data, detected by ground based lightning location networks, to the time of the TGF. Lightning data are obtained from WWLLN (Rodger et al., 2009), and GLD360 provided by Vaisala Inc (Said & Murphy, 2016). Both networks provide time and geolocation by detecting sferics produced by lightning flashes. For each TGF that satisfies the criterion of minimum count statistics we inspect the ground-based lightning data. We inspect the full map of lightning activity in a time interval around the TGF, to see whether there are multiple active lightning cells at different locations. If there is one sferic correlated with the TGF within the absolute timing uncertainty, we assume it is the sferic associated with the TGF. If there are sferics from different locations within the uncertainty, we cannot perform spectral analysis on the TGF as we do not know its geolocation. If an independent source location is provided by LED imaging, this will be considered as an independent check of the lightning location.
The third and last criterion is to avoid pile-up and saturation for very high flux levels, which is a common challenge for all gamma-ray instruments. For HED such pile-up and saturation will manifest itself as a large number of fast events, which will be removed according to the safety time criteria. For LED they will be labeled as multi-hits, which cannot be used for spectral analysis. To assure that the accepted counts from HED and LED used for spectral analysis are representative for the event, we claim that a fraction of fast to normal events should not exceed 20% for HED, and the multi-hits in LED should not be more than 30% of the total. This criterion will bias the TGF sample toward low intensity and longer duration TGFs.

Forward Modeling
To be able to compare modeled TGFs with ASIM measurements, we need to model the TGF source, the photon propagation through the atmosphere, the deposition of counts in LED and HED, and instrumental effects in the detectors. We model TGFs at different altitudes and different source beam opening angles to check which models fit the measured energy spectrum.
The forward modeling is computed in several steps. For a particular TGF produced at an altitude with a certain beaming the following steps are performed: 1. Source assumptions. 2. Photon propagation through the atmosphere from source to satellite altitude. 3. The time-energy distribution from step 2 is used as an input distribution to the ASIM mass model. 4. The output time-energy distribution from step 3 is used as an input distribution to simulate duration and instrumental effects. 5. The energy distribution from step 4 is compared to the measured data.
In step 1 the TGF source is modeled assuming a RREA photon spectrum proportional to the analytical function 1/E ⋅ exp(−E/(7.3 MeV)), with a maximum energy of 40 MeV. This is the only TGF source model spectrum we 5 of 18 assume and we will explore beam width, production altitude, and fluence. In the following, we will refer to the term model as any tested combination of beam width and altitude. We do not consider tilting of the TGF beam to prevent overfitting. After generating TGF photons with a RREA energy spectrum, we perform step 2 where the photons are beamed upward with a Gaussian profile with a σ θ that defines the beaming profile. Figure 1 illustrates a Gaussian beaming profile and shows a simulated TGF in the satellite altitude plane at 408 km with σ θ = 30°. At source this is a two dimensional Gaussian, meaning that 39% of the photons are within one σ θ before they are further scattered in the atmosphere. We elaborate more on how to compare isotropic and Gaussian beaming models in Appendix B. The time profile at source is instantaneous, and the time profile at satellite altitude is due to the transport through the atmosphere in addition to the source duration that will be implemented in step 4.
In step 2, we use the Geant4 Monte Carlo code presented in Sarria et al. (2019) to model the photon propagation through the atmosphere. The atmospheric composition is obtained using the NRLMSIS-00 model (Picone et al., 2002). Geant4 takes into account all the relevant physical processes of photon, electron and positron transportation (Compton scattering, Rayleigh scattering, photo-electric absorption, electrons inelastic and elastic scattering, bremsstrahlung, and positron annihilation). The magnetic gyration of electrons and positrons produced by the gammas are not taken into account as this enhance computation time and are only relevant for Terrestrial Electron Beams (Ebert et al., 2010;Köhn & Ebert, 2015;Sarria et al., 2019Sarria et al., , 2021. The energy, position, momentum and timing information of photons reaching satellite altitude are saved. Note that the energy distribution of the photons from an instantaneous source is dependent on time after propagation through the atmosphere as shown in Figure 6 in Marisaldi et al. (2019).
Step 2 is simulated for TGFs produced at 9, 11, 13, 15, 17, and 19 km altitude, for σ θ = 5°, 10°, 15°, 22° and 30°. In total, this gives 30 time-energy spectra matrices dependent on production altitude and σ θ . For the 30 time-energy spectra (from an instantaneous source) we calculate different distributions depending on the radial distance. Figure 2 shows a TGF produced at an altitude of 13 km and σ θ = 30° observed at different radial distances. We see that the energy spectra get softer at larger radial distances due to a larger fraction of Compton scattering compared to the energy spectra at smaller radial distances. Figures 3a-3c compares two models with σ θ = 15° produced at an altitude of 9 and 19 km, Figure 1. Photons produced by a modeled terrestrial gamma-ray flash (TGF) at latitude = 0° and longitude = 0°, reaching the satellite plane at an altitude of 408 km. The TGF is beamed upward with a Gaussian profile with σ θ = 30° from a source altitude of 15 km. The black ellipse marks the 30° off-nadir angle from the source. 10.1029/2021JD035347 6 of 18 observed at radial distance 0, 150, and 420 km. Figures 3d-3f compares two models produced at an altitude of 13 km with σ θ = 10° and 30°. These plots are further discussed in Section 5.1.
In step 3 the ASIM mass model uses the time-energy spectra from step 2 as input to model the detection of counts. The radial distance is fixed based on the TGF source position obtained by lightning detection data. This produces 60 (HED and LED) new time-energy spectra for the given radial distance. The ASIM mass model is a full scale Geant4 model of the ASIM instrument mounted on the Columbus module on the ISS (Østgaard, Balling, et al., 2019). Note that we removed from the Monte Carlo code the ACES instrument from Figure 11 in Østgaard, Balling, et al. (2019) as it is not present on ISS together with ASIM.
In step 4 the output spectra from step 3 are used to simulate TGF duration and instrumental effects in HED. We generate a very large photon list sampled from the time-energy output spectra from step 3. In order to account for the physical duration at source we add a random number, sampled from a Gaussian time distribution, to each 10.1029/2021JD035347 7 of 18 simulated time tag. The standard deviation of the Gaussian time distribution is fitted such that the modeled time profile matches the duration of the measured TGF. This duration and the time profile after propagation through the atmosphere does not bias the analysis as the safety time criteria are implemented both in the forward modeling and on the observed data. The energy spectrum after step 4 consists of a large number of separate TGFs simulated under the same conditions, including instrumental effects, cumulated to obtain a sample spectrum with large count statistics. We take into account the individual energy range of each of the 12 HED detectors, energy resolution, pulse-pileup, lower energy threshold for fast events and the previously discussed safety criteria in Section 2.2. This gives a modeled spectra that we can compare to the measured spectra. For LED we only apply the energy resolution to the output from the mass model.

Statistical Tests
The output from step 4, which is given for a range of altitudes and beaming angles can be compared with the measured spectra when we know the geolocation of the TGF and therefore the incident angle of photons. To find the most likely altitude and beaming angle we use maximum likelihood estimation (MLE) and χ 2 statistics to compare and quantify the goodness of fit. MLE based on Poisson statistic is the best approach to estimate parameters in counting experiments, especially if we have small total numbers of counts (Hauschild & Jentschel, 2001). We estimate the best fit model using MLE by minimizing − 2log  Mailyan et al., 2016) with where O i and M i are the observed and modeled number of counts in each energy bin i. We will refer to − 2log  as the negative log likelihood (NLL). The modeled energy spectra are fitted to the measured energy spectrum by minimizing NLL. When the difference in NLL between two models is greater than five, the best fit model is preferred above the other model with greater than 99% confidence level (Mailyan et al., 2016;Sarria et al., 2021). This method does not state if the best model is any good in an absolute term, but only states the relative goodness. Therefore we also use the Pearson 2 statistic to quantify the goodness of fit.
The reduced 2 is given by where ν is the number of degrees of freedom and is given by ν = "number of bins" − "number of free parameters" − 1, where the free parameters are production altitude, beaming angle, intensity of HED, and intensity of LED. That gives three free parameters using only HED data and four free parameters if LED data are available.
If the 2 statistic is close to one it is a good fit. If it is below the critical value given by the χ 2 -distribution, the model is considered compatible with the measurement with a 95% significance level. If it is above, the model is rejected. The energy bins in both MLE and χ 2 statistics are chosen by the following criteria: (a) At least five measured counts in each bin, (b) every bin has the same number of counts according to the model, (c) as many bins as possible (James, 2008).
A third statistical test is also implemented to assess the absolute goodness of fit using MLE. This is done by generating a NLL distribution by calculating the NLL for each individual simulated TGF in step 4 to the full energy spectrum of all the simulated TGFs after step 4. The NLL for the observed spectrum compared to the full simulated energy spectrum is compared to the NLL distribution (Lyons, 1986). We chose a 90% confidence level to reject the model, that is, the model is rejected if 90% or more of the NLL distribution have a lower NLL than the observed data. There is a generally good compatibility between the three statistical tests. Therefore, this third method is not discussed further in this study but a comparison between the three statistical tests is shown in Figure S1 in Supporting Information S1.
When we have data from both HED and LED, multi-hits in LED are removed and we scale both simulations to the HED and LED data independently. We do not model the multi-hit counts in LED because they do not take part LINDANGER ET AL.
10.1029/2021JD035347 8 of 18 in the spectral analysis. This is the reason why we keep the normalization between HED and LED independent. The energy range considered for HED is 400 keV to 40 MeV, and the energy range considered for LED is 60-350 keV.

Results
The total number of TGFs detected by ASIM is about 900 TGFs between June 2018 and December 2020. Of these, 17 TGFs meet the TGF sample selection where we have good count statistics, known production geolocation, and minimal instrumental effects. The TGFs are shown in Figure 4 where we plot the number of counts per TGF versus duration. T 90core is defined as the minimum duration where 90% of the counts in the TGF are included. Table 1 shows the observed properties of the TGFs included in the analysis. The TGFs are ordered in increasing radial distance from the ISS footpoint. The datetimes marked with an asterisk are time corrected down to less than 2 ms by correlating MMIA photometer data with ground based lightning detections. The total number of counts per TGF in HED and LED are shown, and we also show the number of counts used in the spectral analysis. The counts used for spectral analysis are counts inside the suitable energy range and not removed by safety time criteria (HED) or multi-hits (LED). If there are no LED data for the event the number of counts is marked with a dash. The radial distance between the subsatellite point and the source geolocation of the TGF foot-point are given together with the latitude and longitude of the source geolocation. The TGF geolocation is estimated from lightning detections associated with the TGF and independently checked with imaging results from LED when available. The duration of each TGF in HED is presented as both T 50core and T 90core . These are the minimum durations where 50% and 90% of the counts in the TGF are included. The standard T 50 and T 90 are calculated by taking the cumulative count distribution of  Note. The durations T 50core and T 90core are calculated for HED data. The latitude and longitude are the geolocation of the sferic associated with the TGF, and the radial distance is the distance between the subsatellite point and the geolocation of the associated sferic. The datetimes marked with an asterisk are time corrected down to less than 2 ms absolute timing accuracy. 10.1029/2021JD035347 9 of 18 the signal and calculate the duration between 25% and 75% (5% and 95%) of the count distribution. However, this method is more sensitive than the core method to where the start and the end of the signal is defined. As it is not always straight forward to define the exact starting-and ending point of a TGF, the core method is preferred.
We present the results of the spectral analysis in Table 2. For each individual TGF the accepted altitude range according to MLE-and χ 2 statistics is presented for different beaming angles σ θ . We also have a column for the best fit model. If LED data are available they are included in the spectral analysis, however the results using only HED are not very different from the results using both HED and LED. Note that TGF 9 and TGF 15 are originating from the same thunderstorm with ∼1 min separation.
For clarity we will present the analysis results of TGF 2 in detail. Figure 5a) shows WWLLN and GLD360 detections within ±5 min of the TGF, and within ±30 s of the TGF. We also show the lightning detection that is closest in time to the TGF within 0-30 ms and within a radial distance of 1,000 km as a triangle. This lightning detection is assumed to be the lightning associated to the TGF. In the center of the map we plot the subsatellite point surrounded by a circle with a radius of 500 km. Figure 5b shows a time-energy scatter plot of the TGF detected in HED and LED. Only the counts with correctly measured energies are shown. Figure 5c shows the light curve of all the counts detected in HED and LED, and the light curve using only counts used in the spectral analysis. Figure 5d shows the results of the MLE statistical test. As only the relative difference in NLL is used to decide which models are a better fit we plot NLL − min(NLL), where min(NLL) is the minimum NLL value, that is, the best fit model. If the model has a value below the critical value of 5, the model is considered accepted by the MLE, for example, for θ σ = 22°, the accepted altitudes are 13-19 km. Figure 5e shows the results of the reduced χ 2 test. If a model has a value below the critical value shown with the dashed line, the model is considered accepted. In general there is a good agreement between the MLE-and χ 2 -tests except for models that are close to the critical values of both tests. Figures 5f and 5g shows the measured energy spectrum together with the best fit model. Figure 5f shows the energy spectrum with counts on the y-axis. As discussed in Section 3.3 the energy bins are selected to have a flat model spectrum. The uncertainty of the data points is ±1 standard deviation assuming Poisson statistics. Figure 5g shows the differential energy spectrum. Figures with the same format as Figure 5 are found in Supporting Information S1 for each TGF analyzed in this study. For four of the TGFs the Note. The table shows the accepted production altitudes in km for different σ θ .
where f o is the observed fluence and f m is the modeled fluence for the given production altitude, beaming, and radial distance, assuming 10 17 photons above 1 MeV at source.

Discussion
As this is the first spectral analysis of TGFs detected by ASIM, we keep the TGF sample as clean as possible selecting only 17 of ∼900 TGFs. The number of suitable TGFs is low because of strict requirements on the reliability of the energy estimate, good count statistics, and on the availability of a reliable geolocation. The sample size would increase if we relaxed the minimum number of counts or the maximum allowed fraction of counts affected by instrumental effects. However, it would not result in a better quality of the scientific interpretation as the uncertainty in the results would increase. In this study, we cannot claim any generalization of the properties of TGFs as we analyze a small sub-sample that is biased by the TGF sample selection discussed in Section 3.1.
All of the 17 TGFs analyzed in this work have a GLD360 lightning match. The two briefest duration TGFs, TGF 4 and 13, also have a WWLLN match at a position compatible with GLD360 within location uncertainties. This is in agreement with Connaughton et al. (2013); Lindanger et al. (2020) that show that brief duration TGFs are more likely to have a WWLLN match.

Source Properties of TGFs
ASIM TGFs are compared to models with source altitude of 9-19 km and Gaussian beaming angle σ θ from 5° to 30°. A larger beaming angle, σ θ = 40°, was simulated but not included in the study as it gave very similar results to σ θ = 30°. There were no cases where all altitudes for σ θ = 30° were accepted while all altitudes for σ θ = 40° were rejected, or opposite. As no TGFs have an accepted model by MLE with σ θ = 5° we may reject σ θ = 5° as a plausible TGF beaming based on the 17 TGFs in this sample. This is consistent with the minimum beam width expected from bremsstrahlung produced by RREA in a perfectly uniform electric field . Figure 2a in Hazelton et al. (2009) shows approximately a Gaussian distribution with FWHM = 18°, which LINDANGER ET AL.

10.1029/2021JD035347
12 of 18 corresponds to σ θ = 7.6°. Therefore, the lack of good fits with σ θ = 5° is a consistency check, because the model, if accepted, would be incompatible with the basic physics of photon scattering in air. Cummer et al. (2014); Pu et al. (2019) estimate TGF production altitudes between 10 and 15 km using LF radio measurements. Heumesser et al. (2021) get similar results by modeling light propagation through clouds assuming typical values of size and density of cloud particles. Mailyan et al. (2016) analyzed Fermi-GBM TGFs and found that 11.6 km models gave frequently a better fit than higher altitude models, but 13.6 km could not be rejected for those cases. Sometimes also TGFs had a best fit at 20.2 km, but lower altitude models could not be rejected. Table 1 in Mailyan et al. (2016) shows which models are accepted and even though the best fit model is highlighted, most of the other models are also accepted for most of the TGFs. This is in agreement with our results, see Table 2, where only a few models are rejected per TGF. The MLE and χ 2 tests state that if a model is below or above a critical value it is accepted or rejected. If a model is accepted it could explain the observation even though there exists a better fit. Therefore the "best fit model" should be handled with caution and not used to draw general conclusions on the source spectrum. Figures 3a-3c show that it is increasingly harder to distinguish 9 and 19 km production altitude with increasing radial distance. Remember that we only have ∼100 counts per TGF and the number of counts at high energies are few. Figures 3d-3f shows that it is very hard to distinguish beaming angles σ θ at a radial distance of 150 km. For a radial distance of 420 km it is easier to distinguish the beaming angle as the observation point is then outside the direct photon beam for small σ θ relative to large σ θ , seeing a larger fraction of Compton scattered photons softening the energy spectrum for small σ θ . These properties of TGF modeling are reflected in the results in Table 2. It is easier to distinguish altitudes for TGFs observed at smaller radial distances, and it is easier to distinguish σ θ at larger radial distances.
The TGF id's are sorted in increasing radial distance. At distances larger than TGF 9 (193 km) only one TGF has an accepted model with σ θ = 10°. As it is easier to distinguish σ θ at larger radial distances, this may indicate that in general the beaming of TGFs has σ θ ≥ 15°. As discussed in the beginning of this subsection, the intrinsic beam width of bremsstrahlung photons from RREA in a perfectly uniform field correspond to σ θ ≈ 7.6°. It is not unreasonable to expect further widening due to non-perfectly uniform electric field in a thundercloud that leads to σ θ > 10°. However, based on the 17 TGFs analyzed we cannot state, based on observations, that in general σ θ ≥ 15°. Also to distinguish σ θ we need to observe TGFs at large radial distances. If we assume that TGFs with smaller σ θ exist, they would not have been included in this analysis due to few observed counts at large radial distances, that is, our TGF sample is biased according to our selection criteria. Figure 6a shows the observed fluence of the detected TGFs for all accepted models according to the MLE test. If all TGFs had the same brightness at source the observed fluence would decrease with increasing radial distance. This is not evident in Figure 6a as TGF 10 and 17 are observed at large radial distances with high observed fluence. This supports the idea that TGFs have a wide range of brightness at source or significant tilting with respect to the vertical axis. Figure 6b shows the estimated number of photons at source with energies above 1 MeV. All the accepted models are shown with different colors representing altitudes, and markers representing σ θ . It is clear from the plot that the uncertainty of the number of photons at source is three orders of magnitude. This is due to high photon absorption in the atmosphere from 9 km to space, compared to 19 km to space. A wide range of altitudes is accepted for each TGF. The best fit models are indicated with a black circle and range from 1.9 ⋅ 10 16 to 1.5 ⋅ 10 20 photons at source as both 9 and 19 km altitude are best fit models dependent on the TGF. This is roughly in agreement with Figure 8, top panel, in Mailyan et al. (2016) that shows the number of relativistic electrons above 1 MeV at source. The number of electrons in Mailyan et al. (2016) varies from 4 ⋅ 10 16 to 3 ⋅ 10 19 . Conversion from relativistic electrons to bremsstrahlung photons is done using Equations 5 and 9 in Dwyer et al. (2017) resulting in a bremsstrahlung photons to relativistic electrons ratio of 0.33 for an electric field strength at sea-level of 400 kV/m. This is the electric field assumed in Mailyan et al. (2016). In order to compare the results in Mailyan et al. (2016) with ours, we convert from electrons to photons and get a range of 1.3 ⋅ 10 16 to 10 19 photons. The one order of magnitude difference for the maximum number of photons at source can well be explained by that the lowest production altitude modeled in this study is 9 km, while in Mailyan et al. (2016) the lowest altitude modeled was 11.6 km. The uncertainties in Figure 6b would be reduced by several orders of magnitude if the source altitude of TGFs were estimated by for example LF radio measurement (Cummer et al., 2014;Pu et al., 2019), lightning mapping arrays (LMA) (Lu et al., 2010), radar measurements (Mailyan et al., 2018), or the electric field antenna part of the TARANIS mission (Lefeuvre et al., 2008). Given the accepted idea that a TGF is produced inside a cloud it is not likely that production altitudes above the tropopause are physical. Even LINDANGER ET AL.

10.1029/2021JD035347
13 of 18 if production altitudes of 17 and 19 km are accepted solutions of the spectral analysis, they may not be realistic as the tropopause height is on average ∼ 16.5 km in the tropical latitude band of the TGFs analyzed (Seidel et al., 2001). One would need a significant overshooting top and TGF production close to the top for a TGF source to be located at 17-19 km altitude. Production altitudes above the tropopause are also not in agreement with typical values from Cummer et al. (2014); Pu et al. (2019), andHeumesser et al. (2021).
The maximum energy count per TGF of the 17 individual TGFs analyzed in this study is in the range 16-28 MeV. The median of the maximum energies is 21 MeV. There is a total of 8 counts with energies above 24 MeV during the 17 TGFs. Note that this is not the photon energy, but the energy of the count in the detector. A larger photon energy is expected, as a partial energy deposit in the detector is likely at these energies. To assess if these 8 counts above 24 MeV originates from the TGFs and not background radiation we use data between −900 ms and −100 ms prior to the TGF detection as background. The average background rate was found to be 482 counts above 24 MeV per second. Given this background rate, the Poisson probability of having 8 counts or more above 24 MeV during the 17 TGFs (4060 µs) is 0.001. Therefore we can conclude that the maximum energy produced by the TGFs is larger than 24 MeV. For comparison RHESSI , AGILE (Marisaldi et al., 2010), and Fermi  teams reported TGF single photon energies of >= 20 , 43, and 38 MeV, respectively. Tavani et al. (2011) reported higher energies detected by AGILE, however, Marisaldi et al. (2019) questioned this after an improved understanding of instrumental effects under high-flux conditions.

Limitations on Spectral Analysis From Space
This study emphasize that the spectral analysis does not constrain well the parameter space of the source models. Therefore we want to investigate what the intrinsic limitations of this method are. To assess the limitations on the spectral analysis of TGFs detected from space, we performed a pseudo spectral analysis by randomly sampling 100 and 1,000 photons from the energy distribution of photons after atmospheric propagation from a production altitude of 11 km and σ θ = 22°. We performed the spectral analysis on these pseudo observations assuming a perfect instrument without taking into account a mass model or instrumental effects. The radial distances considered are 0, 150, and 420 km. The results are shown in Figure 7 where the pseudo observation consist of 100 photon energies on the left side, and 1,000 photon energies on the right side. Note that the results can be slightly different selecting a different random seed for the random sampling. The correct model corresponding to the pseudo observation is always accepted, but it is not always the best fit. Several models can be accepted and in agreement with Hauschild and Jentschel (2001), MLE performs better than χ 2 estimating the correct model. The χ 2 test is unreliable for low count statistics and it should be avoided in such cases, however it can be used as an additional test to give a measure of absolute goodness of fit, not provided by MLE method, if count statistics is large enough. For a radial distance of 420 km, the altitude is nearly impossible to constrain as MLE accepts altitudes between 9 and 17 km even for a pseudo observation of 1,000 photons. It is clear from Figure 7 that for 100 photons the correct model is accompanied by other accepted models. Therefore all accepted models should be taken into account when trying to confine the source properties of TGFs. If independent measurements of the production altitude were provided, the beaming properties may be confined by the spectral analysis, depending on the radial distance between the source and the satellite. This conclusion is only valid assuming no or negligible tilting of the TGF beam. Tilting of the TGF beam is expected both in leader models and models dominated by large-scale electric fields. However, no observations have shown evidence of TGF tilting yet, as several spacecraft need to be observing the same TGF from different locations in order to resolve the degeneracy between tilt and beam opening angle. The results of these simulations provide valuable recommendations for the design phase of future TGF detecting missions: a large effective area, that is, large count statistics, is not enough for the purpose of spectral analysis of individual TGFs, and complementary observations to constrain the production altitude should be planned.

Summary
This study provides the first spectral analysis of TGFs detected by ASIM, comparing individual TGFs to modeling. A sample of 17 TGFs is selected by rigorous selection criteria to keep a clean sample. In agreement with Mailyan et al. (2016) the observed energy spectra are diverse meaning the cumulative spectral analysis of TGFs should be avoided. Monte Carlo modeling of individual TGFs have been compared to observations allowing us to study the possible source altitudes and beaming geometries of TGFs. A careful statistical method using both χ 2 and maximum likelihood estimation is implemented to assess which models fit the observations. A large effort has also been made to properly account for instrumental effects in the high energy detector. For all the TGFs in the sample, several combinations of source altitudes and beaming geometries are accepted by the statistical tests. Tilting of the TGF beam has not been considered in this analysis as adding another free parameter to the statistical tests would lead to overfitting. TGFs may not be centered perfectly in the upward direction, however, only simultaneous TGF observation from a constellation of satellites could solve the degeneracy between beaming angle and tilt.
This work also highlights the limitations of spectral analysis of TGFs from space if no additional measurement can be used to narrow the parameter space, for example by setting a smaller altitude range by radio measurements. The ability to confine source altitude and beaming angle depends on the radial distance in addition to count statistics. All accepted models, according to maximum likelihood estimation, should be considered when trying to confine the source properties of TGFs, not only the best fit. The χ 2 test is unreliable for low count statistics and maximum likelihood estimation should be used instead (Hauschild & Jentschel, 2001).
The analyzed TGFs show diverse observed fluence independent of the distance between the source and ASIM supporting the idea that TGFs have a wide range of brightness at source. The number of photons at source with energies larger than 1 MeV ranges from 10 16 to 10 20 and an independent measure of the altitude, for example by LF-radio (Cummer et al., 2014;Pu et al., 2019) or LMA (Lu et al., 2010;Mailyan et al., 2018), is needed to further constrain the number of photons at source. Based on the 17 analyzed TGFs a lower threshold of the maximum photon energy produced by TGFs is estimated to be 24 MeV.

Appendix A: The HED Safety Time Criteria
For high fluxes and high energy photons the HED instrument suffers from high voltage drops in the PMTs which have to be handled carefully. This high voltage drop is shown in Figure S19 in Supporting Information S1 where flight data from one of the 12 detectors in HED are shown. Figure S19a in Supporting Information S1 shows the first count following a count with energy E 0 = 500 keV ±10%. The fast events are shown in orange color and normal events are shown in blue color. Figure S19b in Supporting Information S1 shows the same data displayed as a 2D histogram. Figure S19c in Supporting Information S1 shows the first count following an energetic count with energy E 0 = 25,000 keV ±10%. Figure S19d in Supporting Information S1 shows the same data displayed as a 2D histogram. The proton peak, which should be at ∼31 MeV, is clearly decreasing toward lower energy channels for smaller dt after an energetic count. This effects is less substantial when the energy of the first count, E 0 , is lower. The white gap for fast events when dt < 10 µs in Figure S19c and S19d in Supporting Information S1 is due to a higher threshold for fast events to be recorded in the data, than for normal events. Note that HED consists of 12 independent BGO-PMT detectors and a voltage drop introduced by an energetic pulse in one of the detectors does not affect the other 11 detectors.
The voltage drop effect was investigated further at ground using a spare BGO-PMT detector used in HED. In the laboratory, three natural radioactive background energy peaks were used together with a weak 60 Co source. The 60 Co source emits gamma-rays with energies 1.17 MeV and 1.33 MeV. From the radioactive background in the laboratory we had 1.46 MeV from 40 K , 2.6 MeV from 208 Tl , and ∼32 MeV from cosmic muons. The energy deposit of muons on the BGO was calculated using the ASIM mass model. A diode, emitting light through the BGO, was used to mimic the energy deposits corresponding to muons sending a new light pulse every 100 µs. Figure S20 in Supporting Information S1 shows the first count after an energetic pulse created by the diode or a muon. We clearly see the muon/diode peak decreasing to a lower energy channel as the time after the energetic pulse gets smaller. Note also that the effect in the lower energy counts in channels 200 and 400 is less significant, but at dt ≈ 5 µs and dt ≈ 7 µs the counts are measured in higher energy channels. The voltage drop in the PMT following an energetic count leads to underestimated pulse heights for high energies and overestimated pulse heights for low energies. We do not know the behavior for counts with energies between 2.6 and 30 MeV. In principle, a correction function dependent on E 0 , dt, and the energy of the count following the count with energy E 0 , could be defined. However, Figure S19c and S19d in Supporting Information S1 show that an energy band of 30 MeV is compressed to an energy band of ∼1 MeV at dt ≈ 2 µs. A correction function going from an energy band of 1-30 MeV must be very accurate to be trusted to restore the original pulse height. We do not have enough data points and calibration lines in the energy band between 1,275 keV and 30 MeV to make such a correction function reliable.
Instead of a correction function we account for this instrumental effect by implementing an energy dependent "safety time". We define the safety time as the time we have to wait after an energetic pulse for the proton peak to be above 24.8 MeV (within 20% of 31 MeV), see dashed magenta line in Figure S19 in Supporting Information S1. As the on-ground experiment ( Figure S20 in Supporting Information S1) showed that the high energy counts at 32 MeV are most affected by the voltage drop, we know that the energy measurement accuracy at the safety time is ≤20% . A plot of the safety time is shown in Figure S21 in Supporting Information S1.
The safety time presumably does not bias the spectral analysis as it is part of the forward modeling of instrumental effects, as well as implemented on the observed TGF data. We do not remove counts dependent on their energy, but we do remove counts dependent on the energy of the previous count in the same detector. If there exists an unknown significant correlation between the energy of consecutive counts in TGFs at source, then the safety time criteria may bias our spectral analysis as the only time dependence we model is the time profile caused by the transport through the atmosphere added to a Gaussian time distribution that is fitted to each individual TGF. The alternative to the safety time criteria is a fixed dead time of ∼ 30 µs which would remove too many counts during TGF detection.

Appendix B: TGF Beaming Type Modeling
In the literature (Dwyer & Smith, 2005;Gjesteland et al., 2010Gjesteland et al., , 2011Hazelton et al., 2009;Mailyan et al., 2016Mailyan et al., , 2019Østgaard et al., 2008), the angular distribution of the TGF beam at the source is usually modeled using two possibilities: 1. An isotropic distribution within a given angle range, that is parameterized as a cone half angle θ cone . 2. A Gaussian distribution that is parameterized with σ θ .
In the Gaussian beaming the direction of the photons are given by x, y, z, where z is upwards and equal to one, and x and y are independently sampled from a Gaussian distribution. Even if the Gaussian beaming could be considered more realistic (see e.g., Hazelton et al., 2009), the half angle is easier to understand (visualize) and leads to quite close results after propagation through the atmosphere, due to the angular scattering of the photons. For any θ cone , there is an equivalent σ θ for which the fluence distribution after atmospheric propagation is quite close. In Figure B1, we present the equivalence between σ θ and θ cone . To produce this plot, the TGF propagation code was run for a series of σ θ and θ cone , and we calculated for which values the resulting fluence (photons/cm 2 ) profile, as a function of radial distance, is as close as possible (using a maximum likelihood evaluation).

Data Availability Statement
ASIM data are available at the ASIM Science Data Center (https://asdc.space.dtu.dk). Additional data for this paper are available at https://doi.org/10.5281/zenodo.4882745. The library of simulated TGFs are available as time-energy matrices and fluence estimations for different altitudes, beaming and radial distances at https://doi. org/10.5281/zenodo.5493579. The TGF propagation code is available at https://doi.org/10.5281/zenodo.5493589. The TGF simulations were performed on resources provided by UNINETT Sigma2-the National Infrastructure for High Performance Computing and Data Storage in Norway, under project no. NN9526K.