Modeling the Morphology of Fast Radio Bursts and Radio Pulsars with fitburst

We present a framework for modeling astrophysical pulses from radio pulsars and fast radio bursts (FRBs). This framework, called fitburst, generates synthetic representations of dynamic spectra that are functions of several physical and heuristic parameters; the heuristic parameters can nonetheless accommodate a vast range of distributions in spectral energy. fitburst is designed to optimize the modeling of features induced by effects that are intrinsic and extrinsic to the emission mechanism, including the magnitude and frequency dependence of pulse dispersion and scatter broadening. fitburst removes intrachannel smearing through two-dimensional upsampling, and can account for phase-wrapping of “folded” signals that are typically acquired during pulsar-timing observations. We demonstrate the effectiveness of fitburst in modeling data containing pulsars and FRBs observed with the Canadian Hydrogen Intensity Mapping Experiment telescope.


INTRODUCTION
One of the key observables in radio-transient astrophysics is the underlying shape of the detected spectral energy over time and electromagnetic frequency, which we hereafter refer to as "morphology."In the context of radio pulsars and fast radio bursts (FRBs), aspects of the observed pulse morphology often serve as tools for probing various physical phenomena that are both intrinsic or extrinsic to the underlying emission mechanism.For example, the morphology of radio pulsars usually exhibits a high degree of constancy that forms a central tenet to pulsar timing; this stability allows for accurate and precise measurements of arrival times, and for the resolution of time-domain effects with broad applications (e.g., Lorimer & Kramer 2012).By contrast, the emerging field of FRBs is providing strong evidence for delineations between distinct types of FRB sources, based purely on morphology, that may reflect underlying differences in progenitor populations and/or emission physics (e.g., Pleunis et al. 2021a).
The primary quantity that characterizes single-pulse morphology is the dynamic spectrum, D(ν, t), which maps source intensity across different values of time (t) and electromagnetic frequency (ν).For radio pulsars, the periodic nature of the emission leads to time-averaged morphology being defined across ν and rotational phase ϕ ≡ ϕ(t ′ ) where t ′ is the timescale of rotation, meaning that the "phase-resolved dynamic spectrum" D ≡ D(ν, ϕ) for pulsars. 1 In either case, D encodes information that is intrinsic and extrinsic to the radiation source, arising from the generation and propagation of radio emission.
Prior studies have characterized the effects that impact observed morphology using frameworks which model a subset of physical effects while marginalizing over other data dimensions and/or avoiding treatment of the underlying spectral energy distribution (SED), both of which can impact parameter estimation (e.g., Ravi 2019;Aggarwal et al. 2021;Jankowski 2022).Moreover, the rise of broadband radio telescope receivers is producing wide-band data that allow for simultaneous evaluation of SEDs and parameters that govern propagation effects.In this work, we describe a framework that generates broadband models of D in terms of physical and phenomenological parameters that encapsulate all aspects of morphology and for bursts of different SED shapes.Our framework, hereafter referred to as fitburst 2 , is a substantial expansion of the initial spectrum model used to study FRB 20110523A (Masui et al. 2015).
The goal of fitburst is to generate synthetic, noiseless representations of D that are functions of several macroscopic quantities which can then be studied through subsequent statistical analyses.We consider this work complementary to recent efforts in constructing broadband "portrait" models of D for pulsar-timing data (e.g., Pennucci et al. 2014;Liu et al. 2014;Lentati et al. 2017;Pennucci 2019), as our methods seek to characterize morphological structures while others assume prior knowledge of the underlying pulse shape.The fitburst algorithm has been under development over the past few years and initial descriptions of the algorithm and results obtained with it have been described across previous works (e.g., Masui et al. 2015;CHIME/FRB Collaboration et al. 2020;Pleunis et al. 2021a;Andersen et al. 2023;Sand et al. 2023).
In this work, we describe fitburst in its fullest generality, publicly release the codebase and, for the first time, present a detailed description of its exact estimation of the covariance matrix.In Section 2, we describe the theoretical basis of our models for D. Afterwards, we outline the optimization procedure for smoothly varying and stochastic SEDs in Section 3. In Section 4, we demonstrate the effectiveness of fitburst when applied to different types of data acquired with the Canadian Hydrogen Intensity Mapping Experiment (CHIME; CHIME Collaboration et al. 2022) telescope.We conclude our work with discussions regarding various technical aspects of fitburst in Section 5.

MODELS OF DYNAMIC SPECTRA
In order to construct a morphological model, we assumed that any pulse within a dynamic spectrum can be intrinsically characterized by its temporal width (σ) and mean time of arrival at the detector (t 0 ); the underlying temporal shape can then be chosen to reflect Gaussian or other forms, and model comparisons can be enacted to discern which of these shapes best describe the data.In what follows, we assumed that any distinct pulse present in D is intrinsically Gaussian in its temporal shape.The frequency dependence of a pulse observed in D(ν, t), however, is subject to greater modeling uncertainty due to observed differences in SEDs between pulsars and FRBs.
At an extrinsic level, pulse dispersion and broadening from multi-path scattering will induce frequency-dependent variations in observed arrival times and morphology of all visible pulse structure (e.g., Cordes 2002).These respective plasma effects are often characterized by a "dispersion measure" (DM) and a timescale of one-sided, exponential scattering-broadening of an intrinsically Gaussian profile (τ ), though deviations from ideal scattering geometries can lead to complex pulse morphology (e.g., Geyer & Karastergiou 2016).The "thin-screen" formalism of scatterbroadening, widely used in pulsar and FRB astronomy, predicts that τ ∝ ν −4 , while propagation through cold plasma leads to a time delay ∆t ∝ ν −2 .
Using these assumptions, we defined the model of a single pulse to consist of three distinct parts: a frequencydependent SED (F k ); a dispersed temporal shape of the pulse (T kn ); and a global amplitude (A).The k and n indices denote labels for discrete frequency channels (ν k ) and time bins (t n ), respectively.For a spectrum with two or more distinct pulses, the full model (M kn ) is a superposition of models for each component labeled with an index l.The analytic expression for T kn,l depends on whether scatter-broadening is measurable (e.g., McKinnon 2014); if scatterbroadening is significant, then the full model of a dynamic spectrum with N scatter-broadened, intrinsically Gaussian pulses is written as , where 1 This distinction is important because pulsar literature often uses the phrase "dynamic spectrum" to mean D = D ′ (ν, t, ϕ)dϕ, i.e., a "true" dynamic spectrum that is marginalized over morphology.For the purposes of this present work, we define D = D ′ (ν, t, ϕ)dt in order to preserve morphology for analysis. 2https://chimefrb.github.io/fitburst/A l = 10 α l , (1) The scatter-broadening timescale τ k and dispersed timeseries t kn can be expressed in terms of a (fixed) reference frequency ν r , where τ r is the scattering timescale at ν r while {δ, ϵ} are the exponents of the frequency dependence for scatterbroadening and dispersion, respectively.The "dispersion constant" k DM = (2.41 × 10 −4 ) −1 pc −1 cm3 GHz 2 s in keeping with community convention. 3We replaced Equation 3 with a Gaussian shape of mean t 0,l and width σ l in the case where scatter-broadening cannot be resolved at the native time resolution of D kn : A summary of symbols defined and used throughout this work is presented in Table 1.

Running power-law SED for Radio Pulses
As indicated by Equation 2, we assumed the SED to follow a "running" power-law (RPL) model, where the spectral index κ l = γ l + β l ln(ν/ν r ) depends on ν. 4 Such a form can be interpreted as a first-order Taylor expansion of a spectral index κ as a function of ln(ν/ν r ), where β l = [dκ l /d(ln[ν/ν r ])]| 0 .In practice, this SED model is desirable because it can adequately describe both broadband and Gaussian-like SEDs as described in later sections.
The distinction between broadband and Gaussian-like SEDs lies in the values of {β, γ}: broadband signals correspond to β ≈ 0 and |γ| ≤ 10, while Gaussian-like SEDs corresponds to large magnitudes for both parameters.The RPL model is thus ideally flexible for analyzing a diverse range of morphology with a single framework.A discussion on the relationship between {β, γ} and the traditional Gaussian mean and variance (for band-limited signals) is provided in Appendix A.

Parameterization of Amplitude
Equation 1 shows that the "amplitude parameter" in fitburst (α) is the base-10 logarithm of the "full" amplitude (A).This choice is largely motivated by the interplay between the {α, β, γ} parameters and the chosen value of ν r , which is discussed in Appendix A. For SEDs with β ≈ 0, A will vary by factors of 1-10 for modern values of receiver bandwidth (∆ν) and values of γ that are typical for radio pulsars in the Milky Way Galaxy.However, for Gaussian-like SEDs, A will be considerably small in order of magnitude if ν r is sufficiently different from the frequency of peak emission.If A was chosen as the variable parameter, then an initial guess for A could easily differ by many orders of magnitude from its "true" value when applied to spectra that have |β| ≫ 0 and if ν r is not carefully chosen.This large difference leads to prolonged or failed searches for the optimal value of A for pulses with Gaussian-like SEDs if ν r is arbitrarily chosen.
One solution to this problem is for the user to choose a value of ν r that is approximately equal to the frequency of peak emission for a pulse with Gaussian-like SED in D kn .However, such a method will be cumbersome when seeking Temporal profile at ν k and tn for burst component l to model a large and morphologically diverse collection of D kn where pulses with Gaussian-like SEDs span different parts of an observed bandwidth.Moreover, it may be useful to obtain a large sample of {β, γ} values that are all referenced to the same ν r for statistical analysis of SED variations.By choosing α as the fit parameter, fitburst can estimate amplitudes of observed spectra in a coordinate phase space where linear steps in α correspond to logarithmic steps in A, and thus efficiently find the optimal value for A regardless of SED shape.

Global and per-pulse Parameters
Equations 1-5 define the model M kn in terms of two parameter subsets.The {DM, τ r , δ, ϵ} parameters govern the extrinsic plasma effects of dispersion and scatter-broadening; we assumed that these "global" parameters apply to all astrophysical features within a given dynamic spectrum and do not vary over the timescales of observed features.The five "per-pulse" parameters {α l , σ l , t 0,l , β l , γ l } are instead uniquely defined for each of the N pulses present in D kn , regardless of the prevalence of scatter-broadening.In all cases, ν r is considered to be a fixed parameter that can be set to any arbitrary, non-zero value.The most general form of M kn , where scatter-broadening is assumed and all variable parameters are chosen as degrees of freedom, therefore depends on N fit = (4+5N ) fit parameters.
In practice, it can be beneficial to generate models where one or more of the parameters are held fixed to predetermined values.For example, and as described above, plasma theory under commonly employed assumptions predicts that δ = −4 and ϵ = −2.Moreover, τ r = 0 ms for features where no scatter-broadening across the total bandwidth is resolvable.One or more of these fixed-parameter combinations define models that are "nested" from the most general Figure 1.Calculation of dispersion smearing in fitburst models.By default, M kn will inherit the original time and frequency resolution of D kn ; this resolution will yield a model for unresolved, highly-dispersed profiles as shown in red in the top-most timeseries.To account for smearing within each channel, fitburst will first produce a high-resolution version of each channelized timeseries, upsampled along both axes by user-specified amounts.fitburst will then evaluate the model across the upsampled grid, where intra-channel dispersion can be better resolved.(The above figure shows a dispersed model of the original timeseries in blue shading that is obtained after upsampling in frequency by a factor of 8, and in time by a factor of 4.) The highresolution model is then downsampled to the original resolution to form a dispersion-smeared pulse, where the signal of a dispersed, intrinsically low-width pulse now appears in more than one time sample.
form of M kn .Several such realizations of M kn can subsequently be compared against observed data through statistical analyses of their goodness of fit, which we discuss in later sections.We nonetheless present our model in its fullest generality, and clarify fixed-parameter assumptions when they are used during analysis.

Treatment of Intra-channel Dispersion Smearing
We used the kn-index notation above to emphasize the discretized nature of observed dynamic spectra.A consequence of discretely sampling raw telescope data is the incoherent removal of pulse dispersion when applying Equation 5 to channelized data: while all channels will be de-dispersed relative to each other, each channel will suffer a smearing of the intrinsic signal evaluated over its finite bandwidth (∆ν k ).The temporal extent of pulse smearing across a channel bandwidth ∆ν k can be estimated to be ∆t smear ≈ ϵk DM DMν ϵ−1 ∆ν k .The full, "observed" width of a pulse (∆t k ) can thus be estimated by adding the intrinsic and broadening terms in quadrature, where ∆t n denotes the limit in resolvable width from the native time resolution.Incoherent de-dispersion will always lead to a non-zero amount of ∆t smear .The combination of dispersion and time/frequency resolution therefore impacts our modeling of the intrinsic pulse width σ.
The magnitude of ∆t smear can nonetheless be minimized through an additional modeling procedure when generating M kn .Our chosen method consists of two steps: first generate M k ′ n ′ at higher frequency and time resolutions than afforded in the observed spectrum D kn ; then downsample M k ′ n ′ , in both dimensions and by appropriate factors, in order to find M kn with dimensions that match those of D kn .A schematic of this dispersion-smearing removal is illustrated in Figure 1.This method of dispersion-smearing removal introduces additional computations in our modeling, but allows for model-independent estimation of both σ and DM.
It is sometimes advantageous to manipulate raw telescope data in order to remove the deleterious effects of dispersion prior to analysis.For example, this manipulation is regularly performed for "coherently dedispersed" timing observations of radio pulsars, where intra-channel smearing for a specified DM is completely removed prior to the formation the total-intensity dynamic spectrum (Hankins & Rickett 1975).The above procedure for smearing removal can nonetheless apply for coherently dedispersed data, so long as the DM parameter instead corresponds to a value offset from the DM used for coherent dedispersion.fitburst can accommodate these differences through proper configuration when modeling data that are coherently or incoherently dedispersed.

FITTING METHODS
We evaluate the goodness of fit between M kn and D kn using a weighted-χ 2 statistic: where w k represents a per-channel "weight" that we define as the standard deviation in each channel, evaluated over a user-specified fraction of time.For our model described in Section 2, the χ 2 statistic is at most a function of the (4+5N) parameters described in Section 2.3 and therefore spans a (4+5N)-dimensional phase space at maximum.The global minimum in the χ 2 phase space represents the optimal fit of M kn to D kn , where the Jacobian vector ∇χ 2 = 0.
Our implementation of the fitburst model and fitting procedure is written in the Python language.We employed a weighted least-squares minimization algorithm in our implementation of fitburst, using Equation 8 to obtain best-fit estimates of M kn from our observed D kn .fitburst can also perform ordinary least-squares fitting, where w k = 1, though we encourage the use of Equation 8 so that arbitrary normalization of D kn can be taken into account.Our implementation of the fit procedure used the least squares routine in the scipy.optimizepackage.In order to efficiently determine the region of best fit, we derived exact expressions for the components of ∇χ 2 and supplied it as a subroutine to least squares.We also derived exact expressions for the components of the Hessian matrix, in order to estimate the covariance matrix and statistical uncertainties of all fitted parameters.For completeness, the exact form of ∇χ 2 is provided in Appendix B.

Amplitude-independent Fitting in Cases of Scintillation
The models described in Section 2 are most appropriate for signals with smoothly-varying SEDs across a given bandwidth.However, a common astrophysical occurrence is the stochastic variation of brightness across frequency due to interstellar scintillation (e.g., Rickett 1970;Cordes et al. 1986).As a consequence, the ISM imparts intensity fluctuations that cannot be adequately modeled using simple analytic functions.Direct modeling of scintillation can nonetheless be performed by estimating an amplitude for each frequency channel, though such an effort will substantially increase the number of degrees of freedom by an amount equal to the number of frequency channels (N ν ), which can often exceed N ν > 10 3 for modern telescope-receiver systems.
As shown by Taylor (1992) and Pennucci et al. (2014), scintillation can be indirectly modeled by modifying Equation 8to be defined for a subset of model parameters.We first note that the quantity A l F k,l that defines M kn represents a per-channel amplitude A k,l .In cases where scintillation is apparent, we are interested in minimizing χ 2 in the subspace where all derivatives with respect to A k,l are already minimized.This condition corresponds to evaluating the derivative of χ 2 with respect to A k,l and requiring that all such quantities be equal to zero; the result is a system of N equations that are linear functions of the A k,l parameters.In the case where N = 1, we find that which is consistent with the form derived by Pennucci et al. (2014).However, the general expression of A k,l for an N -component model is a convoluted function of D kn and T kn,l that defines a system of N linear equations for each channel with index k.This system is most easily summarized with a recursion relation: The model M kn = N l A k,l T kn,l can therefore be fully determined in an amplitude-independent manner using Equation 9for N = 1, or numerically solving Equation 10 for each value of A k,l .While ideal for generating models in the presence of scintillation, this method introduces additional modeling computations and necessarily relinquishes the power to constrain SED parameters.This circumstance leads to the reduction of per-component fit parameters to {t 0,l , σ l }, for a maximum of (4+2N ) parameters.The DM and scatter-broadening parameters can nonetheless be determined uniquely as they define T kn,l , regardless of the assumed SED model.

Assessment of Viability with Simulations
In order to assure the validity of our proposed framework, we simulated a large number of dedispersed radio pulses using the independently developed simpulse library (Merryfield et al. 2023) for analysis with fitburst.For all simulations, we used the simpulse library to produce coherently dedispersed dynamic spectra in the CHIME band (i.e., 400-800 MHz) with N ν = 1024, N = 1, σ = 1 ms, and τ r = 2 ms at ν r = 600 MHz.Moreover, all simpulse simulations assumed a value of γ = 0, corresponding to a flat spectrum; as simpulse does not use an RPL model in its calculations, we expect that all fitburst results derived from these simulations should yield β = 0, i.e., no spectral running.We then created 100 realizations with a time resolution of 40 µs after adding Gaussian noise, setting the variance such that the peak S/N ∼ 200 for each busrt.
A set of results that demonstrate the reliability of fitburst under different noise realizations is presented in Figure 2.For all fitburst fits described below, we used the weighted least-squares fitting procedure described above to optimize all parameters against the "high-S/N" data described above, except {δ = −4, ϵ = −2} which were held fixed; we applied the same initial guess {α = 0, β = 0, γ = 0, σ = 0.5 ms, τ r = 10 ms, DM = 0 pc cm −3 } for each execution of fitburst and for a fixed value of ν r = 600 MHz, leading to N fit = 7 for all fits.We did not perform any upsampling and treated the DM as an offset parameter as all spectra were constructed to be coherently dedispersed.Each panel shows the best-fit value and 68.3%-confidence statistical uncertainty in the denoted parameter as determined by fitburst.
The best-fit estimates of all parameters agree with the simpulse input values to within reasonable factors of the statistical uncertainties.In order to demonstrate this consistency, we computed "uncertainty-weighted deviates" of τ r , i.e., the deviation of each best-fit τ r measurement from the mean value of the 100-sample distribution, divided by its corresponding statistical uncertainty.The resultant distribution of uncertainty-weighted τ r deviates is shown in Figure 2; this distribution appears to reflect a normal distribution, which is expected as the noise of all simulations is derived from a normal distribution.The Anderson-Darling test statistic (A 2 ) for this high-S/N distribution of uncertainty-weighted τ r deviates is A 2 ≈ 0.8, indicating that the null hypothesis -these deviates are drawn from a normal distribution -cannot be rejected at a significance level of at least 2.5%.The uncertainty-weighted deviates for the other fitburst parameter yield similar statistics.
Moreover, this consistency between measurement and noise statistics remains valid for a separate set of 100 simulations that use the same simulation parameters as mentioned above, but possess lower S/N values by a factor of ∼ 10, i.e., S/N ∼ 20.We obtained best-fit models of these "low-S/N" bursts with fitburst to perform a similar statistical analysis as was done to the high-S/N data.The distribution of uncertainty-weighted τ r deviates for this low-S/N data set is also shown in the histogram of Figure 2. We found that the low-S/N distribution of uncertainty-weighted τ r deviates yielded an Anderson-Darling test statistic A 2 ≈ 0.28, indicating consistency with a normal distribution.The differences in A 2 between the low-S/N and high-S/N data sets is likely not significant due to the low number of samples.We reserve a fitburst analysis of bursts with even lower S/N -on the order of S/N ∼ 8-10, corresponding to survey-detection thresholds -for future work.
As a final test, we generated 100 simulations of a low-S/N, two-component burst -i.e., with N = 2 -by adding each low-S/N simpulse realization described above with a time-delayed copy of itself.The time separation between the two components was large enough to ensure that the signal from each component did not overlap with signal from the other component, but was otherwise arbitrarily chosen.We then analyzed each N = 2 simulation with fitburst in order to find optimal values of the same parameters analyzed above, leading to N fit = 12 simultaneously-fit parameters.The uncertainty-weighted deviates for τ r from the N = 2 fitburst fits are also shown in the histogram presented in Figure 2. The Anderson-Darling test statistic for this N = 2, low-S/N sample is A 2 ≈ 0.33, also indicating consistency with a normal distribution despite the substantial increase in fitburst degrees of freedom for an N = 2 model.We are therefore confident that fitburst can produce adequate fit statistics for a wide range of burst brightness and degrees of freedom needed for multi-component fits.

ANALYSIS WITH CHIME DATA
Using the framework described in Sections 2 and 3, we generated models with fitburst using a series of data products obtained with the CHIME telescope and its backends.We primarily analyzed "filterbank" (i.e., channelized timeseries) observations collected by the CHIME/FRB instrument (CHIME/FRB Collaboration et al. 2018).The raw CHIME/FRB data comprise ∼2 seconds of total-intensity timeseries, with a time resolution of 0.98304 ms, evaluated across the 400-800 MHz bandwidth of the CHIME telescope using 16,384 frequency channels.We extracted dynamic spectra for each FRB by incoherently dedispersing the filterbank using the upon-detection value of DM and t 0 determined by the real-time CHIME/FRB detection system, and preserving at least ±80 ms centered on t 0 .
Prior to fitting, all CHIME/FRB dynamic spectra are normalized, baseline-subtracted, and cleaned of radio frequency interference (RFI) by using a variance-threshold criterion appropriate for CHIME/FRB data.5While optional, the application of these or similar external algorithms to uncalibrated data, prior to analysis with fitburst, is strongly encouraged in order to minimize the contamination of D by terrestrial RFI signals.Moreover, the definition of M kn and use of Equation 8 both presume that off-pulse measurements exhibit "white" noise that fluctuate about zero flux density; therefore, baseline subtraction of D is implicitly assumed by fitburst.Nevertheless, as demonstrated in the following subsections, fitburst can adequately model a variety of signal morphology observed in (baseline-subtracted) D even in the presence of minimal RFI.The investigation of fitburst performance on D that contain "red" noise will be the subject of future work.
For all fits discussed in this work, we set ν r = 400.1953125MHz6 and held the {δ = −4, ϵ = −2} parameters fixed.All best-fit parameters are reported below with uncertainties determined from the covariance-matrix formalism discussed in Appendix B; these uncertainties are quoted in parentheses and reflect 68.3% confidence intervals.

Quality of Fits for SEDs
Figure 3 shows examples of best-fit, single-component dynamic spectra of two CHIME/FRB detections.FRB 20220502B displays a power-law SED while FRB 20200928B exhibits a band-limited, Gaussian-like energy distribution within the CHIME band.For each fit, fitburst used the upon-detection estimates of DM and t 0 as initial guesses for the least-squares fitting of the model spectrum; all other parameters were initially set to generic values of {α = 0,  for this FRB, then α = −0.367(11) and γ = −3.95(9).The case where γ < 0 conforms with the "steep-spectrum" SED that is typical of radio pulsars, but the inclusion of β as a variable parameter improves the fit by an amount ∆χ 2 ≈ 600 despite only adding one degree of freedom.The significance of β characterizes the slight turnover in power at low frequencies that deviates from a purely power-law form.
By contrast, FRB 20200928B yields a best-fit model with substantially different magnitudes of the amplitude and SED parameters: α = −24.2(1.2);β = −198(10); and γ = 207(11).No model with β = 0 held fixed can adequately describe the band-limited SED for FRB 20200928B, supporting the notion that β serves as a measure of Gaussianity as described in Appendix A. Values with these magnitudes are common for other "band-limited" signals in the CHIME/FRB data set, especially when the SED is largely encapsulated within the observed bandwidth.

Removal of Dispersion Smearing
As discussed in Section 2.4, fitburst uses an upsampling approach to account for dispersion smearing in channelized, dispersed data.The upsampling factors can be chosen arbitrarily, but should be set to values that are appropriate for the observed signal DM, receiver bandwidth, and time/frequency resolutions of the input data.
An example of impact in upsampling-factor choices for CHIME/FRB analysis is shown in Figure 4.The best-fit values and uncertainties of σ for FRB 20200722B -detected by CHIME/FRB with an observed DM of 2,157.572(4)pc cm −3 -clearly decrease for increasing amounts of frequency/time upsampling.This variation demonstrates that dispersion smearing is detectable in the CHIME band at such high DMs.The difference in the two sets of σ measurements illustrates the interplay between upsampling-factor choices and resolution of the CHIME/FRB total-intensity data: if there is no upsampling in time, then more upsampling in frequency is needed to better resolve dispersion smearing and thus estimate σ.
Recent CHIME/FRB results have used fitburst to find that intrinsic temporal widths are significantly different between confirmed repeating-FRB sources and those that have not been observed to repeat (CHIME/FRB Collaboration et al. 2019;Pleunis et al. 2021a).This difference has increased in significance as the CHIME/FRB catalog continues to grow (Fonseca et al. 2020;CHIME/FRB Collaboration et al. 2021, 2023).In practice, and unlike the no-scattering model, the scatter-broadened fitburst model is sensitive to initial-guess values for {DM, t 0,l , σ, τ r } due to their dependence in the error function.We found that fits for scatter-broadening are aided by following a two-step procedure: first by fitting a model where τ r = 0 ms is held fixed and all other parameters are optimized; and then fitting a separate, scatter-broadened model where τ r is an additional degree of freedom and the initial guess corresponds to the results of the no-scattering model.

Quality of Fits for Scatter-Broadening
The result of this two-step procedure is shown in Figure 5 for the CHIME/FRB detection of FRB 20210124B.We extracted a dynamic spectrum for FRB 20210124B from CHIME/FRB filterbank data using the DM and t 0 values determined from the no-scattering fitburst model.We then applied a N = 1 model where all parameters except {δ = −4, ϵ = −2} were simultaneously fitted, yielding a best-fit τ r = 80(4) ms at ν r = 400.1953125MHz; the rightmost panel of Figure 5 shows the best-fit difference between the model and data, hereafter referred to as "residuals", which indicates that the model optimally describes the raw data.The inclusion of τ r as a fit parameter improves the optimization of the fitburst model by an amount ∆χ 2 ≈ 700, despite only adding one degree of freedom between the no-scattering and scatter-broadened models.
Figure 5 also demonstrates the covariance between DM and τ r , and the important of simultaneously optimizing both parameters.If scatter-broadening is neglected as a model feature, then the scattered signal is absorbed by the no-scattering model as dispersion smearing of the pulse at an incorrect DM.The correct DM is recovered only when incorporating τ r as a degree of freedom; for FRB 20210124B, the difference between the optimized no-scattering and scatter-broadened estimates of DM is ∆DM = 1.16(6) pc cm −3 .

Multi-component Fitting with fitburst
FRBs often show multi-component structure in their dynamic spectra, typically consisting of distinct bursts that occur at different times and span different frequency ranges.For example, repeating-FRB sources often emit multicomponent bursts with "downward-drifting" substructure, where successive band-limited components occur at decreasing radio frequencies.We therefore designed fitburst to accommodate multi-component models in its optimization against input data.
Figure 6 shows the data, optimized fitburst model and the residuals for FRB 20211229A.For the initial-guess parameters, we used the upon-detection estimate of DM but set each arrival time to the value determined by a peakfinding algorithm applied to the band-averaged timeseries. 7We also arbitrarily set ν r = 400.1953125MHz for all components.All other per-component parameters were initially set to generic values of {α l = 0, β l = 0, γ l = 0, σ l = 0.5 ms}.The {δ = −4, ϵ = −2, τ r = 0 ms} global parameters were held fixed to their initial values; all other model parameters were subject to least-squares minimization, leading to N fit = 31.
The best-fit model for FRB 20211229A adequately describes each burst component, as indicated by the noisedominated residuals shown in Figure 6.Moreover, the fitburst framework obtained a robust model despite the large number of fit parameters and use of sub-optimal guesses for most quantities.However, the adequacy of the initial guess nonetheless relied on a peak-finding algorithm for establishing the number and arrival-time estimates of burst components.This peak-finding method assumes that the initial guess for DM can sufficiently resolve individual pulses within the band-averaged timeseries; a sufficiently incorrect value of DM will smear pulses and lead to sub-optimal determination of burst components.4.5.Fitting to "Folded" and/or Coherently-Dedispersed Data Radio pulsars are often observed using a real-time "folding" algorithm that averages consecutive pulses into a single representation of improved statistical significance, based on an existing model of neutron-star rotation (e.g., Lorimer & Kramer 2012).Moreover, modern real-time pulsar instruments usually perform "coherent dedispersion" (Hankins & Rickett 1975) in order to remove the smearing effects of pulse dispersion at the complex-voltage level.The union of these real-time processing steps form the regular observing mode of the "CHIME/Pulsar" telescope backend (CHIME/Pulsar Collaboration et al. 2021), which produces coherently dedispersed timing data for hundreds of pulsars a day, evaluated across 1,024 frequency channels and with time resolution as high as 2.56 µs.
We incorporated features into fitburst that enable direct fitting of folded and/or coherently dedispersed versions of D in pulsar data.For coherently-dedispersed data, we constructed fitburst to interpret the DM as an offset parameter measured relative to the value used for coherent dedispersion.For folded data, fitburst will first generate a realization of T kn,l over a timespan equal to twice the folding period, and then average the two periods into a single realization.An example of these features applied to CHIME/Pulsar timing data is shown in Figure 7 for PSR J2108+4516, a pulsar recently discovered to be orbiting a high-mass OBe star (Andersen et al. 2023).The two-realization step allows for wrapping of the scattered pulse and thus robust modeling of the folded spectrum.

Fitting in the Presence of Scintillation
A significant number of observed pulsars and FRBs exhibit scintillation in the 400-800-MHz band at which CHIME operates.Figure 8 shows an example of scintillation in FRB 20230702B as detected by CHIME/FRB.A statistically adequate model for D can be obtained when using the RPL SED, though clear residual structure remains.These features do not reflect the smoothly varying decay in intensity from scatter-broadening, and instead indicate incomplete modeling of the SED due to scintillation.These features were nonetheless accounted for when using the amplitudeindependent formalism presented in Section 3.1, as shown by the flat residuals in Figure 8.
The best-fit values of {t 0 , σ, DM} between the RFL-SED and amplitude-independent fitburst models for FRB 20230702B are statistically consistent, indicating that moderate degrees of scintillation do not negatively impact optimization of M kn if an RPL model of the SED is used.The amplitude-independent optimization of M kn is nonetheless visually superior and indicates that scintillation is present in D kn .The per-channel amplitudes afforded by the technique described in Section 3.1 can be further interpreted to derive scintillation parameters such as decorrelation bandwidths and timescales.We reserve such efforts to future work.

DISCUSSION & CONCLUSIONS
The CHIME data set provides an ideal opportunity for developing robust models of dynamic spectra due to its operation at low radio frequencies and apparent diversity in realizations of D. However, the framework described in Section 2 is generalized to work with data from any radio telescope.We therefore encourage interested users to apply fitburst to their raw data products, in order to obtain a census of morphological parameters for different source classes and at different parts of the electromagnetic spectrum.We conclude this work by discussing various aspects of the fitburst framework, both in terms of applications and infrastructure development.

Modularity and Runtime of fitburst
fitburst is designed to separate distinct operations that, when combined, can contribute to the modeling of D kn .This modularity within fitburst can therefore be used to simulate realizations of M kn and/or compare existing fitburst models without the need for accessing real data or executing a full fitburst pipeline.Moreover, different SED models can be incorporated into the fitburst framework, so long as exact expressions for the Jacobian-vector and Hessian-matrix components are derived and included as routines within fitburst.
The runtime of the fitburst optimization algorithm depends on the number of floating point operations needed to estimate Equation 8.This number should largely scale as N ν N t so long as the model-upsampling factors discussed in Section 2.4 are modest in comparison to N ν or N t .In order to verify this scaling, we timed the least-squares fitting  Andersen et al. (2023).fitburst can adequately account for wrapping of signal in rotational phase, which is important for pulsars with large duty cycle and/or significant scatteringbroadening as seen above.
algorithm within fitburst when executed on FRB 20219124B, the CHIME/FRB detection first discussed in Section 4.3.We found that the least-squares optimization algorithm needed ∼ 260 s to optimize all parameters except {δ, ϵ} when no downsampling is performed, i.e., when the algorithm is executed on a version of D kn with size (N ν , N t ) = (16384, 162).The runtime decreases linearly when downsampling in frequency by increasing integer multiples of 2, with the algorithm needing only ∼ 3 s when executed on a version of D kn with size (N ν , N t ) = (265, 162).These runtimes apply for fits that upsample M kn in frequency by a factor of 8 and in time by a factor of 4.

Comparison between fitburst and PulsePortraiture
The fitburst codebase complements PulsePortaiture such that the former constructs a synthetic representation of D, while the latter assumes knowledge of D (e.g., σ l ) to extract an accurate estimate of t 0 for high-precision pulsar timing.The two frameworks therefore differ in parameterization of M kn , with the variable parameters in PulsePortraiture consisting only of {t 0 , DM, and τ r }.8Moreover, PulsePortraiture executes a two-dimensional, template-matching algorithm in the Fourier domain for optimization of t 0 while fitburst currently optimizes M kn in the time domain.
Despite the differences in frameworks, estimates of {DM, τ r } from fitburst can be compared with those determined by PulsePortraiture to ensure consistency.Such a comparison was performed by Andersen et al. (2023) in their Figure 8.The CHIME/FRB detection of FRB 20230702B (left).The middle panels corresponds to residuals obtained when fitting a "full", N = 1 fitburst model as described in Section 2. The right-most panel shows residuals for the same data and model, but computed using the amplitude-independent modeling procedure described in Section 3.1.
pulsar-timing analysis of PSR J2108+4516; this pulsar exhibits persistent and extreme variations in dispersion and scatter-broadening as it orbits its OBe companion star, thus lending itself to simultaneous analysis of morphological and timing properties.The best-fit estimates of τ r between fitburst and PulsePortraiture were found to be statistically consistent.By contrast, the best-fit estimates of DM were found to be globally offset between the two frameworks by ∼0.15 pc cm −3 .Such an overall offset has been noted to arise from measurable differences in the underlying morphology estimated or assumed between the two algorithms (e.g., Pennucci et al. 2014).The variations in DM for PSR J2108+4516 were nonetheless observed to be the same between fitburst and PulsePortraiture, indicating that fitburst can reliably describe extrinsic effects while simultaneously modeling morphological properties.

Comparisons of Best-fit fitburst Models
As described above, there are various scenarios where it is desirable to derive and compare several fitburst models for the same D.An example was shown in Section 4.3, where two models were generated and compared for FRB 20210124B that fitted or ignored morphological parameters that quantify scatter-broadening.In this case, we argued that scatter-broadening was statistically significant since ∆χ 2 = χ 2 2 − χ 2 1 ≫ ∆N fit = 1 between the two models, where χ 2 2 refers to the goodness-of-fit statistic for the model that directly fits for τ r .
A rigorous comparison between nested models that differ by ∆N fit ≥ 1 can be performed using statistical tests.A commonly used scheme is the F-test, which assumes that the statistic  A simpler model comparison can be performed by comparing a best-fit fitburst model with one where all component amplitudes A l = 0, i.e., a model of zero amplitude.In this case, the value of ∆χ 2 represents a measure of model significance over the noise in D. We therefore defined a measure of signal-to-noise ratio to be S/N = ∆χ 2 when comparing best-fit and no-amplitude models.Figure 9 shows S/N measurements made for signals published in the first CHIME/FRB catalog; all S/N values from fitburst are greater or equal to those estimate from the real-time CHIME/FRB detection system, indicating that direct modeling with fitburst can extract more information than real-time detection pipelines.

Translating RPL and Gaussian Parameters
The analysis presented in Appendix A presents a framework for connecting the RPL parameters to the parameters of a corresponding Gaussian SED, the central emission frequency ν c and bandwidth σ ν .Figure 10 displays the same time-integrated spectrum of FRB 20200928B as shown in Figure 3, along with Gaussian SEDs derived from best-fit spectral and amplitude parameters that each use different fixed values of ν r .
While all RPL-based fitburst models yield comparable fit statistics, the Gaussian SEDs derived from RPL values clearly vary in accuracy: RPL models with ν r closer to the "true" value of ν c produce Gaussian SEDs that better reflect the data.This discrepancy arises due to the approximations made in deriving the relationships presented in Appendix A, based on the assumption that the band-limited signal is defined over a range of frequencies ν/ν r ∼ 1; Figure 10.Derived Gaussian SEDs from best-fit fitburst models of FRB 20200928B using the RPL SED, where each fit used a different value of νr.While all models yielded statistically consistent fits, the choice in νr ultimately impacts the accuracy in translating between RPL and Gaussian parameters for reasons described in Section 5.5.
given the low frequencies at which CHIME operates, the condition that ν/ν r ∼ 1 is only met for ν r ≈ 675 MHz, which is closest to the true value.This issue is less prominent for observations over similar bandwidths at higher radio frequencies, where the relative difference between frequency channels is smaller.
We note that results from multi-component fitting with fitburst can be used to derive statistics that describe spectral behavior seemingly unique to repeating-FRB sources.A notable example is the rate of change in the SED centroid (e.g., dν c /dt) often observed in mult-component bursts from repeating FRBs, which was first shown to exhibit a negative, linear ramp in FRB 20121102A (e.g., Hessels et al. 2019).Estimates of the drift rate can be computed from best-fit estimates of the SED parameters and their translation to estimates of {ν c , σ c } at distinct times.For low-frequency observations, we urge caution in first ensuring accuracy of the translation to {ν c , σ c } noted above and in Figure 10, prior to estimating drift rates from fitburst data.More generally, additional care must be taken to account for covariances in all fitburst parameters when estimating, e.g., dν c /dt, as incorrect dedispersion or suboptimal fits of scattering-broadening can dramatically effect constraints on the drift rate (Gopinath et al. 2024).

Use of Exact vs. Approximate Gradients
Most modern scientific-computing packages offer optimization routines that numerically estimate ∇χ 2 through finitedifference methods.While convenient, we found that using the exact form of ∇χ 2 and the Hessian matrix offers at least three important advantages: 1. when executed on CHIME/FRB data, fitburst fitting yields a runtime that is ∼200% faster than fits that use the built-in numerical methods of scipy.optimize.leastsquares; 2. for a fixed set of initial guesses, fitburst is considerably more successful in obtaining a robust model when using the exact form of ∇χ 2 ; 3. statistical uncertainties that are derived from the exact Hessian matrix can be at least ∼80% larger than those derived from the Gauss-Newton approximation of the Hessian, and thus provide conservative estimates of the underlying uncertainties in the fit parameters.
These circumstances are most consequential for efficient characterization of FRB morphology, e.g., for real-time systems like CHIME/FRB that detect several FRBs per day.As next-generation instruments work towards maximizing The inner-to-outermost red contours contain 68.3%, 95.4%, and 99.5% of all probability density, respectively.In the marginalized-PDF panels, the red lines denote the median value (solid) and 68.3% credible interval (dashed) derived directly from the PDF.The vertical gray lines denote the same measures but instead determined from computation of the exact Hessian matrix as described in Appendix B. The overlap of these statistics indicates that the exact-Hessian method yields robust measures of uncertainty and covariance.
detection rates and issue near-real-time alerts, fitburst can serve as a reliable tool for producing FRB-parameter data sets with robust statistics.
In order to confirm the robustness of our methods, we performed a Bayesian analysis of {DM, τ r } for the CHIME/FRB detection of FRB 20210124B that was discussed in Section 4.3.We first computed a two-dimensional grid of best-fit χ 2 values with fitburst for different, fixed combinations of {DM, τ r } while allowing all other fit parameters to be optimized during each fit; this χ 2 map was then converted to a two-dimensional likelihood density function p(data|{DM, τ r }) = 0.5 exp(−0.5∆χ 2 ), where ∆χ 2 = χ 2 − min(χ 2 ).We then used Bayes' theorem and assumed no prior knowledge of either parameter, meaning that the posterior distribution p({DM, τ r }|data) ∝ p(data|{DM, τ r }).The posterior probability distribution function (PDF) for either parameter can then be computed by integrating p({DM, τ r }|data) over the appropriate coordinate.
Figure 11 presents the results of this Bayesian analysis of {DM, τ r } for FRB 20210124B, including marginalized posterior PDFs for each parameter.Evaluation of the χ 2 grid for FRB 20210124B required ∼3 core-days due to the large volume of the CHIME/FRB data set.The median and 68.3% credible interval for each parameter are essentially identical to the best-fit values and uncertainties computed directly from the exact-Hessian formulation.This consistency confirms that our exact-Hessian method for computing statistical uncertainties produces reasonable estimates.

Future Prospects for fitburst
Several forthcoming works are underway that generate large-scale catalogs of morphological parameters for pulsars and FRBs observed using fitburst and CHIME data.One such study comprises the second CHIME/FRB catalog, based on total-intensity filterbank acquired over a 5-yr observing span, that presents measurements for ∼4,000 FRBs.Another effort is the first CHIME/FRB catalog of raw-voltage data acquired for ∼120 FRBs, which allow for morphological analyses down to timescales as low as 2.56 µs.A third project consists of a census of single-pulse morphology from several Galactic "rotating radio transients" observed with the CHIME/Pulsar instrument.A fourth project presents a morphological analysis of bursts from SGR 1935+2154, a Galactic magnetar recently observed to emit pulses with brightnesses comparable to nearby FRB sources (CHIME/FRB Collaboration et al. 2020;Bochenek et al. 2020), detected by CHIME/FRB over the past three years (Giri et al. 2023).
As shown above and in prior works, fitburst has been developed against CHIME data.However, fitburst is generalized to model pulses observed across any part of the radio spectrum, including the lowest frequencies where smearing effects from dispersion are most prominent (e.g., Parent et al. 2020;Pleunis et al. 2021b;Gopinath et al. 2024).fitburst is modularly designed to enable use of different underlying models of the SED and/or intrinsic shape, though such new feature will require calculations of their partial derivatives for estimation of the covariance matrix.
With this work, we have also made fitburst open-source in order to encourage community-based, algorithmic development of features that either improve existing functionality or introduce new methods for characterizing the morphology of radio-transient phenomena.One forthcoming feature will be the incorporation of a Markov Chain Monte Carlo sampling algorithm the employs the fitburst model, originally developed by Giri et al. (2023).
B9 constitute the components of the Hessian matrix, which quantifies the local curvature of the phase space about the global minimum and thus encodes information on the statistical uncertainties.
We determined exact expressions for the Jacobian and Hessian matrices using our definition of model spectrum.For a general fit parameter q i in a N -component model M kn = N l M kn,l , the partial derivatives of χ 2 needed to evaluate Equation B9 are: The computation of ∇χ 2 therefore requires the knowledge of the model's first-derivatives with respect to each fit parameter (∂M kn /∂q i ), while the Hessian matrix requires knowledge of both first-and mixed-partial derivatives.The use of per-component and global parameters, as discussed in Section 2.3, lead to slight differences in the number of terms that contribute to ∂M kn /∂q i , since global parameters span all models while per-component parameters do not.We found the first-order partial derivatives for a scatter-broadened, N -component model to be: In cases where scatter-broadening is insignificant, Equations B15-B20 are replaced with appropriate derivatives of a Gaussian temporal function since τ r and δ are not defined.
We then computed exact expressions for all possible mixed-partial derivatives ∂ 2 χ 2 /∂q j ∂q k using Equations B12-B20.An important subtlety with the Hessian calculation is that all first derivatives are functions of M kn,l , and thus any second derivative with respect to a per-component parameter requires knowledge of ∂M kn,l /∂q l .The fitburst codebase accounts for this step when computing all second derivatives for models with N > 1.Once calculated with Equation B11, the Hessian matrix is inverted to yield the covariance matrix.Figure 12 shows the best-fit correlation matrix for FRB 20220502B computed using the exact Hessian components provided by fitburst.

Figure 2 .
Figure 2. (Left.)Best-fit measurements of fitburst parameters derived from data simulated with simpulse, using an initial guess described in Section 3.2.The statistical uncertainty for each point corresponds to the 68.3% confidence interval computed by fitburst, and black-dashed lines correspond to simpulse input values.(Right.)Distributions of uncertainty-weighted deviations of τr from the mean value derived from "high-S/N" (∼ 200) and "low-S/N" (∼ 20) simulated bursts, as well as one-component (N = 1) and two-component (N = 2) bursts, each consisting of 100 simulations.All sets of simulations yield estimates of uncertainty-weighted τr deviates that largely reflect the normal distribution of the simulation noise (black line).

Figure 3 .
Figure 3.The CHIME/FRB detections of FRBs 20220502B (left) and 20200928B (right), with marginalized timeseries and SEDs shown in the outer panels for each FRB; the orange lines correspond to the best-fit M kn projected onto each axis.The best-fit models for both FRBs were optimized against raw CHIME/FRB data, though the products shown here are downsampled along the frequency axis to 256 channels.

Figure 4 .
Figure 4. Measurements of σ when using different frequency/time upsampling factors to model FRB 20200722B.Each point and uncertainty in this figure corresponds to a best-fit estimate of σ made with fitburst for a specific amount of time/frequency upsampling of input data while simultaneously optimizing all other morphological parameters, based on the method described in Section 2.4.The significant differences in σ and its statistical uncertainty indicate the presence of dispersion smearing in FRB 20200722B.

Figure 5 .
Figure 5.The CHIME/FRB detection of FRB 20210124B (left), optimized scatter-broadened model determined by fitburst (middle) and the corresponding residuals (right).The best-fit model was optimized against raw CHIME/FRB data, though each panel shows data that are downsampled along the frequency axis to 256 channels.All parameters except for {δ = −4, ϵ = −2} were simultaneously fitted using an initial guess described in Section 4.3.

Figure 6 .
Figure 6.The CHIME/FRB detection of FRB 20211229A (left) and best-fit, N = 6 model determined by fitburst.The best-fit model was optimized against raw CHIME/FRB data, though each panel shows data that are downsampled along the frequency axis to 256 channels.All parameters except for {δ = −4, ϵ = −2, τr = 0 ms} were simultaneously fitted using an initial guess described in Section 4.4.

Figure 7 .
Figure 7.The coherently dedispersed CHIME/Pulsar observation of PSR J2108+4504 on MJD 59438 and best-fit model determined by fitburst, whose results were first published byAndersen et al. (2023).fitburst can adequately account for wrapping of signal in rotational phase, which is important for pulsars with large duty cycle and/or significant scatteringbroadening as seen above.

Figure 9 .
Figure9.S/N for FRBs published in the first CHIME/FRB catalog, with S/N from best-fit fitburst models derived using the method discussed in Section 5.4.Over 95% of all fitburst models yield S/N greater than the value derived by the real-time CHIME/FRB detection pipeline.
an F -distribution for parameters {∆N dof , N dof,2 } under the null hypothesis.We included an initial set of routines within fitburst that computes the p-value for Equation 11 and uses the cumulative distribution function of the Fdistribution; when applied to the aforementioned models of FRB 20210124B in Section 4.3, we found that the p-value ∼ 10 −10 , indicating that the improvement in fit quality of the scatter-broadened model is unlikely due to chance.Such low p-values can be used as an indicator that scatter-broadened models are statistically favored over no-scattering models.We encourage the inclusion of additional tests for discerning statistically superior fitburst models.5.4.A fit-based Estimate of S/N for fitburst Models

Figure 11 .
Figure11.The two-dimensional posterior probability density function (blue shading) and marginalized PDFs (blue curves in outer panels) of {DM, τr} for FRB 20210124B.The inner-to-outermost red contours contain 68.3%, 95.4%, and 99.5% of all probability density, respectively.In the marginalized-PDF panels, the red lines denote the median value (solid) and 68.3% credible interval (dashed) derived directly from the PDF.The vertical gray lines denote the same measures but instead determined from computation of the exact Hessian matrix as described in Appendix B. The overlap of these statistics indicates that the exact-Hessian method yields robust measures of uncertainty and covariance.

k
DM DM ν ϵ k ln ν k − ν ϵ r ln ν r τ k M kn,l + 2A l F k,l

Figure 12 .
Figure 12.The correlation matrix of the best-fit model for FRB 20210124B (right) presented in Figure 5.

Table 1 .
Variables defined and used throughout this work.