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

TL;DR

This paper presents a new wideband radio imaging procedure that corrects residual interferometer errors, enabling high-fidelity images of distant blazars like NRAO 530 with unprecedented detail, revealing complex jet structures.

Contribution

The paper introduces a novel data reduction procedure for wideband radio observations that significantly improves image dynamic range and fidelity, applicable to broadband continuum data.

Findings

Achieved dynamic range exceeding 1,000,000:1 in radio images.

Revealed detailed jet bending and structure in NRAO 530.

Identified brightness asymmetry suggesting mildly relativistic jets.

Abstract

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.

Figures

Tables

Table 1. Table 1: Log of dataset

Project (PI)	Observing date	Array	Band†	$\nu$	$\Delta \nu$	HA range	Time	Sources
				(GHz)	(GHz)		on source (h)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)
12A-037 (Zhao)	2012-03-29	C	C‡	5.5	2	$-3^h.0$ — $+0^h.8$	0.94	J1733$-$1304
12A-037 (Zhao)	2012-04-22	C	C‡	5.5	2	$-0^h.5$ — $+2^h.9$	0.82	J1733$-$1304
12A-037 (Zhao)	2012-07-24	B	C‡	5.5	2	$-0^h.9$ — $+3^h.0$	0.63	J1733$-$1304
12A-037 (Zhao)	2012-07-27	B	C‡	5.5	2	$+0^h.4$ — $+3^h.9$	0.71	J1733$-$1304
14A-346 (Morris)	2014-05-17	A	C‡	5.5	2	$-3^h.2$ — $+0^h.7$	0.31	J1733$-$1304
14A-346 (Morris)	2014-05-26	A	C‡	5.5	2	$-0^h.4$ — $+3^h.5$	0.47	J1733$-$1304
14A-231 (Yusef-Zadeh)	2014-03-01	A	X‡	9.0	2	$-3^h.0$ — $+3^h.2$	0.27	J1733$-$1304
						$-3^h.0$ — $+3^h.5$	1.18	J1744$-$3116
14A-232 (Yusef-Zadeh)	2014-04-17	A	X‡	9.0	2	$-3^h.3$ — $+2^h.0$	0.19	J1733$-$1304
						$-3^h.4$ — $+3^h.4$	1.12	J1744$-$3116
SF-0853 (Haggard)	2014-04-28	A	X‡	9.0	2	$-3^h.0$ — $+3^h.2$	0.26	J1733$-$1304
						$-3^h.0$ — $+3^h.5$	1.05	J1744$-$3116
15A-293 (Yusef-Zadeh)	2015-09-11	A	Ka♯	33	8	$-2^h.6$ — $+0^h.0$	0.17	J1733$-$1304

Table 2. Table 2: Compact components in NRAO 530

ID	$S_p\pm \sigma$	$T_B$	$\nu$	$\alpha$	$S_t\pm \sigma$	$\Delta \theta \pm \sigma$	PAθ	${\Theta }_{\mathrm{maj}}\pm \sigma$	${\Theta }_{\min }\pm \sigma$	PAΘ
	(mJy beam-1)	(K)	(GHz)	(${\nu }^{\alpha }$)	(mJy)	(arcsec)	(deg)	(arcsec)	(arcsec)	(deg)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)
Core	$\frac{3400\pm 1}{\dots }$	$\frac{4.6\times {10}^5}{\dots }$	$\frac{33}{\dots }$	…	$\frac{3422\pm 1}{\dots }$	$\frac{0.000}{\dots }$	$\frac{0.0}{\dots }$	$\frac{0.011\pm 0.001}{\dots }$	$\frac{0.003\pm 0.001}{\dots }$	$\frac{176\pm 2}{\dots }$
	$\frac{4915\pm 2}{4715\pm 2}$	$\frac{1.4\times {10}^6}{3.7\times {10}^6}$	$\frac{9}{5.5}$	$\frac{0.08\pm 0.01}{\dots }$	$\frac{4922\pm 2}{\dots }$	$\frac{0.000}{\dots }$	$\frac{0.0}{\dots }$	$\frac{0.016\pm 0.001}{\dots }$	$\frac{0.007\pm 0.001}{\dots }$	$\frac{174\pm 4}{\dots }$
W1	$\frac{0.23\pm 3}{\dots }$	$\frac{31}{\dots }$	$\frac{33}{\dots }$	…	$\frac{0.57\pm 0.11}{\dots }$	$\frac{0.31\pm 0.02}{\dots }$	$\frac{-27\pm 2}{\dots }$	$\frac{0.21\pm 0.05}{\dots }$	$\frac{0.05\pm 0.03}{\dots }$	$\frac{146\pm 7}{\dots }$
W2	$\frac{0.31\pm 0.02}{\dots }$	$\frac{42}{\dots }$	$\frac{33}{\dots }$	…	$\frac{1.95\pm 0.12}{\dots }$	$\frac{0.69\pm 0.01}{\dots }$	$\frac{-42\pm 1}{\dots }$	$\frac{0.31\pm 0.02}{\dots }$	$\frac{0.14\pm 0.01}{\dots }$	$\frac{111\pm 3}{\dots }$
	$\frac{3.07\pm 0.24}{4.70\pm 1.70}$	$\frac{904}{3700}$	$\frac{9}{5.5}$	$\frac{-0.86\pm 0.18}{\dots }$	$\frac{8.38\pm 0.89}{\dots }$	$\frac{0.69\pm 0.01}{\dots }$	$\frac{-42\pm 1}{\dots }$	$\frac{0.49\pm 0.05}{\dots }$	$\frac{0.30\pm 0.02}{\dots }$	$\frac{125\pm 6}{\dots }$
W3	$\frac{0.022\pm 0.005}{\dots }$	$\frac{3}{\dots }$	$\frac{33}{\dots }$	…	$\frac{0.18\pm 0.05}{\dots }$	$\frac{1.23\pm 0.04}{\dots }$	$\frac{-52\pm 1.3}{\dots }$	$\frac{0.35\pm 0.10}{\dots }$	$\frac{0.16\pm 0.05}{\dots }$	$\frac{128\pm 18}{\dots }$
W4	$\frac{0.12\pm 0.01}{0.30\pm 0.02}$	$\frac{35}{240}$	$\frac{9}{5.5}$	$\frac{-1.8\pm 0.17}{\dots }$	$\frac{0.71\pm 0.03}{\dots }$	$\frac{2.35\pm 0.01}{\dots }$	$\frac{-58\pm 0.3}{\dots }$	$\frac{0.65\pm 0.03}{\dots }$	$\frac{0.49\pm 0.02}{\dots }$	$\frac{126\pm 5}{\dots }$
W5	$\frac{0.29\pm 0.01}{0.46\pm 0.02}$	$\frac{85}{360}$	$\frac{9}{5.5}$	$\frac{-0.94\pm 0.11}{\dots }$	$\frac{1.52\pm 0.08}{\dots }$	$\frac{3.74\pm 0.01}{\dots }$	$\frac{-63\pm 0.2}{\dots }$	$\frac{0.64\pm 0.03}{\dots }$	$\frac{0.45\pm 0.02}{\dots }$	$\frac{106\pm 4}{\dots }$
W6	$\frac{0.21\pm 0.01}{0.32\pm 0.01}$	$\frac{62}{250}$	$\frac{9}{5.5}$	$\frac{-0.86\pm 0.12}{\dots }$	$\frac{1.20\pm 0.06}{\dots }$	$\frac{4.28\pm 0.01}{\dots }$	$\frac{-64\pm 0.2}{\dots }$	$\frac{0.74\pm 0.03}{\dots }$	$\frac{0.41\pm 0.02}{\dots }$	$\frac{112\pm 2}{\dots }$
W7	$\frac{0.08\pm 0.01}{0.13\pm 0.01}$	$\frac{24}{100}$	$\frac{9}{5.5}$	$\frac{-0.99\pm 0.30}{\dots }$	$\frac{1.20\pm 0.09}{\dots }$	$\frac{6.06\pm 0.03}{\dots }$	$\frac{-69\pm 0.2}{\dots }$	$\frac{1.17\pm 0.08}{\dots }$	$\frac{0.75\pm 0.05}{\dots }$	$\frac{80\pm 6}{\dots }$
W8	$\frac{0.17\pm 0.01}{0.32\pm 0.02}$	$\frac{50}{250}$	$\frac{9}{5.5}$	$\frac{-1.3\pm 0.17}{\dots }$	$\frac{0.87\pm 0.06}{\dots }$	$\frac{7.61\pm 0.01}{\dots }$	$\frac{-73\pm 0.1}{\dots }$	$\frac{0.55\pm 0.03}{\dots }$	$\frac{0.49\pm 0.03}{\dots }$	$\frac{98\pm 22}{\dots }$
W9	$\frac{0.10\pm 0.01}{0.25\pm 0.01}$	$\frac{29}{200}$	$\frac{9}{5.5}$	$\frac{-1.8\pm 0.21}{\dots }$	$\frac{1.74\pm 0.09}{\dots }$	$\frac{8.86\pm 0.03}{\dots }$	$\frac{-76\pm 0.3}{\dots }$	$\frac{1.34\pm 0.07}{\dots }$	$\frac{0.81\pm 0.04}{\dots }$	$\frac{70\pm 4}{\dots }$
W10			$\frac{9}{5.5}$	$\frac{>-1.9}{\dots }$	$\frac{0.67\pm 0.06}{\dots }$	$\frac{9.40\pm 0.04}{\dots }$	$\frac{-79\pm 0.3}{\dots }$	$\frac{0.74\pm 0.06}{\dots }$	$\frac{0.59\pm 0.05}{\dots }$	$\frac{24\pm 14}{\dots }$
W11	$\frac{0.06\pm 0.01}{0.13\pm 0.01}$	$\frac{18}{100}$	$\frac{9}{5.5}$	$\frac{-1.6\pm 0.37}{\dots }$	$\frac{0.62\pm 0.06}{\dots }$	$\frac{10.03\pm 0.04}{\dots }$	$\frac{-84\pm 0.2}{\dots }$	$\frac{1.07\pm 0.11}{\dots }$	$\frac{0.55\pm 0.05}{\dots }$	$\frac{56\pm 5}{\dots }$
W12	$\frac{5.00\pm 0.16}{7.72\pm 0.32}$	$\frac{1500}{6100}$	$\frac{9}{5.5}$	$\frac{-0.88\pm 0.11}{\dots }$	$\frac{38.1\pm 1.3}{\dots }$	$\frac{11.20\pm 0.01}{\dots }$	$\frac{-87\pm 0.1}{\dots }$	$\frac{0.76\pm 0.03}{\dots }$	$\frac{0.54\pm 0.02}{\dots }$	$\frac{63\pm 4}{\dots }$
W13	$\frac{0.05\pm 0.01}{0.08\pm 0.02}$	$\frac{15}{63}$	$\frac{9}{5.5}$	$\frac{-0.96\pm 0.65}{\dots }$	$\frac{1.35\pm 0.06}{\dots }$	$\frac{11.58\pm 0.02}{\dots }$	$\frac{-97\pm 0.1}{\dots }$	$\frac{1.32\pm 0.06}{\dots }$	$\frac{1.09\pm 0.05}{\dots }$	$\frac{70\pm 9}{\dots }$
W14†	$\dots \frac{}{0.124\pm 0.002}$	$\dots$	$\frac{9}{5.5}$	$\dots$	$\dots$	$\dots$	$\dots$	$\dots$	$\dots$	$\dots$
E1	$\frac{0.12\pm 0.02}{\dots }$	$\frac{35}{\dots }$	$\frac{9}{5.5}$	…	$\frac{0.20\pm 0.04}{\dots }$	$\frac{0.62\pm 0.02}{\dots }$	$\frac{97\pm 2}{\dots }$	$\frac{0.37\pm 0.05}{\dots }$	$\frac{0.24\pm 0.03}{\dots }$	$\frac{148\pm 10}{\dots }$
E2†	$\dots \frac{}{0.10\pm 0.4}$	$\dots$	$\frac{9}{5.5}$	…	$\dots$	$\dots$	$\dots$	$\dots$	$\dots$	$\dots$
E3	$\frac{0.02\pm 0.01}{0.03\pm 0.01}$	$\frac{6}{24}$	$\frac{9}{5.5}$	…	$\frac{0.62\pm 0.06}{\dots }$	$\frac{6.5\pm 0.1}{\dots }$	$\frac{89\pm 2}{\dots }$	$\frac{1.27\pm 0.11}{\dots }$	$\frac{1.11\pm 0.09}{\dots }$	$\frac{89\pm 26}{\dots }$
E4†	$\dots$	$\dots$	$\frac{9}{5.5}$	…	$\dots$	$\dots$	$\dots$	$\dots$	$\dots$	$\dots$
E5	$\frac{0.03\pm 0.01}{0.06\pm 0.01}$	$\frac{9}{47}$	$\frac{9}{5.5}$	$\frac{-1.4\pm 0.7}{\dots }$	$\frac{5.4\pm 0.2}{\dots }$	$\frac{13.2\pm 0.1}{\dots }$	$\frac{73\pm 1}{\dots }$	$\frac{3.4\pm 0.2}{\dots }$	$\frac{2.7\pm 0.2}{\dots }$	$\frac{173\pm 6}{\dots }$

Table 3. Table 3: Extended components in NRAO 530

ID	St	Size	$\frac{S_{\mathrm{core}}}{S_t}$	$P_{\nu }$
	(mJy)	(arcsec)		(${10}^{27}$ W Hz-1)
(1)	(2)	(3)	(4)	(5)
5.5 GHz†
Core	4450$\pm$5	$\dots$	1	11
W-lobe	96$\pm$10	19	46	0.36
E-lobe	42$\pm$5	18	106	0.16
N-extension	7$\pm$1	15	636	0.03
S-extension	3$\pm$1	13	1480	0.01
Extended	148$\pm$13	37	30	0.56
Total	4598$\pm 14$	…	0.968	11.6
1.46 GHz‡
Core	6130$\pm$6	$\dots$	1	23
Extended	518$\pm$46	37	12	2
Total	6648$\pm$47	…	0.922	25
Sp. index♯	Core♠		Extended
${\alpha }_{1.46/5.5}$	$-0.23\pm 0.12$		$-0.94\pm 0.13$

Table 4. Table 4: Total flux density of NRAO 530

$\nu$ (GHz)	Sν (Jy)	For fitting	Duration	Ref.		$\nu$ (GHz)	Sν (Jy)	For fitting	Duration	Ref.
(1)	(2)	(3)	(4)	(5)		(6)	(7)	(8)	(9)	(10)
0.076	9.93$\pm$ 0.20	Yes	…	[18]		0.080	7.00	Yes	…	[1]
0.084	9.49$\pm$0.15	Yes	…	[18]		0.092	9.52$\pm$ 0.14	Yes	…	[18]
0.099	9.09$\pm$ 0.12	Yes	…	[18]		0.107	8.71$\pm$ 0.12	Yes	…	[18]
0.115	8.68$\pm$ 0.09	Yes	…	[18]		0.122	8.60$\pm$ 0.08	Yes	…	[18]
0.130	8.23$\pm$ 0.07	Yes	…	[18]		0.143	8.09$\pm$ 0.07	Yes	…	[18]
0.150	7.65$\pm$ 0.77	Yes	…	[17]		0.151	7.96$\pm$ 0.06	Yes	…	[18]
0.158	7.65$\pm$ 0.05	Yes	…	[18]		0.166	7.61$\pm$ 0.05	Yes	…	[18]
0.174	7.92$\pm$ 0.06	Yes	…	[18]		0.181	7.63$\pm$ 0.06	Yes	…	[18]
0.189	7.61$\pm$ 0.07	Yes	…	[18]		0.197	7.66$\pm$ 0.08	Yes	…	[18]
0.200	7.49$\pm$ 0.04	Yes	…	[18]		0.204	6.83$\pm$ 0.14	Yes	…	[18]
0.212	6.72$\pm$ 0.09	Yes	…	[18]		0.220	6.93$\pm$ 0.15	Yes	…	[18]
0.227	7.26$\pm$ 0.18	Yes	…	[18]		0.327	7.42$\pm$ 0.34	Yes	…	[19]
0.365	8.44$\pm$ 0.35	Yes	…	[22]		0.408	6.58	Yes	…	[1]
0.750	6.73$\pm$0.21	Yes	…	[2]		1.400	5.99$\pm$0.14	…	…	[2]
1.41	5.2	Yes	…	[1]		1.46	6.648$\pm$0.047	…	…	[13][15]
1.5	5.2	Yes	…	[16]		2.64	3.79 $-$ 5.41	Yes	8yr	[20][21]
2.7	4.3	Yes	…	[1]		4.8	4 $-$9	Yes	46yr	[9]
4.85	6.991$\pm 0.099$	…	…	[3]		5.0	4.1	Yes	…	[1]
5.0	5.0	…	…	[16]		5.5	4.74$\pm 0.13$	…	…	[15]
8.0	3.5 $-$ 12.5	Yes	46yr	[9]		8.1	10.5	…	…	[16]
8.4	6.2	…	…	[1]		9.0	5.002$\pm$0.050	…	…	[15]
14.5	2.5 $-$ 15	Yes	46yr	[9]		15	11.0	…	…	[16]
22	8.22 $-$ 11.78	…	0.3yr	[10]		33	3.426$\pm$0.017	…	…	[15]
37	2.25 $-$ 15.55	Yes	25yr	[4]		37	10.51 $-$ 12.20	…	0.1yr	[10]
43	8.98$\pm$1.87	…	…	[16]		43	16.6$\pm$0.5	…	…	[6]
86	11.2$\pm$0.3	…	…	[7]		86	2.10 $-$ 6.11	Yes	8yr	[21]
88	7.8$\pm$0.5	…	…	[8]		90	3 $-$ 15	…	11yr	[9]
95	12.6$\pm$0.7	…	…	[6]		107	13.2$\pm$0.4	…	…	[8]
150	7.09$\pm$0.35	…	…	[11]		215	6.2$\pm$1.1	…	…	[7]
230	2.4 $-$ 13.0	…	1yr	[12]		230	0.857 $-$ 3.868	Yes	16yr	[14]
270	5.66$\pm 0.30$	…	…	[11]		345	0.604 $-$ 2.41	Yes	15yr	[14]
350	0.73 $-$ 4.08	Yes	8yr	[5]		375	4.67$\pm$0.47	…	…	[11]

Equations

$${\mathrm{D}}_{\mathrm{Jet}}^i=\sum _{k=1}^{k=i}\left|{\vec{J}}_{k+1}-{\vec{J}}_k\right|,$$

$${\mathrm{PA}}_{\mathrm{Jet}}=-40.4°{\mathrm{Log}}_{10}\left(\frac{{\mathrm{D}}_{\mathrm{Jet}}}{\mathrm{arcsec}}\right)-52.5°,$$

$$r=a{e}^{\theta /b},$$

$$S_{\mathrm{c}}=S_0{\left(\frac{\nu }{{\nu }_0}\right)}^{5/2}\left(1-\exp \left[-{\left(\frac{\nu }{{\nu }_0}\right)}^{{\alpha }_{\mathrm{c}}-5/2}\right]\right),$$

$$S_{\mathrm{ext}}=S_1{\nu }^{{\alpha }_{\mathrm{e}}}.$$

$$V(\mathrm{u},\mathrm{v},\nu ,t)=V^{\mathrm{S}}(\mathrm{u},\mathrm{v},\nu ,t)e^{2\pi i\nu {\tau }_{\mathrm{A}}(\mathrm{t})}{G}_A(\nu ,t),$$

$$V^{\mathrm{S}}(\mathrm{u},\mathrm{v},\nu ,t)=\sum _{i=0}^{m-1}{V}_i^{\mathrm{Ext}}(\mathrm{u},\mathrm{v},\nu )+\sum _{i=0}^{n-1}{V}_i^{\mathrm{C}}(\mathrm{u},\mathrm{v},\nu ,t)$$

$$V^{\mathrm{S}}(\mathrm{u},\mathrm{v},\nu ,t)=V^{\mathrm{C}}(\mathrm{u},\mathrm{v},\nu ,t)+V^{\mathrm{Ext}}(\mathrm{u},\mathrm{v},\nu ).$$

$$V^{\mathrm{S}}(\mathrm{u},\mathrm{v},\nu )=V^{\mathrm{C}}(\mathrm{u},\mathrm{v},\nu )+V^{\mathrm{Ext}}(\mathrm{u},\mathrm{v},\nu )$$

$$V^{\mathrm{P}}(\mathrm{u},\mathrm{v},\nu ,t)=V(\mathrm{u},\mathrm{v},\nu ,t)/V^{\mathrm{S}}(\mathrm{u},\mathrm{v},\nu ),$$

$$V^{\mathrm{C}}(\mathrm{u},\mathrm{v},\nu ,t)=V_0^C(\mathrm{u},\mathrm{v},\nu )+\Delta V^C(\mathrm{u},\mathrm{v},\nu ,t).$$

$$V^{\mathrm{P}}(\mathrm{u},\mathrm{v},\nu ,t)=\left[1+\frac{\Delta I(t)}{I_0+V^{\mathrm{Ext}}(\mathrm{u},\mathrm{v},\nu )}\right]e^{2\pi i\nu {\tau }_{\mathrm{A}}}{G}_A(\nu ,t).$$

$$V^{\mathrm{P}}(\mathrm{u},\mathrm{v},\nu ,t)=e^{2\pi i\nu {\tau }_{\mathrm{A}}(\mathrm{t})}{G}_{\mathrm{A}}(\nu ,t).$$

$${\tau }_{\mathrm{A}}(t)={\tau }_0+\delta {\tau }_{\mathrm{A}}(t)$$

$${\color[rgb]{0,0,0}\mathscr{G}_{\rm A}(\nu,t)=G_{\rm A}(t)B_{\rm A}(\nu)e^{2\pi i{\rm\nu\delta\tau_{\rm A}}(t)},}$$

$$G_{\mathrm{A}}(t)=\mathcal{F}{g}_{\mathrm{A}}(t)e^{i{\varphi }_{g_{\mathrm{A}}}(t)},$$

$$B_{\mathrm{A}}(\nu )=b_{\mathrm{A}}(\nu )e^{i{\varphi }_{b_{\mathrm{A}}(\nu )}}.$$

$$V^{\mathrm{P}}(\mathrm{u},\mathrm{v},\nu ,t)\approx e^{2\pi i\nu \delta {\tau }_{\mathrm{A}}(\mathrm{t})+i{ϵ}_{\varphi }}\left(1+O\left[g_{\mathrm{A}}(t)\right]+O\left[b_{\mathrm{A}}(\nu )\right]\right),$$

$$V^{\mathrm{P}}(u,v,\nu ,t)\approx e^{2\pi i\nu \delta {\tau }_{\mathrm{A}}(\mathrm{t})}\left(1+O\left[g_{\mathrm{A}}(t)\right]\right).$$

Full text

ISM: individual objects (NRAO 530) — BL Lacertae objects: general — galaxies: active — galaxies: jets — radio continuum: galaxies — techniques: image processing — techniques: interferometric

††slugcomment: Draft:

1 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 residual 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$\sim$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_{\rm L}=5.8$ Gpc (1″= 7.8 kpc)111Calculated using the cosmology calculator of Ned Wright, UCLA http://www.astro.ucla.edu/ wright/CosmoCalc.html assuming $H_{0}=71$ km s*-1* Mpc*-1*, $\Omega_{\rm M}=0.27$, $\Omega_{\rm vac}=0.73$ (Komatsu et al., 2009).. 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., 2016, 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 MOJAVE222Monitoring Of Jets Active galactic nuclei with VLBA Experiments: http://www.physics.purdue.edu/MOJAVE/. 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 $\sim 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.

2 Data sets and imaging

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).

2.1 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.

2.2 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, 2013, 2016; Morris, Zhao & Goss, 2017). NRAO 530 was used as a gain calibrator. Using CASA333the Common Astronomy Software Applications package, is being developed with the primary goal of supporting the data post-processing needs of the next generation of radio astronomical telescopes such as ALMA and VLA. https://casa.nrao.edu/, we followed the standard calibration for a continuum source. In order to achieve a high dynamic range (HDR) in a high-resolution image (A-array), in addition to the VLA standard calibrations, two further corrections must be taken: 1) correction for residual delays444Small phase errors, present as slopes across sub-bands with a typical value of $\sim$2° per 128 MHz, remain in the data after applying JVLA standard calibration. The slopes appear to be caused by residual delays with typical values of $\sim$0.04 nsec. The residual delays are antenna-based and vary with time but with unpredictable trend. 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 $\mu$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 by combining the A- and B-array data. An rms noise of 3.5 $\mu$Jy beam*-1* is achieved in the region without source emission ($>$15” from the image center) while the rms noise is 4.5 $\mu$Jy beam*-1* near the core ($\sim 4"$)555Note: the regions used for statistics may contain faint extended emission near the core., comparable to the expected thermal noise of $\sim$3 $\mu$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 flow-chart (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 ${\rm 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 $\mu$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 C-array 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 applied to the combined data with OUTERTAPER 666OUTERTAPER is an input parameter of the CASA task clean, which sets the boundary for long baselines of the sampled data. $=300$ k$\lambda$ 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”$\times$0.69” (4°). The rms noise is $\sigma=5\mu$Jy beam*-1*. The top-right panel in Fig. 2 shows the central 6”$\times$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$\sigma$. 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).

2.2.2 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 Appendix 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 $\mu$Jy beam*-1*, giving DR of 1,300,000:1 with the synthesized beam of 0.27”$\times$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\sigma$ (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 $\mu$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 hour-angle 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 UV coverage in the three observations.

2.2.3 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$\mu$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.

3 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 well-confined 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 E- and W-lobes, respectively. At 5.5 GHz, the specific power of the core is about 1.1$\times 10^{28}$ Watt Hz*-1*, dominating the power of $\sim 5.6\times 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).

3.1.1 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_{\rm 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”$\times$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$\sigma$ limit ($S_{\rm p}=24~{}\mu$Jy beam*-1*) of the A-array image at 33 GHz is consistent with E1 being a newly 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.

3.1.2 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_{\rm p}=$ 30 and 60 $\mu$Jy beam*-1* as determined from the images convolved to a commom beam (0.27”$\times$0.19”, PA=$-17$°). The intensity contrast to the core is $S_{\rm p}^{\rm E5}$: $S_{\rm p}^{\rm Core}$ = 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\pm 0.7$.

3.1.3 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_{\rm 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 $\alpha=-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$\sigma$. The rest of the components in the western jet appear to be resolved out and not detected at 33 GHz.

3.1.4 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).

3.1.5 North-south extension to the core

In the C-array image with low angular resolution of $\sim$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-, 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$\pm$1 mJy, about 5 percent of the western lobe, including the jet, and the flux density of the southern extension is 3$\pm$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 $\sim$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 $\vec{J}$ – the vector of jet axis

Fig. 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°$\pm$7° for W1 at an angular offset $\Delta\theta=0.31$” from the core and PA=111°$\pm$3° for W2 at $\Delta\theta=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 $\vec{J}$ as the vector that connects the core to a jet component. The vector $\vec{J}$ can be quantitatively described by its position angle (PAθ) and its angular offset from the core ($\Delta\theta$). 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 $\Delta\theta$ increases. Table 2 lists $\Delta\theta$ and PAθ in columns 7 and 8, respectively, for all components.

Fig. 7 (top) plots the position angle PAθ as function of $\Delta\theta$. The western (eastern) jet components are represented in green (red); the quantities PAθ and $\Delta\theta$ 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 $\vec{J}$ is essentially to the north on the scale of VLBA component ($<50$ mas) and gradually changes to west as $\Delta\theta$ increase, the angular distance. The transition where PAθ changes rapidly occurs in the range between 0.05” and 1” in $\Delta\theta$. The northern VBLI jet and the western JVLA jet belong to one astrophysical entity, hereafter referred as the northwestern jet. (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 $\sim 90$° over the scale of $\sim$10” (or $\sim$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.

3.2.2 $\vec{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 $\vec{T~{}}$ defined as the direction change of the vector $\vec{J~{}}$ between consecutive jet components: at the $i$-th jet component $\vec{T_{i}}=\displaystyle{\vec{J}_{i+1}-\vec{J_{i}}\over|\vec{J}_{i+1}-\vec{J_{i}}|}$, where $\vec{J~{}}$ is the vector pointing from the core to the jet component, defined in the previous section. The direction of $\vec{T_{i}}$ can be described by PAJet, the position angle of $\vec{T_{i}}$. For a curved jet, PAJet is a function of D${}_{\rm Jet}^{i}$, the distance measured along the projected jet trajectory from the core to the location of the $i$-th component. The quantity D${}_{\rm Jet}^{i}$ 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:

[TABLE]

where, $i=1,2,...,n$ (see the caption of Fig. 7 for more details). Fig. 7 (bottom) plots PAJet as a function of DJet with a logarithmic scale in the DJet axis to facilitate visualization of the PAJet changes in the transition from the northern jet on the VLBI scale to the western jet on the VLA scale. Within a scale of 25 mas, PAJet varies between $-6^{\rm o}$ and 59o based on the six brighter VLBI components d, e, f, g, i and j. The function PAJet(DJet) 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 PAJet 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${}_{\rm Jet}=$7.8” (PA${}_{\rm 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 PAJet, 30 deg arcsec*-1*. The rapid change in PAJet 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; Hardee & Norman, 1988; Norman & Hardee, 1988; Zhao et al., 1992a, b, for example).

In short, from the diagram of PA${}_{\rm Jet}-$DJet, the change in PAJet seems to scale with jet kinematics, i.e. larger change in PAJet occurs in the inner region of the core, where the helical motion of the jet material dominates the PAJet change while the PAJet change appear to be caused by a slow precession of jet axis over the VLA scales. On the VLBI scale, the change rate in PAJet drops from 83 deg mas*-1* at D${}_{\rm Jet}\sim$1 mas to 2 deg mas*-1* at 25 mas. The rapid change in PAJet 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 PAJet changes about 90°. We made a linear fit to the plotted data and found that,

[TABLE]

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 (${\rm PA}_{\rm jet}$, ${\rm D}_{\rm Jet}$) is substituted by ($\theta$, $r$):

[TABLE]

where the variables $\theta$ and $r$ are in the units of degree and arcsec and the coefficients $a=0.0502^{\prime\prime}$ and $b=-17.5\degree$.

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.

3.3 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; An et al., 2013, and the F-GAMMA monitoring program777https://www3.mpifr-bonn.mpg.de/div/vlbi/fgamma/fgamma.html). 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$\pm$0.97c and median 13.74$\pm$0.40c based on the inner six components by the MOJAVE project888https://www.physics.purdue.edu/MOJAVE/sourcepages/1730-130.shtml (Lister et al., 2013, 2016, 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 $\nu>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,

[TABLE]

where $\nu_{0}$ is the turnover frequency at which the synchrotron optical depth becomes unity, $S_{0}$ is the flux density at $\nu_{0}$ multiplied by a factor of $\displaystyle{e\over e-1}=1.58$, and $\alpha_{\rm 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 $\alpha_{\rm e}$ and the flux density of the extended emission at 1 GHz of $S_{1}$ for $\nu$ in units of GHz,

[TABLE]

There are five parameters ($S_{0}$, $\nu_{0}$, $\alpha_{\rm c}$, $S_{1}$ and $\alpha_{e}$) to fit. The spectral index $\alpha_{e}$ can be determined from the extended emission. With $\alpha_{\rm 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 $\nu_{0}=0.16$ GHz, $S_{0}=6.2$ Jy or $S_{\rm c}=3.9$ Jy at $\nu=\nu_{0}$, and $\alpha_{\rm 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 $\alpha_{\rm 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 $\nu=142$ MHz, the contributions in flux density from the core and extended emission are comparable. At $\nu<$ 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 $\nu>>$ 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\times 10^{28}$ and $4.5\times 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 at high frequencies is relativistically Doppler boosted (e.g. Lu, Krichbaum & Zensus, 2011).

4 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 edge-brightened 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_{\rm 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_{\rm radio}\propto~{}L_{\rm opt}^{1.8}$. Based on their nearby sample ($z<0.5$), the FR I/II division appears to occur over a range in $P_{\rm 1.4~{}GHz}\sim 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_{\rm R}$ of the host galaxy to the central black hole mass $M_{\rm BH}$ (McLure & Dunlop, 2001) and the radio power at 1.4 GHz, $P_{\rm 1.4GHz}$, to the kinetic luminosity of the jet, $L_{\rm 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, $\displaystyle{L_{\rm jet}/M_{\rm 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).

4.1 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 flux-density-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_{\rm 1.4GHz}=1.0\times 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_{\rm ext,1.4GHz}=2.1\times 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., 2013, 2016), becoming highly variable. Taking the median value determined from the inner core components with the VLBA, the apparent transverse velocity va in units of $c$ (the speed of light) is $\beta_{\rm a}=14$, giving a lower limit on the Lorentz factor of $\gamma_{\rm min}=\sqrt{1+\beta_{\rm a}^{2}}\approx 14$ , assuming the jet to be viewed at a critical angle $\theta_{\rm c}$. Then, the corresponding view angle $\theta_{\rm c}$, the angle between the direction of jet velocity and the line of sight to the observer, can be determined as ${\theta_{\rm c}=\displaystyle{{\rm arcsin}\left(1\over\gamma_{\rm min}\right)}\approx 4}$° (Lu, Krichbaum & Zensus, 2011). Thus, the minimum Doppler factor, ${\displaystyle\delta\equiv{1\over\gamma}{1\over\left(1-\beta{\rm cos}\theta\right)}}$, would be $\delta_{\rm min}=\gamma_{\rm min}\approx 14$.

The enhanced factor of observed luminosity for a relativisitic jet is: $\delta^{p}$ with $p=3-\alpha_{\rm c}$ for a spherical core and $p=2-\alpha_{\rm 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 $\gamma^{\alpha_{\rm c}-3}(1-\beta{\rm cos}\theta)^{\alpha_{\rm c}-3}$ (Kellermann & Owen, 1988; Cawthorne, 1991; Urry & Padovani, 1995). In the case of $\gamma=\gamma_{\rm min}=14$, the flux densities from the inner components of the core can be boosted by $\sim\gamma_{\rm min}^{3}$, about 3 orders of magnitude, assuming $\alpha_{\rm 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 $\sim\gamma_{\rm min}^{-3}$ 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 1027 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 MOJAVE 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 jet-flow. 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.

4.2 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 narrow-angle-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 mirror-symmetric 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_{\rm hspt}$= 150 is owing to the relativistic Doppler effect, then we can estimate the bulk Lorentz factor $\gamma$ (e.g. Laing & Bridle, 2015). The emission from the western jet, plausibly moving towards the observer, is enhanced by $\delta^{2.94}$ for a coninuous jet with spectral index $\alpha_{\rm Jet}\approx\alpha_{\rm 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 $\displaystyle{\left[\gamma_{\rm min}(1-\beta_{\rm a}^{2})\right]^{-2.94}}$ for the approaching western jet while the receding eastern jet is suppressed by a factor of $\displaystyle{\left[\gamma_{\rm min}(1+\beta_{\rm a}^{2})\right]^{-2.94}}$. From the intensity ratio $R_{\rm hspt}$ for W12 and E5, the bulk Lorentz factor of $\gamma_{\rm min}\sim\sqrt{{R_{\rm hspt}^{1/2.94}+1\over 2}}=1.8$ can be inferred. Then, the view angle of the hotspot W12 can be estimated as $\theta_{\rm c}={\rm sin}^{-1}\left(\displaystyle{1\over\gamma_{\rm min}}\right)=34$°. Compared to $\theta_{\rm 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).

5 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.

@firstsectionfalse

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 point-like source at an angular resolution of $\sim$1 arcsec when the DR is below 10,000:1. In addition, at sub-arcsec resolution, residual delays ($\delta\tau_{\rm 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 multiple-scale (MS) imaging (Rau & Cornwell, 2011). The effect becomes aggravated and significantly degrades the image quality in the field near a strong compact source. For example, Fig A2 shows residual bandpass solutions in phase, solved from JVLA A-array data of NRAO 530 at 9 GHz after the standard delay corrections derived from 3C 286 are applied. A phase slope of 2° across a sub-band of BW 128 MHz is present, suggesting that a typical residual delay of 0.04 nsec or 1/3 $\lambda$ remains in the data. A time variation of the residual delays is also evident from the time-dependent bandpass solutions. A larger phase slope of 4° per 128 MHz is present occasionally on some antennas (see Antenna=’ea08’ Time = ’11:12:35.1’ of Fig. A2, bottom-left). As a consequence, the presence of residual delays in deep image has several implications.

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 ($\Delta\phi$) caused by residual delays ($\delta\tau_{\rm A}$) across a sub-band (${\rm BW}=\Delta\nu$) 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 ($\Delta\theta$) on the radio core can be related to the quality of data described by the residual delays ($\delta\tau_{\rm A}$) remaining after the calibration process, $\Delta\theta={\displaystyle{c\delta\tau_{\rm A}\over 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 $\sim$0.07″, about a quarter of the synthesized beamwidth, $\theta_{\rm 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 $\delta\tau_{\rm A}$, the positional error $\Delta\theta$ 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, $\theta_{\rm syn}={\displaystyle{\lambda\over B}}$, where $\lambda$ is the observing wavelength. The baseline term is cancelled in the ratio of the positional error to the size of the sythesized beam, namely, ${\displaystyle{\Delta\theta\over\theta_{\rm syn}}}={\displaystyle{c\delta\tau_{\rm A}\over\lambda}}$. 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 $\delta\tau_{\rm 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 vector-averaging technique for sub-bands is that the phase variation is effectively turned into signal noise, producing a coherence loss. For a visibility function $V(\nu)=A(\nu){\rm exp}[i\phi(\nu)]$, the ratio of the vector-averaged amplitude to the actual amplitude drops with the amount of phase changes in phase across the sub-band ($\Delta\phi$) due to a residual delay: ${\displaystyle{<A>\over A}=1-{\Delta\phi^{2}\over 24}}$. 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 $\sim$0.06% of the core flux-density ($F_{\rm core}$). The uncertainty of the vector-averaged amplitude is ${\displaystyle\sigma=A{\Delta\phi\over 2\sqrt{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 $\displaystyle\left[{\rm sin}(\pi\Delta\nu_{\rm BW}\tau_{rs})\over\pi\Delta\nu_{\rm BW}\tau_{\rm rs}\right]$, the delay beam or the bandwidth pattern (Thompson, 1986). The residual delay here ($\tau_{\rm rs}$) is considered from two aspects: (1) The residual delay $\delta\tau_{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, $\delta\tau_{\rm pnt}=\displaystyle{B\over c}\Delta\theta_{\rm pnt}$, for a source offset by $\Delta\theta_{\rm pnt}$ from the phase tracking center or the array pointing center. Hereafter, we refer to the geometrical delay, $\delta\tau_{\rm pnt}$, as the source delay. Both delay terms, $\delta\tau_{A}$ and $\delta\tau_{\rm pnt}$, are time variable. Given the unclear origins, the behavior of $\delta\tau_{A}$ appears unpredictable in time. The term $\delta\tau_{\rm 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 $\Delta\nu_{\rm 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 bandwidth-smearing effect. At the radial distance of $\Delta\theta_{\rm pnt}$ away from the phase center, a point source will be smeared to an angular extent of $\Delta\theta_{\rm pnt}\displaystyle{\left(\Delta\nu_{\rm BW}\over\nu\right)}$ (Thompson, 1986). For example, at $\Delta\theta_{\rm pnt}=10$” from the phase center, the angular distortion of a point source can be $\sim$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, $\delta\tau_{\rm A}$, varies with time in an unpredictable way. The unclear origin of $\delta\tau_{\rm A}$ causes difficulties in eliminating such frequency-dependent phase errors from the data of target sources. The phase issue caused by $\delta\tau_{\rm 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 ($\sim\nu^{\alpha}$) with $\alpha=-1$, the fractional amplitude changes by $\Delta S/S=\Delta\nu/\nu\sim 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 $\left[u,v\right]$-plane of a radio source sampled with a wide-band interferometer at observing time $t$ and frequency $\nu$ can be written as:

[TABLE]

where $\tau_{A}(t)$ and $G_{A}(\nu,t)$ correspond to time delays and time-dependent complex gains, respectively. Both terms are antenna-based, containing instrumental errors and atmospheric effects. The source function, $V^{\rm S}(u,v,\nu,t)$, composed of numerous emission components, can be described as:

[TABLE]

where $V_{i}^{\rm Ext}({\rm u},{\rm v},\nu)$ is the source visibility function of the $i$-th extended component, and $V_{i}^{\rm C}({\rm u},{\rm v},\nu,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^{\rm C}({\rm u},{\rm v},\nu,t)$ and an arbitrary extended emission $V^{\rm Ext}({\rm u},{\rm v},\nu)$, so that,

[TABLE]

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^{\rm S}({\rm u},{\rm v},\nu,t)\sim V^{\rm C}({\rm u},{\rm v},\nu,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:

[TABLE]

A.2 Building a point source model

Dividing the observed complex visibility $V({\rm u},{\rm v},\nu,t)$ of Eq(A1) by the source model $V^{\rm S}({\rm u},{\rm v},\nu)$ of Eq(A4), one can approach to point-like visibility function:

[TABLE]

where the normalized source function $V^{\rm P}({\rm u},{\rm v},\nu,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, $V_{0}^{C}({\rm u},{\rm v},\nu)=V^{\rm C}({\rm u},{\rm v},\nu,t_{0})$, the flux-density at the beginning of the observation, $t_{0}$, and a time-variable term giving the deviations from that initial value, $\Delta V^{C}({\rm u},{\rm v},\nu,t)$:

[TABLE]

For calibration concerns, the radio-dominant core is placed at the phase tracking center of the interferometer. Given a variable intensity $I_{0}+\Delta I$ of the radiation from its core, the complex finction (A6) of the compact source can be simplified as a real function, $V^{\rm C}({\rm u},{\rm v},\nu,t)=I_{0}+\Delta I(t)$. If the extended component $V^{\rm Ext}({\rm u},{\rm v},\nu)$ is not time variable, then the normalized source function of Eq(A5) can be expressed as,

[TABLE]

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^{\rm P}({\rm u},{\rm v},\nu,t)$ is essentially a complex function tracing the antenna gains multiplied by the antenna-based delay term, i.e.,

[TABLE]

the right of Eq(A8) can be defined as the overall antenna-based complex gain $\mathscr{G}_{\rm A}({\rm u},{\rm v},\nu,t)\equiv e^{2\pi i{\rm\nu\tau_{\rm A}(t)}}G_{\rm A}(\nu,t).$

Furthermore, the antenna-based delay can be further separated into two terms, a constant delay $\tau_{0}$ and a time-dependent residual delay:

[TABLE]

For data from the JVLA, $\tau_{0}$ can be removed by the standard delay correction in the beginning of data calibration or in a pipeline process, and $\delta\tau_{\rm 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_{\rm A}(\nu,t)=G_{\rm A}(t)B_{\rm A}(\nu)$, the overall antenna-based gains becomes

[TABLE]

where $G_{A}(t)$ is the antenna-based complex gain, a function of time; with a scaling factor of the amplitude $\mathscr{F}$, the antenna gain can be expressed as,

[TABLE]

and $B_{\rm A}(\nu)$ is the antenna-based, instrumental complex bandpass, a function of frequency,

[TABLE]

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 pipeline999see https://science.nrao.edu/facilities/vla/data-processing/ are essentially to correct for the parameter terms $\tau_{0}$, $b_{\rm A}(\nu)$, $\phi_{b_{\rm A}}(\nu)$, $g_{\rm A}(t)$, $\phi_{g_{\rm A}}(t)$ and $\mathscr{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_{\rm A}(t)$ and $B_{\rm A}(\nu)$, or Eq(A11) and Eq(A12), can be assessed by the normalized gains, $g_{\rm A}(t)=1+O[g_{\rm A}(t)]$ and $b_{\rm A}(\nu)=1+O[b_{\rm A}(\nu)]$ for amplitude, as well as $e^{i\epsilon_{\phi_{g_{\rm A}}}(t)}$ and $e^{i\epsilon_{\phi_{b_{\rm A}}}(\nu)}$ for phase. The residual errors are small, i.e., the residual amplitude errors $O[g_{\rm A}(t)]<<1$ and $O[b_{\rm A}(\nu)]<<1$ and the residual phase error $\epsilon_{\phi}\equiv\epsilon_{\phi_{g_{\rm A}}}(t)+\epsilon_{\phi_{b_{\rm A}}}(\nu)$ fluctuates around zero as a function of time and frequency. The point source model as described in Eq(A8) is reduced to be a function of residual errors:

[TABLE]

neglecting the terms with high order of residual errors in amplitude. If the residual errors $O\left[b_{\rm A}(\nu)\right]$ in frequency-dependent gains are much smaller than the residual errors $O\left[g_{\rm A}(t)\right]$ in the time-dependent gains and the phase term owing to residual delay dominates frequency-dependent residual errors in phase, ignoring $O\left[b_{\rm A}(\nu)\right]$ and $\epsilon_{\phi}$ in Eq(A13), the first order of the residual errors is reduced to:

[TABLE]

The remaining residual errors can be separated, namely the residual delay $\delta\tau_{\rm A}(t)$ and the residual gain $O\left[g_{\rm A}(t)\right]$. Both are time dependent. The residual gains can be corrected with self-calibration, the well-known technique (Schwab, 1980) that has been successfully applied in continuum imaging for radio interferometer arrays. The residual delay $\delta\tau_{\rm A}(t)$ can be fit to time-dependent 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 $\delta\tau_{\rm 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 RE-correcting 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 time-variable 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 Ka-band 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 (time-dependent) 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 (step 3), the improved image is constructed, showing an image DR $\sim$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 ${\displaystyle 2\left(S_{max}-S_{min}\over S_{max}+S_{min}\right)=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 flux-density 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 errors101010The closure errors referred here are baseline-based errors from the possible origins as discussed in Cornwell (1986) and Clark (1981) for the former VLA. Unlike the antenna-based errors that are vanished in the well-known closure relationships, e.g. Eqn. (9.104) and Eqn. (10.44) of Thompson, Moran & Swenson (2017) for phase and amplitude, and can be corrected with self-calibration technique (Schwab, 1980; Cornwell, 1986), the baseline-based errors remain in the closure relationships and the self-calibration will fail to correct them. Baseline-based corrections for the closure errors can be done, however (Perley, 1986). 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 $\sim$10 $\mu$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 bandpass 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 sub-bands, 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 lower-level 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. M8 - a new imaging model created from the data corrected for antenna-based error. This model is used to correct for baseline-based 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.