Bursts from Space: MeerKAT - The first citizen science project dedicated to commensal radio transients

The newest generation of radio telescopes are able to survey large areas with high sensitivity and cadence, producing data volumes that require new methods to better understand the transient sky. Here we describe the results from the first citizen science project dedicated to commensal radio transients, using data from the MeerKAT telescope with weekly cadence. Bursts from Space: MeerKAT was launched late in 2021 and received ~89000 classifications from over 1000 volunteers in 3 months. Our volunteers discovered 142 new variable sources which, along with the known transients in our fields, allowed us to estimate that at least 2.1 per cent of radio sources are varying at 1.28 GHz at the sampled cadence and sensitivity, in line with previous work. We provide the full catalogue of these sources, the largest of candidate radio variables to date. Transient sources found with archival counterparts include a pulsar (B1845-01) and an OH maser star (OH 30.1-0.7), in addition to the recovery of known stellar flares and X-ray binary jets in our observations. Data from the MeerLICHT optical telescope, along with estimates of long time-scale variability induced by scintillation, imply that the majority of the new variables are active galactic nuclei. This tells us that citizen scientists can discover phenomena varying on time-scales from weeks to several years. The success both in terms of volunteer engagement and scientific merit warrants the continued development of the project, whilst we use the classifications from volunteers to develop machine learning techniques for finding transients.


INTRODUCTION
The latest generation of radio telescopes provide 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), whilst 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), whilst 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) allow 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 aretransient 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).Whilst 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 can cause short time-scale, large amplitude variations in observations of pulsars, whilst 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 maximise the science yield of our observations.
ThunderKAT 1 (Fender et al. 2016) is a large survey project dedicated to monitoring radio transients with MeerKAT.The Thun-derKAT team regularly observes known transients such as X-ray binaries, cataclysmic variables and gamma-ray bursts (XRBs, CVs and GRBs respectively).The large field of view of MeerKAT (∼ 1 square degree 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), Driessen et al. (2022a) and Andersson et al. (2022) describe the first commensal transients found with MeerKAT, 1 The Hunt for Dynamic and Explosive Radio transients 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) make use of the best-sampled Thun-derKAT 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. (submitted), wherein there are many newly described variables, of which most are likely scintillating AGN.The ThunderKAT survey also makes use of Meer-LICHT (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.Firstly, the volume of data to analyse is far greater than any one person can achieve by eye on reasonable time-scales.In the 4 years 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 minute 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) andMHONGOOSE (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 ∼ 1Jy (Heywood et al. 2022).These data overload issues are only exacerbated when imaging on shorter time-scales, including down to the 8s integration time of MeerKAT, as is currently being tested within ThunderKAT (e.g. Caleb et al. 2022;Chastain et al. submitted;Fĳma et al. in prep).
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 dataset, 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 Zooniverse2 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 six years since launch3 .Similarly, Citizen ASAS-SN users have discovered >10,000 new variable sources that are not present in the existing star catalogs 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 star formation rates (Cardamone et al. 2009) and 'Hanny's Voorwerp', an extended region of gas ionised 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 (Christy et al. 2022;Zevin et al. 2017).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 allowing 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.

THUNDERKAT OBSERVATIONS AND PRE-PROCESSING
Our observations consist of a subset of ThunderKAT XRB images, based on which datasets 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 minute images.Each observing block consists of first scanning a primary calibrator, then bookending source observations with phase calibrator observations.A table of the 11 datasets used in this study, with the number of observations in each, is given in Table 1.The varying number of epochs between datasets 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, whilst 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 dataset, evaluated in the central 8th 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 datasets.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 is typically reduced using O KAT (Heywood 2020), a semi-automatic set of scripts that perform calibration, flagging and imaging of MeerKAT data.O KAT makes use of several existing radio astronomy packages including (Mcmullin et al. 2007) for tasks such as gain and bandpass calibration, self calibration and flagging, C (Kenyon et al. 2018) for further self calibration procedures, for further flagging (Hugo et al. 2022) and (Offringa et al. 2014) for imaging.These steps are broken into 1st 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 O KAT, however these still follow the same basic reduction of flagging using -(Offringa 2010), bandpass, phase calibration and flux scaling in , and imaging with , DDF (Tasse et al. 2018) or .The commensal nature of this work means that, due to different science requirements and observational conditions for each dataset, the resultant images are heterogeneous in their properties, although mostly homogeneous within a particular field.

Pre-processing and subject generation
Each set of images was processed using the Transients Pipeline (T P Swinbank et al. 2015), first designed for detecting transients with LOFAR.Below is a brief description of how the pipeline works and some of the key parameters used.The T P 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 database as new observations are added.Most T P parameters are kept at their default values4 .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 dataset, meaning wherever a source had been found, a light curve with datapoints 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 = T .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 T P outputs, the most relevant are the light curves and two statistics computed based on the time series.For a light curve consisting of  data points of flux density   ±   , observed at frequency , the two variability statistics are defined as and where ,  and   * denote standard deviation, mean and weighted average respectively.Generally speaking we can use  and  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,  and  should be anti-correlated 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 T P, we generate figures of the light curves and local sky around every runningcatalog source -that is, each of the 8874 unique entries to the database.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.

CITIZEN SCIENCE PLATFORM
The Bursts from Space: MeerKAT (hereafter BfS:MKT) Zooniverse project5 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 Figure 1.Volunteers are given a tutorial to familiarise 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 visualisation 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.Citizen scientists were given five classes to which they can assign a source, examples of which can be seen in Figure 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 -making use of 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 classeseither due to visualisation 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 under-or over-prediction 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.We require 10 volunteers to classify a subject before we consider it classified, resulting in a total of 88740 classifications.These were classified over 90 days, or an average of 1 submission per 1-2 minutes.The histogram of all classifications for the project to date can be seen in Figure 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.
Simple aggregation is performed for this one question workflow, using the standard aggregation scripts6 , 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 make use of 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. (2019) to characterise 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  values, as seen in Figure 5.For the majority (131) of our new variables we see that the predicted modulation is equal to or greater than that measured by T P 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  values.This is likely caused by the heteregeneous 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  0 = 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.

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, 9 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 Figure 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 datasets).Indeed, SAX J1808 received 1 'Unsure' vote which, had it instead been for 'Transient/Variable' would have pushed this subject above our classification threshold.
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 T P as they become resolved and move away from the core of the jet.In total there are 6 moving, transient jet components that are detected by the T P as unique, of which 3 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.

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 datapoints 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, whilst 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 datapoints and an additional detection and one with only two detections and two upper limits (see Figure 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.

COUNTERPARTS AND ASSOCIATIONS
We make use of 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 minute exposures of a given field, alternating between q-band (440-720nm) and each of 5 Sloan bands u, g, r, i and z.Whilst 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 database at a 2 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 the uncertaintyweighted mean position of each of our candidates as returned by the TraP.Running this crossmatch returns 25 counterparts in the Meer-LICHT database and 143 3 upper limits.This low rate of optical counterparts is not surprising as all of our XRB fields are within 10 degrees of the Galactic plane and so many optical counterparts may be heavily extincted.
The radio-optical plane for our variables can be seen in Figure 7, with comparison source types from Stewart et al. (2018)  7 .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 andAndersson 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.

Highlights
We also search for counterparts at other wavelengths and with preexisting classifications, with the aid of pre-existing code available in Driessen (2021).Crossmatching makes use of the astroquery package to search the SIMBAD 8 database (Wenger et al. 2000) and several catalogues within Vizier 9 (Ochsenbein et al. 2000).We again search at a 2 radius to our radio sources.We searched for X-ray and gamma-ray counterparts to our radio variables.To do this we crossmatched with catalogs 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).Below we make note of a few interesting objects returned in our search with IR or radio detections, or otherwise known counterparts.

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 7 Code for reproducing this plot is available at https://github.com/4pisky/radio-optical-transients-plot 8 The Set of Identifications, Measurements and Bibliography for Astronomical Data, available online at http://simbad.cds.unistra.fr/simbad/ 9https://vizier.cds.unistra.fr/viz-bin/VizieRback almost 50 years (Evans et al. 1976).This asymptotic giant branch (AGB) star is heavily dust-obscured at optical wavelengths, but very bright at IR wavelengths -1 = 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 Figure 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 −5  yr −1 (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.
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 ∼1730d (Samus' et al. 2017), with more recent estimates from the WISE W1 and W2 bands at 1950 ± 50 and 1520 ± 20 days 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 onto which material is shed (Allen 1984).Radio emission has been seen from symbiotic binaries (Seaquist et al. 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 ionised 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 ionisation (e.g.HeII or [OIII]) we cannot claim that this is a symbiotic system.There is overlap between maser systems and Mira-type symbiotic binaries (Seaquist et al. 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 4 non-detections before the system reached ∼3mJy 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 ionisation).

Pulsar -BfS 20
Pulsar PSR B1845-01 (J1848-0123), whose light curve can be found in Figure 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  ∝   -of −1.3 ± 0.3 (McEwen et al. 2020).
The observed variability is consistent with RISS (see section 4), with a predicted modulation greater than the measured  ∼ 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 −3 (McEwen et al. 2020; compared to typical DMs of order 10s of pc cm −3 for pulsars with || > 25 • ; Manchester et al. 2005), which is known to be linked to long time-scale, low amplitude refractive scintillation (Stinebring et al. 1990(Stinebring et al. , 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 minute epochs consisting of 8s 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).

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 board10 .This source is in the field of MAXI J1820 but outside the 0.5 degree radius set by R22.Its light curve can be seen to vary very smoothly in Figure 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 There were 4 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.
density of 1.9 ± 0.3mJy.There are no counterparts to this source in any higher energy band, despite this field being our furthest from the Galactic plane ( ∼ 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  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.

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 T P 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 Figure 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 3 votes as a transient/variable and 1 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, whilst 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 dataset: 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 emphasising the use of nondetections 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 T P 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-hour 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 realtime 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, whilst 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 datasets 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 prioritise sources of interest and optimise 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 optimise the volunteer classification of optical galaxies (Walmsley et al. 2022a).We are currently applying Lochner & Bassett (2021)'s A 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 prep).

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 made use of 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.
Table A1: Table of transients and variables found by volunteers during BfS:MKT, with their associated positions, median flux densities ( med ) and variability parameters  and .The observed date given here is that in which the source was detected at highest S/N by the T P, whilst the distance recorded is how far an entry is from the pointing centre of that observation.The Transient Fraction (TF) is the fraction of classifications as a transiet/variable, from 10 volunteers.Note that their Zooniverse subject ID is a unique identifier from which one can find the image and light curve online and freely accessible at https://www.zooniverse.org/projects/alex-andersson/bursts-from-space-meerkat/talk/subjects/ID, replacing ID with the numeric value below.This table will be available in machine-readable form online.

Name
Subject

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

Figure 2 .Figure 3 .
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.MNRAS 000, 1-19 (2022)

Figure 4 .
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, whilst marginal distributions of  and  can also be seen.Known variable sources (e.g.XRBs) in our fields are circled.Imaging artefacts appear at low  and high  , whilst flux calibration uncertainties can produce high , low  sources (due to lack of systematic uncertainty in equation 1.Most known transients are found by citizen scientists, whilst many new sources are identified and show a wide spread of values in this parameter space.

Figure 5 .
Figure 5.The ratio between observed  and the variability predicted by Hancock et al. (2019)'s model for Refractive Interstellar Scintillation (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.

Figure 6 .
Figure 6.The T P 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 whilst the others were not.We believe these to be caused by bad figure generation and lack of data points respectively.

Figure 7 .
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).Black crosses denote counterparts within the MeerLICHT database whilst grey triangles are upper limits.Diagonal lines denote a constant ratio between radio and optical flux density, whilst the   marker indicates the horizontal displacement caused by 5 magnitudes of optical extinction.The majority of our radio sources are likely extragalactic as they overlap in parameter space with quasars and GRBs.

Figure 8 .
Figure 8. Upper: Light curve of OH Maser 30.1-0.7 picked up by volunteers.There were 4 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.

Figure 9 .
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.

Figure 10 .
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.

Table 1 .
Properties of the 11 ThunderKAT datasets used in this work.Each field's approximate Galactic latitude is given for relevance in section 4 and Figure5.The number of sources and central RMS values are calculated by the T P (see section 2).