All-Sky Monitor

The EPIC imager forms an image of the whole sky above the horizon every 83.6 ms in each of its 128 channels. The all-sky monitor (skymon, in src/python/skymon of the LWA_EPIC repository) searches these images for dispersed radio bursts at every position above an elevation of 10° in real time.

Each imaging node runs its own monitor for its own band of 3.06 MHz. A monitor depends on no other node, on neither the database nor the shared file system for its search, and on no other monitor. The loss of a node therefore removes only that band.

Data path

The imager writes each image into the next slot of the ring and never waits for a reader. A monitor that falls behind loses the oldest images and records the loss; a restart of either program does not affect the other. The pixels of the published images are in the order of the saved FITS images, so the monitor uses the imager's celestial coordinate system (epic_geometry.apply_wcs) unchanged.

Sky-fixed pixels

A source drifts across the image by one pixel every few minutes. A dispersed burst at low frequency can take longer than that to cross the band. The time series of a fixed image pixel would therefore follow different parts of the sky over the duration of one burst.

The monitor instead follows HEALPix pixels (J2000, RING order). The tracked set is recomputed every 10 s:

• a pixel that rises above the elevation limit is assigned a slot in the state arrays;
• the slot of a pixel that sets is cleared and returned to the pool.

Between updates, the position of each tracked pixel on the image is interpolated linearly in time, so its sample of the sky moves continuously with the rotation of the Earth. Each image is interpolated at these positions by bicubic convolution (Keys 1981). With about three image pixels across the beam, bilinear interpolation would miss the peak of a source by up to 6% depending on its position between pixels, and the error would ripple at the same time in every channel as the sky turns; bicubic interpolation reduces it to about 2%.

The HEALPix resolution is the coarsest at which a pixel is no larger than half of the synthesised beam $\lambda / D$ ($D \approx 100$ m), and never finer than the image pixel. This gives:

• $N_\mathrm{side} = 64$ ($0.92°$) above 40 MHz, about 20 300 pixels tracked at once;
• $N_\mathrm{side} = 32$ near 35 MHz;
• $N_\mathrm{side} = 16$ at 16–21 MHz.

Normalisation and excision of interference

The value of each pixel and channel is normalised by a running mean with a time constant of 2 min and a running variance with a time constant of 5 min. The variance is estimated as half the mean square difference between consecutive images, which a slow change of brightness, such as that of a source rising through the primary beam, does not inflate. Deviations beyond five standard deviations are limited before they enter these statistics, so that interference and bursts do not raise the noise estimate.

For each channel and image, the median over the tracked sky is then subtracted. Interference received through the sidelobes and changes of receiver gain affect the whole image, whereas a celestial burst occupies one beam, so the median removes the former and leaves the latter unchanged.

After the median, the spread of each channel over the sky (1.4826 times the median absolute deviation) is close to 1 for noise. A narrowband carrier received over much of the sky, such as the transmitter of an aircraft or a distant transmitter reflected by a meteor trail, widens it; a burst, confined to one beam, does not. A channel whose spread exceeds 1.6 times that of the band, the median spread over the channels or 1 if smaller, is set to zero for every pixel of that image. A rise of noise over the whole band, as from a thunderstorm or a solar burst, widens every channel alike; it is not flagged, but searched and recorded in the conditions of the band. A channel flagged in more than half of the images over 30 s is masked until its rate of flags falls below 0.2.

A single sample beyond 12 standard deviations is set to zero. A burst of S/N $s$ adds about $s/\sqrt{N_c\,n}$ to one sample, where $n$ is the number of samples per channel in its sweep, which is below 9 for $s = 100$ and 128 channels; the limit therefore removes only values no burst produces. Flagged channels and samples set to zero do not update the running statistics, except that a sample withheld for 30 s in succession updates them again, with its deviation limited to five standard deviations, so that a lasting change of level is eventually followed.

Hierarchical dedispersion

The dispersion delay across the band between $f_\mathrm{lo}$ and $f_\mathrm{hi}$ is

At 57 MHz the delay is 0.135 s per unit of DM, so a burst at DM 3000 takes 6.7 min to cross the band. At 16 MHz the same burst takes 5.2 h.

The smearing within a channel equals the sweep across the band divided by the number of channels $N_c$. The search exploits this in tiers:

• Tier $j$ averages $2^j$ images.
• Each tier computes the fast dispersion measure transform (FDMT; Zackay & Ofek 2017) over delays of $0$ to $D = N_c$ samples at its own resolution.
• Because the delay is proportional to DM, every tier uses the same delay table, and tier $j$ covers DMs $2^j$ times those of tier 0.
• A tier above tier 0 reports only delays above $D/2$, so the tiers cover the DM range in octaves.
• At the upper end of each tier, a pulse spans about one sample per channel.

The history a tier needs is $D$ samples plus the block being searched, independent of DM. For 128 channels the search state of all tiers takes about 0.7 MB per pixel, or 14 GB for the whole sky at $N_\mathrm{side} = 64$; the 6-hour spectrogram accumulation takes a further 12 GB. Storing the full-resolution series instead would take 222 000 samples per channel at 16 MHz and DM 3000, or about 1 TB for the sky.

Whether a burst can be detected at a given DM depends on how long its position stays above the elevation limit. A burst whose sweep ends after its position has set is not completed. Above 30 MHz nearly every position in view remains up long enough for DM 3000. At 16–21 MHz, positions near the western horizon limit the reach to a few hundred.

Detection and classification

Before each transform, a straight line fitted over the window is subtracted from every channel of every pixel, which removes the remaining slow trends; a burst, which occupies about one sample per channel, loses about $1/W$ of its signal for a window of $W$ samples. Above tier 0 the mean over channels at each time is also subtracted (a zero-DM filter; Eatough et al. 2009). This removes any change common to all channels, such as a gain variation or the residual interpolation ripple of a bright source, while a burst in these tiers occupies one or two channels per sample and loses one or two parts in $N_c$. Tier 0 is not filtered in this way, so that bursts and pulsars of low DM keep their full signal. Each tier searches boxcars of 1, 2, 4 and 8 samples at every delay. The noise of every delay and boxcar width is measured from the data themselves, as the clipped root mean square of the previous blocks, so the signal-to-noise ratio does not assume white noise. Values above S/N 8 are retained. They are then:

• clustered over the sky (twice the beam width), over DM and over arrival time;
• each cluster becomes a candidate at its brightest pixel.

A candidate is classified as interference in any of the following cases:

• its cluster extends over more than six beam areas;
• more than 1% of the tracked sky triggers within its time span;
• clusters at four or more well-separated positions coincide with it;
• its DM is below 2 pc cm⁻³;
• its signal is confined to a few channels.

The last test separates a dispersed burst, which spreads over most of the band, from a carrier of an aircraft or of meteor scatter, which occupies one or two channels of 24 kHz. Along the dispersion track, the samples of each channel are summed over the boxcar at that channel's delay. From these contributions $c_i$ (positive part), the effective number of channels is $(\sum c_i)^2/\sum c_i^2$. A candidate with fewer than 6 effective channels, or with more than 60% of the sum in two channels, is classified as narrowband interference, and its data are kept.

A candidate at the position of a known source is stored with its data in the class known, with the kind, name and separation of the source:

Known candidates do not open events. Those not archived are removed after 14 days.

A burst detected again by another tier or block is stored once, as the detection of highest S/N.

Data products

Every candidate is written to <out>/candidates/YYYYMMDD/ as a JSON record. The record holds:

• position and pixel;
• DM, S/N and width;
• arrival times at infinite frequency and at both band edges;
• elevation and azimuth;
• tier;
• the pixels of its cluster.

Every candidate not classified as interference also has an .npz file holding:

• the dynamic spectrum of its pixel at the resolution of the tier;
• the DM–time plane of that pixel;
• the accumulated spectrum of the 30 min before the burst;
• the noise table of the search.

Every 10 s of data, each monitor records the conditions of its band:

• the noise relative to its baseline (the median over pixels and channels of the running variance, divided by its average over 1 h);
• the rates of interference and of all candidates over the last minute;
• the fraction of the tracked sky that triggered;
• the fraction of images missing.

A band is marked disturbed when the noise exceeds 1.5 times its baseline, interference exceeds 2 candidates per minute, or more than 1% of the sky triggers. Thunderstorms, solar activity and episodes of interference produce such periods, lasting minutes to hours. The samples of the last 6 h are served by the monitor and written to its local database. Each candidate carries the conditions of the minute before its arrival, and each decided event those of the bands that detected or measured it.

Calibration and flux

When the imager is calibrated by the streaming calibrator (epic++ --gain_table), each image carries the version of the gain table it was made with. The monitor counts the changes of version in the conditions of its band. Each candidate records the version of its images and whether the version changed during its sweep, in particular whether calibration started or stopped then, which changes the images more than a routine update.

The gain tables are relative: in each channel their median amplitude is one, and the imager applies no primary beam, so the images carry no absolute flux scale. With every solution the calibrator therefore publishes a record of it:

• the absolute scale removed from the table, per channel and polarisation;
• the compact sources of the solution, with their model spectra, the flux scales solved for them, their positions and the response of the dipoles towards them;
• a table of the dipole response over the sky;
• the quality of the solution.

Every minute the monitor measures each calibrator above 20° in its own images: the excess of the calibrator's pixel over a ring two to four beams around it, in the spectrum averaged over 10 s. Divided by the calibrator's flux and the response towards it, this gives the image units per jansky in every channel, independently of the normalisation of the imager. The fluence of a candidate in each channel then follows from its signal along the dispersion track, the noise of its pixel in image units and the response towards its position:

where $x_c$ is the normalised signal in tier $j$, $\sigma_c$ the pixel's noise per image, $\Delta t$ the image length and $k_c$ the measured image units per jansky. Each candidate stores its fluence over the band and per channel, its mean flux density over the boxcar, and the calibrators used. Its data products hold the noise of its pixel, the calibration record and the calibrator measurement, so that the fluence can be recomputed. On a simulated sky with a record in the simulator's units, the measurement of Cyg A gave 0.98 of the true scale, and a burst of known fluence was measured at 0.82 of its fluence, close to the loss in its S/N (34 found for 40 injected): the burst lies off the centre of its sky pixel, and the filters of the search remove part of its signal.

The spectrum of every tracked pixel, averaged over 10 s, is retained for 6 h. It serves the All-sky monitor page of the Waterfalls site. The page shows the sky of a chosen monitor, the mean of the last 5 min of its accumulation on its HEALPix grid, beside the state of every monitor. Any pixel of that sky can be selected for its spectrogram, either the latest search window at the image cadence or the 10 s averages of up to 6 h. The candidates and the events across bands are listed below.

Test injections

At random times, with a mean interval of one hour by default, a monitor adds a simulated burst to an incoming image before interpolation. The burst has:

• a position that stays above the elevation limit for the whole sweep;
• a DM drawn log-uniformly over the range that position can still complete;
• a target S/N between 10 and 30, set from the measured noise of the pixel.

A candidate that matches an injection is classified as a test injection. Every injection is recorded when it is scheduled and again when it is recovered or missed (<out>/injections.jsonl), so the detection efficiency of the running system is measured continuously. The injections exercise every stage after the imager.

Verification with a simulated sky

The monitor includes a simulator of the sky in the imager's format. It produces images at the imager's rate, with:

• diffuse Galactic emission;
• the four brightest radio sources;
• a primary beam;
• radiometer noise correlated on the scale of the beam;
• whole-sky and narrow-band interference;
• FRBs at random positions, DMs, widths and S/N, recorded in a truth file.

The simulator writes into a shared-memory ring as the imager does, and the monitor reads it unchanged.

The run summarised below covered 30 min of the 55.88–58.94 MHz band. It used:

• 128 channels and 83.6 ms images;
• $N_\mathrm{side} = 64$, with 20 310 pixels tracked;
• DM up to 3000 pc cm⁻³;
• eight threads on cores 24–47 of an imaging node (AMD Threadripper PRO 5965WX);
• 40 FRBs per hour, with DM 20–3000 pc cm⁻³ and target S/N 8–30;
• RFI flashes at 6 per hour;
• a live test injection every 5 min.

A second run lasted 60 min with the same settings, except that the FRBs were drawn at 30 per hour with DM 300–3000 pc cm⁻³ and target S/N 10–30, and live injections came every 10 min.

A tier begins reporting once its window has filled twice: after about 21 s for tier 0 and 8.5 min for tier 5. This happens once, when a monitor starts.

Consensus across bands

Each monitor decides alone whether a candidate exceeds its threshold. A burst from the sky must, in addition, appear at the same position and DM in every band that can see that position, at the times given by the dispersion law:

Across the bands of LWA-SV these times differ by up to hours: at DM 500 the burst reaches the bottom of the 74–77 MHz band 6 min after $t_\infty$ and that of the 33–36 MHz band 31 min after it. Interference at the station does not follow this law. The combiner (combiner.py) uses this difference.

The combiner runs as a single process on one imaging node. Every 10 s it reads the candidates of every monitor. A new candidate opens an event, or joins an open event whose position, DM and predicted arrival it matches within the uncertainty of the DM extrapolated across the bands.

A candidate at the position of an event whose arrival falls within the event's sweep through its band, but at a clearly different DM, is listed with the event as another DM at this position, with a note. Such a candidate may be the same burst found at a neighbouring trial DM, or a second burst from the same place, such as from a repeating source. Its data remain linked to the event and are archived with it, but it neither replaces the band's detection nor changes the fitted DM; the reviewer decides between the two cases.

For every other band the combiner:

1. predicts the arrival at the top and bottom of the band, with a time tolerance equal to the extrapolated DM uncertainty;
2. determines whether the position is above that monitor's elevation limit throughout the sweep;
3. sends that monitor a watch: a request to measure the S/N at the position, DM and time, whatever its value.

The monitor measures a watch when the tier covering the DM searches the predicted time. It uses the pixel nearest the position and its neighbours, the delay rows within the DM uncertainty, and the times within the tolerance, and stores the result with the pixel's data. A burst seen first in the high bands is therefore measured in the lower bands at a threshold of S/N 6, a level the blind search over the whole sky could not use.

With detections in two or more bands, the DM and $t_\infty$ are fitted to all arrival times together. The lever arm between the bands determines the DM far more closely than a single band of 3 MHz.

A burst need not be detected in every band that can see its position. Its spectrum may be steep, or confined to a few MHz, as is common for repeating FRBs, and scattering, which grows as $\nu^{-4}$, can spread the pulse at the lowest frequencies far beyond the search's boxcars, so that the S/N falls although the fluence does not. Each band of an event therefore contributes a value rather than a vote:

• the S/N of a detection, or of a followup measurement, whatever its value, with an uncertainty of 1;
• an upper limit at the band's threshold, where the burst had passed the band before the event was known and only the band's own search covered it.

These values are compared (spectrum.py) with a power law of spectral index between −10 and +10, with emission confined to a Gaussian 5–30% wide in frequency, and with either suppressed by scattering of any strength. The fit states the spectra consistent with the measurements, for example a power law of index −2 to +3, or emission confined near 55–60 MHz.

An event is decided when every band is resolved:

With the bands of LWA-SV several MHz apart, almost any burst seen in one band fits emission confined near that band or just beyond the observed range. An implausible spectrum therefore indicates a contradiction between bands that cover nearly the same frequencies.

The combined S/N is the sum of the S/N measured in every band, divided by the square root of the number of bands; the largest single-band S/N is also given, since a band-limited burst is bright in one band only. Events are stored in the local databases of the imaging nodes (see Storage) and served on port 8096. A detection arriving after an event is decided reopens it.

Terrestrial coincidences

Impulsive interference of terrestrial origin, such as lightning sferics, meteor afterglows and electrical discharges, is not dispersed: it reaches every band at the same moment, whereas a burst from the sky reaches the bands seconds to hours apart. The combiner retains the candidates and interference of every band for 15 min. Two candidates from different bands are simultaneous when their arrivals at the bottom of their bands differ by less than 1.5 s plus their sampling intervals and widths. A simultaneous pair is undispersed when the dispersion law, at the DM of either candidate, predicts a gap between the two bands of more than 3 s plus the same tolerances. Positions are not compared, since interference in the near field is imaged at different positions in different bands. For bands close in frequency, or at low DM, the predicted gap is too small to distinguish the two cases, and such pairs are left to the consensus.

The candidates of all bands involved form one event of the kind terrestrial, which needs no review. A coincident candidate opens no event from the sky. An open event whose detections are all coincident is given the status terrestrial, with a note, and its record is kept.

External context

Solar radio bursts and nearby thunderstorms are the main natural sources of false events. When an event is decided or reopened, the combiner attaches for the reviewer:

• the elevation of the Sun at LWA-SV at the first arrival, and its angle from the event's position;
• the reports of NOAA SWPC within 15 min of the arrival in any band: radio sweeps (types II, III, IV and others) overlapping the band in frequency, and X-ray flares, each with its class and its times of beginning, maximum and end;
• lightning strikes within 10 min of the arrivals, as the number within 100 km and the distance of the nearest, when a source of lightning data is configured.

Each source is polled in its own thread every 15 min, and its last data are retained on disk. A failure of a source is recorded once and does not affect the combiner. The context of an event decided within the last 24 h is updated when new reports arrive, since the reports of SWPC are revised hours after the events.

Storage

Nothing in the real-time path depends on the shared file system, whose NFS mounts on the imaging nodes have been lost repeatedly, nor on the head node that serves it. Each imaging node keeps on its local disk:

• the monitor's code and a copy of its Python environment;
• the data products of every candidate;
• a PostgreSQL cluster of its own (port 5433), separate from the node's Citus worker, whose data are on the shared file system.

The local database holds the record of every candidate and serves as the node's mailbox. The combiner reads new candidates from it and writes its requests into it (watches and coordinated injections), which the monitor picks up within seconds and answers in the same tables. A request therefore survives a restart of the monitor. If the local database is unavailable, the monitor continues to search and writes its records to a spool directory on local disk, from which they are inserted when the database returns.

The combiner copies every event into the database of every node it can reach, with a version number. A combiner started on another node rebuilds its state from the newest copies. The combiner holds no state elsewhere.

The shared file system is the archive. On each node an archiver, a process separate from the monitor, copies to it the records and data of the node's candidates that belong to events worth keeping:

• confirmed events;
• events reviewed as an FRB or as unclear;
• coordinated injections.

It also archives the events themselves. Each copy is made in a child process with a time limit, compared byte for byte with the original, and only then marked as archived. When the archive is not mounted the archiver waits, and the data remain on the node.

Local disk is freed only by fixed rules:

Data that might be a burst and have not been archived are never removed.

Failures

Retention of data and review

No result causes data to be discarded. The candidate and its data products are kept in every band:

• the detection with its dynamic spectrum and DM–time plane;
• each followup measurement with the same products;
• for a band the burst had already passed when the event was known, the accumulated spectrum of the position for 30 min either side of the predicted arrival.

Each monitor also keeps the data of candidates classified as RFI, unless the RFI is certain: a DM below 2 pc cm⁻³ (no dispersion) or a moment at which more than 1% of the sky triggers at once. Every candidate carries measures for its review within its own band: the S/N against trial DM, the ratio of the S/N at the lowest trial DM to the best S/N, the S/N in each quarter of the band along the dispersion track, the effective number of channels, the fraction of channels flagged along its sweep, and the conditions of its band.

Archived candidates and events are kept on the shared file system for five years. An event awaiting review, and the data of its candidates, are never removed.

Every event from the sky is marked for review, confirmed events first and events with an implausible spectrum last. A reviewer logged in to the Waterfalls site records a verdict (FRB, RFI or unclear) with an optional note on the All-sky monitor page; the verdict, the reviewer and the time are stored with the event.

At random times, every 6 h on average, the combiner requests the same burst from every monitor: one position, DM and $t_\infty$, each band setting the amplitude from its own noise. These coordinated injections measure the efficiency of the whole chain, including the consensus.

The monitors do not depend on the combiner. If the combiner or its node stops, each band continues to search and store its own candidates. The combiner keeps its state in the local databases of the imaging nodes, so Nomad restarts it on any other imaging node, where it reopens the events still open.

The consensus was verified with simulators of several bands that received the same list of FRBs, each band running its monitor in real time on its own imaging node and the combiner on another.

In the three-band run, the low-DM interference of one band (DM 3–8) remained an event of one band. In the four-band run, the events matching no simulated FRB were mostly detections of bright bursts at a neighbouring trial DM, which the combiner lists with the burst's own event as other DMs when they fall within its sweep at the same position.

Operation

The monitor runs as the Nomad system job epic-skymon, one allocation per imaging node. The archiver runs as the system job epic-skymon-archiver, also one allocation per node. The combiner runs as the service job epic-skymon-combiner, with a single allocation on one of the imaging nodes. The local databases, directories and software are installed by the playbook 25_Skymon_local_storage.yml. It runs at reduced priority on processor cores separate from the imager, with eight compute threads, and serves its HTTP interface on port 8095. The imager publishes its images when started with --sky_shm epic_sky0; until then the monitor waits for the ring.

PYTHONPATH=LWA_EPIC/src/python python -m skymon run --shm epic_sky0 --out skymon_data --port 8095 --threads 8
PYTHONPATH=LWA_EPIC/src/python python -m skymon simtest --realtime --duration 1800 --out simtest