The Celestial Reference Frame at K Band: Imaging. I. The First 28 Epochs

We present K-band (24 GHz) images of 731 compact extragalactic radio sources with submilliarcsecond resolution, based on radio interferometric observations made with the Very Long Baseline Array of 10 telescopes during 29 day long sessions spanning from 2015 to 2018 and recorded at 2048 Mbps. Many of these sources are imaged with submilliarcsecond resolution for the first time at frequencies above X band (8 GHz). From each of the K-band images, we derive the following source properties: peak brightness, core and total flux density, the ratio of peak and core to total flux (compactness measure), radial source extent, structure index, source size, and jet direction. The vast majority of sources are imaged at multiple epochs, providing insights into their temporal behavior. The use of K band was motivated by the fact that the sources are generally intrinsically more compact at higher frequencies, as well as by the factor of 3 improvement in interferometer resolution relative to the historically standard S/X band (2.3/8.4 GHz) used for a large amount of reference frame and calibrator work. Lastly, as most of the sources imaged here are in the K-band component of the third International Celestial Reference Frame, these images serve to characterize the objects used in that International Astronomical Union standard.


Introduction
For millennia, the positions of objects in the sky have been measured against what were previously called "fixed" stars. However, by the end of the 1980s, it became clear that the socalled "fixed" stars had significant motions that were corrupting the reference frames based on those stars. In response, the International Astronomical Union (IAU) resolved (IAU resolution C1, 1988 6 ) to start using active galactic nuclei (AGNs) at cosmological distances (redshifts of the order of z = 1) in order to access sources with negligible parallaxes and proper motions, where any errors in the position due to parallax and proper motion were much smaller than one part per billion (1 nanoradian or 206 μas). The first realization of this new strategy was the International Celestial Reference Frame (ICRF; Ma et al. 1998), hereafter ICRF-1, followed by its successor, the ICRF-2 (Ma et al. 2009), both of which are based on dual-frequency S/X-band (2.3/8.4 GHz) observations that have been conducted since the end of the 1970s. The third and most recent realization, the ICRF-3 (Charlot et al. 2020) included data up until 2018 from the original S/X frequencies as well as observations at K band (24 GHz) and at X/Ka band (8.4/32 GHz), extending the ICRF up from 8 GHz for the first time. Various efforts are now underway to continue the improvement of the ICRF, thereby answering the IAU's call for maintenance of the ICRF-3 (IAU Resolution B2, 2018 7 ). This paper documents our efforts to improve and better understand the K-band component of the ICRF-3 and its successors.
Why were K-band observations added to the ICRF-3? K band was added for several reasons: (1) to allow data from multiple independent data sets to be compared in order to better assess accuracy, (2) to circumvent weaknesses inherent in S/X observations, and (3) to realize potential advantages of higher frequencies. The ideal extragalactic radio source would be compact and core-dominated on very long baseline interferometry (VLBI) scales (nanoradians), be relatively bright at the frequency of observation, have no or little detectable proper motion, and should have accurately known positions, stable at the submilliarcsecond level. AGNs, most of which are at cosmological distances, do not exhibit any measurable proper motion or parallax, making them ideal Celestial Reference Frame (CRF) sources.
Unfortunately, at the standard S/X frequencies, many ICRF sources exhibit spatially extended intrinsic structure that may vary with time and frequency (e.g., Beasley et al. 2002;Gordon et al. 2016;Hunt et al. 2021). The degree to which such variations are resolvable varies with baseline length and observing frequency. Such structure can introduce significant errors in the VLBI measurements thereby contributing to the uncertainty in the measured astrometric source positions (Charlot 1990). Extended source structures, if not corrected for, create an astrometric noise floor limiting the ability to further improve the accuracy of VLBI source positions. Although it is possible in principle to correct for the effects of such structure using high-resolution images of the sources, the experience we gained in the course of the present work suggests that it is not realistic to do so for the thousands of sources in the ICRF at the near monthly cadence needed at the submilliarcsecond accuracy of the ICRF.
Furthermore, it is becoming more difficult to obtain reliable S/X-band VLBI data. The increasing radio frequency interference, especially at S band (e.g., WiFi), is significantly degrading the ability to collect clean S/X data. In addition, spacecraft tracking sites that have historically contributed to VLBI observations are phasing out S band due to the low telemetry rates relative to X and Ka bands.
The advantages of K band over the historically standard S/Xband frequencies are: (1) generally more compact structure (e.g., Charlot et al. 2010), (2) ∼3× higher resolution, and (3) reduced core shift effect-the frequency-dependent shift of the VLBI peak brightness or core position (e.g., Bietenholz et al. 2004;Sokolovsky et al. 2011). Frequency-dependent opacity, in particular synchrotron self-absorption, is the main reason why the apparent source structure is a function of frequency. Typically the core components (or base) of the jets are optically thick, and have flat or inverted radio spectra, while the jets are optically thin and have normal or steep radio spectra. In particular, VLBI images of AGNs show that the observed position of the peak brightness point moves closer to the central black hole as the frequency increases (e.g., Sokolovsky et al. 2011). Thus, astrometric VLBI observations at higher radio frequencies should allow for more accurate and well-defined CRFs, including improved ties of the VLBI frame to optical frames such as Gaia (Gaia Collaboration et al. 2016. For example, Liu et al. (2021) compared the ICRF-3ʼs multifrequency radio positions and the Gaia Early Data Release 3 and concluded that the K band to Gaia position offset is statistically smaller than the S/X-band to Gaia offset for sources that show extended structures. For all of these reasons, we have been motivated to build a K-band CRF.
K-band astrometric and imaging observations commenced in 2002 with Lanyi et al. (2010) and Charlot et al. (2010), resulting in a K-band catalog of 268 sources, but with weak coverage from 0°to −45°decl. and no coverage below −45°, as well as several large regions north of −45°with no sources. Starting from this foundation, in 2014, a new K-band collaboration was formed in order to produce full sky coverage, improve spatial density, and improve astrometric accuracy (e.g., de Witt et al. 2019). As of early 2022, the K-band CRF consists of 1035 relatively uniformly distributed sources derived from more than 1.6 million observations collected in 122 observing sessions (D. Gordon et al. 2023, in preparation). These observations comprise sessions from the Very Long Baseline Array (VLBA) in the United States and singlebaseline sessions between the Hartebeesthoek Radio Astronomy Observatory 26 m antenna in South Africa and the Hobart 26 m antenna in Tasmania, Australia.
Because extragalactic sources are known to vary in both flux density and structure on timescales of months to years, and because there is also the potential for sudden flaring events, it is important to directly image CRF sources. While the CRF sources are in general expected to appear more compact at K band than at X band, they can still exhibit measurable extended emission at K band. It is therefore desirable that these sources be monitored on a regular basis, even at higher radio frequencies, through periodic VLBI observations. Accordingly, we have used our VLBA observations to create sensitive, highresolution, multiepoch images of 731 sources at K band from which we derive quantitative estimates of source strength, size, compactness, and the temporal stability of these quantities. These images aid in estimating the accuracy of the K-band CRF as well as characterizing the sources for optimal use in various applications. For many sources, images at frequencies above X band are published for the first time, thus enabling the first direct comparisons of source structures at X versus K band for many of the sources common to both bands. In summary, since 2014, significant progress has been made toward an improved CRF at K band and a roadmap for further improvements has been defined .
Our work is motivated by the benefits of high-accuracy CRFs for many applications such as determining the Earthʼs orientation in space (e.g., Eubanks 1993;Krásná et al. 2019), studying the motion of tectonic plates (e.g., Carter & Robertson 1993), the alignment of radio and optical reference frames (e.g., Bourda et al. 2010;Gaia Collaboration et al. 2018;Liu et al. 2021), studies of Galactic aberration (e.g., MacMillan et al. 2019), satellite tracking, orbit determination, deep-space navigation (e.g., Thornton & Border 2000), alignment of the planetary ephemerides (e.g., Park et al. 2021), and as phase reference calibrators for VLBI imaging of weak and extended sources and measurements of parallaxes and proper motions (e.g., Reid et al. 2019). The ICRF also contributes toward the Global Geodetic Reference Frame, 8 which is the subject of a recently adopted United Nations resolution on a global geodetic reference frame for sustainable development. The move to higher radio frequencies is also essential to the National Aeronautics and Space Administration (NASA), the European Space Agency (ESA), and the Japan Aerospace Exploration Agency (JAXA), where deep-space missions are using Kaband (32 GHz) CRFs (e.g., Jacobs et al. 2019) to navigate closer to the Sun while also achieving higher telemetry data rates. In addition to the improved source compactness and enabling observations closer to the Sun, Kband also allows observations closer to the Galactic Plane. Differential VLBI astrometry on water masers, which improve our understanding of the structure of the Milky Way galaxy (e.g., the BESSEL project, Reid et al. 2019), is only possible at K band where the maser water line resides, thus requiring well-characterized K-band reference sources.
The discussion of the imaging results obtained from our VLBA multiepoch K-band CRF observations is organized as follows: In Sections 2 and 3 the details of the observations and data reduction are given. The main results from the imaging are given in Section 4. We use an analysis technique similar to that of Charlot et al. (2010) to characterize the source structure and evaluate the suitability of these sources for high-frequency CRF work. The image analysis techniques and results are given in Section 5. In addition to the imaging, estimates of the source characteristics were obtained from model fitting, and the model-fitting approach and results are discussed in Section 6. Short notes on sources with significant structure are provided in Section 7.

Observations
The imaging observations were taken with the VLBA, which consists of ten 25 m Cassegrain antennas (Napier 1995), using the 1.2 cm feed. The VLBA projects BJ083 and UD001 provided the series of observing sessions used for the images presented in this paper. The ten VLBA antennas and their station codes are listed in Table 1, together with typical system equivalent flux density (SEFD) values at 24 GHz for the elevation angles >40°achieved for the BJ083 and UD001 sessions.
Sessions were nominally 24 hr in duration with scans of individual sources of 90-120 s duration. Most sources were observed for 1-3 scans per session. The BJ083 series consisted of five sessions between 2015 July and 2016 June at a center frequency of 24.568 GHz for sessions BJ083A and B, and 23.968 GHz for BJ083A1, C and D. The UD001 series consists of 24 sessions (UD001A to X), observed between 2017 August and 2018 July, at a central frequency of 23.568 GHz. We recorded right circular polarization (IEEE convention) with a total data rate of 2048 Mbps using the ROACH Digital Back End in the Polyphase Filter Bank (PFB) mode and Mark-5C recorders (Whitney et al. 2010). The data were sampled with a 2-bit resolution yielding a total sky bandwidth of 512 MHz, correlated with 32 spectral points in each of 16 intermediate frequency (IF) bands. IF16 was corrupted due to a combination of edge effects and the PFB algorithm that was used, and thus was edited out. From the EVN Sensitivity Calculator, 9 using a bandwidth of 480 MHz and a 90 s scan time, we calculate a theoretical baseline sensitivity of 3 mJy (1σ), allowing for a 10σ detection of a source with 30 mJy of correlated flux density, and a 1σ image sensitivity of 0.5 mJy beam −1 . Table 2 gives a detailed list of the sessions, observing dates, the number of sources observed and imaged in each session, and notes on significant outages.
The set of 1035 sources in the current K-band CRF (D. Gordon et al. 2023, in preparation) was taken from a number of inputs. First we took the 268 sources of the KQ project (Lanyi et al. 2010), which were generally chosen to be strong, flatspectrum sources to ensure detection at frequencies higher than X band. Then we added sources from the K-band Galactic survey (Petrov et al. 2011) to include the region near the Galactic Plane, and sources from the X/Ka celestial frame (Jacobs et al. 2012) to maximize overlap with the X/Ka celestial frame to increase the usefulness for spacecraft navigation. We also included additional flat-spectrum sources from the ICRF S/X-band sample to increase the overlap with the S/X celestial frame and maximize the number of ICRF-3 defining sources in the K-band frame. Of the 1035 sources, 731 were observed with the VLBA as part of the BJ083 and UD001 sessions. J1745-2900, or Sgr Aå is not an extragalactic source, and although it was observed as part of the BJ083 and UD001 series and included in the imaging analysis, it is not part of the K-band CRF astrometric solution. Results of the position and proper motion of J1745−2900 derived from the K-band VLBA astrometric observations will be detailed in a forthcoming paper by Gordon et al. (2023). The 732 sources observed in the BJ083 and UD001 sessions are listed in Table 3, together with their J2000 coordinates, optical identifications, magnitudes, redshifts, and the number of BJ083 and UD001 sessions that each source was observed in.  Note. Session UD001P was considered a failed session and none of the images were considered for analysis. Note that from 2017 September 17 to 2018 March 3, SC was out of the array due to hurricane damage. The loss of SC as the eastern most station during this period degraded the u,v-coverage significantly. Table columns: (1) session name, (2) session start date, (3) number of sources observed in each session, (4) number of successful images produced from each session, and (5) comments on outages and flagged antennas.
The experiment schedules were generated using the NRAO scheduling software SCHED. 10 All of the sessions were scheduled dynamically to allow for flexibility in the start date and time, mainly to avoid bad weather at many of the stations. The schedules were generated using a sequence that minimizes the slewing time under several constraints to optimize the u,vcoverage for absolute astrometry and imaging. The hour angle based optimization mode OPTMODE = HAS, first used by the VLBA Imaging and Polarization Survey (Helmboldt et al. 2007), was used to schedule source scans well separated in hour angle over time. A good distribution of source hour angle is essential to estimate good source positions for astrometry, to build up good u,v-coverage for imaging, and to obtain a good spread of elevations to help separate the troposphere from other parameters. Since separating the effect of vertical station locations from that of the troposphere depends mostly on obtaining observations at the extremes of elevation, all sessions were scheduled with a 5°minimum elevation.
The source list for each schedule was chosen in such a way as to provide a good geometric distribution of sources in both R.A. and decl. with time. A good spread in R.A. and decl. is needed to optimize the number of sources and scans observed in each session, while a wide range of decl. is essential to cover the widest range in elevations to get a well-calibrated troposphere. The number of requested scans for each session was set lower for sources farther south to help match the optimization to what is practical. Sources below a decl. of −30°, were mostly scheduled with one scan only, and sources between −30°and −15°with two scans. Sources at higher decl. were scheduled with three scans per session to build up good coverage in the u,v-plane. The final number of scans per source, limited by the number of requested scans (one to three scans as described above), is determined by the dynamic schedule and optimization constraints at the time of observation. In some of the earlier sessions, a few sources were mistakenly scheduled twice per session, which led to five sources that were observed in four scans, four sources that were observed in five scans, and three sources that were observed in six scans.
The BJ083 and earlier UD001 sessions used a scan length of 120 s, which permits around 170 unique sources per session. From UD001I onwards, the scan length was decreased to 90 s, which increased the number of unique sources observed in each session to around 250. As a result, the source list was divided into four sets that were rotated with each new session. A core set of 40 sources was scheduled in each session as a unifying tie among all sessions to keep the astrometric solution stable, and to use as fringe finders for correlation and phase calibrators for imaging. This core set was made up from ICRF-2 defining sources that are also in the Lanyi et al. (2010) catalog, and are highlighted with an asterisk ( * ) next to the source name in Table 3.

Data Reduction
The data were correlated with the VLBA DiFX (Deller et al. 2011) software correlator at the Array Operations Center in Socorro, New Mexico, using 32 spectral points (each 1 MHz wide) and an integration time of 1 s. The correlated visibility data were calibrated with the NRAO's Astronomical Imaging Processing System (AIPS; Greisen 2003), using a semiautomated approach. The data calibration, for the most part, followed the VLBA calibration pipeline described in Appendix C of the AIPS Cookbook, 11 using the standard utilities available in AIPS. Data inspection and initial editing were done in the standard manner using information from observing logs and correlator reports. This was followed by amplitude gain corrections due to sampler voltage offsets and calibration of the instrumental  Note. This table is published in its entirety in machine-readable format. A portion is shown here for guidance regarding its form and content. Table columns: (1), (2) J2000 and IVS source names. The core set of 40 sources are listed with an asterisk ( * ) next to the J2000 source name, (3) D means ICRF-3 defining source, (4), (5) J2000 coordinates from the ICRF-3 K-band astrometric catalog, rounded to the nearest mas. For state-of-the-art positions, see Charlot et al. (2020). The coordinates for Sgr Aå are taken from Gordon et al. (2023), (6) single-band delays using a manual phase-calibration approach. The task APCAL in AIPS was used to perform an initial amplitude calibration for each IF, using the measurements of the antenna gain and T sys collected during the observations, as well as correcting for atmospheric opacity using the available weather data. Before the standard global fringe-fitting was performed, an additional nonstandard calibration step for removing the tropospheric and clock errors was added (described in Mioduszewski 2009). Multiband or group delays were obtained from fringe-fitting in AIPS. The task DELZN was then used to estimate the zenith atmospheric delay and clock error for each antenna from the multiband delays by fitting a polynomial function. The task DELZN can take a maximum of 30 sources as input, which were selected from the core set of 40 sources. After removal of the tropospheric and clock errors, the standard global fringe-fitting was done using a point-source model with a solution interval equal to the scan length.
The amplitude calibration for each session was further improved by self-calibrating ∼22 of the brightest and most compact sources (which have a strong signal on all baselines), spread out over the observing run. This step serves to improve the initial images for each of the sources, and thus the reliability of the per-source self-calibration done subsequently (described below). A CLEANed image of each of these ∼22 sources was used as a model to determine the gain correction factors needed to make the visibility data of the source/model conform as closely as possible to the CLEAN component model. The derived gain corrections were combined into one solution table and examined. Amplitude versus u,v-distance plots from the self-calibrated data for the ∼22 selected sources were also examined. Antenna gain corrections were applied only if a consistent trend or clear pattern was visible across different sources, e.g., the gain or amplitude being consistently lower or higher for one antenna. Twelve of the sessions observed between BJ083A and UD001K showed consistently low amplitudes for PT over the entire observing run, and large manual amplitude gain corrections of about 30% on average, were needed. Session UD001O had three antennas, BR, HN, and NL, that showed consistently low amplitudes over the entire observing run that required manual gain corrections of more than 50%. The data quality of UD001P was deemed to be so poor that it was discarded.
After inspection and editing of obvious bad data points, the data were separated into single-source files and written to UVFITS file format. All subsequent processing of the singlesource visibility files was carried out using the Caltech Difference Mapping software (DIFMAP; Shepherd 1997). An automated pipeline was used for self-calibration, imaging (Fourier inversion), and deconvolution (using the standard CLEAN algorithm of Högbom 1974) in DIFMAP. A point-source model was used as a starting model for all sources. The data were initially phase self-calibrated and CLEANed using uniform weighting, and then after several iterations of phase selfcalibration, we switched to natural weighting, giving more weight to the shorter baselines that may reveal additional extended emission present in the field. The iterative process of self-calibration and CLEANing was taken to have converged when the peak in the residual image became less than six times the rms brightness level of the residual image from the previous iteration. The data were then amplitude self-calibrated using first a single, time-independent, amplitude gain correction factor for each antenna, followed by a time-dependent amplitude selfcalibration using a solution interval of 60 minutes. These additional amplitude gain corrections were typically <6%, and amplitude gain corrections for each antenna averaged over all sources and sessions are listed in Table 4. The gain corrections for BR, NL, and PT were typically larger than those for other stations, with BR requiring gain corrections of >20% for session UD001A.
After the imaging and self-calibration, estimates of the source characteristics were obtained through fitting models directly to the visibility data by least squares. A simple twocomponent model was fitted to the calibrated visibilities for each source using the MODELFIT task in DIFMAP. In addition to the model fitting in DIFMAP, a line was fit through the locations of the image CLEAN components to estimate the source characteristics over the total extent of the source structure and demonstrate the robustness of the model fitting in DIFMAP. The model-fit process and the results from the model fitting are detailed in Section 6.

Imaging Results
Images for 731 of the 732 sources were produced from the phase and amplitude self-calibrated data using DIFMAP. J2207 −0002 (IVS B2205−002, ICRF J220755.2−000215) was not detected in any of the sessions. A total of 5879 images from a possible 6225 were produced from the BJ083 and UD001 sessions, giving a failure rate of only 5.6%, with the failures mainly a result of poor u,v-coverage for some of the observations and the loss of session UD001P. There are sources such as J0432+4138 (3C 119, ICRF J043236.5 +413828) that appear to have a very complex structure that was too difficult to fully recover using the automated DIFMAP pipeline, for which better images could probably be obtained using a more careful interactive imaging approach. However, due to the sheer volume of this data set, no attempt was made to redo any of the images by hand. Figure 1 shows the sky distribution of the 731 sources imaged in the BJ083 and UD001 sessions, color-coded by the number of sessions from which images were produced for each source. From the total of 731 sources that were imaged, 26 were imaged at a single epoch only, while 705 sources were imaged at two or more epochs, with 63 sources imaged at 2 epochs, 147 sources imaged at 3-5 epochs, 303 imaged at 6-9 epochs, and 192 sources imaged at 10-28 epochs. For the majority of the images where the resolution was close to the full resolution of the array, inspection showed that most of the emission is within ∼4 mas of the brightest pixel. Weak apparent emission at a distance of more than ∼10 synthesized beams (which corresponds to ∼6 mas considering the average value of ∼0.61 mas measured for the geometric mean of the full width beam major and minor axis over all of the images) is in many cases, not real, with for example, emission showing up at different locations in different epochs, or disappearing from one epoch to the next. Such artifacts are likely the result of imaging errors due to residual calibration errors or poor u,v-coverage. Figure 2 shows the images, together with the corresponding u,v-coverage and amplitude versus u,v-distance plots for a representative sample of four sources. The amplitude versus u, v-distance plots show the scan-averaged calibrated visibilities. The two plots were generated using auxiliary tools from the VLBI processing processing software PIMA 12 (Petrov et al. 2011).
The four sources in Figure 2 were chosen to represent a wide range of decls. (δ = − 40°to δ = + 70°), which leads to a wide range of u,v-coverages. In particular, sources north of δ ≈ 60°, such as J1459+7140 (3C 309.1, ICRF J145907.5+714019), will be above the horizon for all telescopes in the array for all times (circumpolar) and thus have the potential for excellent u, v-coverage from the rotation of the baselines through the 24 hr session. On the other hand, given that the VLBA is an all northern array, sources near the equator, such as J0433+0521 (3C 120, ICRF J043311.0+052115), and sources farther south, such as J0106−4034 (IVS B0104−408, ICRF J010645.1 −403419), will be visible for a rather limited range of hour angles and thus the u,v-coverage, in general, will be poorer. A source at mid-northern decl., J0005+3820 (IVS B0003+380, ICRF J000557.1+382015), was also included to give a wide sampling of the expected quality of u,v-coverages. The complete figure set (731 figure sets containing 5879 images) is available. Table 5 lists the relevant parameters and quality measures obtained from the final K-band images for each of the sources: the source name, the observing session, the number of participating antennas (N), the number of scans (M), the peak brightness (S p ), the total CLEAN flux density (S cln ), and the median correlated flux density over short (S s ) and long (S l ) baselines, the background rms brightness level over the entire residual image (s im ), the image signal-to-noise ratio (S/N im ) defined as s S p im , the quality of the fit between the observed and model visibilities after self-calibration (σ), the maximum absolute brightness value in the residual map (S r ), the CLEAN beam minor and major axes FWHM and position angle (q min , θ maj , bpa), and the estimated residual rms phase calibration error (ò).
The median correlated flux density was computed (using the Auxiliary tools from the PIMA software) over two ranges of projected baseline length: <70 Mλ (<900 km), and >400 Mλ (>5000 km) with the latter corresponding to the flux density on angular scales of <0.52 mas. Of the 5879 K-band images, two had no data on projected baselines <70 Mλ, and 586 had no data on baselines >400 Mλ. Figure 3 shows the distribution of the source peak brightness, the total CLEAN flux density, and the median correlated flux densities, taken from the most recent epoch that each of the 731 sources was observed at. The mean (median) values are 0.43 (0.24) Jy/beam −1 for the peak brightness and 0.56 (0.31) Jy for the total CLEAN flux density. The mean (median) of the correlated flux densities are 0.53 (0.29) Jy over the shortest baselines, and 0.35 (0.20) Jy over the longest baselines.
A similar approach to that of Charlot et al. (2010) was used to analyze the fractional variation of the peak brightness and the total CLEAN flux density from one session to the next. The average peak brightness and the average total CLEAN flux density were computed over each of the sessions (excluding the sessions BJ083A and UD001P, both of which had poor data quality), for a sample of 25 sources that were observed in all of    Notes. This table is published in its entirety in machine-readable format. A portion is shown here for guidance regarding its form and content. Table columns: (1), (2) source name and session that it was observed in, (3), (4) number of participating antennas and number of scans, (5) number of visibilities, (6) image peak flux density, (7) total flux density defined by the sum of all image CLEAN components, (8) correlated flux density measured on baselines <900 km (70 Mλ), (9) correlated flux density measured on baselines >5000 km (400 Mλ), (10) image background rms brightness level, (11) image S/N, (12) rms of the squared differences between the observed and model visibilities (after self-calibration) divided by the variances implied by the visibility weights, (13) maximum absolute brightness value in the residual map, (14), (15) minor and major axes FWHM of the restoring beam, (16) position angle of the beam measured from north going toward east, and (17) estimated upper limit of the residual rms phase calibration error. a The number of visibilities is counted individually for each IF and then added together to give N vis . b This is equivalent to the square root of the reduced χ 2 , except that no account is taken of the number of degrees of freedom implied by the number of gains being varied.
(This table is available in its entirety in machine-readable form.) the sessions. Taking the ratio of the standard deviation and the average over this distribution gives an estimate of the fractional variation of the peak brightness and the total CLEAN flux density of 9.66% and 9.82%, respectively, across all of the sessions. Good diagnostic measures of the image quality are the background rms brightness level, s im , hereafter called image rms, and the image signal-to-noise ratio, S/N im , which as mentioned earlier is defined as s S p im . Table 5 lists the values of s im and S/N im for each of the images, and the median and standard deviation computed over each session are listed in Table 6. For comparison, the median and standard deviation of the theoretical thermal noise (σ th ) and the expected K-band signal-to-noise ratio (S/N th ), computed over each session, are also listed. We define S/N th as the image peak flux density divided by the expected theoretical thermal noise (S p /σ th ). The expected value of σ th , as calculated from the EVN calculator 13 is 0.224 mJy beam −1 (1σ), assuming a total integration time of 360 s (120 s × 3 scans per source), a total bandwidth of 512 MHz (taken over all 16 IFs), and assuming a nominal zenith SEFD for all 10 antennas. The expected value of σ th for each source was scaled from this value to account for the number of visibility records remaining after calibration, and actual bandwidth (accounting for the discarded IF16).
The mean (median) values of s im and σ th over all of the Kband images are 1.62 (0.70) and 0.40 (0.38) mJy beam −1 , respectively. The mean (median) values of S/N im and S/N th over all of the K-band images are 375:1 (374:1) and 1528:1 (695:1), respectively. While the image rms values for some of the weaker sources agree well with the expected theoretical thermal noise, the mean image rms is ∼4 times higher, and the median image rms is ∼1.8 times higher, than the theoretical thermal noise. Similarly, the mean S/N im is ∼4 times smaller, and the median S/N im is ∼1.8 times smaller, than S/N th . For brighter sources, systematic errors such as errors in the residual gain phase or amplitude are the dominant factors limiting the S/N, rather than the thermal noise, and thus the image rms values in brighter sources are expected to be larger than the thermal noise. On the other hand, sources with a very complex structure may show a much larger image rms because of the sometimes limited u,v-coverage and the automated imaging approach. The maximum deviation from the expected theoretical thermal noise is seen for J2253+1608 (3C 454.3, ICRF J225357.7+160853), a ∼14 Jy source with a complex structure. Perley (1999) derived the relationship between the image dynamic range and the residual calibration errors, in the absence of other errors. In particular for errors in the gain phase only, with rms f in radians, they found: where D is the image dynamic range, and M, N are the number of scans and number of antennas, respectively, where it is assumed the calibration is per scan. Perley (1999) also showed that the same value D is obtained with errors in gain amplitude only when the rms of the fractional amplitude gain error ò is substituted for f in Equation (1). In the presence of both phase and amplitude gain errors, the image dynamic range would be lower than that given by Equation (1). D as calculated from Equation (1) therefore gives an upper limit on the dynamic range achievable in the presence of gain phase errors with rms f. If any other errors are present, the image dynamic range would be lower than D.
We invert Equation (1) to estimate an upper limit on the residual gain phase errors (f), from the observed image dynamic range, S N im . More precisely, f represents the rms errors in gain phase (in radians) in the case that there were no other errors or noise affecting the image. The median and standard deviation of f, averaged over each session, are listed in Table 6, and show that the upper limit on the residual error is typically less than ∼0.08 radians, which is equivalent to ∼5°.

Image Analysis
For a source to be suitable as a VLBI calibrator or reference source (e.g., for astronomical phase referencing or spacecraft navigation), it should be relatively bright at the frequency of observation and should be detectable with a good S/N (ideally ?10) on all baselines and within the coherence time (typically less than 2 minutes at K band). A calibrator or reference source BJ083A Note. Table columns: (1) session name, (2) minimum and maximum number of participating antennas for each session, (3) median (standard deviation) of the image rms for each session, (4) median (standard deviation) of the image signal-to-noise ratio (S/Nim) for each session, (5) median (standard deviation) of the expected theoretical rms noise for each session, (6) median (standard deviation) of the expected signal-to-noise ratio (S/N th ) for each session, and (7) median (standard deviation) of the estimated upper limit of the residual phase errors.
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. For most calibrator uses, amplitude variability within the observing time (24 hr) is bad, while amplitude variability over longer timescales is only a problem if a source becomes so faint as to fall below S/N thresholds. The image parameters and CLEAN component models of the 731 K-band sources were used to determine the suitability of each as a calibrator or reference source. The quality of each image as well as the flux density variability of each source over time, and the amount of source structure and its variation over time, were determined using several different measures or quantities. These measures are described in detail in the following subsections.

Image Quality
In order to make informed decisions about the suitability of each source as a calibrator, and reliably determine any variations in time of this suitability, it is important to take into consideration the quality of each of the images. Since the scheduling of these K-band observations was mainly optimized 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 is also affected by the availability of antennas, weather conditions, and unexpected outages and errors, which will vary with observing session. In order to evaluate and compare the quality of the images, the following parameters were considered: 1. The visibility fraction, V f , measures the ratio of the number of visibility measurements to the total expected number of visibility measurements and is defined as: where N vis is the number of visibilities from the final calibrated visibility data of each image (listed in Table 5), and N vis_std is the expected number of visibilities for a source observed in a standard way, i.e., N vis_std =194,400 considering the standard three scans of 90 s each, and including all 45 baselines of the 10 VLBA antennas and all 16 IFs. The visibility fraction, V f , is mostly <1. However, for a few sources that were observed in more than three scans, or with scans >90 s in length, V f can be >1. The actual values of V f range from 0.05-1.49, with a mean (median) value of 0.50 (0.47). 2. The resolution fraction, R f , provides a measure of the image resolution and is the ratio of the minimum calculated full interferometric resolution to the geometric mean of the synthesized beam FWHM minor and major axes. The resolution fraction is defined as: where q min and θ maj are the minor and major axes of the synthesized beam for each image, and θ calc is the angular scale of 0.30 mas corresponding to the longest VLBA baseline of 8611 km at a frequency of 23.5 GHz. The actual values of R f range from 0.10-0.86, with a mean (median) value of 0.54 (0.57).
3. The beam fraction, B f , provides a measure of the elongation of the synthesized beam and is defined as: where q min and θ maj are the minor and major axes of the synthesized beam. The values of B f can range from 0-1, where 1 is for a circular beam. The actual values of B f range from 0.06 to near unity, with a mean (median) value of 0.46 (0.41). 4. The error fraction, E f , is based on Equation (1), which provides an estimate of the magnitude of the residual calibration errors in each image, provided that calibration errors are in fact the dominant source of image errors. The error fraction is defined as: im f E f can range from <0 to 1, with 1 representing an image with no residual calibration errors. The actual values of E f for our images range from 0.16 to near unity, with a mean (median) value of 0.96 (0.97).
Finally, we take the overall quality of each image to be the average over all of the quality parameters. The image quality factor is then defined as: The image quality factor, Q, ranges from <0 for the worst possible images to ∼1 for our best images, with the actual values ranging from 0.33-0.97. Figure 4 shows the distribution of Q taken over all 5879 images, with the mean and median value of this distribution both being 0.61. Table 7 lists the mean and median value of Q for each of the BJ083 and UD001 sessions. The degradation in the image quality seen for sessions UD001I to UD001S corresponds to the period when the SC antenna was out of the array due to hurricane damage. The loss of SC, the easternmost station, significantly degrades the u,v- coverage and resolution of the array, and therefore also the image quality. The image quality factor for each image is listed in Table 8.

Flux-density Variability
The flux-density variability index is a measure of the degree of variability of the flux density. It is defined as the standard deviation divided by the mean of the flux density over all sessions that the source was imaged. A variability index of 0 would indicate a source that does not vary. Thus, the variability index of S cln is defined ass S S cln cln , where s S cln is the standard deviation andS cln is the mean value measured over all sessions in which a source was imaged. The variability index of S s and S l for each source is defined ass S S s s , ands S S l l , respectively. Similarly, the variability index of S p is defined ass S S p p . However, S p is manifestly dependent on the beam for resolved sources, and so the S p variability may only be telling us that the beam is variable.
We also define the core flux density, S c , as the sum of the flux densities of all of the CLEAN components within an angle of 0.27 mas of the brightest pixel in the image. The 0.27 mas limit was chosen because that is the maximum angular resolution of a global VLBI array with maximum baseline lengths of 10,000 km at 23.5 GHz, and therefore the size below which the component would be unresolved by any Earth-based 23.5 GHz observations. The S c value for each image is listed in Table 8, and the S c variability index is defined ass S S c c .
The distributions of the variability indices of S p , S c , S cln , S s , and S l for each of the 705 sources that were observed at two or more epochs are shown in Figure 5. The mean (median) values of the variability indices shown in Figure 5 are 0.19 (0.17) for both S p and S c , and 0.17 (0.16) for both S cln and S s , while the mean (median) values for the S l variability indices are slightly higher at 0.26 (0.24). The peak brightness and flux density variability indices for each source are listed in Table 9.

Source Compactness and Variability
As mentioned, (angular) compactness is an important quality of good calibrator sources. Compactness is strongly correlated with the degree to which a source is dominated by an unresolved core. Since the source structure can be timevariable, the degree of compactness and core domination can also be time-variable. There are several measures of source compactness and core dominance. We consider the following: 1. The peak fraction, C 1 , is the ratio of S p to S cln , where an error-free image of a completely unresolved or pointlike source would have C 1 = 1. Due to imperfect deconvolution and the occurrence of negative CLEAN components, C 1 values of >1 are possible. The values of C 1 for each of the 5879 images are listed in Table 8. Figure 6 shows the distribution of C 1 over all 731 sources, taken from the last epoch that each source was observed at. The mean (median) of the distribution in Figure 6 is 0.80 (0.83), which indicates that on average the K-band CRF sources are relatively compact.
The standard deviation over all sessions for which a source was imaged (s C 1 ) gives a measure of the variability of C 1 . We list the per-source mean values of C 1 ,C 1 and the s C 1 in Table 9. The distribution of s C 1 for the 705 sources that were observed at two or more epochs is shown in Figure 6. The mean (median) values of s C 1 are relatively low, at 0.06 (0.05), implying that the variability of C 1 is relatively low. 2. The core fraction, C 2 , is the ratio of S c to S cln , where a value of 1 represents a compact source where all of the flux density is within the core (i.e., within an angular radius of 0.27 mas). As with C 1 , values >1 are possible.
The values of C 2 are listed in Table 8 for each of the 5879 images. Figure 6 shows the distribution of C 2 over all 731 sources, taken from the last epoch that each source was observed at. The computed mean (median) values from the distribution in Figure 6 are 0.84 (0.87), which indicates that on average the sources are relatively compact, with 87% of the flux density in the core. The standard deviation over all sessions for which a source was imaged (s C 2 ) gives a measure of the variability of C 2 . We list the per-source mean values of C 2 ,C 2 and the s C 2 in Table 9. The distribution of s C 2 for the 705 sources observed at more than one epoch is shown in Figure 6. The mean (median) of s C 2 is 0.06 (0.05), indicating that on average the variability of the flux density in the core compared to the total CLEAN flux density is low. 3. The correlated flux density fraction, C 3 , is yet another measure of how resolved the source is, and is defined as the ratio of S l to S s . For perfect images, C 3 could range from 0-1; however, in the presence of calibration errors or very poor u,v-sampling, values slightly >1 are  Table 8. Figure 6 shows the distribution of C 3 over all of the sources, taken from the last epoch that each source was observed at. The mean (median) values of C 3 are 0.71 (0.73), which indicates that on average the source flux density remains relatively constant on all baselines. The standard deviation over all sessions for which a source was imaged (s C 3 ) gives a measure of the variability of C 3 . We list the per-source mean values of C 3 ,C 3 and the s C 3 in Table 9. The distribution of s C 3 for the 705 sources observed at two or more epochs is shown in Figure 6. The mean (median) of the s C 3 is 0.11 (0.10).

Radial Extent and Variability
The radial extent, R, as defined by Ojha et al. (2004), is an indicator of the extent of the source structure and is a measure of the overall angular size of a source. R is defined as the radius around the peak-brightness point that contains 95% of the total CLEAN flux density. Recalling that, for our sources, most of the flux density is within ∼4 mas from the brightest pixel, we also define a radial extent, R 4 , which is the radius that includes 95% of the flux density that is within 4 mas of the brightest pixel. This serves to exclude low-brightness extended emission beyond 4 mas, which has less effect on the suitability for astrometry. The values of R and R 4 are listed in Table 8. Figure 7 shows the distribution of both R and R 4 , taken from the last epoch that each source was observed at. The mean (median) values are 1.38 (0.64) mas for R and 0.79 (0.56) mas for R 4 , which indicates that, for our sources, most (95%) of the emission is generally within ∼1 mas of the brightest pixel in the image.
The standard deviation over all sessions for which a source was imaged (σ R and s R 4 ) gives a measure of the variability of R and R 4 . We list the per-source mean values of R and R 4 (R and R 4 ) and the σ R and s R 4 in Table 9. The distributions of σ R and s R 4 are shown in Figure 7. The mean (median) values are 1.03 (0.31) mas for σ R and 0.31 (0.23) mas for s R 4 . The high mean value of σ R (compared to s R 4 ) is probably an artifact of the distant components being the most poorly recovered in the CLEAN images, and the possibility that the most distant components, for at least some of the poorer-quality images, are spurious as well. This supports our reasoning for defining R 4 in the first place.

Flux-density-weighted Radial Extent and Variability
The flux-density-weighted radial extent, E, as defined by Ojha et al. (2004), provides another measure of the source compactness by considering the extent of the source structure weighted by flux density. E is defined as: where r i is the radius from the brightest pixel in the image at which the ith CLEAN component lies, and S i is its flux density.
The flux-density-weighted radial extent, considering only CLEAN components within 4 mas from the brightest pixel, E 4 , was also computed. The values of E and E 4 are listed in Table 8. Figure 8 shows the distributions of E and E 4 , taken from the last epoch that each source was observed at. The mean (median) values are 0.25 (0.17) mas for E, and somewhat smaller at 0.16 (0.12) mas for E 4 . The median value of 120 μas measured for E 4 is of a similar magnitude to the astrometric agreement between the ICRF-3 S/X and K-band frames, which has a weighted rms value of ∼111 μas for the decl. differences 14 .
The standard deviation over all sessions for which a source was imaged (σ E and s E 4 ) gives a measure of the variability of E and E 4 . We list the per-source mean values of E and E 4 (Ē and E 4 ) and the σ E and s E 4 in Table 9. The distributions of σ E and Table 8 Results from the Image Quality and Source Structure Analysis, for All of the Sources Observed in the BJ083 and UD001 Sessions

Compactness
Radial Extent Note. This table is published in its entirety in machine-readable format. A portion is shown here for guidance regarding its form and content. Table columns: (1) source name, (2) session name, (3) image quality factor, (4) core flux density, (5)-(7) source compactness ratios; the peak fraction, the core fraction, and the correlated flux fraction, (8)-(11) the radial extent, the radial extent within 4 mas, the flux-density weighted radial extent, and the flux-density weighted radial extent within 4 mas, (12) the median value of the structure delay corrections (τ median ), and (13) source structure index.
(This table is available in its entirety in machine-readable form.) 14 This is still only a crude measure given that we have not measured the Xband weighted radial extent, which is presumably a bit worse and is not yet accounted for.
s E 4 are shown in Figure 8. The mean (median) values are 0.16 (0.09) mas for σ E and 0.05 (0.04) mas for s E 4 , showing that on average the variation in the flux-density-weighted radial extent, per source, is relatively low.

Structure Index and Variability
The CLEAN component models were also used to calculate the source structure index, SI, defined by Fey & Charlot (1997) and based on the analysis by Charlot (1990), to provide an estimate of the astrometric quality of the sources. Charlot (1990) demonstrated the significant contribution that intrinsic source structure can have on a VLBI bandwidth synthesis delay measurement. The effect of source structure on the delay measurements can be estimated by calculating corrections to the bandwidth synthesis delay based upon the observed source structure for a range of u,v-plane coordinates. We used the code developed by and described in Shabala et al. (2015) to calculate the additional phase term due to source structure for each CLEAN component (combination of brightness ratio and position) considered, 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 K-band IF, and the resultant average slope of phase against frequency then yields the group delay, τ, due to the source structure (structure delay correction). The SI is then defined as: where a value of 1 or 2 indicates a source with a compact structure or weak extended emission, while a value of 3 indicates a source with significant structure, and a value of 4 a source with strong extended emission and/or a complex structure.   Note. This table is published in its entirety in machine-readable format. A portion is shown here for guidance regarding its form and content. Table columns: (1) source name, (2) peak brightness mean (variability index), (3) core flux density mean (variability index), (4) clean flux density mean (variability index), (5) correlated flux densities over short baselines mean (variability index), (6) correlated flux densities over long baselines mean (variability index), (7) peak fraction mean (standard deviation), (8) core fraction mean (standard deviation), (9) correlated flux fraction mean (standard deviation), (10) radial extent mean (standard deviation), (11) radial extent within 4 mas mean (standard deviation), (12) flux-density-weighted radial extent mean (standard deviation), (13) flux-density-weighted radial extent within 4 mas mean (standard deviation), and (14) structure index mean (standard deviation).
(This table is available in its entirety in machine-readable form.) The values of τ median and SI are listed in Table 8. The distributions of τ median and SI, taken from the last epoch that each source was observed at, are shown in Figure 9. The mean (median) values are 4.35 (2.15) ps for τ median , and 1.66 (1.67) for SI, showing that on average, the source structure has only a small effect on delay measurements.
The standard deviation over all sessions for which a source was imaged (σ SI ) gives a measure of the variability of SI. We list the per-source mean values of SI,SI, and the σ SI in Table 9. The mean (median) values of σ SI are 0.55 (0.49), showing that on average, per source, the SI does not vary over more than one index (i.e., 1, 2, 3, or 4). When considering only the discreet values as defined in Equation (9), the difference between the minimum and maximum SI discrete is either 0 or 1 in most cases; with 239 sources with a difference of 0, and 334 sources with a difference of 1. There are 122 sources where the difference is 2, and only 10 sources with a difference of 3.

Summary of Image Analysis
In order to evaluate the quality of a source for use as a calibrator or astrometric reference, the source structure quantities, C 1 , C 2 , C 3 , R, R 4 , E, E 4 , and SI, and the variability of each over time, should all be considered. However, these measures should be used with caution, especially the variability indicators, as the quality of the images for a particular source may differ significantly, leading to spurious apparent variation in the image-based quantities. 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. We do not have an uncertainty for the source structure quantities, and it is therefore advisable to take into consideration the image quality factor, Q, when evaluating the source structure and its variability over time. One should be particularly careful when interpreting the variation in the image-based quantities or comparing the variability between different sources. Apart from poor image quality that may lead to spurious variation, the measured variability is also affected by the sampling size, range, and spread (i.e., the number of observations, the time frame of the observations, and how spread out the observations are over that time frame), and the timescale of the intrinsic variation of each source. Figure 10 shows the time-series plots of the peak brightness, flux densities, and each of the source structure quantities, with the print version having just two example sources, while the plots for all 731 sources are available. The two sources in Figure 10 were chosen to represent both a source with compact structure and one with a complex structure, and to illustrate the effect that poor image quality may have on our evaluation of the suitability of a source as a calibrator or astrometric reference. The first example source in Figure 10, J0423−0120 (IVS B0420−014, ICRF J042315.8−012033), is one with very high variability in flux density over time. The values of S p , S c , S cln , S s , and S l all have a variability index of close to 40% over the 28 epochs of observation. However, based on the source structure quantities, J0423-0120 is a relatively compact source, over all of the epochs, withC 1 ,C 2 , andC 3 values of 0.82, and s C 1 , s C 2 , and s C 3 values of 0.06, and withĒ andĒ 4 values of 0.25 mas and σ E and s E 4 values of 0.08 mas. TheSI over all of the epochs is 1.30, with a σ SI value of 0.34. From the plots in Figure 10, and visual inspection of the images, it can be seen that the core of J0423-0120 fades as the structure becomes more extended (e.g., Figures 11(a) and (b)), but the extended emission is weak and the majority of the flux remains in the core over time. For this particular source, the poorer-quality images, where Q < 0.6 (e.g., Figure 11(c)), do not give significantly different results for the source structure quantities when compared to images of higher quality.
For the second example source in Figure 10, J2230+6946 (IVS B2229+695, ICRF J223036.4+694628), the images (e.g., Figures 12(a) and (b)) reveal a close, equal double with one component fading over time while the other component brightens, with the brighter of the two components always at the center of the map (as a result of self-calibration starting with an initial point-source model). Visual inspection of the 18 images of J2230+6946 shows that only a few of them have a high enough resolution to distinguish between the two components, and most of the images show a single component (e.g., Figure 12(c)), sometimes with a one-sided bulge. The double nature of this source is most evident from the epoch 2-8 images, where the image quality is highest and close to 0.8. The source structure quantities also differ significantly between the different images, suggesting a fairly compact source for those images where the image quality is lowest. In conclusion, we find that in some cases, poor-quality images may suggest a compact source, when in fact it is a double.

Model Fitting
The imaging results presented in the previous sections provide valuable information about the source structure and the suitability of each source as a reference or calibrator source. However, when comparing two images taken at different epochs, it can be difficult to determine whether the apparent changes from one image to the next are due to real changes in the source or due to differences in the sampling of the u,vplane. A more straightforward approach would be to compare the measured visibilities directly, such as comparing the amplitude versus u,v-distance plots between different epochs (as shown in Figure 2). Since the u,v-coverage is different for each epoch, comparing the visibilities directly is not feasible. Instead we fitted models directly to the calibrated visibilities by least squares (without imaging or CLEANing). Doing so allows us to estimate the flux density and FWHM angular size of the core and the second brightest component directly from the visibilities, and is complementary to the imaging results. Estimates of the overall extent and direction of the source structure were also obtained by fitting a line through the locations of the image CLEAN components.
These estimates of the source characteristics can be used to study the effect of source structure on the astrometric and geodetic VLBI observables and to either filter out the sources with the least pointlike structure or to apply structure corrections (due to the added delay) to the delay observables (e.g., Sovers et al. 2002). Both of these techniques have been used to improve astrometric and geodetic VLBI observations, although filtering out sources with significantly resolved structure could come at the cost of losing a significant fraction of the sources in the reference frame. Estimates of the component flux density and FWHM angular size would also provide some indication of whether these sources are bright and compact enough for K-band observations on longer baselines and for selecting possible reference sources for higher radio frequency frames, such as Ka band (32 GHz), Q band (43 GHz), and W band (86 GHz). The model fitting works well for sources with bright second components or extended emission. However, because of the automated model-fitting approach, care should be taken when interpreting the results from noisier observations that suffer from poor u,v-coverage and for sources with a core-dominated structure, in particular those observations where the source structure only appears to be core-dominated due to poor u,vcoverage.

Model-fitting Approach
A simple two-component model to estimate the source characteristics was fit to the calibrated visibilities for each source using the MODELFIT task in DIFMAP. Since we expect most of any extended structure to be resolved out and the sources to be relatively compact and core-dominated at K band, the simplest approach to determine the suitability of each source as a reference or calibrator source would be to fit a single-component circular Gaussian model to the calibrated visibility data to get an estimate of the FWHM angular size for each source. However, the images show that a number of sources have complex structure, not well modeled with a single circular Gaussian. One of the simplest extensions beyond the single-component model is a model consisting of two circular Gaussian components. For this reason, a second circular Gaussian component (i.e., a two-component model) was fit to all of the sources to get an estimate of the offset and direction of the second brightest component or to distinguish between possible blended components in the core.
Initial tests using the BJ083C session showed that the fits converged more reliably for circular than for elliptical Gaussian components. For this reason, circular Gaussian models were used to get a basic idea of the scale and direction of the source structure. Visual inspection of the visibility plots and the brightness distribution of the CLEANed images for session BJ083C showed that the circular Gaussians generally provide adequate fits to the visibilities for all sources. A circular model is a fairly simple approximation and caution should be used when drawing conclusions about the actual source geometry from the results of the fits, since even a model with two circular Gaussians will be only a very approximate representation of the true source structure.
Because of the large volume of data, a two-staged automated pipeline was developed for the model fitting. During the first stage, a single circular Gaussian model was fit to the core (or main component) of each source by minimizing the χ 2 . The free parameters, in this case, are the component flux density and the source FWHM angular size. Due to the self-calibration in phase, the peak brightness point is very near the phase-center position; thus, we do not fit the position of the first component. We take the total CLEAN flux density as the starting value for the component flux density. Because the angular size is not well known, we started the fits with a range of FWHM sizes, and then picked the one resulting in the smallest χ 2 as the best fit.
The second stage again involves the fitting of only a single circular Gaussian, but this time the starting values of the component flux density and FWHM size are obtained from the stage one trial model that gave the smallest reduced χ 2 . We then Fourier-transformed the residual visibilities, i.e., after the subtraction of the best-fit single Gaussian component, to form a residual map. The peak brightness and position of the brightest point in this residual map are then used as starting values for fitting a second component. The first and second component are then fit simultaneously, and the six free parameters are then the flux density and the FWHM angular size of the first and second component and the position of the second component with respect to the phase center. We again started the fits with a range of FWHM sizes for the second component and selected, as the final best fit, the set of fitted values with the lowest χ 2 .

Model-fitting Results
Table 10 lists the parameters obtained from the model fitting for each of the sources and each session in which it was imaged: the reduced χ 2 from the two-component model fitting, the flux density and the FWHM angular size of the main or core  Table 10 only record values for a second component if it satisfies the following criteria: the flux density of the second component is at least ś 5 im (where s im is the image rms), the radius of the second component from phase center is greater than q 2 min (where q min is the minor axis of the synthesized beam), and it is smaller than or equal to 10× , the image quality (bottom left), the radial extent (top right), the flux-density-weighted radial extent (middle right), and the structure index (bottom right). The gray dashed horizontal line in the bottom-left plots, of each source, helps to distinguish between image quality factors above and below 0.6, and in the bottom-right plots, to distinguish between structure indices above and below 2. For each of the time-series plots, the mean value is shown together with the corresponding variability index (for peak brightness and flux density) or standard deviation (for structure quantities) in parentheses. The complete figure set (731 plots) is available. the geometric mean of the FWHM beam size. From the 5879 images produced from the BJ083 and UD001 sessions, 587 had no significant second component from the MODELFIT results, i.e., either the fit did not converge for the second component or the criteria mentioned above were not satisfied.
When considering the estimated FWHM angular size of the components, it is important to note that it is not possible to reliably determine sizes much smaller than the resolution and that the minimum measurable size depends also on the S/N of the observations. An approximate expression for the minimum measurable size in the visibility plane is given by Lee et al. (2008), to be where q _ min size is the minimum measurable size in units of mas, a and b are the minor and major axes of the restoring beam, and β is the weighting function, which is 0 for natural weighting and 2 for uniform weighting. Consequently, if the FWHM angular size of a component is q < _ min size , then only an upper limit equal to q _ min size is listed in Table 10. We used β = 0 for natural weighting, which gives a minimum measurable size of 0.03 mas for an image with = S N 100 im at the maximum resolution of the VLBA at 24 GHz, i.e., a = b = 0.30 mas. Figure 13 shows the distribution of the MODELFIT estimates of the FWHM angular sizes of the core and second component, as well as the distribution of the radius or distance of the second component from the phase center, all measured in units of milliarcseconds. The mean (median) values of the estimated FWHM angular sizes are 0.14 (0.11) mas for the core components and 0.29 (0.17) mas for the second components. On average, the core components are relatively compact compared to the maximum resolution of the VLBA at 24 GHz, which is 0.30 mas. The mean value of the FWHM size for the second components is somewhat larger than the median value, suggesting a tail to the distribution, with a few cases of values much larger than the median. However, there are some sources that have significant extended emission on opposite sides of the core component, and in a few cases, the emission on both sides of the core was modeled as a single second component in DIFMAP (i.e., the second component is large enough to extend on either side of the core), explaining some of the large outliers obtained for the FWHM size. The mean (median) values of the estimated radius or distance from the phase center for the second components are 0.91 (0.62) mas, showing that on average the second components are within 1 mas from the core.

Error Estimation
To estimate realistic uncertainties for the parameters derived from model fitted VLBI data is a difficult task. The statistical uncertainties obtained from the model fitting are calculated under the assumptions that the errors in the visibilities are independent and Gaussian-distributed and that the visibilities are perfectly calibrated, none of which are generally true. The statistical uncertainties, therefore, are almost guaranteed to be too small. Lee et al. (2008) estimated uncertainties based on the S/N of detection of a given model-fit component by adopting the approximations given by Fomalont (1999), but these estimates do not consider the systematic errors, and are also likely too small. More realistic uncertainties, in particular, estimating the contribution due to any residual miscalibration, can be obtained using Monte Carlo simulations (e.g., de Witt et al. 2016 andPetrov et al. 2019 used Monte Carlo simulations randomizing the antenna gain values to obtain uncertainties of the fitted sizes). However, this method does not work well for multiple components (e.g., if the second component is weak, the uncertainty would also depend on the ratio between the flux of the core and the second component and the separation between them) or parameters other than the estimated FWHM size. Since it would be too labor intensive to carefully extract uncertainties for each image and each component, and because the purpose of the model fitting was only to provide some estimate of the size and radius for astrometric and geodetic VLBI applications, we calculate only an average measure of the uncertainties to provide an order-of-magnitude guide.
To obtain an average measure of the uncertainties for the estimated parameters of the core component, the mean of the estimated flux densities and the mean of the estimated FWHM angular sizes were computed for each of the sessions for a sample of 25 sources that were observed in all of the sessions. The standard deviation over the mean values obtained from each of the sessions for the sample of 25 sources gives an estimate of the uncertainty. Taking the ratio of the standard deviation and its mean gives an estimated average relative uncertainty of 9% for the flux density of the main components, and taking the standard deviation gives an estimated average uncertainty of 0.02 mas for the FWHM size of the main components. Since there could be intrinsic variation over the observing sessions, these estimates of the uncertainty are upper limits. We expect, however, that the contribution due to intrinsic variation is small compared to the observed variation. A similar approach was followed to estimate the average uncertainties for the estimated parameters of the second component, except that a subset of only 19 of the 25 sources, which had a clearly distinguishable second component, was used. The average uncertainties for the estimated parameters of the second component are 14% for the flux density, 0.09 mas for the FWHM size, 0.20 mas for the radius, and 5°.4 for the angle. Note that larger uncertainties are expected if the second component is weak compared to the core.

Model-fitting Approach
In addition to the two-component model fitting in DIFMAP, estimates of the principal angle of extension of the total emission region (a.k.a. "jet direction") were obtained by fitting a line through the locations of the CLEAN components. Fitting CLEAN component locations is less resource-intensive than fitting to visibilities and provides a simple measure of the source's elongation, and the robustness of the model fitting in DIFMAP.
The set of CLEAN components that models the source brightness distribution has three dimensions: the flux density and two coordinates, X, Y (relative R.A. and decl.), with the coordinates being relative to the map center. Thus, by fitting a line through the locations of the CLEAN components, we are actually fitting a linear locus of brightest emission to simplify the inherently two-dimensional brightness distribution. We fit the line by minimizing the sum of the square distances between each CLEAN component and the locus line. This is different to a normal least-squares fit, where there is a set of independent variables, in dimension X, associated with a set of observed, dependent ones in dimension Y, and the line of best fit is the one that minimizes the squared sum of the Y-distances between the fitted line and observed Y-values.
The normal method of least-squares fitting does not solve our particular problem, where we have to minimize the sum of the square distances between each CLEAN component and the fitted locus line. One method of solving this problem is to first fit the data assuming the X-coordinates as independent variables and the Y-coordinates as dependent variables and the fit residuals being represented by the vertical (Y-direction) distance between the observed data points and the fitted curve and then, second, to fit the data assuming the Y-coordinates as independent variables and the X-coordinates as dependent variables, with the fit residuals being represented by the vertical (X-direction) distance between the observed data points and the fitted curve. The best fit to the data is then taken as the arithmetic mean between the two fits. However, the unsatisfactory nature of this simplified approach is well known and was evident from the results obtained from applying this method to the image CLEAN components.
A more sophisticated solution is that of "fitting by rotation," where the points in the XY-plane are incrementally rotated clockwise or counterclockwise through 180°by some small angle δθ about the origin. A normal least-squares fit to the data is done for each rotation, and the best fit is the one with the smallest fit residuals. However, since the goal here is merely to get some estimate of the direction of the overall extended emission, a simplified version of this approach is applied in which the XY-plane is rotated so that the X-axis lies on the straight line between the origin and CLEAN component farthest from the phase center. This method was used to perform both an unweighted and flux-density-weighted fit, going through (0, 0), to the CLEAN component locations. Lastly, since this method would not be resistant to outliers caused by weak emission and/or residual errors, all CLEAN components with a flux density of s <5 im and at distances of >10 times the FWHM beam size (where the FWHM beam size is the geometric mean of the beam major and minor axes), were excluded from the fit. In addition, the fit was only performed if there were at least three CLEAN components. The angles obtained from this method agree reasonably well (within a few degrees) with those obtained from the model fitting in DIFMAP, even for sources that appear to be core-dominated from their images (as shown in Figure 14). Where the angles do not agree, it is mostly for sources with either bent jets or with separate jets at large angles to each other (as shown in Figure 14). Thus, this method also allows for easy identification of sources with unique or complex structure, i.e., when the angles do not agree, since it is more sensitive to the direction of the CLEAN component farthest from the phase center. The three sources in Figure 14 were chosen to represent a wide range of source morphologies: (1) core-dominated with weak extended emission, (b) bright secondary components, and (3) curved jets.
The BJ083A1 image of the source J0334+0800 (IVS B0332 +078, ICRF J033453.3+080014) shows a bright second component at a radius of 0.72 mas from the phase center and at an angle of ∼−54° (Figure 14(a)). The spread in position angles obtained from the unweighted and flux-densityweighted linear fit to the image CLEAN component positions, and the estimated angle obtained from the two-component model fitting in DIFMAP is only 0°.5. The image from the UD001I session (Figure 14(b)) suffers from poor u,v-coverage and resolution, preventing the two components from being distinguished. However, the model fitting in DIFMAP does pick up the two components, and the fitting to the CLEAN component locations gives position angles close to those obtained from the model-fitting of this source at other epochs, showing the value of model fitting for identifying individual components, even if they are not distinguishable in the image.
The BJ083C image of the source J0237+2848 (IVS B0234 +285, ICRF J023752.4+284808) shows a weak extended outflow with a position angle around −30°that appears to curve toward the west and then turns toward the north (Figure 14(b)). A spread of more than 14°is seen between the position angles obtained from the unweighted and fluxdensity-weighted linear fit to the image CLEAN component positions and the estimated angle obtained from the twocomponent model fitting in DIFMAP. For this source, a higherorder polynomial would most probably provide a much better fit to the CLEAN component positions than the linear fit. Similar results are seen in the maps and from the model fitting of this source at other epochs.
J1616+0459 (IVS B1614+051, ICRF J161637.5+045932) is a rare example of a source that appears to have two jets at large angles to each other. Session UD001U shows the brighter of the two at an angle of about −53° (Figure 14(d)), and the weaker more extended jet at an angle of around −155°. Since the model fitting to the CLEAN component positions is more sensitive to the maximum extent of the emission, the line of best fit is in the direction of the weaker jet, while the estimated angle from the model fitting in DIFMAP is in the direction of the Notes. This table is published in its entirety in machine-readable format. A portion is shown here for guidance regarding its form and content. Table columns: (1) J2000 source name, (2) session name, (3) reduced χ 2 from the two-component model fitting, second-brightest component. The two jets are also seen in the maps of these sources at other epochs.

Notes on Individual Sources of Interest
This section provides short notes on sources that have significant structure that may be of interest to astronomers. These sources are not ideal for use as calibrators or reference sources and should be used only with caution.

Triple Components
The following sources all show three distinct, bright components, for example J0403+2600 (IVS B0400+258, ICRF J040305.5+260001) shown in Figure 15 There is only one image of this source, which shows three bright components with the core component being the strongest followed by the slightly weaker third component (see Figure 15(a)). The image has SI = 3.8, C 2 = 0.4, and E = 0.7 mas (recalling that on average our sources are relatively compact, with median values of SI, C 2 , and E being 1.7, 0.9, and 0.2 mas, respectively). J1602+3326 (IVS B1600+335, ICRF J160207.2+332653).
The higher-quality images of this source (eight out of the 11 images) all show three components, with the second and third components brightening and fading over time. These eight images have SI > 3, C 2  0.5, and E  0.8 mas.

Close, Near-equal Doubles
The following sources all show two closely separated, nearequal components, for example J0006-0623 (IVS B0003−066, ICRF J000613.8−062335) shown in Figure 15 Four of the nine images show two components of nearequal size that are separated by 1.1 mas (see, e.g., Figure 15(b)), with the second component fading and brightening over time. These four images have SI > 3.8, C 2 < 0.6, and E  0.5 mas. The remainder of the images are of poorer quality, with SI < 3.6, and the two components are not distinguishable in these images.  Figure 12(c)). However, the model fitting in DIFMAP picks up the two components in all but one of the images. All 18 images have C 2  0.6 and E  0.4 mas, and 13 of the images have SI  3.

Bright Second Components
The following sources all have a relatively bright, wellseparated or clearly distinguishable second component, for example J1242+3720 (IVS B1239+376, ICRF J124209.8 +372005) shown in Figure 15  All four images show a close second component at a distance of ∼0.9 mas that fades and brightens over time.
All eight images show a second component that fades and brightens over time. The SI values range between 2.3 and 3.4. The C 2 values range between 0.6 and 0.8, and the E values between 0.5 and 1.2 mas. J0557+2413 (IVS B0554+242, ICRF J055704.7+241355).
All three images show a close, bright second component at a distance of ∼0.8 mas, with SI ∼ 3.2 and E = 0.4 mas for two of the images, with C 2 = 0.5 for the one and C 2 = 0.6 for the other. The third image has SI = 3.7, C 2 = 0.5 and E = 0.4 mas. J0725−0054 (IVS B0723−008, ICRF J072550.6−005456).
All six images show an extended second component that fades and brightens over time. All six images have SI > 3.3, with C 2 ranging between 0.2 and 0.5, and E between 1.1 and 2.2 mas. One of the images has SI = 5.3, with C 2 = 0.2, and E = 1.1 mas. J1024+1912 (IVS B1022+194, ICRF J102444.8+191220).
This source is an ICRF-3 defining source. All four images show a bright second component at a distance of ∼0.9 mas. All of the images have an SI of ∼3.4, and C 2 of ∼0.7, with E ranging between 0.3 and 0.4 mas. J1217+5835 (IVS B1214+588, ICRF J121711.0+583526).
Five of the eight images clearly show a second component that fades over time, while a third component that brightens over time, emerges on the opposite side of the core. The other three images are of much lower quality, and it appears that the source structure was not fully recovered in these images. The SI ranges between 2.7 and 5, with C 2 ranging between 0.4 and 0.7, and E between 0.7 and 1.0 mas. J1242+3720 (IVS B1239+376, ICRF J124209.8+372005).
This source is an ICRF-3 defining source. The images show very weak extended emission that brightens significantly over time into a well-separated second component. One of the seven images is of much lower quality, and it appears that the source structure was not fully recovered in this image. The SI ranges between 2.5 and 4.0 over all seven images, with C 2 ranging between 0.4 and 0.7, and E between 0.3 and 0.8 mas. J1448−1620 (IVS B1445−161, ICRF J144815.0−162024).
All but one of the seven images show a well-separated, second component that fades slightly over time. The SI ranges between 3.1 and 4, with C 2 ranging between 0.5 and 0.7, and E between 0.8 and 1.4 mas. J1500+4751 (IVS B1459+480, ICRF J150048.6+475115).
All but one of the six images show a close, relatively bright second component at a distance of ∼1.1 mas. The SI ranges between 2.8 and 3.3, with C 2 ranging between 0.5 and 0.7, and E between 0.4 and 0.6 mas. J1713−3226 (IVS B1710−323, ICRF J171350.7−322612). A second component is visible in three of the images, with SI ranging between 3.0 and 3.6, C 2 between 0.5 and 0.6, and E between 0.6 and 0.7 mas. This second component is not visible in the fourth lower-quality image, which shows a compact core with SI = 1.1, C 2 = 0.9, and E = 0.1 mas. J2001+3323 (IVS B1959+332, ICRF J200142.2+332344).
Six of the eight images show a strong second component that brightens over time. The SI over these six images ranges between 3.6 and 4.6, while C 2 ranges between 0.4 and 0.5, and E between 0.6 and 0.9 mas. The other two images are of poorer quality, with one showing a compact core with a very weak second component, with SI = 3, C 2 = 0.6, and E = 0.4 mas. The other shows an extended core component with SI = 4.7, C 2 = 0.7, and E = 1.5 mas. The very high SI value of 4.7 is likely the result of artifacts in the image due to the poor u,v-coverage. J2038+5119 (3C 418, ICRF J203837.0+511912). The 26 images show bright emission protruding from the core that first brightens and later fades as it moves farther away from the core into a well-distinguished second component. The second component is clearly visible in all but two of the images that are of poorer quality. The SI ranges between 2.2 and 3.5 over the 24 higher-quality images, with C 2 ranging between 0.5 and 0.7, and E between and 0.3 and 1.1 mas. J2131−1207 (IVS B2128−123, ICRF J213135.2−120704).
Twenty of the 27 images show a close second component at a distance of ∼1.2 mas, with SI ranging between 2.6 and 3.3, C 2 between 0.4 and 0.5, and E between 0.4 and 1.1 mas. The other seven images are mostly of poorer quality and show either a compact core or a core with a slight bulge in the direction of the second component. J2136+0041 (IVS B2134+00, ICRF J213638.5+004154). A bright, well-separated second component is seen in 14 of the 15 images. One of the images is of poorer quality and shows only a compact core. The SI ranges between 2.4 and 4.2, with C 2 between 0.4 and 0.7, and E between 1 and 1.1 mas. J2354+4553 (IVS B2351+456, ICRF J235421.6+455304).
The 13 images show very weak extended emission that brightens over time to reveal a bright second component. The SI ranges between 1.7 and 3.7, with C 2 ranging between 0.6 and 0.8, and E between 0.4 and 1.1 mas.

Bright Extended Emission
The following sources all show bright extended emission that appears connected to or blended with the core, for example J1550+0527 (IVS B1548+056, ICRF J155035.2+052710) shown in Figure 15 Very weak extended emission is seen in the first two epoch images, while later images (apart from two images of poorer quality) show extended emission protruding from the core and growing brighter over time. The SI ranges between 2.7 and 3.4 over all nine images, with C 2 ranging between 0.4 and 0.8, and E between 0.2 and 0.5 mas. J0433+0521 (3C 120, ICRF J043311.0+052115). The 24 images show very weak extended emission connected to the core that brightens significantly in the fourth and fifth epoch image, and then fades again. The SI ranges between 1.4 and 3.3, with C 2 ranging between 0.6 and almost unity, and E between 0.1 and 0.5 mas. J0750+4814 (IVS B0746+483, ICRF J075020.4+481453).
Eleven of the 13 images show relatively bright, extended emission, connected to the core. Two of the images, which are also of poorer quality, show only a compact core. The SI ranges between 2.5 and 4 over all 13 images, while C 2 values range between 0.5 and 0.7, and E between 0.4 and 1.9 mas.
J0927+3902 (4C 39.25, ICRF J092703.0+390220). The images all show bright extended emission protruding from the core. Thirteen of the 17 images have SI > 3. The SI ranges between 2.5 and 4 over all 17 images. The four images where SI < 3 are all of poorer quality, and the extended emission appears more blended with the core component. The values of C 2 range between 0.4 and 0.6 over all 17 images, with E ranging between 0.4 and 0.6 mas, except for one of the poorer-quality images where E = 0.2. J1550+0527 (IVS B1548+056, ICRF J155035.2+052710).
Twenty-four of the 27 images show a long, narrow tail of bright, extended emission protruding from the core that appears to fade over time (see, e.g., Figure 15(d)). The SI over these 24 images ranges between 1.6 and 3.6, with C 2 ranging between 0.5 and 0.7, and E between 0.6 and 0.9 mas. The extended emission is not visible in three of the images that are also of poorer quality and that show only a compact core. These three images have SI = 0.8, 1.2, and 2.7, respectively. The higher SI value of 2.7 is likely the result of artifacts in the image due to the poor u, v-coverage. J2203+3145 (IVS B2201+315, ICRF J220314.9+314538).
The images show bright extended emission connected to the core that almost completely fades over time. The SI ranges between 2.8 and 3.3 over the six images, with C 2 ranging between 0.6 and 0.8, and E between 0.4 and 0.6 mas.

Complex Structure
The following sources appear to have a complex structure that is difficult to image with snapshots only and the automated imaging pipeline in DIFMAP, for example J0432+4138 (3C 119, ICRF J043236.5+413828) shown in Figure 15 Most of the images show a poorly modeled, unresolved core with s im as high as ∼2 mJy beam −1 for some of the images (recalling that the median s im is 0.7 mJy beam −1 for our images). The estimated residual calibration error, f, ranges between 0.15 and 0.34 radians over the six images (recalling that f is typically less than 0.08 radians for our sources). The values of SI range between 1.4 and 4.2. J0121+1149 (IVS B0119+115, ICRF J012141.5+114950).
Three of the four images show poorly modeled, bright, extended emission with artifacts clearly visible in the maps (e.g., components surrounded by a "bowl" of apparently negative flux density). The values of s im range between 2 and 7 mJy beam −1 over these three images, with f ranging between 0.12 and 0.31 radians. The fourth image has s = 1.6 im mJy beam −1 and f = 0.05. The SI over all four images ranges between 3.8 and 4.5. J0432+4138 (3C 119, ICRF J043236.5+413828). All four images show poorly modeled emission with s im ranging between 3 and 5 mJy beam −1 , and f values ranging from 0.35 to as high as 0.84 radians. The SI ranges between 2.3 and 2.9. It is clear that the structure is complex and too difficult to recover with snapshots only (see, e.g., Figure 15(e)). J0644+3914 (IVS B0641+392, ICRF J064453.7+391447).
This source is an ICRF-3 defining source. Although f and s im are not excessively large, with f < 0.13 radians and s < 1.1 im mJy beam −1 , there seems to be no consistent trend in the emission features seen from one image to the next, with artifacts clearly visible in some of the images. The SI values range between 2.5 and 3.7 over the ten images. J0731−2341 (IVS B0728−235, ICRF J073106.6−234147).
Four of the five images show poorly modeled extended emission with artifacts clearly visible in some of the maps.
This source is an ICRF-3 defining source. A very bright second component, at a distance of 0.9 mas, is visible in one of the four images, with SI = 3.9, C 2 = 0.4, and E = 0.7 mas. The other three images are of much lower quality, and it appears that the source structure was not fully recovered in these images, with f ranging between 0.17 and 0.27 radians, and s im ranging between 5 and 7 mJy beam −1 . The SI over these three images ranges between 3 and 3.6. J1159−2148 (IVS B1157−215, ICRF J115951.9−214853).
Two of the seven images show a second component that is relatively bright in the one image, with SI = 4.5, C 2 = 0.6, and E = 1.3 mas, and much weaker in the other, with SI = 3.5, C 2 = 0.8, and E = 1.3 mas. The other five images are of much lower quality, and it appears that the source structure was not fully recovered in these images, with f ranging between 0.05 and 0.17 radians, and s im ranging between 1 and 2 mJy beam −1 . The SI over these five images ranges between 1.4 and 4.4. J1832−1035 (IVS B1829−106, ICRF J183220.8−103511).
Only one of the six images has data on baselines >400 Mλ, and shows an extended core component with an FWHM size of ∼0.9 mas. Artifacts are clearly visible in the map, with f = 0.3 radians, and s = 3 im mJy beam −1 , and SI = 3.6. The source has galactic latitude of −0°.63, and we suspect that the source is broadened sufficiently by the interstellar medium (ISM), so that much of the flux is resolved away on longer baselines.

Sources Near the Galactic Center
The following sources suffer from broadening by the ISM from the Galactic Plane, limiting our detections on baselines longer than about 200 Mλ, for example, J1744−3116 (IVS B1741−312, ICRF J174423.5−311636) shown in Figure 15 The visibility amplitude versus baseline-length plots for this source all show the visibility amplitude dropping rapidly with increasing baseline length (see Figure 15(f)). The source is resolved on baselines longer than ∼200 Mλ. The estimated FWHM angular size of the core component is ∼1 mas (recalling that the median value of the estimated FWHM angular sizes for our sources is 0.14 mas for the core components). J1745−2900 (Sgr Aå). Sgr Aå is believed to be a supermassive black hole at the center of the Galaxy (Ghez et al. 2008;Gravity Collaboration et al. 2018 Figure 16(b)) have similar resolution and are in good agreement with the 22 GHz observations of Hada et al. (2020). However, our 2017 February 23 observations (see Figure 16(a)) are of much higher resolution, and the core-jet structure is more clearly resolved, showing a bright core and a weaker second component at a distance of ∼1.3 mas, consistent with the ∼1.5 mas separation measured by Hada et al. (2020) from their 43 GHz images. J1833−210A (IVS B1830−211, ICRF J183339.8−210340). J1833−210A is the primary image of a strongly lensed quasar, located close to the direction of the Galactic Center. It consists of two main lensed images, A and B, about 1 1 apart, and embedded in a full Einstein ring. Millimeter observations by Muller et al. (2020) identified the much weaker, pointlike source, located about 0 3 away from image B, as the third lensed image of the system. The foreground lens is a nearly face-on, spiral galaxy. Earlier VLBI images of J1833−210A at 2.3, 8.4, 15, and 43 GHz reveal a highly variable, complex structure, with significant differences between the two main lensed images, A and B (e.g., Garrett et al. 1997;Guirado et al. 1999). In contrast, more recent 43 GHz VLBA observations, of eight epochs taken over 14 weeks, show a bright, compact core and a weak extended jet for both A and B (Jin et al. 2003). However, large variations in the angular separation between the core positions of A and B, of up to 87 μas between successive epochs, are seen. Our 24 GHz images of J1833−210A (lensed image A) all suffer from poor u,v-coverage, showing only a single, compact component with an FWHM angular size of ∼0.20 mas, although the amplitude versus u,v-distance plots suggest a more complex structure. The source is detected on the longest baselines and does not appear to be suffering from ISM scattering.
The above list is by no means a comprehensive list of sources that may be of interest to the astronomical community, but rather a list of those sources with the least pointlike structure. There are many other sources in the K-band catalog that show interesting behavior, such as the appearance and disappearance of jet components over time, and sources that show curved jets and jets that change orientation over time. Also included in our sample is the source J0906+6930 (QSO J0906+6930, ICRF J090630.7+693030), which is at z = 5.47 and is one of the highest-redshift blazars known  (Romani et al. 2004). The model fitting in DIFMAP estimates the FWHM angular size of the source to be ∼0.18 mas or less, and some of the images show very weak extended emission to the southwest of the core.
It should be noted that a number of sources observed in the BJ083 and UD001 sessions, in particular sources south of −30°d ecl., suffer from very poor u,v-coverage for all of the images. For most of these, neither the images nor their structure quantities provide adequate or reliable information to assess the quality of the source structure.

Conclusions
The K-band VLBA reference frame program has so far provided 5879 images for 731 sources-the vast majority at multiple epochs. These K-band VLBA sessions have substantially increased the number of calibrators with VLBI images. A detailed analysis of the images has allowed us to determine several quantities that provide useful indicators of the quality of each image and the suitability of each source as a calibrator or reference source, including: an image quality factor, the source flux density and flux density variability, the source compactness and compactness variability, the radial extent of the source structure and radial extent variability, and the structure index and structure index variability. In addition, model fitting has allowed us to determine the FWHM angular size, radial extent and position angle of the core, and a second component in each image. The results of our imaging and analysis are summarized as follows: 1. Analysis of the image parameters indicates that the data are well calibrated, as evidenced by residual phase calibration errors that are typically less than ∼0.08 radians. We find that on average the fractional variation of the flux density from one session to the next is ∼10% over the 28 sessions. On average, the image rms over all of the K-band images is 1.62 mJy beam −1 , with a median value of 0.70 mJy beam −1 . The average S/N over all of the images is 459:1. 2. Our image quality factor, which depends on the magnitude of the residual calibration errors in each image, the image resolution, the elongation of the synthesized beam, and the number of visibility measurements, ranges between 33% and 95% for our sources, with an average value of 61%. We find that the lowerquality images (Q < 60%) are mostly a result of poor u,vcoverage, from observations of sources south of −30°d ecl., and from sessions UD001I to UD001S, which correspond to the period when the SC antenna was out of the array due to hurricane damage. 3. The distribution of the source flux densities shows a mean value of 560 mJy for the total CLEAN flux density and 350 mJy for the correlated flux densities on longer baselines (>5000 km). Variations over the different observing epochs in the individual source flux densities are 17% on average for the total CLEAN flux density and 26% for the correlated flux densities on longer baselines. These results suggest that the sources are sufficiently bright for astrometric use, likely even at frequencies higher than K band. 4. Our measures of source compactness indicate that on average our sources are relatively compact at K band, with 87% of the flux density in the core, and that the variation in the compactness, per source, is relatively low. Measures of the radial extent indicate that, for our sources, most (95%) of the emission is generally within ∼1 mas of the brightest pixel in the image. The median value of 120 μas measured for the flux-densityweighted radial extent is similar in size to the astrometric agreement between the ICRF-3 S/X and K-band frames. On average, the per-source temporal variation in the flux-density-weighted radial extent is relatively low. The distribution of the structure index shows a mean value of 1.7, indicating that on average, the source structure has only a small effect on delay measurements. We also found that the variability of the structure index for a given source over different observing epochs was low. 5. Results from the model fitting show that on average, the core components are relatively compact, with a mean value of 0.14 mas for the estimated FWHM angular size. Where a second component is significant, it is generally within 1 mas from the core. 6. Care should be taken when evaluating the source structure and image-based quantities and their variability over time. We find that in some cases, poor-quality images may suggest a compact source with a low structure index, when in fact it is a double or has some extended emission, and that images with large residual calibration errors may have spuriously large structure index values. It is therefore recommended to not only consider the source structure quantities when evaluating the astrometric suitability of a source, but also to visually inspect the images and visibility amplitude versus baseline-length plots over all epochs, and in particular to take into consideration the image quality factor of each image. Similarly, care should be taken when interpreting the model-fitting results from noisier observations that suffer from poor u,v-coverage. 7. We identified 45 sources (6% of the 731 sources) that should be handled with extreme care when used as a calibrator or reference source, including: two gravitational lenses, two sources near the Galactic Center that suffer from ISM broadening, two sources with bright triple components, four sources with two closely separated, near-equal components, 20 sources with bright second components, seven sources with bright extended emission connected to or blended with the core, and eight sources with a complex structure that appears to be too difficult to image with snapshots only. Six of these sources are also ICRF-3 defining sources, including: J0644+3914 (IVS B0641+392), J0745−0044 (IVS B0743−006), J1024+1912 (IVS B1022+194), J1419 +5423 (IVS B1418+546), J2216+3518 (IVS B2214 +350), and J2230+6946 (IVS B2229+695).
Future Prospects: The work of characterizing calibrator sources has continued since 2018 July (last session covered in this paper) up to the time of publication. These new sessions have and will continue to contribute to the monitoring of the evolution of source strength, morphology, and angular coordinates. In addition, since mid-2018, the source list has been expanded by more than 200 sources. We are seeking to do more than just extend our observation history by implementing a roadmap to continually improve the quality of our observations (de Witt et al. 2019). Our goals are to: (1) improve sensitivity, (2) improve u,v-coverage, and (3) improve calibration. Sensitivity was increased starting in 2019 November by doubling our data rate from 2 to 4 Gbps by changing from single to dual (RCP, LCP) polarization. Moving to dual polarization will allow our imaging to extend from simple intensity to full Stokes parameters. The next increase in sensitivity being pursued is to augment the VLBA with the 50 m Large Millimeter Telescope (LMT) in Mexico (Kurtz et al. 2020). If successful, baselines from the LMT to the VLBA will have approximately double the sensitivity of current VLBA baselines. These new baselines would also enhance u,vcoverage. In parallel, we are seeking collaborators to improve coverage with North-South baselines from East Asia to Hobart, Tasmania and Yebes, Spain to Hartebeesthoek, South Africa. We have already obtained fringes on this last baseline, and it is approved for full 24 hr sessions in mid-2022. In the longer term, we are exploring a next-generation broadband receiver for the VLBA, which would cover from 8-36 GHz (Fung et al. 2020). A prototype is scheduled for testing at the Owens Valley VLBA site in late 2023. If this prototype becomes operational and is deployed across the VLBA, it would enable up to four simultaneous bands (X, Ku, K, Ka) each band up to 4 GHz in width and having dual polarization. The increase from 0.5 to 4 GHz would provide almost a factor of 3 increase in sensitivity. The capability for multiple simultaneous bands would allow for direct ionospheric calibrations, thus removing one of our major error sources . In summary, the success of the work described in this paper has prompted continuing observations and ongoing improvements in data quality through enhanced sensitivity, u,v-coverage, and calibration. The images presented in this paper, when combined with future images from our continually improving observation program, should provide a well-documented set of calibrators as well as a rich source of astrophysical insights into AGNs. For all of these reasons, we see a bright future for K-band CRF work.