© ESO 2018

1. Introduction

Radio source counts above mJy flux densities are dominated by radio galaxies and quasars powered by active galactic nuclei (AGN). Below mJy flux densities, there is an observed upturn far in excess of those predicted by extrapolating source counts of high luminosity radio galaxies and quasars. This upturn is found to comprise an increasing fraction of active star forming galaxies and faint “non-jetted” or radio-quiet AGN plus a decreasing fraction of classical radio-loud sources (see Prandoni et al. 2001; Huynh et al. 2015; Padovani 2016, and references therein). The majority of extragalactic radio surveys are carried out at arc-second resolutions (corresponding to galactic/∼10’s kpc physical scales at z≥ 0.1) where it can be difficult to distinguish between the sub-kpc scale AGN activity and the kpc star-formation related emission based purely on their radio morphologies. This is particularly important if we are to characterise the properties of radio-quiet AGN whose radio emission in local systems are confined within the host galaxy (see Orienti et al. 2015, and references therein). As a result, these surveys rely on multi-wavelength diagnostics, such as radio-excess, SED fitting, X-ray emission etc., in order to identify any AGN activity (e.g. Bonzini et al. 2013; Smolčić et al. 2017b). These diagnostics are often incomplete with dust masking the signatures of AGN activity. For example, X-rays often do not detect Compton-thick AGN which are estimated to account for over a third of the total AGN population (Mateos et al. 2017).

These hidden AGN can be found using high resolution, dust-independent radio observations. Indeed, surveys using e-MERLIN, such as the e-MERlin Galaxy Evolution (e-MERGE) survey (Muxlow et al. 2005, and in prep.), and Very Long Baseline Interferometry (VLBI; e.g. Middelberg et al. 2011, 2013; Herrera Ruiz et al. 2017) have shown that deep, sub-arcsecond and sub-kpc observations can effectively isolate AGN activity from compact star-forming related emission in distant galaxies.

VLBI observations detect bright, compact objects with brightness temperatures in excess of 10⁵ K. In nearby galaxies, these brightness temperatures can be typically reached by either AGN, supernovae (SNe) and their remnants (SNRs). However, in more distant galaxies (z > 0.1), these brightness temperatures can typically only be attained by AGN-related emission processes (e.g. Kewley et al. 2000), thus making VLBI a unique and invaluable tool to survey distant galaxies for AGN activity. However, until the last decade, there have been many factors preventing VLBI from being used as a survey instrument.

Conventional wide-field VLBI observations mapped a significant proportion of the primary beam by using a single correlation pass at a ultra-fine temporal and frequency resolution in order to limit time and bandwidth smearing towards the edge of the primary beam (Garrett et al. 2001). As a result, the observer would receive a single large and unwieldy (often ∼TB size) data set. With the ever increasing number of VLBI-ready telescopes along with widening bandwidths, the bit rates of modern VLBI arrays are rapidly increasing and this method of correlation has become computationally infeasible. Software correlators established the “multiple simultaneous phase centre observing” approach to correlation (Deller et al. 2011; Morgan et al. 2011; Keimpema et al. 2015) which substantially reduces the computational load. Here, the observer defines a number of sub-fields (also referred to as phase centres) which can be either sources of interest or can be arranged to cover the entire primary beam. When correlated, these data are split and each sub-section is correlated at the ultra high temporal and frequency resolution required to restrain smearing. It is then copied and phase shifted to the various sub-fields of interest and averaged to a small field-of-view (typically 30–60″). The result is a small (∼GB) dataset per sub-field which is easily manageable and parallelisable when calibrating and imaging.

By combining multi-phase centre correlation with advanced calibration techniques such as in-beam phase referencing (Garrett et al. 2001, 2005; Lenc et al. 2008) and multi-source self-calibration (Middelberg et al. 2013; Radcliffe et al. 2016), wide-field VLBI surveys of milliarcsecond scale extragalactic radio sources to μJy flux densities have become increasing feasible (e.g. Middelberg et al. 2011, 2013; Chi et al. 2013; Morgan et al. 2013; Cao et al. 2014; Deller & Middelberg 2014; Rampadarath et al. 2015; Herrera Ruiz et al. 2017).

We here present a new wide-field VLBI survey targeting the well studied Great Observatories Origin Deep survey North (GOODS-N) field using the European VLBI Network (EVN). The GOODS-N field covers 160 arcmin² with complementary deep multi-wavelength data including Chandra, Spitzer, Herschel, UBVRIJHK photometry and spectroscopy.

Previous wide-field VLBI observations targeted the Hubble Deep Field (HDF) and Flanking Fields (HFF) for which the GOODS-N field encompasses. Garrett et al. (2001) used the EVN to target MERLIN sources within a 3.5″ radius from the EVN pointing centre to rms sensitivities of 33 μJy beam⁻¹.

This resulted in the detection of 3 sources. These observations were substantially expanded upon by Chi et al. (2013) who used Global VLBI to target the 92 VLA-MERLIN sources of Muxlow et al. (2005) within a 10″ × 10″ field to rms sensitivities of 7.3 μJy beam⁻¹. This resulted in 12 compact radio source detections (including the 3 detected by Garrett et al. 2001) thus beginning the characterisation of the faint compact radio population in GOODS-N. However, these surveys were invariably limited because computational limitations prevented imaging of the entire primary beam at that time. Our survey aims to substantially expand upon this sample, encompassing and surpassing the field-of-view and sensitivities of previous VLBI surveys in GOODS-N by targeting sources within a 30″ × 30″ area to 1σ central sensitivity of ∼2−3 μJy beam⁻¹ with the completion of this survey.

In this paper, we present our initial catalogue of the 31 compact sources detected in the first data release to a 1σ sensitivity of ∼9 μJy beam⁻¹ (corresponding to ∼17.5 h on source) along with derived radio properties of these objects using complementary 1–2 GHz VLA data. In paper II, we compare our VLBI-selected population to other AGN detection diagnostics used in other wavebands. A future publication, paper III, will describe the final data release which will include an additional 48 h of observations which comprise of the first wide-field VLBI observations using a combined eMERLIN-EVN array.

For this paper, we adopt a spatially-flat 6-parameter ΛCDM cosmologywithH₀ = 67.8 ± 0.9kms⁻¹ Mpc, Ω_m = 0.308±0.012 and Ω_Λ = 0.692 ± 0.012 (Planck Collaboration XIII 2016). We assume S_ν ∝ ν^α throughout, where S_ν is the radio integrated flux density and α is the intrinsic source spectral index.

The paper is organised as follows. Section 2 outlines our observations, source selection strategy, calibration and source detection methodology. Section 3 details the primary beam correction method used for the EVN. Section 4 describes the VLBI catalogue accompanying this paper while a formatted version is presented in Tables 2 and 3. Section 5 presents our results and associated discussion, including redshifts, astrometry, comparisons with other VLBI surveys and the radio properties of the VLBI-selected population. We conclude our findings in Sect. 6.

2. Observations and data reduction

2.1. EVN observations

The EVN observed the GOODS-N field at 1.6 GHz for 24 h in total on the 5–6th June 2014 (EVN code EG078B). The pointing centre used was the Hubble Deep Field-North (HDF-N), J2000 12:36:50.0 62:12:58.0. Ten telescopes, including the 100 m Effelsberg and the 76 m Lovell (Jodrell Bank) telescopes, were involved in the observation. In order to attain a uniform sensitivity profile, the Effelsberg and Lovell telescopes were nodded between 5 different pointing centres over the course of the observation while smaller telescopes remained pointed at the HDF-N centre. The observing strategy and participating telescopes are summarised in Table 1.

These data were recorded at a bit rate of 1024 Mbits s⁻¹ (8 × 16 MHz bands) in both right and left hand circular polarizations. The fringe finders used were 3C345 and DA193. The observations were made using the standard phase referencing mode. Two phase calibrators were used; a strong, ∼0.4 Jy, primary calibrator J1241+602 lying approximately 2^◦ from the target centre, and a weaker, 17 mJy, secondary calibrator J1234+619 lying 23.5″ from the target centre. The primary calibrator was observed for 1.5 min on source every ∼27 min. To permit more accurate phase corrections, the secondary calibrator was observed more frequently (1.5 min on source every ∼7.5 min). In total, the on-source integration time on the GOODS-N field was approximately 17.5 h.

2.1.1. Source selection and correlation

These data were correlated using the SFXC correlator (Keimpema et al. 2015) at the Joint Institute for VLBI ERIC, Dwingeloo, the Netherlands. The correlation implemented the multiple simultaneous phase centre observing technique (see Deller et al. 2007, 2011; Keimpema et al. 2015) to target 699 sub-fields/phase centres. Two source selection strategies were implemented and the criteria are as follows:

• (a)

  A survey designed to completely cover the central 7′̣5 radius area with VLBI coverage whilst restraining bandwidth and time smearing to <10%. This is designed to complement the 1.5 GHz eMERGE survey (Muxlow et al., in prep.) for use in integrated imaging. This comprises of 582 phase centres split into three main categories:

  • (i)

    339 1.5 GHz eMERLIN-VLA sources complete to ∼11 μJy beam⁻¹ from the upcoming eMERGE survey (Muxlow et al., in prep.).

  • (ii)

    151 SCUBA-2 850 μm sources (Smail, priv. comm.).

  • (iii)

    92 other positions, denoted “legacy”, which cover gaps in coverage across the field.

• (b)

  A targeted survey of 117 1.4 GHz VLA radio-bright sources outside the central 7′̣5 radius area with integrated VLA flux densities, S_{i, VLA} > 0.1 mJy (Morrison et al. 2010). These are split into:

  • (i)

    91 sources with 0.1 < S_{i, VLA} < 1 mJy targeted to a radius of 12″ from the central pointing centre.

  • (ii)

    26 sources with S_{i, VLA} > 1 mJy targeted to the edge of the Morrison et al. (2010) VLA survey. We note that 4/26 of these sources are within the 15″ and were accidental double entries in the correlation catalogue. These are kept in for clarity in-case these data are re-reduced in the future.

The source positions targeted are shown in Fig. 1. With source positions determined, correlation proceeded as follows. Short sub-integrations of data were correlated at the required high spectral and temporal resolutions in order to reduce time and bandwidth smearing. In this observation, each sub-integration had a frequency resolution of 1.953 kHz and a time resolution of 13.056 ms in order to restrain time and bandwidth smearing to below 1% on the longest baseline (∼8400 km) at 5″ from the pointing centre. At the end of each sub-integration, the visibilities were phase shifted to every desired source position to create a separate data set per position. Each data set was averaged to a temporal resolution of 3 s and a frequency spacing of 0.5 MHz (corresponding to a 10% time and bandwidth smearing at 30″ from the assigned source position) and then added to previous sub-integrations until the entire data set was correlated. As a result, 699 separate, narrow field-of-view (FoV) data sets were produced, one per source position. Attached to one data set, containing source J123462+621331, were the scans of the phase calibrators, J1241+602 and J1234+619, and the fringe finders, 3C345 and DA193 used for calibration. The phase referencing calibration and flagging tables derived for this data set can then be easily copied to the other data sets. Despite the total size being 3.79 TB, calibration of this data set is easily parallelised and computationally less intensive than previous wide-field VLBI projects (e.g. Chi et al. 2013). For further clarification, we will refer to the target field as the GOODS-N field as a whole, while the sub-fields are the small FoV phase centres within the GOODS-N field whose coordinates were correlated upon.

2.1.2. Data reduction

These data were reduced using the Astronomical image processing (AIPS) software developed by NRAO¹ (Greisen et al. 2003), and its Python interface, Parseltongue (Kettenis et al. 2006).

Before describing the data reduction, we note that there was an error found in the position of the secondary phase calibrator (J1234+619) when we tested phase referencing from the primary calibrator (J1241+602) to the secondary phase calibrator. This incorrect position originated from the Chi et al. (2013) observations. The correct position was found to be J2000 12:34:11.7413(57) +61:58:32.478(07). Independent 5 GHz e-MERLIN observations of J1234+619 using multiple phase calibrators verified that this new position is correct (see Appendix A). The model of J1234+619 derived from initial phase referencing tests were then used when fringe fitting in order to ensure the correct position is used.

With this issue established and solved, these data were calibrated as follows. Gains were calibrated using the system temperature, T_sys, measurements from each antenna and the data were edited to remove any Radio Frequency Interference (RFI) using the AIPS tasks SPFLG and CLIP. Instrumental phase offsets between the spectral windows² were then removed by solving for the phase and delays on a two minute integration of 3C345, using the task FRING. This allowed the spectral windows to be combined when the data is fringe fitted. We note that the dispersive delays were not corrected for, however we are confident that this is a minimal contribution as the fully calibrated data shows no phase deviations across the frequency band in excess of 10–15 degrees on all baselines.

The group delays for the phase calibrators and fringe finders were calibrated using FRING (using a model of J1234+619 created when investigating the positional offsets), edited (using SNEDT) and smoothed (using SNSMO) to remove noisy and spurious solutions. With the delays calibrated, the phase and rates were then calibrated, edited and smoothed and applied to the data. We note that fringe fitting was conducted on both phase calibrators. This is because the Lovell Telescope (Jb) did not observe the primary phase calibrator due to a restriction on the number source changes per hour. With fringe fitting complete, the bandpass response was calibrated using AIPS task BPASS. 3C345 was used for bandpass calibration on all telescopes

The primary phase calibrator (and furthest from the target), J1241+602, underwent three rounds of phase only self calibration (with solution intervals of five, one, and one minute(s) respectively) and one round of amplitude and phase self-calibration with a five minute solution interval. These solutions were applied to the closer, secondary phase calibrator, J1234+619. Three rounds of phase only self-calibration (with solution intervals of five, four, and two minutes respectively) were conducted on J1234+619 and these solutions were then applied to the sub-field containing J123462+621331. Spectral windows were not combined when self-calibration was performed on J1241+602, or during the first round of self-calibration on J1234+619. This would effectively correct any residual dispersive delay errors caused by a variable ionosphere by approximating the true phase correction (smoothly variable with frequency) with one that is a stepwise constant (one value per subband per solution interval).

The calibration solutions and flagging tables derived and applied to J123642+621331 were then applied to the other 698 sub-fields. All sub-fields were then imaged using AIPS task IMAGR using both natural (UVWTFN=“NA”) and uniform weighting schemes and these images were searched for emission. A detection threshold of 6σ was used to reduce the chance of false positives (see Radcliffe et al. 2016, for more in depth discussion).

To reduce residual phase errors arising from atmospheric inhomogeneities between the phase calibrator and target field, we utilised the Multi-source Self-calibration (MSSC) technique developed by Middelberg et al. (2013) and Radcliffe et al. (2016). The nine brightest sources were used in MSSC. These sources were detected when imaged with both uniform and natural weighting schemes. If a source was detected in both images, it is highly suggestive that the source can be detected on all baselines. As a conservative precaution, we excluded sources outside the primary beam of the large telescopes (Effelsberg and Lovell) because we would expect considerable phase and gain errors to be induced by the attenuation of the primary beam. These errors would not be the same for each sub-field and will simply add noise into the solutions for MSSC when the individual sub-fields are combined. By performing this we ensure that the dominant error is from differing atmospheric paths between the phase calibrator and the target field. Three rounds of phase-only MSSC were performed using a solution interval of 2 min and the corrections derived were applied to all sub-fields. A primary beam correction scheme (outlined in Sect. 3) was then applied to the central 7′̣5 radius field and the sub-fields were searched for emission. The following section describes the method used to detect sources once MSSC and the primary beam correction was applied.

2.2. Source detection methodology

To determine accurate peak brightnesses and integrated flux densities, we tested multiple source detection algorithms namely AIPS task SAD, BLOBCAT (Hales et al. 2012) and PYBDSF (Mohan & Rafferty 2015). It has been noted that Gaussian fitting routines, namely SAD and PYBDSF, were found to routinely over estimate integrated flux densities in the low signal-to-noise (S/N) regime where noise fluctuations across the extent of a source can induce sub-optimal fitting (see Middelberg et al. 2013). In the low S/N regime (S/N∼ 6−10), the measured integrated flux densities were on average ∼16% and ∼22% higher than the BLOBCAT measured values when using SAD and PYBDSF respectively. In the high S/N regime, the effects are less pronounced with both SAD and PYBDSF reporting fluxes only 4–5% larger than BLOBCAT.

We therefore used BLOBCAT to measure peak brightnesses and integrated flux densities of our sources. Initially, BANE (Hancock et al. 2012) was used to generate a rms map of each field which is in turn input to BLOBCAT. All parameters were set to default apart from the surface brightness error which was assumed to be ∼10% which is caused by amplitude calibration errors (--pasbe=0.1), the minimum S/N detection threshold (--dSNR=6) and, as our point spread function (psf or synthesised beam) is vastly oversampled, the peak brightness pixellation error was set to 1% (--ppe=0.01; see Appendix A of Hales et al. (2012) for more information). The surface brightness error included an additional error factor which is proportional to the value of the primary beam correction. Note that BLOBCAT does not provide any size information therefore source sizes were measured using PYBDSM³.

For the central 7′̣5 radius field, each sub-field was imaged using natural weighting only (to optimise sensitivity) and then the method outlined above was used to search for detections. Detections were then imaged with uniform weighting (whose rms is 1.6 × the naturally weighted rms) in order to obtain the highest resolution image possible.

To optimise the number of detections, we used a different strategy to image sub-fields outside the central 7′̣5 radius field which are beyond the half-power beam widths (HPBWs) of the large telescopes (Effelsberg and Lovell). Two images were made for each sub-field. For the first image, the large telescopes on all pointings apart from the closest pointing centre to the sub-field were flagged. This was performed because while the large telescopes retain significant sensitivity well beyond the HPBW of their primary beam, the retention of the more distant pointing centres would induce significant amplitude errors which would out-weigh any sensitivity gains. This method produced an additional 6 detections, the majority of which (4/6) are within a 12″ radius of the pointing centre. For the second image, all of the large telescopes are flagged, so that sources up to the HPBW of the smaller 32 m and 25 m telescopes could be detected without being affected by amplitude errors from the large telescopes still present in these data. This method produced just one additional detection and none of the 6 sources detected with the larger telescopes included were detected with this method. This is likely due to the significant sensitivity reduction when the large telescopes are removed. Additionally, we note that the primary beam models are poorly constrained outside the HPBW, therefore these sub-fields do not have primary beam correction applied.

Once detections were identified, each sub-field was re-imaged with both uniform weighting (∼5.3 × 4.5 mas) and natural weighting (∼16 × 16 mas) schemes and re-catalogued resulting in a total of 31 detections (24 from the central field and 7 from the targeted survey beyond the 7′̣5 radius field). These detections are shown in Fig. 3 which is an update on Fig. 1 from Chi et al. (2013). We note that Radcliffe et al. (2016) only stated an initial 20 sources, however this study invariably missed detections because the majority of the annulus sub-fields were not included. The derived peak brightnesses, flux densities and positions of our objects are described in Table 2.

2.3. VLA observations

In addition to the EVN observations, archival L-band Karl. G. Jansky Very Large Array (VLA) A-array data (PI F. Owen) were reduced to provide a zero-spacing flux density for our VLBI sources and forms part of the eMERGE survey data (Muxlow et al., in prep.). This section briefly describes the data reduction process. The VLA observed the GOODS-N field between the 7th August–11th September 2011 for a total of 38 h, in the A-array configuration. These data were flagged using the AOFlagger software (Offringa et al. 2012) and calibrated using the VLA CASA calibration pipeline⁴ (packaged with CASA version 4.3.1). These data were transferred to AIPS and two sources (J123452+620236 and J123538+621932) were peeled.

Due to the large fractional bandwidth (∼68%) and large data size, postage stamps centred on each VLBI detection were produced using the multi-term multi frequency algorithm within CASA task tclean (Rau & Cornwell 2011). These images were primary beam corrected using the CASA routine widebandpbcor, which also corrects for the induced spectral index caused by the varying primary beam attenuation across the bandwidth. The resulting images have a rms of ∼2−5 μJy beam⁻¹ with a restoring beam of 1.54″ × 1.34″. Flux densities were extracted using BLOBCAT (Hales et al. 2012) and we conservatively assign a standard 10% surface brightness error originating from calibration.

3. EVN primary beam correction

For these observations, we used and developed one of first primary beam models European VLBI Network (Keimpema et al., in prep.). We followed a similar prescription to primary beam modelling as described in Strom et al. (2004) and Cao et al. (2014). Due to the lack of accurate primary beam models for many EVN telescopes, the primary beam power response of each telescope can be approximated by using a normalised, symmetric, 2D Gaussian of the form,

(1)

where P(θ, φ) is the relative primary beam power response. θ and φ are the respective azimuthal and polar angular distances from the antennas’ pointing centres. The azimuthal and polar coordinates of the telescope’s pointing centres are defined by θ₀ and φ₀ respectively. The standard deviation, σ, can be related to the FWHM of the primary beam, θ_1/2, through the expression,

(2)

where the FWHM of the primary beam is defined as,

(3)

Here, λ_c, is the observing wavelength, and D is the aperture diameter (Keimpema 2015, unpublished). A small correction factor, K = 1.05, was used to take into account any aperture blockages (Wrigley et al., in prep.). For some telescopes (namely Effelsberg, Jodrell Bank and Shenzhen in these observations) accurate beam models are available and therefore replace D/K in Eq. (3) with the fitted HPBW of these telescopes. We note that there are some uncertainties associated with the derived HPBW of Jodrell Bank because this is derived from a modelled aperture distribution rather than more accurate methods such as holographic scans (Wrigley 2016). The fitted aperture diameters are summarised in Table 1.

P(θ, φ)^−1/2, that is the inverse primary beam voltage response, was calculated for each telescope at every time integration step, each spectral window and each sub-field. These were recorded into an AIPS SN table (one per sub-field) which were then be applied to the corresponding uv data set using AIPS task CLCAL. The application of this calibration table multiplies the visibility amplitudes of each baseline by a correction factor which is the inverse product of the primary beam voltage responses, (P_i(θ, φ)^−1/2 × Pj(θ, φ)^−1/2) of the two telescopes i, j that form the baseline. The application of this table also adjusts the weights by the inverse of this correction factor. The simultaneous adjustment of weights means this primary beam correction will also correctly weight multiple pointings, thus permitting mosaicking using the EVN array.

We note that the use of a Gaussian model for our primary beam model implies heavy tapering of EVN telescopes, however these telescopes were designed for single dish observations so we would expect a large amount of tapering. In addition, we tested multiple models (1D Gaussian, 2D Gaussian, Airy disk, and a polynomial) when fitting to these beam models and found that differences between the models are only significant towards the primary beam null and these models are indistinguishable within the HPBW.

However, because of the lack of available primary beam models for many EVN telescopes (most notably the Lovell telescope for this observation), there are considerable uncertainties on the derived peak brightnesses and integrated fluxes. Incorporating this model to include higher order corrections, such as beam squint and more physically motivated models would be simple. Nevertheless, if we are to significantly improve EVN primary beam correction, beam measurements for all individual EVN stations are required, because uncertainties towards the edge of the primary beam are dominated by the lack of information on the primary beams’ HPBWs.

As a result of this, sub-fields which lie outside the central 7′̣5 radius area were not primary beam corrected because, the Effelsberg and Lovell telescope corrections are beyond the FWHM of the Gaussian models derived, and errors due to incorrect beam models will rapidly increase. Figure 2 shows the rms sensitivity of our observations after primary beam correction using natural weighting. The central rms is approximately 9 μJy beam⁻¹.

We note that this primary beam correction is constantly under-development and will be updated with the latest EVN beam models. The code is publicly available⁵.

4. Catalogue description

In this section we describe the VLBI catalogue of 31 compact radio sources which accompanies this paper (see Tables 2 and 3). The designated column numbers correspond to the associated machine-readable version of the catalogue. The column descriptors are as follows:

Column 1: Source ID. Radio name adopted in this paper which is of the form Jhhmmss+ddmmss based upon the J2000 Right Ascension (in hours) and Declination (in degrees). Note that some source identifiers are slightly different to that of Morrison et al. (2010) because of the improved astrometric accuracy.

2–4: z. The redshifts for the 31 VLBI detected sources can be found in Col. 2. The 68% lower and upper confidence intervals are in Cols. 3–4. A description of how these redshifts were compiled can be found in Sect. 5.1.

5,6: z type/ref. Column 5 describes how these redshifts were determined. Spectroscopic redshifts are denoted with “S” in the z type column while photometric redshifts are denoted with “P”. Column 6 contains the reference for which the redshift was acquired

7: RA (J2000) Right Ascension (J2000) in hours with the format hh:mm:ss.ssss.

8: Dec (J2000) Declination (J2000) in degrees with the format dd:mm:ss.sss.

9–10: VLBI P. Peak brightness (Col. 9) and associated error in μJy beam⁻¹ (Col. 10). This is determined using the flood-filling algorithm BLOBCAT (Hales et al. 2012). Errors derive from source fitting, calibration, and primary beam correction. Peak brightnesses of those sources without primary beam correction are merely lower limits (denoted by −99.0 in the error column in the machine-readable table).

11–12: VLBI I. Integrated flux density (Col. 11) and associated error in μJy (Col. 12). This is determined when deriving the peak brightness. Errors originate from source fitting, calibration, and primary beam correction. Integrated flux densities of those sources without primary beam correction are merely lower limits (denoted by −99.0 in the error column in the machine-readable table).

13: S/N. Signal-to-noise ratio.

14–16: Beam. Major axes (Col. 14), minor axes (Col. 15) and position angle (Col. 16) of the CLEAN restoring beam in milliarcseconds. The two restoring beams of ∼16 × 16 mas and ∼5.3 × 4.5 mas correspond to natural weighting and uniform weighting schemes respectively. The VLBI peak brightnesses and integrated flux densities presented in Cols. 9–10 and 11–12 respectively have been derived using the beam sizes specified in these columns.

17–18: VLA P. Peak brightness (Col. 17) and associated errors (Col. 18) of the 1.5 GHz VLA observations described in Sect. 2.3. Errors are determined from calibration (∼10%) and source fitting.

19–20: VLA I. Integrated flux density (Col. 19) and associated errors (Col. 20) of the 1.5 GHz VLA observations. Errors originate from calibration (∼10%) and source fitting.

21: α. Spectral index between the 5.5 GHz VLA observations of Guidetti et al. (2017) and 1.5 GHz VLA observations presented in this paper.

22–23: L_1.5GHz. Monochromatic rest-frame radio luminosity in WHz⁻¹ (Col. 22) and its associated error (Col. 23). Values were derived using the 1.5 GHz VLA integrated flux densities and k-corrected using the spectral index. Associated errors originate primarily from those sources with photometric redshift errors. The median spectral index of −0.56 were used to determine the radio luminosities of sources without spectral index information.

24–26: T_b. Brightness temperatures, in K, calculated using Eq. (5) (Col. 24). Column 25 is a flag to denote lower limits (set as 1 in machine readable version to denote lower limits) and Col. 26 corresponds to whether the brightness temperature was calculated using uniform or natural weighting. This is denoted as U or N respectively in the machine-readable table and in Table 3 brightness temperatures derived using natural weighting are italicised. See Sect. 5.4.3 for further details.

27–30: Angular sizes. Columns 27 and 29 describe the major and minor axes of the deconvolved projected angular size of the VLBI source in milliarcseconds. Columns 28 and 30 contain flags to denote upper limits (hence unresolved source sizes) for the major and minor axes respectively (set as 1 in machine readable version to denote upper limits). Sizes were fitted using PYBDSF (Mohan & Rafferty 2015) and see Sect. 5.4.3 for further details. 31–34: Linear size. Columns 31 and 33 describe the major and minor axes of the deconvolved projected linear size of the VLBI source in parsecs and Cols. 32 and 34 contain flags to denote upper limits for the major and minor axes respectively (set as 1 in machine readable version to denote upper limits).

5. Results and discussion

5.1. Redshifts

The VLBI positions were matched to the many spectroscopic and photometric catalogues to within a radius of one arcsecond from the VLBI position. In order to prevent misidentifications potential matches were visually compared to the HST optical/near-IR images of Skelton et al. (2014) to ensure the correct redshift is assigned. Redshift information was found for all 31 objects. This comprises of 20 spectroscopic redshifts (Cowie et al. 2004; Smail et al. 2004; Chapman et al. 2005; Barger et al. 2008; Murphy et al. 2017; Cowie, priv. comm.) and 11 photometric redshifts (Skelton et al. 2014; Yang et al. 2014; Cowie et al. 2017). The redshift distribution is shown in Fig. 4. The median redshift is 1.146 and ranges from 0.11 to 3.44.

The redshift distribution shows an abundance of sources around ∼2 which could be an imprint from the previously identified z = 1.99 proto-cluster in GOODS-N (Blain et al. 2004; Chapman et al. 2009). This proto-cluster has approximate redshift bounds of 1.982 < z < 2.010, and two VLBI sources (J123618+625541 and J123621+621708) are associated this proto-cluster (Casey 2016). It is expected that this structure is extended spatially beyond the limited field-of-view of the GOODS-N survey as the deep spectra does not extend much outside of the HST coverage. It has been suggested that J123642+621331 (z = 2.018) could also member of this proto-cluster (Murphy et al. 2017), and it is also possible that J123721+621130 (z = 2.02) could be associated, however an accurate spectroscopic redshift would need to be acquired. This raises the possibility that, with deep VLBI surveys and improved number densities of sources, over-densities of VLBI-detected AGN could act as a tracer of proto-clusters for which there is evidence of enhanced AGN activity in multiple wavebands including the radio (e.g. Wylezalek et al. 2013; Krishnan et al. 2017).

5.2. Astrometry

In order to check the astrometry of these observations, we compared these VLBI observations to the MERLIN-VLA observations from Muxlow et al. (2005, hereafter M05), and the reprocessed wide-field images of the same data by Wrigley (2016). The positions were not compared to the global VLBI observations of Chi et al. (2013) due to the known positional uncertainty on the phase calibrator used in those observations (see Appendix A), plus there is a larger number of concordant sources between these EVN observations and M05. Due to computational limitations at the time, M05 targeted 92 radio sources within a 10″ × 10″ square field with VLA-only flux densities >40 μJy in the HDF-N field to a central rms of 3.3 μJy beam⁻¹.

These data were re-processed by Wrigley (2016) and a primary beam corrected 18″×16″ image was used for the subsequent analyses. We used BLOBCAT (with a detection threshold of 6σ) to generate a catalogue of 155 MERLIN-VLA detected sources (with a restoring beam of 0.2″). These were then cross matched with the EVN positions to within 1 arcsec. A total of 25/31 sources were matched to an MERLIN-VLA detected source with the remaining unmatched sources beyond the sky area considered by Wrigley (2016). We estimate a conservative 5 mas error for the astrometry of these new observations, due to uncertainties on the position of the secondary phase calibrator J1234+619, and a 10 mas error on the VLA-MERLIN data arising from calibration and source fitting errors.

As Fig. 5 shows, there is a small systematic offset of 5.5 mas in RA and 0.4 mas in Dec Note that all sources apart from wide-angle tail FR-I source, J123726+621129, are within 60 mas of the MERLIN-VLA positions. This source was excluded from the derivation of the systematic astrometric offset and Fig. 5. The systematic offsets can originate from core-jet blending of the radio emission, standard errors associated with source fitting and calibration, and also errors originating from repeated use of the AIPS task UVFIX when peeling bright sources from the VLA-MERLIN data (e.g. see Morgan et al. 2011). The astrometric scatter is expected to be dominated by blending in the MERLIN-VLA data as substructure such as AGN jets will blend with AGN core emission thus causing an offset in position of peak brightness. These errors should have orientations which are randomly distributed, hence the median 5 mas astrometry offset indicates that we are in fairly good agreement with the MERLIN-VLA positions.

5.3. Comparison to other VLBI surveys

As previously stated, the GOODS-N field has been targeted by two previous wide-field VLBI surveys by Garrett et al. (2001) and Chi et al. (2013). It is worth noting that surveys had restricted field-of-views which are encompassed by the field-of-view of these new EVN observations, therefore we would expect that we should be able to detect all previous identified sources. This survey recovers 11/12 of the Chi et al. (2013) detections and all three of the Garrett et al. (2001) detections. The missing source, J123642+621545, illustrates significant radio variability and has an e-MERLIN integrated flux density of only 60 μJy during the period of these observations, whereas Chi et al. (2013) detects an integrated flux density of 343 μJy.

We compared our results with other wide-field VLBI surveys to ensure that our observations and detection fractions are consistent. In order to calculate the detection fraction, we used the VLA-A array observations provided by Morrison et al. (2010) and cross-matched these with to our VLBI survey. For this analysis, we only consider the central area where we have contiguous imaging and where our primary beam model is most reliable. Within this region, our EVN data were used to image the locations of the 314 VLA sources. Of these sources, 94 were theoretically detectable assuming that they are unresolved on VLBI angular scales with a flux density greater than our 6σ VLBI detection threshold (based upon the VLBI rms noise distribution, as shown in Fig. 2). Of this sample 24 sources were detected with VLBI, thus giving a detection threshold of (24/94). Errors were determined using the Bayesian binomial estimator of Cameron (2011). This is consistent with previous wide-field VLBI surveys. For example, the Chandra Deep Field-South (55 μJy beam⁻¹ rms, Middelberg et al. 2011) has a detection fraction of , the Lockman Hole/XMM (24 μJy beam⁻¹ rms, Middelberg et al. 2013) detects 30±3% and the COSMOS survey (10 μJy beam⁻¹ rms, Herrera Ruiz et al. 2017) detects 20 ± 1%. The mJIVE survey has a detection fraction of 20 ± 1% to a an rms of 60 μJy beam⁻¹ (Deller & Middelberg 2014).

Note that there are some caveats because our VLBI sample is surface brightness limited due to the resolution and flux sensitivity of our EVN observations. For the fainter sources in our targeted sample (e.g. S ≲100 μJy), we can only detect those VLA sources with relatively large VLBI-VLA flux density ratios. For example, as the median flux density ratio is ∼0.6 for this sample (see Sect. 5.4.1), most VLA sources with integrated flux densities ≲100 μJy will go undetected. This motivates deeper VLBI observations in the future as the improved rms noise levels will recover these sources.

5.4. Radio properties of the VLBI selected population

5.4.1. VLBI-VLA flux densities

We compared the 1.5 GHz VLA flux densities to our 1.6 GHz EVN observations in order to investigate the origin of the radio emission. To do this we can use the VLBI to VLA flux density ratio (R) to establish whether a source is dominated by milliarcsecond-scale emission from AGN cores or arc-second scale emission from AGN jets/lobes or star-formation related processes. Note that for this comparison, we only use the VLBI sources which have been primary beam corrected (24/31). As both the VLBI and VLA observations have complete 24 h uv coverage, we use the integrated flux density measurements to define the VLA-VLBI flux density ratio (R≡ S_{i, VLBI}/S_{i, VLA}).

We find that the 18/24 (66.7%) VLBI sources have over 50% of their radio emission originating from a milli-arcsecond scale component. 2/24 (8%) has R > 1 which is most likely due to AGN variability because the observation times of VLA (2011) and VLBI (2014) data vary by a few years. Excluding the two known variable sources with R > 1, we find that the median VLBI-VLA ratio of our observations is 0.625. This is largely in agreement with the COSMOS VLBI survey which find a median VLBI-VLA ratio of 0.6 (Herrera Ruiz et al. 2017). We note that, at low flux densities, our VLBI observations are expected to preferentially detect core-dominated systems, with the majority of VLA arcsecond-scale emission confined to a high brightness temperature core which is detectable by VLBI observations. This is consistent as only a small fraction (∼10%) of our VLBI-detected sources exhibit large-scale radio jets or lobes, while the remaining sources are compact.

An evolution towards more core-dominated systems as we approach μJy flux densities has been hinted at in other VLBI surveys, most notably the mJIVE-20 and COSMOS VLBI surveys (Deller & Middelberg 2014; Herrera Ruiz et al. 2017). There is some evidence that suggests that this evolution may be true. It has been shown that a population of radio sources with core fraction of ∼0.3 below a 1.4 GHz luminosity of 10²⁵ WHz⁻¹ are required by empirical simulations in order to accurately extrapolate the established populations from low-frequency (<5 GHz) surveys to the >10 GHz source populations (Whittam et al. 2017). This could be equivalent to the postulated population of “FR0” sources in the local universe which have core dominated, compact radio emission extending to at most just 3 kpc (Baldi et al. 2015). These wide-field VLBI surveys could be beginning to detect the high-z analogues to this population of radio sources.

5.4.2. Luminosities

The radio luminosity of an object can be used to determine the probability that the radio emission of a source is dominated by AGN activity or star-formation. The monochromatic 1.5 GHz radio power can be calculated using:

(4)

where α is the source spectral index, z is the redshift, d_L is the luminosity distance and S_obs is the measured source flux density (ideally, the zero spacing flux density). Integrated flux densities for all 31 sources were derived using the VLA-A array data outlined in Sect. 2.3. Spectral indices for 24/31 sources were derived using the 5.5 GHz integrated flux densities from the VLA 5.5 GHz survey of Guidetti et al. (2017). For the seven remaining sources, we used the median spectral index of −0.56 from the sources with 5.5 GHz detections, but we note in passing that these redshift based k-corrections only contribute a small factor to the resulting luminosities. As Fig. 6 shows, we sample a large range of radio powers from ∼10²²−10²⁶ WHz⁻¹ which have a median luminosity of 2.7 × 10²⁴ WHz⁻¹.

As expected, in the low redshift range (0 < z < 1), we detect lower luminosity AGN of the order 10²²−10²⁴ WHz⁻¹ which is typical of objects such as Seyfert galaxies. Higher luminosity sources are not detected at low redshift due to the combination of a low density of high power sources plus a smaller cosmic volume surveyed due to the restricted field-of-view. At higher redshifts, these observations preferentially detect higher luminosity objects of the order 10²⁴−10²⁶ WHz⁻¹ which is a consequence of the sensitivity (illustrated by the VLBI sensitivity limit plotted in Fig. 6). This corresponds to radio-loud AGN systems such as FR-I, radio galaxies and quasars (Owen & Ledlow 1994; Smolčić et al. 2017a). Only three of these sources are lobe-dominated (J123644+621133, J123726+621128 and J123636+615659) with extended morphologies in the lower resolution VLA data while the remaining objects are core dominated, with any jets unresolved or of low luminosity undetectable by the VLA.

We used the selection criteria of Magliocchetti et al. (2018) to illustrate the radio populations that this survey is probing. Their criteria defines the crossover point, P_cross, as where AGN related emission is dominant over star-formation related emission in a radio-selected population. At z≤ 1.8, their selection is based upon the radio luminosity functions of McAlpine et al. (2013). In this regime, P_cross approximately scales with redshift as 10^{log10(P_{0, cross})+z} where P_{0, cross} = 10^21.7 WHz⁻¹ sr⁻¹ which is the crossover point from the local universe. Above z = 1.8, the radio luminosity function for star forming galaxies drops off rapidly and P_cross is kept constant at 10^23.5 WHz⁻¹ sr⁻¹. This selection criteria was found to keep contaminants from star-forming galaxies to below 10% at z > 1.8 (see appendix of Magliocchetti et al. 2018).

As Fig. 6 shows, these VLBI observations clearly probe the AGN dominated regime with all radio luminosities at least 3 × P_cross. For VLBI surveys to truly detect statistically significant samples of hybrid systems with both AGN and star-formation related emission, and provide valuable information regarding feedback, either ultra-wide surveys should be used to probe the tail of star-forming galaxies within the AGN dominated luminosity regime (e.g. Herrera Ruiz et al. (2016) investigated radio-quiet AGN using the 2 square degree VLBA survey of COSMOS), or ultra-deep surveys should to used to probe the star-formation dominated luminosity regime which is potentially achievable using the EVN and SKA-VLBI. Nevertheless, we undoubtedly have uncovered hybrid AGN-starburst systems, as we will show in Paper II.

5.4.3. Brightness temperatures

Brightness temperatures were calculated for the VLBI detected objects which were primary beam corrected (24/31). For these measurements, we used PYBDSF (Mohan & Rafferty 2015) to fit an elliptical Gaussian model to each source. We adopted a different selection criterion compared to our detection methodology. In this criterion, referred to as the size measurement detection threshold, sources would be excluded if their S/N ratio were less than ten. This S/N cut-off was chosen because, when using randomly distributed model sources injected onto a noise field from our VLBI data, PYBDSF does not detect all of the injected sources when their S/N < 10. In addition, below this threshold, the variance of fitted sizes is often larger than 20% of the full-width half-maximum of the psf.

As we are concerned with only compact emission, we attempted to exclude emission from radio jets or compact star formation by fitting to the uniformly weighted images (with a restoring beam ∼5.3 × 4.5 mas) if possible. If the S/N of the uniformly weighted image was less than 10, then the naturally weighted images (with restoring beam ∼16 × 16 mas) would be used to calculate the brightness temperatures. Sources with S/N < 10 in both weighting regimes would be excluded completely. Using this selection criteria, 18 sources were selected, 9 using uniform weighting and 9 with natural weighting.

For the calculation, we assume the brightness temperature distribution of a source at redshift z is can be modelled as an elliptical Gaussian radio emission region with major axis θ_maj, minor axis, θ_min, and position angle of the major axis φ. Here in our case, θ_maj, θ_min and φ are the deconvolved major and minor axes and their position angle of the deconvolved elliptical Gaussian model. The brightness temperature, T_b can then approximated using,

(5)

where S ν is the observed flux density and ν is the observing frequency (Condon et al. 1982; Ulvestad et al. 2005). In addition, we took into account the resolution limit for both θ_maj and θ_min using the prescription described in Lobanov (2005). The following expression for the minimum resolvable size, θ_lim,ψ, along each axes of the fitted Gaussian can be calculated using the following equation,

(6)

where b_ψ is the measured FWHM of the psf along the direction of the axis for which the resolution limit is being calculated. S/N is the signal to noise of the image and β is a constant that takes into account the weighting of the visibilities ( β = 0 for uniform weighting and β = 2 for natural weighting). If θ_maj or θ_min were lower than the corresponding minimum resolvable size along each axes, then θ_lim,ψ was used instead to calculate the source frame brightness temperature. A source is classified as unresolved if either axes was below the minimum resolution limit. In this case, the size of the radio emitting region cannot be constrained and only lower limits to the brightness temperature can be derived.

Figure 7 shows the brightness temperature distribution of our VLBI selected sources. Lower limits are derived for those sources which are unresolved or are unresolved in one axis. Sources detected with a 5.3 × 4.5 mas restoring beam have brightness temperatures of the order 10⁷−10⁹ K and, as expected, sources detected with only a 16 × 16 mas restoring beam have lower brightness temperatures in the range (10⁵−10⁶ K). In both weighting regimes, it is highly unlikely that the radio emission is caused by star-formation related processes as typical star-burst galaxies exhibit brightness temperatures of <10⁵ K (Condon et al. 1991).

Note that the majority of the uniformly weighted sources (8/9) are unresolved, hence emission in these objects come from a compact component. In the naturally weighted sources, all of these are classified as resolved. This is most likely caused by the improved surface brightness sensitivity as a result of the increased weighting of the sensitive, but crucially shorter, central European baselines (especially the Effelsberg to Lovell baseline). As the uniformly weighted images have approximately 1.5 × the noise of the naturally weighted images, sources detected in natural weighting will most likely have a compact component with higher brightness temperatures, however this component is currently below the size measurement detection threshold (10σ) with uniform weighting. Indeed 6/9 sources do have a compact component in the uniformly weighted images which is above the source detection threshold of 6σ. The remaining three sources do not have a compact component above 6σ indicating that some flux may be resolved between the two weighting schemes.

6. Conclusions

We present a catalogue of 31 faint VLBI detected sources in GOODS-N to ∼9 μJy beam⁻¹ 1σ rms noise levels and radio luminosities of the order 10²² WHz⁻¹. This has substantially increased the number of VLBI detected sources over previous GOODS-N surveys by Chi et al. (2013) and Garrett et al. (2001), providing a valuable addition to the understanding of the AGN content in this well studied field. We also present a primary beam correction scheme developed for the European VLBI Network. This is publicly available and will be constantly updated.

Additional ancillary information about the radio properties of these objects were derived using VLA 1.5 and 5 GHz data. The detected objects have radio luminosities of the order 10²²−10²⁶ WHz⁻¹ and brightness temperatures in the range 10⁶−10⁹ K. The high fraction of compact emission (as defined by the VLBI/VLA flux density ratio) may be hinting at the existence of radio-loud but core-dominated systems at high redshift which may be equivalent to the population of FR0 objects in the local universe (Baldi et al. 2015, 2018).

With previous VLBI surveys effectively characterising the >50 μJy radio sky, the scheduled 48 h of remaining observations will enable us to reach limiting rms sensitivities of ∼2 μJy beam⁻¹, thus providing the next step in analysing the faint radio-selected AGN population in GOODS-N. The final data release will be able to test whether the increasing fraction of star-forming galaxies at these low flux densities are influenced by AGN activity as postulated in lower resolution surveys. This will be presented in a future publication. In addition, the last 24 h of these observations will include the first ever wide-field VLBI observations using an integrated EVN-eMERLIN array, providing valuable surface brightness sensitivity improvements with the addition of intermediate length baselines.

---

Acknowledgments

We thank the anonymous referee for their comments which helped improve this manuscript. The authors gratefully acknowledge Len Cowie and Amy Barger for allowing us to use unpublished data and for their useful comments. We would also like to thank P. Padovani and R. Windhorst for their useful discussions. This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018). The research leading to these results has received funding from the European Commission Seventh Framework Programme (FP/2007-2013) under grant agreement No 283393 (RadioNet3). The European VLBI Network is a joint facility of European, Chinese, South African, and other radio astronomy institutes funded by their national research councils. e-MERLIN is a National Facility operated by the University of Manchester at Jodrell Bank Observatory on behalf of STFC. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

References

Appendix A Astrometry of J1234+619

We acquired 12 h of e-MERLIN time at 5 GHz in order to independently confirm the position of J1234+619. In order to effectively characterise changes in the delay path induced by tropospheric variations, we implemented a multi-source phase referencing strategy using three VLBA calibrators surrounding J1234+619. These sources have astrometric accuracies <2.6 mas. Whilst not as sophisticated as astrometric techniques such as MultiView calibration (Rioja et al. 2017), multi-source phase referencing has been used to obtain astrometric precisions of ∼20 μas with the VLBA (Fomalont & Kopeikin 2002). We observed three calibrators and then the target with a total cycle time, t_cyc, of ∼10 min (summarised in Fig. A.1 and Table A.1). The use of multiple phase calibrators in different relative positions allows the spatial dependence of the phase to be resolved and a short cycle time is desired in order to map the atmosphere-induced phase variations across all sources more frequently.

In order to estimate the uncertainties in the astrometric measurements, we consider the approach of Rioja et al. (2017). Errors arise from three main sources: (1) Thermal noise which originates from the uncertainty in the position measurement due to random noise. This is characterised as σ_t ∼ θ_B/2(S/N) where θ_B is the restoring beam and S/N is the signal-to-noise. The target field has a S/N of 182, hence the thermal noise error is approximately 0.11 mas (2) Fractional flux recovery. This is the difference between the peak brightness in the phase referenced maps and the self-calibrated maps. (3) Accuracy. The observed phase reference sources have astrometric errors <2.6 mas with 2/3having errors <0.5 mas. Positional errors associated with effects such as the core-shift effect, the use of phase delays over geodetic group delays (Porcas 2009) and differential structure blending effects between the EVN VLBI observations at 1.6 GHz and these observations at 5 GHz is expected to be a few milliarcseconds. In total, errors on these observations is expected to be dominated by unknown effects such as core-jet blending. We therefore assign a conservative estimate of ≤5 mas astrometric accuracy.

With the phase referencing from the three phase calibrators complete, these solutions were transferred to J1234+619 and positions on the resultant image were measured using AIPS task JMFIT. The results of these observations are presented in Table A.2. We find that the eMERLIN 5 GHz positions of J1234+619 are more consistent with those derived as part of these new EVN observations, rather than than the position of Chi et al. (2013). However, due to the unknown errors, such as core-jet blending and structure blending associated with the differing resolutions, we will continue to use the positions derived in these new VLBI observations, as they are more consistent with the 5 GHz eMERLIN positions.

All Tables

All Figures