The ASKAP Variables and Slow Transients (VAST) Pilot Survey

Abstract

The Variables and Slow Transients Survey (VAST) on the Australian Square Kilometre Array Pathfinder (ASKAP) is designed to detect highly variable and transient radio sources on timescales from 5 s to $\sim\!5$ yr. In this paper, we present the survey description, observation strategy and initial results from the VAST Phase I Pilot Survey. This pilot survey consists of $\sim\!162$ h of observations conducted at a central frequency of 888 MHz between 2019 August and 2020 August, with a typical rms sensitivity of $0.24\ \mathrm{mJy\ beam}^{-1}$ and angular resolution of $12-20$ arcseconds. There are 113 fields, each of which was observed for 12 min integration time, with between 5 and 13 repeats, with cadences between 1 day and 8 months. The total area of the pilot survey footprint is 5 131 square degrees, covering six distinct regions of the sky. An initial search of two of these regions, totalling 1 646 square degrees, revealed 28 highly variable and/or transient sources. Seven of these are known pulsars, including the millisecond pulsar J2039–5617. Another seven are stars, four of which have no previously reported radio detection (SCR J0533–4257, LEHPM 2-783, UCAC3 89–412162 and 2MASS J22414436–6119311). Of the remaining 14 sources, two are active galactic nuclei, six are associated with galaxies and the other six have no multi-wavelength counterparts and are yet to be identified.

1. Introduction

Radio variability is produced by a wide range of astronomical sources and phenomena, including the Sun, planets, stellar systems, neutron stars, supernovae, gamma ray bursts (GRBs) and active galactic nuclei (AGN) (Cordes et al. 2004). Among the most exciting and anticipated transient sources from extragalactic imaging surveys are neutron star mergers, tidal disruption events, and orphan GRB afterglows (Metzger et al. 2015). These can help answer questions regarding the true rates of these events (independent of beaming), their energetics, and the properties of their surrounding media. However, because the source of the emission is similar in most cases (shock-driven synchrotron emission), distinguishing these sources from other transients can be hard, especially with small numbers of detections well after the sources have faded. Therefore prompt identification and multi-wavelength comparisons are crucial (Ofek et al. 2011; Stewart et al. 2018).

At the same time, we must distinguish between extragalactic transients and Galactic foreground sources, primarily low-mass stars and pulsars (analogous to the situation in searches for fast optical transients, where M star flares are the dominant foreground; see e.g., Cowperthwaite & Berger 2015; Ho et al. 2018; Webb et al. 2020). Again, multi-wavelength identification can be crucial, but in some cases the properties of the radio emission itself can be used to identify high brightness temperatures and/or coherent emission that rule out extragalactic sources (e.g., Kaplan et al. 2008; Roy et al. 2010; Pritchard et al. 2021).

Early untargeted surveys generally used archival data, and consisted of either many repeats of relatively small sky areas (e.g., Carilli et al. 2003; Bower et al. 2007); a few repeats of large sky areas, such as the comparison between the NRAO VLA Sky Survey (NVSS; Condon 1997) and Faint Images of the Radio Sky at Twenty-centimetres (FIRST; White et al. 1997) surveys (Gal-Yam et al. 2006); or inconsistently sampled repeats over a range of sky areas (e.g., Bannister et al. 2011a). These surveys typically resulted in a small number of detections, but the historical nature of these detections, and the lack of multi-wavelength data made it difficult to classify or confirm these objects (e.g., Bower et al. 2007; Thyagarajan et al. 2011).

Despite the limitations of these surveys there have been a number of interesting discoveries. For example, Hyman et al. (2005) and subsequent papers identified several interesting sources towards the Galactic Centre whose origins are still unknown (see also Hyman et al. 2009), and Ofek et al. (2011) used near real-time reduction and multi-wavelength follow-up to identify a single transient that may be consistent with a neutron star merger.

More recently, the advent of Square Kilometre Array (SKA) pathfinder telescopes and upgrades of established telescopes have allowed a new wave of radio transients surveys (e.g., Driessen 2020). For example, the Caltech-NRAO Stripe 82 survey (CNSS; Mooley et al. 2016) has used the upgraded Karl G. Jansky Very Large Array (VLA) to repeatedly survey an equatorial strip and identified a number of sources, including a tidal disruption event (TDE; Anderson et al. 2020). Stewart et al. (2016) used repeated observations of the north celestial pole at 60 MHz with the Low Frequency Array (LOFAR; van Haarlem et al. 2013) and found a single transient event of unknown nature. The ongoing VLA Sky Survey (VLASS; Lacy et al. 2020) has led to the identification of several sources through comparison with archival surveys (Law et al. 2018; Ravi et al. 2021). The ThunderKAT large survey project (Fender et al. 2016) on the MeerKAT telescope (Jonas & MeerKAT Team 2016) will both conduct its own observations and analyse commensal data from other projects, and has already identified new sources including a chromospherically active K-type sub-giant (Driessen et al. 2020). Increasingly these detections are being made closer to real-time (within days of observation), aided by improved radio surveys (Lacy et al. 2020; McConnell et al. 2020) to allow rapid identification of new (or missing; e.g., Ravi et al. 2021) sources.

This progress at radio wavelengths is coupled with significant improvements (in depth, wavelength coverage, resolution, and temporal sampling) in multi-wavelength coverage, such as the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS; Chambers et al. 2016), the Vista Hemisphere Survey (VHS; McMahon et al. 2013), the VISTA Variables in the Via Lactea survey (VVV; Minniti et al. 2010), SkyMapper (Keller et al. 2007), the Dark Energy Survey (DES; Abbott et al. 2018), DESI Legacy Imaging Surveys (especially Dark Energy Camera Legacy Survey, or DECaLS; Dey et al. 2019), the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010) and, especially for stellar science, the Gaia mission (Gaia Collaboration et al. 2016). These surveys improve cross-matching and classification for a variety of sources (Stewart et al. 2018).

The Australian Square Kilometre Array Pathfinder (ASKAP; Johnston et al. 2007; Hotan et al. 2021) is an array of $36{\times}12$ -metre prime-focus antennas located at the Murchison Radio-astronomy Observatory in Western Australia. It operates over a frequency range of 0.7–1.8 GHz with an instantaneous bandwidth of 288 MHz and a nominal field of view of 31 square degrees. ASKAP was designed as a rapid wide-field survey instrument. The ASKAP Variables and Slow Transients Survey (VASTa; Murphy et al. 2013) is one of the key survey science projects (SSPs)b approved for ASKAP time allocation.

VAST was designed to detect astronomical phenomena that vary on timescales accessible in the ASKAP imaging mode (from $\sim\!5$ s to several years). Since the original survey proposal described by Murphy et al. (2013), the scope of the science case has been extended to incorporate detection of radio emission from gravitational wave events (in particular binary neutron star mergers) and fast radio bursts and their host galaxies (in the image domain). More rapidly varying sources (on timescales of less than seconds) are being explored as part of the Commensal Real-Time ASKAP Fast-Transients Survey (CRAFT; Macquart et al. 2010).

In the Early Science period, we have demonstrated ASKAP’s capability to detect time-variable phenomena in imaging mode. We have reported: an initial search for transients and variables in a single $30\,\mathrm{deg}^2$ field using 8 daily epochs and a subset of the final ASKAP array (Bhandari et al. 2018, also see Heywood et al. 2016); the discovery of a new millisecond pulsar J1431 $-$ 6328 (Kaplan et al. 2019); a sample of 23 previously undetected flaring stars (Pritchard et al. 2021); and radio bursts from UV Ceti (Zic et al. 2019) and Proxima Centauri (Zic et al. 2020). In addition, ASKAP was used to follow-up several gravitational wave events from LIGO/Virgo, notably GW190814 (Dobie et al. 2019; Abbott et al. 2020).

During 2019–2020, each of the SSPs was allocated 100 h of observing time to conduct a Phase I Pilot Survey with ASKAP. In 2021, another 100 h have been allocated to conduct a Phase II Pilot Survey. The purpose of these pilot surveys is to test the observing strategy and implementation, data processing capabilities, and scientific analysis in preparation for the full surveys, which are expected to begin in 2022.

In this paper we describe the VAST Phase I Pilot Survey, present some initial results, and discuss the plan for the VAST Phase II Pilot Survey. In Section 2 we discuss our observational specifications and scheduling; in Section 3 we discuss our transient detection pipeline software; in Section 4 we present results from an untargeted variability search on two of the VAST Phase I Pilot Survey regions (3 and 4); and in Section 5 we discuss our future plans.

2. Observations and survey strategy

2.1. The rapid ASKAP continuum survey

In addition to our dedicated VAST observations, we incorporated data from the Rapid ASKAP Continuum Survey (RACS; McConnell et al. 2020) as our first epoch. RACS is the first large-area survey to be conducted with the full 36-antenna ASKAP telescope. When complete, RACS will cover the sky south of declination $+51^\circ$ (a total sky area of $36\,656\ \mathrm{deg}^2$ ), across three ASKAP bands (centred on 888, 1296 and 1656 MHz). The angular resolution is $15{-}25\,\mathrm{arcsec}$ in the low band and $8-$ 27 arcsec in the mid band, depending on the declination.

The first data release of RACS consists of 903 images covering the sky south of declination $+41^\circ$ , observed in the low band (centred on 888 MHz). Each field was observed for $\sim\!15\,\mathrm{min}$ , achieving a typical rms noise of $0.25\,\mathrm{mJy\,beam}^{-1}$ . All four instrumental polarisation products (XX, XY, YX, and YY) were recorded to allow images to be made in Stokes parameters I, Q, U, and V. The data were processed using the ASKAPSoft pipeline (Cornwell et al. 2012) and are available through the CSIRO ASKAP Science Data Archive (CASDAc; Chapman et al. 2017) under project code AS110.

We used the RACS low-band survey as ‘epoch 0’ of the VAST Phase I Pilot Survey (VAST-P1) which was conducted entirely in the lowest frequency band. We will use the RACS mid-band survey as ‘epoch 0’ of the VAST Phase II Pilot Survey (VAST-P2), which will be conducted in both the low and mid-bands. This allows us to take advantage of the RACS tiling pattern and survey design, which had very similar parameters to those required for VAST. An important note is that the publicly-released RACS images have had an additional flux correction applied, as described by McConnell et al. (2020), so any epoch 0 flux densities reported in this work may differ from the (improved) published RACS catalogue (Hale et al., submitted).

In in this paper, we use the following terminology to describe the different components of our survey and others:

1. beam footprint The arrangement of phased-array feed (PAF) beams used to cover a single field. We use the square_6x6 and closepack36 beam footprints (see Hotan et al. 2021).

2. field A single observation or pointing.

3. tiling pattern How individual fields are placed to cover a larger area.

4. region A grouping of pointings (usually but not always covering adjacent) meant to cover a particular area of the sky or type of source.

5. survey footprint The total sky area covered by a survey (for instance, DES) or one epoch of a survey (as in VAST-P1).

2.2. Phase I pilot survey observations

VAST-P1 consists of 113 fields taken from the RACS low-band tiling pattern. The survey coverage of VAST-P1 is shown in Figure 1. The survey consists of six regions, chosen to test some of the different VAST science goals. These are summarised below.

Figure 1. The VAST-P1 survey footprint, showing the number of observations of each field. The green region shows the planned survey footprint of the VAST-P2 mid-band observations. VAST-P2 low-band will cover the same survey footprint as VAST-P1. The sky map is plotted with J2000 equatorial coordinates in the Mollweide projection and the background diffuse Galactic emission at 887.5 MHz is modelled from Zheng et al. (2017) using Price (2016).

Region 1: An extragalactic region of 40 fields ( $1840\,\mathrm{deg}^2$ ) along an equatorial strip centred at RA of $0^{\rm h}$ . This was chosen to overlap with the CNSS, which itself had considerable overlap with FIRST. In addition, this overlaps with part of the DES/DECaLS coverage.

Region 2: A second extragalactic region of 28 fields ( $1300\,\mathrm{deg}^2$ ) along an equatorial strip centred at RA of $12^{\rm h}$ . The equatorial regions have the advantage of greater multi-wavelength coverage from telescopes in the northern hemisphere. This region in particular has good overlap with FIRST and DECaLS.

Region 3: An extragalactic region of 19 fields ( $820\,\mathrm{deg}^2$ ) chosen within the larger DES survey footprint, centred at ${\rm RA}=4.5\,\mathrm{hr}, {\rm Dec}=-50\deg$ .

Table 1. VAST observing parameters. Note the number of epochs for VAST-P2 (as well as the minimum and maximum spacing) are planned estimates, and may change when the survey is conducted. The image rms and total area for VAST-P2 (mid) are estimated from early RACS-mid observations; these estimated values are marked in italics. See Hotan et al. (2021) for details about the beam footprints.

Region 4: A second region of 19 fields ( $830\,\mathrm{deg}^2$ ) within the DES survey footprint, centred at ${\rm RA}=22\,\mathrm{hr}$ , ${\rm Dec}=-50\deg$ . Region 4 also overlaps with the Evolutionary Map of the Universe (EMU; Norris et al. 2011) Phase I Pilot Survey and the Polarisation Sky Survey of the Universe’s Magnetism (POSSUM; Gaensler et al. 2010) Phase I Pilot Survey.

Region 5: A small 5-field ( $265\,\mathrm{deg}^2$ ) region covering the Galactic Centre.

Region 6: A small 2-field ( $135\,\mathrm{deg}^2$ ) region covering the Small and Large Magellanic Clouds.

Note that the region area above is not the number of fields times $31\,\mathrm{deg}^2$ for two reasons. First, the instantaneous field-of-view is typically more than $31\,\mathrm{deg}^2$ as it includes sky sampled by the PAF beams which is outside the region over which the PAF can point (which is the origin of the $31\,\mathrm{deg}^2$ number). The field-of-view is frequency dependentd but can be up to $64\,\mathrm{deg}^2$ at the low-band frequency, including lower-sensitivity wings of the PAF beams. Second, when fields are adjoining the lower-sensitivity wings overlap, improving sensitivity but reducing total sky coverage.

The ASKAP observing specifications in the context of VAST-P1 are summarised in Table 1, and the full set of VAST-P1 observations is summarised in Table 2. Our survey strategy was to observe five epochs of the VAST-P1 survey footprint, with each epoch separated by months, and one pair of epochs separated by a single day. However, there were also extra test observations conducted, and some epochs in which extra fields were observed to accommodate observations where the schedule initialisation failed (or was excessively delayed). Hence, in addition to the six full epochs (epochs 1, 2, 8, 9, 12 and 13 in Table 2), there are a number of epochs that have partial coverage of the VAST-P1 survey footprint. These epochs are labelled with an ‘x’ in their name in Table 2 (3x, 4x, 5x, 6x, 7x, 10x and 11x). The data for epoch 12 (marked with as asterisk in Table 2) only became available after submission of this paper. We have included it in the table for completeness but the data is not included in any of the analysis presented in this paper.

The total combined area of a full VAST-P1 epoch is $5131\,\mathrm{deg}^2$ . Each field was observed for 12 min integration time (rather than 15 min for RACS). Excluding RACS (as epoch 0) we observed 813 fields, for a total observing time of $\sim\!162\,\mathrm{h}$ . In this paper we present results from two of the regions (3 and 4). The rest of the VAST-P1 survey footprint will be analysed in subsequent papers.

2.3. Phase II pilot survey strategy

The VAST Phase II Pilot Survey will commence in 2021. The aims of VAST-P2 are: to build on the scientific analysis from VAST-P1; to test faster turn-around of imaging and transient detection; and to demonstrate commensality with ASKAP mid-band surveys by observing some epochs at 1296 MHz.

Table 2. Summary of VAST-P1 observations, giving the number of fields in each epoch, the start and end dates for the epoch, and the total sky area. Epochs with an ‘x’ in the name only have partial sky coverage, as discussed in the text. Epoch 0 is the RACS survey. Epoch $12^*$ was only available after submission of this paper and is included here for completeness but is not included in any of the analysis.

We plan to observe five epochs spread over the six months of the pilot programme. We will observe two epochs at low-band, using the same survey footprint as VAST-P1 and the square_6x6 beam footprint. The remaining three epochs will be observed in the ASKAP mid-band. This will use the tiling pattern from the RACS mid-band survey and the closepack36 beam footprint (Hotan et al. 2021).

The size of the primary beams scales proportionally with wavelength. Since we tile the sky with beam footprints that provide nearly uniform sensitivity across the observed area, our mid-band survey will require a different tiling pattern where each beam footprint will cover less area than the low-band beam footprints. Just as we have done for our low-band survey, we select our mid-band tiling pattern from the RACS mid-band survey tiling pattern. Note that since the mid-band beam footprints are smaller, the low- and mid-band field centres are not aligned.

Subtracting the duration of two low-band epochs from the total VAST-P2 observing time will leave enough time to observe 3 epochs of 91 mid-band fields. This will not be enough to cover the entire low-band survey area, so we have iteratively selected mid-band fields from RACS to optimise the overlap with the low-band fields from the following VAST regions in order: 5 (Galactic Centre), 6 (Magellanic Clouds), 1 (equatorial strip), then 4 (DES/EMU/POSSUM). This results in 88 mid-band fields. The remaining 3 mid-band fields have been selected from region 3 (DES) as it is adjacent to the already selected region 6 fields that cover the LMC. The full set of selected mid-band fields is shown in Figure 1.

2.4. Data products

VAST-P1 data was processed using the same calibration and imaging strategy as that used by the RACS survey (McConnell et al. 2020), the only difference was that VAST used a standard Gaussian primary beam correction whereas RACS applied a position-dependent correction based on holography data to correct for slight non-Gaussianity of the PAF beams. The use of holographic data to correct for primary beam attenuation was made available for general processing in the ASKAPsoft pipeline after most VAST processing for VAST-P1 was completed, however, it was used in the processing of epoch 12. There are two main image data products:

1. 1. The individual field (or ‘tile’) images produced by ASKAPSoft. These were created by combining the 36 individual beam images. We produced multi-frequency synthesis (MFS) images for all four Stokes parameters I, Q, U and V. Two Taylor-terms are used in Stokes I imaging to consider the spectral slope of field sources across the observed band.

2. 2. Combined images made by mosaicing the individual field images together (these are not produced by ASKAPSoft) across each region of VAST-P1. This means that images at different times within a given epoch are combined, so they can only be used for epoch-to-epoch comparisons.

A typical VAST-P1 image is shown in Figure 2. When imaging, we removed short ( $< 100$ m) baselines in order to; (i) to minimise solar interference; and (ii) as there is insufficient uv-coverage to reliably recover extended structure in short VAST-P1 observations.

Figure 2. A VAST image within Region 4 (field VAST_2131 $-$ 56A) from epoch 12 with two cutouts. Cutout A: a $1\deg$ image centerd on (J2000) $\alpha$ = 21:49:26.5, $\delta$ = –55:17:52.87 containing several bright sources, including the large radio galaxy 2MASX J21512991–5520124. Cutout B: a $0.3\deg$ image centerd on (J2000) $\alpha$ = 21:36:18.9, $\delta$ = –58:00:12.68 containing a range of source morphologies.

A histogram of the rms noise in each individual VAST image in regions 3 and 4 is shown in Figure 3, where we separately show Stokes I and Stokes V. The latter is dominated by thermal noise in the images, while the former also has significantly more contributions from source confusion, sidelobe confusion, and deconvolution artefacts. The median rms is $0.24\,{\rm mJy\,beam}^{-1}$ in Stokes I, while it is $0.20\,{\rm mJy\,beam}^{-1}$ in Stokes V.

Figure 3. Distribution of median image rms values (computed over the central half of each image) for each field in each epoch of regions 3 and 4. We plot the rms values from Stokes I (blue) and Stokes V (orange).

In addition to images, the ASKAPSoft pipeline also produces a catalogue of groups of pixels of contiguous emission (‘islands’), a catalogue of Gaussian components fit to the islands (‘components’, where each island may be fit with one or more components), and an rms map and a median map corresponding to each field. These are created by the Selavy source finding software (Whiting 2012) incorporated into ASKAPSoft. We ran Selavy separately on our combined images, and used these additional data products as inputs to our VAST transient detection pipeline, as described in Section 3. We used the default parameter settings in ASKAPSoft, so sources were considered ‘detected’ if they exceeded 5 times the local rms (roughly 1.2 mJy in Stokes I for an unresolved source).

Two of the VAST-P1 epochs (8 and 9) were observed on consecutive days, and were matched in local sidereal time (LST). This was done to explore whether LST-matching improved the ability to search for transients using image subtraction. The observations for the matched epochs were taken over a period of two weeks, but within that period each individual field was LST matched and observed on consecutive days.

2.5. Data quality analysis

Before images were released in CASDA, we performed some quality control checks of the astrometric accuracy and flux density scale (both relative and absolute). In addition, we visually inspected each image to check the overall data quality (on the basis of this, a single field was excluded: VAST_0534–43A in epoch 7x, which contains Pictor A). In this section, we summarise the results for the two VAST-P1 regions (3 and 4) presented in this paper (see Section 4.2).

2.5.1. Astrometric accuracy

To evaluate the astrometric accuracy of our images, we extracted bright (SNR $>=7$ ), compacte sources using Selavy that are isolated by a minimum of 150 arcsec from other sources, and then crossmatched them with the International Celestial Reference Frame (ICRF) catalogue (Charlot et al. 2020). The ICRF consists of 4536 radio sources with accurate positions measured using very long baseline interferometry (VLBI).

The left panel of Figure 4 shows the offsets for VAST-P1 sources associated with the 41 ICRF sources in regions 3 and 4 (a total of 331 comparisons across all epochs). The median and standard deviation of the positional offsets is $-0.19 \pm 0.53\,\mathrm{arcsec}$ in right ascension and $0.07 \pm 0.48\,\mathrm{arcsec}$ in declination, with a standard error of $29\,\mathrm{milliarcsec}$ in both coordinates. The median positions of each epoch of VAST-P1 are shown as coloured markers.

Figure 4. Astrometric accuracy for compact sources in regions 3 and 4 of VAST-P1 compared to Left: sources from the ICRF catalogue and Right: sources from the RACS catalogue, as described in Section 2.5.1. The image pixel size of $2.5\times 2.5$ arcsec is shown as a red dashed box. The median offset for each individual epoch is shown with coloured markers. The solid lines show the overall median offsets, and the dashed lines are the median $\pm 1$ standard deviation. For the comparison with RACS, the background is a 2D histogram of the source counts where the colour scale represents the number of sources per bin.

We also compared the VAST-P1 sources with the positions in the published RACS catalogue (Hale et al., submitted). There were $378\,823$ unique measurements of $32\,769$ compact, isolated RACS sources across all epochs within regions 3 and 4 with SNR $>=7$ . The results of this comparison are shown in the right panel of Figure 4. The median and standard deviation of the positional offsets are $0.58\,\ \pm\, 1.01\,\mathrm{arcsec}$ in right ascension and $0.14 \pm 0.99\,\mathrm{arcsec}$ in declination, with a standard error of $1.6\,\mathrm{milliarcsec}$ in both coordinates.

2.5.2. Flux density scale

In addition to the absolute flux density scale, for variability analysis it is important that the flux density scale from epoch to epoch is consistent. To evaluate both of these, we crossmatched compact, isolated sources with $\mathrm{SNR}>=7$ in each epoch of VAST-P1 with sources in the published RACS catalogue. Figure 5 shows the flux density ratio for these $378\,823$ sources.

Figure 5. Absolute and relative flux density scale comparison for $378\,823$ sources in regions 3 and 4, comparing those measured in each VAST-P1 epoch to RACS (as discussed in Section 2.5.2). The median ratio for each epoch is shown by the dotted lines. The background is a 2D histogram of the source counts where the colour scale represents the number of sources per bin. The overall median ratio is 0.98 with a scatter of 0.15.

Nine of the twelve epochs have a median VAST-P1/RACS flux density ratio between 0.97 and 1.03 (within 3%). Epoch 10x has a ratio of 0.93 and Epochs 05x and 06x have ratios of 0.94 and 0.95 respectively. The overall ratio is 0.98 with a $1\sigma$ scatter of 0.15. This scatter is probably the result of a few differences in the way that RACS and VAST-P1 data were processed. The most significant of these is that the published RACS images had a holography correction applied that improved the flux density scale, particularly in the outer regions of images. These discrepancies will be corrected in future analyses, but for the current work our conservative variability thresholds (see below) mean that they are not significant.

2.5.3. Known issues

As discussed above, the observations presented here are part of a pilot study for a full survey. Therefore some of the observations and analysis used were less than optimal, and evolved as we discovered deficiencies. These lessons will inform the full VAST survey. Nonetheless, the data are useful for scientific analysis, but we note the various issues below.

From epoch 4x onwards, we made all observations within $\pm 1\,\mathrm{h}$ of meridian transit to ensure a consistent synthesised beam size for all fields. In epochs 1–3x, this did not occur and some fields were observed at large zenith angles, resulting in very extended synthesised beams.

From epoch 5x onwards, an intermittent system issue at the initialisation of a new field observation caused random phase delays on a few antennas on a small subset of field observations. These were identified after processing, but could not be corrected in all cases. Furthermore, even for corrected fields, this occasionally resulted in increased astrometric and flux scale errors. Fields that were affected by this issue are highlighted in the release notes in CASDA. Other more subtle changes in the bandpass and beam-shape may also have occurred in early observations as beam-forming and use of the on-dish calibrator improved over the course of the pilot survey. As a result, more recent data tend to be of higher quality than earlier data.

The process used to obtain an initial deconvolution model for use in self-calibration in the ASKAP pipeline can introduce a slight positional offset in each beam image. This is because self-calibration will tend to snap to a dominant field source and place it at a pixel centre. This not only causes astrometric errors, up to half a pixel (1.25 arcsec) in extent, but can also cause apparent extension of a source if the shift is in opposite directions for a source in two overlapping beams.

The source finding for the VAST pipeline input was performed on final mosaiced images. In early pilot survey images, it was assumed that the point spread function (PSF) of each beam would be consistent in each field. In practice, the PSF can vary from beam to beam as a result of differences in flagging and elevation. As only a single PSF is recorded in the image header, typically the PSF associated with beam 0, this can affect the accuracy of the flux scale across the entire mosaic.

Also, note that the RACS flux density scale correction (as described in Appendix A of McConnell et al. 2020) has not been applied to the VAST data. Hence we also used the pre–flux-corrected RACS data in epoch 0, for consistency. For future work, we plan to convolve all fields to a consistent PSF and apply these corrections to ensure consistency between RACS-low, RACS-mid and both the VAST pilot surveys.

2.6. Data availability

The raw data and individual field images for the entire VAST-P1 survey are publicly available through CASDA under project code AS110 for RACS and AS107 for VAST. Epoch 4x of VAST-P1 was conducted as a test observation and is available under project code AS113. Higher order data products such as variability catalogues will be released once full analysis is complete.

3. Transient detection pipeline

The primary technical goal of the VAST Pilot Surveys is to demonstrate the capability to detect highly variable and transient sources. One of the main open-source software packages for image-based radio transient detection is the LOFAR Transients Pipeline (TraP; Swinbank et al. 2015). TraP can be applied to radio image data from any facility, and has been used successfully in a number of transient projects (e.g., Stewart et al. 2016; Driessen et al. 2020; Sarbadhicary et al. 2020).

Testing TraP on the early ASKAP datasets revealed some scalability issues. After investigating a number of options, we decided to implement a new pipeline to search ASKAP data in an efficient way. Most of the key logic of this pipeline follows the approach of TraP, but with some implementation differences that significantly improve the performance.

The main VAST Pipeline technology stack is: Python 3.7+, PostgreSQL f $12+$ , Astropy g $4+$ , Django h $3+$ , Pandas i $1.2+$ , Dask j $2+$ and Bootstrap k 4. A more complete technical description of the pipeline implementation is given by Pintaldi et al. (2021). In this section we discuss the scientific functionality and resulting data products in the context of the results presented here.

3.1. Pipeline processing method

The VAST pipeline takes as input a set of images, rms maps and mean (background) maps (all as FITS files) from Selavy as well as the Selavy source finding component list. These data products are all default outputs from the ASKAPSoft continuum imaging process, which means the VAST pipeline can be run immediately after the ASKAPSoft data reduction is completed.

For the work presented in this paper, we used the combined images (created by mosaicing individual fields together) as the input. This means that only epoch-to-epoch timescales could be explored. The entire Selavy source list was read in, but any duplicate sources from overlapping fields within the same epoch were filtered out using a tight crossmatching radius of 2.5 arcsec. In principle it would also be possible to search for variability within epochs in the regions where fields overlap: this will be the subject of future analysis. Investigating faster variability (within a single observation, for example, Wang et al. 2021) will also be explored in future work.

We ran the VAST pipeline on these combined images, following the steps below:

1. 1. Image ingest: Image metadata and Selavy catalogues were ingested for each image.

2. 2. Uncertainties on flux densities and positions for every source component were calculated following Condon (1997).l

3. 3. Source association was done by crossmatching each epoch to the next epoch using the de Ruiter method (Scheers (2011), based on de Ruiter et al. 1977). We used the default beamwidth limit of 1.5 and the search radius of 5.68 (these are both unitless values). The beam size that is used for the epoch mode association should be equal to the largest beam contained in the epoch images.

4. 4. Forced extractions (where a flux density is measured at a specified position) were performed for sources that were not detected in every epochm, so that a complete light curve could be built for each source. Where there were multiple observations of a given source within an epoch (e.g., due to image overlap), the forced extractions were taken from the image whose centre was nearest to the source position.

5. 5. New sources that had not appeared in previous epochs were identified.

6. 6. Variability metrics were calculated for all sources identified by the pipeline (see Section 3.2). The forced extraction measurements were incorporated when calculating all variability metrics.

7. 7. Results were saved in a database and output to parquet n files, an open source columnar data format that allows flexible custom post-processing (Vohra 2016).

The resulting data products are a light curve for every source detected in one or more epochs, and a set of standard variability metrics. Figure 8 shows light curves for the highly variable sources presented in this paper.

The pipeline results are available to users via either an interactive web application (including some basic multi-wavelength information), or through the parquet files that can be read by a variety of data-analysis software packages. The VAST Tools python package, described in Section 3.3, offers a streamlined method to explore the pipeline results in a Jupyter notebook environment. For more detailed information about the VAST pipeline, see the publicly available software and documentationo.

3.2. Variability metrics

The VAST pipeline identifies highly variable and transient sources using the two key light curve parameters identified by Rowlinson et al. (2019). These are: the modulation index V, a measure of the variability in a light curve; and the reduced $\chi^2$ relative to a constant model $\eta$ , a measure of the significance of the variability. These parameters are defined as:

(1) \begin{equation}V = \frac{1}{\overline{S}}\sqrt{\frac{N}{N-1} (\overline{S^2} - \overline{S}^2)},\end{equation}

(2) \begin{equation}\eta = \frac{N}{N-1}\left( \overline{w S^2} - \frac{\overline{w S}^2}{\overline{w}}\right),\end{equation}

where N is the number of data points, the flux density and uncertainty on the i-th epoch are $S_i$ and $\sigma_i$ , means are denoted by overbars (i.e., $\overline{S} \equiv \frac{1}{N}\sum_i S_i$ ) and in calculating $\eta$ the flux density measurements are weighted by the uncertainties according to $w_i=1/\sigma_i^2$ (also see Swinbank et al. 2015).

The parameter V selects sources with high variability without assessing the statistical significance of that variability. In contrast, the parameter $\eta$ acts as a measure of the statistical significance of that variability without knowledge of its absolute amplitude. Used in combination, these parameters select a sample of statistically significant, highly variable objects.

The pipeline also allows for two-epoch variability searches (where the variability is calculated based on pairs of observations), but this functionality was not used in the results presented here.

3.3. VAST tools software

VAST Tools p is a python package containing tools for accessing and exploring VAST data products. The current functionality at the time of latest release (Version 2.0.0) is outlined below. VAST Tools can be used via a command line interface or in Jupyter notebooks. VAST Tools can also be used as an interface to analyse results from a VAST Pipeline run.

3.3.1. Survey information

The VAST Tools package includes the complete details of the spatial and temporal coverage of VAST-P1. The beam centres for all fields are listed in comma-separated-values (CSV) files on a per-epoch basis, while Hierarchical Equal Area isoLatitude Pixelisation (HEALPix; Górski et al. 1999) Multi-Order Coverage (MOCq; Fernique et al. 2019) maps contain the spatial coverage of each survey field, and each epoch. We have also created a space-time MOC for the full pilot survey coverage. The MOCs can be accessed though a mocpy r wrapper included in the VAST Tools package, which also allows users to load MOCs and also query VizieR tables for sources within the MOC coverage.

3.3.2. Querying coordinates

The Query class allows users to search the VAST pilot for observations of a given location. Three main query types exist: (1) finding the fields containing the specified coordinates; (2) searching for radio emission at the specified coordinates; and (3) finding all sources within a given radius of the specified coordinates. Query types 2 and 3 require VAST tools to be run on a machine with access to the relevant FITS files and Selavy catalogues.

VAST tools can crossmatch coordinates with the Selavy catalogues, and produce light curves, FITS cutouts and postage stamp images with various overlays (e.g., Figure 9). Forced fitting is available, and there are options to produce light curves with upper limits (based on Selavy rms maps), forced fits for non-detections, or forced fits for all epochs.

3.3.3. Pipeline interface

Users are able to explore the results from the VAST Pipeline in an interactive manner using the Pipeline class, which contains built-in methods to perform tasks such as transient analyses, plotting lightcurves and querying external resources. This is achieved by using the parquet files that are provided as output by the pipeline, and using data exploration libraries such as Pandas and vaex s.

4. Untargeted variability search

Some of the early results from VAST focus on the time-domain behaviour of samples of known objects selected in radio or other wavebands, for example the search for radio afterglows from Swift GRBs that resulted in a detection of the late-time afterglow from GRB 171205A (Leung et al. 2021).

In contrast, in this section we present initial results from an untargeted search for highly variable sources and transients in regions 3 and 4 of the VAST Phase I Pilot Survey. This demonstrates the functionality of the VAST pipeline, and shows some of the diversity of variable sources that the full VAST survey will detect. Analysis of the other survey regions will be presented in subsequent papers.

4.1. Data analysis

We ran the VAST pipeline on the fields from regions 3 and 4, following the process outlined in Section 3. We then selected compact, isolated sources by excluding sources that satisfied any of the following criteria:

1. • fewer than two measurements (forced or Selavy);

2. • neighbouring sources within 30 arcseconds;

3. • extended (a ratio of average integrated flux to average peak flux greater than 1.5);

4. • a SNR of $<\!7$ ;

5. • negative modulation index (caused by a negative mean flux density).

After applying these criteria, our remaining sample of compact, isolated sources consisted of 155 071 unique sources each detected in at least one epoch. This is a very conservative selection. For instance, requiring no nearby neighbours helps to eliminate imaging artefacts but will also exclude sources in a galaxy with a detectable nucleus or sources in crowded regions. Therefore the results below should be taken as a lower limit, and future searches should identify more objects even in the same regions of the sky.

Figure 6 shows the distribution of sources in the two key variability metrics V and $\eta$ (defined in Equations 1 and 2). Using $2\sigma$ cutoffs in both of these parameters ( $V>0.51$ , $\eta > 5.53$ ) we identified 171 sources as highly variable. These appear in the shaded top right-hand quadrant of Figure 6.

Figure 6. A plot of the two key variability metrics, V and $\eta$ . Sources that appear in the shaded top-right quadrant are variable candidates as they exceed the $2\sigma$ thresholds on V and $\eta$ calculated by fitting a Gaussian function to the sigma-clipped distributions of each metric. The thresholds are $V > 0.51$ and $\eta > 5.53$ . Sources that have been classified as variables after manual inspection are marked in red.

We manually inspected the radio light curves, images and multi-wavelength data for each of these sources, and identified 108 as artefacts near bright sources. Another 14 were in regions of poor data quality, and 21 were detections of marginal significance. This left 28 reliable astronomical transients and variables, which are marked as red in Figure 6. We discuss these variable sources in the next section.

Figure 7. Predicted modulation index (for interstellar scintillation) V versus 888-MHz flux density for pulsars. We show all of the pulsars in the ATNF pulsar catalogue (Manchester et al. 2016) as grey points, with the flux densities computed from the catalogued values at 1.4 GHz or 430 MHz assuming a spectral index of $-1.6$ (Jankowski et al. 2018). Predicted modulation index is calculated using the Galactic electron-density model of Yao et al. (2017) and diffractive scintillation following Lorimer & Kramer (2012) and Cordes & Lazio (1991). The pulsars contained but not detected in VAST-P1 regions 3 and 4 are black outlined circles. The highly variable pulsars identified here are blue circles (based on catalogue values, with the addition of Corongiu et al. 2021 for PSR J2039 $-$ 5617 and Bhattacharyya et al. 2019 for PSR J2144 $-$ 5237) that are connected to their observed properties in VAST-P1 (blue squares). Rather than mean flux density we show the maximum flux density, as that governs detectability. Our thresholds of flux density $\geq\!1.2\,$ mJy (based on the median image rms and a $5\sigma$ detection threshold) and modulation index $\geq\!0.51$ are the dashed lines. The black open circle at the right near $300\,$ mJy is PSR J0437 $-$ 4715: it was not identified in our sample because it is below our V threshold (it has $V=0.4$ ) and it has an unrelated source within the 30 arcsec neighbour limit.

Figure 8. Lightcurves of the variables identified in Regions 3 and 4 of VAST-P1. Blue round points are peak flux density measurements from Selavy. Grey triangles are the $3\sigma$ upper-limits of the forced-fitted flux density for images where there was no Selavy detection.

Figure 9. Images of the variable sources identified in this paper. The left panel shows the VAST Stokes I image for the epoch with the maximum flux density. The middle panel shows the Stokes V image for the same epoch where positive flux density corresponds to right handed circular polarisation and negative to left handed. The ellipse in the lower left corner of each radio image shows the FWHM of the restoring beam. The right panels show Stokes I contours at 30%, 60%, and 90% of the peak Stokes I flux density overlaid on an RGB image of optical data from either DES DR1 (DES) or SkyMapper (SM; Onken et al. 2019) with red=i-band, green=r-band, blue=g-band. All images have been centred on a frame aligned with the position of the radio source. Optical data have been astrometrically corrected to the VAST epoch according to the proper motion of the target when available.

4.2. Results

Our 28 highly variable and transient sources are listed in Table 4 with light curves in Figure 8 and image cutouts in Figure 9. Seven of these sources are known pulsars and seven are radio stars. We discuss each of these classes below. The remaining sources are discussed in Section 4.2.3.

4.2.1. Pulsars

Seven of our highly variable sources are known pulsars: PSRs J0255 $-$ 5304, J0418 $-$ 4154, J0600 $-$ 5756, J2039 $-$ 5617, J2144 $-$ 5237, J2155 $-$ 5641, and J2236 $-$ 5527. Pulsars can be variable for intrinsic reasons (nulling, intermittency, e.g., Backer 1970; Kramer et al. 2006; Wang et al. 2020), external reasons (interstellar scintillation; e.g., Rickett 1970; Stinebring & Condon 1990; Bell et al. 2016; Kumamoto et al. 2021), or because of absorption by gas in binary systems (e.g., Broderick et al. 2016; Polzin et al. 2020; Kudale et al. 2020). In fact this variability can help separate pulsars from the background of continuum sources (e.g., Dai et al. 2017).

Three of the sources are in binary systems: PSRs J2039 $-$ 5617, J2144 $-$ 5237, and J2236 $-$ 5527. Of these, PSR J2039 $-$ 5617 does have eclipses that cause radio pulsations to disappear for half the orbit (Corongiu et al. 2021). That source also has significant variability in the pulsed amplitude, which Corongiu et al. (2021) largely ascribe to scintillation but also argue could be caused by more ‘exotic’ mechanisms (such as absorption by intra-binary gas). PSR J2144 $-$ 5237 is a binary with a low-mass companion (Bhattacharyya et al. 2019) that, based on its mass alone, could be a ‘redback’ system. Such systems are often eclipsing (Roberts 2011), but there is no sign of eclipses in the timing data and so a helium-core white dwarf seems more likely. PSR J2236 $-$ 5527 is similarly presumed to have a helium-core white dwarf companion (Burgay et al. 2013). Therefore binary companions do not seem to be causing the variability in these systems, although PSR J2039 $-$ 5617 may merit further analysis.

For the other four pulsars, there are no exotic phenomena noted in the literature (e.g., Newton et al. 1981; Xie et al. 2019; Johnston et al. 2021). Therefore, in the absence of any other reason we would expect these detections to be due to variability via refractive interstellar scintillation on longer timescales (e.g., Kumamoto et al. 2021), plus diffractive interstellar scintillation on shorter timescales (Rickett 1990; Cordes & Rickett 1998) although since not all pulsars have been observed extensively we cannot rule out other effects. To examine this in more detail, in Figure 7 we plot for all known pulsars the predicted modulation index for interstellar scintillation versus the flux density, and highlight those identified here. The predictions are far from perfect. We have assumed a single spectral index for our flux density estimates, and a rough electron density model (there is more than an order of magnitude variation in scattering timescale as a function of dispersion measure; Kumamoto et al. 2021) along with a single source velocity for our scintillation estimates. However, the sources identified through variability are among those expected to be brighter and more heavily modulated. We also highlight the difference between the predicted properties (which will reflect the mean flux densities) and the observed properties, which will be biased to high flux densities and modulation levels because of the nature of these detections.

The maximum dispersion measure (DM) in regions 3 and 4 is about $30\,{\rm pc\,cm}^{-3}$ (Yao et al. 2017), which is also consistent with the DMs of the pulsars we identified (they span 14– $30\,{\rm pc\,cm}^{-3}$ ). As the scintillation timescale and bandwidth are largely dependent on the DM (e.g., Kumamoto et al. 2021), the DM range can be used as a proxy to help identify which other sources will be detected as variable in the remainder of our survey.

4.2.2. Stars

Seven of our highly variable sources are identified with known stars. Two of them (CD–44 1173 and UPM J0409–4435) were found in the circular polarisation search conducted by Pritchard et al. (2021), and are discussed there. Our new detections (AB Dor A/C, SCR J0533–4257, UCAC3 89–412162, LEHPM 2-783 and 2MASS J22414436 $-$ 6119311) are discussed below and the images are shown in Figure 9.

VAST J052644.9 $-$ 652650 is identified as AB Dor, a quadruple system composed of a K1V pre-main sequence star + M5.5V and M3.5 + M4.5 pair of binary systems (Close et al. 2007). The pre-main sequence star AB Dor A, is by far the brightest star in the system; it has a rotation period of 0.514 d, and shows optical super-flares (Schmitt et al. 2019) and long-term X-ray variability (Lalitha & Schmitt 2013). VLBI radio observations at 1.4, 8.4, and 22.3 have revealed milliarcsecond radio structure (Climent et al. 2020). The source is detected in VAST-P1 in 11 epochs in Stokes I with a maximum flux density of $7.8\,\mathrm{mJy\,beam}^{-1}$ , and in three epochs in Stokes V showing both right- and left-handed circular polarisation with circular polarisation fraction ranging from 52% to 60%. We speculate that the three events with high circular polarisation are associated with known prominence formation in young active stars such as AB Dor, HK Aqr, BO Mic and PZ Tel (Doyle & Collier Cameron 1990; Barnes et al. 2001; Jardine et al. 2001). AB Dor lies in the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) southern continuous viewing zone, allowing near continuous optical photometry at 2 min cadence extending over a year (Ioannidis & Schmitt 2020).

VAST J053328.0 $-$ 425719 is identified as SCR J0533–4257, an M4.5 dwarf (Riaz et al. 2006). Ultraviolet activity (Schneider & Shkolnik 2018) and X-ray flares (Fuhrmeister & Schmitt 2003) have been observed from this star, though there are no previous radio detections reported in the literature. The source is detected in three epochs. The maximum flux density is $2.83\pm0.24\,\mathrm{mJy\,beam}^{-1}$ . The circular polarisation fraction ranges between 33% and 57%. It was observed using TESS (TIC 302962949) in 2 min cadence during sectors 5 and 6 and in 20 sec cadence in sector 33 (total duration 67.9 d). The TESS data show evidence for modulations at periods of 1.25 and 0.46 d, hinting it may be a binary. Optical flares are frequently seen, most of which have a peak amplitude of a few $\sim\!0.01\,\mathrm{mag}$ but up to 0.2 mag on one occasion.

VAST J201949.8 $-$ 581650 is identified as LEHPM 2–783, an M6 dwarf (Riaz et al. 2006) with high proper motion ( $\mu_{\alpha}=-23.4$ , $\mu_{\delta}=-341.1\ \mathrm{mas\ yr}^{-1}$ ) (Gaia Collaboration et al. 2018). There are no previous radio detections reported in the literature. The source is only detected in a single epoch with a flux density of $2.50\pm0.27\,\mathrm{mJy\,beam}^{-1}$ and a circular polarisation fraction of 77%. It was observed using TESS (TIC 387220832) in 2 min cadence during sectors 13 and 27 (total duration 50.4 d). The TESS sector 13 data shows evidence for a period at 6.7 d whilst the sector 27 data suggest a longer period of 11.6 d. It is flare-active, with several optical flares having amplitudes greater than 0.5 mag. It has a location on the Gaia colour-magnitude diagram red-ward of the main sequence, indicating it maybe a binary or a young star. The high circulation polarisation of 77% could point towards auroral magnetospheric activity, similar to that seen by Zic et al. (2019) in UV Ceti, a star of similar spectral type.

VAST J212217.5 $-$ 454631 is identified as UCAC3 89–412162, which has a distance of 31.6 pc and with a Gaia colour $(\mathrm{BP}-\mathrm{RP}) =2.49$ and absolute magnitude $\mathrm{MG} =10.0$ (Gaia Collaboration et al. 2018) is consistent with an M3V spectral type (Pecaut & Mamajek 2013). There are no previous radio detections reported in the literature. This source is also detected in a single epoch with a flux density of $2.31\pm0.25\,\mathrm{mJy\,beam}^{-1}$ and a circular polarisation fraction of 87%. It was observed using TESS (TIC 200645542) in 2 min cadence during sectors 27 and 28 (total duration 41.1 d). There is no evidence for a periodic modulation in the optical light curve or flaring activity. The lack of flaring activity in this object is consistent with reduced or minimum spot coverage. This of course does not mean a lack of flaring, e.g., Schmitt et al. (2019) suggested that large flares only have $\sim\!14\%$ of their energy in the TESS passband. The fact that we detect a radio outburst could suggest that this object may have flares that have a peak flux in the blue, since the contrast over the background emission is much higher at shorter wavelengths and the TESS wavelength band in centred at 7865 Å.

VAST J224144.7 $-$ 611933 is associated with 2MASS J22414436 $-$ 6119311. It has a high proper motion ( $\mu_{\alpha}=150.2$ , $\mu_{\delta}=-87.9\ \mathrm{mas\ yr}^{-1}$ , from Gaia), a colour ( $\mathrm{BP}-\mathrm{RP}=2.81$ ) consistent with a M3.5V spectral type and has a distance of $28.42\pm 0.06\ \mathrm{pc}$ (Gaia Collaboration et al. 2021). It was observed in TESS sectors 1, 27 and 28 (TIC 232064182) in 2 min cadence and shows evidence for a periodic modulation on a period of 0.721 d. Short duration optical flares are seen with amplitudes typically of a few percent, but on one occasion with an amplitude of $\sim\!60\%$ .

All of the stars in our sample show radio emission with high fractional circular polarisation ( $>\!30\%$ ), making this a useful tool for distinguishing these sources from extragalactic transients. They have optical magnitudes between V=17.17 and V=10.9, and span a wide range in brightness and colour. With the exception of SCR J0533 $-$ 4257, these stars have high proper motions ( $>\!100$ mas/yr) suggesting that in the absence of an optical identification, VLBI proper motion measurements are another way in which they can be distinguished from extragalactic sources.

4.2.3. Other variables

We found 14 highly variable sources that were not associated with a known pulsar or star. For these we searched for optical or infrared counterparts in several WISE catalogs, Gaia EDR3 (Gaia Collaboration et al. 2021) and the Dark Energy Survey DR1. Proper motions from Gaia were applied if available. Photometric redshifts were retrieved from the Photometric Redshifts for the Legacy Surveys catalogue (Zhou et al. 2021) when available. For sources that have unclear WISE counterparts due to resolution limits, we used forced-photometry measurements from the unWISE images (Lang 2014; Meisner et al. 2016; Meisner et al. 2017, for W1 and W2) and WISE images (for W3) provided in the DESI Legacy Imaging Surveys DR8 catalogue (Dey et al. 2019).

None of these 14 sources have a detection in circular polarisation within the VAST-P1 data, although many of the limits are not very constraining, with $3\sigma$ Stokes V/I fractional polarisation limits between 11% and 56%. The multi-wavelength properties of these objects are summarised in Table 3. Sources with a WISE counterpart (including some of the stars discussed in Section 4.2.2) are plotted on a WISE colour-colour plot in Figure 10.

Table 3. Multi-wavelength properties of the 14 variable sources not identified as known stars or pulsars; see Section 4.2.3. We provide the WISE cross-ID, or for sources with unclear WISE counterparts, the unique LS object ID from the DESI Legacy Imaging Surveys DR8 catalogue. The median of the photometric redshift probability distributions are provided when available from Zhou et al. (2021) along with the limits of the 95% confidence interval.

Table 4. Highly variable sources identified in the VAST-P1 regions 3 and 4. The coordinate of each source is given as the weighted average of all Selavy detections, where the weight is the inverse square of the positional error. $\sigma_{\text{pos}}$ is the averaged positional uncertainty. $\eta$ and V are the variability parameters described in the text. nE gives the number of epochs (observations) that cover the source location. nD gives the number of detections. $|{\rm V}|/{\rm I}$ is the ratio of Stokes V to Stokes I flux density measured in the epoch for which this is a maximum, or the most constraining 3 sigma upper limit in the case of non-detections in Stokes V.

Figure 10. WISE colours of the variable sources identified in this paper plotted on top of the classification regions from Wright et al. (2010). Sources classified as stars are shown in orange, and those classified as other variables are shown in blue. Sources that did not have a clear counterpart in WISE but had unWISE forced-photometry available in the DESI Legacy Imaging Surveys DR8 catalogue are shown with open markers.

Two of the objects are associated with AGN. VAST J031648.1–625825 is coincident with the IR source WISEA J031648.22–625825.5 (and Gaia source EDR3 472209598679248729). The WISE colours of this object (see Figure 10) suggest it is a QSO with a photometric redshiftt of $z = 1.16\:(0.82, 1.50)$ . There are only two detections in VAST-P1. VAST J212623.6-461124 is coincident with the Gaia source EDR3 657552493938965580, and unWISE colours suggest this is a Seyfert galaxy. There are four detections of this source in VAST-P1. Both of these Gaia sources have parallaxes and proper motions consistent with zero. In the absence of further information it is not possible to determine the cause of the observed variability, but it can likely be explained by scintillation or intrinsic AGN variability.

Six of the sources are associated with galaxies, based on their WISE or DES DR1 identifications. None of these are offset from the central positions of their optical or infrared counterparts and so could be due to nuclear emission, but the limited angular resolution means we cannot rule out other sorts of transients (supernovae, GRBs). Therefore, we also searched for possible associations with known transients at other wavelengths. We checked against catalogues of GRBs detected by the International Gamma-Ray Astrophysics Laboratory (INTEGRAL; Mereghetti et al. 2003), the Swift Burst Alert Telescope (BAT; Lien et al. 2016) and the Fermi Large Area Telescope (LAT; Ajello et al. 2019) and Gamma-ray Burst Monitor (GBM; von Kienlin et al. 2020). Possible associations with known supernovae were checked against the Open Supernova Catalogu (Guillochon et al. 2017) and supernovae reported in the Transient Name Serverv. None of the sources are associated with known GRBs or supernovae. Further monitoring during VAST-P2 may provide more insight into the origin of the variability in these objects.

The remaining six sources have no obvious optical/IR counterparts. None are detected in circular polarisation (although again the limiting ratios are not very constraining, with $3\sigma$ Stokes V/I ratios ranging between 24% and 56%). Again, none of these sources are associated with known GRBs or supernovae. With the exception of VAST J204421.6 $-$ 495935, none exhibit a light curve with a temporal rise or decay that is physically consistent with an uncatalogued GRB afterglow or radio supernova. Further observations of these sources are ongoing.

We used the lightcurves and redshifts (where available) to compute brightness temperatures for all of the sources in Table 3 (as in Pietka et al. 2015, Stewart et al. 2016, Bhandari et al. 2018, Bell et al. 2019, and others). For a timescale, we took the minimum timescale over which the two-point modulation index (e.g., Mooley et al. 2016) exceeded 0.5. Given the sparsity of our temporal sampling, these timescales do not have good resolution, but range from $\sim\!1-2\,\mathrm{d}$ (e.g., VAST J051452.1 $-$ 440838) to 266 d (VAST J212857.7 $-$ 395953). Using the photometric redshifts (or $z=1$ when not available) we found brightness temperatures of $10^{13}\,$ K for the most slowly-varying sources to $>10^{18}\,$ K for the most quickly-varying. However, we do not believe those timescales to be intrinsic to the sources, as they would violate the Compton catastrophe limit (Kellermann & Pauliny-Toth 1969) or the stricter equipartition brightness temperature limit (Readhead 1994). Instead we look at three possibilities: extrinsic variability, Doppler boosting, or coherent emission.

For extrinsic variability, we used NE2001 (Cordes & Lazio 2002) to calculate the scintillation transition frequency, $\nu_0$ for extragalactic sources along the lines of sight contained by regions 3 and 4. Both regions have similar values with the values in region 3 spanning $7.8 {-} 12.3\,\mathrm{GHz}$ and the values in region 4 spanning $8.1 - 14.1\,\mathrm{GHz}$ with median values of 9.2 GHz and 9.8 GHz respectively. This quantity can then be used to estimate the expected level of variability caused by refractive scintillation using the scaling relations of Walker (1998). We find that compact sources in both regions will exhibit $\sim\!25\%$ variability at 888 MHz due to refractive scintillation on timescales of $\sim\!10$ days (the all-sky model for refractive scintillation from (Hancock et al. 2019) predicts a similar timescale, but a mean modulation index about a factor of two higher than that predicted from NE2001). In most cases this is reasonably consistent with what we observe, although as we are selecting the most variable sources we are likely probing the tail of the modulation-index distribution.

However, there are sources that seem inconsistent with this interpretation. For instance, VAST J051452.1 $-$ 440838 and VAST J212758.6 $-$ 470528, which vary by a factor of $>2$ on timescales as fast as 1 d. If such variability were to be intrinsic and not violate the Compton catastrophe limit, then both sources would be nearby ( $\lesssim$ 5 Mpc); one may have then expected that a counterpart should be detectable at other wavelengths. Instead, in these cases, diffractive scintillation (which has higher modulation indices and faster variability; e.g., Walker 1998 and (Macquart & de Bruyn 2006)) may be operating if the source size is very compact. Or, we may be seeing a relativistic jet beamed towards us with a Doppler factor of $\sim\!100$ (e.g., Readhead 1994; Lähteenmäki & Valtaoja 1999). Finally, some of the the sources without optical counterparts may be coherently-emitting Galactic objects such as pulsars, where diffractive scintillation is expected and intrinsic variability can be very fast. For most of regions 3 and 4 the most sensitive pulsar survey is the high-latitude portion of the High Time Resolution Universe (HTRU) survey (Keith et al. 2010), which is sensitive to flux densities as low as 0.2 mJy at 1.4 GHz for low DMs. This is a factor of 10 below the flux densities of our sources, but especially if they are scintillating strongly or otherwise varying, previous surveys may have missed them. Further observations to probe the spectral energy distribution, small-scale structure, and precise position of the varying source will help distinguish between these scenarios, as can optical spectroscopy to search for signatures of AGN (in those classified as ‘Galaxy’ in Table 3) and pulsar searches (for those classified as ‘?’).

Figure 11. Two-epoch transient source surface density limits. Black wedges denote upper limits from previous searches for transients on week-month timescales (Bower & Saul 2011; Dobie et al. 2019; Thyagarajan et al. 2011), while black markers show rates from searches with detected transients (Bannister et al. 2011a; Anderson et al. 2020) with associated two-sided 95% confidence interval (Gehrels 1986). Cyan wedges show the phase space that the proposed VAST Deep single field (VAST-DS), Deep multi-field (VAST-DM) and VAST wide surveys will be sensitive to (note that these predictions come from Murphy et al. (2013) and will change once the full VAST survey specifications are confirmed). The result of this work is shown in red. The predicted rates of neutron star mergers (NSM), magnetars, long and short gamma ray bursts (LGRB, SGRB) and tidal disruption events (TDE) from Metzger et al. (2015) are shown with dashed lines.

4.3. Variability analysis and transient source density

We found 28 highly variable or transient sources out of 155 071 compact sources detected in at least one epoch of the VAST-P1 survey. This is a lower limit as we used relatively strict variability thresholds (Figure 6) and also exclude any source with neighbours within $30\,$ arcsec (Section 4.1). In particular, about 2% of the surveyed area was not searched for ‘transients’ because it was within 30 arcsec of an existing source, but as many as 8% of sources were not searched for variability as they had neighbours. Future analyses will relax these requirements and so should discover more transient and variable sources. The variability metrics we used in this analysis correspond to a variability of $\sim\!50\%$ . This implies that only an extremely small percentage ( $0.02\%$ ) of sources above $\sim\!1.75$ mJy are variable at this level on timescales of a few months. This is consistent with (although much lower than) upper limits from previous studies at 1.4 GHz that showed that the fraction of variables on timescales between minutes and years, and flux densities above 0.1 mJy is less than $\sim\!1\%$ (Bannister et al. 2011a; Bannister et al. 2011b; Thyagarajan et al. 2011; Mooley et al. 2013) as summarised by Mooley et al. (2016).

Of our 28 highly variable sources, nine were detected in only a single epoch of our survey (see the sources with ${\rm nD} = 1$ in Table 4) and hence would have been considered ‘transients’ in previous works that calculated transient source densities and rates (e.g., Mooley et al. 2016). Of these nine single-epoch transients, seven are known pulsars or stars, suggesting that a majority of candidates detected in searches for extragalactic transients at this sensitivity will be foreground Galactic sources, even at Galactic latitudes of $|b|\approx 45\deg$ sampled here. These are easy to identify through cross-matching with existing catalogues, or, in some cases by their circularly polarised emission.

The other 19 sources were detected in multiple (2–7) epochs, but some of them could still have been considered ‘transients’ in previous surveys due to the range of temporal scales probed. For instance VAST J041908.3 $-$ 472233 has 4 detections (Figure 8) but they all occur within a 15 d period, with constraining non-detections before and after. So for a typical synchrotron transient with timescale $>30\,$ d at these frequencies (e.g., Metzger et al. 2015), it could represent a single event. Even some sources that appear multiple times could be considered ‘transients’ in the context of a 2-epoch transient rate. VAST J031648.1 $-$ 625825 was visible, disappeared, and then reappeared 480 d later. For a 2-epoch survey metric it could then represent 2 independent transients.

For the purposes of comparison with previous work, we have calculated the two-epoch transient source density for extragalactic synchrotron sources, on timescales of $>30\,$ d (Metzger et al. 2015; Carbone et al. 2017). Future works will explore longer and shorter timescales as well as more detailed rates calculations. We considered only the unclassified sources from Table 4 and computed two numbers of ‘transients’: total source pairs based on a constraining non-detection (upper limit lower than the lowest detection) combined with a detection, and unique sources, where in both cases all observations had been binned to 30 d samples. The former more accurately represents a 2-epoch snapshot rate, while the latter accounts for multiple detections of individual objects. We found 20 total ‘transients‘ comprising 6 unique sources. The equivalent sky area probed was $41\,034.6\,\mathrm{deg}^2$ , based on unique pairs of observations binned to 30 d samples (Carbone et al. 2016). This then gives a density of $(4.9\pm1.1) \times 10^{-4}\,\mathrm{deg}^{-2}$ (total transients) or $(1.5\pm0.6) \times 10^{-4}\,\mathrm{deg}^{-2}$ (unique sources). These are at a flux density of $>\!1.2\,\mathrm{mJy}$ (based on our median sensitivity and a 5 $\sigma$ detection threshold).

Figure 11 shows this two-epoch transient source density compared to other results (and constraining limits) from the literaturew.

5. Summary

We have introduced the VAST Pilot Survey and described the survey strategy and observations to date. We outlined the VAST transient detection pipeline, and as a demonstration of its functionality, used it to analyse two regions (3 and 4) of the Phase I Pilot Survey.

In this analysis we found 28 highly variable and transient sources out of a total of 155 071 sources detected in these regions. These sources included seven known pulsars, seven stars (three of which have no previous radio detections), two AGN, six galaxies and six sources yet to be identified.

We found a transient source density of $(1.5\pm0.6) \times 10^{-4}\,\mathrm{deg}^{-2}$ for a flux density limit of $>\!1.2\,\mathrm{mJy}$ and a $5\sigma$ detection threshold. This is consistent with rates measured in other surveys.

Future work on the VAST pipeline will incorporate better methods for reducing the false positive rate, in particular for excluding artefacts due to nearby bright sources. We are also developing a pipeline to enable faster continuum imaging of visibility datasets (following a similar method as that of Wang et al. 2021) so that we can process data as soon as possible after the observations.

In the coming year, we will conduct the Phase II Pilot Survey, incorporating both low-band and mid-band observations. The results presented here are a significant milestone in demonstrating ASKAP’s capability to detect transient sources once VAST begins full survey operations in 2022.