Could the “Wow” signal have originated from a stochastic repeating beacon?

Searches for strict periodicity can be readily achieved using classic periodogram techniques. Stochastic repeaters, in contrast, may evade detection using such techniques. Accordingly, we here seek to test the specific hypothesis of a stochastic repeater.

In what follows, we assume that the source maintains a uniform rate of signal production over any fixed interval of time within the baseline of the Big Ear observations. This essentially corresponds to a kind of Copernican Principle in time - we assume that the Big Ear observations do not occupy a special point in time when the signal source changed its emission mode (for example). Indeed, this assumption demands that the mean number of signals produced per unit time (in our direction) is the same when Big Ear began observing the source, as when it finished. This defines a Poisson process, for which the interval between individual signals will follow an exponential distribution. Accordingly, the $(k+1)^{\mathrm{th}}$ signal will occur at a time

where $\mathcal{E}[\lambda]$ is an exponential distribution of mean $\lambda$. For a given choice of $\lambda$, we thus begin our simulation by iteratively populating a list of signal mid-times using Equation 1, where we start with $t_{\mathrm{signal},0}=0$ but ultimately discard this initial entry.

The signal is also characterised by a duration. In what follows, we assume that the signal is a top-hat persisting for a time $T$ (with no variability), since this is both consistent with the Wow signal and represents the simplest possible model. If $t_{\mathrm{signal},k+1}-t_{\mathrm{signal},k}<T$, the two signals are effectively merged into a longer continuous one.

The observations are modelled using the observing logs of the Wow field from Big Ear itself. In September 1985, co-author Robert Gray visited OSU where the logs were reviewed and noted. Print outs were pulled from the OSU archives by Marc Abel, who located runs from 122 days spanning August 1977 to March 1984. Data near the Wow locale was examined for each day, with the locale defined as near RA 19 17 in 1977, and RA 19 22 and 19 25 subsequently, since a horn squint correction had presumably been applied by 1978.

Excluding observations more than 20’ north or south of the Wow declination, or days with bad data (e.g. clock off, printer failure, power failure), we find 90 useful visits, which are listed in Table 1.

These dates are used to define a set of representative observations in our simulation, with the first observation occurring at $t_{\mathrm{obs},1}=0.5$ days, and subsequent observations offset in time from this reference point according to their observation date. Each observation in fact a pair of measurements, since the Big Ear horns are physically separated and thus a sky source passes over each. Accordingly, we define the first of these two occurring at a mid-time of $t_{\mathrm{obs},l}$ (for the $l^{\mathrm{th}}$ observation) and the next occurring at $t_{\mathrm{obs},l}+\Delta t$, where $\Delta t=172.37$ seconds represents the temporal offset.

Each of the pair of observations is modelled as lasting for 72 seconds. For a “detection” to occur, we define that a signal must be throughout the entire duration of this 72-second window. We thus cycle through each observation ($90\times 2$ in total) and query whether any of the generated signals lead to a detection. After cycling through all observations, we tally up the total number of detections. Illustrative examples of both Wow-like and non Wow-like signals are depicted in Figure 1.

The process of generating a set of injected signals, comparing them to our simulated observations and tallying up the number of detections defines a single realisation. Since the signal generation process is stochastic, a statistical interpretation requires repeating many realisations for a given choice of $T$ and $\lambda$. In practice, we approach this by creating a “while” loop that generates many realisations and continues until exactly one detection is made - since this is what was found for the Wow signal. The greater the number of realisations necessary, the less likely it is that the particular choice of $T$ and $\lambda$ could reproduce Wow. Naively, one might reason that if it took $C$ realisations to get a Wow-like signal, then the probability of $T$ and $\lambda$ producing Wow is $1/C$.

For computational efficiency, we cap the maximum number of realisations in our while loop at 1000 and flag such cases appropriately. Although this single while loop provides a first estimate of $1/C$ probability, it is inherently uncertain given the fact only a single loop was executed. Thus, we wish to repeat the loop $W$ times, in order to improve our confidence in the estimated probability.

Since our while loop represents a binary outcome of fixed probability, then $C$ represents the number of Bernoulli trials needed to get one success, defining a geometric distribution i.e. $C\sim\mathcal{G}[p]$, where $\mathcal{G}$ represents the geometric distribution and $p$ is the probability value for this choice of $T$ and $\lambda$. We now wish to infer $p$ from a sequence of $C$ evaluations emerging from our while loops, which form a vector with indices $C_{w}$ of length $W$. To accomplish this, we note that

where $C_{\mathrm{tot}}\equiv\sum_{w=1}^{W}C_{w}$. To infer $\mathrm{Pr}(p|\{C_{1},C_{2},...,C_{W}\})$ (the posterior distribution of $p$), we need to use Bayes theorem. Adopting a uniform prior on $p$ for simplicity, the posterior is simply the normalised likelihood given by

Differentiating, this has a mode (i.e. maximum a-posteriori, MAP) value of

which is close to our original intuition of $1/C_{\mathrm{tot}}$ for $W\sim 1$. However, this process also allows us to define the variance to be

Using these results we can say that after $W$ while loops, each of which took $C_{w}$ realisations to obtain a Wow-like signal, we infer the MAP likelihood of the choice of $T$ and $\lambda$ to be $\bar{p}$ with a precision (fractional error) of $\bar{p}/\sigma_{p}$. We thus continue through the while loops until the precision falls below 10% - which we qualify as a “good” estimate of $\bar{p}$.

From before, $\bar{p}$ represents the maximum a-posteriori probability of obtaining a Wow-like signal with the 90 days of Big Ear data, given a choice of $\lambda$ and $T$. In moving towards inferring $\lambda$ and $T$, $\bar{p}$ can be considered to be the (most probable) likelihood of obtaining the data $\mathcal{D}$ (here defined as a single detection in the 90 pairs of observations) given $\lambda$ and $T$. Thus, we define

To make progress, we need to be able to vary $\lambda$ and $T$ and evaluate likelihoods quickly, to exploit Bayesian inference methods. Given the computational expensive nature of even a single likelihood evaluation (sometimes taking many hours), we elected to build a grid of $\lambda$ and $T$ values across which we could pre-calculate the likelihoods and then use them to build a likelihood emulator at intermediate values. To this end, we set up a $M=100$ square grid of plausible values for each. The grid was chosen to be log-uniformly spaced from $\lambda_{\mathrm{min}}$ to $\lambda_{\mathrm{max}}$ and $T_{\mathrm{min}}$ to $T_{\mathrm{max}}$, such that

where $i$ and $j$ are the integer indices running from 1 to $M$. We set $T_{\mathrm{min}}=72$ seconds, which defines the minimum allowed duration given that the signal shape appears consistent with a continuous emission. For $T_{\mathrm{max}}$, we chose a maximum at least two orders of magnitude larger and specifically adopted $T_{\mathrm{max}}=e^{-1}$ days. For $\lambda_{\mathrm{min}}$ and $\lambda_{\mathrm{max}}$, we chose to bracket them three orders of magnitude either side of $1/T_{\mathrm{max}}$, giving $\lambda_{\mathrm{min}}=e^{-4}$ days^-1 and $\lambda_{\mathrm{max}}=e^{2}$ days^-1.

Our later parameter explorations found that this initial grid did not support certain regions of high likelihood space and needed to be expanded. To accomplish this, we simply increased the indices $i$ and $j$ beyond 100. Specifically, we ran up to $i=150$ and $j=160$. This means that in practice the maximum allowed $T$ value is not in fact $T_{\mathrm{max}}$, but rather $T_{\mathrm{max}}^{\prime}$, and similarly for $\lambda$, where

which evaluates to $T_{\mathrm{max}}^{\prime}=7.97$ days and $\lambda_{\mathrm{max}}^{\prime}=280.4$ days^-1.

Since our grid is regular and evenly spaced (in log space), it is amenable to interpolation. To this end, we use a bi-cubic spline interpolation of the grid. As shown in Figure 2, this leads to an accurate and quick emulator across the high likelihood region of interest.

Having established a likelihood emulator supported over the range $T_{\mathrm{min}}\leq T<T_{\mathrm{max}}^{\prime}$ and $\lambda_{\mathrm{min}}\leq\lambda<\lambda_{\mathrm{max}}^{\prime}$, we now turn to inferring the joint posterior distribution of $T$ and $\lambda$, given the data $\mathcal{D}$. This can be achieved through Bayes theorem:

For the prior, $\mathrm{Pr}(T,\lambda)$, we assume that the priors are independent and follow log-uniform distributions over their supported range. This is most easily achieved by simply walking in $\log T$ and $\log\lambda$ space directly, which is also convenient as this is the native sampling the grid used for the likelihood emulator. The limits on the walkers are set to that of our likelihood grid, since outside this region we cannot reliably emulate the likelihood.

We sample the joint posterior distribution using Markov Chain Monte Carlo algorithm and Metropolis walkers, calculating the posterior probability as a product of our likelihood emulator and the prior (which is uniform in this parameterisation). Each walker completes 110,000 accepted steps, burning off the first 10,000. The starting point for the chain is randomly drawn from the joint prior. We repeat for 100 walkers and combine the final results, with a thinning of a factor of 10 to leave us with $10^{6}$ posterior samples. Chains and posterior distributions were inspected to ensure convergence and good mixing.

The joint posterior distribution is depicted in Figure 3. We find that the marginalised duration distribution peaks close to the lower boundary condition. Since MCMC samplers cannot make jumps below the prior minimum threshold, a slightly positively skewed posterior is expected even when the result is consistent with the minimum. Accordingly, we here conclude that there is no clear lower limit on the signal duration, besides from the constraint that it appears consistent with being fully-on during the window in which it was caught (i.e. $T>72$ seconds). However, we find that long duration signals are disfavoured, which can be understood by the fact that they are increasingly unlikely to have only been detected once during the ensemble of Big Ear observations. Indeed, from this, we can place an upper limit that $T<77$ minutes to $2$ $\sigma$ confidence (95.45%).

For $\lambda$, the mean rate of signal emission, we find a broad but peaked posterior with a maximum at $\lambda=2.01$ days^-1 (i.e. twice per day). We can see that $\lambda$ is inversely correlated to $T$ but compresses in a narrow region at low $T$. This “fine-tuning” reduces the marginalised posterior density in this region, and ultimately gets truncated by the lower boundary on $T$. Very frequent repeat times are only compatible with the data at extremely short durations, in order to explain the lack of additional detections. Together, these constraints imply that $\lambda$ is bound within the interval $0.047$ days^-1 to $64.2$ days^-1 to 2 $\sigma$.

Of course, besides from the original Big Ear observations, there have been several efforts to re-observe the field and look for repetitions, without success. In Gray (1994), the field was observed using the ultranarrowband Harvard/Smithsonian-META system (Horowitz et al., 1986), tracking for up to 4 hour at a time. Gray & Marvel (2001) performed the second attempt with the Very Large Array, tracking the field for no more than 22 minutes a time. Gray & Ellingsen (2002) performed far longer observations with the University of Tasmania’s Hobart 26 m radio telescope, tracking the field for 14 hours continuously on six separate days. Finally, the Allan Telescope Array was used by Harp et al. (2020), accumulating over 100 hours on target for 30 to 180 minutes each time. In all cases, we note that the observations were sufficiently sensitive to have recovered a repeat of the Wow signal to high confidence.

Amongst the four follow-up efforts highlighted above, the Gray & Ellingsen (2002) data with Hobart is most amenable to incorporating into our analysis, since the observing window is approximately the same each time (14 hours) and only six sets were taken. The dates of the observations are provided in Gray & Ellingsen (2002) and thus we could easily imagine directly extending our likelihood formalism to include these dates. However, a much simpler approach would be to ignore the dates and rather just assume that each 14-hour window is an independent and fair sample over all time within which no detection was made. For a $F=14$ hour block, the probability of a Poisson process not yielding a single success equals $e^{-F\lambda}$, and to evade detection in six independent draws would be $(e^{-F\lambda})^{6}=e^{-6F\lambda}$. This serves as an additional likelihood term which we can append to the Big Ear likelihood emulator. With this simple change, we repeated the MCMC analysis using an otherwise identical setup.

The updated chain found that the maximum likelihood dropped from 32.3% using Big Ear alone, to 3.27%, thus placing tension on the stochastic repeater hypothesis at the 2.1 $\sigma$ level. This peak occurs at $\hat{\lambda}=0.233$ days^-1 and $\hat{T}=485$ s. Evaluating the likelihood of this position using the Big Ear likelihood emulator alone (ignoring the Hobart data) yields 7.40%. Thus, the likelihood shifts from 7.40% to 3.27% - a factor of 0.442 of the original value, which of course simply equals $e^{-6F\hat{\lambda}}$. Using this point of highest likelihood, we can test the accuracy of our assumption of independence.

To this end, we used the dates of the Hobart observations and generated 1000 Poisson processes with a mean rate $\hat{\lambda}$ and duration $\hat{T}$ and counted how many detections Hobart would have seen. From this, we measure 439 zero-count cases - consistent with the real observations of Gray & Ellingsen (2002). Thus, we find that a more realistic calculation that accounts for the specific timings of the observations yields a Hobart-likelihood at $\hat{\lambda}$ and $\hat{T}$ of $(0.439\pm 0.016)$ - a value fully consistent with the much simpler independence assumption. On this basis, we argue that our simplifying assumption is justified and produces an accurate final posterior.

From our revised joint posterior, we obtain credible intervals of $\lambda=[0.021,0.441]$ d^-1 and $T=[130\,\mathrm{s},2.2\,\mathrm{d}]$.

Our approach for incorporating the Hobart observations can be extended to other data sets too. Recall that for Hobart, the additional likelihood term was $(e^{-F\lambda})^{6}=e^{-6F\lambda}$, since 6 sets of $F$ baselines were taken. This is equivalent to adding a single penalty of $6F$ duration (i.e. 84 hours). Accordingly, even if the other observations did not use the same observational window each, we can still include them by simply counting up the total observing time of each (on the Wow field).

From Gray (1994), we find 8 hours on each of the two possible Wow positions. Gray & Marvel (2001) dwelled no longer than 22 minutes on the field and used a variety of observing modes, and thus we elected to not include these data as their temporal constraining power is small compared to Hobart. In contrast, the ATA 100 hour campaign from Harp et al. (2020) more than doubles the baseline and thus provides a key constraint. Together then, we add a third likelihood term equal to $e^{-G\lambda}$, where $G=108$ hours.

With this change, the final maximum likelihood is 1.78% (2.37 $\sigma$), occurring at $\hat{T}=659$ s and $\hat{\lambda}=0.121$ days^-1. The final joint posterior is plotted in Figure 4, where we report a revised 2 $\sigma$ credible interval of $T=[130\,\mathrm{s},1.53\,\mathrm{d}]$ and $\lambda=[0.021,0.397]$ d^-1. In practice, the joint posterior is not significantly changed by the inclusion of the other observational constraints.