An insight into the extragalactic transient and variable microJy radio sky across multiple decades
Jack F. Radcliffe, Robert. J. Beswick, A. P. Thomson, Michael A. Garrett, Peter D. Barthel, Thomas W. B. Muxlow
MMNRAS , 1–18 (2018) Preprint 30 September 2019 Compiled using MNRAS L A TEX style file v3.0
An insight into the extragalactic transient and variablemicroJy radio sky across multiple decades
Jack F. Radcliffe, , , , (cid:63) Robert. J. Beswick, A. P. Thomson, Michael A. Garrett, , Peter D. Barthel, and Thomas W. B. Muxlow Kapteyn Astronomical Institute, University of Groningen, 9747 AD Groningen, The Netherlands Jodrell Bank Centre for Astrophysics, School of Physics & Astronomy, The University of Manchester, Alan Turing Building, Oxford Road, Manchester M13 9PL, UK ASTRON, the Netherlands Institute for Radio Astronomy, Postbus 2, 7990 AA Dwingeloo, The Netherlands Department of Physics, University of Pretoria, Lynnwood Road, Hatfield, Pretoria, 0083, South Africa Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
The mJy variable extragalactic radio sky is known to be broadly non-changing withapproximately of persistent radio sources exhibiting variability which is largelyAGN-related. In the faint ( < mJy) flux density regime, it is widely accepted that theradio source population begins to change from AGN dominated to star-formationdominated, together with an emergent radio-quiet AGN component. Very little isknown about the variable source component in this sub-mJy regime. In this paper, weprovide the first insight into the µ Jy variable sky by performing a careful analysis usingdeep VLA data in the well studied GOODS-N field. Using five epochs spread across 22years, we investigate approximately 480 radio sources finding 10 which show signs ofvariability. We attribute this variability to the presence of an AGN in these systems. Weconfirm and extend the results of previous surveys, finding that variability in the faintradio sky is rather modest with only ≤
2% of sources exhibiting significant variabilitybetween any two epochs. We find that 70% of variable sources show variability ontimescales of a few days whilst on longer decadal time-scales, the fraction of variablesources decreases to < . This suggests that the radio variability peaks on shortertimescales as suggested by other studies. We find that 80% of variable sources haveVLBI counterparts, and we use multi-wavelength data to infer that these may wellbe core-dominated FR-I sources as postulated by wide-field VLBI surveys and semi-empirical simulations. Key words: radio continuum: transients – galaxies: active – techniques: interfero-metric
Extremely deep radio observations of extragalactic fieldshave now revealed the nature of the persistent sub-100 µ Jy radio source population. It is generally agreed that at µ Jy flux densities the radio source population is dominated bymoderate to high redshift star-forming galaxies and radio-quiet active galactic nuclei (AGN). This is in contrast to theradio source counts at flux densities of ∼ mJy and above,which are dominated by radio galaxies and quasars (e.g.Muxlow et al. 2005; Padovani et al. 2015; Smolˇci´c et al.2017, and references therein).Whilst deep field studies are providing a view of the µ Jy radio sky in the spatial and spectral domains, at present (cid:63) E-mail: jack.radcliff[email protected] comparatively little is known about their intrinsic variabil-ity in this faint flux density regime. Previous studies haveshown that the variable and transient radio sky consists ofa range of interesting astrophysical phenomena. Within ourown galaxy, these can range from flaring stars and X-raybinaries (e.g. Thyagarajan et al. 2011; Williams et al. 2013)to novae and pulsars (e.g. Healy et al. 2017) whilst the tran-sient and variable extragalactic radio sky includes AGN,tidal disruption events, supernovae, gamma ray bursts, fastradio bursts and the recently observed binary neutron starmergers (e.g. Berger et al. 2003; Fong et al. 2014; Lorimeret al. 2007; Hallinan et al. 2017). Of the slower transients(on timescales of minutes or longer), AGN are by far themost common extragalactic radio variable sources ( (cid:38) perdeg above 0.3 mJy) and are found to vary at all frequencieson time-scales from minutes to years (e.g. Dennett-Thorpe © a r X i v : . [ a s t r o - ph . GA ] S e p J. F. Radcliffe et al. & de Bruyn 2002; Lovell et al. 2008; Ofek et al. 2011; Hodgeet al. 2013; Mooley et al. 2013; Pietka et al. 2015; Mooleyet al. 2016).Many transient surveys are focused on relatively short-duration events, such as radio afterglows of gamma-raybursts, however there have been few studies on the slowly-varying radio sky. From the few longer-term variability sur-veys that exist, the faint ( > . ) radio sky appears to bebroadly non-variable, with only ∼ . of sources showingvariations over week to 1.5 year time-scales (Mooley et al.2016). Over still longer time-scales (7-20 years), it was foundthat only ∼
1% of extragalactic sources in the mJy regime(excluding galactic transients) are variable (Croft et al. 2010;Becker et al. 2010; Bannister et al. 2011b).However, these previous studies have a number of lim-itations. Long term (several years) variability surveys arereliant on archival data such as the FIRST/NVSS surveywith limited sensitivity (1 σ , 0.2 mJy) and thus are only ableto probe the brighter end of the transient radio luminosityfunction. The next generation of radio interferometers, suchas the Karl G. Jansky Array (JVLA; Perley et al. 2011),e-MERLIN, ASKAP (Johnston et al. 2008), LOFAR (vanHaarlem et al. 2013), MeerKAT (Booth et al. 2009), Aper-tif/WSRT (Oosterloo et al. 2010), the Square Kilometer Ar-ray (SKA; Dewdney et al. 2009) and the Next GenerationVLA (e.g. Carilli et al. 2015), with their much improved sur-vey speeds and snapshot sensitivity will be able to study thefaint transient radio population across large areas of the sky.The GOODS-N field has had multiple µ Jy sensitivityobservations at 1.4 GHz spanning over 22 years using boththe Very Large Array (VLA; Richards 2000; Biggs & Ivison2006; Morrison et al. 2010) and the recently upgraded JVLA(Owen 2018, Radcliffe et al. in prep.) with cadence timesfrom a few days to decades . In this paper, our goal is touse this unique data-set to investigate the µ Jy variable skyover a 22-year period, probing the variable properties of over400 hundred faint extragalactic sources.This paper is organised as follows. In Section 2 wepresent the observations, data reduction and imaging. Weoutline the steps used to define our variable sample in Sec-tion 3 and our results are presented in Section 4. We discussthe implications of our results and compare to other vari-ability surveys in Section 5. Concluding remarks are givenin Section 6.Throughout this paper we adopt a spatially-flat 6-parameter ΛCDM cosmology with H = . ± . − Mpc , Ω m = . ± . and Ω Λ = . ± . (Planck Collabo-ration et al. 2016). For spectral index measurements, we usethe convention S ν ∝ ν α throughout, where S ν is the radiointegrated flux density and α is the intrinsic source spectralindex. In November 1996, a total of 50 hours of VLA (A-configuration; VLA project code AR0368) 1.4 GHz observa-tions, centred within the GOODS-North region at J2000 RA h m . s DEC + ◦ (cid:48) . (cid:48)(cid:48) × S . = . ) was used as aamplitude, phase and bandpass calibrator and was observedevery 40 minutes. 3C286 was used as the flux density cal-ibrator where the assumed flux densities of 15.328 Jy and14.947 Jy at 1.385 and 1.435 GHz respectively are derivedfrom the Perley & Butler (2017) absolute flux density scale.We discuss how differing flux scales between various obser-vations are a vital consideration for any variability study inSection 3.4.All calibration was conducted using the Common As-tronomy Software Applications (CASA; McMullin et al.2007). The raw VLA data sets for the entire 50 hours wereimported into a single u v measurement set using importvla .Due to known issues with the old VLA MODCOMP con-trol computers, the task fixvis was used to correctly re-calculate the u , v and w visibility coordinates. These datawere then inspected and Radio Frequency Interference (RFI)was removed using a combination of CASA tasks flagdata , viewer and plotms along with the rfigui viewer packagedwith the AOFlagger software (Offringa et al. 2012). Oncethe strongest RFI was excised, a gain curve, describing theforward gain of each antenna, was derived using gencal .The absolute flux scale was set using an L-band model of3C286 using setjy (Perley & Butler 2017). Bandpass cal-ibration was conducted on a per scan basis using calibra-tor, J1313+6735, and a complex gain calibration (ampli-tude and phase) was conducted on J1313+6735 and 3C286with five minute averaging interval using gaincal . This cal-ibration table was used to bootstrap the flux density scalefrom 3C286 to J1313+6735. The integrated flux density ofJ1313+675 was found to be be 2.396 Jy at 1.40 GHz whichis within 0.2% of the quoted value of 2.40 Jy. These calibration solutions were applied to theGOODS-N target field, which was split from the originaldata set. The data was reweighed to obtain the ideal sensi-tivity (using statwt ), and any remaining RFI was excisedfrom the data. Confusing sources were removed from thedata following the procedures outlined in Section 2.3 andthen self-calibration was performed on the target field (cal- Note that in this paper, we use observations from before andafter the VLA electronics upgrade hence the term JVLA will cor-respond to post-upgrade and the term VLA to pre-upgrade. https://science.nrao.edu/facilities/vla/observing/callist MNRAS000
1% of extragalactic sources in the mJy regime(excluding galactic transients) are variable (Croft et al. 2010;Becker et al. 2010; Bannister et al. 2011b).However, these previous studies have a number of lim-itations. Long term (several years) variability surveys arereliant on archival data such as the FIRST/NVSS surveywith limited sensitivity (1 σ , 0.2 mJy) and thus are only ableto probe the brighter end of the transient radio luminosityfunction. The next generation of radio interferometers, suchas the Karl G. Jansky Array (JVLA; Perley et al. 2011),e-MERLIN, ASKAP (Johnston et al. 2008), LOFAR (vanHaarlem et al. 2013), MeerKAT (Booth et al. 2009), Aper-tif/WSRT (Oosterloo et al. 2010), the Square Kilometer Ar-ray (SKA; Dewdney et al. 2009) and the Next GenerationVLA (e.g. Carilli et al. 2015), with their much improved sur-vey speeds and snapshot sensitivity will be able to study thefaint transient radio population across large areas of the sky.The GOODS-N field has had multiple µ Jy sensitivityobservations at 1.4 GHz spanning over 22 years using boththe Very Large Array (VLA; Richards 2000; Biggs & Ivison2006; Morrison et al. 2010) and the recently upgraded JVLA(Owen 2018, Radcliffe et al. in prep.) with cadence timesfrom a few days to decades . In this paper, our goal is touse this unique data-set to investigate the µ Jy variable skyover a 22-year period, probing the variable properties of over400 hundred faint extragalactic sources.This paper is organised as follows. In Section 2 wepresent the observations, data reduction and imaging. Weoutline the steps used to define our variable sample in Sec-tion 3 and our results are presented in Section 4. We discussthe implications of our results and compare to other vari-ability surveys in Section 5. Concluding remarks are givenin Section 6.Throughout this paper we adopt a spatially-flat 6-parameter ΛCDM cosmology with H = . ± . − Mpc , Ω m = . ± . and Ω Λ = . ± . (Planck Collabo-ration et al. 2016). For spectral index measurements, we usethe convention S ν ∝ ν α throughout, where S ν is the radiointegrated flux density and α is the intrinsic source spectralindex. In November 1996, a total of 50 hours of VLA (A-configuration; VLA project code AR0368) 1.4 GHz observa-tions, centred within the GOODS-North region at J2000 RA h m . s DEC + ◦ (cid:48) . (cid:48)(cid:48) × S . = . ) was used as aamplitude, phase and bandpass calibrator and was observedevery 40 minutes. 3C286 was used as the flux density cal-ibrator where the assumed flux densities of 15.328 Jy and14.947 Jy at 1.385 and 1.435 GHz respectively are derivedfrom the Perley & Butler (2017) absolute flux density scale.We discuss how differing flux scales between various obser-vations are a vital consideration for any variability study inSection 3.4.All calibration was conducted using the Common As-tronomy Software Applications (CASA; McMullin et al.2007). The raw VLA data sets for the entire 50 hours wereimported into a single u v measurement set using importvla .Due to known issues with the old VLA MODCOMP con-trol computers, the task fixvis was used to correctly re-calculate the u , v and w visibility coordinates. These datawere then inspected and Radio Frequency Interference (RFI)was removed using a combination of CASA tasks flagdata , viewer and plotms along with the rfigui viewer packagedwith the AOFlagger software (Offringa et al. 2012). Oncethe strongest RFI was excised, a gain curve, describing theforward gain of each antenna, was derived using gencal .The absolute flux scale was set using an L-band model of3C286 using setjy (Perley & Butler 2017). Bandpass cal-ibration was conducted on a per scan basis using calibra-tor, J1313+6735, and a complex gain calibration (ampli-tude and phase) was conducted on J1313+6735 and 3C286with five minute averaging interval using gaincal . This cal-ibration table was used to bootstrap the flux density scalefrom 3C286 to J1313+6735. The integrated flux density ofJ1313+675 was found to be be 2.396 Jy at 1.40 GHz whichis within 0.2% of the quoted value of 2.40 Jy. These calibration solutions were applied to theGOODS-N target field, which was split from the originaldata set. The data was reweighed to obtain the ideal sensi-tivity (using statwt ), and any remaining RFI was excisedfrom the data. Confusing sources were removed from thedata following the procedures outlined in Section 2.3 andthen self-calibration was performed on the target field (cal- Note that in this paper, we use observations from before andafter the VLA electronics upgrade hence the term JVLA will cor-respond to post-upgrade and the term VLA to pre-upgrade. https://science.nrao.edu/facilities/vla/observing/callist MNRAS000 , 1–18 (2018) he variable extragalactic radio sky Table 1.
A summary of the individual VLA/JVLA epochs utilised in this paper. The epochs correspond to the data that is comparedon long timescales, typically 7-22 years whereas sub-epochs are those data used to establish any short term variability typically within afew days. This epoch and sub-epoch naming conventions used here are adopted throughout the paper. Note that the VLA observationswere comprised of 24 observations spread over one month where the total observing time was approximately 47 hours on source.Telescope VLA Project Code Epoch Sub-epoch Date (UT) Frequency rms[GHz] [ µ Jy beam − ]VLA AR0368 VLA1996 1996 Nov 01 - 1996 Dec 01 1.35-1.45 5.6JVLA TLOW0001 JVLA2011 JVLA2011 ep1 2011 Sep 09 12:57:47 - 17:56:57 1.00-2.03 2.802011 Sep 09 18:00:25 - 22:59:35JVLA2011 ep2 2011 Sep 11 11:49:53 - 15:49:012011 Sep 11 15:49:11 - 20:48:202011 Sep 11-12 20:48:26 - 01:47:36JVLA 18A-392 JVLA2018 JVLA2018 ep1 2018 May 09 06:27:27 - 08:25:28 1.00-2.03 3.982018 May 12 06:21:33 - 08:19:39JVLA2018 ep2 2018 May 15 03:24:20 - 05:22:262018 May 16 01:14:41 - 03:12:432018 May 19 01:23:49 - 03:34:10 Figure 1. ◦ × ◦ image of the 1996 VLA data illustrating the source removal routine used outlined in Section 2.3. The yellow dashedcircle corresponds to the HPBW of a VLA antenna at 1.4 GHz ( ∼ (cid:48) . ) while the dotted circles correspond to the VLA HPBW at1 GHz and 2 GHz ( ∼ (cid:48) . and ∼ (cid:48) . respectively). The symmetric logarithmic colour-scale is the same for both panels and ranges from − µ Jy beam − to . − Left Panel:
The data before the source removal process. Ripples from off-axis sources severely affect thecentral target field. Inset is a cut-out of S1 (J2000 12h35m38.044s +62d19m32.097s) illustrating the un-deconvolvable sidelobes presentin the 1996 VLA data which are directed towards the phase centre.
Right Panel:
The data post the removal of confusing sources revealingsources located in the centre of the field. ibrating each polarisation separately) using multiple tar-get sources, contained within a central 10 (cid:48) diameter area,as a model. Due to the complexity of the sky model, self-calibration used an iterative strategy with solution intervalsbeginning at 20 mins and culminating with 2 min solution in-tervals for the final self-calibration cycle. No amplitude selfcalibration was conducted.
The JVLA (A-configuration) observed the GOODS-N fieldfor a total of 38 hrs in between August and September 2011(VLA project code TLOW0001; PI: F. Owen) and a further10 hrs in May 2018 (VLA project code 18A-392; PI: J. Rad-cliffe). For the purposes of this study we only use 20 hoursof the JVLA2011 data which provides two epochs with sim-ilar observing times (see Table 1). The observations have a
MNRAS , 1–18 (2018)
J. F. Radcliffe et al. total bandwidth of 1024 MHz comprising of 16 spw with 64channels per spw.The source J1313+6735 was used as a phase and band-pass calibrator and 3C286 was used as the flux calibrator.Initial flagging was conducted using
AOFlagger (Offringaet al. 2012). The VLA CASA pipeline was then used toconduct the flux density scaling and phase referencing. Someepochs required additional flagging (and re-calibration) afterthe initial pipeline calibration, especially if strong RFI hadnot been flagged adequately before or during calibration. Intotal, approximately 20-40% of the data in each epoch waslost due to RFI. As has been noted in previous publications dealing with theGOODS-N field, there are significant issues with bright, con-fusing sources which limit the dynamic range of the centralimage. In order to deal with these sources, we were requiredto remove a number of sources from these data. Confusingsources were identified by generating a large, ◦ × ◦ , un-deconvolved image using wsclean (Offringa et al. 2014).The influence of off-axis confusing sources was more se-vere in the old VLA observations. The pre-WIDAR, VLAcorrelator had significant multiplicative off-axis errors whenthe continuum mode was used. These errors are correlatedwith the visibility strength resulting in un-deconvolvableside-lobes around the brightest sources which are orientatedtowards the phase centre (see the inset in Fig. 1). This hasbeen previously noted by Richards (2000), Biggs & Ivison(2006) and Morrison et al. (2010). As a result, sources asfar as 2 ◦ from the phase centre had to be carefully sub-tracted. For the JVLA data, these errors have been fixedwith the upgrade to the WIDAR correlator. As a result,fewer sources needed to be peeled and the main source ofconfusion was due to the primary beam attenuation induc-ing direction-dependent errors upon bright sources near theprimary beam half-power.With confusing sources identified, each was removed in-dividually using the following procedure:(i) A frequency-dependent model of the source was gen-erated using CASA task tclean and the multi-term multifrequency synthesis (MT-MFS) algorithm (Rau & Cornwell2011). Separate models were generated for each of the cir-cular polarisations due of the significant beam squint of theVLA antennas.(ii) The model generated was then Fourier transformedand gridded in order to generate model visibilities (usingtask ft ).(iii) Calibration corrections, derived by comparing the ob-served visibilities to the model visibilities and solving perantenna, were generated using gaincal . The observed visi-bilities were adjusted by these corrections using task apply-cal .(iv) Steps (i)-(iii) were repeated until the side-lobe struc-ture around the confusing source was reduced. https://science.nrao.edu/facilities/vla/data-processing/pipeline (v) The final self-calibration model visibilities were sub-tracted from the observed visibilities using task uvsub . Thisremoved the source and the side-lobe structure from the ob-served visibilities.(vi) The observed visibilities now have the correct phaseand amplitude solutions at the location of the confusingsource rather than the pointing centre. Therefore, we mustrevert our complex gains to the centre of the target field.To do this, the amplitude and phase corrections containedin final calibration table were inverted using the contributedCASA task invgain (Hales 2016).(vii) These inverted corrections were applied (using task applycal ) to adjust the observed visibilities to be correctfor the point centre again.(viii) Finally, steps (i)-(vii) were repeated for the otherconfusing sources.The results of this process are presented in Fig. 1. Thisremoval method is identical to the ‘peeling’ method com-monly used in radio interferometric data processing (e.g.Owen & Morrison 2008; Intema et al. 2009) with the excep-tion that we do not reinsert the self-calibrated source modelinto the visibilities after step (viii). This meant that we werenot required to deconvolve these sources when generating thefinal images, thus reducing the image size required. A large (cid:48) × (cid:48) image was generated for each data set usingCASA task tclean . Imaging for both the JVLA and VLAused the wprojectft algorithm (Cornwell et al. 2008). Thisalgorithm accounts for the non-coplanar array term (knownas the w term) thus minimising the induced smearing awayfrom the phase centre. For these VLA data, deconvolutionwas performed using the clarkstokes algorithm. A primarybeam model for the VLA correction was derived using theCASA recipe makePB and then was divided through the im-age using impbcor .For these JVLA data, the imaging method was slightlydifferent. Due to the large fractional bandwidth ( ∼ ),deconvolution was performed using the multi-term multi-frequency synthesis (MT-MFS) technique which permitsmore accurate deconvolution by taking into account the fre-quency dependence of the sky model (Rau & Cornwell 2011).In addition, this method allows the correction of the differ-ence in reference frequencies utilising the in-band spectralindices to correct fluxes to that of the same frequency asthe VLA data. This technique is tested and outlined in Sec-tion 3.3. The images were corrected for the primary beamattenuation and primary beam induced spectral curvatureusing the widebandpbcor task. The final data sets, epochsand sensitivities of the various epochs are described in Ta-ble 1. MNRAS000
AOFlagger (Offringaet al. 2012). The VLA CASA pipeline was then used toconduct the flux density scaling and phase referencing. Someepochs required additional flagging (and re-calibration) afterthe initial pipeline calibration, especially if strong RFI hadnot been flagged adequately before or during calibration. Intotal, approximately 20-40% of the data in each epoch waslost due to RFI. As has been noted in previous publications dealing with theGOODS-N field, there are significant issues with bright, con-fusing sources which limit the dynamic range of the centralimage. In order to deal with these sources, we were requiredto remove a number of sources from these data. Confusingsources were identified by generating a large, ◦ × ◦ , un-deconvolved image using wsclean (Offringa et al. 2014).The influence of off-axis confusing sources was more se-vere in the old VLA observations. The pre-WIDAR, VLAcorrelator had significant multiplicative off-axis errors whenthe continuum mode was used. These errors are correlatedwith the visibility strength resulting in un-deconvolvableside-lobes around the brightest sources which are orientatedtowards the phase centre (see the inset in Fig. 1). This hasbeen previously noted by Richards (2000), Biggs & Ivison(2006) and Morrison et al. (2010). As a result, sources asfar as 2 ◦ from the phase centre had to be carefully sub-tracted. For the JVLA data, these errors have been fixedwith the upgrade to the WIDAR correlator. As a result,fewer sources needed to be peeled and the main source ofconfusion was due to the primary beam attenuation induc-ing direction-dependent errors upon bright sources near theprimary beam half-power.With confusing sources identified, each was removed in-dividually using the following procedure:(i) A frequency-dependent model of the source was gen-erated using CASA task tclean and the multi-term multifrequency synthesis (MT-MFS) algorithm (Rau & Cornwell2011). Separate models were generated for each of the cir-cular polarisations due of the significant beam squint of theVLA antennas.(ii) The model generated was then Fourier transformedand gridded in order to generate model visibilities (usingtask ft ).(iii) Calibration corrections, derived by comparing the ob-served visibilities to the model visibilities and solving perantenna, were generated using gaincal . The observed visi-bilities were adjusted by these corrections using task apply-cal .(iv) Steps (i)-(iii) were repeated until the side-lobe struc-ture around the confusing source was reduced. https://science.nrao.edu/facilities/vla/data-processing/pipeline (v) The final self-calibration model visibilities were sub-tracted from the observed visibilities using task uvsub . Thisremoved the source and the side-lobe structure from the ob-served visibilities.(vi) The observed visibilities now have the correct phaseand amplitude solutions at the location of the confusingsource rather than the pointing centre. Therefore, we mustrevert our complex gains to the centre of the target field.To do this, the amplitude and phase corrections containedin final calibration table were inverted using the contributedCASA task invgain (Hales 2016).(vii) These inverted corrections were applied (using task applycal ) to adjust the observed visibilities to be correctfor the point centre again.(viii) Finally, steps (i)-(vii) were repeated for the otherconfusing sources.The results of this process are presented in Fig. 1. Thisremoval method is identical to the ‘peeling’ method com-monly used in radio interferometric data processing (e.g.Owen & Morrison 2008; Intema et al. 2009) with the excep-tion that we do not reinsert the self-calibrated source modelinto the visibilities after step (viii). This meant that we werenot required to deconvolve these sources when generating thefinal images, thus reducing the image size required. A large (cid:48) × (cid:48) image was generated for each data set usingCASA task tclean . Imaging for both the JVLA and VLAused the wprojectft algorithm (Cornwell et al. 2008). Thisalgorithm accounts for the non-coplanar array term (knownas the w term) thus minimising the induced smearing awayfrom the phase centre. For these VLA data, deconvolutionwas performed using the clarkstokes algorithm. A primarybeam model for the VLA correction was derived using theCASA recipe makePB and then was divided through the im-age using impbcor .For these JVLA data, the imaging method was slightlydifferent. Due to the large fractional bandwidth ( ∼ ),deconvolution was performed using the multi-term multi-frequency synthesis (MT-MFS) technique which permitsmore accurate deconvolution by taking into account the fre-quency dependence of the sky model (Rau & Cornwell 2011).In addition, this method allows the correction of the differ-ence in reference frequencies utilising the in-band spectralindices to correct fluxes to that of the same frequency asthe VLA data. This technique is tested and outlined in Sec-tion 3.3. The images were corrected for the primary beamattenuation and primary beam induced spectral curvatureusing the widebandpbcor task. The final data sets, epochsand sensitivities of the various epochs are described in Ta-ble 1. MNRAS000 , 1–18 (2018) he variable extragalactic radio sky . . . . . . P u r i t y ( P r ) S/N detection threshold . . . ∇ P r . . . h S PYBDSF / S m o d e l i S p S int VLAJVLA6080100 % r ec o v e r e d s o u r ce s S/N model . . . . . σ S PYBDSF / S m o d e l JVLA2011JVLA2018VLA1996
Figure 2.
Summary of the various source extraction tests performed using
PYBDSF . Left Panel:
The purity, P r parameter against S/Ndetection threshold. The number of false detections in all epochs ( > ) is large at S / N < . The parameter asymptotes for all epochsat around S/N thresholds of 5.5 as illustrated by the gradient of the purity parameter, ∇ P r , in the bottom panel. Right Panel:
Theresults from the flux recovery simulations using
PYBDSF . The top panel illustrates the percentage of model sources recovered by thesource extractions software. The middle panel shows the average flux recovery performance while the bottom panel shows the typical 1 σ variance on the average recovered flux densities. The 7 σ detection threshold adopted in these analyses provides good peak brightnessrecovery with a typical error of around and recovers the majority of sources inserted ∼ . With the images generated, we can now proceed onto thevariability analysis. This section aims to summarise our se-lection method for defining the variable population whileoutlining and addressing the various systematic offsets whichcould induce artificial variability between the various epochs.These effects include the source extraction software, the ref-erence frequency, the differing u v coverage and the absoluteflux scale. We shall address in turn each of these issues. Be-cause we expect from previous studies that the number ofvariable sources will be low, systematics which can affectindividual source flux densities (e.g. source extraction) aremuch more important to address than those which affect thesource flux densities as a whole (e.g. absolute flux densityscaling) as these can be modelled and corrected. One of the most important systematics to address is thevarying performance of the chosen source extraction routine.There are two main effects to consider, the first being thefalse detection rate and the second being the recovered fluxdensities of the extraction software. To identify sources wewill use
PYBDSF (Mohan & Rafferty 2015). The reference frequency is the centre of the band and, whenusing MT-MFS, is where the wide-band image flux density shouldequal the source flux density (see Section 3.3)
In order to determine the appropriate signal-to-noise (S/N)threshold at which to allow
PYBDSF to identify and cata-logue sources in our maps, we adopt the empirical ‘purity’parameter of Stach et al. (2018). We ran
PYBDSF on thereal and inverted images (i.e. the real maps multiplied by − ) for each epoch with S/N thresholds between - , andwe then evaluated the ‘purity parameter’, P r = N p − N n N p , (1)for each image, where N p is the number of sources detectedabove a given S/N limit in the real image and N n is thenumber of sources detected above the same S/N limit in theinverted image. For an idealised image comprising of Gaus-sian noise and a real source population, at an arbitrarilyhigh S/N threshold the purity parameter should asymptotetowards unity (i.e. N p > and N n = ). However, at lowerS/N thresholds, we expect to detect a combination of gen-uine sources and spurious noise peaks in the real image, andonly spurious noise peaks in the inverted image.As Figure 2 illustrates, the purity parameter asymptotestowards unity around a S/N threshold of 5 for epochs ob-served with JVLA, whereas for the pre-upgrade VLA mapsthe purity parameter asymptotes to ∼ . . We determinethat the VLA false-positive rate asymptotes to ∼ for S / N (cid:38) . . This is indicative of non-Gaussian noise proper-ties in the VLA maps due to a combination of multiplicativebaseline errors (which cause visible corrugations in the maps MNRAS , 1–18 (2018)
J. F. Radcliffe et al. towards the pointing centre), and the more limited u v cov-erage of the old array with respect to the post-2010 JVLA.Visual inspection of the inverted maps confirms this andthese highly localised areas are excluded during subsequentcataloguing. In addition to ensuring that the S/N threshold reduces thenumber of false positives to a minimum, for variability stud-ies, we also have to evaluate the source-extraction techniqueand any systematic biases to the recovered flux densities.One of these biases is the well-noted ‘flux-boosting’ effect(Jauncey 1968; Vernstrom et al. 2016). This causes the fittedflux densities, especially in the low S/N regime, to be system-atically over-estimated due to a combination of confusionand instrumental noise. For deep radio surveys, this mani-fests itself as an over-estimation of the faint source counts.In the case of variability studies, this effect can be significantif the sources between the two epochs being compared havediffering S/N ratios. The flux boosting effect is a function ofS/N, therefore a source will be flux boosted in the less sensi-tive image resulting in artificial variability between epochs.The severity of this effect is also related to the choice of thesource extraction routine.In these tests, we performed a hybrid approach whichuses real observations (including residual calibration errors)with simulated sources injected. This was done separatelyfor the VLA and JVLA due to the differing u v coverageand therefore different noise characteristics. Before simu-lated sources were injected, a sky model was generated using wsclean , utilising the auto-masking routine to ensure thatonly real sources are removed and the noise characteristicsare not changed. This model is Fourier transformed and sub-tracted from the visibilities and then the u v -data is imagedagain to generate blank × pixel noise fields. These im-ages are used to generate an r.m.s. map using PYBDSF .These model sources comprise of simple delta functionswith a range of S/N ratios. These are convolved with therestoring beam and added into the noise fields. The posi-tions of these sources are randomised across the field of viewwith the restriction that a source cannot be 10 pixels fromanother source in order to prevent blending complications.In addition, the periphery of the image ( ∼ of the totalimage size) is avoided in order to reduce complications dueto CLEAN aliasing. Each image is then catalogued using PYBDSF using a 5.5 σ detection threshold (as motivated bythe purity tests).The results of these simulations are summarised in Fig-ure 2. The peak brightnesses of the simulated sources aregenerally recovered by PYBDSF at around a S/N of 7. Be-low this threshold, the detected source flux distribution iscensored where some sources which have their flux densitiesreduced by natural noise fluctuations are not detected. Asa result, we see the flux-boosting effect which increases theaverage recovered flux densities in the low S/N regime. It isworth noting here that at intermediate S/N ratios ( S / N = - ) the peak brightnesses are, on average, underestimatedby ∼ σ for our detection thresholdwhich provides a high reliability (as shown by the purity andsources recovered metrics) and adequate peak flux recovery(albeit with ∼ errors) in this S/N regime. u v coverage One potential source of induced variability could be causedby the differing u v -coverage of the various epochs. The multi-frequency synthesis imaging approach improves the sensi-tivity and fidelity of interferometric images by utilising thelarge bandwidths to fill the u v plane. The large differencesbetween the bandwidths of the VLA and JVLA observationstherefore causes a mismatch in the u v coverage. This effect ismore significant on complex structures with differing radiopower on multiple spatial scales, however, the majority ofthe radio sources in this field are not extended and subtendonly 1-2 VLA beam widths. To address this issue, we ex-cluded large extended sources (with classical radio jet struc-tures) by masking them before cataloguing, and we did notconsider those sources that required models comprising ofmultiple Gaussian components by PYBDSF in order to modeltheir flux distribution. In addition, we only consider the peakbrightnesses of sources as these radio emission mechanismswill have the same power on all spatial scales and are thusless sensitive to differences in the u v coverage. It is worthnoting that these cuts will inevitably bias our result as wemay miss variable sources associated with extended objects.However such variability would be difficult to disentanglefrom artificial variability induced by the differing u v cover-age between the epochs. The reference frequency is the frequency in which the wide-band flux densities in interferometric images equal the nar-row band flux densities. This frequency ν ref is simply definedas the centre of the observing band, ν ref = ν high − ν low (2)where ν low and ν high corresponds to the low and high edgesof the observing bandwidth . The intensity of an image issuch that the wideband flux density is equal to the narrow-band source flux density at the reference frequency. Withouttaking into account the frequency dependence of the skybrightness distribution, this intensity is simply the weightedmean of the flux density across the bandwidth. This makesthe returned fluxes susceptible to bias from the flagging oflarge frequency ranges. However, if MTMFS is used, the fluxdensities are fitted across the band and so the returned wide-band flux density will be closer to the real narrow-band fluxdensity (Emonts private communication). This makes thereturned flux densities more robust to any flagged data. see https://casa.nrao.edu/casadocs/casa-5.6.0/imaging/synthesis-imaging/wide-band-imaging for more informationMNRAS000
J. F. Radcliffe et al. towards the pointing centre), and the more limited u v cov-erage of the old array with respect to the post-2010 JVLA.Visual inspection of the inverted maps confirms this andthese highly localised areas are excluded during subsequentcataloguing. In addition to ensuring that the S/N threshold reduces thenumber of false positives to a minimum, for variability stud-ies, we also have to evaluate the source-extraction techniqueand any systematic biases to the recovered flux densities.One of these biases is the well-noted ‘flux-boosting’ effect(Jauncey 1968; Vernstrom et al. 2016). This causes the fittedflux densities, especially in the low S/N regime, to be system-atically over-estimated due to a combination of confusionand instrumental noise. For deep radio surveys, this mani-fests itself as an over-estimation of the faint source counts.In the case of variability studies, this effect can be significantif the sources between the two epochs being compared havediffering S/N ratios. The flux boosting effect is a function ofS/N, therefore a source will be flux boosted in the less sensi-tive image resulting in artificial variability between epochs.The severity of this effect is also related to the choice of thesource extraction routine.In these tests, we performed a hybrid approach whichuses real observations (including residual calibration errors)with simulated sources injected. This was done separatelyfor the VLA and JVLA due to the differing u v coverageand therefore different noise characteristics. Before simu-lated sources were injected, a sky model was generated using wsclean , utilising the auto-masking routine to ensure thatonly real sources are removed and the noise characteristicsare not changed. This model is Fourier transformed and sub-tracted from the visibilities and then the u v -data is imagedagain to generate blank × pixel noise fields. These im-ages are used to generate an r.m.s. map using PYBDSF .These model sources comprise of simple delta functionswith a range of S/N ratios. These are convolved with therestoring beam and added into the noise fields. The posi-tions of these sources are randomised across the field of viewwith the restriction that a source cannot be 10 pixels fromanother source in order to prevent blending complications.In addition, the periphery of the image ( ∼ of the totalimage size) is avoided in order to reduce complications dueto CLEAN aliasing. Each image is then catalogued using PYBDSF using a 5.5 σ detection threshold (as motivated bythe purity tests).The results of these simulations are summarised in Fig-ure 2. The peak brightnesses of the simulated sources aregenerally recovered by PYBDSF at around a S/N of 7. Be-low this threshold, the detected source flux distribution iscensored where some sources which have their flux densitiesreduced by natural noise fluctuations are not detected. Asa result, we see the flux-boosting effect which increases theaverage recovered flux densities in the low S/N regime. It isworth noting here that at intermediate S/N ratios ( S / N = - ) the peak brightnesses are, on average, underestimatedby ∼ σ for our detection thresholdwhich provides a high reliability (as shown by the purity andsources recovered metrics) and adequate peak flux recovery(albeit with ∼ errors) in this S/N regime. u v coverage One potential source of induced variability could be causedby the differing u v -coverage of the various epochs. The multi-frequency synthesis imaging approach improves the sensi-tivity and fidelity of interferometric images by utilising thelarge bandwidths to fill the u v plane. The large differencesbetween the bandwidths of the VLA and JVLA observationstherefore causes a mismatch in the u v coverage. This effect ismore significant on complex structures with differing radiopower on multiple spatial scales, however, the majority ofthe radio sources in this field are not extended and subtendonly 1-2 VLA beam widths. To address this issue, we ex-cluded large extended sources (with classical radio jet struc-tures) by masking them before cataloguing, and we did notconsider those sources that required models comprising ofmultiple Gaussian components by PYBDSF in order to modeltheir flux distribution. In addition, we only consider the peakbrightnesses of sources as these radio emission mechanismswill have the same power on all spatial scales and are thusless sensitive to differences in the u v coverage. It is worthnoting that these cuts will inevitably bias our result as wemay miss variable sources associated with extended objects.However such variability would be difficult to disentanglefrom artificial variability induced by the differing u v cover-age between the epochs. The reference frequency is the frequency in which the wide-band flux densities in interferometric images equal the nar-row band flux densities. This frequency ν ref is simply definedas the centre of the observing band, ν ref = ν high − ν low (2)where ν low and ν high corresponds to the low and high edgesof the observing bandwidth . The intensity of an image issuch that the wideband flux density is equal to the narrow-band source flux density at the reference frequency. Withouttaking into account the frequency dependence of the skybrightness distribution, this intensity is simply the weightedmean of the flux density across the bandwidth. This makesthe returned fluxes susceptible to bias from the flagging oflarge frequency ranges. However, if MTMFS is used, the fluxdensities are fitted across the band and so the returned wide-band flux density will be closer to the real narrow-band fluxdensity (Emonts private communication). This makes thereturned flux densities more robust to any flagged data. see https://casa.nrao.edu/casadocs/casa-5.6.0/imaging/synthesis-imaging/wide-band-imaging for more informationMNRAS000 , 1–18 (2018) he variable extragalactic radio sky .
80 0 .
85 0 .
90 0 .
95 1 .
00 1 .
05 1 .
10 1 . S . /S . S peak ˜ S int S peak S int Figure 3.
Comparison between the 1.4 GHz and 1.52 GHz fluxesof the JVLA data of epoch 5 (observed in 2018). The indigofilled and black edged histograms correspond to the ratio of peakbrightnesses and integrated flux densities respectively. The corre-sponding dashed lines indicated the median of each distribution.The large scatter of both distributions are dominated by the in-trinsic scatter in the individual source spectral indices.
One possible source of induced variability is the differ-ing reference frequencies between the VLA (1.4 GHz) andJVLA (1.52 GHz) epochs. Over the entire radio source sam-ple, the differences between the two reference frequencies willmanifest as a systematic flux scaling issue around the me-dian spectral index of the source population. Indeed, Owen(2018) finds a ∼ systematic flux density decrease be-tween the JVLA data used in these analysis and the Morri-son et al. (2010) VLA observations which consistent with thedifference in flux density expected between 1.4 and 1.52 GHzwith a mean spectral index in the expected range of − . to − . . While the net effect of this can be fitted and removed,the large intrinsic scatter in the spectral index distribution( α = − . ± . , Smolˇci´c et al. 2017) means that a smallfraction of sources with extreme deviations could signifi-cantly change the flux densities between 1.4 and 1.52 GHz.For studies of the entire population as a whole, this isnot significant, but in the case of variability where we expecta small number of variable sources, this effect could be signif-icant on the final variable population identified. To correctfor this, we utilise the MT-MFS algorithm when imaging theJVLA epochs to adjust the reference frequency to 1.4 GHz.This algorithm uses the in-band source spectral index to re-scale each source to the appropriate flux density at 1.4 GHz.To test that this was performing as expected, we used asingle JVLA data set ( JVLA2018, observed in 2018). Twoimages covering an (cid:48) × (cid:48) were generated using CASA task tclean implementing the MT-MFS (with nterms=2 ) mode.For one of the images, the reference frequency for where theevaluation of the Taylor expansion, which takes into accountthe frequency dependence of the sky model during deconvo-lution, was set to 1.4 GHz. Both of these images were pri-mary beam corrected using the widebandpbcor task. These were then identically catalogued using PYBDSF with an arbi-trary σ detection threshold (to reduce scatter induced bythe fitting routine). As Fig. 3 shows, the median differencebetween the peak brightnesses and integrated flux densitiesare 5.2% and 4.7% respectively. This corresponds to an me-dian spectral index of − . which is consistent with the lit-erature (e.g. Condon 1984; Kimball & Ivezi´c 2008; Smolˇci´cet al. 2017), thus validating this correction method in ac-counting for the reference frequency differences.The scatter in the distribution indicates that there is anon-negligible number of sources where this difference cancause 10-20% amplitude adjustments. While this is not suffi-cient for outright classification using the metric used in thispaper , the combination with other factors, such as resid-ual calibration errors, could induce artificial variability in asmall number of sources.Another possible artificial source of variability due tothe reference frequency shifting correction could be due tothe induced source spectral index causes by differing primarybeam attenuation across the bandwidth. The spectral index,used by MTMFS to re-scale the flux densities to 1.4 GHz,would be a linear combination of the source-intrinsic spec-tral index and a primary beam induced spectral index, thuscausing an incorrect rescaling of the flux densities. The extraprimary beam spectral index term, α E , can be approximatedby : α E = − ( ) (cid:18) θθ (cid:19) (cid:18) νν ref (cid:19) (3)where θ is the primary beam FWHM at the reference fre-quency, θ is the distance from the pointing centre and ν isthe observing frequency. Therefore, a flat spectrum source( α = ) located 10 (cid:48) from the pointing centre would have aadditional primary beam induced spectral index of ∼ − . .This would correspond to a shift in flux density of ∼ be-tween the frequencies of 1.52 GHz and 1.4 GHz. To solve this,the wide-band primary beam correction, applied with CASAtask widebandpbcor , models and removes the primary beaminduced spectral index, α E . Due to the use of the VLA pipeline, a different flux scalewas used for the VLA and JVLA data during calibration.However, there is little difference between these flux den-sity scales and any discrepancies are within the 3-5% error(Perley & Butler 2017). Indeed, the phase calibrator fluxdensities of all epochs agree fairly well and are all within5% of each other at 1.4 GHz. These small offsets should beproportional to the source flux density and can therefore canbe modelled and correction factors derived.To do this, we use robust linear regression (implementedin the
ODR scipy
Python package) to fit the peak fluxes be-tween each pair of epochs in order to derive single epoch cor-rection factors. Only sources with peak flux densities larger A 20% flux difference would only give a V s (Eqn. 5) of 1.414assuming a typical 10% error due to calibration and fitting see https://casa.nrao.edu/casadocs/casa-5.6.0/imaging/synthesis-imaging/wide-band-imaging MNRAS , 1–18 (2018)
J. F. Radcliffe et al. than 55 µ Jy beam − were considered because this correspondsto σ in the least sensitive epoch (VLA 1996). The moresensitive JVLA 2011 epoch 2 is used as the reference epoch.The scaling factors derived and applied to the peak bright-nesses were 1.01, 1.04, 0.96 for the JVLA 2011 epoch 1, VLA1996 and JVLA 2018 epochs respectively. In order to identify variable and transient sources we adoptthe formulation developed by Mooley et al. (2016). We as-sume that our radio interferometric noise is Gaussian dis-tributed (Condon et al. 1998) and therefore we can comparethe flux densities of a source, S , across two epochs (1,2) withthe following equation: ∆ S σ = S − S (cid:113) σ + σ . (4)where σ is the measurement error. This comprises of the PYBDSF fitting error , the in-band spectral index error ( < ), and a surface brightness error to quantify slight differ-ences in the u v coverage, pointing errors, and primary beamcorrection. We assume a flux density error for the JVLA andVLA data of between 3-5% (Morrison et al. 2010; Perley &Butler 2013, 2017). The quantity, ∆ S / σ , is expected to followthe two-sided Student-t distribution and, utilising the samenotation as Mooley et al. (2016), we define the variabilityt-statistic, V s as: V s = (cid:12)(cid:12)(cid:12)(cid:12) ∆ S σ (cid:12)(cid:12)(cid:12)(cid:12) (5)The extent of a source’s variability between two epochs canbe quantified using the ‘modulation index’, m , which is de-fined as, m = S − S S + S = ∆ S (cid:104) S (cid:105) , (6)where (cid:104) S (cid:105) is the mean flux density of the source. Due tothe expected low number of variable sources, we are moreconcerned with reliability over completeness. Therefore weadopt the rather conservative criterion for a variable sourcewhich is identical to that used by Mooley et al. (2016). Wedefine a source as variable if the variable t-statistic, V s , isbeyond the 95% confidence interval ( V s ≥ . ) and its peakbrightness changes by ≥ (corresponding to | m | > . ). This error is based upon the Condon (1997) formulation anddoes not require the Monte-Carlo approach used in
PYBDSF toevaluate errors as these are only needed for extended sources thatare not considered here For two degrees of freedom, the 95% confidence interval in theStudent-t distribution corresponds to a Gaussian probability ofmore than ± σ . For the Gaussian distribution, σ correspondsto a probability of about 1/16,000 (0.00625%). The number ofindividual measurements in our variability analysis is 3464 whichcorresponds to 0.2165 false-positive variable sources. See Mooleyet al. (2016) and p. 65-67 and Table C.8 of Bevington & Robinson(2003) for further details. To summarise, the definitions and systematics outlinedabove leads to the following strategy being adopted in thispaper, in order to identify variable and transient sources,(i) Each epoch was searched for extended sources andthese were identically masked in all epochs.(ii) For each epoch, sources were extracted using
PYBDSF using a S/N threshold of 7 σ based upon the false detectionand flux recovery tests outlined in Section 3.1. The totalnumber of radio sources considered is around 250 to 480 dueto the differing sensitivities of the epochs.(iii) Flux scaling issues due to small errors in the absoluteflux scaling and differences in restoring beams are correctedusing the routine outlined in Section 3.4.(iv) These catalogues were cross-matched and combinedusing a 1 (cid:48)(cid:48) search radius. Sources not present in some epochshave upper limits derived which correspond to 7 × the localr.m.s.(v) Using the definitions outlined in Section 3.5, V s and m were derived for all two epoch combinations and variablecandidates are identified if V s ≥ . and | m | ≥ . .(vi) Sources only present in one of the two epochs beingcompared have upper limits derived which correspond tothe 7 σ detection threshold at the source position where σ is derived using the r.m.s. maps generated by PYBDSF . If thedetection threshold is smaller than the peak brightness then V s and m are calculated, using the detection threshold as thepeak brightness, otherwise this source is not considered inthe two epoch comparison.(vii) Next we separate those sources into transient andvariable sources. A source is classed as being variable if thesource satisfies the variability criterion and is detected inmultiple epochs while a transient source must satisfy thevariability criterion but is only detected in a single epoch.(viii) Finally, variable and transient candidates arechecked manually to see if other factors such as nearby cal-ibration errors or poor source extraction fitting could haveartificially induced the variability. The aim of this paper is to investigate the long-term µ Jy radio sky. In order to achieve this, we split our JVLAepochs into two sub-epochs separated by short day-weektimescales. This is done to establish whether any long termvariability behaviour is artificially induced by any shortterm variability. In light of this, we establish two definitionscorresponding to those sources with long and short termvariability respectively. The long-term variable candidatesare those which satisfy the variability criteria in any two-epoch comparison between the VLA1996, JVLA2011, andJVLA2018 epochs, while the short-term variables are identi-fied as those which satisfy the variability criteria for any twosub-epoch comparisons (i.e. JVLA2011 ep1 and JVLA2011ep2, JVLA2018 ep1 and JVLA2018 ep2). In Figure 4 weshow the variability index, V s against the modulation index, m , for all two-epoch and sub-epoch comparisons consideredhere. Sources which are located in the gold shaded areas MNRAS000
PYBDSF using a S/N threshold of 7 σ based upon the false detectionand flux recovery tests outlined in Section 3.1. The totalnumber of radio sources considered is around 250 to 480 dueto the differing sensitivities of the epochs.(iii) Flux scaling issues due to small errors in the absoluteflux scaling and differences in restoring beams are correctedusing the routine outlined in Section 3.4.(iv) These catalogues were cross-matched and combinedusing a 1 (cid:48)(cid:48) search radius. Sources not present in some epochshave upper limits derived which correspond to 7 × the localr.m.s.(v) Using the definitions outlined in Section 3.5, V s and m were derived for all two epoch combinations and variablecandidates are identified if V s ≥ . and | m | ≥ . .(vi) Sources only present in one of the two epochs beingcompared have upper limits derived which correspond tothe 7 σ detection threshold at the source position where σ is derived using the r.m.s. maps generated by PYBDSF . If thedetection threshold is smaller than the peak brightness then V s and m are calculated, using the detection threshold as thepeak brightness, otherwise this source is not considered inthe two epoch comparison.(vii) Next we separate those sources into transient andvariable sources. A source is classed as being variable if thesource satisfies the variability criterion and is detected inmultiple epochs while a transient source must satisfy thevariability criterion but is only detected in a single epoch.(viii) Finally, variable and transient candidates arechecked manually to see if other factors such as nearby cal-ibration errors or poor source extraction fitting could haveartificially induced the variability. The aim of this paper is to investigate the long-term µ Jy radio sky. In order to achieve this, we split our JVLAepochs into two sub-epochs separated by short day-weektimescales. This is done to establish whether any long termvariability behaviour is artificially induced by any shortterm variability. In light of this, we establish two definitionscorresponding to those sources with long and short termvariability respectively. The long-term variable candidatesare those which satisfy the variability criteria in any two-epoch comparison between the VLA1996, JVLA2011, andJVLA2018 epochs, while the short-term variables are identi-fied as those which satisfy the variability criteria for any twosub-epoch comparisons (i.e. JVLA2011 ep1 and JVLA2011ep2, JVLA2018 ep1 and JVLA2018 ep2). In Figure 4 weshow the variability index, V s against the modulation index, m , for all two-epoch and sub-epoch comparisons consideredhere. Sources which are located in the gold shaded areas MNRAS000 , 1–18 (2018) he variable extragalactic radio sky − . − . − . − . . . . . . VLA1996 - JVLA2018 (21.59 yr)% var = 3 /
295 (1 . var = 5 /
290 (1 . . . . . . . . JVLA2018 - JVLA2011 (6.68 yr)% var = 5 /
471 (1 . . . . . . . . − . − . − . − . . . . . . JVLA2011 ep1 - JVLA2011 ep2 (2 days)% var = 5 /
479 (1 . . . . . . . . JVLA2018 ep1 - JVLA2018 ep2 (5 days)% var = 2 /
197 (1 . V s (Peak) m ( P e a k ) V s > . Figure 4.
The two epoch comparison comparing the variability statistic V s to the modulation index, m . The black markers correspondto sources which are present in both epochs whilst the blue triangles are those sources which have upper limits in one epoch. The markersizes are proportional to the mean peak brightness, (cid:104) S p (cid:105) . The gold region corresponds to the variability criteria where the variabilitymust be statistically significant i.e. V s > . and it must have adjusted flux by 30% (i.e. | m | > . ). The total number of variable sourcesis higher in this figure as some sources are classed as variables in more than one two epoch comparison. satisfy the variability criteria ( V s > . and | m | > . ). Thetop row corresponds to the long-term variables while thebottom row corresponds to the short-term variables. In to-tal, there are around 479 unique sources being consideredin these analyses. The number of sources varies between ap-proximately 200 and 480 due to the differing sensitivities ofthe epochs.Initial investigations revealed 15 unique variable sourcesof which 5 were excluded after manual inspection. Thesewere excluded due to a combination of poor source ex-traction, extended structure, residual calibration errors andsmearing effects. In total, therefore we found 10 uniquesources based upon the epoch and sub-epoch comparisons.The epoch-averaged peak flux densities ranging from 86 to419 µ Jy beam − with an average flux of 173 µ Jy beam − . Thepaucity of low flux density sources < µ Jy beam − may beindicative of the change in the radio properties of galax-ies being primarily driven by AGN to those driven by star-formation processes, however we note that the variabilitystatistic V s is biased towards large S/N, because the fittingerror increases with decreasing S/N.For the candidate long-term variables, around 1-1.8%of the persistent radio sources exhibit significant variabil-ity. Note that for the rest of this paper, the term ‘vari- able fraction’ indicates the percentage of sources which ex-hibits variability within the total persistent radio popula-tion. There is no clear correlation between the cadence timeand the variable fraction, however it is worth noting that theJVLA2018 and JVLA2011 two-epoch comparison is moresensitive ( ∼ . µ Jy beam − detection threshold) hence thenumber of persistent sources compared differs. If we restrictthe JVLA2011 and JVLA2018 epochs to the VLA1996 de-tection threshold, then we achieve a similar variable fractionof 1.4%. In total, we find eight sources which are classed aslong-term variables. Five of the eight sources are classed asvariable on multiple two-epoch comparisons, while two showvariability between all epochs.We classify ∼ of our sample as candidate short-termvariables, based upon measurements of their flux densitiesbetween the JVLA2011 sub-epochs and/or the JVLA2018sub-epochs. It is worth noting that the JVLA2011 sub-epochs are more sensitive than the JVLA2018 sub-epochs.If we match the JVLA2011 sub-epochs detection thresholdto the JVLA2018 sub-epoch detection thresholds, we find ahigher variable fraction of ∼ . for the JVLA2011 sub-epochs. In total, we find that 7 sources are classed as short-term variables and only one these sources show variabilitybetween both sub-epochs. MNRAS , 1–18 (2018) J. F. Radcliffe et al.
Based upon our definitions of long and short-term vari-ability, we find that 5/10 variable sources are classed asboth long and short term variables. In these sources, theshort-term variability causes the flux to change sufficientlywhen the sub-epochs are combined such that the sourcesare also classed as long-term variables. We highlight thiseffect in the light curves presented in Figure 5, where thebottom panel shows the influence of the short term variabil-ity has upon the long term variability measurements. Wecan partially mitigate this in our sample by isolating thosesources identified as long term variables only. This resultsin only three sources. However, it is worth stating that thesub-epoch cadence scales sampled here ( < week) cannot ruleout that the long term variability in this sample is beinginduced by > week timescale variability. We would expectthat the physical mechanism causing variability on shorttimescales to be different to the longer timescales. For exam-ple, Ofek et al. (2011) suggest that most short term variabil-ity ( ∼ week timescales) is extrinsic to the source, whilst longterm variability could be intrinsic (e.g. quiescent, AGN flar-ing). Therefore, in order to mitigate this properly and trulyseparate the physical mechanisms causing variability on dif-ferent timescales, we would require multiple epochs whichhave cadences from day to decade timescales such that anyshort term variability that could induce long term variabilitycan be identified.To summarise, we find a total of 10 unique sources ex-hibiting variability which corresponds to ∼ of the totalpersistent radio source population comprising of 479 sources.Of these 10, we find that 8 are classed as long-term variableswhile 7 are short-term variables. There are 5 sources whichexhibit both short and long term variability for which theflux-changes in the short-term influences the long-term vari-ability. This leaves 3 long term variable candidates, howeverit is worth noting that these sources could still be influencedby > week scale variability.This leads to a few conclusions. Firstly, the 1.4 GHz longterm variable radio sky is relatively sedate, continuing thetrend of previous variability surveys with shorter cadences(e.g. Mooley et al. 2016; Hancock et al. 2016). Secondly, atthese flux densities, long cadence variables seem to containless variable sources than the day to week timescale epochswith the fraction of purely long term variables being lessthan 1% of the persistent radio population. However, we notethat with the relatively small numbers of variable sourcesidentified in this study, strong statistical conclusions cannotbe made. For the purposes of this paper, we define a transient sourceas a source which only appears in a single epoch and satisfiesthe variability criteria (based upon flux density upper limitsderived from epochs where the source is not detected). Wefind no source which satisfies this criterion. We discuss theimplications of this result in Section 5.4. It is worth notinghowever, that by using this definition, a transient sourcecould simply be a variable source at a flux density below thedetection threshold which simply increases in flux during asingle epoch.
J123742.3+621518.2
J123644.4+620121.0
J123623.5+621642.7
J123615.0+620613.0
J123513.6+621350.5 P e a k B r i g h t n e ss [ µ J y b e a m − ] Long-term variability
J123751.2+621919.0
J123719.9+620902.7
J123709.4+620837.5
J123701.1+622109.5
J123642.2+621545.4
J123622.5+620653.8
J123620.2+620844.1
Observation time P e a k B r i g h t n e ss [ µ J y b e a m − ] Short-term variability
Figure 5.
Light curves for the variable candidates. The top panelcorresponds to those which are classified as long-term variablecandidates and the bottom panel corresponds to short-term vari-able candidates. The entries with asterisks are both long and shortcadence candidates. The open markers show the combined fluxesfor the JVLA epochs illustrating how short-term variability canartificially induce long-cadence variability.
While this survey does not have the advantage of an ultrawide field of view which is advantageous for finding tran-sients and variables (see Appendix A of Mooley et al. 2016),it does have the considerable advantage of probing deeperthan other surveys across extremely long decadal time scales.Other surveys at these frequencies suggests that the variablefraction is around 1% or less on minute to year timescales(Croft et al. 2010; Thyagarajan et al. 2011; Mooley et al.2013; Bannister et al. 2011a). However, deeper surveys sug-gest that the variable fraction is maybe higher. Carilli et al.(2003) studied variability in a similar manner to this studyusing a deep single VLA pointing at 1.4 GHz in the Lock-man Hole region. They derived an upper limit of ∼ ofsources above a detection threshold of µ Jy beam − vary-ing by more than 50% over 19 days and 17 month cadences.This is largely consistent with the results presented here.In our short timescale study (3 day cadence) we find no MNRAS000
While this survey does not have the advantage of an ultrawide field of view which is advantageous for finding tran-sients and variables (see Appendix A of Mooley et al. 2016),it does have the considerable advantage of probing deeperthan other surveys across extremely long decadal time scales.Other surveys at these frequencies suggests that the variablefraction is around 1% or less on minute to year timescales(Croft et al. 2010; Thyagarajan et al. 2011; Mooley et al.2013; Bannister et al. 2011a). However, deeper surveys sug-gest that the variable fraction is maybe higher. Carilli et al.(2003) studied variability in a similar manner to this studyusing a deep single VLA pointing at 1.4 GHz in the Lock-man Hole region. They derived an upper limit of ∼ ofsources above a detection threshold of µ Jy beam − vary-ing by more than 50% over 19 days and 17 month cadences.This is largely consistent with the results presented here.In our short timescale study (3 day cadence) we find no MNRAS000 , 1–18 (2018) he variable extragalactic radio sky variable sources above µ Jy beam − have a fractional vari-ability >
50% ( | m | > . ).There is evidence that the number of variable sourcesdetected is a function of sampling cadence. Ofek et al. (2011)carried out a 5 GHz survey across 2.66 square degrees of thesky at low Galactic latitudes which aimed to investigate thevariable population on short ( < day) timescales. For sourcesabove flux densities of 1.8 mJy, they found a high fractionof persistent sources ( ∼ ) were variable at the > σ con-fidence level across days to weeks cadences. This is signifi-cantly larger than the results in other surveys, which theyattribute to other blind surveys being insensitive to fast vari-ability (e.g. scintillation), caused by averaging across mul-tiple observations. In addition, this survey was conductedat 5 GHz, which has a selection bias towards flat-spectrumradio AGN. To this end, Ofek et al. (2011) presented a struc-ture function of the variability (defined as ( S −(cid:104) S (cid:105))/(cid:104) S (cid:105) ) whichshowed that the amount of variability is high on day to weektimescales and is constant above around 10 days (albeit withno constraints on longer timescales). Our results supportthis, because we find that 70% of our variable sources arevariable the shortest timescales (days) which could corre-spond to the high levels of variability suggested by the Ofeket al. (2011) structure function.Mooley et al. (2016) extended upon this using the 3 GHzCNSS 50 deg pilot survey covering around 3500 radiosources. They find that 2.6% of sources vary on 1.5 yeartimescales which is considerably more that the weekly (1%)and monthly (0.8%) cadences which is suggestive that thenumber of variable sources is a function of time up to theyearly cadences. In total, they find that . + . − . % of sourcesare variable on timescales < . years. Ignoring the differencein frequency, our results for the shortest timescales (3 days)are consistent with ∼ . of the sources showing variability.On the longest timescales, Hodge et al. (2013) comparedthe Stripe 82 VLA observation with FIRST and found that ∼
6% of sources in the mJy flux density regime show frac-tional variability above 30% on 7-22 year timescales. Thisis considerably higher than seen in other studies which findthat the fraction of variables is of the order a few percent orless (Croft et al. 2010; Becker et al. 2010; Bannister et al.2011b). Our results presented here seem to agree with thelatter studies as we find a long term variable fraction of ≤ .Bannister et al. (2011a) found some evidence that thereis a peak in radio variability on timescales of between 2000and 3000 days. This was reinforced by Hancock et al. (2016)using a 1.4 GHz study of variability in the Phoenix deepfield. Our results suggest that variability reduces in the longdecadal cadences and as such implies that the peak in vari-ability suggested by Bannister et al. (2011a) and Hancocket al. (2016) may be viable. However, it is worth notingthat this survey is probing a different flux density regime(0.03-2 mJy) to the majority of surveys before (typically > . mJy). In the standard picture, the radio source pop-ulation at the flux densities probed by this paper are tran-sitioning from an AGN to a star-forming galaxy dominatedregime, thus the smaller fraction of variable sources may alsobe an imprint of the change in the nature of the underlyingsource population. MNRAS , 1–18 (2018) J . F . R ad c l i ff ee t a l . Table 2.
A summary of the multi-wavelength properties and AGN classification for the variable candidates. Long and short term variable candidates are denoted by L and S respectivelyin the Var. type column. Checkmarks or bold font indicates that the source is classed as an AGN. Redshifts with 68% upper and lower bounds are photometric. The redshift referenceabbreviations are as follows: W04 - Wirth et al. (2004), S04 - (Smail et al. 2004), B08 - Barger et al. (2008), S14 - Skelton et al. (2014), Y14 - Yang et al. (2014). The star formationrates (SFR) were compiled from Whitaker et al. (2014) who used a combination of IR+UV to derive SFRs. The VLBI checkmarks in brackets correspond to those with - σ VLBIcounterparts coincident with the e-MERLIN positions and not reported in Chi et al. (2013) or Radcliffe et al. (2018). The Donley column corresponds to the IR-AGN classificationtechnique proposed by Donley et al. (2012) and the WISE column corresponds to the Stern et al. (2012) AGN classification technique. The crosses indicate that the sources had the bandsnecessary to evaluate these metrics but were not classified as AGN. Entries with hyphens indicate that the measurement could not be taken as the source was out of the field-of-view ofthe required bands while empty entries indicate that the source was not-detected in the required bands for classification.Source ID Var. type z z ref. Morph. SFR (cid:104) S p (cid:105) q Donley WISE X-rays VLBI e-MERLIN[ M (cid:12) yr − ] [ µ Jy beam − ] [ erg s − ]J123742.33+621518.27 L . B08 2 0.4 103 0.92 ± × . × ( (cid:88) ) (cid:88) J123623.55+621642.73 L . S04 1 419 − . ± . × (cid:88) (cid:88) J123615.02+620613.06 L . + . − . Y14 1 - 86 - (cid:88)
J123751.23+621919.00 LS . + . − . S14 1 5.6 163 × (cid:88) (cid:88) J123709.43+620837.55 LS . W04 1 11.4 134 − . ± . × . × (cid:88) (cid:88) J123644.48+620121.09 LS . + . − . Y14 1 - 151 - - × - -J123622.51+620653.90 LS . + . − . S14 3 55.0 161.0 − . ± . × (cid:88) (cid:88) J123620.27+620844.15 LS . B08 1 3.5 122 × × (cid:88) (cid:88)
J123701.11+622109.55 S . B08 1 6.2 262 × × . × (cid:88) J123642.21+621545.42 S . B08 2 85.7 137 . ± . × × . × (cid:88) (cid:88) M N R A S , ( ) he variable extragalactic radio sky To understand the nature of the variability in these sources,we compiled multi-wavelength data for each of these sourcesutilising the extensive photometry of objects in the GOODS-N field. These data are summarised in Table 2. Redshiftswere compiled for the 10 sources from various catalogues(Wirth et al. 2004; Smail et al. 2004; Barger et al. 2008;Skelton et al. 2014; Yang et al. 2014; Momcheva et al. 2016).Of these 12 redshifts, 6/10 are spectroscopic redshifts andthe rest are photometric redshifts (of which one is uncer-tain). The redshifts range from 0.07 to 1.94 with a medianredshift of 0.96.Optical and near infra-red (NIR) Hubble Space Tele-scope (HST) data (Giavalisco et al. 2004; Skelton et al. 2014)plus CFHT K s band imaging (Wang et al. 2010) were used todetermine the host galaxy morphologies. For the morpholo-gies, we visually classified the sources into four categories: 1 -Early-type / bulge dominated, 2 - late-type / spiral galaxies,3 - irregular or 4 - unclassified.The early type / bulge dominated group are those cir-cular/elliptical extended objects whose surface brightnessdistribution drops towards the edge, while late type galax-ies must have clear spiral features. The irregular categoryencompasses those with clumpy surface brightness distribu-tions and sources which are unclassified are those for whicha morphology cannot be attained due to being faint. Thesesources often have undetected low surface brightness areashence they often appear ‘point-like’.We find that the majority of all variable sources (70%)have early-type morphologies which are typical of radio-loudAGN (e.g. Hickox et al. 2009; Griffith & Stern 2010). Theremainder comprise of two late-type galaxies and one ir-regular galaxy. To establish whether the radio emission isconstrained to the nuclear region we compared these opti-cal images to the VLA, e-MERLIN (Muxlow et al. in prep.)and VLBI positions (Chi et al. 2013; Radcliffe et al. 2018).Note that the e-MERLIN data forms part of the upcominge-MERGE survey (Muxlow et al. in prep.; Thomson et al.in prep.)The VLBI observations alone provide the most com-pelling evidence that the radio emission is caused by anAGN. For extragalactic sources, a VLBI detection impliesbrightness temperatures in excess of K which cannotbe reliably attributed to star-formation processes alone andmust instead be AGN related (Condon 1992; Kewley et al.2001). We find that 80% (8/10) of the variable candidatesexhibit VLBI emission. This includes one source with a 6-7 σ VLBI detection coincident with the e-MERLIN emission,which was not formally reported in Radcliffe et al. (2018)due to their adopted 7 σ detection threshold. Interestingly,J123642+621545 was measured to have a total flux den-sity of 343 µ Jy in the 2004 Global VLBI observations ofChi et al. (2013), almost 2.5 × the VLA flux density. Thishas subsequently reduced in flux density significantly andwas not detected in the 2014 EVN observations of Radcliffeet al. (2018). All of the VLBI positions acquired are consis-tent with nuclear activity, apart from J123742.33+621518.30which seems to be offset by around 0 . (cid:48)(cid:48)
24 from the opticalmaximum. We discuss this source in Section 5.3.For the remaining 2 sources which are not detected byVLBI observations, we use the more sensitive e-MERLIN ob- servations ( ∼ µ Jy beam − ; Muxlow et al. in prep.) to searchfor evidence of nuclear emission in the remaining sources.This resulted in one more source with radio emission coin-cident with the nuclear regions while the remaining sourcewas out of the field-of-view of the e-MERLIN observations.We searched for other indications of AGN activity inthese sources in the infra-red (IR) and X-ray bands. Addi-tional Spitzer
IRAC (3.6-8 µ m Wang et al. 2010; Ashby et al.2013; Yang et al. 2014) and MIPS (24, 70 µ m Dickinson et al.2003; Magnelli et al. 2011) plus
WISE (3-22 µ m , Cutri & etal. 2012) photometry was compiled. The X-ray data was de-rived from the 0.5-7 keV, Chandra exposure with thecatalogue provided by Xue et al. (2016). A 1 (cid:48)(cid:48) search radiuswas used and any counterparts were visually evaluated usingthe higher resolution radio and optical HST images.These multi-wavelength data can be used to further dis-tinguish between AGN or star-formation related emissionmechanisms. For
Spitzer
IRAC NIR bands, we used the IRcolour-colour diagnostics of Donley et al. (2012) to establishwhether the source had AGN signatures. Here the presenceof a dusty AGN torus causes the spectral energy distribu-tion (SED) to represent a power law in the region betweenthe . µ m stellar bump and the 25-50 K emission from colddust heated by star-formation. The Donley et al. (2012) cri-terion classifies a source as an AGN if it satisfies all of thefollowing: x ≥ . y ≥ . y ≥ ( . × x ) − . y ≤ ( . × x ) + . S . µ m > S . µ m ∧ S . µ m > S . µ m ∧ S . > S . where x = log ( S . µ m / S . µ m ) and y = log ( S . µ m / S . µ m ) .In total we found that 8/10 sources have NIR counterparts inthe four Spitzer
IRAC bands with S/N ratios larger than σ in each band. As Fig 6 shows, the majority of these sourcesshow little or no signs of AGN activity in the infra-red withtheir NIR characteristics closely following star-formationdominated emission, as shown by the close correlation be-tween the host morphology and the observed NIR colours.For those sources outside of the Spitzer
IRAC coverage, wecross matched to the WISE all-sky catalogue (Cutri & et al.2012) which yielded one additional match. Due to the poorsensitivity in the W3 and W4 bands (only one detection inband 3 and none in band 4), we used the Stern et al. (2012)criterion ( W1 − W2 ≤ . ) in order to identify AGN. Again,this source showed no IR-AGN signatures.In addition, we use the Spitzer
MIPS 24 µ m flux densi-ties to identify those sources which deviate from the IR-radiocorrelation. Despite there being some contamination effects(e.g. from high redshift AGN, Del Moro et al. 2013), theratio between the 24 µ m , S µ m and 1.4 GHz flux densities, S . can be used effectively distinguish between AGN andstar-formation related activity where systems with an AGNpresent will produce more radio emission than what is ex-pected from star-formation alone (e.g. Appleton et al. 2004;Garn & Alexander 2009; Chi et al. 2013). In addition, thisband was chosen to alleviate the effects of confusion which MNRAS , 1–18 (2018) J. F. Radcliffe et al. plagues the longer wavelength observations. For each sourcewith µ m counterparts (5/10), we calculated q using, q = log ( S µ m /(cid:104) S . (cid:105)) . (7)For those sources with µ m counterparts, we find thatthree show clear radio-excess signatures from AGN activ-ity ( q < ; Del Moro et al. 2013). We also cross-matchedthe sample with the 2Ms Chandra
X-ray observations pre-sented in Xue et al. (2016). In total we find only 4 sourceshave X-ray counterparts all of which are classified as AGNby their selection criteria (see Section 2.3.5 of Xue et al.2016).In light of this, these commonly used multi-wavelengthdiagnostics are not able to identify all these variable sourcesas AGN and three sources rely on the identification of acompact radio component. Indeed based upon the star for-mation rate (SFR) measurements and IR colours, some ofthese systems (e.g. J123642.21+621545.42) only show evi-dence of AGN activity by the existence of a high brightnesstemperature core revealed by the VLBI and e-MERLIN ob-servations. This reinforces our belief that high resolutionradio observations are integral to obtaining a full consensusof AGN activity across cosmic time.In total, we find compelling evidence for AGNactivity in 8/10 sources. Of the remaining sources,J123742.33+621518.27 is a possible supernovae (see Sec. 5.3)whilst J123644.48+620121.08 is outside of the e-MERLINand VLBI field-of-view. However it is worth noting that itis hosted by a typical QSO host, an elliptical galaxy.Almost all of the variable objects are detected withVLBI, illustrating that there is a milliarcsecond scale com-pact radio source producing the variability seen here. Re-cent VLBI surveys have shown that VLBI-detected sourcesare becoming more compact at lower flux densities (Deller &Middelberg 2014; Radcliffe et al. 2018). This is supported bysemi-empirical simulations which require FR-I type sourcesto become core-dominated at low flux densities in order toreconcile the high-frequency ( >
10 GHz ) number counts withthe low-frequency number counts ( < ) (Whittam et al.2017). It could be that the variable population at low fluxdensities traces this population because variable sources willinherently contain a compact radio component. Indeed, thelack of mid-IR AGN signatures (which require a dusty AGNtorus to re-radiate UV photons into the mid-IR bands) inany of these objects further adds to this argument, as FR-Igalaxies are known to not have dusty tori (e.g. van der Wolket al. 2010).For the short term variable objects, the physical originof the variability could be scintillation. This has been pre-viously been suggested as the primary cause of most short-term variability at frequencies below 5 GHz. (e.g. Qian et al.1995; Frail et al. 2000; Ofek et al. 2011). For the long-termvariable sample, the source of variability is most likely to beintrinsic to the source and could be caused by changes inthe AGN jet power which has timescales of years to decades(e.g. Aller et al. 1999; Arshakian et al. 2012) or shocks in thejet (e.g. Woo & Urry 2002; Hovatta et al. 2008). It is worthnoting that, without continual monitoring of the source, theorigin of the variable emission cannot be distinguished eas-ily. For a sparsely sampled cadence of observations as pre-sented here, the physical processes which cause variability − . − . . . . . . (S . µ m / S . µ m ) − . − . . . . . . . l og ( S . µ m / S . µ m ) α = − . α = − . AGN
Donley+2012EarlyLateIrr.M82Mrk231S0Sc . . . . . . . R e d s h i f t Figure 6.
Spitzer
IRAC AGN selection criteria for the vari-able candidate sample. The solid line corresponds to the revisedAGN selection wedge by Donley et al. (2012). The solid line withblack markers corresponds to the IR power law locus in the range − . ≤ α ≤ − . . Overlaid are the predicted SED colours of thestarburst galaxy M82, AGN Mrk231, an Sc type and an S0 typegalaxies from the SWIRE library (Polletta et al. 2006). The tracksare color-coded across a redshift range of 0-3.4 with the perpendic-ular bars corresponding to integer redshift intervals. The markercolors of the variable sample correspond to their redshift and arematched to the track colors while the marker shapes correspondto the host morphology. The AGN selected objects using thesecriteria are those in the top-right of the figure for all selectionmethods on month to year timescales can masquerade as processes onmuch longer timescales. To solve this, more epochs spreadover short and long cadences are needed to distinguish be-tween the two. This source was noted to be unusual due to the offset be-tween the radio position (provided by e-MERLIN) and op-tical emission ( ∼ . (cid:48)(cid:48) ) along with low redshift ( z = . ) ofthe host galaxy. As shown in Figure 7, the VLA1996 datadoes not show the compact component which has appearedin the JVLA2011 epoch. The peak radio luminosity from the2012 JVLA 5.5 GHz observations of Guidetti et al. (2017) is . × W Hz − and the spectral index of − . . This is atypical luminosity for core-collapse radio supernovae, whichcan range between × and . × W Hz − (Weileret al. 2002). The q radio excess measure shows no radio-excess emission which would be indicative of an AGN andthere are no signs of AGN activity in the IR diagnostics.Prior to the emergence of this new compact componentthe VLA1996 observations detected diffuse radio emissionwhich was below the formal 7 σ detection threshold used MNRAS000
IRAC AGN selection criteria for the vari-able candidate sample. The solid line corresponds to the revisedAGN selection wedge by Donley et al. (2012). The solid line withblack markers corresponds to the IR power law locus in the range − . ≤ α ≤ − . . Overlaid are the predicted SED colours of thestarburst galaxy M82, AGN Mrk231, an Sc type and an S0 typegalaxies from the SWIRE library (Polletta et al. 2006). The tracksare color-coded across a redshift range of 0-3.4 with the perpendic-ular bars corresponding to integer redshift intervals. The markercolors of the variable sample correspond to their redshift and arematched to the track colors while the marker shapes correspondto the host morphology. The AGN selected objects using thesecriteria are those in the top-right of the figure for all selectionmethods on month to year timescales can masquerade as processes onmuch longer timescales. To solve this, more epochs spreadover short and long cadences are needed to distinguish be-tween the two. This source was noted to be unusual due to the offset be-tween the radio position (provided by e-MERLIN) and op-tical emission ( ∼ . (cid:48)(cid:48) ) along with low redshift ( z = . ) ofthe host galaxy. As shown in Figure 7, the VLA1996 datadoes not show the compact component which has appearedin the JVLA2011 epoch. The peak radio luminosity from the2012 JVLA 5.5 GHz observations of Guidetti et al. (2017) is . × W Hz − and the spectral index of − . . This is atypical luminosity for core-collapse radio supernovae, whichcan range between × and . × W Hz − (Weileret al. 2002). The q radio excess measure shows no radio-excess emission which would be indicative of an AGN andthere are no signs of AGN activity in the IR diagnostics.Prior to the emergence of this new compact componentthe VLA1996 observations detected diffuse radio emissionwhich was below the formal 7 σ detection threshold used MNRAS000 , 1–18 (2018) he variable extragalactic radio sky VLA1996 JVLA2011 ep2 JVLA2018 -21124681115
HST125W + e-MERLIN (2013) . × .
21” ( − . ◦ ) EVN (VLBI) × . ◦ ) Figure 7.
Type II supernovae / LLAGN candidate. The top row shows the evolution of this source from the VLA1996 to the JVLA2018epochs. The colour scale is identical for each VLA epoch and are in units of µ Jy beam − . The bottom row illustrates the e-MERLIN offset(Muxlow et al. in prep.) from the astrometry corrected HST125W optical image (Grogin et al. 2011; Koekemoer et al. 2011; Skeltonet al. 2014) with the 6.7 σ VLBI detection in the bottom right panel. in this variability study. We hypothesis that this diffuse off-nuclear emission may originate from on going star-formationprocesses in the galaxy hence the diffuse radio emitting HIIregion. Low resolution WSRT imaging (observed in May1999; Garrett et al. 2000) recorded an integrated flux den-sity of ∼ µ Jy , which is far in excess of the flux densi-ties recorded by the later JVLA2011 and 2018 epochs (seeFigure 5). These JVLA epochs show that a new compactradio component dominates over the diffuse emission seenin the VLA1996 epoch. e-MERLIN imaging (observed in2013) (Figure. 7) shows this diffuse emission across the faceof the optical host along with an offset compact compo-nent. This diffuse emission could be the same star-formingregion detected in the VLA1996 data or could be from alow-luminosity AGN jet. The compact source has a 6.7 σ EVN detection (observed in 2014) which corresponds to anestimated brightness temperature of ∼ . × K . This isslightly larger than expected from star-formation processesalone (Condon et al. 1991).The physical origins of this off-nuclear compact radiosource remain unclear but the luminosity, spectral indexand light-curve timescales are broadly consistent with thecharacteristics of a powerful core-collapse radio supernovaewhich has occurred between 1996 and 1999 which now dom-inates the radio emission from the HII region. If this newvariable source is a core-collapse radio supernovae, it wouldbe the most distant ever discovered. Unfortunately, with the information available, it is notpossible to categorically determine if this off-nuclear variablesource is a distant radio supernovae or a variable off-nuclearaccretion powered source such as one component of a wideseparation binary AGN system. However this remains an in-triguing variable source which will be monitored by futureVLBI observations to investigate any expansion of the com-pact radio emitting region. In this survey, we detected no transient sources i.e. thosesources which only appear in one epoch alone. Despite thesmall FoV and therefore limiting constraints that this placeson the transient fraction we estimated an upper limit to thetransient fraction was calculated following the prescriptionof Ofek et al. (2011). If we assume an arbitrary luminosityfunction, a uniform density of transient sources in a Eu-clidean universe and a primary beam that can be reliablyapproximated by a Gaussian, then the surface density, κ , oftransients brighter than flux density, S , can be calculatedusing, κ ( > S ) = κ (cid:0) S / S min , (cid:1) − / , (8) κ = N b ln 22 π r (cid:16) − exp (cid:16) − r ln ( )/ r (cid:17)(cid:17) − (9) MNRAS , 1–18 (2018) J. F. Radcliffe et al. where N b is the number of transients, r HP is the primarybeam HPBW, r max is the maximum radius considered, S min , is the detection limit at r = i.e. the primary beam max-ima. For the case considered here (no transients), the σ limit is three events ( N b = , Gehrels 1986). The tight-est constraints require the largest number of independentepochs. In this condition the limiting detection flux den-sity is constrained by the JVLA2018 sub-epochs which hasa central σ detection limit of 54.4 µ Jy beam − . The primarybeam has a r HP = .
258 deg and transients were searchedto a maximum radius of r max = .
233 deg . Substitutioninto equation 8 and dividing by the number of indepen-dent epochs (5) gives a 2 σ upper limit on the areal den-sity κ ( > . µ Jy beam − ) = .
21 deg − . This can be repeatedto lower flux densities using different independent epochcombinations. These new constraints are shown Figure 8which shows the differential source counts and the limitsexcluded by previous radio surveys and the results shownhere. We have included the differential source counts of per-sistent sources from observational surveys (White et al. 1997;Bondi et al. 2008; Vernstrom et al. 2016; Smolˇci´c et al. 2017)and the Tiered Radio Extragalactic Continuum Simulation(TRECS; Bonaldi et al. 2019) for comparison. Our resultsillustrate that there is no upturn in transient events in the µ Jy flux density regime. This is entirely expected, howeverit is worth stating here that the majority of these surveyshave very different cadence scales, therefore each surveywill have a different sensitivity to very-short lived transientevents. One conclusion that could be made is that there isno significant population below our detection threshold thatare extremely variable on long timescales. This follows ontothe generally accepted paradigm that the µ Jy radio sourcepopulation becomes increasingly dominated by star-forminggalaxies which would not be expected to exhibit strong vari-ability on decadal timescales. We have conducted a study on the variability in theGOODS-N field using 5 epochs of JVLA and VLA dataspread across 22 years. We find a total of 10 variable sources,of which 8 show variability on long (year-decade) timescaleswhile 7 show variability on very short (day) timescales. How-ever, we find that the variability in five of the long term vari-ables is maybe influenced by short term variability, leavingonly 3 sources exhibiting true long-term variability. Nearlyall variable sources (9/10) have peak brightnesses above100 µ Jy beam − , which could be an imprint of transition inthe radio populations from AGN to star-forming galaxies atthis flux density regime.Multi-wavelength and ultra-high resolution radio obser-vations reveal that the variability in almost all of thesesources can be reliably attributed to AGN activity. Thelack of IR-AGN signatures, suggests that the variable sourcepopulation in this regime could be the population of core-dominated FR-I galaxies, without dusty tori, that have beenpostulated in other wide-field surveys and semi-empiricalsimulations (Deller & Middelberg 2014; Whittam et al.2017). The only source without an AGN signature may bea Type-II supernovae which, at z = . , would make it oneof the most distant ever discovered. − − − − − − S [Jy] − − − − S . d N / d S [ J y . s r − ] G OO D S - N M o + C r + M o + B e + T h + G Y + B a + C a + B o + L e + . d e g − . d e g − . d e g − . d e g − . d e g − . d e g − White+97Bondi+08Vernstrom+16Smolˇci´c+17TRECS (Bonaldi+19) 0.8 GHz1.4 GHz3.0 GHz5.5 GHz
Figure 8.
The transient normalised areal densities and upper lim-its with respect to flux density. For comparison the normalised dif-ferential number counts for persistent sources are overlaid (Whiteet al. 1997; Bondi et al. 2008; Vernstrom et al. 2016; Smolˇci´c et al.2017). Over-plotted are the expected total extragalactic differen-tial number counts from the TRECS simulation (Bonaldi et al.2019). The wedges correspond to the various 95% upper limits forsurveys without transient detection and their color corresponds tothe survey frequency. The shaded area shows the phase space ex-cluded by the current surveys and the abbreviations correspondsto the following references: Ca+03 - Carilli et al. (2003), Bo+07 -Bower et al. (2007), Be+15 - Bell et al. (2015) , Mo+13 - Mooleyet al. (2013), Mo+16 - Mooley et al. (2016), GY+06 - Gal-Yamet al. (2006), Ba+11 - Bannister et al. (2011a), Cr+10 - Croftet al. (2010), Le+02 - Levinson et al. (2002).
In light of this, transient and variable studies are a func-tion of many different parameters such as cadence time, fluxdensity, observing frequency, observing time and resolution.This paper highlights how important the range of cadenceshave upon the interpretation of the result. In the future, theseparation of variability as a function of the aforementionedparameters will require multi-epoch wide-field observations.While such surveys have been conducted at the mJy andJy regime (e.g. Bower et al. 2010; Croft et al. 2013; Moo-ley et al. 2016), characterising the µ Jy regime requires im-proved snapshot sensitivity and u v coverage which will beprovided by the next generation of interferometers, such asASKAP, MeerKAT, SKA and ngVLA, and their associatedtransient projects such as VAST (Murphy et al. 2013) andThunderKAT (Fender et al. 2017). ACKNOWLEDGEMENTS
The authors would like to acknowledge the referee JimCondon and the scientific editor for their useful commentsthat have greatly helped focus this paper. J.F.R. acknowl-edges support by the Science and Technologies FacilitiesCouncil (STFC), and the Ubbo Emmius Scholarship of theUniversity of Groningen and the South African Radio As-tronomy Observatory (SARAO) fellowship that funded this
MNRAS000
MNRAS000 , 1–18 (2018) he variable extragalactic radio sky research. A.P.T. acknowledges support from STFC grant(ST/P000649/1). The National Radio Astronomy Observa-tory is a facility of the National Science Foundation oper-ated under cooperative agreement by Associated Universi-ties, Inc. e-MERLIN is a National Facility operated by theUniversity of Manchester at Jodrell Bank Observatory onbehalf of STFC. The European VLBI Network is a joint fa-cility of independent European, African, Asian, and NorthAmerican radio astronomy institutes. Scientific results fromdata presented in this publication are derived from the fol-lowing EVN project codes: EG078, GG053. REFERENCES
Aller M. F., Aller H. D., Hughes P. A., Latimer G. E., 1999, ApJ,512, 601Appleton P. N., et al., 2004, ApJS, 154, 147Arshakian T. G., Le´on-Tavares J., B¨ottcher M., Torrealba J.,Chavushyan V. H., Lister M. L., Ros E., Zensus J. A., 2012,A&A, 537, A32Ashby M. L. N., et al., 2013, ApJ, 769, 80Bannister K. W., Murphy T., Gaensler B. M., Hunstead R. W.,Chatterjee S., 2011a, MNRAS, 412, 634Bannister K. W., Murphy T., Gaensler B. M., Hunstead R. W.,Chatterjee S., 2011b, MNRAS, 418, 2813Barger A. J., Cowie L. L., Wang W.-H., 2008, ApJ, 689, 687Becker R. H., Helfand D. J., White R. L., Proctor D. D., 2010,AJ, 140, 157Bell M. E., Huynh M. T., Hancock P., Murphy T., Gaensler B. M.,Burlon D., Trott C., Bannister K., 2015, MNRAS, 450, 4221Berger E., Kulkarni S. R., Frail D. A., Soderberg A. M., 2003,ApJ, 599, 408Bevington P. R., Robinson D. K., 2003, Data reduction and erroranalysis for the physical sciences; 3rd ed.. McGraw-Hill, NewYork, NY, https://cds.cern.ch/record/1305448
Biggs A. D., Ivison R. J., 2006, MNRAS, 371, 963Bonaldi A., Bonato M., Galluzzi V., Harrison I., Massardi M.,Kay S., De Zotti G., Brown M. L., 2019, MNRAS, 482, 2Bondi M., Ciliegi P., Schinnerer E., Smolˇci´c V., Jahnke K., CarilliC., Zamorani G., 2008, ApJ, 681, 1129Booth R. S., de Blok W. J. G., Jonas J. L., Fanaroff B., 2009,preprint, ( arXiv:0910.2935 )Bower G. C., Saul D., Bloom J. S., Bolatto A., Filippenko A. V.,Foley R. J., Perley D., 2007, ApJ, 666, 346Bower G. C., et al., 2010, ApJ, 725, 1792Carilli C. L., Ivison R. J., Frail D. A., 2003, ApJ, 590, 192Carilli C. L., et al., 2015, preprint, ( arXiv:1510.06438 )Chi S., Barthel P. D., Garrett M. A., 2013, A&A, 550, A68Condon J. J., 1984, ApJ, 287, 461Condon J. J., 1992, ARA&A, 30, 575Condon J. J., 1997, PASP, 109, 166Condon J. J., Huang Z.-P., Yin Q. F., Thuan T. X., 1991, ApJ,378, 65Condon J. J., Cotton W. D., Greisen E. W., Yin Q. F., PerleyR. A., Taylor G. B., Broderick J. J., 1998, AJ, 115, 1693Cornwell T. J., Golap K., Bhatnagar S., 2008, IEEE Journal ofSelected Topics in Signal Processing, 2, 647Croft S., et al., 2010, ApJ, 719, 45Croft S., Bower G. C., Whysong D., 2013, ApJ, 762, 93Cutri R. M., et al. 2012, VizieR Online Data Catalog, 2311Del Moro A., et al., 2013, A&A, 549, A59Deller A. T., Middelberg E., 2014, AJ, 147, 14Dennett-Thorpe J., de Bruyn A. G., 2002, Nature, 415, 57Dewdney P. E., Hall P. J., Schilizzi R. T., Lazio T. J. L. W., 2009,Proceedings of the IEEE, 97, 1482 Dickinson M., Giavalisco M., GOODS Team 2003, in Ben-der R., Renzini A., eds, The Mass of Galaxies at Lowand High Redshift. p. 324 ( arXiv:astro-ph/0204213 ),doi:10.1007/10899892 78Donley J., et al., 2012, The Astrophysical Journal, 748, 142Fender R., et al., 2017, preprint, ( arXiv:1711.04132 )Fong W., et al., 2014, ApJ, 780, 118Frail D. A., Waxman E., Kulkarni S. R., 2000, ApJ, 537, 191Gal-Yam A., et al., 2006, ApJ, 639, 331Garn T., Alexander P., 2009, MNRAS, 394, 105Garrett M. A., de Bruyn A. G., Giroletti M., Baan W. A., SchilizziR. T., 2000, A&A, 361, L41Gehrels N., 1986, ApJ, 303, 336Giavalisco M., et al., 2004, ApJ, 600, L93Griffith R. L., Stern D., 2010, AJ, 140, 533Grogin N. A., et al., 2011, ApJS, 197, 35Guidetti D., et al., 2017, Monthly Notices of the Royal Astronom-ical Society, 471, 210Hales C. A., 2016, invgain v1.3, doi:10.5281/zenodo.163000, https://doi.org/10.5281/zenodo.163000
Hallinan G., et al., 2017, Science, 358, 1579Hancock P. J., Drury J. A., Bell M. E., Murphy T., GaenslerB. M., 2016, MNRAS, 461, 3314Healy F., O’Brien T. J., Beswick R., Avison A., Argo M. K., 2017,MNRAS, 469, 3976Hickox R. C., et al., 2009, ApJ, 696, 891Hodge J. A., Becker R. H., White R. L., Richards G. T., 2013,ApJ, 769, 125Hovatta T., Nieppola E., Tornikoski M., Valtaoja E., Aller M. F.,Aller H. D., 2008, A&A, 485, 51Intema H. T., van der Tol S., Cotton W. D., Cohen A. S., vanBemmel I. M., R¨ottgering H. J. A., 2009, A&A, 501, 1185Jauncey D. L., 1968, ApJ, 152, 647Johnston S., et al., 2008, Experimental Astronomy, 22, 151Kewley L. J., Dopita M. A., Sutherland R. S., Heisler C. A.,Trevena J., 2001, ApJ, 556, 121Kimball A. E., Ivezi´c ˇZ., 2008, AJ, 136, 684Koekemoer A. M., et al., 2011, ApJS, 197, 36Levinson A., Ofek E. O., Waxman E., Gal-Yam A., 2002, ApJ,576, 923Lorimer D. R., Bailes M., McLaughlin M. A., Narkevic D. J.,Crawford F., 2007, Science, 318, 777Lovell J. E. J., et al., 2008, ApJ, 689, 108Magnelli B., Elbaz D., Chary R. R., Dickinson M., Le Borgne D.,Frayer D. T., Willmer C. N. A., 2011, Astronomy & Astro-physics, 528, A35McMullin J. P., Waters B., Schiebel D., Young W., Golap K.,2007, in Shaw R. A., Hill F., Bell D. J., eds, AstronomicalSociety of the Pacific Conference Series Vol. 376, AstronomicalData Analysis Software and Systems XVI. p. 127Mohan N., Rafferty D., 2015, PyBDSF: Python Blob Detec-tion and Source Finder, Astrophysics Source Code Library(ascl:1502.007)Momcheva I. G., et al., 2016, ApJS, 225, 27Mooley K. P., Frail D. A., Ofek E. O., Miller N. A., KulkarniS. R., Horesh A., 2013, ApJ, 768, 165Mooley K. P., et al., 2016, ApJ, 818, 105Morrison G. E., Owen F. N., Dickinson M., Ivison R. J., Ibar E.,2010, ApJS, 188, 178Murphy T., et al., 2013, Publ. Astron. Soc. Australia, 30, e006Muxlow T. W. B., et al., 2005, MNRAS, 358, 1159Ofek E. O., Frail D. A., Breslauer B., Kulkarni S. R., ChandraP., Gal-Yam A., Kasliwal M. M., Gehrels N., 2011, ApJ, 740,65Offringa A. R., van de Gronde J. J., Roerdink J. B. T. M., 2012,A&A, 539, A95Offringa A. R., et al., 2014, MNRAS, 444, 606MNRAS , 1–18 (2018) J. F. Radcliffe et al.
Oosterloo T., Verheijen M., van Cappellen W., 2010, inISKAF2010 Science Meeting. p. 43 ( arXiv:1007.5141 )Owen F. N., 2018, ApJS, 235, 34Owen F. N., Morrison G. E., 2008, AJ, 136, 1889Padovani P., Bonzini M., Kellermann K. I., Miller N., MainieriV., Tozzi P., 2015, MNRAS, 452, 1263Perley R. A., Butler B. J., 2013, ApJS, 204, 19Perley R. A., Butler B. J., 2017, ApJS, 230, 7Perley R. A., Chandler C. J., Butler B. J., Wrobel J. M., 2011,ApJ, 739, L1Pietka M., Fender R. P., Keane E. F., 2015, MNRAS, 446, 3687Planck Collaboration et al., 2016, A&A, 594, A13Polletta M. d. C., et al., 2006, ApJ, 642, 673Qian S. J., Britzen S., Witzel A., Krichbaum T. P., Wegner R.,Waltman E., 1995, A&A, 295, 47Radcliffe J. F., et al., 2018, A&A, 619, A48Rau U., Cornwell T. J., 2011, A&A, 532, A71Richards E. A., 2000, ApJ, 533, 611Skelton R. E., et al., 2014, ApJS, 214, 24Smail I., Chapman S. C., Blain A. W., Ivison R. J., 2004, ApJ,616, 71Smolˇci´c V., et al., 2017, A&A, 602, A1Stach S. M., et al., 2018, ApJ, 860, 161Stern D., et al., 2012, ApJ, 753, 30Thyagarajan N., Helfand D. J., White R. L., Becker R. H., 2011,ApJ, 742, 49Vernstrom T., Scott D., Wall J. V., Condon J. J., Cotton W. D.,Perley R. A., 2016, MNRAS, 461, 2879Wang W.-H., Cowie L. L., Barger A. J., Keenan R. C., Ting H.-C.,2010, ApJS, 187, 251Weiler K. W., Panagia N., Montes M. J., Sramek R. A., 2002,ARA&A, 40, 387Whitaker K. E., et al., 2014, ApJ, 795, 104White R. L., Becker R. H., Helfand D. J., Gregg M. D., 1997,ApJ, 475, 479Whittam I. H., Jarvis M. J., Green D. A., Heywood I., Riley J. M.,2017, MNRAS, 471, 908Williams P. K. G., Bower G. C., Croft S., Keating G. K., LawC. J., Wright M. C. H., 2013, ApJ, 762, 85Wirth G. D., et al., 2004, AJ, 127, 3121Woo J.-H., Urry C. M., 2002, ApJ, 579, 530Xue Y. Q., Luo B., Brandt W. N., Alexander D. M., Bauer F. E.,Lehmer B. D., Yang G., 2016, ApJS, 224, 15Yang G., et al., 2014, ApJS, 215, 27van Haarlem M. P., et al., 2013, A&A, 556, A2van der Wolk G., Barthel P. D., Peletier R. F., Pel J. W., 2010,A&A, 511, A64This paper has been typeset from a TEX/L A TEX file prepared bythe author. MNRAS000