Sampling

Animation of nested sampling converging on channel 3 of BATSE trigger 8099 with a FRED pulse fit.

The equation in the nested sampling animation should read

\[S(t|A,\Delta,\tau,\xi) = A \exp \left[ - \xi \left( \frac{t - \Delta}{\tau} + \frac{\tau}{t-\Delta} \right) -2 \right]\]

Photon Counting

Photon emission is a stochastic process. The mean number of photons emitted during time interval \(t\) is

\[\left< N(t) \right> = \lambda t,\]

where \(\left< N(t) \right>\) is the expectation value for \(N\) over the integration time, with Poisson rate \(\lambda\). The probability distribution for the random variable \(N\) is

\[\Pr(N(t) = n) = \frac{(\lambda t)^n e^{-\lambda t}}{n!}\]

For photons \(\vec{\gamma}=\{\gamma_1, \gamma_2, \ldots, \gamma_n\}\) emitted via a homogeneous Poisson process (\(\lambda(t)=\lambda\)), then the inter-arrival time between photon \(\gamma_i\) and \(\gamma_{i+1}\) is exponentially distributed. The inter-arrival times between photon \(\gamma_i\) and \(\gamma_{i+k}\) follow a gamma distribution with

\[X \sim \Gamma(k, \lambda^{-1}),\]

where \(k=1\) recovers the exponential distribution

\[X \sim \exp (\lambda^{-1}).\]

High energy detectors like BATSE accumulate photons at discrete multiples of their clock cycles (sampling frequencies). BATSE has a sampling frequency of 500 kHz (a clock cycle of 2 \({\mu s}\)). These detectors record each photon arrival time (time-tagged event, tte) to the nearest integer multiple of the clock cycle. For BATSE, hardware limitations restricted this to the first 32,768 Large Area Detector (LAD) events, inclusive of all detectors. Only the shortest, moderately bright \({\gamma}\)-ray bursts are contained completely within the tte data. After the tte period has been exhausted, BATSE counts are collected in 64ms intervals (Discriminator Science Data, discsc) for \({\sim 240}\) seconds after the trigger.

BATSE trigger 2611. This is a short bright \({\gamma}\)-ray burst for which the tte data cuts out before the event is over.

Notes

Pre-trigger tte data is available for all 8 LADs. Post-trigger tte data is available for the triggered detectors only.

BATSE sums the counts of the triggered detectors aboard the satellite. These are then converted to rates in units of counts per second by dividing by the integration time (bin width). To use Poisson statistics, PyGRB converts rates back to counts by multiplying by the bin widths, and forcing them to be integer. This does not take into account the detector deadtime (see below). Using rates rather than counts would underestimate the uncertainty in the radiation field, since for a Poisson distribution

\[\sigma_\text{Poisson}\sim \sqrt{N}.\]

A true Poisson process occurs at each detector. The photon arrival times are convolved with the response of the detector, resulting in a Poisson distributed count spectrum. Due to the onboard summing across triggered detectors we are forced to apply Poisson statistics to the summed counts, rather than the counts at each detector.

For \({\gamma}\)-ray bursts which are completely resolved in tte data, it is possible to analyse the counts at each detector. However, since PyGRB’s main focus is the analysis of discsc data (prebinned), which is summed over the triggered detectors, this has not yet been implemented. It is entirely possible to run the program independently over each triggered detector for tte data. There is not yet a unified joint-likelihood framework which will consider a single radiation field across these detectors.

Dead Time

BATSE has a small dead time after each photon count of approximately one clock cycle. This dead time is proportional to the energy of the incident photon (or particle event) which triggered the count [1].

\[\tau \sim \alpha \ln \frac{E_{\gamma}}{E_0}\]

Where \(E_{\gamma}\) is the energy of the incident photon, \(E_0= 5.5\) keV is the reset level of the detector, and \({\alpha}=0.75\) \({\mu s}\) is the signal decay time [2]. This means that photon counting is not a true Poisson process when the count rate approaches the sampling frequency. Rather, it follows a truncated Poisson distribution.

Poisson rate

The rate passed to the likelihood function is the sum of the individual pulses, specified in Pulse Types. An example rate for two FRED pulses would be:

\[\begin{split}\begin{split} S(t|A_1,\Delta_1,\tau_1,\xi_1,A_2,\Delta_2,\tau_2,\xi_2) = &A_1 \exp \left[ - \xi_1 \left( \frac{t - \Delta_1}{\tau_1} + \frac{\tau_1}{t-\Delta_1} \right) -2 \right] \\ + &A_2 \exp \left[ - \xi_2 \left( \frac{t - \Delta_2}{\tau_2} + \frac{\tau_2}{t-\Delta_2} \right) -2 \right] \end{split}\end{split}\]

Background

The background is by default modelled as a constant. This is sufficient for all but the longest gamma-ray bursts, excepting periods of unusually high background variability. More complex background models, such as a polynomial, can be included by specifying them as a rate function and including relevant priors.

The complete rate is then:

\[\begin{split}\begin{split} S(t|A_1,\Delta_1,\tau_1,\xi_1,A_2,\Delta_2,\tau_2,\xi_2) = B + &A_1 \exp \left[ - \xi_1 \left( \frac{t - \Delta_1}{\tau_1} + \frac{\tau_1}{t-\Delta_1} \right) -2 \right] \\ + &A_2 \exp \left[ - \xi_2 \left( \frac{t - \Delta_2}{\tau_2} + \frac{\tau_2}{t-\Delta_2} \right) -2 \right] \end{split}\end{split}\]

Likelihood

The rate function is then passed into the Poisson likelihood, which is a sum of the specified rates.

Priors

The default priors are

parameter	minimum	maximum	type	units
\(\Delta_i\)	–	–	uniform	seconds
\(\Delta_{i+1}\)	\(\Delta_i\)	–	uniform	seconds
\(B\)	\(10^{-1}\)	\(10^{3}\)	log-uniform	counts / bin
\(A\)	\(10^{0}\)	\(10^{5}\)	log-uniform	counts / bin
\(\tau\)	\(10^{-3}\)	\(10^{3}\)	log-uniform	seconds
\(\xi\)	\(10^{-3}\)	\(10^{3}\)	log-uniform	–
\(\gamma\)	\(10^{-1}\)	\(10^{1}\)	log-uniform	–
\(\nu\)	\(10^{-1}\)	\(10^{1}\)	log-uniform	–
\(\Delta_\text{res}\)	–	–	uniform	seconds
\(A_\text{res}\)	\(10^{0}\)	\(10^{3}\)	log-uniform	counts / bin
\(\tau_\text{res}\)	\(10^{-3}\)	\(10^{3}\)	log-uniform	counts / bin
\(\omega\)	\(10^{-3}\)	\(10^{3}\)	log-uniform	–
\(\varphi\)	\(-\pi\)	\(\pi\)	uniform	radians

The priors on \(\Delta_i\) and \(\Delta_\text{res}\) are determined by the length of the light-curve passed to PulseFitter.

Further reading

For more on nested sampling, the reader is referred to the following links.

https://lscsoft.docs.ligo.org/bilby/basics-of-parameter-estimation.html

https://dynesty.readthedocs.io/en/latest/overview.html

https://dynesty.readthedocs.io/en/latest/dynamic.html

References

BATSE Appendix G

1

T. Gjesteland, N. Østgaard, P. H. Connell, J. Stadsnes, and G. J. Fishman. Effects of dead time losses on terrestrial gamma ray flash measurements with the Burst and Transient Source Experiment. Journal of Geophysical Research (Space Physics), 115(3):A00E21, May 2010. doi:10.1029/2009JA014578.

2

B. W. Grefenstette, D. M. Smith, J. R. Dwyer, and G. J. Fishman. Time evolution of terrestrial gamma ray flashes. Geophysical Research Letters, 35(6):L06802, March 2008. doi:10.1029/2007GL032922.