Search and identification of transient and variable radio sources using MeerKAT observations: a case study on the MAXI J1820+070 field Open Access

ABSTRACT

Many transient and variable sources detected at multiple wavelengths are also observed to vary at radio frequencies. However, these samples are typically biased towards sources that are initially detected in wide-field optical, X-ray, or gamma-ray surveys. Many sources that are insufficiently bright at higher frequencies are therefore missed, leading to potential gaps in our knowledge of these sources and missing populations that are not detectable in optical, X-rays, or gamma-rays. Taking advantage of new state-of-the-art radio facilities that provide high-quality wide-field images with fast survey speeds, we can now conduct unbiased surveys for transient and variable sources at radio frequencies. In this paper, we present an unbiased survey using observations obtained by MeerKAT, a mid-frequency (∼GHz) radio array in South Africa’s Karoo Desert. The observations used were obtained as part of a weekly monitoring campaign for X-ray binaries (XRBs) and we focus on the field of MAXI J1820+070. We develop methods to efficiently filter transient and variable candidates that can be directly applied to other data sets. In addition to MAXI J1820+070, we identify four likely active galactic nuclei, one source that could be a Galactic source (pulsar or quiescent XRB) or an AGN, and one variable pulsar. No transient sources, defined as being undetected in deep images, were identified leading to a transient surface density of <3.7 × 10⁻² deg⁻² at a sensitivity of 1 mJy on time-scales of 1 week at 1.4 GHz.

1 INTRODUCTION

The past decade has seen a renaissance of the radio transient sky. While a number of transient and variable radio sources were known for many years from targeted searches of sources discovered at other observing frequencies, for example X-ray binaries (XRBs), active galactic nuclei (AGNs), and gamma-ray burst (GRB) afterglows, the typical radio transient sky was not well probed. The rapid development of new instrumentation has enabled us to conduct large-scale surveys to systematically explore the radio transient sky over a range of time-scales. For example, at high time resolution, typically <1 s, this led to the discovery of a new category of radio transient sources referred to as fast radio bursts (FRBs; Lorimer et al. 2007) that are pushing to the extremes of physics and enabling us to probe the contents of the Universe (see Petroff, Hessels & Lorimer 2019, for a recent review of FRBs).

At a lower time resolutions, typically >1 s, new wide-field imaging software (e.g. wsclean; Offringa et al. 2014), telescopes with good instantaneous uv-coverage, and significant advances in computational power have enabled large sky areas to be imaged on multiple time-scales. This has led to a number of large surveys for transient and variable sources ranging from seconds (e.g. Kuiack et al. 2021) to tens of years (e.g. Bhandari et al. 2018; Radcliffe et al. 2019) at a wide range of observing frequencies in the radio spectrum. Initially, these surveys were only capable of identifying bright sources in a large sky area (e.g. Bower & Saul 2011) or faint sources in a much smaller sky area (e.g. Frail et al. 2012). New facilities have now come online that can quickly survey large sky areas to high sensitivity. While transient radio astronomy focuses on a wide range of observing frequencies (30 MHz–150 GHz; e.g. Whitehorn et al. 2016; Varghese et al. 2019), in the following work we shall consider observations at ∼1.4 GHz to enable direct comparison to our observations. At mid-frequencies, the MeerKAT (Meer Karoo Array Telescope; Camilo 2018) and ASKAP (Australian SKA Precursor; Hotan et al. 2021) radio telescopes have recently started observing with unprecedented sensitivity and rapid survey speeds, enabling deep probes of the radio transient sky on a wide range of time-scales, leading to the discovery of transient and variable sources (e.g. Driessen et al. 2020; Wang et al. 2021).

Transient sources are typically considered to be those sources which appear and then disappear due to a cataclysmic event, while variable sources are those whose flux density varies between some minimum and maximum value. However, the observational division between transient and variable sources is more complex. For instance, a source that appears and disappears during an observation could be truly transient but, on the other hand, it could be a variable source with a minimum flux density that is below the detection threshold in the image. In this work, we define transient sources as those appearing and disappearing during the observations with the caveat that they may be previously unknown variable sources. We note that this bimodal classification and the definition of transient versus variable sources are subject to debate within the literature. As we now have a growing population of these sources, we will likely move away from this bimodal classification and towards a classification system based upon specific object types such as supernovae, AGNs, etc.

While many unbiased transient surveys have been conducted at 1.4 GHz (e.g. Bower et al. 2010; Bell et al. 2011; Hodge et al. 2013; Mooley et al. 2016), only a small number of transient sources have been identified at 1.4 GHz (e.g. Levinson et al. 2002; Thyagarajan et al. 2011; Aoki et al. 2014, see Table B1 for more details). Levinson et al. (2002) identified 25 transient candidates on time-scales of ∼1 yr and they are believed to be orphan afterglows of radio afterglows from GRBs, radio supernovae, and radio-loud AGNs. Thyagarajan et al. (2011) found 71 transient candidates on a wide range of time-scales from 3 min up to 1 yr using all the images obtained as part of the FIRST survey (Becker, White & Helfand 1995). Many of the transient candidates have no known counterpart or identification, with the majority being consistent with galaxies and quasi-stellar objects. Finally, Aoki et al. (2014) conducted a shallow survey for transients with flux densities greater than 3 Jy covering a very large area of sky and found one highly significant transient event on the time-scale of 1 d. Thus, there are transient sources detected across a range of time-scales at 1.4 GHz. Further observations are required to determine the nature of these sources and their populations.

Variable sources have proven to be more prolific at ∼1.4 GHz, with known sources being part of targeted monitoring campaigns (e.g. Bright et al. 2020) and the discovery of variable sources in wide-field searches of survey data sets (e.g. Murphy et al. 2021; Driessen et al. 2022a). The variable sources discovered in wide-field surveys vary on a wide range of time-scales and include flaring stars (e.g. Driessen et al. 2020, 2022b), pulsars (e.g. Murphy et al. 2021), tidal disruption events (e.g. Anderson et al. 2020), and AGNs (e.g. Murphy et al. 2021). A summary of expected variable and transient sources on specific time-scales is given in Pietka, Fender & Keane (2015). Wide-field surveys with multiple snapshots are ideal for identifying variable sources at 1.4 GHz, leading to the discovery of new variable sources and examples of rare extreme variability for known sources.

In this paper, we present a transient and variability search using observations obtained as part of the ThunderKAT XRB monitoring campaign (Fender et al. 2016). Specifically, we used observations of the field of MAXI J1820+070 that were attained on a roughly weekly cadence (Bright et al. 2020). In Section 2, we outline the observational data used, the processing strategy, and image quality control. In Section 3, we outline the strategy and results for two optimized searches for transient and variable sources. Section 4 considers the progenitors of the identified transient and variable sources from this survey.

Throughout this work, we adopt a cosmology with H₀ = 71 km s⁻¹ Mpc⁻¹, Ω_m = 0.27, and Ω_λ = 0.73.

2 MEERKAT OBSERVATIONS

The observations used in this publication were obtained as part of the ThunderKAT weekly monitoring programme for XRBs. The target source, MAXI J1820+070, was observed weekly during its active phases (Bright et al. 2020). 64 observations were obtained at a central frequency of 1.28 GHz, with a bandwidth of 0.86 GHz, and for a total integration time of 15 min per epoch. The observation dates are provided in Table A1. These data were processed using the oxkat¹ pipeline with standard parameters unless otherwise stated (Heywood 2020). Heywood et al. (2022) describe the oxkat pipeline in depth and it has three key stages. First, initial automatic flagging is conducted on the visibilities of the calibrator field using tricolour² (Hugo et al. 2022), followed by calibration of the data using both a primary and secondary calibrator. The primary calibrators used for these data, J1934-638 and PKS B1934-638, provided delay, bandpass, and frequency-independent amplitude and phase solutions. The gain solutions are then transferred to the secondary calibrator, J1733-1304, which is used to calculate the complex, frequency-dependent gain solutions. These solutions are then applied to the target data. The second phase of the oxkat pipeline conducts automated flagging of the target using tricolour and creates an initial image of the field and generates a mask. The third phase of oxkat conducts a masked deconvolution of the target data and predicts the model visibilities that are used in a final self-calibration and imaging step. A single self-calibration cycle was conducted using cubical (Kenyon et al. 2018) to obtain both phase and delay solutions. All imaging is conducted using wsclean (Offringa et al. 2014). The images per time-step are created at four frequency bands (0.96, 1.18, 1.39, and 1.61 GHz), together with a single combined multifrequency image (created using a fourth-order spectral polynomial fit to the data, for more details see Offringa & Smirnov 2017). Each image is 10 222 × 10 222 pixels, covering a field of view of 2.9 × 2.9 deg². Observations that failed to be processed by this automated pipeline were excluded from the sample. Finally, a primary beam correction was applied to the images using the oxkat tool pbcor_katbeam.py³ with image masking beyond the primary beam power point of 1 per cent.

We note that there may be slight flux density variations between the images caused by offsets in the absolute image flux density scaling. To correct for this, we need to assume that, on average, the sources detected in the images do not have significant intrinsic variability. We chose the first successful image (Observation 6) at the highest observing frequency (1.61 GHz; as it has the highest resolution) and extracted the sources in that image using pyse (Carbone et al. 2018). We carefully choose sources to correct for the absolute flux scale offset using the following criteria.

• We want to exclude any variations in the outer regions of the image, where the scaling may be affected by primary beam uncertainties. Thus, we extract all sources within 0.5 deg of the image centre, minimizing the primary beam uncertainties to 5 per cent.

• Sources must be bright and significantly above the image noise. We use all sources detected above 20σ.

• Sources must be point sources, as extended sources can lead to flux density variations caused by source extraction methodology. The sourcefinder used in this analysis is pyse, which is purposefully designed to handle point sources efficiently and is known to model extended sources less well. We chose the point sources by visual inspection of the sources in the image used to extract the sources.

Following these selection criteria, we obtain 13 point sources that can be used to correct for any systematic offsets in the absolute flux density scale of the image. We input the list of sources into the LOFAR Transients Pipeline (trap; Swinbank et al. 2015) using the monitoring list capability and ran the pipeline on all the images. The detection threshold was chosen such that no other sources were detected in the images. We take the average flux density of a source across the observations as the reference value to compare to the individual extractions of that source. For each image, we plot the extracted flux density of each source against its average flux density and fit a linear regression model through the origin with the numpy (Harris et al. 2020) least squares fitting algorithm. The gradient of the fitted linear model gives the flux density correction factor. The absolute flux density scale of the image is then corrected by dividing each pixel value by the flux density correction factor. This method is conducted for each observing frequency separately. The typical flux density correction factor is ∼4 per cent.

To assess image quality, we measure the rms noise variation in the inner eighth of the images. We fit a Gaussian distribution to the rms noise values to obtain the average rms value and the 1σ uncertainty. In Table 1, we give the average rms noise for the sample of images at each observing frequency. Images with rms noise values that are >3σ deviant from the average rms value are rejected. Additionally, images where the restoring beam ellipticity is >2 are rejected as this can cause issues in source association within the trap. Finally, images are checked by visual inspection. The images at 1.18 and 1.61 GHz from Observation 29, on 2019 May 18, are rejected due to having a very poor resolution. The images rejected are given in Table A1. All other images are deemed as being of suitable quality for transient and variability hunts.

3 TRANSIENT AND VARIABLE HUNT

We input the images attained in Section 2 directly into trap using default settings except for the parameters given in Table 2. The trap takes a time series of images and conducts basic quality control on these images prior to further processing. After the quality control step, trap searches for all sources in the images using the pyse source finder. Sources are then associated across time and frequency to produce multiwavelength light curves. If a source is not detected, a constrained fit of the flux density is obtained at the location of the source. Finally, variability metrics are calculated for the output light curves (Swinbank et al. 2015).

The detection_threshold parameter determines the signal-to-noise ratio used by the source finder to detect sources within the images; in this analysis we use different values for the transient and variability analysis in order to optimize those searches (see the following subsections). The extraction_radius parameter is the radius (in pixels) out to which sources are found by the source finder. The force_beam parameter, when set to true, assumes all sources are point sources as expected for transient sources on the time-scales and resolutions probed by this analysis. The new_source_sigma_margin parameter can be used to raise the detection threshold so that new sources need to be more significantly detected before they are labelled as being transient. In this analysis, we choose to turn-off the new_source_sigma_margin and instead filter transient candidates via other methods (see the following subsections). The elliptical_x parameter is used to filter images of lower quality; in our analysis any images where the restoring beam is significantly elliptical (i.e. B_major/B_minor > 2) are rejected from the analysis. Finally, the beamwidths_limit parameter is used to control how far, in units of the restoring beam, from a newly detected source the source association algorithm can search for associations.

Following the quality control steps outlined in the previous sections, 3–27 per cent of images were rejected. All of the rejected images were due to having high image noise and none were rejected due to the ellipticity of the restoring beam. In Table 1, we provide the number of images per observing frequency that remain in the final sample.

3.1 Transient hunt

Transient sources are defined as those sources that are newly detected during the time-scale of the observations. To determine a reliable detection threshold for new sources, we require to detect less than 1 false positive source in all the images we process. Due to correlated noise in radio interferometers, features in the images take the shape of the restoring beam. Thus, we need to calculate the false positive rate by using the number of independent beam elements in the image rather than the number of pixels. The probability of this occurring is given by |$P(X\le 1) = 1 - \frac{1}{N}$| (assuming a binomial distribution and taking the maximum likelihood of there being one false detection), where the number of trials N is given by the total number of pixels divided by the number of pixels in one restoring beam. We extract the pixel values for all the pixels contained within the extraction_radius used by trap (given in Table 2) and plot them in a histogram. In Fig. 1, we show an example histogram for the images obtained at a frequency of 0.96 GHz. The data set contains 2.5 × 10⁷ pixels per image, corresponding to ∼1.6 × 10⁹ pixels per frequency band.

To calculate the optimal detection threshold for each observing frequency, we assume that the noise pixels follow a Gaussian distribution. We note that the distribution does not follow a Gaussian distribution; this is due to an excess of positive and negative sidelobes around bright sources in the outer regions of the image. Direction-dependent calibration can significantly reduce these sidelobes; however, we chose not to do this as it could not be fully automated, and takes significantly longer to process and the sidelobes can be easily distinguished from transient candidates (see Section 3.1). Even though these data do not follow a Gaussian distribution, this assumption is good enough to set a threshold for the initial search for transients. The pixel values also contain real sources detected in the images, thus we conduct clipping of the pixel data to remove those sources. We calculate the median value of the pixels and then remove all pixels that are outside of two standard deviations from the median. After clipping, the remaining noise pixels are fitted with a Gaussian distribution and we show the parameters in Table 3. Using the Gaussian distributions, we calculate the sigma detection threshold required to provide P(X ≤ 1). To be conservative in our transient identification, we use the maximum detection threshold from the four observing frequencies. We find that the detection threshold providing less than 1 false positive for all observing frequencies should be 5.56σ.

To search for new sources, we run trap with a 4σ detection threshold on all of the images at once in order to obtain a list of all new sources detected throughout the data set. For this low detection threshold, the source association becomes increasingly complex and can fail. Therefore, we significantly simplify the source association to sources within 1 restoring beamwidth of each other in the image being processed. This simplification may lead to missed source associations for sources within the image but this does not hinder the transient search due to the filtering strategy outlined below. The number of new sources (those not detected in the first time-step processed) detected is 42 090, across the four observing bands. We note that this number does not exclude multiple detections of the same source across the frequency bands. These sources are filtered using the following criteria.

• Extract only sources from the new source list that are detected with a signal-to-noise ratio above the calculated detection threshold of 5.56σ. This reduces the number of new sources by 99 per cent to 334 remaining.

• Reject any sources that can be associated with the existing sources in the data base. These are likely caused by source association issues within the pipeline. This results in a reduction of new sources by 31 per cent, bringing the candidate list to 231 sources.

• Reject any sources close to the source extraction radius. The trap measures the noise within the extraction radius by subdividing the area into smaller blocks and extrapolating the noise between these blocks (see Swinbank et al. 2015, for further details). This can lead to an underestimate of the noise around the source extraction radius and thus a number of false positive transient detections (e.g. Rowlinson et al. 2016). We reject any new sources identified within 3 arcmin of the source extraction radius. This removes 9 per cent of candidates remaining from the second filtering step, leading to 210 sources remaining on the candidate list.

• Reject any candidates only detected in one frequency band. Often in radio images there are noise features surrounding bright sources, known as sidelobes, which can be mistaken for new sources. Additionally, correlated noise features can also lead to false positive source detections. These noise features are an artefact of the images being made via an incomplete Fourier transformation on sparsely sampled data. By conducting self-calibration, as with the data presented in this paper, it is possible to significantly reduce these artefacts but some remain. The locations of these noise features are highly dependent on the observing frequency, meaning that one of these false positive results will not have a corresponding source at the same time and location at a different observing frequency. Thus, we require that all new transient candidates are detected in at least two frequency bands in the same time-step. At least one detection must exceed the detection threshold of 5.56σ and at least two detections must exceed the 4σ detection threshold used in trap. This reduces the number of candidates remaining from the previous step by 98 per cent, leading to a total of four candidates.

• Confirm that candidates are not underlying faint sources near the detection threshold used in this analysis. We create deep images for a given observing frequency by summing all the pixels in each image obtained and dividing the resultant pixel values by the number of images used (note that this is not as good as re-imaging all of the calibrated visibilities together but is sufficient for this analysis). To confirm the flux density scale in these averaged images is reasonable; we plot the observed flux density in the mean image against the average flux density measured by the trap for all sources with a detection significance above 8σ, and there is no significant systematic offset between the fluxes. We then reject any candidates detected in the deep image for the detection frequency where the candidate maximum flux density is more than 5σ deviant from the flux density measured in the deep image. Sources not associated with a source in the deep image are also retained as transient candidates. This removes none of the candidates from the previous filtering step, leading to four sources remaining on the candidate list.

• Conduct visual inspection of all remaining candidates using both the individual images and the deep image of the field. This leads to the rejection of one candidate as it is comparable to noise features in the nearby region.

• Run trap with the remaining three candidates in a monitoring list to obtain a full light curve and the variability parameters. Their properties are given in Table 4.

In summary, three candidate sources pass these criteria, although all of them are associated with sources in the deep images of the field and they are defined to be variable sources. Their light curves are shown in Figs 2–4. The method outlined above shows an efficient method to automatically identify reliable transient candidates in a fully automated manner using the standard outputs of the trap.

3.1.1 Transient surface density

The sources detected in this analysis are all associated with faint sources in the deep images and are observed in multiple epochs. These are identified as variable sources rather than the transient sources typically considered in standard transient searches. Therefore, we do not include these sources in comparisons to other transient surveys. In order to compare this transient search to other published searches, we need to determine the transient surface density at a single observing frequency. We use the data presented in Table B1 and obtained from Mooley et al. (2016).⁴ We choose to use the images observed at 1.39 GHz as these images were the best quality and we are comparing this result to other surveys conducted at 1.4 GHz. The transient surface density is defined as being the number of transients detected per square degree surveyed and is typically compared to the sensitivity of the observations (the faintest transient detectable in the survey).

The images typically have structure in their noise properties; for instance the quality decreases with increasing distance from the image centre due to the response of the primary beam. Thus, we sample a small area to high sensitivity and a significantly larger area to a lower sensitivity limit.

Following a similar strategy to that outlined by Kuiack et al. (2021), we divide each image into 100 annuli, determine the sky area covered by each annulus, and then calculate the rms noise in each annulus (following a simple sigma clipping, with a 3σ threshold, to remove bright sources from the data). Following this step, we determine the maximum and minimum rms values and create a range of values between this limit with 100 bins. The sky area surveyed by each unique annulus within a given rms bin is then summed up. Finally, for each rms bin we sum up the total sky area surveyed for all rms bins less than or equal to this bin. Hence, the final bin, at the highest rms value, has the total area surveyed in all the images. As we search out to a radius of 0.9 deg, the area surveyed in one image is equal to Ω = 0.9²π deg². Thus, the total area surveyed is 0.9²π × (N − 1) deg⁻², where N is the number of images surveyed. The ‘−1’ in this expression is because one of the images is needed as a comparison image and hence not included in the total area. For large data sets, (N − 1) tends towards N, therefore the summation method used above will tend towards the total area surveyed.

To convert this to a transient surface density ρ, the number of transients T found is divided by the total area surveyed, i.e. |$\rho = \frac{T}{\Omega (N-1)} = \frac{T}{A}$|⁠, where A is the area calculated via the binning strategy outlined above. Following previous studies, we calculate the 95 per cent confidence limit that no transient sources were detected assuming Poissonian statistics (e.g. Bell et al. 2014; Rowlinson et al. 2016), and this equation approximates to |$\rho \simeq \frac{3}{A}$|⁠. The transient surface density is then calculated for each rms bin. To determine the sensitivity, or the faintest transient detectable, we multiply the rms value by the detection threshold used, i.e. σ × rms, where σ = 5.56 for this work.

Finally, we plot the transient surface density as a function of the sensitivity in the right-hand panel of Fig. 5 using the thick red line. The diamond at the end of this line represents the faintest transient detectable in the total area surveyed (i.e. the lowest transient surface density probed). In the left-hand panel of Fig. 5, we plot the lowest value of the transient surface density probed versus the median minimum separation time-scale of the images in days (7 d for these data). We find a transient surface density of <3.7 × 10⁻² deg⁻² for transients brighter than 0.1 mJy on time-scales of 1 week.

3.2 Variable hunt

In this analysis, variable sources are defined as those sources being detected in the first observation of the field and showing significant variation in their flux densities throughout the observing campaign. As variability is difficult to quantify for faint sources, we use a higher detection threshold than in the previous section, i.e. 8σ. Additionally, as we are concerned with variability of sources, it is important to not miss source associations. Using a higher detection threshold reduces the complexity of the source association procedure and we can use the default value of 3 for the trap beam widths limit to prevent missed source associations. Variable sources are typically identified using the two key variability parameters calculated by the trap: the reduced weighted χ² of a fit assuming a constant brightness, η, and the coefficient of variation, |$V = \frac{s}{\overline{I}}$| (where s is the standard deviation of the flux density measurements and |$\overline{I}$| is the average flux density of the source, as defined in Swinbank et al. 2015). These variability parameters are calculated for each unique source detected in the images. As we handle newly detected sources separately (see Section 3.1), in this section we only consider sources detected in the first observation. We process each observing frequency separately using trap, as the data quality between the observing frequency differs and this can impact the observed variability. In Fig. 6, we plot the variability parameters for each of the sources at the four observing frequencies. We can model the distribution of stable sources as having a Gaussian distribution in both η and V, with outliers being candidate variable sources (Swinbank et al. 2015; Rowlinson et al. 2019). We define variable candidates as being >1.5σ offset from the η and V distributions. In Table 5, we provide the 1.5σ thresholds used for each observing frequency and they are shown as the dashed black lines in Fig. 6. This leads to four sources being identified as variable sources, including the XRB MAXI J1820+070 at the centre of the image as expected (highlighted with the red circle in Fig. 6). The variable candidates are listed in the bottom part of Table 4 with their variability parameters for each observing frequency. Variability is confirmed by visual inspection and comparison to nearby sources. The multiwavelength light curves of these sources are shown in Figs 7–9. Finally, we searched for associations within the NVSS 1.4 GHz radio catalogue (Condon et al. 1998). Each of these variable sources are associated with known sources and their associations are given in Table 5.

In addition to the variable sources identified by the trap in this analysis, we also consider the variable candidates found in the transient search conducted in Section 3.1 given in the top part of Table 4. We plot their variability parameters on Fig. 6 using red ‘+’ symbols and note that they do not pass the variability threshold on η as outlined in Table 5. This is likely due to the majority of flux density measurements being non-detections and hence the η parameter is instead being dominated by the low variability in the noise in the images.

In summary, we have identified three variable sources in addition to the target of these observations, MAXI J1820+070 (which is considered in a separate publication; Bright et al. 2020). In the following section, we will discuss these sources in more depth.

4 ANALYSIS OF NEW VARIABLE SOURCES

In this paper, we have identified three new variable sources via the transient search outlined in Section 3.1, and three new variable sources via the standard variability search given in Section 3.2. Here, we consider their multiwavelength properties to determine the likely progenitor source. We use key diagnostic plots from Stewart et al. (2018)⁵ and Pietka et al. (2015)⁶ to compare the properties of our sources to the known populations of transient and variable sources.

4.1 PSR J1822+0705

We searched the existing multiwavelength catalogues by querying the location in simbad (Wenger et al. 2000), and it is associated with PSR J1822+0705. PSR J1822+0705 is a regular pulsar with a spin period of 1.36 s, a flux density of 3.8 ± 0.1 mJy at 400 MHz, and a dispersion measure of 62.2 pc cm⁻³ (from the ATNF pulsar catalogue; Manchester et al. 2005).⁷ This pulsar has not been extensively studied and was not previously known to be variable.

Pulsars are known to vary on short time-scales, due to intrinsic effects, and on long time-scales, due to both interstellar scintillation and intrinsic effects. To check if the variation observed is also seen on short time-scales, we take the observation that recorded the brightest flux density for this pulsar (observation id 1543743065 on 2018-12-02 at 0.96 GHz) and re-imaged on the shortest time-scales of 8 s, producing 112 snapshot images. In Fig. 10, we show three consecutive 8 s images showing that the variation appears to be present on the shortest time-scales we can image. The time-scale of this variation is much shorter than both diffractive interstellar scintillation (minutes to hours; Hewish et al. 1968) and refractive interstellar scintillation (days to years; Sieber 1982). Therefore, this observed variation is likely to be intrinsic to the pulsar.

4.2 MKT J182015.5+071455

This newly identified source is faint, with a maximum flux density of 0.2 ± 0.1 mJy (at 1.39 GHz). It is visible in the deep images of this field, with a flux density of 0.089 ± 0.013 mJy (at 1.39 GHz), showing that the source is likely persistent and has increased in brightness by a factor of ∼2 at its maximum. At this flux density, it is nearly an order-of-magnitude fainter than typically assumed limits for radio surveys like NVSS. There are no optical sources associated with this source in the PanSTARRS archive (Flewelling et al. 2020), with a lower limit of 22.7 mag in the i band. The deep radio image and optical image are shown in Fig. C1. This source is estimated to have a variability time-scale of ∼1 week as it was typically observed in just one epoch with non-detections in the previous and following epochs.

From Fig. 11, we note this source lies in a region of the plot where the optical survey used as a reference for AGN (SDSS) is not complete. While it is consistent with GRB afterglows in Fig. 11, this is unlikely given that the source appears to be persistent. Although poorly sampled, there are pulsars and stellar sources in this region. While there is potential of optical extinction towards the source, we note that the Galactic colour excess towards the source is E(g-r) = 0.25 mag, suggesting A_I ≲ 0.5 for Galactic sources (Green et al. 2019). The most likely conclusion is that this source is consistent with a fainter AGN population than was sampled by SDSS.

In Fig. 12, we consider the variability time-scale of ∼2 weeks in comparison to other transient and variable populations. As the distance to this source is unknown, we take two typical examples: a Galactic origin with a distance of 5 kpc, and an extra-galactic origin with a distance 100 Mpc. If this source is Galactic, we determine that it is most likely an RS CVn or an XRB; however, the typical optical luminosities of both classes argue against this conclusion. Toet et al. (2021) have detected RS CVn at 144 MHz with flux densities of ∼0.5 mJy), but these sources were only at 15–300 pc (with the sources above 100 pc always having at least one subgiant or giant star). At 300 pc, even low-mass dwarf stars should have been detected by PanSTARRs. In Fig. 6, the optical fluxes of typical XRBs in outburst were 100 times brighter than the optical limit. While we cannot rule out quiescent accretion on to an XRB, only the quiescent radio flux of the black hole XRB V404 Cyg (∼2.4 kpc; Miller-Jones et al. 2009) has been observed at ∼0.1 mJy (e.g. Rana et al. 2016). If it is extra-galactic, it may be consistent with a supernova; however, the persistent nature of this source discounts this possibility, so it would then most likely be an AGN. For distances beyond the assumed 100 Mpc, the source becomes more consistent with the faint end of AGN, albeit with a shorter variability time-scale. We note that the latter can be compressed by relativistic effects, but that is usually more severe for the brighter blazars. Although we cannot confidently conclude which type of astrophysical variable this source is, the leading contenders are a pulsar, an atypical quiescently accreting XRB, or an AGN at a distance beyond 100 Mpc with a faster than average variability time-scale. Given that a large fraction of variable radio sources are AGNs, the latter may be the most likely of these conclusions.

4.3 MKT J182001.4+072945

As with the previous source, this newly detected source is faint with a flux density of 0.3 ± 0.1 mJy at 1.39 GHz. This source is detected in the deep image with a flux density of 0.13 ± 0.02 mJy (1.39 GHz), showing that the source is likely persistent and has increased in brightness by a factor of ∼3 at its maximum. We find this source is associated with the PanSTARRS source PSO J275.0057+07.4958, which is detected in all observing filters except the g band. The observed magnitudes of this optical sources are r = 21.26 ± 0.07, i = 20.58 ± 0.05, z = 20.27 ± 0.08, and y = 19.87 ± 0.05 (Flewelling et al. 2020). The deep radio image and optical image are shown in Fig. C2. This source is estimated to have a variability time-scale of ∼2 weeks as it was typically observed in two consecutive epochs with non-detections in the previous and following epochs.

The maximum radio flux density of 0.5 mJy at 0.96 GHz and variability time-scale of ∼2 weeks is similar to MKT J182015.5+071455 above. Thus, the conclusions discussed in Section 4.2 based just on radio flux density and time-scale are largely the same. However, with an optical brightness of ∼20 mag, some of the Galactic scenarios dismissed need to be re-evaluated. At distances of 300 pc and 5 kpc, the absolute optical magnitude is roughly ∼12 and ∼6.5 mag, respectively. The former still tends to rule out most dwarf stars and thus the RS CVn category. However, the latter is consistent with a late K dwarf and suggests that a somewhat nearby XRB, particularly one in quiescence, cannot be ruled out. Although we cannot confidently conclude which type of astrophysical variable this source is, the leading contenders are a pulsar, an XRB, or an AGN at a distance beyond 100 Mpc with a faster than average variability time-scale. Again, an AGN is the most likely of these conclusions given the typical source densities of such sources.

4.4 NVSS J181849+062843

First, we searched for associations within the existing radio catalogues at 1.4 GHz, specifically NVSS (Condon et al. 1998) and VLASS (Gordon et al. 2021). We find that this source is associated with NVSS J181849+062843, and has flux densities of 2.2 ± 0.4 mJy (1.4 GHz; NVSS) and 13.8 ± 0.3 mJy (2–4 GHz; VLASS), respectively. We note that the average flux density at 1.39 GHz in these observations is five times higher than the observed NVSS flux density, confirming this source is variable on long time-scales.

We searched the existing multiwavelength catalogues by querying the location in simbad (Wenger et al. 2000) and no associated sources are found. To search for a fainter optical source association, we queried the PanSTARRS catalogue (Flewelling et al. 2020). We find that NVSS J181849+062843 is associated with the PanSTARRS source PSO J274.7073+06.4791, which was detected in just two of the PanSTARRS observing filters: i = 21.80 ± 0.20 mag and z = 20.88 ± 0.02 mag. In Fig. C3, we show the MeerKAT and PanSTARRS images of this source.

With an optical brightness of ∼21 mag and an average radio flux density of 10 mJy at 1.39 GHz, this source is most likely an AGN or quasar, as shown in Fig. 11. With a variability time-scale of ∼1–3 months (determined by visual inspection of peaks in the light curve), this is typically very short for observed AGN intrinsic variability (on time-scales >days to years; e.g. Hovatta et al. 2008) but not inconsistent as shown in Fig. 12. Thus, we conclude that this source is most likely an AGN.

4.5 NVSS J181752+064638

This bright variable source is associated with NVSS 181752+064638 with an NVSS flux density of 100.8 ± 3.1 mJy at 1.4 GHz (Condon et al. 1998). This source is also detected in the VLASS survey with a flux density of 73.0 ± 0.3 mJy at 2–4 GHz (Gordon et al. 2021). We find this source is on average 10 per cent brighter in our MeerKAT observations. The VLASS flux density at 2–4 GHz is ∼30 per cent fainter, signifying either spectral behaviour or variability.

We searched the existing multiwavelength catalogues by querying the location in simbad (Wenger et al. 2000), and NVSS J181752+064638 has a flat spectrum and has been identified as a blazar (Healey et al. 2007). Additionally, we find this source is associated with the PanSTARRS source PSO J274.4693+06.7773, which is detected in all observing filters except g. The observed magnitudes of the optical counterpart are r = 21.80 ± 0.02, i = 21.45 ± 0.03, z = 21.09 ± 0.21, and y = 20.39 ± 0.19 (Flewelling et al. 2020). In Fig. C4, we show the MeerKAT and PanSTARRS images of this source.

With an optical brightness of ∼21 mag and an average radio flux density of 113 mJy at 1.39 GHz, this source is most likely an AGN, as shown in Fig. 11, further supporting the blazar identification. A variability time-scale of ∼3–6 months (as determined by visual inspection of peaks in the light curve) is relatively short for an AGN but is not inconsistent as shown in Fig. 12. Thus, we conclude that this source is also most likely an AGN.

4.6 NVSS J182029+063419

This source was also associated with sources in NVSS (NVSS J182029+063419) and VLASS, with flux densities of 5.8 ± 0.4 mJy (1.4 GHz) and 13.4 ± 0.3 mJy (2–4 GHz), respectively (Condon et al. 1998; Gordon et al. 2021). Our MeerKAT observations show that the flux density has approximately doubled since the NVSS observations.

We searched the existing multiwavelength catalogues by querying the location in simbad (Wenger et al. 2000) and no associated sources were found. We queried the PanSTARRS catalogue and find that the bright point source PSO J275.1254+06.5721 is associated. The optical source is detected in all of the PanSTARRS filters with observed magnitudes of g = 18.89 ± 0.02, r = 18.27 ± 0.01, i = 18.18 ± 0.01, z = 17.91 ± 0.02, and y = 17.78 ± 0.03 (Flewelling et al. 2020). In Fig. C5, we show the MeerKAT and PanSTARRS images of this source.

With an optical brightness of ∼18 mag and an average radio flux density of 13 mJy at 1.39 GHz, this source could be either a quasar or an XRB, as shown in Fig. 11. With an observed variability time-scale of ∼1 yr (determined by visual inspection of peaks in the light curve), this variable source could be an XRB, but is more likely an AGN. We show this source in Fig. 12, assuming it is extragalactic at a distance of 100 Mpc. We note that, with a Galactic distance of 5 kpc, this source could still be associated with XRBs or novae.

Therefore, we determine that this source is either an AGN, a nova, or an XRB. Further, radio observations or an optical spectrum will be able to confirm this identification.

5 DISCUSSION AND CONCLUSIONS

In Section 3.1, we presented a transient hunt on 1 week time-scales for 62 epochs and a search radius of 0.9 deg, giving a total sky area surveyed of 155 deg² to a limiting sensitivity of 1 mJy. A known pulsar was found via this search, and also two transient candidates that were reclassified as variables, as they were associated with persistent emission in deep images. The 1 week time-scales for transient sources at ∼1.4 GHz is a relatively unexplored parameter space, as shown in the left-hand panel of Fig. 5, with only one other survey on comparable time-scales by Bell et al. (2011). The transient surface density limit presented in this work is similar to that attained by Bell et al. (2011) using 5037 epochs of observations from the Very Large Array; however, our survey is an order-of-magnitude more sensitive than that survey and using two orders of magnitude fewer observations. Thus, with a relatively small sample of MeerKAT data, we are able to place highly competitive limits on the presence of transients on time-scales of ∼1 week. Significantly more observations on these time-scales have been obtained by monitoring XRBs using MeerKAT; when processed, these will be able to determine if there is a rare population of transient sources on these time-scales.

The known pulsar PSR J1822+0705, which was detected as part of the transient hunt, was previously poorly studied and was not known to show variable emission. We identified it as varying on weekly time-scales and further investigation of the brightest epoch demonstrated that it is varying on time-scales of seconds. This short time-scale variability indicates that it is intrinsic to the pulsar and is not caused by interstellar scintillation. From studies of known variable pulsars (e.g. Parent et al. 2022), there are a number of possible causes for intrinsic variation including giant pulses (e.g. Lundgren et al. 1995), nulling pulsars (e.g. Backer 1970), intermittent pulsars (e.g. Kramer et al. 2006), and variable flux density emission states (e.g. Young et al. 2014). Further detailed analysis is required to determine the origin of the intrinsic variability of PSR J1822+0705. Imaging surveys are an excellent method to detect extreme pulsar behaviour and a method to discover new highly variable pulsars. This is because large areas of sky can be monitored over long time-scales, which can be challenging with traditional pulsar timing methods (e.g. Hobbs et al. 2016).

In addition to PSR J1822+0705, we detected two other likely variable sources. Based on their variability time-scales and radio luminosities, we determined that they are most likely AGNs though one source could also be a Galactic pulsar or a quiescent XRB. Though it is unexpected that AGNs vary on time-scales ∼1–2 weeks, this may be subject to observational biases in AGN surveys. Alternatively, the observed variability time-scale may only show the peak of the emission due to the sources being near to the detection threshold in this survey and the ‘true’ variability time-scale could be much longer.

In addition to a transient search, we also conducted a search for variable sources in the field in Section 3.2. We identified three likely AGNs/blazars that are varying on time-scales from ∼1 month up to ∼1 yr. Although not unexpected, the number of AGNs known to vary on these time-scales is relatively small (e.g. Pietka et al. 2015), though the known population is growing rapidly (Sarbadhicary et al. 2021; Driessen et al. 2022a). By identifying the variability of these sources via imaging surveys, we are able to significantly increase the known population of variable AGNs, especially on short time-scales, leading to a greater understanding of these sources. Following identification, dedicated AGN studies can determine if the variability is an intrinsic effect caused by, for example, variable accretion rates on to the central black hole or if the variability is caused by extrinsic effects such as scintillation.

In summary, this paper has presented an optimized transient and variability search method for MeerKAT observations. The transient limits attained are competitive despite this data set being comparatively small. Finally, we identified four AGNs varying on relatively short time-scales, one source that may be a pulsar or a quiescent XRB (though, based on statistical arguments, also most likely to be an AGN), and one pulsar demonstrating extreme variability on both short (8 s) and long (weeks) time-scales.

ACKNOWLEDGEMENTS

We thank Ben Stappers for useful discussions regarding this work. AR acknowledges funding from the NWO Aspasia grant (number: 015.016.033). PAW acknowledges financial support from the University of Cape Town and the National Research Foundation. GRS is supported by NSERC Discovery Grants RGPIN-2016-06569 and RGPIN-2021-0400.

We acknowledge use of the Inter-University Institute for Data Intensive Astronomy (IDIA) data intensive research cloud for data processing. IDIA is a South African university partnership involving the University of Cape Town, the University of Pretoria, and the University of the Western Cape.

The MeerKAT telescope 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. We would like to thank the operators, SARAO staff, and ThunderKAT Large Survey Project team.

This research used astropy, a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018). We also use numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), pandas (The pandas development team 2020), sqlalchemy,⁸ and matplotlib (Hunter 2007). This research has used the SIMBAD data base, operated at CDS, Strasbourg, France.

The Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, the Queen’s University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation Grant No. AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation.

DATA AVAILABILITY

The data and scripts underlying this article are available in Zenodo, at https://zenodo.org/record/6826588#.Ys60VC8Rpqs. The MeerKAT data used are available in the MeerKAT data archive at https://archive.sarao.ac.za/. The PanSTARRS fits images were obtained from https://panstarrs.stsci.edu.

Footnotes

REFERENCES

APPENDIX A: OBSERVATIONS

APPENDIX B: PREVIOUS TRANSIENT SURVEYS

 

APPENDIX C: RADIO AND OPTICAL IMAGES OF VARIABLE SOURCES