# A VLA Search for the Geminga Pulsar: A Bayesian Limit on a Scintillating Source

## 1 Introduction

Geminga has long been known as a strong X-ray and gamma-ray continuum source (Fichtel et al. (1975)) and, in 1987, was identified with an optical star (Bignami et al. (1987)). In 1992, Geminga was detected at X-ray (Halpern & Holt (1992)) and gamma-ray (Bertsch et al. (1992)) energies as a 237-ms pulsar with a period derivative of $`10.95\times 10^{15}`$ s s<sup>-1</sup>, indicating a characteristic age of $`3.4\times 10^5`$ yr and a magnetic field of $`1.6\times 10^{12}`$ Gauss. A proper motion of $`170\pm 6`$ mas yr<sup>-1</sup> and a parallax of $`6.36\pm 1.7`$ mas have been measured for Geminga’s optical counterpart (Caraveo et al. (1996)). From these measurements, a distance of 157$`{}_{34}{}^{}{}_{}{}^{+59}`$ pc and a transverse velocity of $`140\pm 15`$ km s<sup>-1</sup> have been calculated.

Since its discovery as a gamma-ray source in 1975, many observers have unsuccessfully attempted radio detection of Geminga both as a continuum source (Spelstra & Hermsen (1984)) and as a pulsar (Seiradakis (1992)). Recently, however, three independent groups (Malofeev & Malov (1997); Kuz’min & Losovskii (1997); Shitov & Pugachev (1997)) have reported successful detections of pulsed radio emission from Geminga. All 3 detections were made at the Pushchino Radio Astronomy Observatory at a radio frequency of 102 MHz, a lower frequency than those used for previous searches. All three groups calculate a dispersion measure (DM) of roughly 3 pc cm<sup>-3</sup>. The measured fluxes vary greatly, perhaps due to interstellar scintillations, and range from 5-500 mJy. Reported pulse widths range from 10 to 120 ms. Recently, Shitov & Malofeev (1997) reported a 41 MHz detection of Geminga with mean flux density 300 mJy. Comparing this result with Malofeev and Malov’s quoted 102 MHz mean flux of 60 mJy, we may estimate a spectral index of $`\alpha 1.8`$, where $`F(\nu )\nu ^\alpha `$. Comparing Malofeev and Malov’s quoted mean flux of 60 mJy with the previous 1 mJy upper limit at 1.4 GHz (Seiradakis (1992)), we calculate a lower bound on Geminga’s spectral index of $`\alpha 1.6`$.

Since Geminga is the only known gamma-ray pulsar which is not also a strong radio source, confirming its existence as a radio pulsar and determining the reasons for its radio weakness are essential to understanding the relationship between pulsar radio and gamma-ray emission. Geminga is also important as it may be a prototype for a class of radio-weak or radio-quiet high-energy pulsars which could account for some of the 80 EGRET sources not yet identified with known pulsars or active galactic nuclei.

We have therefore revisited the search for Geminga as a radio pulsar with a long, low frequency, fast-sampled VLA observation. In this paper, we present our results, introducing Bayesian methods for determining the pulsed flux in a signal of unknown width and phase and for estimating the intrinsic flux density of a source, accounting for interstellar scintillations.

## 2 Data

We observed Geminga (PSR J0633+1746) on 2 February 1998 with the D-configured VLA in phased array mode. The data spanned a bandwidth of 25 MHz, centered on 317.5 MHz, with 14 (not necessarily consecutive) 1-MHz channels (situated to avoid known interference frequencies). Full polarization data were recorded with the High Time Resolution Processor (Moffett (1997)) using a sampling frequency of 1152 Hz. We obtained 7 one-hour observations (scans) of Geminga which spanned a total time of 9 hours. Two other known pulsars (B0329+54 and B0950+08) were observed for 10 minutes each with the same setup as that used for Geminga. Each scan was started on a 10 second tick to ensure accurate pulse phase referencing between scans. This 10-second tick was tied to the VLA’s hydrogen maser, which is compared to Universal Time through the GPS. The array was phased, with rms phase errors less than 1, between consecutive scans. We observed a strong VLA calibrator source (B0813+482) with known flux density ($`S`$ = 45 Jy) to calculate the gain in each of the 56 channels (14 frequencies and 4 polarizations). As we expect only a 7-mJy uncertainty in our flux density calibration due to radiometer noise, our main source of uncertainty is the intrinsic source flux variation. By comparing the VLA flux measurement with those from the Westerbork and Texas surveys (Rengelink et al. (1997); Douglas et al. (1996)), we estimate a maximum intrinsic flux density variation of 10%. We therefore expect errors in our final reported flux densities to be less than 10%.

## 3 Analysis

For each of the 7 Geminga data sets, we summed polarizations to calculate the total intensity and dedispersed all 14 frequency channels using a DM of 3.2 pc cm<sup>-3</sup>, the midpoint of the DM range predicted by the TC (Taylor & Cordes (1993)) model for Galactic electron density given Geminga’s measured distance and direction. Our estimated DM is consistent with those reported for the 102 MHz Geminga detections. We note than even a 1 pc cm<sup>-3</sup> error in our assumed DM would produce only 7 ms of pulse smearing, negligible for Geminga’s period of 237 ms.

Each one-hour dedispersed time series was folded using the gamma-ray ephemeris of Mattox et al. (1998) (P = $`0.23709574610`$ s, $`\dot{\mathrm{P}}=10.974012\times 10^{15}`$ s s<sup>-1</sup>, epoch = 2446600). We used a single period to fold each 1 hour data set, but updated the period between scans, using the gamma-ray ephemeris and TEMPO (Taylor & Weisberg (1989)). We note that using a single period for each scan results in only 0.2 ms of pulse smearing. Table 1 lists the modified Julian date of the midpoint of each scan, the corresponding period, and the total length of the scan. The first 6 scans are 1 hr in length, while the last is 50 minutes, yielding a total integration time of 6 hr 50 min.

Figure 1 shows the folded profiles for scans 1 through 7, in addition to the composite profile formed by calculating the phase at the start of each scan, and shifting and summing the profiles accordingly. We have verified that our phase-shifting algorithm is correct by applying it to the observatory generated 19.2 Hz waveguide switch signal. There is structure in some of the profiles of Figure 1. However, no consistent pulse shape and/or phase is apparent. Furthermore, coherently summing all profiles dramatically decreases the signal-to-noise level of any features. We have also calculated Fourier transforms of all 7 dedispersed time series. We have not detected any harmonics of Geminga’s 237-ms period.

Figure 2 shows the folded pulse profiles for our test pulsars, B0329+54 (DM = 26.8 pc cm<sup>-3</sup>) and B0950+08 (DM = 3.0 pc cm<sup>-3</sup>). These data were processed with identical dedispersion and folding routines to the Geminga data. We have subtracted the off-pulse mean to calculate $`F`$, the phase-averaged pulsed flux density. For PSR B0329+54, $`F`$ = 1071 mJy, and for PSR B0950+08, $`F`$ = 8.9 mJy. We note that we also detected both pulsars with high signal-to-noise in Fourier transforms of their time series.

Because Geminga may be too weak to be detected through its time-averaged flux, we searched for Crab-like ‘giant’ pulses, or individual pulses with amplitudes much greater than the mean pulse amplitude (Lundgren et al. (1995)). We search for these aperiodic, dispersed pulses by dedispersing each data set over a range of trial DMs and, for each dedispersed time series, recording those pulses with amplitudes above some signal-to-noise threshold. We enchanced our sensitivity to broadened pulses by repeating this thresholding with different levels of time series smoothing. We found no evidence for isolated pulses from Geminga. In Figure 3, we show the results of this analysis for Geminga, PSR B0950+08, and PSR B0329+54, from which we detect a large number of individual pulses.

### 3.1 Bayesian Pulsed Flux Estimation

In Appendix A, we derive a Bayesian procedure for estimating the pulsed flux density given a folded profile with a pulse of unknown amplitude, width, and phase. Assuming the additive noise is Gaussian, we are able to calculate a probability density function (PDF) for the phase-averaged pulsed flux density independent of pulse width, pulse amplitude, pulse phase, off-pulse mean, and off-pulse rms. Because no prior assumptions are necessary to calculate an upper limit, this method is an improvement over standard methods used to test for pulse shape uniformity and to calculate upper limits. The $`\chi ^2`$ statistic, the method most commonly used, returns only the deviation of a binned pulse profile from a uniform distribution. More sophisticated tests, such as the Z<sup>2</sup> and H tests (DeJager (1994)), developed to analyze sparse X-ray and gamma-ray profiles, are less dependent on binning and more sensitive to a wide variety of pulse shapes than the $`\chi ^2`$ statistic, but still yield only the probability that a distribution differs from uniformity. To calculate an actual upper limit, one must still assume a pulse width and phase. Furthermore, these methods do not provide a simple mechanism for calculating a PDF for the pulsed flux density. In following papers, we will use our method to calculate more accurate X-ray and gamma-ray upper limits for known radio pulsars and to search for new high-energy pulsars. The method can also be adapted to a wide range of pulse shapes, in addition to the simple square pulse presented in Appendix A.

We apply the method of Appendix A to calculate PDFs for the pulsed flux density for the folded Geminga pulse profiles shown in Figure 1. The resultant PDFs are also shown in Figure 1 and are discussed further in Section 3.2.2.

### 3.2 Interstellar Scintillations of the Geminga Pulsar

Since Geminga is a nearby (d $``$ 150 pc) pulsar, we expect its flux density to be strongly modulated by diffractive interstellar scintillations (DISS). We therefore present a method which allows us to calculate a PDF for a source’s flux density given several measurements of the scintillated flux density. We restrict our analysis to the strong scattering regime (a good assumption at low radio frequencies), where $`\varphi _{\mathrm{rms}}`$, the rms phase perturbation due to electron-density variations, is much greater than 1. In this case, the amplitude fluctuations are ‘saturated’ because the electric field has been sufficiently randomized to have complex Gaussian statistics.

We describe DISS through a characteristic timescale $`\mathrm{\Delta }t_\mathrm{d}`$ and characteristic bandwidth $`\mathrm{\Delta }\nu _\mathrm{d}`$. These quantities depend upon the scattering measure (SM), the integrated turbulence strength (assuming a power spectrum for electron density irregularities) along the line of sight to the pulsar. In the strong scattering regime, and assuming a uniform distribution of scattering material and a Kolmogorov power spectrum, $`\mathrm{\Delta }t_\mathrm{d}`$ and $`\mathrm{\Delta }\nu _\mathrm{d}`$ scale with SM, observing frequency, distance, and pulsar velocity as (Cordes & Lazio (1991); Cordes & Rickett (1998))

$`\mathrm{\Delta }t_\mathrm{d}=10^{2.52}\nu ^{1.2}\mathrm{SM}^{0.6}V^1\mathrm{s},`$ (1)   
$`\mathrm{\Delta }\nu _\mathrm{d}=10^{0.775}\nu ^{4.4}\mathrm{SM}^{1.2}D^1\mathrm{kHz}`$ (2)   
where $`\nu `$ is the radio frequency in GHz, SM is the scattering measure in kpc m<sup>-20/3</sup>, $`V`$ is the pulsar’s transverse velocity in km s<sup>-1</sup>, and D is its distance in kpc. The scaling is such that a distant source with high velocity observed at a low frequency has a shorter characteristic bandwidth and timescale than a nearby, low velocity source observed at high frequency.

Once we have characterized the DISS properties of a source (through a measurement or an estimation using Eqs. 1 and 2), we calculate the posterior PDF of the intrinsic source flux density given several measurements of the scintillated flux density using the method described in Appendix B. The source’s intrinsic flux density, $`S`$, is modulated by interstellar scintillations by a factor $`g`$, so that the measured flux density, $`F`$, is

$`F=gS+N,`$ (3)   
where $`N`$, the additive noise, may include both radiometer noise and radio frequency interference. We note that this additive noise model holds only for the case where the time-bandwidth ($`TB`$) product of the measurements (ie. integration time $`T`$ times bandwidth $`B`$) is large, $`TB1`$. In the strong scattering regime, the scintillation time-bandwidth product, $`\mathrm{\Delta }t_\mathrm{d}\mathrm{\Delta }\nu _\mathrm{d}`$, which is the characteristic size of a “scintle” in the $`\nu t`$ plane, may or may not be much larger than $`TB`$ for the measurement. When the time-bandwidth product of the measurement is smaller than that of one DISS ‘scintle,’ the PDF for $`g`$ is a one sided exponential:

$`f_g(g)=e^g,g0`$ (4)   
with CDF

$`F_g=P\{<g\}=1e^g.`$ (5)   
When $`TB\mathrm{\Delta }t_\mathrm{d}\mathrm{\Delta }\nu _\mathrm{d}`$, the measurements average over many scintles, increasing the number of degrees of freedom from 2 (for unquenched DISS) to 2$`N_{\mathrm{ISS}}`$, where

$`N_{\mathrm{ISS}}\left(1+0.2{\displaystyle \frac{B}{\mathrm{\Delta }\nu _\mathrm{d}}}\right)\left(1+0.2{\displaystyle \frac{T}{\mathrm{\Delta }t_\mathrm{d}}}\right).`$ (6)   
Then, $`g`$ is a $`\chi ^2`$ random variable with $`2N_{\mathrm{ISS}}`$ degrees of freedom whose PDF is

$`f_g(g,N_{\mathrm{ISS}})={\displaystyle \frac{(gN_{\mathrm{ISS}})^{N_{\mathrm{ISS}}}}{g\mathrm{\Gamma }\left(N_{\mathrm{ISS}}\right)}}e^{gN_{\mathrm{ISS}}}U(g),`$ (7)   
where $`U(g)`$ is the Heaviside function and $`\mathrm{\Gamma }`$ is the gamma function. For pulsars at large distances or those observed at low frequencies, the quantity $`N_{\mathrm{ISS}}`$ will approach infinity and $`f_g`$ will tend toward a delta function, $`\delta (g1)`$. Figure 4 illustrates the dependence of $`f_g`$ on $`N_{\mathrm{ISS}}`$.

#### 3.2.1 Results for Known Pulsars

For the known pulsars B0329+54 and B0950+08, $`\mathrm{\Delta }t_\mathrm{d}`$ and $`\mathrm{\Delta }\nu _\mathrm{d}`$ have been measured (Stinebring et al. (1996); Phillips & Clegg (1992)). Applying the known frequency scaling of these quantities (see Eqs. 12), we estimate $`\mathrm{\Delta }t_\mathrm{d}=142`$ s, $`\mathrm{\Delta }\nu _\mathrm{d}=25`$ kHz for B0329+54 and $`\mathrm{\Delta }t_\mathrm{d}=2757`$ s, $`\mathrm{\Delta }\nu _\mathrm{d}=102`$ MHz for B0950+08 at 317 MHz. Given these values and our observing bandwidth ($`B=25`$ MHz) and time ($`T=600`$ s), we may calculate $`N_{\mathrm{ISS}}`$ from Eq. 6. For B0329+54, $`N_{\mathrm{ISS}}=371`$, $`TB\mathrm{\Delta }\nu _\mathrm{d}\mathrm{\Delta }t_\mathrm{d}`$, and $`f_g`$ will approach a delta function. B0950+08 is much closer to Earth, so its time-bandwidth product $`\mathrm{\Delta }\nu _\mathrm{d}\mathrm{\Delta }t_\mathrm{d}>TB`$ and $`N_{\mathrm{ISS}}=1.095`$. Although this pulsar is likely in the transition regime from weak to strong scattering, $`f_g`$ will approximate an exponential and will have a similarly long tail.

Using the method of Appendix B, we calculate the PDFs for the intrinsic source fluxes of these pulsars given by our measurements. In Figure 5, we plot $`f_S(S)`$, where $`S`$ is the intrinsic phase-averaged flux density. For B0329+54, the PDF of the intrinsic source flux density is strongly peaked at the measured flux density, and follows a roughly Gaussian distribution away from the peak. For B0950+08, the PDF for the intrinsic source flux density again peaks at the measured flux, but is characterized by a broad exponential tail. We can calculate confidence intervals for the flux densities of these pulsars from their PDFs. For B0329+54, the 95% confidence interval is 968 mJy $`<S<`$ 1190 mJy, while for B0950+08, the 95% confidence interval is 5.4 mJy $`<S<`$ 880.1 mJy. We note that these pulsars’ predicted 317 MHz fluxes from the Princeton Pulsar Catalog (Taylor et al. (1993)) 400 MHz fluxes and measured spectral indices (Lorimer et al. (1995)) are 1890 mJy for PSR B0329+54 and 590 mJy for PSR B0950+08. While the predicted flux for PSR B0329+54 lies outside of our confidence interval, this discrepancy is consistent with the $`40\%`$ modulations due to refractive interstellar scintillation that have been measured for this pulsar (Stinebring et al. (1996)). Refractive ISS is probably also important in modulating the flux of PSR B0950+08, although its predicted flux does fall within our confidence interval.

#### 3.2.2 Results for Geminga

Because the scintillation bandwidth and timescale have not been measured for Geminga, we use Eqs. 1 and 2 to predict the characteristic timescale and bandwidth given our observing frequency, Geminga’s measured distance, and the SM predicted from the TC model. For Geminga, we estimate SM $`10^{4.4}`$ kpc m<sup>-20/3</sup> and, from Eq. 2, $`\mathrm{\Delta }\nu _\mathrm{d}`$ 1.5 MHz. Given Geminga’s transverse velocity of $`V_{}140`$ km s<sup>-1</sup>, Eq. 1 yields $`\mathrm{\Delta }t_\mathrm{d}`$ 275 s. Therefore, $`TB\mathrm{\Delta }\nu _\mathrm{d}\mathrm{\Delta }t_\mathrm{d}`$, $`N_{\mathrm{ISS}}16`$, and we expect statistical independence between measurements. Applying Eq. B9 to the distributions for $`f_F(F)`$ (see Figure 1), we calculate $`f_S(S)`$ for each of the 7 observations individually. These PDFs are plotted in Figure 6. Using Eqs. B8 and B9, we may use all 7 Geminga observations to calculate a composite PDF, shown in Figure 7. This PDF does peak at a non-zero value of $`S`$. However, as shown in Figure 10 (see Appendix A), the PDFs of simulated noise often peak at non-zero values. We note that the rms noise levels in the Geminga pulse profiles are roughly equivalent to those of the simulated profiles in this figure. In addition, we expect more structure in the Geminga profiles due to the large amount of radio frequency interference at 317 MHz. Furthermore, the marginalized PDFs for pulse width and phase for each of the 7 Geminga scans do not peak at a consistent value of these quantities, as expected for a real pulsar. We therefore treat our result as an upper limit and calculate a 95% confidence upper limit to Geminga’s pulsed flux of 3.0 mJy.

The results so far have assumed that the pulse width is completely unknown. We have also calculated upper limits assuming reasonable values for Geminga’s pulse width. If we assume Geminga’s radio pulse has the same width (roughly 180<sup>o</sup>) as its gamma-ray pulse, we calculate a 95% confidence upper limit to the pulsed flux density of 4.0 mJy. However, except in the case of the Crab, the radio pulse widths of the known EGRET pulsars are much narrower than the gamma-ray pulse widths. For this reason, we also calculate an upper limit for a pulse width of 25<sup>o</sup> (calculated assuming the $`wP^{1/2}`$ scaling of pulse width $`w`$ with period $`P`$ given by Biggs (1990)) of 1.6 mJy.

We note that the distribution of material along the line of sight will most likely not be uniform, and $`\mathrm{\Delta }t_\mathrm{d}`$ and $`\mathrm{\Delta }\nu _\mathrm{d}`$ may be somewhat different from those predicted by Eqs. 1 and 2. We therefore explore what happens in the limiting case where the scintillation timescale is longer than the 9-hr total timescale of our observations. For this case, we use Eq. B7 to calculate the PDF for intrinsic source flux density given all 7 observations, again integrating over the PDFs of Figure 1. The resultant PDF for this limiting case, also plotted in Figure 7, is much broader than that derived assuming statistically independent scintillations. The 95% confidence upper limit on Geminga’s pulsed flux for this case is 13.9 mJy. However, for the scintillation timescale to be as large as 9 hrs, the scattering measure would have to be 10<sup>-4</sup> times smaller than that estimated using the TC model. Alternatively, a timescale this long could also be explained by a thin screen located only 0.021 pc from the observer (Cordes & Rickett (1998)). Obviously, both of these scenarios are very improbable, and the true 95% upper limit is likely very close to 3 mJy.

## 4 Discussion and Conclusions

Comparing our 3 mJy upper limit to Geminga’s flux density at 317 MHz with the mean flux density of 60 mJy reported by Malofeev & Malov at 102 MHz, we calculate a lower limit on the spectral index of Geminga of $`\alpha 2.7`$ ($`F(\nu )\nu ^\alpha `$), comparable to the spectral index of the Crab. If the spectral index of Geminga is indeed this high, it would be extreme, as the Crab and Geminga pulsars have quite different ages, spin periods, and spin-down rates.

In Figure 8, we plot our upper limit along with previously measured upper limits and detections. Unfortunately, while there have been many radio searches for Geminga, there are few published upper limits. This is partly due to Geminga’s migrating gamma-ray positional error box. For example, Boriakoff et al. (1984) searched for radio pulsars in Geminga’s error box, determined from COS B data, but did not cover the current, well-determined position. Of special interest is the recent non-detection of Geminga at both 35 and 327 MHz by Ramachandran et al. (1998). While their quoted 327 MHz upper limit of 0.3 mJy is much lower than ours, the rms sensitivities of our 7-hr VLA observation and their 6-hr Ooty observation are very similiar. For consistency, we have plotted the upper limit obtained through applying our method to their data, taking into account scintillation and the unknown width and phase of Geminga’s pulse.

The reason for Geminga’s radio weakness is still debatable. One explanation is that some gamma-ray pulsars simply do not emit at radio energies. Halpern & Ruderman (1993) suggest that Geminga’s radio emission is quenched by the copious pair production in the inner magnetosphere. However, the transient dips in Geminga’s soft X-ray profile (Halpern & Wang (1997)) suggest that the processes responsible for supplying $`e^\pm `$ to the inner magnetosphere may be variable. Because the radio silence will be occasionally broken as this $`e^\pm `$ quenching plasma clears away, we may expect Geminga to be a transient radio emitter. If Geminga continues to be the only high-energy pulsar found to exhibit this transient phenomenon, studying its X-ray properties may be useful for determining if and why radio emission is suppressed.

A more likely explanation for Geminga’s radio weakness is simply that there is some misalignment between the gamma-ray and radio pulsar beams. Geminga’s pulse shape in gamma-rays (like that of the Crab and Vela pulsars) is broad and double-peaked, suggesting an origin in the “outer gaps” of the pulsar magnetosphere near the light cylinder. Because radio emission is expected to be associated with the open field line region centered on the magnetic poles, it is possible that our line-of-sight intersects the gamma-ray beam only. Romani & Yadigaroglu (1995), using the outer gap model for pulsar gamma-ray emission, have modeled the orientation and size of the radio and gamma-ray pulsar beams given radio and gamma-ray data on known pulsars. They find that, because the gamma-ray beam is much wider than the radio beam, 45% of young pulsars will be detected only at high energies, while only 19% of young pulsars will be detected at both radio and gamma-ray energies. Scaling from the 5 EGRET detections of radio pulsars, Romani & Yadigaroglu expect 12 nonradio pulsars to be visible in gamma-rays at flux levels comparable to the radio selected objects. This implies that most of the $``$ 20 unidentified Galactic EGRET sources are young radio-quiet or radio-weak pulsars like Geminga. This is consistent with other studies which find that the properties, such as flux variability (McLaughlin et al. (1996)) and spectral index (Merck et al. (1996)), of many unidentified sources are similar to those of the known EGRET pulsars. It has also been shown (Yadigaroglu & Romani (1997)) that many of the unidentified sources are associated with regions of massive star formation and death, which are expected to be breeding grounds for pulsars.

We thank J. Mattox and K. Xilouris for helpful discussions. M. A. McLaughlin acknowledges support from an NSF fellowship. The work was also supported by the National Astronomy and Ionosphere Center, which is operated by Cornell University under cooperative agreement with the National Science Foundation (NSF). The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This research was partially funded by NSF grants AST93-15285 and AST95-28394.     

## Appendix A Bayesian Pulsed Flux Estimation

Pulsar flux densities are usually reported as phase-averaged quantities. Therefore, to calculate a pulsed flux upper limit, one must assume a pulse width, pulse phase, off-pulse mean, and off-pulse rms. To avoid trial and error estimates, we have developed a method to calculate a PDF for the pulsed flux density independent of these quantities.

We assume an $`M`$-bin folded pulse profile, where the amplitude of each bin is described by $`d_i`$ ($`i1,M`$). Our model for $`\widehat{d_i}`$, the expected amplitude within each bin, is parametrized by   
w \- width of square pulse   
m \- first bin of pulse   
F \- phase-averaged pulsed flux   
$`\mu `$ \- off-pulse mean   
$`\sigma `$ \- off-pulse rms

Given this model, we may express $`\widehat{d_i}`$ as

$`\widehat{d_i}=\{\begin{array}{cc}FM/w+N_i\hfill & \text{}\text{}\text{m, m+w}\text{1}\hfill \\ N_i\hfill & \text{otherwise}\hfill \end{array}`$ (A3)   
where $`N_i`$ is the amplitude of the noise in bin $`i`$. The likelihood function of the data is given by the joint probability density function (PDF) of the noise samples. Assuming Gaussian noise with mean $`\mu `$ and rms $`\sigma `$, and statistical independence of the bins, we may write this as

$`=(2\pi \sigma ^2)^{M/2}\mathrm{exp}\left[{\displaystyle \frac{1}{2\sigma ^2}}{\displaystyle \underset{i=1}{\overset{M}{}}}(N_i\mu )^2\right].`$ (A4)   
To find the best model parameters, we use Bayes’ Theorem to write the posterior probability $`p(F,w,m,\mu ,\sigma |𝐝,𝐈)`$ of a model parametrized by $`F,w,m,\mu ,`$ & $`\sigma `$ given the data d and any prior information I as

$`p(F,w,m,\mu ,\sigma |𝐝,𝐈)={\displaystyle \frac{p(𝐝|F,w,m,\mu ,\sigma ,𝐈)p(F,w,m,\mu ,\sigma |𝐈)}{p(𝐝|𝐈)}}.`$ (A5)   
The factor $`p(𝐝|F,m,w,\mu ,\sigma ,𝐈)`$ is simply the likelihood function given in Eq. A4. Because we expect all parameter values to be equally likely a priori, we assume flat priors for $`F,w,m,`$ & $`\mu `$ (i.e. $`p(F|𝐈)`$ = constant) and a prior $`p(\sigma |𝐈)1/\sigma `$ (i.e. $`p(\mathrm{log}(\sigma )|𝐈)`$ = constant). We choose a prior that is uniform in $`\mathrm{log}\sigma `$ to express our ignorance of the scale of the noise variation. Then, we can express the posterior probability as simply

$`p(F,m,w,\mu ,\sigma |𝐝,𝐈){\displaystyle \frac{}{\sigma }}.`$ (A6)   
Since we would like to calculate the PDF of the phase-averaged flux $`F`$ independent of the other parameters, we must marginalize the posterior probability over the “nuisance parameters” $`w,m,\mu ,`$ & $`\sigma `$. We write the PDF for $`F`$ as

$`f_F(F){\displaystyle \frac{dw}{w}𝑑m𝑑\mu 𝑑\sigma p(F,w,m,\mu ,\sigma |𝐝,𝐈)}.`$ (A7)   
Substituting Eq. A4 into Eq. A7, we find that the integrals over $`\mu `$ and $`\sigma `$ may be done analytically. After some algebra, we find

$`f_F(F){\displaystyle \frac{dw}{w}𝑑m\left(D_2+\frac{F^2M^2}{w}\frac{2FMD_p}{w}\frac{(FMD_1)^2}{M}\right)^{(1M)/2}},`$ (A8)   
where constant factors have been eliminated and

$`D_1={\displaystyle \underset{i=1}{\overset{M}{}}}d_i,D_2={\displaystyle \underset{i=1}{\overset{M}{}}}d_i^2,\&D_p={\displaystyle \underset{i=m}{\overset{m+w1}{}}}d_i.`$ (A9)   
In the case that the width and/or phase of the pulse is known, the marginalization over $`m`$ and/or $`w`$ in Eq. A8 can be removed to further constrain the PDF for $`F`$.

To test our method, we have created many simulated profiles consisting of a pulse superimposed on Gaussian noise, and have confirmed that the PDF of Eq. A8 does maximize at the correct value of $`F`$. Some example simulation results for square wave pulses are shown in Figure 9. Figure 10 shows the results of applying our algorithm to profiles consisting only of Gaussian random noise.

## Appendix B Bayesian Flux Estimation of a Scintillating Source

We assume that a source’s intrinsic flux density is modulated by interstellar scintillations by a factor $`g`$, so that the measured flux density is

$`F=gS+N,`$ (B1)   
where $`N`$ is the additive noise and $`S`$ is the intrinsic flux density. To calculate a PDF for $`S`$ we assume we have $`m`$ measurements of $`F`$ and $`m`$ separate estimates of $`N`$, corresponding to making on and off-pulse flux measurements. We assume that the noise PDF is Gaussian (a good assumption for radio observations with large time-bandwidth product), and that we have estimated $`\mu `$, the average value of the off-pulse mean, and $`\sigma `$, the standard deviation of this quantity, for each of $`m`$ profiles. The PDF of the phase-averaged noise, $`f_N(N)`$, is a Gaussian with mean $`\mu `$ and standard deviation $`\sigma `$, and the PDF for an individual measurement of $`F`$ is

$`f_F(F)={\displaystyle 𝑑Nf_N(N)f_g\left(\frac{FN}{S}\right)},`$ (B2)   
where $`f_g`$ is the PDF of $`g`$. The likelihood function is a combination of all the data using PDF’s of this form. So, for $`m`$ statistically independent observations, the likelihood function for $`S`$ is

$`(S)={\displaystyle \underset{j=1}{\overset{m}{}}}f_F(F_j)`$ (B3)   
We cannot use this scheme for all cases because generally we need to use the joint PDF for the scintillation gain $`g`$ for the different measurements. Because there is no closed expression for this, we consider only two limiting cases.

### B.1 Case I: Statistically Independent ISS

The first case we consider is when the interstellar scintillations are statistically independent between measurements. For this case, we may write

$`_S(S)={\displaystyle \underset{j=1}{\overset{m}{}}}{\displaystyle 𝑑Nf_N(N)f_g\left(\frac{F_jN}{S}\right)}.`$ (B4)   
where $`f_N(N)`$ is a Gaussian described by $`\mu _j`$ and $`\sigma _j`$.

For $`\mathrm{\Delta }t_\mathrm{d}T`$ (i.e. the ISS time scale $`\mathrm{\Delta }t_\mathrm{d}`$ is much less than the averaging time $`T`$) and/or $`\mathrm{\Delta }\nu _\mathrm{d}B`$ (the ISS bandwidth is much less than the receiver bandwidth), the ISS is quenched and has a PDF that is related to a $`\chi ^2`$ PDF with the number of degrees of freedom given by $`2N_{\mathrm{ISS}}`$, where $`N_{\mathrm{ISS}}`$ is the number of scintles averaged over. In this case, statistical independence between observations is a good assumption.

### B.2 Case II: Perfectly Correlated Scintillations

The second case we consider is when the time scale for interstellar scintillations is much longer than the total data acquisition time (e.g. $`mT`$ for contiguous measurements). In this case, the scintillation gain $`g`$ is identical for all measurements, and has some unknown value with PDF described by Eq. 4. In this case, individual measurements are statistically independent for the noise but completely identical with respect to $`g`$. The likelihood function is therefore simply the product of the noise PDFs, and can be expressed as

$`_S(S)={\displaystyle \underset{j=1}{\overset{m}{}}}f_F(F_j)={\displaystyle \underset{j=1}{\overset{m}{}}}f_N(F_jS).`$ (B5)   
We can use this equation to get the likelihood function for the product $`S^{}gS`$ and then integrate over $`f_g(g)`$ to get the likelihood for $`S`$:

$`_S(S)`$ $`=`$ $`{\displaystyle 𝑑gf_g(g)_S^{}(gS)}`$ (B6)   
$`\mathrm{where}_S^{}(S^{})`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{m}{}}}f_N(F_jS^{}).`$ (B7)   
Once we have chosen the scintillation regime of interest, and calculated $`_S(S)`$ with Eq. B4 or Eq. B7, we may calculate the PDF of $`S`$ as the normalized likelihood function

$`f_S(S)={\displaystyle \frac{_S(S)}{_0^{\mathrm{}}𝑑S_S(S)}}.`$ (B8)   
When no source contributes to the ‘on-source’ measurement, we expect $`f_S(S)`$ to maximize at $`S=0`$ and have a width determined by $`\sigma `$ and the number of measurements. An upper bound on $`S`$ can be calculated by choosing a probability level $`1ϵ`$ and calculating the value of $`S`$ such that the area of $`f_S(S)`$ above that value is $`ϵ`$. If $`f_S(S)`$ maximizes at a non-zero value, then a confidence interval for $`S`$ can be similarly calculated.

When, as in the calculation of an upper limit using the method of Appendix A, we do not have a single measured value for $`F_j`$, but instead a PDF, we must replace $`F_j`$ by an integral over its PDF. In this case, Eq. B4 becomes

$`_S(S)={\displaystyle \underset{j=1}{\overset{m}{}}}{\displaystyle f_F(F_j)𝑑F_j𝑑Nf_N(N)f_g\left(\frac{F_jN}{S}\right)}.`$ (B9)