Bursts from Space: MeerKAT – the first citizen science project dedicated to commensal radio transients

ABSTRACT

1 INTRODUCTION

The latest generation of radio telescopes provides us with unprecedented detail of the radio sky. Regular, wide-field images from highly sensitive telescopes, including Square Kilometre Array (SKA) pathfinders such as MeerKAT (Jonas & MeerKAT Team 2016) and the Australian SKA Pathfinder (ASKAP; Hotan et al. 2021), allow us to probe a wide range of physics at novel time-scales and depths. For example, stellar radio emission can provide insight into magnetic re-connection (Rigney et al. 2022) and has implications for orbiting planet habitability (Airapetian et al. 2017; Günther et al. 2020), while the radio afterglows of gamma-ray bursts constrain jet physics and kinetic feedback of the most violent eruptions in the Universe (e.g. Rhodes et al. 2021). Similarly, fast radio bursts can probe the baryonic content of the Universe (Macquart et al. 2020), while the afterglows of neutron star mergers can provide key insights into the structure of relativistic jets (Margutti & Chornock 2021). The combination of sensitivity, regular cadence, and (crucially) wide field of view (FoV) allows for commensal, untargetted strategies in order to search for these known transient phenomena and new classes of objects as yet undiscovered.

Previous investigations from both MeerKAT and other instruments have found that only a few per cent of point sources are transient or variable above sensitivity limits at 1.4 GHz (Ofek et al. 2011 and references therein), with source classes spanning a wide range of time-scales and physical processes. The majority of radio variables found are active galactic nuclei (AGN; see e.g. Thyagarajan et al. 2011), whose variations can be attributed to a combination of refractive scintillation (Rickett 1990) and shock-induced flaring in their jets (Mooley et al. 2016). While these AGN dominate samples of variables, active or flaring stars have been found in untargetted radio surveys (Mooley et al. 2016; Driessen et al. 2020; Andersson et al. 2022), as have supernovae and GRB orphan afterglow candidates (Levinson et al. 2002; Gal-Yam et al. 2006). Pulsars can vary intrinsically in the image plane – indeed some of the slowest known pulsars are discovered in imaging data (Tan et al. 2018; Caleb et al. 2022). Futhermore, diffractive scintillation through the interstellar medium (ISM) can cause short time-scale, large amplitude variations in observations of pulsars, while refractive scintillation can produce lower amplitude variability occurring on time-scales of hours to years for point-like sources (Rickett 2001; Hancock et al. 2019). There are also numerous accounts of radio transients being discovered without clear progenitor systems or multiwavelength counterparts (Bower et al. 2007; Stewart et al. 2016; Murphy et al. 2017). These include the elusive sources near the Galactic centre (Davies et al. 1976; Zhao et al. 1992; Hyman et al. 2005; Chiti et al. 2016) – including the newest such transient found by ASKAP (Wang et al. 2021). The serendipitous discoveries, elusive nature, and broad physics at play in this zoo of radio transients all point towards the need for new searches and the development of novel methods to maximize the science yield of our observations.

ThunderKAT¹ (Fender et al. 2016) is a large survey project dedicated to monitoring radio transients with MeerKAT. The ThunderKAT team regularly observes known transients such as X-ray binaries, cataclysmic variables, and gaama-ray bursts (XRBs, CVs and GRBs respectively). The large FoV of MeerKAT (∼1 square deg at 1.28 GHz) sampled at approximately weekly cadences with high sensitivity also allows for unprecedented commensal searches for transients, variables, and other ancilliary science. Driessen et al. (2020, 2022a) and Andersson et al. (2022) describe the first commensal transients found with MeerKAT, detailing flaring and quiescent behaviour from stellar systems. Similarly, Driessen et al. (2022b) and Rowlinson et al. (2022, hereafter D22 and R22, respectively) use the best-sampled ThunderKAT fields surrounding XRBs GX339−4 and MAXI J1820+070 to discover new radio variables including pulsars and variable AGN. Images of short GRB fields have been searched for transients at both fast and slower time-scales in Chastain et al. (in preparation), wherein there are many newly described variables, of which most are likely scintillating AGN. The ThunderKAT survey also makes use of MeerLICHT (Bloemen et al. 2016), a robotic facility (65cm primary mirror) whose goals include shadowing MeerKAT observations, providing spatial and temporal coverage of the optical sky to complement the radio data.

Despite these methods of searching for radio transients bearing fruit, they are not optimal. First, the volume of data to analyse is far greater than any one person can achieve by eye on reasonable time-scales. In the 4 yr since operations began, ThunderKAT has observed over 30 XRBs in total at weekly cadence, typically following each source for over a month. This results in over 100 TB of raw data to reduce and analyse, producing over 500 final images. Each image then contains of order several hundred sources, from just a single 15 min observation. ThunderKAT also has memoranda of understanding with many of the other Large Survey Projects on MeerKAT, such as LADUMA (Blyth et al. 2016), MIGHTEE (Jarvis et al. 2016), and MHONGOOSE (de Blok et al. 2016), to use their data commensally. As a result there are many hundreds of observations in the growing archive, in which transients may reside, probing right down to 1σ sensitivity limits ∼1 μJy (Heywood et al. 2022). These data overload issues are only exacerbated when imaging on shorter time-scales, including down to the 8 s integration time of MeerKAT, as is currently being tested within ThunderKAT (e.g. Caleb et al. 2022; Chastain et al. in preparation; Fijma et al. in preparation).

Radio observations are not free from false positives. Two main causes of these false positives are the non-Gaussian artefacts that typically occur around bright sources in radio images, and the changes in the point spread function (PSF, or restoring beam) caused by differing elevations over a set of observations, which induces non-intrinsic variability in the measurements of resolved objects. As these issues might plague only one observation in a data set, they can lead to measurements easily confused for bona fide transients by automated methods.

One method to search for radio transients is by harnessing the power of citizen science. Citizen science projects hosted on the Zooniverse² have been highly successful in transient astrophysics and astronomy more generally, starting with the original Galaxy Zoo (Lintott et al. 2008). Since then the hundreds of public projects have received over 700 million classifications from 2.5 million users (taken from the website’s live tracker). In transient astronomy specifically, Wright et al. (2017)’s Supernova Hunters combined a neural network with human classifications to outperform either classifier alone and is still discovering supernovae, over 6 yr since launch.³ Similarly, Citizen ASAS-SN users have discovered >10 000 new variable sources that are not present in the existing star catalogues of the Southern hemisphere (Christy et al. 2022). The Zooniverse’s Talk feature (project-specific forums where users, moderators, and experts can discuss individual subjects, classifications, and so on) provides room for novel discovery space – classic Galaxy Zoo examples include the ‘Green Peas’, a class of compact galaxies with extremely high specific star formation rates (Cardamone et al. 2009) and ‘Hanny’s Voorwerp’, an extended region of gas ionized by the now-faded AGN of IC 2497 (Lintott et al. 2009). Similar finds from other projects include unusual variable stars and new classes of systematic noise in LIGO/Virgo detectors (Zevin et al. 2017; Christy et al. 2022). In this work we will describe the first citizen science project dedicated to radio transients. The aims of this project are to discover new transients, eliminate spurious false positives, and provide complementary analysis to other commensal search methods, as well as allow us to assess the viability of further citizen science work.

In Section 2 we detail the observations and processing prior to the Zooniverse project launch discussed in Section 3, the results of which are found in Section 4. We search for counterparts to radio variables in Section 5 before the discussion and conclusions found in Sections 6 and 7.

2 THUNDERKAT OBSERVATIONS AND PRE-PROCESSING

Our observations consist of a subset of ThunderKAT XRB images, based on which data sets were available at the time of research and contained more than a few epochs. The observations used in this work were taken between mid-2018 and late 2021. Generally, the observing strategy is determined by reports from X-ray facilities of activity from an XRB, which is then observed at weekly cadence by ThunderKAT in 15 min images. Each observing block consists of first scanning a primary calibrator, then bookending source observations with phase calibrator observations. A table of the 11 data sets used in this study, with the number of observations in each, is given in Table 1. The varying number of epochs between data sets is a direct result of the radio activity of the central XRB i.e. if a source is seen to fade below detection it is removed from the weekly scheduling block, while the central root mean square flux density (RMS) varies due to baseline coverage and presence of diffuse structures. The values presented are the median values of the RMS calculated across all images of a given data set, evaluated in the central eighth of an image. It is worth mentioning that the GX339–4 field has been observed every week since ThunderKAT observations began, in contrast to the one or two outburst cycles followed for all other data sets. It is important to remember that the commensal nature of this work constrains us to whatever observational cadence was used for monitoring the XRB.

ThunderKAT data are typically reduced using OxKAT (Heywood 2020), a semi-automatic set of scripts that perform calibration, flagging, and imaging of MeerKAT data. OxKAT uses several existing radio astronomy packages including casa (Mcmullin et al. 2007) for tasks such as gain and bandpass calibration, self-calibration, and flagging, Cubical (Kenyon et al. 2018) for further self-calibration procedures, tricolour for further flagging (Hugo et al. 2022), and wsclean (Offringa et al. 2014) for imaging. These steps are broken into first-generation calibration (1GC; direction-independent effects), flagging and self-calibration (2GC), with optional 3GC steps to account for direction-dependent effects. Some of the earlier observations were reduced prior to the release of OxKAT, however these still follow the same basic reduction of flagging using aoflagger (Offringa 2010), bandpass, phase calibration, and flux scaling in casa, and imaging with wsclean, DDFacet (Tasse et al. 2018), or casa. The commensal nature of this work means that, due to different science requirements and observational conditions for each data set, the resultant images are heterogeneous in their properties, although mostly homogeneous within a particular field.

2.1 Pre-processing and subject generation

Each set of images was processed using the Transients Pipeline (TraP Swinbank et al. 2015), first designed for detecting transients with LOFAR. Next is a brief description of how the pipeline works and some of the key parameters used. The TraP finds sources above a user-defined threshold in a set of astronomical images, creating light curves of each unique source and calculating statistics for each source. This is done by fitting a Gaussian component to each source in every epoch and associating it with those found at that position in all previous images, updating the data base as new observations are added. Most TraP parameters are kept at their default values.⁴ The detection_threshold, the signal-to-noise ratio (S/N) above which sources are detected, was fixed at 8 throughout as a trade-off between detecting false positives and missing genuine sources of interest. Once a source has been found, a Gaussian component is fit from its peak down to 3σ above the noise (analysis_threshold  = 3). The expiration, i.e. the number of force fits to a position where a source was found in a previous epoch, was always kept at greater than the total number of observations in a data set, meaning wherever a source had been found, a light curve with data points for all remaining time-steps was created. We are interested in unresolved, point sources and their variability, so we set all source fits to be fixed at the size of the PSF via force_beam  = True. For extended regions of emission, the change in size and position angle of the PSF between observations can lead to non-intrinsic variability measurements, as discussed in Section 3. To allow for deblending we set deblend_nthresh  = 10, which accounts for overlaps between nearby sources such as double-lobed radio galaxies. Finally, the extraction_radius_pix, describing how far ‘out’ in the image to search for sources, was always kept to approximately 1.5× the main lobe of the primary beam, in accordance with Sarbadhicary et al. (2021).

Of the TraP outputs, the most relevant are the light curves and two statistics computed based on the time series. For a light curve consisting of N data points of flux density F_i ± σ_i, observed at frequency ν, the two variability statistics are defined as

and

where s, |$\overline{F}$|⁠, and |$\overline{F_\nu }^{*}$| denote standard deviation, mean, and weighted average, respectively. Generally speaking we can use η and V to determine the statistical significance and the amplitude of variation, respectively. We expect η to correlate with average flux density; the brightest sources will have the smallest statistical uncertainties. Similarly, V and |$\overline{F}$| should be anticorrelated as we are only able to measure small variations for the brightest sources. We also note that flux calibration errors can lead to overinflated statistics for bright sources – e.g. systematic differences between epochs with small statistical uncertainties producing very large values of η, as equation (1) does not include systematic uncertainties in its calculation.

Once all observations have been processed by the TraP, we generate figures of the light curves and local sky around every runningcatalog source – that is, each of the 8874 unique entries to the data base. The local sky figure is a square arcminute image centred on the weighted mean right ascension and declination (RA, Dec.) of each source, using the image where the source was detected at the highest S/N. Finally, these figures are uploaded to the Zooniverse platform along with some basic metadata (RA, Dec., median flux density, and the time stamp of the highest S/N observation), creating the subjects for citizen scientists to classify.

3 CITIZEN SCIENCE PLATFORM

The Bursts from Space: MeerKAT (hereafter BfS:MKT) Zooniverse project⁵ launched on 7th December 2021 with classifications concluding by early March 2022. During this time 1038 volunteers classified our sources using the workflow seen in Fig. 1. Volunteers are given a tutorial to familiarize them with the data and describe the figures shown in the project, as well as accounting for common pitfalls due to figure processing – namely there are a few visualization issues that make classifications more difficult, discussed further in Section 6. There are also description pages for the project detailing the team, the telescope and the kinds of objects we are searching for so that citizen scientists can learn more about astronomy and the work we do.

Figure 1.

Classification workflow for BfS:MKT, showing the light curve and local sky figures for GX339–4.

Open in new tabDownload slide

Citizen scientists were given five classes to which they can assign a source, examples of which can be seen in Fig. 2. Stable sources are unresolved, point sources whose light curves are judged to be (within uncertainties) consistent with flat. The Extended Blob classification is intended to catch the resolved, extended sources, regardless of variability. We know that the changing size and angle of the PSF between epochs introduces non-intrinsic variability and so we instructed all volunteers to classify subjects they deem to be resolved as Extended Blobs – using the local sky figure, comparing the source size to that of the PSF in the lower left corner. The Artefact classification was implemented to account for any spurious, non-astrophysical sources that may be present in our images. Transient/Variable classifications are those we are searching for, which are point sources with variable light curves. Finally, if volunteers are uncertain if a subject fits into any of these classes – either due to visualization issues, their own interpretation or anything other reason – they can say they are Unsure. We included Unsure to assess the confidence of volunteers – if a subject does not clearly fit into one class this will be seen quantitatively (not just in e.g. the Talk board). Also, without an ‘unsure’ option, volunteers may have settled for classifying as either stable or transient, leading to an underprediction or overprediction of interesting sources. There is also a Field Guide (accessible on the main project webpage) with examples to demonstrate the type of subjects intended for each class, as well as some help text describing the rough thought process behind each source. If volunteers feel a subject is particularly interesting, or they have questions on a particular source, they can create individual threads on the dedicated Talk forum for the project, where (citizen and project) scientists can discuss.

Figure 2.

Examples of the four observational classes within our workflow on the Zooniverse, showing both the light curve and an image of each. From top to bottom these are: Stable – no variation in the light curve given the error bars and a point source; Extended – variations caused by changes to the PSF and a source that is larger than the beam (lower left); Transient – a clear variable light curve for a source the same size as the PSF; Artefact – a spuriously transient light curve and a faint source on the outskirts of a very bright object, with non-Gaussian noise structure. The final class, Unsure, by definition has no archetypal characteristic so we show no figure here.

Open in new tabDownload slide

We require 10 volunteers to classify a subject before we consider it classified, resulting in a total of 88 740 classifications. These were classified over 90 d, or an average of 1 submission per 1–2 min. The histogram of all classifications for the project to date can be seen in Fig. 3, showing the expected steep power law distribution of votes (Spiers et al. 2019), as well as a number of ‘super users’ who have classified thousands of sources each. The median, mean, and standard deviation of user classifications are 4.5, 86, and 490, respectively. The chosen retirement value of 10 is enough to mitigate outlier votes, but not so high that it would take many months for a single subject to be fully classified. We note here that if further iterations of the project gain as much traction as this first batch of subjects we will be able to determine what the optimal trade-off could be.

Figure 3.

Histogram of the number of classifications made per user on the BfS:MKT project. As might be expected, many users perform only a few classifications, while a few volunteers dedicate thousands of votes to the project.

Open in new tabDownload slide

Simple aggregation is performed for this one question workflow, using the standard aggregation scripts,⁶ where we take the Boolean values for each classification and sum over every vote for each subject. This gives us 10 votes for a given source, from which we can calculate fractional classifications and determine how many subjects are deemed to be transient/variable by some number of citizen scientists (TF  = transient fraction).

We set a threshold of 4/10 volunteers classifying a subject as transient/variable, reducing our sample from 8874 to 381 sources. This 0.4 threshold was chosen as a trade-off between having many sources to vet and missing some low vote fraction variables. These 381 subjects were visually inspected by project scientists to confirm or reject each source as one where both volunteers and experts agree. This final vetting reduces our number of transient candidates to 168 sources (i.e. a true positive rate of 168/381 = 0.44). Reasons for disagreement between citizen and project scientists are numerous and include how much volunteers use error bars, PSFs, and other parts of the figures, as well as the existence of systematic variability that is present in earlier MeerKAT observations (see e.g. the appendix of D22) that would only be noticeable to experts who have compared many light curves in a given field. We will discuss what information we get from the false positives in Section 6.

4 RESULTS

Using citizen scientists to inspect a wealth of data from the MeerKAT telescope produces 168 variable and transient radio sources which have also been vetted by project scientists. Of these 168 variables and transients found by volunteers, 142 are not detailed in previous commensal or XRB work. This constitutes one of the largest samples of radio variable candidates to date, and their positions are listed in Table A1. This table provides the TF, median 1.28 GHz flux density, and η and V statistics described in Section 2.1 for each of our sources, as well as the date on which they were detected at highest S/N. We strongly encourage follow-up of sources of interest to the community.

In order to characterize these sources, we can examine the variability plane defined by η and V (equations 1 and 2) seen in Fig. 4. The structure seen in this figure consists of many sources at low significance and amplitude, some spurious artefacts at high V and generally bright objects detected at high significance (η) but low variability amplitudes. Systematic differences between data sets due to their heterogeneous sampling and imaging are also present (e.g. the sources not voted as transients at log(V) ∼−0.6). The most important thing to note is that there are many variable objects whose parameters are ‘normal’ i.e. non-anomalous. This means that, were we instead to take outliers above some Nσ threshold in η and V, we could have missed sources that, upon individual inspection, appear variable or transient. So citizen scientists can find interesting variable radio sources that could have been missed by other techniques or without detailed analysis (e.g. D22, R22 and Chastain et al. in preparation). The previously known transient sources (i.e. the circled sources of Fig. 4) with high η and V are almost all recovered, as are many of those found in previous studies (see Sections 4.1 and 4.2 for more detail), meaning that volunteers are able to recover or discover interesting sources across a range in statistical parameters. So our volunteers have been able to analyse large data volumes in just a few months and produce the largest sample to-date of variables from a radio telescope.

Figure 4.

Variability plane for the 168 sources found by citizen scientists to be variable, along with those they find to not be so in grey. The colourbar denotes the fraction of classifications as a transient/variable source, while marginal distributions of η and V can also be seen. Known variable sources (e.g. XRBs) in our fields are circled. Imaging artefacts appear at low η and high V, while flux calibration uncertainties can produce high η, low V sources (due to lack of systematic uncertainty in equation (1). Most known transients are found by citizen scientists, while many new sources are identified and show a wide spread of values in this parameter space.

Open in new tabDownload slide

We can use these discovered (142) and recovered (26) variables, along with the other known transients (19) in our fields, to estimate that at least 2.1 per cent of radio sources seen at 1.28 GHz with our approximately weekly observations are transient or variable. This rate is in line with previous work on radio transients; see Ofek et al. (2011) for a review, as well as D22 and Sarbadhicary et al. (2021). We also note that these different studies, including this work, probe slightly different time-scales, though all of them ‘long’ (≥ a week). This means that our volunteers and project scientists can find comparable transient rates as previous work. This 2.1 per cent is a lower limit on the number of variables and we can estimate what fraction of variables/transient are missed by assessing which known transients are not recovered. The fraction of previously known variables not recovered is 19/45, implying the true amount of variables in our fields could be as high as ∼5 per cent. However, this estimate is only valid under the assumption that new variables are recovered at the same rate as the known transients in our field. In Section 6 we discuss the selection effects that are evident in the different kinds of light curves selected by volunteers when compared to the known transients in our field.

The majority of these 142 variables that our volunteers discover show long-term variability, with light curves showing variations over weeks to months. One cause of variability for extragalactic sources is scintillation through the Galactic ISM. All of our observations are within latitude |b| ≲ 10°of the Galactic plane (see Table 1) so this may be a large contributor to the variability seen in our sample. We can use the model from Hancock et al. (2019) to characterize the effect of refractive interstellar scintillation (RISS) for our set of variables, using an underlying H α map from Finkbeiner (2003) to quantify the electron scattering along a given line of sight through the Galaxy. This model predicts the level of variability for a given radio frequency. We can directly compare this predicted maximum RISS-induced variability to our measured V values, as seen in Fig. 5. For the majority (131) of our new variables we see that the predicted modulation is equal to or greater than that measured by TraP i.e. the observed variability is consistent with (though not exclusively explained by) the scintillation of light from AGN. By contrast, the known Galactic XRBs and their jets – whose variability is caused by shocks and particle acceleration – show variability that is much greater than what would be expected due to refractive scintillation. The predicted RISS variability is in some cases much higher than our observed V values. This is likely caused by the heterogeneous sampling of our fields, the coarseness of the H α model grid, and/or the assumptions of the model. Finally, we can calculate a weighted average time-scale of variation for our sources, which we find to be t₀ = 8 ± 4 months, where the weights used are propagated through from the uncertainties in the underlying map and the quoted uncertainty is the standard deviation. This range of time-scales of variation at 1.28 GHz matches well with the length of typical observations for our XRB fields. Both these matching amplitudes and time-scales of variability provide evidence for the majority of our transients being scintillating AGN or other point-like extragalactic radio sources.

Figure 5.

The ratio between observed V and the variability predicted by Hancock et al. (2019)’s model for RISS, for our sample of transients and variables. Most variability can be explained by this model of a scintillating, extragalactic source, apart from the known XRBs and jetted systems.

Open in new tabDownload slide

4.1 Comparison to target sources

Of all the 8874 sources classified in this project, there are 45 known variable/transient objects published in the literature, including the 11 XRBs listed in Table 1. Of the 11 XRBs, nine are classified as transient by citizen scientists. The only two that are missed are Swift J1858 and SAX J1808, whose light curves can be seen along with that of EXO 1846 for contrast in Fig. 6 – exactly as citizen scientists would have seen them. We can explain why Swift J1858 was not classified as transient due to a combination of there being only one significantly bright data point, as well as the figure generation creating an overly large legend. Similarly, SAX J1808’s misclassification can be understood as stemming from uncertainty surrounding so few data points (especially compared to other data sets). Indeed, SAX J1808 received one ‘Unsure’ vote which, had it instead been for ‘Transient/Variable’, would have pushed this subject above our classification threshold.

Figure 6.

The TraP light curves of EXO 1846–031, Swift J1858.6–0814, and SAX J1808.4–3658 (upper left, upper right, and lower, respectively). The foremost was classified as transient by citizen scientists, while the others were not. We believe these to be caused by bad figure generation and lack of data points, respectively.

Open in new tabDownload slide

In addition to these central 11 sources, several of the XRBs also display discrete, resolved jet components (MAXI J1820+070, MAXI J1348–630, MAXI J1848–015, and 4U1543–47; see Table 1) that are counted as unique sources by the TraP as they become resolved and move away from the core of the jet. In total there are six moving, transient jet components that are detected by the TraP as unique, of which three are classified as variable by citizen scientists. Given that neither the software pipeline nor project classifications were designed to pick up moving point sources, it is interesting to note that there is room for unexpected discoveries even in such simple workflows.

4.2 Previous commensal studies

One of the largest works on variable sources in MeerKAT images to date is that of D22. In that work, 21 new long-term variables (LTVs) are detailed, along with GX339–4, the first MeerKAT transient MKT J170456.2–482100 and the known mode-changing pulsar PSR J1703–4851 (both described in Driessen et al. 2020). Of these 24, excluding the XRB discussed earlier, MKT J170456.2–482100, PSR J1703–4851, and 10 of the LTVs are missed by classifications i.e. only 11 are marked as transient/variable by volunteers. One reason for this is that in D22, each light curve is binned by a factor of 10 (i.e. every 10 data points are represented by their uncertainty-weighted mean) to make long-term variability more apparent, something not done here. Furthermore, here we made no use of a deep, stacked observation in order to follow sources through every epoch (see Section 2.1’s details on the expiration parameter), something that will result in different light curve shape and in one case resulted in a source not being detected at all in our set of images. In general, the LTVs not identified by volunteers were generally fainter and have less smooth light curve evolution. This may give some insight into a bias of our method; it is easier to find light curves with long-term, clear evolution as opposed to more stochastic variability.

Similarly, when comparing to the work of R22 we see that the three variables found therein – NVSS sources J181849+062843, J181752+064638, and J182029+063419 – were recovered by citizen scientists. However, the three sources found in the ‘transient hunt’ (where therein transient is defined as sources not detected in a deep observation of the field, see their section 3.1) were not identified by volunteers. As above we note that the variables recovered have smooth light curve evolution over long time-scales, while those not identified as transient are much fainter, with larger uncertainties and figures that are more ‘cluttered’.

Finally, the M dwarf SCR 1746–3214 is a radio transient that exhibits flares, serendipitously seen in early ThunderKAT data (Andersson et al. 2022). In order to assess what light curves citizen scientists were most comfortable with classifying, we provided two light curves of this source – one with more data points and an additional detection and one with only two detections and two upper limits (see fig. 3 in the above article). Interestingly, when provided with the shorter light curve of the source, the transient source was mostly classified as ‘Unsure’ or ‘Stable’. However, when given the full light curve volunteers correctly identify the subject as transient. We can perhaps use this to infer that citizen scientists were least unsure when classifying sources with more data points and less reliant on upper limits, as is the case for many of the new variables found in this study.

5 COUNTERPARTS AND ASSOCIATIONS

We use the MeerLICHT optical telescope to gain physical insight into possible source classes of our radio variables and transients. We use MeerLICHT for this due to its position in South Africa and its mission to follow the radio observations of MeerKAT, resulting in highly complementary spatial and temporal coverage of our radio sources (if observed at night). A typical observing schedule consists of 1 min exposures of a given field, alternating between q band (440–720 nm) and each of five Sloan bands u, g, r, i, and z. While MeerLICHT operates in these 6 bands, here we use only the q band due to its highest sampling rate and its broad wavelength coverage. We crossmatch with the MeerLICHT data base at a 2 arcsec search radius – large enough to partially account for MeerKAT astrometric uncertainties (see D22) and proper motion (e.g. of nearby stars) but not so large as to include many false matches – using the uncertainty-weighted mean position of each of our candidates as returned by the TraP. Running this crossmatch returns 25 counterparts in the MeerLICHT data base and 143 3σ upper limits. This low rate of optical counterparts is not surprising as all of our XRB fields are within 10 deg of the Galactic plane and so many optical counterparts may be heavily extinct.

The radio-optical plane for our variables can be seen in Fig. 7, with comparison source types from Stewart et al. (2018).⁷ The first thing to note is that the majority of these cross-matches exist in the region where extragalactic sources have been detected – either quasars or GRBs. If most of our candidates are extragalactic in nature this agrees with previous studies, who find the vast majority of radio variables are extragalactic (Thyagarajan et al. 2011; Sarbadhicary et al. 2021). There are a few sources overlapping the ‘stellar’ region of the parameter space – these include the transients already reported by ThunderKAT (in Driessen et al. 2020 and Andersson et al. 2022). It is important to note that by comparing to archival data in this way does not leave room for unknown astrophysical classes, however this still gives an indication of the overall distribution of multiwavelength counterparts.

Figure 7.

The mean optical and radio flux densities of our sample of radio variables, atop an underlying distribution of astrophysical classes (Stewart et al. 2018). The black crosses denote counterparts within the MeerLICHT data base, while the grey triangles are upper limits. The diagonal lines denote a constant ratio between radio and optical flux density, while the A_R marker indicates the horizontal displacement caused by 5 mag of optical extinction. The majority of our radio sources are likely extragalactic as they overlap in parameter space with quasars and GRBs.

Open in new tabDownload slide

5.1 Highlights

We also search for counterparts at other wavelengths and with pre-existing classifications, with the aid of pre-existing code available in Driessen (2021). Crossmatching uses the astroquery package to search the SIMBAD⁸ data base (Wenger et al. 2000) and several catalogues within Vizier⁹ (Ochsenbein, Bauer & Marcout 2000). We again search at a 2 arcsec radius to our radio sources. We searched for X-ray and gamma-ray counterparts to our radio variables. To do this we crossmatched with catalogues from the Fermi (Schinzel et al. 2014), Chandra, (Evans et al. 2010), Swift (Evans et al. 2020), and XMM–Newton (Traulsen et al. 2020) facilities. There were no counterparts for any of our sample, aside from the known transients of interest (i.e. the XRBs). Next, we make note of a few interesting objects returned in our search with IR or radio detections, or otherwise known counterparts.

5.1.1 OH Maser – BfS 80

We detected transient emission from a source in the EXO 1846 field which volunteers confidently classified as transient. Cross-matching reveals that the source in question is a known maser star OH 30.1–0.7 otherwise known as V1362 Aql, with radio observations stretching back almost 50 yr (Evans et al. 1976). This asymptotic giant branch (AGB) star is heavily dust-obscured at optical wavelengths, but very bright at IR wavelengths – W1 = 6.9 ± 0.1 mag, or 279 Jy at 25 μm (Cutri et al. 2013; Gonidakis et al. 2014). A comparison of a MeerKAT radio detection (contours) and Spitzer Glimpse imaging can be seen in Fig. 8, alongside its light curve. AGB stars are post-main-sequence systems, whose low surface temperatures, radii of several hundred times the solar radius, and stellar pulsations give rise to strong winds, expelling as much as ∼10⁻⁵ M_⊙ yr⁻¹ (Höfner & Olofsson 2018). These winds create an oxygen-rich, dusty circumstellar environment that generate masers at 1612 and 1667 MHz as infrared photons pump OH molecules formed through photodissociation of water by interstellar radiation.

Figure 8.

Upper: Light curve of OH Maser 30.1–0.7 picked up by volunteers. There were four epochs in which this source was not detected prior to this observed variability. Lower: Overlay of MeerKAT radio contours over Spitzer imaging from the GLIMPSE survey (Benjamin et al. 2003). Contours are spaced linearly in 0.5 mJy increments from −0.5 to 3.5 mJy.

Open in new tabDownload slide

The cause of the variability of OH 30.1–0.7 is not clear. Perhaps the OH maser emission is varying, due to stellar pulsations. The derived stellar period from the General Catalog of Variable Stars is ∼1730 d (Samus’ et al. 2017), with more recent estimates from the WISE W1 and W2 bands at 1950 ± 50 and 1520 ± 20 d, respectively (Groenewegen 2022). These time-scales are much longer than the variability seen in the radio light curve here (of order a few months), so it is not clear if stellar pulsations are responsible for variable maser emission (averaged across the entire L band from Jy to mJy levels). There could also be inhomogeneities at the site of the maser emission. The second cause for variability could be due to binary interactions – the system likely has a companion, as inferred from ALMA CO data in Decin et al. (2019). This, combined with the lack of any optical counterpart implies the system could be a dusty (D-type) symbiotic binary – all D-type symbiotics host Mira stars (Whitelock 1987). These binaries consist of a windy red giant and a smaller companion on to which material is shed (Allen 1984). Radio emission has been seen from symbiotic binaries (Seaquist, Taylor & Button 1984) and is known to vary in some sources, typically interpreted as optically thick emission from an inhomogeneous region in the red giant’s wind that is ionized by its companion (Seaquist 1988). The nature of OH 30.1–0.7’s companion is as yet unknown due to the heavy extinction in this region and without evidence of high ionization (e.g. He ii or [O iii]) we cannot claim that this is a symbiotic system. There is overlap between maser systems and Mira-type symbiotic binaries (Seaquist, Ivison & Hall 1995) and so the observed variability could be due to a combination of emission mechanisms discussed, or perhaps something else entirely. One thing to note is that there were four observations of this field prior to the initial data point of the light curve – i.e. there are four non-detections before the system reached ∼3 mJy as seen. Future radio observations could help determine what the nature of the variability is (including the initial non-detections), perhaps combined with spectroscopic searches for the nature of the companion (e.g. if it is a white dwarf, or for evidence for ionization).

5.1.2 Pulsar – BfS 20

Pulsar PSR B1845-01 (J1848–0123), whose light curve can be found in Fig. 9, was seen in the field surrounding XRB MAXI J1848. In our MeerKAT observations we measure a mean flux density of ∼15 mJy, in close agreement with that measured by the Parkes radio telescope recently (15.2 ± 0.9 mJy; Johnston & Kerr 2018). However, reports of this pulsar’s flux density have been as low as 8.9 ± 0.9 mJy at 1400 MHz (Hobbs et al. 2004), indicating long-term changes in the received brightness from the source. The pulsar is brighter at lower frequencies (e.g. measured to be 79 ± 6 mJy at 408 MHz by Lorimer et al. 1995), with a spectral index α – defined at frequency ν as flux density F ∝ ν^α - of −1.3 ± 0.3 (McEwen et al. 2020).

Figure 9.

Light curve of PSR B1845–01 with MeerKAT. 7 of 10 volunteers voted for this as a transient/variable source. The low amplitude variability observed over a time-scale of 100s of days is consistent with RISS.

Open in new tabDownload slide

The observed variability is consistent with RISS (see Section 4), with a predicted modulation greater than the measured V ∼ 0.02 and the observed time-scale of variation matching the estimated 11 ± 2 months. Furthermore, this pulsar has a relatively high dispersion measure (DM) of 159.1 ± 0.2 pc cm⁻³ (McEwen et al. 2020; compared to typical DMs of order 10s of pc cm⁻³ for pulsars with |b| > 25°; Manchester et al. 2005), which is known to be linked to long time-scale, low amplitude refractive scintillation (Stinebring et al. 1990, 2000). Similarly, the pulsar’s location in the Galactic plane is in keeping with its radio emission traversing a large free electron content, hence the high DM and clear scintillation. We note that B1845-01’s spin period of ∼0.65943s is much shorter than that of our observations (typically 15 min epochs consisting of 8 s correlator sampling) so this cannot be contributing to the observed variability.

This pulsar adds to the diverse range of behaviours seen in pulsars spotted by MeerKAT in imaging data. Similar examples include the mode changing pulsar observed in the GX339–4 field (see D22) and, in the most extreme case, one of the slowest pulsars discovered (Caleb et al. 2022).

5.1.3 VLASS1 J181955.28+074418.7 – BfS 146

Radio source VLASS1 J181955.28+074418.7 is an object spotted by our volunteers not only in our aggregated results (scoring 7/10 votes as a transient/variable) but also in our Talk board.¹⁰ This source is in the field of MAXI J1820 but outside the 0.5 deg radius set by R22. Its light curve can be seen to vary very smoothly in Fig. 10, where it is worth noting that the non-detections are all due to drops in data quality in those images, which should be filtered out in future work. This variation is not similar to any nearby source in our data, nor is it correlated with PSF size or shape. There is a counterpart to this source in the Very Large Array’s Sky Survey (VLASS QuickLook Epoch 1; Lacy et al. 2020; Gordon et al. 2021), with a 3 GHz flux density of 1.9 ± 0.3 mJy. There are no counterparts to this source in any higher energy band, despite this field being our furthest from the Galactic plane (b ∼ 10°). This kind of source is very typical of our sample – scant extra data but an intriguing radio light curve. Further information at other wavelengths would help determine source type, as would simultaneous radio data – e.g. to determine a precise spectral index α and see if this points towards the source being an AGN or a pulsar. Our RISS analysis (see Section 4) produces a ratio of 1.17 between the observed V and predicted scintillation amplitude. This could indicate the source has some intrinsic variation; however, the scintillation parameters were not exhaustively tested, so this ratio could be explained by an incorrectly assumed distance to, or relative motion between, the model scattering screen and the observer.

Figure 10.

Light curve of source VLASS1 J181955.28+074418.7, showing smooth, near sinusoidal variations. The epochs denoted upper limits are in error and due to not filtering out low-quality images in our pre-processing.

Open in new tabDownload slide

6 DISCUSSION

Using citizen science to discover transients is a fruitful endeavour, as shown by both the uptake of our Zooniverse project and the novel results produced. For our volunteers, the project provides new experiences with a branch of observations perhaps less familiar to their perception of astronomy (e.g. compared to galaxy morphologies or solar physics). The uptake of our project demonstrates a clear appetite for further development of citizen science for radio transients, as evidenced by >1000 citizen scientists applying their time to our science case. The time taken to individually check several years worth of data is much greater than the 3 months taken by our volunteers, so citizen scientists help project scientists analyse their observations more efficiently. In this study we have been able to recover or discover a broad range of astrophysical transients that occur over several orders of magnitude in time-scale, so the science results discussed here also give merit to developing BfS:MKT further. These results include flaring from stellar systems, the discrete and compact jets of XRBs, maser emission, pulsars, the long-term variation in likely AGN and potentially new source classes. Given these successes, we plan to launch a second wave of classifications to the Zooniverse site shortly, using different data from surveys on MeerKAT to explore the range of parameter space for both citizen and project scientists. Of course this method is not restricted to MeerKAT data and any such set of (radio) images and light curves could be analysed in this way.

However, citizen science is not without its challenges and we can understand some of these by comparing our findings to those of previous work from the ThunderKAT team – namely D22 and R22, as well as the published XRB work (see Table 1). The first thing to note is that not all of the transients from D22 and R22 were recovered here. Some of this can be explained by the difference in pre-processing of images and TraP light curves (discussed in 4.2). However, several of the transients missed are due to how we have produced subject images for the Zooniverse project. For example, the two XRBs not classified as transients can be explained due to the lack of clarity in Fig. 6 (see Section 4.1). Similarly, the image of source MKT J182015.5+071455 (R22’s source 2) provided to volunteers was automatically scaled to the brightest pixel in the image, not be the central source of interest. As such, the source appears not to be present and was classified mostly as an artefact. We will use these issues to improve our procedure for future batches of data e.g. by manually setting the pixel scale in images, altering legend sizes etc.

Aside from these pre-processing issues, there are still some known transients not recovered, including the pulsars in D22 and R22. When we compare the light curves of these pulsars to variables that are recovered we see a clear trend – recovered transients show long-term, smooth variations, and are typically brighter, with smaller uncertainties, resulting in clearer patterns on display to scientists. By contrast, the faint, transient pulsar light curves show very ‘noisy’ light curves, with less eye-catching patterns, despite being precisely the kinds of transients we want to discover. R22’s pulsar received three votes as a transient/variable and one unsure classification – this could be due to the heterogeneous nature of our classifiers, or it could reflect the uncertainty surrounding a less clear pattern in a light curve. As mentioned in Section 4.2, we provided two light curves of Andersson et al. (2022)’s M dwarf to volunteers, one with force fit measurements (more data points and an additional detection) and one without. The light curve with more data points passed our threshold of 0.4, while the latter did not, perhaps indicating that classifiers are more comfortable with longer light curves. These unrecalled transients give insight into the limitations of this data set: our sample of transients and variables is likely biased towards brighter, slower evolving objects that occur in our most sampled fields. We hope to alleviate some of this bias in future Zooniverse runs by emphasizing the use of non-detections and by encouraging volunteers to label things as transient. We could also implement the transients and variables discussed in this work into the Field Guide for volunteers as more examples of the types of sources for which we are looking.

We can quantify how scalable our method is by comparing our observations to the MeerKAT MIGHTEE survey. MIGHTEE’s observations produce ∼6000 sources per square degree on the sky (Heywood et al. 2022), for which our TraP processing would produce a light curve. For comparison, in this work we produce ∼1000 subjects from the ThunderKAT fields that are devoid of large diffuse structures or exceedingly bright sources. If we assume that only the same 1000 volunteers contribute to all future data releases as with this study, at the same classification rate, then it would take approximately 60 days to receive 10 classifications on every source at MIGHTEE’s sensitivity (1σ RMS noise ∼ 1 μJy), per square degree. Multiplying this across the sky coverage of just the ThunderKAT fields used in this study (∼1.5 square degrees over 11 fields) results in a required volunteer classification time of 990 days. This is far greater than the time taken for the observations (e.g. of order tens of 8-h epochs) and which is needed to image the data and process them to form light curves. If the volunteer results were only analysed when all classifications were finished then this would also be too late for real-time follow-up of transients. To bring the classification time to that of this work (90 days) would require an order of magnitude more volunteers, which is achievable for Zooniverse projects, particularly when disseminated widely. For example, Galaxy Zoo variants receive many 10s of thousands of volunteers, while Gravity Spy has had over 30 000 participants. However, for any considerable survey area the overall time required would again balloon to far larger time-scales than reasonable, particularly if we process data sets and classifications in batches once observations are complete. Finally, observations at the 8s integration time for MeerKAT produce far fewer sources to classify compared to deeper images, however the imaging and processing time prior to volunteer classification increases hugely so the overall time-scale remains long. So this kind of analysis will not scale easily to the most sensitive observations of MeerKAT fields (at 1.28 GHz), let alone those expected from the SKA, ngVLA (Hallinan et al. 2019) or DSA-2000 (Selina et al. 2018).

One way to alleviate this data deluge might be to develop machine learning methods to remove ‘bogus’ and ‘boring’ sources in favour of the rarer variables and transients for which we are searching, as has been done for e.g. supernovae and galaxy morphologies (Wright et al. 2017; Walmsley et al. 2022b). Active learning, where humans feed back to machine learning techniques in order to prioritize sources of interest and optimize precious human attention, has been shown to uncover unique light curves (Ishida et al. 2021) and radio morphologies (Lochner et al. 2023), and are able to optimize the volunteer classification of optical galaxies (Walmsley et al. 2022a). We are currently applying Lochner & Bassett (2021)’s Astronomaly active learning framework to the data presented in this work in the hopes that we can use the combined power of human classifiers and machine processing to extract the most science from our wealth of data (Andersson et al. in preparation).

7 CONCLUSIONS

In this work we have presented the first citizen science project for finding transients in image plane radio surveys. The uptake of the project was very strong, with >1000 volunteers taking part, demonstrating a healthy appetite for further Zooniverse data releases. We were also able to use the known transients in our fields to understand some reasons why interesting sources may be missed and will fold this learning through to future iterations of the project. Citizen and project scientists uncovered a large sample of interesting transient and variable sources, some of which we may have not uncovered were it not for our volunteers’ dedication. We provide the full catalogue of 168 radio transients and variables, the largest catalogue of candidate radio variables to date. This includes links to images and light curves, and we encourage others to follow up these sources, and additional future catalogues that this project will deliver. We used archival multiwavelength data, including the MeerLICHT telescope, to help classify the systems found. The sources found span a broad range of physical phenomena including pulsars, radio-loud stars, XRBs, and a large set of AGN – likely varying due to scintillation. These results demonstrate the wealth of science possible with new radio facilities. Finally, we hope to use volunteer classifications to develop anomaly detection algorithms, with an eye towards current and future surveys such as the SKA.

ACKNOWLEDGEMENTS

The authors would like to thank our reviewer Matthew Bailes for his helpful comments that have improved this study. We also thank Patricia Whitelock for her helpful discussion points around OH masers, Mira variables, and symbiotic variables. AA acknowledges the support given by the Science and Technology Facilities Council through an STFC studentship. CJL acknowledges support from the Alfred P. Sloan Foundation. FC acknowledges support from the Royal Society through the Newton International Fellowship programme (NIF/R1/211296). LND acknowledges support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 694745). PJG acknowledges support from NRF SARChI Grant 111692. MeerKAT is operated by the South African Radio Astronomy Observatory (SARAO), which is a facility of the National Research Foundation, an agency of the Department of Science and Innovation.

This publication uses data generated via the Zooniverse.org platform, development of which is funded by generous support, including a Global Impact Award from Google, and by a grant from the Alfred P. Sloan Foundation.

We thank the SARAO staff involved in obtaining the MeerKAT observations. We acknowledge the use of the ilifu cloud computing facility – www.ilifu.ac.za, a partnership between the University of Cape Town, the University of the Western Cape, the University of Stellenbosch, Sol Plaatje University, the Cape Peninsula University of Technology, and the South African Radio Astronomy Observatory. The ilifu facility is supported by contributions from the Inter-University Institute for Data Intensive Astronomy (IDIA – a partnership between the University of Cape Town, the University of Pretoria, the University of the Western Cape, and the South African Radio astronomy Observatory), the Computational Biology division at UCT and the Data Intensive Research Initiative of South Africa (DIRISA). The MeerLICHT consortium is a partnership between Radboud University, the University of Cape Town, the Netherlands Organisation for Scientific Research (NWO), the South African Astronomical Observatory (SAAO), the University of Oxford, the University of Manchester and the University of Amsterdam, in association with and, partly supported by, the South African Radio Astronomy Observatory (SARAO), the European Research Council and the Netherlands Research School for Astronomy (NOVA).

This research has made use of the SIMBAD data base, operated at CDS, Strasbourg, France (Wenger et al. 2000), as well as the ATNF Pulsar catalogue (Manchester et al. 2005) available at https://www.atnf.csiro.au/research/pulsar/psrcat/.

This research has used the astropy package, an Astronomy-based, community-developed Python package (Astropy Collaboration 2013, 2018, 2022).

DATA AVAILABILITY

ThunderKAT raw data are available on the SARAO archive (https://archive.sarao.ac.za/). The data for each light curve mentioned in Table A1 and throughout the text will be available at https://github.com/AnderssonAstro/BfS-MKT-Analysis and we encourage readers interested in particular sources to investigate them in further detail, giving credit to this work. The data used to produce our figures will also be in said repository, excluding personal data from volunteers, which can be shared on reasonable request to the authors.

Footnotes

References

APPENDIX A: TABLE OF RESULTS