A procedure for making high dynamic-range radio images: Deep imaging of the kiloparsec-scale radio structures of a distant blazar, NRAO 530, with JVLA data

Using JVLA data of Sgr A* and NRAO 530 (as a calibrator) at 5.5, 9 and 33 GHz, we developed a procedure for the reduction of wideband data. We have demonstrated that, correcting for residual interferometer errors such as residual delays, astronomers can now achieve high-fidelity radio images with a dynamic range exceeding 1,000,000:1. We outline the detailed procedure, noting that it can have broad application to the analysis of broadband continuum observations. We apply this procedure to observations of a distant blazar, NRAO 530, revealing its structures in unprecedented detail. Our 5.5-GHz image shows that the structure of NRAO 530 is prominently characterized by a moderately curved western jet terminating at a hot spot. Close to the radio core, an abrupt bending of the jet is revealed in the high-resolution (<100 mas) images at 33 GHz, showing an evolution of the position angle of the jet from the north at the VLBI scale (50 mas, or ~400 pc projected), increasing toward the west at the larger VLA scales (1 arcsec, or ~10 kpc). The continuation of the jet axis drift forms the curved western jet extending out to 20 arcsec. In contrast, a faint and broad counter-jet is present on the eastern side with a curvature antisymmetric to the western jet. The eastern jet terminates at a bright hotspot, forming an edge-brightened diffuse lobe. The observed contrast in brightness between the western and eastern jets suggests that the jets on the VLA scale are mildly relativistic. The radiation from the western jet is boosted while the radiation from the receding eastern jet is plausibly suppressed owing to the relativistic Doppler effect.


INTRODUCTION
As a number of radio interferometer systems increase sensitivity by employing large increases in bandwidth, several important corrections must be applied due to time variations of various instrumental parameters during the observation. An example is antenna-based residual delays that vary across the band and that change with time. For the wideband data distorted by these time-variable issues, it becomes difficult to restore the true information in the process of imaging the continuum structure of a radio source. Consequently, without corrections for such residual effects, the dynamic range of the resulting images is limited.
In previous narrow-band observations, the dynamic range of continuum observations was typically sensitivity-limited, so it was challenging at best to learn about the full spatial structure of the brightest radio sources such as blazars, distant radio galaxies and QSOs characterized with dominant bright radio cores. Such objects have often served the purpose of being calibrator sources for interferometer data inasmuch as they appear largely as bright point sources. Historically, the emission from their extended environs, as well as background radio sources, were neglected, except at the highest frequencies, where the bright emission components might be partially resolved.
When extended-source signals from such high-contrast regions are sampled with recent wideband techniques, the resid-ual errors (RE) left over from standard calibration procedures from previous eras challenge the imaging reliability, particularly for the faint features in the vicinity of bright radio sources, features that are often confused with artifacts resulting from the uncorrected RE. The dynamic range (DR) of the image is therefore limited, and the image fidelity is degraded owing to the residual errors in phase and amplitude. DR is a good indicator of image fidelity -the difference between any produced image and the true image (Perley 1999), so DR can be used as a measure of image quality. In the presence of bright radio sources, DR is particularly important for reliably capturing the distribution of much fainter extended emission. For example, NRAO 530 was listed in the VLA calibrator catalog accompanied by images showing little more than a dominant radio core at 14.9 and 8.46 GHz (see Fig. A1 of this paper). Faint extended structure was revealed by later VLA observations at 1.46 GHz with the A-array, showing east-west radio lobes and jets (Kharb et al. 2010). With an angular resolution of 1.5" and noise level of 0.34 mJy beam −1 , this observation achieved a DR of∼18,000:1.
Motivated by our JVLA observations of the Galactic center source, Sgr A (Zhao, Morris & Goss 2016;Morris, Zhao & Goss 2017), we developed a procedure for restoring the emission structure of the calibrators used for our programs to a substantially greater depth than had previously been accomplished. The image DR was improved by iteratively converging the reference model to the true structure by applying various nonstandard corrections, in addition to eliminating radio frequency interference (RFI). In particular, we corrected for the time-variable residual delays.
In this paper, we illustrate our detailed procedure by applying it to the radio structure of the complex-gain (phase and amplitude) calibrator NRAO 530 -a distant radio blazar -using the JVLA data obtained in support of the observations of Sgr A. Using the best DR image model of the calibrator, the solutions for correcting RE were derived and applied to the Sgr A data (Morris, Zhao & Goss 2017). As by-products, an unprecedentedly deep image of NRAO 530 was also obtained. So while we use this source to demonstrate the wideband capability for RE correction to maximize the DR of the image, we also report a new scientific perspective on the blazar itself.
NRAO 530 is a gamma-ray blazar, associated with an optically violent variable source at z=0.902 (Healey et al. 2008) or D L = 5.8 Gpc (1 = 7.8 kpc) 1 . The radio source was first discovered with the 300-foot telescope of the NRAO (Pauliny-Toth, Wade & Heeschen 1966). The VLBI observations show a jet structure to the north of its compact core, suggesting a prominent example of a periodic oscillation of the jet axis (Lister et al. 2013) or a helical motion (Lu, Krichbaum & Zensus 2011) that the jet traces on scales of milli-arcsec or parsec. The northern VLBI jet is superluminal (Lu, Krichbaum & Zensus 2011;Lister et al. 2016Lister et al. , 2018. On kpc scales, NRAO 530 is characterized as having double radio lobes in the E-W direction. The western lobe shows relatively higher intensity, showing a curved jet connecting the western hot spot to the core. NRAO 530 is one of the MOJAVE 2 blazar sample mapped with VLA data by Kharb et al. (2010). Based on their kpc-scale radio study of the MOJAVE blazar sample, those authors have questioned the simple radio-loud unified scheme linking BL Lac objects (i.e. blazars) to FR-I galaxies and quasars to FR II sources (Fanaroff & Riley 1974;Urry & Padovani 1995) and have also pointed out that a significant fraction of MOJAVE blazars show parsec-to-kiloparsec scale jet misalignment. NRAO 530 seems to fall into the FR I/II category because of the ∼ 90 • misalignment between the VLBI jet and the VLA structure (Kharb et al. 2010).
In general, given the bright cores of blazars and typically large distances, the relatively faint extended structure of their jets, radio lobes and halos of radio cores are hidden under the detection limits of the old-generation narrow-band telescopes. As an illustrative example for a distant blazar, deep imaging of NRAO 530 with the JVLA endowed with wideband capability is necessary to reveal the puzzles concerning the connection of the morphology on the VLA scale and the powerful core observed with VLBI.
The rest of the paper is outlined as follows. The details of the procedure for correcting RE, and achieving high-DR imaging for wideband radio interferometer array data are given and discussed in Appendix A. In the main text, section 2 describes the data and new images of NRAO 530. Our findings from the blazar NRAO 530 are presented in section 3, and the astrophysical implications are discussed in section 4. Section 5 summarizes our astrophysical results.
Along with the dramatic hardware improvements of the JVLA, data reduction software packages are being developed emphasizing the wideband capability that allows astronomers to make images much deeper than have previously been possible. However, many successive corrections for possible data issues are needed in order to construct VLA images approaching theoretical image fidelity of the equipped broadband instruments (see Appendix A).

Observations of NRAO 530 and datasets
In the recent study of Sgr A with the JVLA, we developed procedures to demonstrate the possibility of correcting for the data errors and minimizing the instrumental and atmospheric effects. In the following text, we discuss the data that were used for the images of J1733-1304 (hereafter, NRAO 530) and J1744-3116, both used as calibrators for observations of Sgr A. The relevant C-band, X-band and Ka-band data sets used in this paper are summarized in Table 1. We have applied the technique of RE-correcting and DR-imaging to the radio study of the Galactic center; e.g., Morris, Zhao & Goss (2017) have published a paper on deep imaging the non-thermal radio filament (SgrAWF) near the Galactic black hole. More papers on the nonthermal radio filaments and faint compact radio sources in the radio bright zone (RBZ) at the Galactic center will follow.

Images of NRAO 530
2.2.1. C-band images at 5.5 GHz Our C-band observations are described in a sequence of papers (Zhao, Morris & Goss 2013Morris, Zhao & Goss 2017). NRAO 530 was used as a gain calibrator. Using CASA 3 , we followed the standard calibration for a continuum source. In order to achieve a high dynamic range (HDR) in a highresolution image (A-array), in addition to the VLA standard calibrations, two further corrections must be taken: 1) correction for residual delays 4 across each of the sub-bands and 2) removal of flux-density variability (Zhao et al. 1991;Zhao, Morris & Goss 2016;Morris, Zhao & Goss 2017). At radio wavelengths, the former is related to antenna-based instrumental issues of unclear origin (see Appendix A.4. for further discussions), and the latter is intrinsic to AGNs.
We processed the C-band data taken in the A-and B-array at 5.5 GHz with the JVLA, following the procedure described in the flow-chart exhibited in Fig. A4 of the Appendix. After the first five steps with antenna-based corrections, the rms noise of the image with robustness weighting R=0 is 4.5 µJy beam −1 in a region far away from the core (> 15 ), giving a DR exceeding 1,000,000:1. However, within the central 1 region, the rms noise is about a factor of 10 higher, indicating residual errors remaining in the data, including minor baseline-based closure issues. Thus with the model produced at step 5, we performed a baseline-based correction with a solution interval of 30 seconds. Fig. 1 shows the final image made with natural weighting 12A-037 (Zhao) 2012-03-29 (1) The JVLA program code and PI name.
(2) The calendar dates of the observations. (3) The Array configurations. (4) The JVLA band codes † . (5) The frequencies at the observing band center. (6) The total bandwidth of the observations. (7) The coverage of hour angle (HA) range for the data.
(8) The time on sources. (9) The source name of gain calibrators included in the data. † The band codes "C", "X" and "Ka" used for the JVLA correspond to the receiver bands in the frequency ranges of 4.0 -8. by combining the A-and B-array data. An rms noise of 3.5 µJy beam −1 is achieved in the region without source emission (>15" from the image center) while the rms noise is 4.5 µJy beam −1 near the core (∼ 4") 5 , comparable to the expected thermal noise of ∼3 µJy beam −1 . Given an observing frequency, the positional errors caused by residual delays are proportional to the size of the synthesized beam (see Appendix). For the C-array data at 5.5 GHz, contaminated with the mixed issues of larger uncertainties of the core position owing to the residual delays, additional flux density from extended emission and lower-level RFI, we initially had a difficulty to converging the RE-correction and high-DR imaging simply following the procedure outlined in the flowchart (Fig. A4). A special procedure is developed with more detailed sub-steps to complement each of the major correction cycles. A detailed description for the reduction of the C-array data is given in Appendix A.5. At step M8, a decent image was obtained from the antenna-based correction process, showing a N-S extension. At this stage, it is not clear whether the weak N-S feature is a pattern of possible faint residual sidelobes likely caused by the postulated baseline-based residual errors. At the final step of the baseline-based correction with the imaging model of the radio core and E-W lobes and jets excluding any N-S structure, we solved baseline-based solutions with intervals of 120 s, 60 s and 30 s and applied them separately to the data. This final step can effectively eliminate the significant basedline-based errors in comparison to the noise levels corresponding to the solution intervals. Therefore, the final image model is guatanteed to be produced from antenna-based signals. The three trials with different solution intervals all give the same source structure, revealing the faint N-S extended emission. This necessary test excludes the possibility that the N-S pattern is a residual sidelobe caused by the baseline-based closure issues on some individual baselines in certain periods, e.g., temporary issues of some correlator chip. Therefore, we conclude that the N-S extension is most likely the true weak emission extended from the core. Fig. 2 (top-left panel) shows the C-array image at 5.5 GHz, revealing a lower brightness feature elongated north-south, labelled as N-extension and S-extension, with an rms noise of 6 µJy beam −1 in the region no contamination from source emission (20" away from the core).
Then with the point-source at the core, we aligned the Carray data with the A-and B-array data. The flux density of the core was reduced down to a common value of 0.18 Jy for all six data sets by subtracting a point source model from the phase center of the data to below the contamination level spread from the wing of the clean beam. A Gaussian taper function was FIG. 1.-The image of NRAO 530 made with natural weighting, i.e. weighting all the visibility data equally, by combining four sets of data observed at 5.5 GHz using the JVLA in the A-and B-arrays in July of 2014 and March/April of 2012, respectively (see Table 1). The peak intensity and rms are 4.715 Jy beam −1 and 3.4 µJy beam −1 , respectively, corresponding to a dynamic range of 1,300,000:1. The wedge (right) shows the color scale (0 − 0.3 mJy beam −1 ) in the display of the extended emission. The synthesized beam is 0.57"x0.41" (−2.4 • ); the ellipse at the left-bottom indicates FWHM of the beam. applied to the combined data with OUTERTAPER 6 = 300 kλ to weight down the long-baseline visibilities. The dirty image was cleaned with the multi-scale and multi-frequency-synthesis (MS-MFS) algorithm (Rau & Cornwell 2011) of CASA and the clean image (Fig. 2, middle) was convolved to a synthesized beam 0.81"×0.69" (4 • ). The rms noise is σ = 5µJy beam −1 . The top-right panel in Fig. 2 shows the central 6"×8" region, revealing a weak extended emission structure stretched out from the core to the north up to 4 synthesized beams (3 arcsec) at an intensity level of 4σ. The combined A+B+C-array image at 5.5 GHz shows evidently the detection of the weak component near the bright core although the brightness of the southern extension is below the detection limit. Our results confirm the N-S extension observed in the A-array image at 1.4 GHz (Fig.  2, bottom) taken from Kharb et al. (2010).

X-band images at 9 GHz
Two gain calibrators, NRAO 530 and J1744-3116, were included in the X-band observations of Sgr A* with the JVLA in the A-array during March and April in 2014 (see Table 1). J1744-3116 is 2.5 • away from Sgr A while NRAO 530 is 16 • away. Applying the calibration solutions derived from closer gain calibrators can more effectively eliminate the residual errors from the data of the target sources, in particular for the time-variable residual delays (see Morris, Zhao & Goss (2017) and the Appendix below). Following the five-step procedure for the RE-correcting and DR-imaging process discussed in Ap-pendix A.4, we performed further corrections for the residual gains, residual delays and flux-density variations after the standard calibration for JVLA data, obtaining deep image models for both calibrators. The final image model of J1744-3116 is present in the Appendix (Fig. A3, right), and is utilized to discuss the consequence of flux-density variation of the core on imaging and is not discussed in the main text. The images of NRAO 530 are discussed as follows. Fig. 3 shows the image of NRAO 530 made with robustness weight R=0 by combining the three data sets taken on 2014-3-1, 2014-4-17 and 2014-4-28. The peak of the image is 4.9 Jy beam −1 and the rms noise in the outer region is 3.6 µJy beam −1 , giving DR of 1,300,000:1 with the synthesized beam of 0.27"×0.19" (−17 • ). Twelve compact components in the western jet (W2-W13) and two components in the eastern jet (E1 and E5) are detected at the level of > 8σ (see Fig. 3). These components are fitted with 2D Gaussian function and the parameters resulting from the fitting are summarized in Table 2. Fig. 4 shows uniform weighted X-band images at three individual epochs, corresponding to the observations at 9 GHz on March 1, April 17 and April 28 in 2014. A comparable rms noise level in the range between 11 -15 µJy beam −1 is achieved. Each of the epoch data sets has been processed independently by following the RE-correcting and DR-imaging procedure (Fig. A4). The resulting images agree with each other. However, the UV data were sampled in different hourangle ranges at the three epochs (Table 1), which resulted in different FWHM of the synthesized beams (see Fig. 4 caption). The detailed difference between the epochal images, e.g., the apparent peak intensity of W4, is attributable to the different The rms noise in the region ∼10" from the core is 6 µJy beam −1 . The synthesized beam is 4.3" × 2.8", PA = 11 • . The lowest positive contour corresponds to brightness temperature T B = 0.7 K. Middle: The image of NRAO 530 made with six data sets observed at 5.5 GHz in the A-, B-and C-arrays. A gaussian taper function was applied to the data, weighting down the long-baseline visibilities with OUTERTAPER = 300 kλ (see footnote 4 in Section 2.1). The final image is convolved to a synthesized beam 0.81"×0.69" (4 • ), with rms = 5 µJy beam −1 . The inner 6"×8" of the core, indicated by a box in the image of the whole source, is enlarged at top-right. The contours are Sp × −10 −4 , 10 −4 × 2 n , where n = 0, 1, 2, ..., 13 and Sp = 0.18 Jy beam −1 . The lowest positive contour corresponds to brightness temperature T B = 1.3K. Bottom: Image of NRAO 530 made with the VLA data observed in the A array at 1.46 GHz with the synthesized beam size of ∼1.5" (Kharb et al. 2010). The peak intensity and rms noise are 6.11 Jy beam −1 and 0.34 mJy beam −1 , respectively, corresponding to a dynamic range of 18,000:1.
UV coverage in the three observations.

Ka-band images at 33 GHz
The VLA Ka-band data set observed on September 11, 2015, was obtained from the NRAO archive (15A-293). The correlator setup for the dataset is described in the footnote of Table  1. The resulting total bandwidth is 8 GHz. The central band frequency is 33 GHz. The data were first calibrated following the standard JVLA procedure. The corrections for the residual errors follow our new method. This data set was used to test and demontrate the procedure for RE-correction and high-DR imaging discussed in the Appendix (A.4). Fig. 5 shows the image at 33 GHz with natural weighting. An rms noise of 8µJy beam −1 is achieved with a resolution better than 0.1", revealing the bending transition from the northern VLBI jet to the western VLA jet.

DESCRIPTION OF THE STRUCTURE OF NRAO 530
3.1. Structure on kiloparsec scales Fig. 1 shows a remarkably detailed image of the radio emission structure related to this blazar. NRAO 530 appears to share the common radio emission structure observed in radio galaxies and QSOs. From the high-resolution (0.5") image at 5.5 GHz made with combined A-and B-array data from the JVLA observations ( Fig. 1) as well as the Gaussian-tapered image at 0.8" resolution with combined A-, B-, and C-array data (Fig. 2,middle), both the dominant bright core and the extended radio structure with an angular size of 10 arcsec show a sequence of slightly resolved emission knots aligned in the wellconfined but curved western jet (W-jet). This structure bridges the hotspot (W12) in the western lobe (W-lobe) with the radio core. The eastern jet (E-jet) is delineated by a faint radiation blob (E2) and a plume of slightly curved faint emission (E3 and E4) connecting to the hotspot (E5) in the eastern lobe (E-lobe) (see the top panel in Fig. 6 for labelling, and also see Fig. 1 and the middle panel in Fig. 2 for the detailed structure). In contrast to the western jet, in which the emission is well confined within a train of bright spots, the eastern jet (or jet path) is represented by a more diffuse emission ridge with no bright spots. The radio morphology is similar to those of radio galaxies and QSOs (Bridle & Perley 1984, for example); the radio core of NRAO 530 dominates the power, with intensity ratios of 610:1 and 79,000:1 with respect to the hotspots in the Eand W-lobes, respectively. At 5.5 GHz, the specific power of the core is about 1.1×10 28 Watt Hz −1 , dominating the power of ∼ 5.6 × 10 27 Watt Hz −1 obtained by integrating over the extended components as revealed in the C-array image (see Fig.  2 and Table 3).

Eastern jet
The 5.5 GHz images (Fig. 1, the middle panel of Fig. 2 and the top panel of Fig. 6) show the slightly curved eastern jet, or jet path, with no bright spots. Fig. 3 shows an image at 9 GHz with angular resolution of 0.2". At this resolution, two (E2 and E4) of the five features along the eastern jet and lobe (E1, E2, E3, E4 and E5) are detected at 9 GHz. All five features are observed at 5.5 GHz with angular resolution of 0.5" (see the top panel of Fig. 6 for labelling). The brightest component in the eastern jet, E1 (S p = 0.12 mJy beam −1 ), is located 0.6" east of the core, and is detected only at 9 GHz in the A-array observations during the period between March 1 and April 28, 2014 (see Fig. 4). The deconvolved size of E1 0.37"×0.24" (PA=148 • ) has been determined from 2D Gaussian fitting; the source is slightly resolved at the angular resolution of 0.2". The component E1 is not detected at 5.5 GHz with angular resolution of 0.5", or cannot be separated from the core owing to inadequate angular resolution. This component is also not detected at 33 GHz. The 3σ limit (S p = 24 µJy beam −1 ) of the A-array image at 33 GHz is consistent with E1 being a newly FIG. 3.-Image of NRAO 530 made with robustness weighting R=0 by combining three data sets from the A-array observations carried out with the JVLA at 9 GHz on 2014-03-01, 2014-04-17 and 2014-04-28. The synthesized beam is 0.27"x0.19" (−17 • ). The image peak and rms noise are Sp =4.915 Jy beam −1 and σ =3.6 µJy beam −1 , corresponding to a dynamic range of 1,300,000:1. The compact components with peak intensity > 8σ are labelled. The contours are 2σ×(-2, 2 n ), for n = 1, 2, 3, ..., 19. The lower-brightness components W14, E2 and E4 are detected at 5.5 GHz as labelled in Fig. 6, top. ejected jet component with a steep spectral index and an intensity below the detection limit at the angular resolution of 0.1". Furthermore, the component E1 is detected in all three epochs of the A-array observations at 9 GHz on 2014-3-1, 2014-4-17 and 2014-4-28. No significant variation in flux density, position or size is observed during this two-month period.

Hotspot of E-lobe
The component E5 appears to mark the hotspot in the eastern lobe where the curved eastern jet terminates, showing the outer edge-brightened structure of the eastern lobe. The surface brightness of E5 is low, with peak intensities at 9 and 5.5 GHz being S p = 30 and 60 µJy beam −1 as determined from the images convolved to a commom beam (0.27"×0.19", PA=−17 • ). The intensity contrast to the core is S E5 p : S Core p = 1:160,000 at 9 GHz and 1:79,000 at 5.5 GHz. At 1.46 GHz, this ratio is less than 1:1,000, determined from the image of Kharb et al. (2010) (Fig. 2 bottom). In addition, the eastern hotspot E5 is trailed by a fainter emission plume elongated 5" towards the northwest. The spectral index determined from the 9 and 5.5 GHz data is −1.4 ± 0.7.

Western jet
Fourteen compact components in the western jet and lobe are identified as W1 to W14. W1 is only detected in the high resolution image at 33 GHz, and is not distinguishable from the core in the JVLA A-array images at both 9 and 5.5 GHz. The components W1 to W12 are located along a curved jet feature, convex towards north. Table 2 lists emission quantities determined from the high-resolution JVLA images. The components are slightly resolved in the A-array images. The 9 and 5.5 GHz data at angular resolution of 0.2" show that the components W2-W11 have peak intensity S p ranging from 0.06 to 3.07 mJy beam −1 at 9 GHz and 0.13-4.70 mJy beam −1 at 5.5 GHz. The spectral index determined from these components at 9 and 5.5 GHz is about α = −1, consistent with synchrotron radiation. W3, located 0.5" northwest of W2, appears to be nearly resolved out at 33 GHz, showing a marginal detection at the level of 4σ. The rest of the components in the western jet appear to be resolved out and not detected at 33 GHz.

Hotspot of W-lobe
After the core, the second brightest component in NRAO 530 is W12, with peak intensities of 5.0 and 7.72 mJy beam −1 at 9 and 5.5 GHz, respectively, marking the hotspot in the western jet/lobe structure. A hotspot usually represents the termination of a jet flow (Blandford & Rees 1974, for example), where the jet flow impacts the ambient medium at the outer edge of the lobe. Similar to the northwestern tail of the eastern hotspot E5, a southwestern extension (W13 and W14) appears to be associated with W12 (see Figs. 1, 2 and 6). We also point out a detailed difference of the W-lobe from the E-lobe as follows. The northwestern tail of the hotspot E5 and the hotspot itself, as well the eastern jet that forms the edge-brightened large lobe all combine to give the appearance that the jet material is flowing back toward the core, at leat in projection. The E-lobe appears to be consistent with the typical FR II morphology. The western jet, on the other hand, turns 90 • to the south after passing through the hotspot W12, and then stretches further out to the west. There appears no indication of a similar continuation of the western jet after the hotspot W12 pointing back towards the core. The diffuse western tail of the W-lobe appears to be more like the morphology as often observed in FR I type jets, such as the mirror-symmetric jets in 3C 449 (Perley et al. 1979).

North-south extension to the core
In the C-array image with low angular resolution of ∼4" (Fig.  2), a lower surface brightness structure (plume) extending 15" north and 13" south of the core is detected at an intensity level in the range between 0.2 and 2 mJy beam −1 . This N-S feature is obvious in the 1.46 GHz A-array image at a resolution of 1.5" from Kharb et al. (2010) (see the bottom panel of Fig. 2). In the higher angular resolution (0.8") image with combined A-,  Col. 1 is component ID. Col. 2 is peak flux densities at 9 GHz / 5.5 GHz. Col. 3 is the surface brightness temperature. Col. 4 is the observing frequency. Col. 5 is spectral index derived from the values listed in col. 2. Col. 6 is total flux density. Col. 7 is the angular offset from the core, RA (J2000) = 17:33:02.706, Dec (J2000) = -13:04:49.55. Col. 8 is the position angle of the offset vector pointing from the core towards the jet components. Cols. 9, 10 and 11 are the sizes along the major and minor axis, and the position angle from 2D Gaussian fitting. The values listed in column 2 and columns 5 to 10 corresponding to 9 GHz were determined the A-array image with a synthesized beam of 0.27"×0.19" (−17 • ) (Fig. 3). The intensities at 5.5 GHz listed in col. 2 were determined from an image with uniform weighting, convolved to the same synthesized beam as the 9 GHz image. The rms noise of 9 µJy beam −1 in the 5.5 GHz image is larger by a factor of 2.5 than that of the 9 GHz image. † The quantities for these components are determined from the natural weighted A+B image at 5.5 GHz with the synthesized beam 0.57"×0.41" (−2.4 • ).
B-, and C-array data at 5.5 GHz, an emission plume with intensity in the range between 0.02 and 0.3 mJy beam −1 is located a few arcsec north of the core. The lower surface brightness emission feature of the southern extension appears to be below the detection limit. Table 3 lists the emission quantities of the north-south extension for comparison with the core and eastern and western lobes. The 5.5 GHz flux density integrated over the northern extension is 7±1 mJy, about 5 percent of the west-ern lobe, including the jet, and the flux density of the southern extension is 3±1 mJy. The flux density integrated over the full N-S extension is about 10% of the total flux density integrated over the overall extended emission, excluding the central core emission. However, the total flux-density from the extended emssion in the north-south extension is only ∼3% of the core flux density at 5.5 GHz and about 10% at 1.4 GHz.
3.2. The northern VLBI jet and western JVLA jet 3.2.1. J -the vector of jet axis Col. 4 is the ratio of core flux density to the integrated flux density for each component. Col. 5 is radio power Pν = St 4πD 2 L /(1 + z) 1+α , assuming D L = 5.8Gpc. † The values are determined from the C-array image at 5.5 GHz (Fig. 2).  6 shows the high-resolution images obtained from the JVLA wideband data at 5.5, 9 and 33 GHz as compared to the VLBA image at 15 GHz. The JVLA images trace the western jet on scales from 200 kpc (5.5 GHz) to 700 pc (33 GHz). At the resolution of 0.1", the components W1 and W2 detected at 33 GHz show an elongation to the NW. This structure begins at the core, linking to the large-scale western jet. The position angles of the major axis (from 2D Gaussian fitting) indicate PA=146 • ±7 • for W1 at an angular offset ∆θ = 0.31" from the core and PA=111 • ±3 • for W2 at ∆θ = 0.69" (see Table 2). The change in PA between W1 and W2 suggests that the projected jet axis bends towards the west as the scale increases from a few 100 pc to a few kpc. At the size scale of the VLBA structure of < 200 pc, the axis of the jet is oriented toward the north (Lu, Krichbaum & Zensus 2011, for example).
To quantitatively investigate the relationship between the jet on the VLA scale and the VLBI jet, we defined a jet vector J as the vector that connects the core to a jet component. The vector J can be quantitatively described by its position angle (PA θ ) and its angular offset from the core (∆θ). For a jet with linear motion or no precession involved, its trajectory is a straight line, i.e., no changes in PA θ are expected as ∆θ increases. Table  2 lists ∆θ and PA θ in columns 7 and 8, respectively, for all components. Fig. 7 (top) plots the position angle PA θ as function of ∆θ. The western (eastern) jet components are represented in green (red); the quantities PA θ and ∆θ on angular scales between 0.1 to 20 arcsec, are determined from the JVLA data (see Table  3). The northern jet components observed with the VLBA (Lu, Krichbaum & Zensus 2011) are in blue. From this plot, the following obvious conclusions can be drawn: (1) The JVLA western jet (green) is the continuation of the VLBI northern jet (blue). The PA of the jet axis J is essentially to the north on the scale of VLBA component (< 50 mas) and gradually changes to west as ∆θ increase, the angular distance. The transition where PA θ changes rapidly occurs in the range between 0.05" and 1" in ∆θ. The northern VBLI jet and the western JVLA jet belong to one astrophysical entity, hereafter referred as the northwestern jet.   (2011) (2) On scales < 50 mas, the PA θ changes by > 50 • . According to proper motion measurements based on multiple VLBA observations by Lu, Krichbaum & Zensus (2011), those authors conclude that the motion on VLBA scales is dominated by the E-W swing of the jet components at a rate of 3.4 • yr −1 . They suggest that the observed angular motion is evidence for helical motion of the jet component on scales from a few mas to a few tens of mas. Such an E-W swing in which PA θ changes by > 50 • is not observed in the JVLA data. Plausibly, the rapid swing motion occurs only at small angular scales < 50 mas and the angular resolution of the JVLA is not adequate to detect such small-scale angular motions, such as jet component rotation around the jet axis or the inferred helical motion. We note that the helical motion of jet material occurred on the VLBI scale appears to be not the same thing as a precession of the jet axis characterized by the jet curvature on VLA scales. Slow precession of jet axis can also result in observed VBLI and VLA jet misalignments.
(3) It is obvious that over life time of the northwestern jet, the PA θ changes about ∼ 90 • over the scale of ∼10" (or ∼100 kpc).
Unlike the western jet in which the bright emission spots clearly delineate the projected jet trajectory, the eastern jet components are determined from the enhanced emission peaks along the relatively diffuse, elongated ridge that likely delineates the eastern jet path. The uncertainties in the parameters of the eastern components could be larger than the errors quoted in Table  2 concerning the issue in identification of jet components, given the strong likelihood that the eastern jet is receding and the radiation from the jet components is therefore suppressed due to the relativistic Doppler effect, assuming implicitly that the jet is relativistic given that the VLBI apparent motions are superluminal.

T -the unit vector of the jet transverse velocity
The apparent projected direction of motion of the northwestern jet can also be described quantitatively. We can determine the direction of the transverse velocity with a unit vector T defined as the direction change of the vector J between consecutive jet components: at the i-th jet component where J is the vector pointing from the core to the jet component, defined in the previous section. The direction of T i can be described by PA Jet , the position angle of T i . For a curved jet, PA Jet is a function of D i Jet , the distance measured along the projected jet trajectory from the core to the location of the i-th component. The quantity D i Jet can be computed approximately as the sum of the lengths of the connecting lines between consecutive jet components from the core to the location of the i-th component: where, i = 1, 2, ..., n (see the caption of Fig. 7 for more details). indicates that the overall VLBI jet displays east-west wiggles as it moves northward. The result is consistent with VLBI observations (Lu, Krichbaum & Zensus 2011). In fact, based on the VLBA observations of NRAO 530 at twenty-seven epochs between 1994 and 2011, the innermost jet PA can be fit to a periodic oscillation with amplitude of 20 • and a period of 9.4 yr (Lister et al. 2013). Then, on the scale transiting from VLBI to VLA, the abrupt northwestward bending of the jet occurs as PA Jet changes −7 • (j) to −54 • (W2) on angular scales from 25 mas to 700 mas. The corresponding change rate is 70 deg arcsec −1 . Then, after the following 7" along the jet route, the jet turns westward at D Jet =7.8" (PA Jet = −90 • at W8) with a rate of 5 deg arcsec −1 . Advancing the next 3.1", at W10 before encountering the hotspot region (W12), the jet turns southwestward at a rate of 14 deg arcsec −1 . Then, the jet appears to terminate at W12 followed by a rapid change in PA Jet , 30 deg arcsec −1 . The rapid change in PA Jet indicates that the jet loses its stiffness (or collimation for straight line jets), likely owing to the impact onto the ambient medium at the hotspot. The jet disruption for FR I source could be also due to an unstable growing wave mode, such as a Kelvin-Helmholtz instability (Blandford & Pringle 1976;Ferrari et al. 1978;Hardee 1979;Zhao et al. 1992a,b, for example).
In short, from the diagram of PA Jet −D Jet , the change in PA Jet seems to scale with jet kinematics, i.e. larger change in PA Jet occurs in the inner region of the core, where the helical motion of the jet material dominates the PA Jet change while the PA Jet change appear to be caused by a slow precession of jet axis over the VLA scales. On the VLBI scale, the change rate in PA Jet drops from 83 deg mas −1 at D Jet ∼1 mas to 2 deg mas −1 at 25 mas. The rapid change in PA Jet on the VLBI mas scale is consistent with the jet undergoing helical motion (Lu, Krichbaum & Zensus 2011) or periodic swing (Lister et al. 2013). However, as compared to the rate on the VLBA scale, the rate on the VLA scale is about three orders of magnitude smaller. In the course of the jet moving from the core to the hotspot W12 over its route length of 10", overall the PA Jet changes about 90 • . We made a linear fit to the plotted data and found that, where the uncertainties of the coefficients in the linear fitting are equivalent to 1.2 • , the rms of the residuals (see Fig. 7). The fitted straight line, a logarithmic spiral, can be expressed in polar coordinates if (PA jet , D Jet ) is substituted by (θ, r): where the variables θ and r are in the units of degree and arcsec and the coefficients a = 0.0502 and b = −17.5 • . A similar logarithmic spiral function can also be fit to the projected trajectory that may serve as evidence for a slow precessing jet from the core in NRAO 530.  Table  4. For the data from the monitoring programs (Aller et al. 1985;Teräsranta 1992;Bower et al. 1997;Robson et al. 2001;Hovatta et al. 2008;Jenness et al. 2010;Aller 2012, and the F-GAMMA monitoring programm 7 ), the range of flux densities from minimum to maximum at the observing frequency is scaled with the vertical dotted line. For individual observations and measurements, the mean values are marked with solid dots with the vertical bars for 1 σ errors. The green solid curve is the synchrotron self-absorption model component (S C ) and the blue line is the optically thin synchrotron emission model component (Sext). The red curve is the resulant sum of the two model components (S C +Sext), derived from a best fitting to the flux-density minimal values as marked by the open circles.

Spectral fitting to the low flux-density state
The blazar NRAO 530 exhibits a remarkable degree of flux density variations over a wide frequency range from 5 GHz to 375 GHz (Aller et al. 1985;Teräsranta 1992;Bower et al. 1997;Robson et al. 2001;Feng et al 2006;Hovatta et al. 2008;Jenness et al. 2010;Aller 2012;, and the F-GAMMA monitoring program 7 ). Superluminal motions of the northern jet components within the VLBI core have been detected with apparent transverse velocities ranging from 2.3c to 26.5c (Lu, Krichbaum & Zensus 2011) which agree with the more accurate determination of maximum 27.37±0.97c and median 13.74±0.40c based on the inner six components by the MOJAVE project 8 (Lister et al. 2013(Lister et al. , 2016(Lister et al. , 2018. Thus, the core-dominated emission is highly enhanced via the relativistic Doppler-boosting effect. Table 4 lists the total flux densities of NRAO 530. Fig. 8 plots the spectrum of NRAO 530, implying two discernible states: (1) the high state characterized by highly fluctuating flux densities at the wavelengths from short centimeter to sub-millimeter; (2) the low state that appears to be delineated by the minimal values of the flux densities determined at each of the observing frequencies ν > 1 GHz. The minimal flux density in each of monitoring programs appears to depend on the length of monitoring time. A model composed of a synchrotron self-absorption component and an optically thin component is used to fit the low-state spectrum. Allowing for the emission from the core, the synchrotron self-absorption can be described as, where ν 0 is the turnover frequency at which the synchrotron optical depth becomes unity, S 0 is the flux density at ν 0 multiplied by a factor of e e − 1 = 1.58, and α c is the spectral index when the core emission becomes optically thin.
The emission from the extended lobes and jets is assumed to be optically thin synchrotron radiation with spectral index α e and the flux density of the extended emission at 1 GHz of S 1 for ν in units of GHz, S ext = S 1 ν αe .
There are five parameters (S 0 , ν 0 , α c , S 1 and α e ) to fit. The spectral index α e can be determined from the extended emission. With α e = −0.94, determined from the VLA images at 1.46 GHz (Kharb et al. 2010) and 5.5 GHz (this paper, see Table 3), the number of free parameters is then reduced to four. We fit the low-state flux-density data in Fig. 8 (open circles) to derive the parameters. The turnover frequency is ν 0 = 0.16 GHz, S 0 = 6.2 Jy or S c = 3.9 Jy at ν = ν 0 , and α c = −0.19 which is consistent with the value determined from the core. The extended flux density at 1 GHz, S 1 = 0.739 Jy, is also consistent with the value extrapolated from the VLA measurement for the extended flux density at 1.46 GHz (Kharb et al. 2010). We note that the core in this paper refers to the region within 50 mas from the source center and the emission from the outside of the region is included in the extended component. The spectral index of α c = −0.19 appears to be consistent with the result derived from the VLBI core (Lu, Krichbaum & Zensus 2011). A few notes on the spectral fitting follow: (1) at ν = 142 MHz, the contributions in flux density from the core and extended emission are comparable. At ν < 142 MHz, the extended emission becomes prominent, e.g. at 80 MHz the flux density from the extended emission is 88% of the total flux density. At ν >> 142 MHz, the core becomes dominont, having 88% and 99% of the total flux densities at 1.4 and 90 GHz, respectively. The corresponding total radio power in the low state is 1.0 × 10 28 and 4.5 × 10 27 W Hz −1 at 1.4 and 90 GHz, respectively.
(2) At most times, the core appears to be in the high flux-density state. Higher variability in flux density tends to occur towards high frequencies. This result is consistent with the VLBI determinations of the variation in superluminal velocity at different angular scales, suggesting that the high-state flux density

IMPLICATION AND DISCUSSION
Extragalactic radio sources with jets and lobes can be classified into two morphological classes, FR I and FR II (Fanaroff & Riley 1974). The FR I jets appear to get weak at large distances from the radio core, e.g., 3C 449 (Perley et al. 1979) while the FR II sources, e.g., 3C 47 (Bridle et al. 1994), show edgebrightened lobes that are characterized by a hot spot. These two distinct classes in apparent radio morphology appear to be related to the powers of the radio sources and absolute isophotal magnitude (Ledlow & Owen 1996). In the logarithmic diagram of the radio power at 1.4 GHz (P 1.4 GHz ) and absolute isophotal magnitude (M), the authors find that a non-zero slope line divides FR I from FR II sources, which corresponds to L radio ∝ L 1.8 opt . Based on their nearby sample (z < 0.5), the FR I/II division appears to occur over a range in P 1.4 GHz ∼ 10 24−26.6 Watt Hz −1 .
On the other hand, studies of the correlation between radio power at 1.4 GHz and optical luminosity (Ledlow & Owen 1996) lead to the conclusion that all radio galaxies live in similar environments in that optical luminosity and other properties of the host galaxy are the most important parameters affecting radio source formation and evolution. Relating the absolute optical R-band magnitude M R of the host galaxy to the central black hole mass M BH (McLure & Dunlop 2001) and the radio power at 1.4 GHz, P 1.4GHz , to the kinetic luminosity of the jet, L jet (Willott et al. 1999), Ghisellini & Celotti (2001) find that the dividing line of the radio power at 1.4 GHz between FR I and FR II corresponds to a constant ratio, L jet /M BH . The statistical analysis of the FR I/II dividing line suggests that for radio sources situated on the FR I/II boundary, higher radio powers imply a larger mass of black holes harbored in their cores (e.g. Ghisellini & Celotti 2001;Willott et al. 1999).

FR II or FR I type source?
Based on 1.4 GHz VLA observations of the extended radio structures associated with the blazars in the complete fluxdensity-limited MOJAVE sample, Kharb et al. (2010) find that a substantial fraction of MOJAVE sources fall into FR I/II category according to their luminosities and morphologies. Either the total power of P 1.4GHz = 1.0 × 10 28 W Hz −1 determined from the spectral fitting to the low flux-density state or the power from the extended emission in kpc scale P ext,1.4GHz = 2.1 × 10 27 W Hz −1 suggests that NRAO 530 has a high optical luminosity if it falls on the FR I/II dividing line (Owen & Ledlow 1994;Ledlow & Owen 1996;Ghisellini & Celotti 2001). Therefore, the energetic jets in NRAO 530 might imply a high accretion rate (Ghisellini et al. 2014) onto a very massive black hole (Rawlings & Saunders 1991) at the blazar core.
However, the radio power at 1.4 GHz from the core is a factor of 5 greater than that of the extended emission even in the low state. The flux density from the core is boosted due to the relativistic Doppler effect as superluminal motions have been observed with the VLBA (Lu, Krichbaum & Zensus 2011;Lister et al. 2013Lister et al. , 2016, becoming highly variable. Taking the median value determined from the inner core components with the VLBA, the apparent transverse velocity v a in units of c (the speed of light) is β a = 14, giving a lower limit on the Lorentz factor of γ min = 1 + β 2 a ≈ 14 , assuming the jet to be viewed at a critical angle θ c . Then, the corresponding view angle θ c , the angle between the direction of jet velocity and the line of sight to the observer, can be determined as θ c = arcsin 1 γ min ≈ 4 • (Lu, Krichbaum & Zensus 2011). Thus, the minimum Doppler factor, δ ≡ 1 γ

(− βcosθ)
, would be δ min = γ min ≈ 14. The enhanced factor of observed luminosity for a relativisitic jet is: δ p with p = 3 − α c for a spherical core and p = 2 − α Jet for an optically thin jet (e.g. Cawthorne 1991;Urry & Padovani 1995). Therefore, the radiation from the core is enhanced due to the Doppler boosting by an amount of γ αc−3 (1 − βcosθ) αc−3 (Kellermann & Owen 1988;Cawthorne 1991;Urry & Padovani 1995). In the case of γ = γ min = 14, the flux densities from the inner components of the core can be boosted by ∼ γ 3 min , about 3 orders of magnitude, assuming α c = 0. The Doppler boosting effect explains why only the northern jet components (approaching the observer) are detected if the emission from the receding southern jet is suppressed by ∼ γ −3 min and the bipolar jet is intrinsically symmetric.
Taking account of the Doppler boosting, we find that the actual intrinsic radio power at 1.4 GHz of NRAO 530 is likely on the order of the power from the extended emission, or a few times 10 27 W Hz −1 . Although the edge-brightened eastern lobe with barely visible jet is very consistent with the FR II morphology, the relatively bright western jet with a faint emission lobe stretching further out appears to be more like an FR I-type radio source. Some of the extremely misaligned MO-JAVE blazar jets could be hybrid morphology sources, with an FR I jet on one side and an FR II jet on the other (Kharb et al. 2010). Our deep images may provide further evidence for NRAO 530 falling into the FR I/II category. However, the classification based on edge-brightening appears to be dependent on the viewing direction with respect to the orientation of jetflow. From a certain azimuth and elevation angle, the western lobe of NRAO 530 could also be viewd as an edge brightened lobe, i.e., FR II morphology.

Mildly relativistic jets on the VLA scale?
On the VLBI scale, the jet axis lies essentially north of the core on angular scales 1-25 mas. The jet axis of the northern jet components in the inner core rapidly swings at a rate of 3.4 • yr −1 in change of the position angle (Lu, Krichbaum & Zensus 2011), that can be modelled with a periodic oscillation of the jet position angle for the innermost component (amplitude of 20 • and period of 9.4 yr) (Lister et al. 2013). On VLA scale, over the anglar scale 0.1-10 arcsec, the position angle of the northwestern jet drifts towards the west. In addition, neglecting the contrast in the brightness of the western and eastern jets, the curvature of the eastern jet appears to be antisymmetric with respect to the western jet. The twin jet in NRAO 530 can be characterized as an S-shaped structure.
The S-shaped morphology differs from the U-shaped narrowangle-tailed and wide-angle-tailed twin-jet sources that are often observed in the environment of rich clusters of galaxies (O'Dea & Owen 1986;O'Donoghue et al. 1993). Such mirrorsymmetric jet structures are usually associated with FR I radio sources, e.g., 3C 449 (Perley et al. 1979). The ram-pressure owing to the motion of the jet-parent galaxies through the intracluster medium is suggested to be responsible for the observed, symmetrically curved radio trails (Begelman, Rees & Blandford 1979). Unlike the collinear twin jet observed in many FR II sources in which the view angle of the jet axis keeps constant, the projected S-shaped jet indicates that the jet axis in those cases likely drifts on the sky slowly over the jet lifetime. Consequently, the view angle of the jet axis also changes over time. The S-shaped antisymmetric jet is a feature predicted by precessing jet models, discussed in detail by Gower et al. (1982) via their extensive numerical simulations and used for the interpretation of the VLBI features observed in AGNs, such as 4C+12.50 (Lister et al. 2003).
For the case of NRAO 530, if the mean intensity ratio of the western jet hotspot (W12) to the eastern one (E5) at 5.5 and 9 GHz, R hspt = 150 is owing to the relativistic Doppler effect, then we can estimate the bulk Lorentz factor γ (e.g. Laing & Bridle 2015). The emission from the western jet, plausibly moving towards the observer, is enhanced by δ 2.94 for a coninuous jet with spectral index α Jet ≈ α e = −0.94. If we further assume that the western jet is viewed at a critical angle and the eastern jet is in antisymmetric to the western jet, the enhancement factor is γ min (1 − β 2 a ) −2.94 for the approaching western jet while the receding eastern jet is suppressed by a factor of γ min (1 + β 2 a ) −2.94 . From the intensity ratio R hspt for W12 and E5, the bulk Lorentz factor of γ min ∼ R 1/2.94 hspt +1 2 = 1.8 can be inferred. Then, the view angle of the hotspot W12 can be estimated as θ c = sin −1 1 γ min = 34 • . Compared to θ c = 4 • , determined for the inner jet in the VLBI core, the large view angle of the hot spot and the curvature of the western jet imply that the jet axis has precessed, presumably in concert with the rotation axis of an accretion disk, changing the critical view angle from 34 • at the earlier epoch to 4 • at the recent epoch. The slow precession of the jet axis could occur if the jet nozzle has been subject to a torque as has been postulated in many theoretical models of accreting black holes (e.g. Ghisellini et al. 2014). The simple model adopted in this paper, which assumes an intrinsic bi-directional symmetry in a mildly relativistic bipolar jet coupled with a slow precession of the jet axis, acoounts for the observed antisymmetry of the curved, large-scale twin jet morphology and the intensity ratio of the oppositely directed jets. The observed properties of the northwestern jet seem to also fit into the theoretical model for the hybrid FR I/II sources, in which the jets become transonic in the core and decelerate to mildly relativistic velocities from moderately relativistic or ultrarelativistic jets (Bicknell 1995).

CONCLUSION
With a procedure developed within the CASA software package, we show that the time-variable residual delay in the JVLA data can be effectively eliminated for deep imaging of faint, extended radio structure associated with dominant cores. Applying the technique that is described in the Appendix, sensitivities or dynamic ranges of images made with wideband interferometric array data can be improved by a factor of several tens to over a hundred. As an example for demonstration, the data for NRAO 530, a gain calibrator used in the JVLA observations of Sgr A at 5.5, 9 and 33 GHz during the period between March 29, 2012 and September 11, 2015, have been processed following the procedure of residual-error correction and high dynamic-range imaging discussed in this paper. Deep images of NRAO 530 at 5.5, 9, and 33 GHz at angular resolutions from 0.1" to 4" have been achieved with dynamic ranges in the order of 1,000,000 : 1. The detailed emission structure revealed on scales of 1-100 kpc suggests that the edge-brightened eastern lobe/jet structure appears to fit with the FR II category. The observed morphology of the western jet, with the western extension, seems to more likely fall into the FR I category. The classification for the western jet perhaps is owing to a specific viewing angle. The apparent antisymmetric, curved twin-jet morphology and the observed contrast of their radiation intensity suggest that the twin jets in NRAO 530 on VLA scales is mildly relativistic and that the view-angle of the jet axis likely evolves over time.
We are grateful to Kumar Golap for advice on the use of CASA. We also thank Ken I. Kellermann and Matt L. Lister for their valuable comments. We want to thank the anonymous referee for providing useful input and suggestions, especially for instructive suggestions on the low-frequency spectrum fitting as well as fitting a logarithmic spiral function to the northwestern jet data. The Very Large Array (VLA) is operated by the National Radio Astronomy Observatory (NRAO). The NRAO is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The research has made use of NASA's Astrophysics Data System.

APPENDIX A A PROCEDURE FOR HIGH DYNAMIC RANGE IMAGING
Wideband capability enabled by recent technology advancements has enhanced the power of interferometer arrays at radio, millimeter and submillimeter wavelengths. Sensitivities of radio telescopes have been dramatically improved in recent years. Source structure and variability of core-dominant radio objects such as blazars and quasars have become noticeable issues in calibration and imaging. Less attention has been paid to the effects of these radio properties in handling radio data when the image dynamic range (DR) is limited by the telescope's sensitivity. For example, the image of NRAO 530 observed with the VLA in 1997 (see Fig. A1) indicates that the blazar is a pointlike source at an angular resolution of ∼1 arcsec when the DR is below 10,000:1. In addition, at sub-arcsec resolution, residual delays (δτ A ), varying with time, become a noticeable issue in high-dynamic range imaging with wideband data taken with modern telescopes. The corresponding phase gradient across a wide frequency band causes sidelobe smearing of a dirty beam which is difficult to clean with the algorithms that have been widely used in radio astronomy (Högbom 1974;Clark 1980;Steer, Dewdney & Ito 1984) as well as the newer algorithms developed for multi-frequency-synthesis (MFS) and multiplescale (MS) imaging (Rau & Cornwell 2011). The effect becomes aggravated and significantly degrades the image quality in the field near a strong compact source. 1. Residual phase delays generate errors in the position of a radio core that is supposed to be placed at the phase center or the delay center. For example, phase errors (∆φ) caused by residual delays (δτ A ) across a sub-band (BW = ∆ν) could be due to pointing errors or can propagate to positional uncertainty of the radio core at the sampling channel frequencies in synthesis imaging with MFS algorithms. For example, given fixed baselines, positional errors (∆θ) on the radio core can be related to the quality of data described by the residual delays (δτ A ) remaining after the calibration process, ∆θ = cδτ A B , where c is the speed of light and B is the separation between two antenna elements, or the baseline length. In the case of JVLA A-array observations at 9 GHz, corresponding to the longest baselines, the positional errors of the core caused by residual errors are ∼0.07 , about a quarter of the synthesized beamwidth, θ syn . The main sidelobe pattern in the dirty image will be smeared owing to the small deviation in phase between the channels, corresponding to the positional errors caused by residual delays. The consequence of the residual delays will limit the benefits of deep cleaning in the process of multiple-frequency synthesis imaging. Given a fixed residual delay error δτ A , the positional error ∆θ is inversely proportional to the baseline length. The shorter baselines of the small array configurations produce larger uncertainties in source positions. On the other hand, the sythesized beam of the array is also inversely proportional to the baseline length, θ syn = λ B , where λ is the observing wavelength. The baseline term is cancelled in the ratio of the positional error to the size of the sythesized beam, namely, ∆θ θ syn = cδτ A λ . Therefore, the issue in deep cleaning of wideband imaging caused by dirty beam smearing due to the residual delay is not dependent on the size of the array configuration, given a fixed observing wavelength. However, for the same array configuration, the issues in deep imaging due to the residual delay become worse for the data observed at shorter wavelengths due to the difficulty in achieving a reliable model from the contaminated visibility data.
2. Vector-averaging sub-band data along with self-calibration can eliminate the residual delay issue but only at the cost of increasing the noise in the data. This approach has been used in handling earlier generation, narrow bandwidth VLA continuum data for each of the IFs, or by creating "ch0" data to form the continuum data from spectral line observations. The "ch0" is simply produced from vector-averaging the channels across a spectral band. Indeed, the frequency-dependent phase errors due to δτ A can be eliminated completely. The remaining phase errors related to source position can be further removed by applying additional self-calibration solutions derived from the "ch0" data. However, the trade-off using the vectoraveraging technique for sub-bands is that the phase variation is effectively turned into signal noise, producing a coherence loss. For a visibility function V (ν) = A(ν)exp[iφ(ν)], the ratio of the vector-averaged amplitude to the actual amplitude drops with the amount of phase changes in phase across the sub-band (∆φ) due to a residual delay: . This ratio can be used for quantitative assessment of the coherence loss. The maximum loss in fringe intensity owing to a residual delay of 0.04 nsec is ∼0.06% of the core flux-density (F core ). The uncertainty of the vector-averaged amplitude is σ = A ∆φ 2 √ 6n , where n is the number of independent data channels in a given sub-band. Thus, the process of vector-averaging produces an uncertainty in the calibration modeling and generates additional noise in the data as well.
3. The method of vector-averaging sub-band data is subject to bandwidth smearing if a residual delay is present. The resultant fringes of vector-averaged sub-band data are modulated by a sinc function sin(π∆ν BW τ rs ) π∆ν BW τ rs , the delay beam or the bandwidth pattern (Thompson 1986). The residual delay here (τ rs ) is considered from two aspects: (1) The residual delay δτ A ; as discussed above, this delay term is referred to the time delay from unknown source when the target source is placed at the phase center that is presumed to be accurately tracked. (2) The geometrical delay, δτ pnt = B c ∆θ pnt , for a source offset by ∆θ pnt from the phase tracking center or the array pointing center. Hereafter, we refer to the geometrical delay, δτ pnt , as the source delay. Both delay terms, δτ A and δτ pnt , are time variable. Given the unclear origins, the behavior of δτ A appears unpredictable in time. The term δτ pnt is traceable. The phase compensation of the source delay at each of the spectral frequencies can be computed with a high precision interferometer model if no sub-band averaging is applied. The MFS algorithms have incorporated such corrections in the imaging process.
However, the shortcut with sub-band averaging will lead to a loss of valuable frequency-dependent information and therefore the image fidelity will be substantially degraded. For sub-band width of ∆ν BW = 128 MHz and the maximum baseline of the JVLA, at angular distance of 5" from the delay center or phase center, the fringe amplitude will drop to 75% of the value in the spectral data without sub-band averaging. In addition, the structure of the source will be distorted due to the bandwidthsmearing effect. At the radial distance of ∆θ pnt away from the phase center, a point source will be smeared to an angular extent of ∆θ pnt ∆ν BW ν (Thompson 1986). For example, at ∆θ pnt = 10" from the phase center, the angular distortion of a point source can be ∼0.14", more than half of the A-array synthesized beam, if sub-band averaging is applied to a 128 MHz sub-band data at 9 GHz.
4. The residual delay, δτ A , varies with time in an unpredictable way. The unclear origin of δτ A causes difficulties in eliminating such frequency-dependent phase errors from the data of target sources. The phase issue caused by δτ A can be corrected partially by applying the frequency-dependent phase corrections interpolated from the time-dependent bandpass solutions derived from a reference calibrator close to the target. A more sophisticated technique using modern computing resources is needed that can extract the requisite information from the data of the target source itself. Development of new technique for data calibration and deep imaging with wideband data appears to be imperative.
In addition to the phase errors caused by residual delays, the variation in visibility amplitude as a function of both time and frequency leads to detrimental effects in deep imaging with wideband data if this variability is not accounted for. The source activity as function of time and the distribution of the radiative energy as a function of frequency or source spectrum are not unusual among astronomical sources. The issues related to the amplitude variation in DR imaging with wideband data include the following: 1. Source spectrum. For a given spectral energy distribution, the visibility amplitude varies as a function of frequency across a broad frequency band. For a source having a steep power-law spectrum (∼ ν α ) with α = −1, the fractional amplitude changes by ∆S/S = ∆ν/ν ∼ 36% across the 2-GHz JVLA wideband at 5.5 GHz. Unlike the former VLA continuum data with narrow IF bandwidth, the amplitude variation across such a wide band challenges traditional imaging techniques. In order to cope with the spectral issues, a new algorithm for imaging the first two terms in the Taylor expansion of the visibility spectra has been developed and has been used in the clean program within CASA (Rau & Cornwell 2011). This technique has been successfully employed for deep imaging of the first term, the intensity distribution.
2. Time variation. Fig. A3 shows the issue in deep imaging owing to the time-variation of core flux-density. A variation of 10% in flux density of a 1-Jy core can reduce the dynamic range by a factor of 6-7, limiting deep imaging. In addition, artifacts around the core can be generated in the "cleaned" image if the variability during the observation is not removed from the data (Zhao et al. 1991). A procedure for removal of the time variation of a compact radio core has been described using the software tools in CASA (Zhao, Morris & Goss 2016).
The rest of the content in this appendix concerning the method of RE-correcting and DR-imaging is organized as follow. Section A.1 formulates a general model of the radiation structure and the variable components for given a source. Section A.2 describes a set of fomulas for transforming extended source structure into a point-source model. Section A.3. discusses the JVLA standard data reduction procedure. Based on the modeling discussed in these three sections, we develop a procedure for residual-error (RE) correction and dynamic-range (DR) imaging utilizing CASA software. Section A.4 outlines five major steps involved in the RE-correcting and DR-imaging process, with a demonstration using JVLA Ka-band data from the A-array observation of NRAO 530 at 33 GHz. Section A.5 introduces a case involving a more sophisticated procedure for handling the JVLA C-array data of NRAO 530 observed at 5.5 GHz. Section A.6 discusses the application of residual-delay corrections derived from gain calibrators to target sources. A summary is given in Section A.7. The techniques discussed here can be applied to general interferometer array data, although the suggested procedures are based on JVLA data and CASA algorithms.
A.1. Modeling structure and variability of a radio source The visibility function on the [u, v]-plane of a radio source sampled with a wide-band interferometer at observing time t and frequency ν can be written as: where τ A (t) and G A (ν,t) correspond to time delays and timedependent complex gains, respectively. Both terms are antennabased, containing instrumental errors and atmospheric effects. The source function, V S (u, v, ν,t), composed of numerous emission components, can be described as: is the source visibility function of the i-th extended component, and V C i (u, v, ν,t) is the visibility function of the i-th unresolved compact or point source. The integers m and n are assumed the total numbers of extended and unresolved components in the field of view, respectively. The extended components are presumably stable during the observing period while the flux-densities of the compact components can vary with time. In practice, for a QSO or a blazar in which a dominant radio core is present, the source function can be simplified as a strong point-source function V C (u, v, ν,t) and an arbitrary extended emission V Ext (u, v, ν), so that, (u, v, ν).
(A3) For the MOJAVE sources at 1.4 GHz, the ratio of the integrated flux density from the extended components to that of the core is distributed across a very wide range between a few thousands to less than a thousandth (Kharb et al. 2010). For the two gain calibrators NRAO 530 and J1744-3116 at 5.5 and 9 GHz, the integrated flux density of the extended portion is only a few percent that of the core. In core-dominant sources such as most blazars and QSOs at high frequencies, V S (u, v, ν,t) ∼ V C (u, v, ν,t) can be used in the initial approach.
As interferometer arrays improve in both sensitivity and angular resolution, the previous assumption of a calibrator (a QSO or a small planet) as a point source is no longer suitable. The extended structure and variability of a calibrator appear to generate significant effects, limiting both the calibration precision and the image quality as well as fidelity. Given a celestial object, a time-invariant structure of the source emission can be constructed, as was shown in the analysis of the radio observations and wideband imaging of the complex region Sgr A after fixing the variable component of Sgr A* to a constant value (Zhao, Morris & Goss 2016). Then, a model for the structure of stable emission from a calibrator or a radio source can be expressed as:

. Building a point source model
Dividing the observed complex visibility V (u, v, ν,t) of Eq(A1) by the source model V S (u, v, ν) of Eq(A4), one can approach to point-like visibility function: where the normalized source function V P (u, v, ν,t) is an initial calibration model in which the deviations from a point source are due to the atmospheric effects, such as amplitude attenuation and phase distortion, and the antenna-based instrumental errors in both amplitude and phase, as well as error due to imperfections of the initial source model. The process of calibrations for various effects will eventually converge the normalized source function to a point source function. For example, the linear changes in phase as a function of frequency across the sampling band can be fit by time delays and then be removed from the data. In addition, for QSOs or blazars, time variations of the flux densities from their radio cores often become significant. Special attention is needed in the calibration and imaging processes (Zhao et al. 1991;Zhao, Morris & Goss 2016). Thus, when a variable core is present, the source function can be separated into two components, a component, , the flux-density at the beginning of the observation, t 0 , and a time-variable term giving the deviations from that initial value, ∆V C (u, v, ν,t): For calibration concerns, the radio-dominant core is placed at the phase tracking center of the interferometer. Given a variable intensity I 0 + ∆I of the radiation from its core, the complex finction (A6) of the compact source can be simplified as a real function, V C (u, v, ν,t) = I 0 + ∆I(t). If the extended component V Ext (u, v, ν) is not time variable, then the normalized source function of Eq(A5) can be expressed as, where the spectrum of the core is further assumed to be flat for simplicity. Removing the time-variable term is necessary prior to further calibration of complex gains (Zhao, Morris & Goss 2016). For a calibrator with no time-variations and the source structure is simple enough, the normalized source function V P (u, v, ν,t) is essentially a complex function tracing the antenna gains multiplied by the antenna-based delay term, i.e., the right of Eq(A8) can be defined as the overall antenna-based complex gain G A (u, v, ν,t) ≡ e 2πiντA(t) G A (ν,t). Furthermore, the antenna-based delay can be further separated into two terms, a constant delay τ 0 and a time-dependent residual delay: For data from the JVLA, τ 0 can be removed by the standard delay correction in the beginning of data calibration or in a pipeline process, and δτ A (t) is the residual delay that appears to be dependent on time. Furthermore, the frequency-dependent and time-dependent functions in the complex gain can be decoupled if the instrumental bandpass is stable enough in time, i.e., G A (ν,t) = G A (t)B A (ν), the overall antenna-based gains becomes where G A (t) is the antenna-based complex gain, a function of time; with a scaling factor of the amplitude F , the antenna gain can be expressed as, and B A (ν) is the antenna-based, instrumental complex bandpass, a function of frequency, A.3. The standard data reduction procedure for JVLA data In principle, the calibrations of delay, bandpass, complex gain and flux scale with the JVLA standard data reduction procedure or in VLA pipeline 9 are essentially to correct for the 9 see https://science.nrao.edu/facilities/vla/data-processing/ parameter terms τ 0 , b A (ν), φ bA (ν), g A (t), φ gA (t) and F , respectively (Fig. A4, top box). After going through the standard data reduction procedure, various errors discussed above are corrected, neglecting the extended emission structure. Then, residual errors remained in G A (t) and B A (ν), or Eq(A11) and Eq(A12), can be assessed by the normalized gains, g A (t) = 1 + Left: Image made from a combination of the two data sets which have been self-calibrated with independently determined models, including corrections for residual delay and intra-day variability of the core, but with no correction for the variation from 0.71 Jy on Day 1 to 0.77 Jy on Day 2. A high noise level is present owing to the 10% variation in flux density of the radio core. The rms noise level in the region far from the core is σ far ∼45 µJy beam −1 . Middle: The image made from the same uncorrected data but adding a value of 0.06 Jy at the core to the data of Day 1, establishing a constant flux density of the core in both data sets. The rms noise of σ far is improved to the level of σ far ∼6 µJy beam −1 . Right: The image made with both data sets with further corrections for the residual errors using the improved model. The rms noise σ far is ∼2 µJy beam −1 . The contours are 2σ far ×(−2, 2 n ) with n = 1, 2, 3, until reaching the peak. The synthesized beam is 0.20"x0.39" (-31 • ).
gains are much smaller than the residual errors O [g A (t)] in the time-dependent gains and the phase term owing to residual delay dominates frequency-dependent residual errors in phase, ignoring O [b A (ν)] and φ in Eq(A13), the first order of the residual errors is reduced to: (A14) The remaining residual errors can be separated, namely the residual delay δτ A (t) and the residual gain O [g A (t)]. Both are time dependent. The residual gains can be corrected with selfcalibration, the well-known technique (Schwab 1980) that has been successfully applied in continuum imaging for radio interferometer arrays. The residual delay δτ A (t) can be fit to timedependent antenna-based phase solutions for the bandpass. Then, the residual delay can be essentially removed in the process of constructing a high-precision point-source model.

A.4. Corrections for residual delays
In processing broad bandwidth interferometer data with high sensitivity, both the residual delay and extended emission struction prevent from building a high precision point-source model. In this section, we discuss the procedure of RE-correcting and DR-imaging based on the JVLA data. Iterations of deliberately designed steps in processing wideband data appears to be a practical way on correction for the residual errors.
For the JVLA data at 9 GHz, the typical residual delays δτ A (t) produces a phase change at a level of 2 • across a sub-band of BW 128 MHz. The origin of the residual errors is not ascertained, but they likely arise from a combination of the instrumental issues related to the changes in the telescope backend electronics and the telescope frontend. The later includes pointing drifting owing to weather conditions such as a strong wind as well as the gravitational deformation of telescope beams as function of sky position as the telescope tracking a celestial source. The changes are dependent on the conditions of the individual telescope elements across the array. Based on the analysis discussed in A.1, A.2 and A.3, a procedure for REcorrecting and DR-imaging can be developed with the tasks and modules in the CASA software. Fig. A4 outlines the basic steps in corrections for the residual delays, residual gains and timevariable issues. The method for removing the time-variable portion of the core flux-density has been discussed (Zhao, Morris & Goss 2016). The critical issue of the procedure is to provide a way of approaching a high-precision source model with the steps of corrections for residual errors. We use the JVLA Kaband data of NRAO 530 as an example to illustrate the process. A description of the steps is outlined in the caption of Fig. A4.
The image created after applying JVLA standard calibration (step 1) is a point-like source with DR of 6,800:1 (top image in Fig. A4). With this model, the sub-band or spectral window (spw) based residual errors of antenna gains in phase (timedependent) are corrected with self-calibration. Then, alignment of the time-dependent offsets in both phase and amplitude between the sub-bands are carried out with the solution (time-dependent) derived from the CASA task bandpass by averaging all the channels in each of the sub-bands. Solutions for both gains and bandpass are derived from time intervals equal to each of the observing scans. The bandpass solutions are equivalent to the scan-based offsets in ampltiude and phase between the "ch0" data of each sub-band. Correcting for the residual phase errors from self-calibration with the task gaincal along with applying the bandpass solutions of spectral window alignment for each of the 64 sub-bands (step 2), a better image model is produced. With the updated model, further bandpass solutions are solved in a shorter time interval 30 sec by averaging every 4 channels in each of the sub-bands. The phase slopes across each of the sub-bands present in the bandpass solutions suggests the residual delays that have not been fully corrected in the data. Applying the new bandpass solutions to the data FIG. A4.-A flow chart illustrating the procedure of corrections and calibrations for the instrumental and atmospheric issues in the process of HDR imaging. The imaging models (on the right column) were produced from a 10-mininute test data of NRAO530 observed on 2015-9-11 using JVLA in the A-array at 33 GHz in the Ka-band with a total bandwidth of 8 GHz. The images are scaled in units of arcsec.
Step 1 is the JVLA standard calibrations, including the corrections for instrumental delay τ 0 and bandpass G A (ν), calibrations for complex gains G A (t) of interferometer data, and boostrap flux density scale F from a primary calibrator. Top-right image shows the model image produced from this step with a DR of 6,800:1.
Step 2 is further corrections for phase errors and alignment of the sub-band offsets in both phase and amplitude, utilizing the self-calibration technique and multiple-spw bandpass solutions, improving the imaging model.
Step 3 is to further remove the residual delays in each of sub-bands; an improved model image (Middle-right) with an improved image DR of 57,000:1 is produced from this step.
Step 4 is to correct for the residual gain errors in both phase and amplitude with solution intervals from self-calibration as short as possible. Step5 is to remove the time variation in flux density from the compact core, reaching the final imaging model if no significant residual closures present. Otherwise, corrections for closure errors are needed. The bottom-right image shows the final model image with image DR 340,000:1 limited by sensitivity.
(step 3), the improved image is constructed, showing an image DR ∼57,000:1 (the middle image in Fig. A4). The reference model is converging to the true source model. With the updated image model produced from the previous iteration, the corrections for the residual gains in both phase and amplitude can be solved with gaincal in the integration time, the shortest time interval of the data. Applying the new gain solution (step 4), the further corrected dataset is generated. Then, the variability of the core is inspected with a light curve showing the varibility index in flux-density 2 S max − S min S max + S min = 0.1% within the observing period of 2.7 hr. A model for the variable component is built on the scan-averaged basis with the value derived from the first time scan as the reference flux density. The offsets from the reference flux density at the following time scans are computed. This variability model is then subtracted from the data with the method discussed by Zhao, Morris & Goss (2016). After step 5, the antenna-based issues as well as fluxdensity variation have been corrected. A new model image can be constructed. In the region a few arcsec from the core, the rms noise in the new model image appears to be at the level of the thermal noise. One may go back to step 2 and redo the procedure if the image DR is not limited by telescope sensitivity. If the baseline-based residual closure errors 10 present, the baseline-based corrections may need to apply to the data. For the case of NRAO 530 Ka-band data, we redo the procedure from step 2 and apply baseline-based corrections with blcal.
The resulting image is consistent with the image made in the step 5. The rms noise of the final image made with robustness weight (R=0) is ∼10 µJy beam −1 , reaching the sensitivity limit of the JVLA. The image DR of the final image (the bottom-left image in Fig. A4) is 340,000:1.
A.5. A note to the 5.5-GHz C-array data The above procedure appears to be quickly converging for the high-resolution data sets used in this paper. Owing to large positional uncertainty generated from the residual delays in the shorter baselines or lower-resolution data as well as contamination by low-level RFI, the initial model is poorly determined from our JVLA C-array data taken at 5.5 GHz in the spring of 2012. The process of RE-correcting and DR-imaging is difficult to converge following the steps outlined in Fig. A4. Consequently, a more elaborate procedure has been developed based on the C-array data, adding sub-steps to complement the major correction steps of the procedure. Fig. A5 (top) shows a flow chart of the RE-correcting and DR-imaging procedure. The details appear to be necessary to converge the imaging process while fixing data issues using software in CASA. In the data correcting flow, eight progressive imaging models, marked with a symbol of capital M followed by a serial number (the filled light-blue circles in Fig. A5), are created. The imaging model is needed for two purposes: (1) to be utilized as an input imaging model to generate solutions from the data processed in previous step; (2) to make statistics with the output image as compared with the model in previous step to judge if the process converging; if failed in converging, adjusting the input parameters and redo the step; if converging, then comparing the data with the improved imaging model to further edit and flag the data to eliminate the lower-level RFI and the erroneous data owing to incorrigible telescope issues. Thus, the details of each step might be different for different data sets. For the C-array data at 5.5 GHz used in this paper, we summarize the eight sub-steps marked in the flow chart after the initial process that includes the JVLA standard calibrations and fixing core flux-density to a constant. The output data from the initial calibration are averaged to 15s in integration and the frequency configuration remains as the original, namely 64 spectral channels with a uniform width in each of 16 sub-bands. A brief description of the sub-steps in the process, identified with the serial number of the imaging model, is provided below: M1 -A model created from the data calibrated initially. Using this model, the data are edited; then using the edited data, solutions for bandpass correction are solved by averaging the 64 channels in the basis of sub-band with a time interval of 60 sec; the solutions are applied to data to align the sub-bands in both phase and amplitude; self-calibration is applied to the aligned data for residual gain corrections. M2 -A model produced from the data processed in step M1; the process in this step is similar to step M1 except for the band-pass solutions computed by averaging every 8 channels in each of the sub-bands. A shorter interval of 30 sec is used in the averaging. The bandpass solutions are applied to the data for correcting the residual errors in delay.
M3 -A new model created from the data post the steps M1 and M2. Using this model, the residual gain errors are further corrected with self-calibration in the interval of 15s. Then, light curves for each of the two epochs are produced to fix the residual variations in flux-density with the module composed of add.component, ft and uvsub in CASA.
M4 -This sub-step is similar to the sub-steps M1 and M2, except for shorter time interval (15s) used for solving bandpass solutions in a spectral bin of every two channels. The remaining residual delays are further corrected.
M5-M6-M7 -The three sub-steps to average 64 channels into single channel by vector averaging the data in each of subbands, producing "ch0" data in each of the sub-bands. Every four channels are averaged in each of the three sub-steps; consequently, the S/N of the new channel data in each of the sub-steps increases by a factor of 2, which helps to discern the lowerlevel RFI. Data editing is applied to reject the erroneous data that are hidden in the noise. Further self-calibration is carried out, ensuring the previous gain corrections to be not affected by bad data and lower S/N issues. Post M7, the output data becomes single channel (continuum) in each of the sub-bands. We note that one needs to be cautious for possible distortions of the source structure owing to the bandwidth smearing effect at a large distance from the phase center when applying the vector averaging across each of the sub-bands.
M8a new imaging model created from the data corrected for antenna-based error. This model is used to correct for baselinebased errors. The scan-averaged solutions are applied to the data. After flagging the data with large deviation from the baseline-based fringe curve, the output is the final data. Fig. A5 (bottom) shows a comparison between the image of M1 (left panel) made from the initially input data (red box) and the image (right panel) from the final output data (green box). The process of RE-correcting and DR-imaging is converged. The residual error correction appears to be the key process in the restoration of image from radio interferometer array data in addition to FFT and deconvolution or cleaning process.
A.6. Applying residual-delay corrections to target sources The time-dependant bandpass solutions for the residual delays derived from calibrators can be applied to target sources. In Morris, Zhao & Goss (2017), we demonstrated that the residual delay corrections derived from NRAO 530 at 5.5 GHz and J1744-3116 at 9 GHz are applied to Sgr A* producing HDR images. However, the resultant images of Sgr A* suggest that the corrected data at 9 GHz produces a better image than the one at 5.5 GHz. We note that both NRAO 530 and J1744-3116 were frequently scheduled as a gain calibrator for Sgr A* in the observations at 5.5 GHz and 9 GHz, respectively. However, J1744-3116 is 2.5 • away from Sgr A* while NRAO 530 is 16 • away. The difference between 9 GHz and 5.5 GHz images might suggest that the accuracy in the correction for the residual delay is sensitive to the amount of positional offsets between target and calibrator.

A.7. Summary
In short, implementing the algorithms and methods for processing wideband data in addition to the fundamental technology developed during the former, narrow bandwidth VLA operation, the CASA has successfully provided users a software platform for handling JVLA data. With the procedures as discussed in this Appendix, we have demonstrated that, within CASA, it is possible to achieve HDR images, reaching the telescope sensitivity limits. However, the approach is tedious and time consuming. Various input parameters for the CASA programs used in the procedures require specifically adjusting in the RE-correcting and DR-imaging process. With the modern computing power available, it becomes possible to have a better way using more advanced algorithms to cope with these issues that we have confronted in the operation of wideband telescopes and handling big data. Therefore, improvements in programming are needed to arrange and manipulate the computer units in a more efficient way for processing wideband data in a combination of the advanced statistical algorithms with the knowledge that radio astronomers have developed over the past decades.