A GLOBAL 86 GHZ VLBI SURVEY OF COMPACT RADIO SOURCES (2024)

Very long baseline interferometry (VLBI) at millimeter wavelengths offers the best tool for imaging compact radio structures on scales of several dozens of microarcseconds. The first detection of single-baseline interference fringes in an 89 GHz (3.4 mm) VLBI observation was reported by Readhead et al. (1983), demonstrating the feasibility of 3 mm VLBI. Since then, many VLBI observations at 86 GHz have been made, probing the most compact regions in active galactic nuclei (AGNs). However, the number of objects detected and imaged at 86 GHz remained small, compared with the number of objects imaged with VLBI at lower frequencies.

Sensitive VLBI observations at 86 GHz have been made for several sources, including 3C 111 (Doeleman & Claussen 1997), 3C 454.3 (Krichbaum et al. 1995, 1999; Pagels et al. 2004), NRAO 150 (Agudo et al. 2007), NRAO 530 (Bower et al. 1997), M87 (Krichbaum et al. 2006), 3C 273 and 3C 279 (Attridge 2001). In order to increase the number of objects imaged at 86 GHz, four detection and imaging surveys were conducted during the 1990s, with a total of 124 extragalactic radio sources observed (see Beasley et al. 1997; Lonsdale et al. 1998; Rantakyro et al. 1998; Lobanov et al. 2000). In these surveys, fringes were detected of 44 objects, but only 24 radio sources have been successfully imaged. Table 1 gives an overview of these surveys. The low detection and imaging rates of the previous 86 GHz surveys were caused by the relatively poor baseline sensitivities, small numbers of telescopes, and short observing times.

The results of the survey of a larger number of sources can be used to investigate the innermost region of compact jets and to observationally test inner jet models (Marscher 1995): accelerating and decelerating jet models. In the accelerating jet model, the jet accelerates hydrodynamically from the base of the jet, and as the internal energy of the jet plasma is converted into the kinetic energy of bulk flow, the jet Lorentz factor increases along the jet. In this model, an ultra-relativistic neutral beam is generated from the central engine and then the neutrons decay into protons and electrons which form a relativistically flowing plasma. In the decelerating jet model, the central engine produces a highly collimated beam of ultra-relativistic electron–positron pair plasma that scatters photons produced outside the jet (particle cascade). The scattered photons emit X-rays and γ-rays, which decelerates the beam and so decreases the Lorentz factor along the jet.

The theoretical prediction from the jet models leads to the fact that the intensity profiles along the jet are different from each other and have a distinctive shape in each of these models. The resulting brightness temperature can be used to probe the difference. The intrinsic brightness temperature of a distribution of observed brightness temperatures can be determined from a statistical modeling (Lobanov et al. 2000). By estimating the brightness temperature at several frequencies (e.g., 15 GHz, 43 GHz, and 86 GHz) and determining the intrinsic brightness temperatures, we would be able to constrain the physical conditions (e.g., dynamics and compositions) of the innermost region of the compact sources. Moreover, the dependence of the intrinsic brightness temperatures on the observing frequencies will tell about the feasibility of VLBI at higher frequencies (e.g., 150 GHz, 215 GHz, etc).

A large global 86 GHz VLBI survey of compact radio sources was carried out from 2001 October to 2002 October using the Coordinated Millimeter VLBI Array (CMVA) (Rogers et al. 1995), which is succeeded by the Global Millimeter VLBI Array (GMVA).⁴ The main aim of this VLBI survey is to increase the total number of objects accessible for future 3 mm VLBI imaging by factors of 3–5, and to provide the database for the subsequent statistical modeling in order to test the inner jet models.

2.1.Source Selection

2.2.Observational Strategy

In this section, we describe the post-correlation processing of the 3 mm VLBI survey datasets. Fringes were searched in two steps using the Haystack Observatory Postprocessing System (HOPS) and the NRAO Astronomical Image Processing System (AIPS). In the first step, the HOPS task fourfit was used to precisely determine phase residuals. The first fourfit was run with a wide search window (e.g., a width of 1 μs for single-band delay, 2 μs for multiband delay, and 500 ps s⁻¹ for delay rate) centered at zero in delay. Since the fourfit produces baseline-based fringe solutions, the mean and standard deviation of the detected fringe solutions on each baseline were estimated and served as the offset and width of the search window for the second fourfit. The detected fringe solutions from the second run were used to interpolate the offset of the single-band delay for non-detected scans. In the final run of fourfit, an interpolated search window was used with a width of 0.02 μs for the single-band delay and centered at the offset interpolated for each non-detected scan. After this first step of the fringe search, the total number of fringe detections for the survey data was improved by up to 20%. In the second step of the fringe search, the baseline-based fringe solutions were imported into the AIPS using a modified AIPS task MK4IN (Alef & Graham 2002). We then made an antenna-based fringe fit to the data using the AIPS task FRING. Pico Veleta (PV) was chosen as the reference antenna for most of the data. When PV was not available, Fort Davis (Fd) was selected as an alternative reference antenna. The antenna-based fringe fitting was done with a solution interval of 7 min in order to achieve a higher signal-to-noise ratio (S/N). Fringe solutions for strong sources were used to define coarse search windows for the fringe solutions for nearby weaker sources. With the fringe fit of fourfit and FRING, 121 out of the 127 observed sources have yielded fringe detections with S/N ⩾6. Figure 3 shows the S/N distribution of the fringe detection in the entire survey data. Only six sources (0710+439, 1458+718 (3C 309.1), 1749+701 (4C 09.57), 2021+614, 2030+407 (Cyg X-3), and 2031+405 (MWC 349)) are not detected. The highest S/Ns of 425 is measured on the "Pico Veleta–Plateau de Bure" baseline for 1741-038 and the "Effelsberg–Pico Veleta" baseline for 1633+382.

The fringe-fitted data were amplitude calibrated using regular measurements of the system temperatures and antenna gains and the weather information for each station made during the observations. Where possible, time-dependent factors in the antenna power gains were accounted for by applying atmospheric opacity corrections. The AIPS task APCAL was used to calibrate the amplitudes. As a check of the accuracy and consistency of the amplitude calibration, we have investigated (independently for each of the detected sources) the calibrated visibility amplitudes using the best-fit Gaussian component models obtained from the data (the corresponding models are given in Table 7). For each of the sources, the antenna gains were allowed to be scaled by a constant factor so as to optimize the fit by the Gaussian model.

The obtained corrections are within 20% for most of sources, which is also reflected in the average correction factors listed in Table 5. On average, the gain factors for Pico Veleta, Onsala, North Liberty, Owens Valley, and Los Alamos remained within 10% in all three observing sessions. The average gains for Effelsberg did not change much except for the session C. Fort Davis and Mauna Kea required average corrections by more than 20%. Time-dependent errors may still be present in the calibrated data. Therefore we an expect overall calibration accuracy of ∼20–30%.

From the phase- and amplitude-calibrated data, the images were made using the Caltech DIFMAP software (Shepherd et al. 1994). After averaging in frequency, the uv-data were averaged in time over 30 s and were edited for deviant data points.

The uv-data were then fitted with a simple Gaussian model. First, a single circular Gaussian component was applied to fit the data. In case that the single-component model did not represent the data satisfactorily, a multiple-component model was applied. Self-calibration and CLEAN deconvolution were applied to produce final images of the detected sources.

The noise in the final image can be expressed quantitatively by the quality ξ_r of the residual noise. Suppose that a residual image has an rms σ_r and the maximum absolute flux density |s_r|. For Gaussian noise with a zero mean, the expectation of s_r is

where N_pix is the total number of pixels in the image. The quality of the residual noise is given by

When the residual noise approaches Gaussian noise, ξ_r → 1. If ξ_r > 1, not all the structure has been adequately recovered; if ξ_r < 1, the image model has an excessively large number of degrees of freedom (Lobanov et al. 2006). The values of ξ_r of the images in the survey are presented in Column 14 of Table 6 and the distribution of them is shown in Figure 4, implying that the images adequately represent the structure detected in the visibility data.

The initial CLEAN cycles for the region around the core component are conducted with using natural weighting and without using uv-tapering. Once the CLEAN models satisfactorily fit the visibility at the longest baselines, uv-tapering was applied to the data on long baselines in order to recover faint emission further out from the core component. We did not modify the visibility amplitudes, except for introducing an overall, time-constant gain correction factor wherever it was required for improving the agreement between the CLEAN model and the data. In addition to the check for the antenna gain corrections described above, we investigated the changes of the visibility amplitudes with and without introducing the time-constant gain correction factor, by using the correlated flux densities S_S, S_L, at the shortest and longest baselines B_S, B_L, (presented in Columns 4–7 in Table 6), for each of the sources. The correlated flux densities S_S,L obtained after introducing the antenna gain corrections are compared with the flux densities S'_S and S'_L before the gain corrections. As shown in Figure 5 for the distributions of the ratios R_S = S_S/S'_S and R_L = S_L/S'_L, the visibility amplitudes on the shortest and longest baselines for each of the sources were not changed for most of the sources during the hybrid imaging. For a small number of peculiar sources, the ratios of the visibility amplitudes fall within a range of 0.75–1.25. This analysis shows again that the amplitude calibration error of this survey observations is 20%–30%.

We introduce a zero-baseline flux to recover a faint structure in the extended region by adding a fake visibility at the point of origin in the Fourier plane (uv-plane). Since the shortest baseline of this survey observations is about 50–100 Mλ, this may result in the faint structures of the extended regions appearing to be negative because the flux in every pixel of the map is offset by a small negative amount. The effect may be countered to some extent by adding the fake visibility at the origin of the uv-plane. The measured total flux density of each source S₈₆ (listed in Column 3 of Table 6) is used as the zero-baseline flux.

In order to extract quantitative information from the images, circular Gaussian-component models were used to fit the self-calibrated uv-data, yielding the following parameters: the total and peak flux densities, positions, and sizes of each component. The uncertainties of the models were estimated, based on the S/N of detection of a given model fit component, using an analytical (first order) approximation (Fomalont 1999). The general fit parameters of a component in VLBI images of radio sources are S_tot (total flux density), S_peak (peak flux density), σ_rms (post-fit rms), d (size), r (radial distance, for jet components), θ (position angle, measured for jet components, with respect to the location of the core component). The uncertainties of the fit parameters can be estimated by adopting the approximations given by Fomalont 1999:

where σ_peak, σ_tot, σ_d, σ_r, and σ_θ are the uncertainties of the total flux density, peak flux density, post-fit rms, size, and radial distance of a component, respectively. When the size, d, of a component was determined, the resolution limits (Lobanov 2005) should be taken into account. So, the minimum resolvable size of a component in an image is given by

where a and b are the axes of the restoring beam, S/N is the signal-to-noise ratio, and β is the weighting function, which is 0 for natural weighting or 2 for uniform weighting. When d < d_min, the uncertainties should be estimated with d = d_min.

We use the results of the model fitting to estimate the brightness temperatures of the core and jet components. The rest frame brightness temperature T_b of the emission region represented by a Gaussian component is

where λ is the wavelength of observation, z is the redshift, and k is the Boltzmann constant. Practically, the brightness temperature can be calculated by simplifying (6):

where the total flux density S_tot is measured in Jy, the size of the circular Gaussian component d in mas, and the observing frequency ν in GHz. If d < d_min, then the lower limit of T_b is obtained with d = d_min.

Out of 127 sources, 109 sources have been imaged and model fitted. The first 3 mm VLBI maps for 90 sources were made in this survey, increasing the number of sources ever imaged with 3 mm VLBI observations up to 110.

In Figure 6, we present two plots and one contour map for each source at each epoch. In the left panel, the plot of the visibility amplitudes against uv-radius is shown. The corresponding uv-sampling distribution is given in the inset. The X-axis of the plot of the visibility amplitude represents the uv-radius which is the length of the baseline used to obtain the corresponding visibility point. The uv-radius is given in the units of 10⁶λ, where λ is the observing wavelength. The Y-axis of the plot shows the amplitude of each visibility point (i.e., correlated flux density) in units of Jy. The uv-sampling distribution in the inset of the left panel describes the overall distribution of the visibility in the uv-plane, whose maximum scale equals that of the uv-radius. In the right panel, the contour map of each source is shown, with the X- and Y-axes in units of milliarcseconds. For each source, the source name and the observation data are given in the upper left corner of the map. The lowest contour level is identified in the lower right corner of the map. The shaded ellipse represents the FWHM of the restoring beam in the image. In all of the images, the contours have a logarithmic spacing and they are drawn at −1, 1, 1.4, ..., 1.4ⁿ of the lowest flux density level. For 12 sources (0133+476, 0149+218, 0212+735, 0234+285, 0238 − 084, 0316+413, 0355+508, 0415+379, 0430+052, 0716+714, 1928+738, and 2255 − 282), multi-epoch images are presented. Most sources are centered on the brightest component (VLBI core), but for some sources with a larger structure we have shifted the center to fit the image in the box.

In Table 6, the parameters of the images presented in Figure 6 are summarized. For each image, Table 6 lists the source name, the observing epoch, the total flux density, S₈₆, obtained from the pointing and calibration measurements made at Pico Veleta during the observation, the correlated flux densities, S_S,L, measured on the shortest and longest baselines, B_S,L, the parameters of the restoring beam (the size of the major axis, B_a and the minor axis, B_b, and the position angle of the beam, B_PA), the total flux, S_t, the peak flux density, S_p, the off-source rms, σ, and the quality of the residual noise in the image.

Table 7 lists the parameters of each model-fit component: the total flux, S_tot, peak flux density, S_peak, size, d, radius, r (only for jet components), position angle, θ (the location of the jet component with respect to the core component), and measured brightness temperature, T_b. For sources with multiple components, parameters of the core component are followed by those of jet components. For sources observed at multiple epochs, individual epochs are marked. The estimated uncertainties are given next to each parameter. The upper limits of size, d, and the lower limits of brightness temperature, T_b, are in italic with brackets.

6.1.Source Compactness

6.2.Brightness Temperature T_b

6.3.Intraday Variable Sources

We have conducted the largest global 86 GHz VLBI survey of compact radio sources during three GMVA sessions in 2001 October, 2002 April, and 2002 October. Participation of sensitive European telescopes augmented by the VLBA antennas ensured high baseline and image sensitivities. A total 121 out of 127 sources observed have been detected at least on one baseline and 109 sources have been imaged with a typical dynamic range exceeding 50. The survey observations have resulted in an increase by a factor of 5 of the total number of sources imaged at 86 GHz with VLBI.

We have used two-dimensional, circular Gaussian components to fit the observed visibilities and parameterize the source structure. Using the results of these fits, the source compactness and brightness temperatures have been derived.

We find that almost all of the survey objects are resolved and the cores of about 70% of the imaged sources are resolved. Radio galaxies are less compact than quasars and BL Lac objects. BL Lac objects are similar to quasars in the compactness at sub-milliarcsecond scales.

The distribution of the core brightness temperatures peaks at ∼10¹¹ K and only 1% of the cores have brightness temperatures higher than 10¹² K. This shows apparently lower brightness temperatures than those derived from other VLBI surveys at lower frequencies (e.g., 5 GHz and 15 GHz).

IDV sources in our sample are similar to non-IDV sources in compactness at sub-milliarcseconds. The cores of IDV sources are smaller in angular size and so yield a higher brightness temperature than non-IDV sources, since the core flux densities of both samples are similar to each other.