The following article is Open access

Multiparameter Constraints on Empirical Infrasound Period-Yield Relations for Bolides and Implications for Planetary Defense

, , , , , , , , ,

Published 2025 June 19 © 2025. The Author(s). Published by the American Astronomical Society.
, , Citation Elizabeth A. Silber et al 2025 AJ 170 38DOI 10.3847/1538-3881/add47d

PDF Opens in a new tab.
ePub

You need an eReader or compatible software to experience the benefits of the ePub3 file format.

1538-3881/170/1/38

Abstract

How effective are methods for estimating bolide energies from infrasound signal period-yield relationships? A single global period–energy relation can obscure significant variability introduced by parameters such as the atmospheric Doppler wind profile and the bolide’s energy deposition profile as a function of altitude. Bolide speed, entry angle, burst altitude, and multiepisode fragmentation may all play a role in defining the detected period of the shockwave. By leveraging bolide light-curve data from the Center for Near Earth Object Studies, we re-examined the period–energy relation as a function of these parameters. Through a bootstrap approach, we show that various event subsets can deviate from widely cited period–energy models and we identify which specific conditions most strongly reshape the period–energy scaling. The results define both the fidelity and reliability of period–energy relations when no additional data beyond the infrasound record is available and improve the outcome when supporting data from bolide trajectories and light curves are included. Ultimately, these findings expand the scope of earlier models, providing a nuanced and robust framework for infrasound-only yield estimation under a range of bolide scenarios.

Export citation and abstractBibTeXRIS

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

1. Introduction

The need for simple, practical estimators of bolide energy using infrasound signals stretches back to methods originally developed for monitoring nuclear tests, most notably the Air Force Technical Applications Center (AFTAC) period-yield relations (D. O. Revelle 1997). The AFTAC relations were established based on an analysis of the dominant period of nuclear-explosion signals at source-to-station distances of 1300–8500 km (S. Glasstone & P. J. Dolan 1977; D. O. Revelle 1997).

The general form for period-based energy relations is:

Equation (1)

where E is energy in units of kilotons of TNT equivalent (1 ktTNT = 4.184·1012 J), P is the dominant infrasound signal period (in seconds), typically measured at maximum amplitude, and A and B are the regression coefficients. Some relations employ an averaged signal period ($\mathop{P}\limits^{\unicode{773}}$), calculated by averaging the measured period across all stations that detected a particular high-energy event. Table 1 lists all empirical period-yield relations published to date, including those derived in this study.

Table 1. Compilation of Empirical Period-yield Relation Coefficients

ABP*UnitsSourceNotes
3.30 (±0.14)−1.89 (±0.10) $P$ ktThis study
3.71 (±0.42)−2.07 (±0.21) $\mathop{P}\limits^{\unicode{773}}$ ktThis study
3.44−2.57 $\mathop{P}\limits^{\unicode{773}}$ ktR. W. Whitaker (2023)
3.68−1.99 $P$ ktN. Gi & P. Brown (2017)
3.84−2.21 $\mathop{P}\limits^{\unicode{773}}$ ktN. Gi & P. Brown (2017)
3.750.50 $P$ tT. A. Ens et al. (2012)Original form
3.75−2.50 $P$ ktT. A. Ens et al. (2012)Units corrected to kt
3.280.71 $\mathop{P}\limits^{\unicode{773}}$ tT. A. Ens et al. (2012)Original form
3.28−2.29 $\mathop{P}\limits^{\unicode{773}}$ ktT. A. Ens et al. (2012)Units corrected to kt
3.26−2.40 $P$ ktR. Whitaker & J. Mutschlecner (2006)
3.34−2.28 $P$ ktD. O. Revelle (1997)E < 200 kt
4.14−3.31 $P$ ktD. O. Revelle (1997)E > 80 kt

Note. The original T. A. Ens et al. (2012) relations were given in tons of TNT; here, we also present the corrected expressions in kilotons and address minor typographical errors that have appeared in the literature. P* represents the signal period at each station (P) or the signal period averaged ($\mathop{P}\limits^{\unicode{773}}$) across all stations that detected a given event.

Download table as:  ASCIITypeset image

As applied to bolides, ReVelle’s approach effectively reduced the total atmospheric energy by half, on the assumption that a significant fraction of the energy was radiated rather than contributing to the shock. This is reasonably in line with the energy partition estimates for meteoroid entries (M. F. Romig 1965; E. A. Silber et al. 2018). S. Glasstone & P. J. Dolan (1977) state that this 50% partitioning only applies to bursts with altitudes <12 km and that at higher altitudes the fraction is lower. They also noted that the distribution between thermal radiation and blast changes more rapidly above 30 km, which could affect the energy–period relationship for events at very high altitudes. It is not well known, however, how this might apply to bolides. In Table 1, the first relations published by D. O. Revelle (1997) are modified from their original form where the energy had been halved. In the original form, the coefficient B is 0.3 smaller due to halving of the energy (i.e., the log(2) = 0.3). In this paper, all logs are assumed to be base 10.

R. Whitaker & J. Mutschlecner (2006) used unclassified atmospheric nuclear-explosion test data and derived a new relation, applicable to all energy ranges. This relation was recently updated by including an additional unclassified nuclear-explosion data set (R. W. Whitaker 2023). These relations, because they are based on nuclear tests, assume a point-source explosion, which differs from shock production by asteroids and large meteoroids (a cylindrical line source and/or a quasi-spherical source).

Studies by T. A. Ens et al. (2012) and N. Gi & P. Brown (2017) built on point-source explosions to produce empirical fits calibrated to bolide observations. Both of these studies utilized bolide events listed on the Center for Near-Earth Object Studies (CNEOS) fireball database, maintained by the National Aeronautics and Space Administration’s Jet Propulsion Laboratory. These bolides were detected by US government (USG) sensors from space (E. Tagliaferri et al. 1994), providing details such as location, time, and altitude of peak brightness. In some cases, the velocity vector is also available; this allows a tentative estimation of the heliocentric orbit (E. Peña-Asensio et al. 2022). Importantly, the CNEOS database includes total radiated energy (in joules) and impact energy (in kilotons) derived from flashes produced by bolides in the atmosphere and observed by USG space-based instruments in the silicon bandpass (E. Tagliaferri et al. 1994; P. Brown et al. 2002). The USG-derived energy estimates have been found to be robust (J. Borovička et al. 2015; H. A. R. Devillepoix et al. 2019; K. S. Wisniewski et al. 2024), thereby providing an important benchmark for validating and refining energy estimates derived from empirical infrasound-based relations. Currently, the light curves for the bolides that occurred prior to and including 2022 have been made available.

To derive period-based energy relations, T. A. Ens et al. (2012) compiled 63 CNEOS bolide events detected by 113 infrasound stations globally. They also used period measurements from multistation detections of a single event to show that averaging can reduce the scatter of data points (relative to the line defined by these relations). Their relations are updated here from their original form (Table 1) to reflect yields in kilotons and to correct minor typographical errors present in the published literature. The correction from tons to kilotons was done by reducing the original coefficient B by 3 (i.e., log(1000) = 3). Remarkably, the relation with averaged periods was very similar to that by D. O. Revelle (1997).

Drawing on data from 78 CNEOS bolides detected by 179 stations, N. Gi & P. Brown (2017) established two additional relations, one that applies to all detections, and another that uses averaged periods for multistation detections (Table 1). A portion of their data set included bolides from T. A. Ens et al. (2012). Their fits from this larger data set, however, are notably different, illustrating the need to rigorously evaluate the sensitivity of periods to energy and examine the robustness of period-based energy relations when applied to global bolide events. Although N. Gi & P. Brown (2017) examined a few bolides in greater detail by incorporating their light curves, neither they nor earlier studies integrated light-curve analyses across a large number of events to directly connect energy-deposition modes with inferred energy estimates.

The fits from all these period-based energy relations are shown in Figure 1. As with any empirical regression, these relations differ in slope and intercept due to differing calibration data sets, some focusing on smaller bolide events, others on nuclear explosions, and still others on partially idealized assumptions about the source function. Depending on the signal period length and the selected energy relation, the discrepancy in estimated yields can range from minor to over an order of magnitude. This outcome shows that both the choice of period–energy relation and whether data originate from single or multiple stations play significant roles in determining the final yield.

Figure 1. Refer to the following caption and surrounding text.

Figure 1. Log–log plot of energy vs. signal period comparing published period–energy relations. Each line represents a different empirical fit: D. O. Revelle (1997), R. Whitaker & J. Mutschlecner (2006), R. W. Whitaker (2023), and bolide-based relations from T. A. Ens et al. (2012) and N. Gi & P. Brown (2017), the latter two of which include both single-station and multistation (averaged) data. T. A. Ens et al. (2012) and N. Gi & P. Brown (2017) calibrated the bolide energy relations against the USG sensor data, while the other relations stem from calibration against nuclear-explosion data.

Standard image High-resolution image

Many of the early infrasound-based yield estimations also incorporated amplitude relations, acknowledging that signal amplitude could, in principle, map to explosive source energy. However, amplitude-based fits typically demand substantial propagation corrections, particularly for range and atmospheric conditions, rendering them more sensitive to local winds, turbulence, and attenuation. By contrast, it can be argued that the signal period is more robust against these propagation effects, at least in the far field, and thus has served as the principal focus for most modern bolide energy studies. Accordingly, despite the inherent scatter in station-by-station period measurements, period-based relationships have remained the mainstay for rapid yield estimations.

The amplitude-based formulae are beyond the central scope of our paper. We include them for completeness and transparency in Appendix A (Tables A13). The tables consolidate the sometimes disparate approaches and correct for minor unit inconsistencies or typographical misprints that have appeared in past literature.

Over the last few decades, the empirical period–energy relations have formed the backbone of infrasound-based yield estimation, with the nuclear-explosion-derived relations (D. O. Revelle 1997) historically serving as a simple and widely used benchmark. Because these period–energy fits require no detailed knowledge of event geometry or atmospheric profiles, they have remained attractive, despite recognized shortcomings such as station-dependent period variations (T. A. Ens et al. 2012). While subsequent studies (e.g., T. A. Ens et al. 2012; N. Gi & P. Brown 2017) have expanded the event catalogs and offered more robust multistation period averaging approaches, a universal “global” regression still fails to adequately capture the complexity arising from more nuanced effects from atmospheric propagation and varied bolide deposition profiles. These complexities also reflect a broad range of meteorite properties, as meter-scale objects that generate noticeable infrasound signals frequently arrive at Earth already weakened by prior collisions (J. M. Trigo-Rodriguez & J. Blum 2009; E. Beitz et al. 2016). Fragmentation can lead to increased ablation due to a higher surface/volume ratio.

Bolides can produce infrasound signals whose measured periods differ considerably from station to station, in part because each receiver can sample a distinct segment along the bolide path, each with potentially different effective burst altitudes or local fragmentation episodes (E. A. Silber et al. 2009; E. A. Silber 2025). In principle, this would suggest the need to partition events not only by standard parameters such as velocity or altitude but also by more nuanced factors such as the number of fragmentation episodes in a single event, the density or dynamic strength, and the entry geometry (e.g., altitude, entry angle). Earlier works largely acknowledged this complexity (e.g., T. A. Ens et al. 2012; N. Gi & P. Brown 2017), but typically did not incorporate it into a unified regression framework.

Amplitude-based relations for both artificial explosive sources (e.g., A. D. Pierce et al. 1971; J. Stevens et al. 2002; J. P. Mutschlecner & R. W. Whitaker 2009) and bolides (W. N. Edwards et al. 2006; T. A. Ens et al. 2012) have also been developed, yet they suffer even more strongly from range and atmospheric effects (e.g., J. Stevens et al. 2002) and thus remain less reliable than period-based yield estimations. In light of these challenges, there is an obvious need to broaden the analysis of period–energy correlations using a more advanced framework that accounts for geometry (angle, altitude), physical parameters of bolides (mass, diameter), and light-curve classification (single versus multiple “bursts”).

When a meteoroid or an asteroid enters Earth’s atmosphere at velocities of tens of kilometers per second (11.2–∼72 km s−1), it travels far above the local speed of sound. As it descends into denser atmospheric layers, it becomes luminous and forms a shock wave (V. A. Bronshten 1983; Z. Ceplecha et al. 1998; E. A. Silber et al. 2018). Provided the object largely remains intact, this shock behaves much like a cylindrical line source because the hypersonic bolide continuously deposits energy along its path. In such a regime, the shock front extends roughly perpendicular to the trajectory and the pressure disturbances can be approximated as a cylindrical wave expanding outward from the flight path (M. N. Plooster 1970; D. O. ReVelle 1976). A useful parameter here is the blast radius (R0) (M. Tsikulin 1970), which defines the radial distance from the bolide trajectory within which the atmospheric pressure is significantly perturbed by the shock (A. A. Few 1969). Closer to the meteoroid, the shock is strongly nonlinear (i.e., overpressures exceed ambient pressure by a large factor), whereas beyond the blast radius the wave transitions into a more linear or weakly nonlinear regime (D. O. ReVelle 1976).

As a first approximation, higher energy deposition yields a larger blast radius, shifting the shock signature toward lower frequencies (longer periods). However, in practice, the picture is more complicated. Atmospheric ablation can alter meteoroid flow regimes (M. Moreno-Ibáñez et al. 2018; E. A. Silber et al. 2018) and fragmentation events release energy in short segments or discrete “bursts” (e.g., J. M. Trigo-Rodríguez et al. 2021), which cause the shock wave to locally transition from a cylindrical (line source) to a quasi-spherical (point source) geometry (E. A. Silber & P. Brown 2019). In a continuous flight with no major fragmentation, the shock predominantly resembles a cylindrical wave trailing behind the bolide. In contrast, an abrupt fragmentation (or airburst) produces a dominant near-spherical expansion, and multifragmentation events generate multiple overlapping quasi-spherical shock regions along the trajectory, complicating how the overall energy is distributed in time and space. In some instances, the observed light curve represents the sum of light curves generated by individual “daughter” objects produced during a cascading fragmentation event. We discuss these complexities further in Section 4, in the context of our results.

The bolide’s associated light curve can be diagnostic in identifying these processes because its shape, duration, and peak structure are sensitive to fragmentation events and can be leveraged to infer critical details about energy release and shock formation (E. A. Silber & P. G. Brown 2014; E. A. Silber 2024b; T. C. Wilson et al. 2025). Figure 2 represents an idealized schematic illustrating how a bolide’s fragmentation behavior can lead to characteristic light-curve signatures. We omit a simple, nonfragmenting asteroid case from our diagram. On the left, four conceptual light curves are shown, each representing a distinct breakup pattern: an airburst with a single narrow peak, a single-point fragmentation with one broader peak, a discrete multipoint fragmentation with a series of peaks, and a continuous multipoint fragmentation producing multiple overlapping peaks. To the right, each fragmentation type is shown with both steep and shallow-entry trajectories, emphasizing how the bolide’s path through the atmosphere influences the spatial distribution of fragmentation episodes and the accompanying shock formation. In a steep entry, the energy is deposited within a relatively narrow vertical column, whereas a shallow entry spreads the disruption over a longer atmospheric path. These variations in breakup processes, combined with differences in entry angles, result in the diverse light-curve shapes observed and offer clues about the shock generation and total energy release of the bolide event. Moreover, information gleaned from simultaneous infrasound and light-curve analysis is of extraordinary importance for improving our understanding of the atmospheric entry of rare superbolides and assessing their potential to become hazardous projectiles (J. M. Trigo-Rodríguez 2022).

Figure 2. Refer to the following caption and surrounding text.

Figure 2. Idealized diagram, not drawn to scale and intended only as a rough conceptual representation, illustrates how a bolide’s fragmentation behavior affects the resulting light curve. The left column shows four conceptual light curves, from a single-peaked airburst to continuous multipoint fragmentation, while the right panels depict how different geometry scenarios might look for steep vs. shallow entries.

Standard image High-resolution image

Here, we investigate the impact of these processes in the period-yield relationships using the CNEOS bolide light-curve data. One of the primary objectives in this study was to infer the mode of shock deposition through the analysis of light curves, since they serve as a diagnostic for various physical processes taking place during a bolide entry. We aim to investigate possible links between types of shock, altitudes of fragmentation episodes, and bolide parameters, and feed these into the energy fits to investigate the linkage. To evaluate the robustness of the period-yield relations, we apply a bootstrap model, a powerful statistical technique that repeatedly resamples the data to produce multiple “bootstrapped” data sets (B. Efron 1992; B. Efron & R. J. Tibshirani 1993). Each resampling run estimates the parameters of interest (e.g., slope and intercept in period–energy fits) and the collection of those estimates is used to derive robust confidence intervals. The bootstrap approach offers a more comprehensive view of the variability and reliability of the derived coefficients than a single, static fit. In principle, testing the reliability of period-based energy relations under diverse event conditions will allow an assessment of the parameter ranges in which these relations are most reliable, or where they fall short.

2. Methods

2.1. Bolide Data Set

A bolide infrasound detections database, including measured signal parameters (see Appendix B), was compiled by drawing on multiple published sources that focused on bolides detected by US government sensors and subsequently listed on the CNEOS website (E. A. Silber et al. 2011; T. A. Ens et al. 2012; N. Gi & P. Brown 2017; T. Ott et al. 2019; C. Pilger et al. 2020; T. Ott et al. 2021; P. Hupe et al. 2024; E. Silber et al. 2024; Figure 3). The CNEOS database supplies contextual metadata to reported energy yields, such as event time, location, and frequently a velocity vector and corresponding light curve data (e.g., E. A. Silber 2024b; E. A. Silber & V. Sawal 2025). Where velocity vectors were available (95 events), we derived entry angles, azimuths, and orbital parameters based on the method of E. Peña-Asensio et al. (2022), which then informed about each bolide’s physical properties.

Figure 3. Refer to the following caption and surrounding text.

Figure 3. Global map showing all 138 bolides (red circles) investigated in this work and 59 distinct infrasound stations that recorded detections (gray triangles).

Standard image High-resolution image

All but 10 bolides had infrasound signals detected by the stations of the International Monitoring System (IMS) of the Preparatory Commission for the Comprehensive Nuclear-Test-Ban Treaty Organization, from which signal periods were derived. To ensure consistency and verify data quality, we recalculated source-to-station distances, back azimuths (direction of infrasound wave arrival at a station), and celerities (defined here as the average acoustic propagation speed from source to receiver). We discarded any entries deemed outliers, most notably those with anomalously large back azimuth discrepancies (where the difference between observed and theoretical values is >50°) or unrealistically high celerities (>360 m s−1). In keeping with established energy-relation derivations in prior studies, we also excluded signals that traveled through the thermospheric waveguide (celerity <220 m s−1 (P. T. Negraru et al. 2011)). Therefore, all events have celerities between 220 and 360 m s−1 and back azimuth deviations <50°.

This refinement process yielded 138 distinct bolide events and 362 total detections (Figure 3), constituting the largest single consolidated data set of infrasound-detected CNEOS bolides to date. Among these, 69 bolides were detected by one station globally, whereas the remaining 69 were observed across two or more infrasound stations worldwide. The number of distinct stations that made detections was 59. Out of these, 53 are IMS stations.

For these bolides, we retrieved event-specific atmospheric profiles from the Ground-2-Space model hosted by University of Mississippi (D. P. Drob et al. 2003; C. H. Hetzer 2024) along each source–station path. Earlier work has shown that Doppler shifts due to stratospheric winds can alter the observed period by more than ∼10% (T. A. Ens et al. 2012). However, they also noted that uncertainties in wind velocity can be as large as the mean itself, and therefore did not correct for Doppler shift. In this study, we explore both raw and Doppler-corrected measured periods. We make Doppler shift corrections following the approach outlined in D. O. ReVelle (2010). Since infrasound propagation is assumed to be dominated by the stratospheric waveguide, we consider the altitude span of 40–60 km in a range-dependent fashion, computing mean and standard deviation values of the wind field. We also accounted for temporal variations, which is particularly relevant for long-distance infrasound paths requiring multihour propagation.

Figure 4 provides histograms illustrating the distributions of bolide parameters available in our data set, including (top row) peak-brightness altitude, source-to-station distance, and entry velocity; (middle row) total energy (on a log scale), entry angle, and impactor diameter; and (bottom row) dynamic strength (also log-scaled), impactor density, and the Tisserand parameter. The Tisserand parameter is a near-constant measure in the restricted three-body problem (e.g., Sun–planet–small body) that characterizes how a smaller object’s orbit changes, or fails to change, before and after a close planetary encounter. In practical terms, the Tisserand parameter computed with respect to a given planet helps categorize interplanetary bodies (e.g., comets versus asteroids) by capturing relationships among orbital elements, such as semimajor axis, eccentricity, and inclination. For instance, many Jupiter-family comets have Tisserand values between 2 and 3 (C. D. Murray & S. F. Dermott 1999). Not all events in the database have complete information for these parameters, which is reflected in some histograms having lower counts than others. Nonetheless, these plots demonstrate the range and diversity of bolide characteristics captured in our final data set.

Figure 4. Refer to the following caption and surrounding text.

Figure 4. Histograms showing the distributions of bolide parameters available in our data set: peak-brightness altitude, source-to-station distance, entry velocity, total impact energy (on a log scale), entry angle, impactor diameter, dynamic strength (also log-scaled), impactor density, and the Tisserand parameter.

Standard image High-resolution image

The CNEOS database hosts an inventory of light curves (LCs) associated with most bolide events up to and including most of 2022. Presently, these come in a PDF format, requiring digitization. We gathered all available light curves, and then digitized and processed them following the analysis and classification algorithm developed by E. A. Silber & V. Sawal (2025). This algorithm, Bolide Light-curve Analysis and Discrimination Explorer, uses the Savitzky–Golay filtering, prominence-based peak detection, and gradient analysis to automate the identification of fragmentation events and, therefore, energy release modes. E. A. Silber & V. Sawal (2025) devised five classification schemes based on the light curve features to extract the dominant process of energy deposition (e.g., airburst, discrete fragmentation, and continuous fragmentation). For the purposes of this study, we adopted a simplified light curve classification scheme: events with a single, dominant fragmentation episode (LC class 1) versus those exhibiting multiple fragmentation episodes, whether continuous or discrete (LC class 2). This simplified categorization enabled us to probe whether repeated fragmentation episodes (LC class 2) might broaden the infrasound signal relative to single-burst events (LC class 1), thereby clarifying fragmentation-driven influences on the measured period. Only three events did not have the light curve available.

We analyzed the light curves in two steps. First, we applied our classification algorithm to all available light curves to determine their respective types. Second, for CNEOS events that provided both a velocity vector and a peak-brightness altitude, we anchored the altitude to the peak-brightness point identified in the light curve and extrapolated the altitude range from there. Here, as the first approximation, we assume constant velocity. This approach allowed us to examine the onset or conclusion of fragmentation events more thoroughly, rather than focusing solely on the peak-brightness altitude. A representative example of this step in our analysis is shown in Figure 5. Here, we compare two distinct bolide fragmentation scenarios by plotting light curve intensity (in W/Sr) against both time and altitude. These examples demonstrate our approach to discerned multifragmentation versus a single, abrupt burst by analyzing the light curve shape and its corresponding altitude profile. As noted earlier, the fragmentation can be rather complex, where multiple fragmentation episodes can be discrete or continuous (see Figure 2). In the latter case, this could result in persistent emissions along the wake. However, we do not make that specific differentiation here, but rather focus on single burst versus multibursts, regardless of how they might have occurred.

Figure 5. Refer to the following caption and surrounding text.

Figure 5. Representative examples of different classes of light curves, plotted as a function of the USG-provided time (left) and calculated altitude assuming no deceleration (right). Panels (a) and (b) show LC class 1, depicting a single dominant fragmentation event, where a single abrupt peak unfolds within 2 s and spans only about 4 km in altitude (37–33 km). This single fragmentation event is consistent with a classic “airburst” in which the entire object is abruptly vaporized due to the fragmentation. Panels (c) and (d) show LC class 2. This is a multifragmentation event, with multiple bursts across 4 s and an altitude range from 48 km down to 35 km, indicative of energy deposition over 13 km.

Standard image High-resolution image

2.2. Ordinary Least Squares Partitioning

In their study, T. A. Ens et al. (2012) employed a multivariate approach that factored in parameters such as range and amplitude; however, their analysis did not account for altitude or light curve classification, both of which can critically alter where and how shock energy is deposited along the bolide trajectory. Moreover, amplitude-based fits typically require explicit propagation corrections, making range an even more central parameter. In our study, by contrast, we focus on period rather than amplitude, which lessens the need for direct range corrections and enables us to isolate other physical factors (e.g., altitude, multiepisode fragmentation) that might influence the period–energy relationship.

Fundamentally, we adopt a log–log regression model of the form: $\mathrm{log}\left(P\right)=m\mathrm{log}\left(E\right)+c$. Recognizing that infrasound detections can be influenced by source altitude, entry angle, velocity, and similar parameters, we analyze not only the data set as a whole but also partition it into subsets. For instance, we might isolate only events below a certain altitude threshold or within a specific velocity range; likewise, we can set distance ≤1000 km or altitude <20 km. Each such partition defines a data set subset for which we perform a distinct ordinary least squares (OLS) fit. This step mitigates confounding effects (D. C. Montgomery et al. 2021) by preventing a single, global fit from obscuring systematic differences in event geometry or fragmentation type. Moreover, by initially running OLS on narrowly defined partitions and then applying a bootstrap framework, we can gauge more transparently how each parameter impacts the resulting slope and intercept. As in previous studies, we pursue two approaches in our analysis: (1) treating each station detection independently, and (2) computing a mean signal period for bolide events with multiple stations, by averaging the individual detections across those stations.

Within our database, each event at minimum includes the signal period, the total impact energy reported by CNEOS (in kt), and the source-to-station distance (in kilometers); in all but one case, a peak-brightness altitude (in kilometers) is also available. A subset of 135 bolides includes a light curve classification, and 95 bolides contain additional parameters (i.e., entry angle, velocity, mass, density, and diameter). Rows missing a specific parameter are excluded only from analyses involving that variable, thus preserving the maximum number of valid records for partitions that do not rely on it. For each partition (e.g., distance ≤1000 km, altitude <20 km, velocity <20 km s−1, etc), we perform fitting (D. C. Montgomery et al. 2021) using the statsmodels.api.OLS routine in Python, yielding a slope $\hat{m}$ and intercept $\hat{c}$. We also apply bin permutations of up to five parameters in the analysis. Since many published energy–period formulas use the inverted form: $\mathrm{log}\left(E\right)\,=A\mathrm{log}\left(P\right)+B$, we convert via $A=1/\hat{m}$ and $B\,=1/\hat{c}$, enabling direct comparisons with reference relations in the literature (e.g., T. A. Ens et al. 2012; N. Gi & P. Brown 2017; and D. O. Revelle 1997). A single OLS pass produces one pair $\left(\hat{m},\hat{c}\right)$, and hence one pair $\left(A,B\right)$. However, we also seek to characterize uncertainties, given that small sample sizes or physically diverse events can cause significant variance in slope estimates. Consequently, a bootstrap framework is employed to estimate confidence intervals for ($\hat{m},\hat{c},A,B)$, hereby assessing the reliability of each partition’s fit under repeated sampling.

2.3. Bootstrap Approach for Period–Energy Fitting

To robustly estimate the variability and confidence intervals of our fitted slope/intercept pairs, and their inverted counterparts, we employ a bootstrap sampling framework (B. Efron 1981, 1992; B. Efron & R. J. Tibshirani 1993; T. Hastie et al. 2009). For each partition of the data set, we draw the same number of rows as in the original subset with replacement, repeatedly generating “pseudo-data sets.” We then fit the OLS model $\mathrm{log}\left(P\right)={m}_{{i}}\mathrm{log}\left(E\right)+{c}_{i}$ on each pseudo-data set and invert the resulting slope and intercept. By iterating this process n times (n = 2000), we accumulate sets $\left\{{m}_{i},{c}_{i}\right\}$ and $\left\{{A}_{i}=\frac{1}{{m}_{i}},\,{B}_{i}=\frac{1}{{c}_{i}}\right\}$, corresponding to each bootstrap iteration $i$. From these distributions, we obtain not only the means $\left(\mathop{A}\limits^{\unicode{773}},\,\mathop{B}\limits^{\unicode{773}}\right)$ but also standard deviations and 95% confidence intervals.

A key point is that we transform each pair $\left({m}_{i},{c}_{i}\right)$ individually into $\left({A}_{i},{B}_{i}\right)$. This captures the true variability of the inverted parameters $\mathrm{log}\left(E\right)-\mathrm{log}(P)$, rather than merely inverting the mean slope $\mathop{m}\limits^{\unicode{773}}$. Indeed, $\mathop{A}\limits^{\unicode{773}}\,\ne 1/\hat{m}$ once distributions are nonlinear. After n bootstrap iterations, we collect distributions $\left\{{m}_{i}\right\}$, $\left\{{c}_{i}\right\}$, $\left\{{A}_{i}\right\}$, $\left\{{B}_{i}\right\}$. We then compute the mean values $\mathop{m}\limits^{\unicode{773}}$, $\mathop{c}\limits^{\unicode{773}}$, $\mathop{A}\limits^{\unicode{773}}$, $\mathop{B}\limits^{\unicode{773}}$, along with their standard deviations. Moreover, we derive confidence intervals in two ways: (i) normal approximation ($\mathop{A}\limits^{\unicode{773}}\pm z{\sigma }_{A}$ and $\mathop{B}\limits^{\unicode{773}}\pm z{\sigma }_{B}$ for 95% intervals, where z = 1.96), which is appropriate if $\left\{{A}_{i}\right\}$ and $\left\{{B}_{i}\right\}$ are roughly symmetric and (ii) percentile-based, where we take the empirical 2.5th and 97.5th percentiles of $\left\{{A}_{i}\right\}$ and $\left\{{B}_{i}\right\}$. This is often more robust if the distributions are skewed or if normality is not a good approximation. Comparing these two intervals helps us assess the degree of symmetry or skew in the underlying distributions.

To discern how specific parameters, such as source-to-station distance, peak-brightness altitude, entry angle, impact velocity, or impactor properties, affect the period–energy relation, we subdivide the data set into bins. For instance, distance bins might include (d ≤ 1000, ≤ 3000, ≤ 5000, ≤ 10,000, all) while altitude bins include (zpb < 20, < 30, < 40, all). We then run separate bootstrap-based OLS fits for each bin (or bin combination), storing: (1) the number of data points N in that subset, (2) the mean $\mathop{m}\limits^{\unicode{773}}$, $\mathop{c}\limits^{\unicode{773}}$, $\mathop{A}\limits^{\unicode{773}}$, $\mathop{B}\limits^{\unicode{773}}$ with their standard deviations and confidence intervals, and (3) a diagnostic figure comparing the resulting line to existing references. A similar procedure applies to partitions defined by entry angle (<30°, 30°−60°, >60°), mass bins, density ranges, velocity, diameter thresholds, and light-curve class classifications. This partition-based, multiparameter nested approach systematically tests whether, for instance, shallow-entry versus steep-entry bolides or high-altitude versus low-altitude bursts exhibit distinct slope–intercept behaviors. By combining partitioning with a bootstrap framework, we obtain a clearer, more statistically robust picture of how each subset’s period–energy fit varies, and the degree of confidence we can place in that variation.

3. Results

3.1. Ordinary Least Squares Results

Table 2 presents the results from the OLS analysis for different subset of bolides. N is the number of subset samples. The fit parameters A and B are given from both uncorrected and Doppler-wind corrected data. The classical regression coefficient is given also. In single detections, the coefficients A and B across 31 primary populations, binned by up to two parameters, were 3.47 ± 0.49 and −1.88 ± 0.24, respectively. Specifically, the range was 2.60 ≤ A ≤ 5.03, and −2.26 ≤ B ≤ −1.18. While the absolute differences in the regression coefficients A and B may seem modest in numerical terms, they reflect statistically and physically significant deviations among partitioned populations.

Table 2. Summary of Regression Coefficients A and B Derived Through OLS Across all Detections

     Doppler Corrected 
DistanceNABr2ABr2
(km)       
global3623.29−1.880.593.28−1.880.58
d = all, z < 30 km1632.88−1.520.672.87−1.520.67
d = all, z < 40 km2823.15−1.770.633.13−1.760.62
d 1000 km, z < 30 km113.10−1.600.662.99−1.550.66
d ≤ 1000 km, z < 40 km273.23−1.700.553.18−1.700.52
d ≤ 1000 km, z = all km373.60−1.890.473.60−1.920.43
d ≤ 3000 km, z < 30 km633.55−1.780.573.45−1.720.58
d ≤ 3000 km, z < 40 km1413.88−2.040.473.82−2.010.46
d ≤ 3000 km, z = all1984.16−2.220.384.11−2.190.37
d ≤ 5000 km, z < 30 km953.24−1.680.603.11−1.610.61
d ≤ 5000 km, z < 40 km2063.59−1.960.503.50−1.920.50
d ≤ 5000 km, z = all2793.82−2.110.433.74−2.070.42
d ≤ 10,000 km, z < 30 km1433.01−1.570.652.98−1.560.65
d ≤ 10,000 km, z < 40 km2603.32−1.840.603.29−1.830.59
d ≤ 10,000 km, z = all3403.49−1.970.543.46−1.950.53
θ < 30°, v ≤ 20 km s−1693.00−1.810.722.96−1.760.71
θ < 30°, v = all842.95−1.690.722.92−1.650.72
30°≤θ < 60°, v < 20 km s−1704.09−2.150.333.92−2.150.30
30°≤θ < 60°, v = all884.43−2.230.314.32−2.250.28
θ≥60°, v ≤ 20 km s−11013.59−2.100.573.57−2.070.57
θ≥60°, v = all1073.61−2.110.573.57−2.070.56
θ = all, v ≤ 20 km s−12403.25−1.860.653.26−1.880.62
θ = all, >20 km s−1395.03−2.260.274.91−2.170.30
1e5 ≤ mass < 1e81343.00−1.650.443.05−1.720.42
2000 ≤ dens < 35001713.83−2.190.383.77−2.180.37
dens ≥ 3500832.60−1.180.692.64−1.220.71
1 ≤ diam < 102603.64−2.040.513.64−2.050.49
diam = all2773.27−1.830.623.27−1.840.61
LC = 11563.31−1.920.463.26−1.860.46
LC = 22033.33−1.880.593.30−1.880.58
LC = all3593.29−1.880.593.28−1.870.58

Download table as:  ASCIITypeset image

Table 3 presents the same data but for averaged-period detections. The coefficients A and B across 18 primary populations, binned by up to two parameters, were 3.63 ± 0.37 and −1.92 ± 0.22, respectively. Specifically, the range was 2.85 ≤ A ≤ 4.21, and –2.18 ≤ B ≤ −1.34.

Table 3. Summary of Regression Coefficients A and B Derived Through OLS Across Averaged Period Events

     Doppler Corrected
Distance [km]NABr2ABr2
global1383.65−2.050.453.60−2.010.46
d = all, z < 30 km353.07−1.600.713.05−1.570.71
d = all, z < 40 km943.31−1.780.553.27−1.740.55
θ < 30°, v ≤ 20 km s−1193.80−2.120.543.78−2.100.54
θ < 30°, v = all273.50−1.860.593.51−1.860.59
30° ≤ θ < 60°, v < 20 km s−1364.11−2.050.374.05−2.000.37
30° ≤ θ < 60°, v = all444.21−2.040.374.12−1.990.37
θ ≥ 60°, v ≤ 20 km s−1213.53−1.960.703.51−1.960.68
θ ≥ 60°, v = all243.67−2.100.623.65−2.090.60
θ = all, v ≤ 20 km s−1763.66−1.960.563.62−1.930.55
θ = all, >20 km s−1193.69−1.820.433.69−1.820.44
1e5 ≤ mass < 1e8193.14−1.690.443.13−1.680.43
2000 ≤ dens < 3500674.08−2.160.363.99−2.100.36
dens ≥ 3500152.85−1.340.772.84−1.330.77
1 ≤ diam < 10934.10−2.160.444.05−2.120.44
diam = all963.71−1.970.523.67−1.940.52
LC = 1413.63−2.130.403.61−2.110.40
LC = 2943.95−2.120.413.91−2.090.41
LC = all1353.69−2.050.453.66−2.030.45

Download table as:  ASCIITypeset image

Figure 6 illustrates the distribution of best-fit coefficients A and B obtained from our OLS regressions across a range of bolide subsets. We performed a total of 509 fits for single-station detections, each corresponding to permutations of up to five parameters (e.g., altitude, distance, velocity, fragmentation type, and light curve class), and 126 fits for multistation averages (where distance was excluded). After removing any partition with insufficient data (N < 20 in single detections, N < 10 in averaged periods), we retained 274 and 70 pairs of coefficients, respectively.

Figure 6. Refer to the following caption and surrounding text.

Figure 6. Distributions of OLS-derived coefficients A and B for both single-station (blue points) and multistation averaged (orange points) period fits.

Standard image High-resolution image

Visually, the single-station data (blue points) cluster in a diagonal pattern, reflecting the fact that more negative intercepts (B) occur at higher slopes (A). The orange points from multistation averages appear more tightly grouped, reflecting the often more stable period estimates when multiple stations contribute to a single, averaged period. These results emphasize both the overall trend and the variability inherent in subdividing bolide data by altitude, distance, and other parameters. We also overlaid coefficients from earlier studies alongside our best fits. Across all single-station and averaged-period detections, respectively, our global regression yields:

Equation (2)

Equation (3)

Some bins suggest a strong correlation between period and energy when the shock is observed closer to ground level and is less “smeared” by propagation. Tails at the lower end of the slope distribution arise in scenarios where energy release is more spread out, either by shallow geometry or multifragmentation combined with distant station observation.

Figure 7 illustrates the overall variability and spread across all results in the form of histograms, subdivided by single (Figures 7(a) and (c)) and averaged populations (Figures 7(b) and (d)), and A and B coefficients.

Figure 7. Refer to the following caption and surrounding text.

Figure 7. Histograms showing the spread in coefficients A and B as a function of populations of single and averaged detections. These include all our fits after outliers were removed.

Standard image High-resolution image

Since our data set is substantially larger than in previous studies, we also performed the fitting using Doppler-corrected periods (D. O. ReVelle 1974). Across all single-station and averaged-period detections, respectively, our global regression yields:

Equation (4)

Equation (5)

Compared to the noncorrected data, the resulting coefficients differed by only about 1.2% on average, with a maximum deviation of 4.3% in both coefficients A and B. This outcome is consistent with T. A. Ens et al. (2012), who similarly concluded that Doppler corrections have little effect when deriving period-based relations for a global population. Nevertheless, Doppler effects may still be important for individual-event analyses because the averaging over numerous events can mask or dilute any single-event corrections.

To illustrate this point, we plotted the effect of Doppler shifts on measured infrasound periods (Figure 8). Figure 8(a) shows wind velocities along each source–station path for which wind data was available, with the shaded region highlighting distances of up to 6000 km where wind speeds often exceed ±50 m s−1. Figure 8(b) compares the Doppler-corrected and raw measured periods, demonstrating that most points cluster near the 1:1 line, although some events deviate substantially at longer periods. The absolute difference between the Doppler-corrected and raw periods (Figure 8(c)) reveals that corrections can shift the period by several seconds, especially for greater source–station separations or higher wind speeds. Figure 8(d) “zooms in” to 14 s, and more closely shows that most points lie near the 1:1 line.

Figure 8. Refer to the following caption and surrounding text.

Figure 8. (a) Wind velocity along the source–station path, plotted against propagation distance. (b) Doppler-corrected vs. measured period. (c) Difference between Doppler-corrected and measured periods, showing distance-dependent variations. (d) Zoomed-in version of (b) for shorter periods, again showing how stratospheric winds can alter observed infrasound periods.

Standard image High-resolution image

3.2. Bootstrap Approach Results

Using a bootstrap approach, we derived A and B coefficient pairs with their respective standard deviations and confidence bounds. Similar to the OLS analysis performed prior to this, for single-station detections, we included source-to-station distance as a factor, whereas for multistation average periods we excluded it. Alongside global A and B pairs, we produced over 500 subset fits, each reflecting permutations of up to five parameters in a given subgroup. Before examining the correlations, we discarded any results with N < 20 for single detections and N < 10 for averaged detections. Further to this, we examined coefficient standard deviations and removed poor fits (i.e., large standard deviations), thereby ensuring only robust regressions were carried forward in the analysis. We were left with 250 and 50 pairs of coefficients A and B, respectively, for single detections and averaged detections. The global fits with uncertainties are as follows:

Equation (6)

Equation (7)

Doppler-corrected periods yield coefficients A and B that closely match those in Equations (6) and (7), remaining well within their respective uncertainties.

Figure 9 shows the coefficient A and B pairs for single (navy blue circles) and averaged detections (diamonds). We opted against plotting the error bars because they obscure the plot. For context, we also plotted single detections derived through the OLS method (light green circles). The global fits with error bars are shown with a triangle (Equation (6)) and square (Equation (7)).

Figure 9. Refer to the following caption and surrounding text.

Figure 9. Coefficient A and B pairs for single detections (navy blue circles) and averaged detections (diamonds). Single detections derived with OLS (light green circles) are also shown, but some extreme points are cut off in this zoomed in version. Error bars for individual points are omitted to maintain clarity. The global fits, with error bars, appear as a triangle (Equation (4)) and a square (Equation (5)).

Standard image High-resolution image

For additional context and to compare the outcomes from the two fitting approaches, we plotted the distributions of coefficients A and B under both OLS (hatched bars) and bootstrap (solid bars) approaches (Figure 10). Figures 10(a) and (b) display histograms of coefficients A and B for single detections and Figures 10(c) and (d) for averaged detections. This allows a straightforward visual comparison of how each method influences the coefficient estimates for single versus averaged detections.

Figure 10. Refer to the following caption and surrounding text.

Figure 10. Histograms of coefficients A and B for single detections (panels (a) and (b)) and averaged detections (panels (c) and (d)). The OLS results are shown with hatched bars, while the bootstrap results are shown with solid bars. Note that the axes scale is different in all panels.

Standard image High-resolution image

Figure 11 provides an at-a-glance summary of the coefficient A and B for various subsets of single and averaged detections. Each label indicates a subset’s primary parameter (e.g., altitude, distance, entry angle, light curve class, or yield range), and the corresponding sample size (N). This figure also shows published coefficients from N. Gi & P. Brown (2017), T. A. Ens et al. (2012), and D. O. Revelle (1997), plotted as vertical lines for A and reference lines for B, to place the new results in context. This layout allows a straightforward comparison of how each subset’s best-fit slope (Figure 11(a)) and intercept (Figure 11(b)) compare to previously reported values, motivating an evaluation of how parameters such as altitude, range, or fragmentation behavior might influence the period–energy relationship.

Figure 11. Refer to the following caption and surrounding text.

Figure 11. Coefficients A (a) and B (b) from $\mathrm{log}\left(E\right)=A\mathrm{log}\left(P\right)+B$, plotted for numerous bolide subgroups. Each point represents a distinct partition (e.g., by altitude, angle, distance, or fragmentation type), with error bars indicating fit uncertainties. Reference lines from N. Gi & P. Brown (2017), T. A. Ens et al. (2012), R. Whitaker & J. Mutschlecner (2006), R. W. Whitaker (2023), and D. O. Revelle (1997) are shown for comparison, demonstrating how geometry (e.g., shallow vs. steep angles) and multifragmentation (LC class 1 vs. 2) can shift slopes and intercepts away from classic global values. The circles and diamonds represent fits for single-station and averaged events, respectively.

Standard image High-resolution image

In Figure 12, we show the bolide period-yield regression fits obtained through this study and compared to all earlier studies (Figure 12(a)). Single detection fits and averaged fits are shown in Figures 12(b) and (c), in the context of relevant earlier studies. Finally, the global fits derived through bootstrap method alongside fits for LC1 and LC2 populations are shown in Figure 12(d).

Figure 12. Refer to the following caption and surrounding text.

Figure 12. Log–log plot of total impact energy vs. signal period comparing (a) published period–energy relations and those derived in this study. Fits for ingle station detections (b) and averaged periods (c) are plotted separately alongside earlier studies. In panel (d), we show the global fits alongside LC1 and LC2 (averaged).

Standard image High-resolution image

3.3. Period and Yield versus Multiparameter Variables

We computed Spearman correlation (C. Spearman 1904) matrices to examine how primary parameters relate to one another across our bolide data set. The parameters are period, CNEOS energy, impact velocity, entry angle, altitude of peak brightness, meteoroid density, meteoroid diameter, and meteoroid dynamical strength. The Spearman’s rank correlation coefficient is a nonparametric measure of how two variables covary in a monotonic manner, using rank orders rather than raw values. By forgoing assumptions of linearity and normality, it remains robust in the presence of outliers or skewed data, capturing whether one variable consistently increases or decreases with the other. Figure 13 presents correlation matrices for averaged-period bolides. Figure 13(a) includes the entire data set, while Figures 13(b) and (c) show two primary fragmentation classes, LC class 1 (single-burst or airburst events) and LC class 2 (multifragmentation). Each cell in the matrices is color-coded according to the strength and sign of the correlation, with red indicating a positive correlation and blue a negative correlation.

Figure 13. Refer to the following caption and surrounding text.

Figure 13. Correlation matrices for averaged-period bolides across (a) all events, (b) single-burst (LC class 1) events, and (c) multiburst (LC class 2) events. Each panel shows pairwise Spearman correlation coefficients among various parameters: infrasound period, CNEOS energy, entry angle, velocity, diameter, density, dynamic strength, and peak-brightness altitude.

Standard image High-resolution image

As expected, the highest correlation coefficient is between CNEOS energy and meteoroid diameter. The meteoroid diameter is typically estimated from the reported yield. However, we also see a strong correlation when the diameter is independently measured from cosmogenic nuclides in recovered meteorites (e.g., P. Jenniskens et al. 2009; P. Jenniskens et al. 2021). Another strong correlation exists between dynamical strength and meteoroid density, both parameters are related to the atmospheric altitude at which key phenomena are observed. CNEOS energy and period are less correlated, reflecting the influence of other parameters in the relation.

To explore the relationship between bolide energy and the maximum distance at which the infrasound signal was detected, we plotted the detection range as a function of satellite-estimated yield (Figure 14). Examining this trend provides a straightforward means of assessing how far infrasound from a given bolide might travel before dropping below detectability thresholds. As in previous work (e.g., T. A. Ens et al. 2012), an empirical upper-range curve can be superimposed to represent a typical outer envelope of observed detections, thereby highlighting how wind conditions or other propagation effects can extend (or limit) the observable range for lower-energy events. This approach also offers a simple visualization for distinguishing bolide signals from other atmospheric sources by illustrating whether an observed detection range is consistent with the event’s reported yield. Our expression for detection range, with energy in units of kt and distance in units of km, is:

Equation (8)

Figure 14. Refer to the following caption and surrounding text.

Figure 14. Detection range (km) vs. total impact energy (kt) for all infrasound stations that recorded a positive signal. The red line represents the upper-range envelope (yellow-black circles) derived in this study, while the black line shows the relation from T. A. Ens et al. (2012), $\mathrm{log}\left(R\right)=0.33\mathrm{log}(E)+3.79$. Each blue circle corresponds to a unique detection, illustrating how measured ranges vary with event yield; the superposed curves provide a simple upper bound on the distances over which infrasound from bolides of given energies can be detected.

Standard image High-resolution image

4. Discussion

4.1. OLS versus Bootstrap and Value of Subgroup Binning

OLS provides a straightforward point estimate of the best-fit line under classic linear regression assumptions, while bootstrap resampling quantifies how sensitive those OLS results are to variations in the underlying data. The close agreement between the OLS (Figures 6 and 7) and bootstrap-derived results (Figures 9 and 10) reinforces the stability of our core period–energy relationships, demonstrating that removing or resampling small subsets of events does not drastically alter the estimated slopes (A) or intercepts (B). The bootstrap analysis further indicates where outlier events or small partitions (e.g., bins characterized by shallow angles and low altitudes) can yield slope–intercept pairs that can deviate from the “typical” 3–4 slope region and −2 intercept. While OLS averages over these extremes, the bootstrap approach (see tails in Figures 9 and 10) illuminates their frequency and variability, providing robust confidence measures for each subgroup fit. By applying multiparameter binning, using both OLS and bootstrap, we identify distinct phenomena such as steeper slopes at low altitudes and short ranges, and shallower slopes at shallow angles or distant stations. Our findings (as depicted in Figures 9 and 10) confirm that these variations are physically meaningful and represent distinct shock-physics regimes, thereby supporting the need for subgroup-specific fits rather than relying on a single universal formula.

4.2. Sensitivity of Coefficients A and B

Our expanded data set exhibits significant variability in coefficients A (slope) and B (intercept), reflecting the complex interplay of geometry, fragmentation type, and propagation effects in period–energy relationships. Although our global fits (Equations (2)–(5)) generally conform to earlier literature (e.g., D. O. Revelle 1997; T. A. Ens et al. 2012; N. Gi & P. Brown 2017), partitioning the data by altitude, angle, distance, or multifragmentation yields values of A that may exceed 4.0 or fall to 2.5–3.0, suggesting a wider parameter space than that provided by a single global fit. These variations can be attributed to shock-physics principles. Deeper or steeper trajectories yield stronger coupling between energy and the dominant period, whereas shallower entries or large station distances dilute period sensitivity to energy. Equally important, multifragmentation (LC class 2) tends to amplify period changes relative to single bursts (LC class 1) because overlapping bursts (i.e., increased energy deposition through fragmentation episodes) or superposition of higher-frequency peaks reinforce low-frequency signal components.

Considerable scatter is observed in both the original data sets (often nuclear-test signals) and the subsequent bolide analyses (e.g., E. Mas-Sanz et al. 2020; E. A. Silber 2024a). Such scatter might arise for reasons spanning signals emanating from different parts of the trajectory (E. A. Silber et al. 2009; E. A. Silber 2024a), altitude effects (W. N. Edwards et al. 2006), source-to-station range effects (N. Gi & P. Brown 2017), station geometry (C. Pilger et al. 2015), station noise (J. R. Bowman et al. 2005), wind-driven Doppler shifts (T. A. Ens et al. 2012), and measurement uncertainties (P. Golden & P. Negraru 2011). Though infrasound signal period, compared to signal amplitude, is less sensitive to range, it might still exhibit variations because stations may receive signals from different portions of the bolide path (E. Silber & D. C. Bowman 2025; E. A. Silber 2025). Furthermore, the range of shock configurations (from straightforward cylindrical line sources to localized quasi-spherical explosions) might plausibly explain the significant scatter in simple period–energy relations. Thus, determining the prevailing shock geometry (cylindrical versus quasi-spherical) is critical for accurately linking measured infrasound periods to underlying energy-deposition physics. We further discuss these points in the next sections.

4.3. Light-curve Classification and the Fragmentation Process

Despite their practical utility in distinguishing single-burst (LC class 1) from multiburst (LC class 2) bolides, CNEOS light-curve records have inherent limitations. These records often capture only the brightest segments of a bolide’s luminous path, often omitting earlier or later phases of ablation. Moreover, the absence of satellite sensor specifications (e.g., spectral response, dynamic range, or field of view) and undisclosed platform positioning precludes robust correction for parallax or viewing-angle effects. Although the reported velocities and altitudes may carry inherent uncertainties (H. A. R. Devillepoix et al. 2019; M. Hajduková et al. 2024), these do not impact the present study. Nonetheless, the light-curve data near peak brightness remains invaluable for determining whether a bolide undergoes a single dominant burst or multiple discrete fragmentation episodes. Even with partial coverage and uncertain sensor calibration, the light-curve signal provides essential information regarding the timing and energy release pattern, allowing classification into LC class 1 or LC class 2. Furthermore, when combined with a velocity vector (if available), the data allow for reconstruction of the approximate altitude span of peak brightness and fragmentation (E. A. Silber & V. Sawal 2025), an invaluable constraint when evaluating both the shock-formation details and the potential for multiepisode bursts.

For LC class 1 events, most of the bolide’s energy is released abruptly as a single disruption, typically transitioning from a near-cylindrical shock to a quasi-spherical expansion. Consequently, this configuration often yields moderate period–energy slopes (Figure 12(d)) because the limited fragmentation constrains shock complexity. In contrast, LC class 2 events tend to yield steeper slopes because they exhibit multiple bursts along the bolide’s path, with each releasing a significant fraction of the total energy. From a shock-physics perspective, these overlapping bursts effectively stack quasi-spherical expansions along the flight, broadening the low-frequency infrasound signature. Consequently, when bolide light-curve data are available, the relationships presented in Table 2 offer more reliable yield estimates.

4.4. Effects of Entry Geometry and Fragmentation on Period–Energy Relations

Bolide entry geometry strongly influences shock development and propagation, thereby shaping the log(E)–log(P) relationship from our regression analyses. Fundamentally, a bolide traveling steeply at low altitudes concentrates its energy in a narrow near-cylindrical or quasi-spherical shock region, resulting in sharper period–energy gradients. In contrast, shallow or high-altitude trajectories introduce more complex shock interactions and extended propagation paths that weaken the energy–period correlation. This observation that steeper angles and lower altitudes correspond to higher slopes supports the premise that geometry governs the conversion of kinetic energy into infrasound. When a meteoroid enters at a steep angle (θ > 60°) and penetrates to relatively low altitudes, it deposits substantial energy into a narrow vertical column (see Figure 2). The shock may initially resemble a cylindrical line-source and then transition into a quasi-spherical expansion with abrupt fragmentation, yielding a larger blast radius and stronger coupling between energy and period. Our regression fits under these conditions indicate that moderate energy increases correspond to notably longer acoustic periods. In contrast, shallow-angle entries (θ < 30°) typically disperse energy over a longer horizontal path, with multiple trajectory segments generating overlapping shock fronts, especially when multifragmentation occurs. This elongated geometry leads to partial wave interference and more gradual “smearing” of higher-frequency shock components, thus reducing the observed slope in period–energy space. Altitude further modulates these effects: deeper bursts (<30 km altitude) frequently behave like localized explosions, whereas higher-altitude fragmentation often allows greater shock attenuation.

The distinction between bolides and nuclear explosions is apparent here. Assuming a single point-source explosion for bolides may be misleading in high-energy, shallow-entry cases such as the Chelyabinsk superbolide (O. P. Popova et al. 2013; C. Pilger et al. 2015; E. A. Silber 2024a, 2025). Traditional infrasound energy–period relations based on point-source models do not capture the complex shock geometries resulting from multifragmentation or extended atmospheric paths. Recent work by E. A. Silber (2025) demonstrates that shallow angles can result in significant station-to-station variability in observed back azimuths (exceeding 10° over thousands of kilometers) because acoustic energy emerges from multiple points along the bolide’s trail rather than from a single point. This multisegment emission can lead to substantial discrepancies between observed and modeled back azimuths when a unique peak brightness or point-source location is assumed (E. A. Silber 2024a). Under favorable energetic and geometric conditions, anisotropy in acoustic emissions may occur, as observed in the Chelyabinsk event (C. Pilger et al. 2015).

4.5. Influence of Bolide Physical Parameters

Dense, cohesive bolides and those with high dynamic strength are negatively correlated with altitude (Figure 12), suggesting that they retain structural integrity until lower altitudes. In these conditions, energy is released into a denser atmospheric region, creating a more localized shock that intensifies wave amplitude and period. Velocity shows only moderate correlations with energy or size, suggesting that a wide range of orbital parameters can obscure simple mass or composition effects. Most meteoroid physical properties (e.g., strength, shape, density, and size) control the ablation and fragmentation rates, and therefore the energy-deposition profile that generates the detectable infrasound signal. Furthermore, the expansion characteristics of the ablated vapor, particularly from the noncondensable volatile component, influence the shock radius and period. This is analogous to the Lamb wave coupling efficiency proposed by M. Boslough & V. Titov (2024), which is suggested to be greater for explosive volcanic events associated with volatile-rich magmas, and impact cratering events into volatile-rich target rocks. Overall, our results indicate that no single physical property dominantly governs the infrasound period. Nonlinear processes during hypervelocity entry (e.g., orientation-dependent drag, the response of internal structure to stress, and sequential fragmentation dynamics) contribute to the complex interplay that determines the shock wave characteristics.

4.6. Dependency on Range

One evident pattern from both the bootstrap method and OLS is that shorter ranges (generally ≤3000 km) often yield steeper slopes (A) compared to the global average (Figure 11), confirming that this effect is robust rather than resulting from a few outliers. In principle, more energetic bolides produce lower frequency infrasound components. As a result, these more powerful events can be detected at longer ranges (Figure 13) because low-frequency energy propagates efficiently through the atmosphere without substantial dissipative losses (L. B. Evans et al. 1972; L. C. Sutherland & H. E. Bass 2004). This selectivity means that many far-range detections in large-area data sets often come from bolides with enough energetic output to maintain a discernible infrasound signal over thousands of kilometers. In far-field propagation, the initially complex shock can become effectively smoothed into its dominant low-frequency component. Consequently, once predominantly low-frequency content is all that remains, the measured periods do not vary as sharply with changes in bolide energy as they might at closer ranges. Thus, the observed period at long ranges can exhibit a slightly better correlation with energy than closer stations, where a fuller frequency spectrum is measured. While this ensures detectability for major bolides, it can also mask or average out certain mid- to high-frequency indicators of the event’s energy, such as details linked to multiple bursts or rapid fragmentation. Our detection-range relation Equation (6) diverges somewhat from that of T. A. Ens et al. (2012), likely because our larger, more diverse data set includes additional low-energy bolides not captured in earlier work, resulting in longer average detection ranges for modest-energy events.

4.7. Are Period-based Energy Relations Reliable?

A key outcome from our analysis is that no single global relation fully encapsulates the entire range of bolide infrasound behavior. While a single global line provides a convenient baseline, caution is warranted for nearby, localized events, where shock coupling can be stronger or more variable. Moreover, our findings show that while the classical slopes and intercepts from global fits remain practically useful, a parameter-aware approach is warranted when angles are shallow, altitudes are low, fragmentation is multiepisodic, or station ranges are especially short or long. Single-value relations are effective for initial energy estimates but may not capture extremes in geometric configuration or shock dynamics. Refining these relations through subgroup-specific fits (Tables 2 and 3) can improve the accuracy of infrasound-based energy assessment for both planetary defense applications and scientific investigations. Even if little else is known about a bolide beyond its infrasound record, our results show that period measurements alone can still yield a reasonable first approximation of total energy, and, if multistation records are available, geolocation. Whenever a station network detects sufficiently strong infrasound signals, regardless of uncertainties about flight path or fragmentation mode, one can employ established period–energy formulas to gauge yield. The effective propagation of low-frequency infrasound over large distances renders it a primary method for rapid bolide characterization, especially when optical or radar observations are lacking. Furthermore, refining energy estimates through data partitioning (e.g., by distance or inferred multifragmentation) advances risk assessment when supplementary information on the object’s physical attributes are unavailable.

4.8. Broader Implications for Planetary Defense

Our results have direct implications for planetary defense, where accurate bolide energy estimates are crucial for assessing potential risks from near-Earth asteroids (NEAs). Robust correlations among diameter, density, and dynamic strength with total energy reveal that larger, more cohesive bolides tend to penetrate deeper, resulting in lower-altitude breakups that increase the potential for surface damage and have distinct infrasound signatures. These observations support the use of period measurements as a reliable, rapidly obtained parameter for energy estimates, even though secondary factors (e.g., entry angle and altitude) modulate the shock’s behavior. Clarifying how multifragmentation (LC class 2) intensifies the period–energy link through overlapping quasi-spherical expansions reinforces the need to account for fragmentation in yield predictions, particularly for low-altitude, multiburst events. Recognizing that no single factor dominates the infrasound signal emphasizes the importance of multiparameter models that incorporate density, dynamic strength, angle, and altitude. The correlations derived here provide a framework for planetary defense agencies to assess potential hazards from incoming NEAs by linking fundamental shock physics to practical infrasound diagnostics.

Future work may benefit from advanced models, including Bayesian hierarchical frameworks and high-fidelity wave-propagation codes, to address complexities such as multiepisode fragmentation and shallow-angle entries. Nonetheless, our bootstrap and partition methodology form a practical basis that bridges classical period–energy formulas with the complexity of real bolide observations, meeting the needs of both planetary defense and meteor physics research.

5. Conclusions

We have conducted a comprehensive investigation into empirical period–energy relations for bolides. We compiled consolidated tables of previously published period- and amplitude-based energy relations, correcting minor unit inconsistencies or typographical issues that have persisted in the literature. By incorporating an expanded bolide database that includes USG light curve and atmospheric Doppler-wind profile data, and by applying both classical OLS fits and a bootstrap framework that partitions data by various parameters (e.g., altitude, source-to-station distance, entry angle, bolide mass, and light curve classification), we have moved beyond the limitations of a single, global relation. This methodology provides a more comprehensive, physics-consistent means of refining yield estimates than relying on any single, universal formula.

Our refined period–energy fits demonstrate that while global or averaged fits are acceptable baselines for distant or minimally constrained events, significant variability appears in cases of low altitude, close stations, or shallow trajectories. Under these conditions, the measured infrasound period becomes much more sensitive to energy, explaining why a single universal relation may under- or over-estimate yields for specific geometries. By systematically isolating factors such as peak-brightness altitude, fragmentation type, and station distance, our study reveals the origins of deviations from classical period–energy lines. This partitioned, bootstrap-driven approach, rarely applied so extensively in bolide research, offers a transparent, data-driven method for quantifying uncertainties and testing physical assumptions in shock formation.

From a planetary defense hazard-assessment perspective, our framework allows for adaptable yield estimations tailored to bolide geometry and fragmentation type. Although an averaged global fit remains suitable for first-order estimates, subgroup partitions become crucial for shallow entries, low altitudes, or multiple fragmentation episodes where period–energy scaling significantly deviates from a single baseline. Our findings confirm that period measurements alone can provide a practical first approximation of bolide energy, even without additional observational constraints, and that multistation records enable geolocation. The demonstrated effectiveness of infrasound detection in characterizing large, potentially hazardous bolides across vast distances emphasizes its immense value for planetary defense.

Acknowledgments

Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC (NTESS), a wholly owned subsidiary of Honeywell International Inc., for the US Department of Energy’s National Nuclear Security Administration (DOE/NNSA) under contract DE-NA0003525. This article has been authored by an employee of National Technology & Engineering Solutions of Sandia, LLC under Contract No. DE-NA0003525 with the US Department of Energy (DOE). The employee owns all right, title, and interest in and to the article and is solely responsible for its contents. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this article or allow others to do so, for United States Government purposes. The DOE will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan https://www.energy.gov/downloads/doe-public-access-plan. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the US Department of Energy or the United States Government. Additionally, this article has been authored by an employee or employees of Triad National Security, LLC (Triad) under contract with the US Department of Energy (DOE). Accordingly, the US Government retains an irrevocable, nonexclusive, royalty-free license to publish, translate, reproduce, use, or dispose of the published form of the work and to authorize others to do the same for US Government purposes.

Data Availability

CNEOS fireballs database is available at: https://cneos.jpl.nasa.gov/fireballs/. Light curve data are available on the CNEOS webpage: https://cneos.jpl.nasa.gov/fireballs/lc/. Data tables for all events examined in this study, along with the accompanying bootstrap outputs, can be accessed through Harvard Dataverse, doi:10.7910/DVN/P6XZ67.

Funding

This work, in part, was supported by the Laboratory Directed Research and Development (LDRD) program (project number 229346) at Sandia National Laboratories, a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA0003525. J.M.T.-R. acknowledges financial support from the Spanish research project PID2021-128062NB-I00 funded by MCIN/AEI/10.13039/501100011033. E.P.-A. acknowledges support from the LUMIO project funded by the Agenzia Spaziale Italiana (2024-6-HH.0). I.O. acknowledges the Universities Research Association-Sandia Graduate Student Fellowship Program which supported this research. S.C. acknowledges support provided by NSF Award EAR-1852339 and the NSF SAGE and GAGE facilities, operated by EarthScope Consortium. P.J. is supported by NASA award 80NSSC22K1467 and a grant from the NASA Asteroid Threat Assessment Project.

Declaration of Competing Interests

The authors declare no competing interests.

Appendix A

Appendix. A.1. Empirical Amplitude-range Yield Relations

In this section, we provide a table with empirical amplitude-range yield relations (Table A1). These relations are in the following form:

Equation or symbol description not available

where E is energy in kilotons (or tons) of TNT equivalent, P is infrasound signal period in seconds, and A, B, and C are the regression coefficients, ${A}^{* }$ is the maximum amplitude (${A}_{\max }$) or wind-corrected maximum amplitude (${A}_{w}$) in Pa, and R* is the range from the source to receiver in either km (R) or degrees (Δ). The amplitude correction, derived by R. W. Whitaker (1995), is:

Equation or symbol description not available

Table A1. List of Amplitude-range Relations

ABCA* [Pa]R*Source
1.000.50−1.54 ${A}_{\max }$ R sin ΔA. D. Pierce & J. W. Posey (1971)
2.002.94−1.84 ${A}_{\max }$ ΔD. Clauter & R. Blandford (1998)
1.472.00−4.96 ${A}_{\max }$ RR. W. Whitaker (1995)
3.033.03−9.09 ${A}_{\max }$ RRussian-crosswind, J. Stevens et al. (2002)
3.033.03−10.00 ${A}_{\max }$ RRussian-downwind, J. Stevens et al. (2002)
2.003.52−10.62 ${A}_{\max }$ RE. Blanc et al. (1997)
1.492.00−4.18 ${A}_{w}$ RJ. P. Mutschlecner & R. W. Whitaker (2009)
1.552.00−8.45 ${A}_{w}$ RM. Davidson & R. W. Whitaker (1992)

Note. Range (R*) is listed in km or degrees, while A* refers to either maximum amplitude (Amax) or wind-corrected amplitude (${A}_{w}$) in Pa. J. Stevens et al. (2002) included all arrivals, including thermospheric.

Download table as:  ASCIITypeset image

Appendix. A.2 Amplitude-range (Wind Corrected) Relations

There is a family of amplitude-range relations that incorporate an empirical correction constant ($k$) accounting for the average wind velocity (${v}_{W}$) (in m s−1) along the great-circle path from source to receiver (Tables A2 and A3). In these expressions, the signal amplitude (${A}^{* }$) is given either as maximum amplitude (${A}_{\max }$) or peak-to-peak amplitude (${A}_{P2P}$) in Pa:

Equation or symbol description not available

Table A2. List of Amplitude-range Relations that use Wind Correction

TypeA* [Pa]abck [s/m]UnitsSource
<3.5 kt ${A}_{\max }$ 4.06−0.99−0.43−0.0084tT. A. Ens et al. (2012)
<3.5 kt ${A}_{P2P}$ 4.33−1.00−0.43−0.0084tT. A. Ens et al. (2012)
>7 kt ${A}_{\max }$ 10.75−1.150.780.0014tT. A. Ens et al. (2012)
>7 kt ${A}_{P2P}$ 11.35−1.200.770.0013tT. A. Ens et al. (2012)
all E, all R ${A}_{\max }$ 4.38−1.06−0.47−0.0068tT. A. Ens et al. (2012)
all E, all R ${A}_{P2P}$ 4.71−1.08−0.46−0.0068tT. A. Ens et al. (2012)
<3000 km ${A}_{\max }$ 4.97−1.18−0.45−0.0086tT. A. Ens et al. (2012)
<3000 km ${A}_{P2P}$ 5.30−1.2−0.45−0.0084tT. A. Ens et al. (2012)
<3.5 kt ${A}_{\max }$ 3.21−1.75−0.0174ktW. N. Edwards (2007)
<3.5 kt ${A}_{P2P}$ 3.36−1.74−0.0177ktW. N. Edwards (2007)
>7 kt ${A}_{\max }$ 2.18−1.26−0.0024ktW. N. Edwards (2007)
>7 kt ${A}_{P2P}$ 2.58−1.35−0.0018ktW. N. Edwards (2007)
all ${A}_{\max }$ −0.46ktJ. P. Mutschlecner et al. (1999)
all ${A}_{\max }$ −1.28−0.43ktJ. W. Reed (1972)

Note. Range (R) is listed in km, while A* refers to either maximum amplitude (${A}_{\max }$) or peak-to-peak amplitude (${A}_{P2P}$) in Pa. T. A. Ens et al. (2012) and W. N. Edwards (2007) derived their expressions from bolide observations. J. P. Mutschlecner et al. (1999) worked with nuclear-test data and J. W. Reed (1972) used high-yield explosive tests. Some relations are expressed in tons (t), whereas others are in kilotons (kt), as indicated in the “Units” column.

Download table as:  ASCIITypeset image

Table A3. List of Half-Yield Amplitude-range Relations that use Wind Correction

TypeAmplitude [Pa]abK [s/m]UnitsSource
<3.5 kt ${A}_{\max }$ 3.55−0.92−0.0083tT. A. Ens et al. (2012)
<3.5 kt ${A}_{P2P}$ 3.79−0.93−0.0082tT. A. Ens et al. (2012)
>7 kt ${A}_{\max }$ 5.79−1.39−0.0006tT. A. Ens et al. (2012)
>7 kt ${A}_{P2P}$ 6.24−1.44−0.0008tT. A. Ens et al. (2012)
all E, all R ${A}_{\max }$ 4.04−1.01−0.0068tT. A. Ens et al. (2012)
all E, all R ${A}_{P2P}$ 4.28−1.02−0.0068tT. A. Ens et al. (2012)
<3000 km ${A}_{\max }$ 4.44−1.10−0.0086tT. A. Ens et al. (2012)
<3000 km ${A}_{P2P}$ 4.68−1.11−0.0084tT. A. Ens et al. (2012)
<3.5 kt ${A}_{\max }$ 3.21−1.75−0.0174ktW. N. Edwards (2007)
<3.5 kt ${A}_{P2P}$ 3.36−1.74−0.0177ktW. N. Edwards (2007)
>7 kt ${A}_{\max }$ 2.18−1.26−0.0024ktW. N. Edwards (2007)
>7 kt ${A}_{P2P}$ 2.58−1.35−0.0018ktW. N. Edwards (2007)

Note. Range (R) is listed in km, while amplitude refers to either maximum amplitude (${A}_{\max }$) or peak-to-peak amplitude (${A}_{P2P}$) in Pa. Some relations are expressed in tons (t), whereas others are in kilotons (kt), as indicated in the “Units” column.

Download table as:  ASCIITypeset image

Here, a, b, and c are the regression coefficients, range (R) is in km, and the constant k is expressed in units of (s/m). Half-yield amplitude-range relations are of the following form, where the coefficient c is replaced by 0.5:

Equation or symbol description not available

Appendix B

The infrasound signal periods used in this study were compiled from published data sources (E. A. Silber et al. 2011; T. A. Ens et al. 2012; N. Gi & P. Brown 2017; T. Ott et al. 2019; C. Pilger et al. 2020; T. Ott et al. 2021; P. Hupe et al. 2024; E. Silber et al. 2024). The detailed procedures for measuring infrasound periods are described in T. A. Ens et al. (2012) and E. A. Silber (2024a). In essence, these methods identify the main segment of the waveform corresponding to the peak amplitude, and then measure successive zero crossings (ZCs) to derive one or more fundamental periods. The dominant period is then computed from these zero-crossing intervals, often cross validated by a spectral analysis to confirm the peak frequency (see E. A. Silber 2024a). An illustrative example of this measurement process is shown in Figure B1, demonstrating how the maximum and minimum amplitudes are located, followed by consecutive ZCs that yield the signal period. Readers interested in a comprehensive discussion of the signal processing steps, potential sources of uncertainty, and station-to-station variability are referred to the above references. Because the measurements used in this work were taken directly from the literature, we do not replicate the entire procedure here.

Figure B1. Refer to the following caption and surrounding text.

Figure B1. Schematic representation of the procedure for measuring infrasound signal parameters and quantifying the signal period. Each ZC is marked where the waveform crosses the zero-pressure line. In this example, ZCs labeled as ZC1, ZC2, ZC3, and ZC4 define the segments used to calculate the period. Notably, based on the analyst’s judgment and the appearance of the waveform, ZC1 may be selected either before the maximum amplitude or at its onset. The green and gray sets of ZCs illustrate these possible scenarios. This flexibility in choosing ZC1 ensures that the extracted period accurately reflects the characteristics of the shock signal and accommodates variations in waveform morphology. Full-cycle ZCs define discrete periods. Averaging these discrete periods yields the mean signal period, and the standard deviation provides an uncertainty estimate. The maximum and minimum amplitudes (±A) are labeled, and the peak-to-peak amplitude is computed as (+A)+∣−A∣. This method enables consistency when extracting signal parameters such as amplitude, period, and their statistical variability.

Standard image High-resolution image

Please wait… references are loading.
10.3847/1538-3881/add47d