The correlated visibility data from our K-CRF VLBA observations are calibrated with the NRAO’s Astronomical Imaging Processing System (AIPS), using a semi-automated approach. The data calibration, for the most part, follows the VLBA calibration pipeline described in Appendix C of the AIPS Cookbook, using the standard utilities available in AIPS. An automated pipeline is used for self-calibration, imaging, and deconvolution using the Caltech Difference Mapping software (DIFMAP). We generate the image maps, u,v-coverage plots, and scan-averaged calibrated visibility amplitude versus baseline-length plots through custom Python routines developed in-house. 

More detailed information about the calibration and imaging process is available from de Witt et al., 2023.

We obtain the following parameters from each of the final K-band images: (1) the number of scans and number of visibilities, (2) the peak brightness, total clean flux density, and the weighted average correlated flux density over four baseline length ranges, (3) the background rms brightness level over the entire residual image and the image signal-to-noise-ratio (peak brightness/rms), (4) the quality of the fit between the observed and model visibilities after self-calibration, (5) the maximum absolute brightness value in the residual map, (6) the clean beam minor and major axes FWHM and position angle, and (6) an estimate of the residual rms phase calibration error. 

The image parameters and clean component models are then used to determine (1) the flux density variability of each source over time, (2) the source structure quantities, that gives a measure of the amount and extent of the source structure, and its variation over time, and (3) the quality of each image.

For a source to be suitable as a VLBI calibrator or reference source it should be relatively bright at  the frequency of observation and should be detectable with a good signal-to-noise (ideally >> 10) on all baselines and within the coherence time (typically less than 2 minutes at K-band). A calibrator or reference source should also be compact or dominated by an unresolved or only marginally resolved "core" at the frequency of observation, and should not show any significant variation in source structure over time.

The source structure quantities, and the variability of each over time are used to evaluate the suitability of each source as a calibrator or reference source. The source structure quantities include: (1) a measure of the degree of source (angular) compactness or core domination, (2) the radial extent as an indicator of the extent of the source structure and overall angular size of a source, and (3) the structure index that provides an estimate of the astrometric quality of the source by calculating the median value of the structure delay corrections (the additional phase terms due to source structure for each clean component and each VLBI baseline). 

It is important to take into consideration the quality of each image. Poor image quality means that most of the structure quantities are poorly determined, i.e. they will be noisier for images of poor quality, and  they may vary from epoch to epoch, not because of intrinsic variation, but because they are noisier. Since the scheduling of these K-band observations was mainly optimised for astrometric reference frame work, the images are not all of the same quality, with some sources observed more often and with better u,v-coverage than others. In addition to the scheduling constraints, the quality of the observations and images are also affected by the availability of antennas, weather conditions and unexpected outages and errors, which will vary with observing session.

More detailed information about the image parameters and structure metrics is available from de Witt et al., 2023.

The core flux density (Score) is defined as the sum of the flux densities contributed by all Clean components within a 0.25 mas angle originating from the brightest pixel in the image. This specific threshold was chosen to match the upper limit of angular resolution for a global VLBI array employing maximal baseline lengths of 10,000 km at 24 GHz.

The flux-density variability index serves as a quantifier for the extent of flux density variability, defined as the standard deviation divided by the mean of the flux density across all imaging sessions for the particular source. An index of 0 signifies a non-varying source.

Given the potential time-varying nature of source structure, the levels of compactness and core dominance can likewise exhibit variability over time. We consider several measures of source compactness and core dominance: (1) the peak fraction (C1), which is the ratio of the peak brightness to the clean flux density, where an error-free image of a completely unresolved or point-like source would have have a value of 1, (2) the core fraction (C2), which is the ratio of the core flux density to the Clean flux  density  where a value of 1 indicates that all flux density resides in the core, and (3) the correlated flux density fraction (C3), which is the ratio of average correlated flux density over the longest baselines to that over the shortest. Values around 1 indicate compact sources.

The standard deviation over all sessions for which a source was imaged gives a measure of the variability of the source compactness measures.

In addition to the metrics mentioned, we generate time-series plots for other source structure attributes, such as the radial extent and radial extent within an area < 10 times the geometric mean of the FWHM beam from the brightest pixel (R and R10), the flux-weighted radial extent (E and E10), and the structure index (SI). We also construct time-series plots for the image quality factor (Q) of each image. Factors considered during the calculation of the image quality factor encompass: (1) the fraction of calibrated visibilities compared to expected visibilities, (2) the proportion of residual calibration errors in each image, (3) the elongation of the synthesized beam in each image, and (4) a comparison of image resolution to the minimum anticipated resolution.

For a comprehensive exploration of these concepts, we recommend referring to de Witt et al., 2023.

The time-series plots for NRAO140, yield valuable insights. Notably, they highlight a trend in which the flux density experiences an increase over time, coinciding with a transition toward a slightly more core-dominated source configuration (and conversely). The images themselves portray a fading of weak extended emission as the flux density amplifies. 

Intrinsic source structure can have a significant effect on a VLBI bandwidth synthesis delay measurement. The influence of source structure on delay measurements can be assessed through the computation of corrections applied to the bandwidth synthesis delay. These corrections are deduced from the inherent spatial brightness distribution of the source - the source model that forms the foundation of the VLBI image (in this case the image Clean components) -  for a range of baseline lengths and orientations projected onto the plane of the sky, quantified in terms of wavelength units (u,v-coordinates).

To carry out these calculations, we employ the code developed by and described in Shabala et al. (2015), based on the work published in Fey & Charlot 1997. This code facilitates the computation of an additional phase term due to source structure for each image Clean component considered, and for each VLBI baseline (less than the diameter of  the Earth), in a 512 × 512 grid in the u,v-plane. This is repeated for each intermediate frequency (IF) band. The resultant average slope of phase against frequency then yields the group delay, arising from the source's structural attributes, effectively serving as a correction for structural delay.  

The structural index, SI, is then introduced as a parameter. It is defined as 1 + 2 × log (median group delay), where the (median group delay) is just the median value of the structure delay corrections, quantified in units of picoseconds (ps). Within this framework, an SI value between 0 and 2 indicates a source characterized by compact structure or faint extended emission. A value closer to 3 signifies a source marked by substantial structural features, while a value of 4 or more implies a source boasting pronounced extended emission and/or intricate structural elements.

For a comprehensive understanding of these methodologies, we recommend referring to Fey & Charlot 1997, and Shabala et al. (2015).

We employed a straightforward two-component model to fit the calibrated visibilities of each source image, utilising the MODELFIT task within DIFMAP. Circular Gaussian models were implemented to establish initial estimates of key source characteristics. This approach particularly aided in approximating the flux density and FWHM angular size for both the primary or "core" component (c1) and the second-brightest component (j1). Additionally, it assisted in assessing the offset and orientation of the second brightest component (j1) or in distinguishing potentially blended components within the "core". Given the extensive volume of data, we developed an automated pipeline dedicated to handling the modelfitting process efficiently.

In tandem with the two-component modelfitting in DIFMAP, we derived estimates for the principal angle of extension, commonly referred to as the 'jet direction'. This was accomplished by fitting a line through the locations of the image Clean components. This custom Python routine, developed in-house, provided a less resource-intensive yet effective means of gauging source elongation and validating the modelfitting robustness in DIFMAP. Both an unweighted and flux-density-weighted fit were executed using this method, originating from (0,0) and spanning the clean component locations. Angles obtained from this approach exhibited reasonably close agreement (within a few degrees) with those extracted from DIFMAP modelfitting, even for sources seemingly core-dominated in their images.

For an in-depth exploration of these methodologies, we strongly recommend referring to de Witt et al., 2023.

Furthermore, we generate time-series plots (not displayed here) that illustrate the flux density, FWHM angular size of the primary (c1) and secondary (j1) component, as well as the offset and direction of the secondary component (j1) for each source across all epochs. A time-series plot spanning all 33 epochs of NRAO140 is accessible at this location.