Interactions between the jet and disk wind in a nearby radio intermediate quasar III Zw 2

We study the low-frequency radio structure of III Zw 2 based on the recent observations from the Murchison Widefield Array (MWA, Tingay et al. 2013), the Giant Metrewave Radio Telescope (GMRT, Swarup et al. 1991) and the Australian SKA Pathfinder (ASKAP, Johnston et al. 2007; Hotan et al. 2014).

The MWA is the precursor of the Square Kilometre Array (SKA) low-frequency telescope and is located in Western Australia. We have downloaded the latest unpublished data of III Zw 2 observed by the MWA Phase II (Wayth et al., 2018) on 2018 June 15. The extended array of the MWA phase II comprises 128 tiles (each consisting of 16 crossed-pair dipole antenna), distributed over a maximum baseline of $\sim$6 km, twice longer than that of Phase I. These data are included in the GaLactic and Extragalactic All-sky MWA survey – eXtended (GLEAM-X, Hurley-Walker et al. 2022). The GLEAM-X survey covers the entire sky south of declination $+30^{\circ}$ and at the same frequency range of 72–231 MHz as GLEAM (Hurley-Walker et al., 2017). GLEAM-X achieves a typical sensitivity of $1\sigma\approx 1.3$ mJy beam^-1 in the frequency range of 170–231 MHz, which is eight times more sensitive than GLEAM. We downloaded 40 observation snapshots containing III Zw 2 from the MWA archive ³³3https://asvo.mwatelescope.org, and processed the data at the China SKA Regional Centre (An et al., 2019, 2022) using the pipeline developed by the GLEAM-X team⁴⁴4https://github.com/tjgalvin/GLEAM-X-pipeline maintained by Tim Galvin. For each snapshot, we first flagged the bad tiles or receivers and obtained reliable calibration solutions. Next, we flagged potential radio frequency interferences. The calibration solutions obtained from the calibrators were then applied to the data. A Briggs robust parameter was set to $+0.5$ in imaging to maximize sensitivity. We made a deep cleaned image, followed by source-finding, ionospheric de-warping, and flux density scaling to GLEAM in the clean images. Then the point spread function (PSF) placeholder of each snapshot image was corrected. In the last step, we combined the astrometrically- and primary-beam-corrected snapshot images into a high signal-to-noise image with reduced noise, allowing for detecting fainter sources and diffuse structures that are not visible in individual snapshot images.

Figure 6-a shows the 216-MHz MWA image of III Zw 2, displaying a compact component. The resolution of the image is $63\farcs 1\times 49\farcs 1$, and the root-mean-square (rms) noise is 2.1 mJy beam^-1 (a factor of $\sim 4.8$ lower than that of the GLEAM image). The rms noise of the III Zw 2 image is 1.7 times higher than the mean value of the GLEAM-X images, due to the lower elevation angle of the MWA when observing III Zw 2 at $\delta=+11^{\circ}$. III Zw 2 is marginally resolved and shows an elongation along the north-south direction. The source is unresolved at other lower frequencies. The GLEAM-X data points themselves can be fitted with a power law with a steep spectral index $\alpha=-0.96\pm 0.09$.

Figure 6-b shows the GMRT image at 150 MHz, showing a resolved structure along the NE-SW direction with a total extent of $\sim 48\arcsec$. The GMRT image of III Zw 2 is from the TIFR GMRT Sky Survey Alternative Data Release (TGSS-ADR1 ⁵⁵5https://vo.astron.nl) (Intema et al., 2017). TGSS-ADR1 covers 90 per cent of the total sky from $\delta=-53^{\circ}$ to $\delta=+90^{\circ}$ with a noise level of $\sim$5 mJy beam^-1 and a resolution of $25\arcsec\times 25\arcsec$/cos(Dec$-19^{\circ}$) at 150 MH. The total flux density is $\sim$ 212 mJy, dominated by the central core. Besides the core C, there are extensions toward the northeast (NE) and toward the south (S).

We obtained the ASKAP image of III Zw 2 from the Rapid ASKAP Continuum Survey (RACS), which is a shallow all-sky (covering the entire sky south of declination $\delta=+41^{\circ}$) pilot survey for future multi-year surveys with the full-scale ASKAP (McConnell et al., 2020). The observations covering the III Zw 2 position began on 2019 April 21 and lasted three weeks. The RACS has an instantaneous bandwidth of 288 MHz and is centered at 888 MHz. The angular resolution of the resulting image using the natural weighting is $14\farcs 20\times 12\farcs 93$. The distance between NE and SW lobes is $\sim 33\arcsec$ and the image resolution enables to resolve the jet structure. The rms noise in the image is 0.33 mJy beam^-1, which is consistent with the mean rms of the RACS images. The ASKAP image is shown in Fig. 6-c, which reveals much richer details than the MWA image: the main jet body is elongated along the northeast–southwest (NE–SW) direction and consists of the core C and NE and SW lobes. The ASKAP image shows a very similar structure to the recently published uGMRT image in Silpa et al. (2020), although the two images are observed at different frequencies. As the observed frequency of the ASKAP image is higher than that of the uGMRT image, the emission from the outermost extended components beyond the NE and SW lobes (labeled as N and S in Fig. 6-c) is fainter due to their steep spectrum nature. The separation between N and S lobes is $\sim 55\arcsec$ ($\sim$98 kpc)⁶⁶6Adopting the following cosmological parameters: H₀ = 71 km s^-1 Mpc^-1, $\Omega_{\Lambda}=0.73$ and $\Omega_{m}=0.27$, at $z=0.0893$, 1 mas angular separation corresponds to 1.65 pc linear size in projection on the plane of the sky., larger than the size of the previously detected structure between NE and SW lobes in VLA images at GHz frequencies (Brunthaler et al., 2005). In the 685-MHz uGMRT map, the southern extension (labeled ML in Silpa et al. 2021) is much longer than our Figure 6-c. The total flux density estimated from the ASKAP images is 92.5 mJy, of which the main body (C+NE+SW) accounts for $\sim$86.3 mJy, while the remaining $\lesssim$6.7% (6.2 mJy) of the flux density comes from the sum of the N and S lobes.

The power-law portion at low frequencies originates from large-scale extended emission structures in the outermost N and S lobes which are evident from the MWA image of III Zw 2 and part of the structures are shown in the uGMRT (Silpa et al., 2021) and ASKAP images. These steep-spectrum features correspond to aged structures, which become very faint above 1 GHz (e.g. Callingham et al., 2017).

The flux density from optically thin synchrotron emission is (Rybicki & Lightman, 1986; Ghisellini, 2013)

where $R$ is the size of the extended emission structure, $D_{L}$ is the luminosity distance, $K_{0}$ is a normalization factor, $B$ is the magnetic field strength, $\nu$ is the observing frequency, $p$ is the energy index, and $c1$ and $c2$ are radiative constants (Pacholczyk, 1970). Assuming a similar power law distribution of synchrotron emitting electron energies, $N(E)~{}dE=K_{0}E^{-p}~{}dE$, and energy equipartition between the magnetic fields and particle kinetic energy densities, the normalization $K_{0}$ can be evaluated similarly to eqn. E2. The synchrotron frequency corresponds to the minimum injected Lorentz factor of emitting electrons and can be calculated as

Using eqn. A2 and the normalization in eqn. A1, and $\gamma_{m}=((p-2)/(p-1))(m_{p}/m_{e})\epsilon_{e}$, in which $\epsilon_{e}$ is the fraction of the total energy density in the emitting region in the magnetic fields. The magnetic field strength in this region can be estimated as

assuming $\epsilon_{e}=0.1$ and an index $p=1-2\alpha_{\rm low}=2.94$ based on the inferred spectral index $\alpha_{\rm low}=-0.96$, and is indicative of emission from a weakened decelerating shock produced by Fermi acceleration of electrons (e.g. Blandford & Eichler, 1987; Jones & Ellison, 1991; Frank et al., 2002). The estimated $\mu$G magnetic field strength is consistent with similar estimates for this source from other studies (Silpa et al., 2021) and AGN in cluster environments (Müller et al., 2021). Assuming a spherical volume, the total energy content in the synchrotron emitting extended region is

where, $U_{e}=\epsilon U_{B}=\epsilon B^{2}/(8\pi)$, where $\epsilon=\epsilon_{e}/\epsilon_{B}$ and using $B$ from eqn. A3. After eqn. A4 relating to the total energy $E_{\rm ext}$, from the mean flux density measured in the MWA bands (88 – 215 MHz), $\left<S_{\nu}({\rm ext})\right>=0.246$ mJy, and using a spectral index $\alpha_{\rm low}=-1.15$ for the relics and a central frequency $\nu=151.6$ MHz, the radio luminosity of the extended emission $L_{\rm R,ext}\approx 4\pi D^{2}_{L}\nu\left<S_{\nu}({\rm ext})\right>(1+z)^{-1-\alpha_{\rm low}}=6.18\times 10^{36}$ erg s^-1. Using the empirical relations in eqn. E8 and eqn. E9, the associated luminosity of the extended jet is $L_{\rm ext}=5.19\times 10^{39}$ erg s^-1.

The source structure of III Zw 2 on parsec scales is revealed from the VLBI imaging data, which include the archive data from the Monitoring Of Jets in Active galactic nuclei with VLBA Experiments (MOJAVE) (Lister et al., 2018) program and the archive data obtained from the Astrogeo database⁷⁷7http://astrogeo.org/. Details of the VLBI data are presented in Table 1. The MOJAVE data of III Zw 2 were observed at 15 GHz over 25 sessions from 1995 July 28 to 2013 June 2. The Astrogeo data were observed simultaneously at the dual frequencies of 2.3 and 8.4 GHz, enabling the determination of spectral indices $\alpha_{\rm 2GHz}^{\rm 8GHz}$ for III Zw 2. Since the source is unresolved in 2.3 and 8.4 GHz images, the analysis of the jet kinematics is based only on the 15 GHz VLBI data.

All of these archival VLBI data have been calibrated, so we only made a few iterations of self-calibration in the Difmap software package (Shepherd et al., 1995) to eliminate some residual phase errors and to increase the dynamic range of the images. Model fitting was performed using the MODELFIT program in Difmap. Lister et al. (2019) studied the parsec-scale jet kinematics of 409 bright radio-bright AGN including III Zw 2. They presented the model fitting results of III Zw 2 from epoch 1995 July 28 to epoch 2013 June 2. In their model fitting, (1) the source was detected with only a single core before 2011 January 11; (2) the distance of the fitted jet component from epoch 2000 July 22 to 2006 June 15 is less than 0.32 mas, i.e., less than half of the minor axis of the synthesized beam, therefore, these components are practically unresolvable. These above epochs were abandoned from the jet kinematics analysis. The model fitting parameters are given in Table 2. The uncertainties of the parameters are first estimated from the equations given in Fomalont (1999) which only take into account the fitting errors. In practical observations, the flux density error also includes a certain level of uncertainty in the flux scale calibration through error propagation. For the 15-GHz VLBA, this systematic error is typically about 5%. Both total intensity (Stokes I) and linear polarization images are created from the MOJAVE data. After self-calibration in Difmap, we separately created images using the Stokes I, Q and U data. Then we combined the Stokes Q and U images to produce the linear polarization intensity $p$, which is the root of the sum of squares of Q and U components, i.e., $p=\sqrt{Q^{2}+U^{2}}$. The fractional polarization was calculated as the ratio of the polarized flux density to Stokes I flux density, $p/I$. The derived fractional polarizations are consistent with the values reported on the MOJAVE webpage. Snapshot-mode Astrogeo data have short integration time and sparse (u,v) coverage. Deconvolution and model fitting are affected by strong sidelobes caused by incomplete (u,v) coverage. For this reason, the uncertainty in component size adopts the root of the quadratic sum of its corresponding statistical error and the fitting error. The position error was then estimated as half of the component size error.

The present paper is focused on the 2009–2010 major flare and its associated jet properties, so we only use the MOJAVE data after 2004, including 13 epochs listed in Table 1. The (u,v) coverage, image sensitivity, and resolution of these VLBI data are all in good agreement, therefore the systematic errors between the fitted parameters are small. In six of the 13 epochs, the jet is not distinguishable, and these data were not used in jet kinematics analysis.

On pc scales, VLBI images of III Zw 2 reveal either a single core or “a core + a compact jet” structure. The compact core dominates the total flux density in all images. The jet extends to the west and is relatively weaker but well distinguished from the core in 7 epochs, allowing us to determine its kinematic properties. We double-checked the model fitting reported by the MOJAVE team and only adopted the high-confidence jet models. If the signal-to-noise ratio of a jet component is lower than 3, or the jet is indistinguishable from the core, then it will not be used for the kinematic analysis. Figure 2 shows the variation of the core-jet separation with time and the variation of the jet position angle with time. It is clear that these jet knots do not belong to a single component, but are two separated jet knots, labeled J1 and J2 in the plot (corresponding to components #1 and #4 in Lister et al. 2019), each of which follows a ballistic trajectory along a different position angle. We made linear regression fitting to the position-time correlation of J1 and J2 and obtained the jet proper motions as $0.24\pm 0.02$ mas yr^-1 and $0.22\pm 0.01$ mas yr^-1, respectively. These convert to jet transverse speeds of $1.35\pm 0.13\,c$ (J1) and $1.21\pm 0.07\,c$ (J2). The derived jet proper motions are in good agreement with the previous studies (Brunthaler et al., 2000) based on the 43-GHz VLBI observations in 1998. To distinguish it from our jet components, we name the jet detected by Brunthaler et al. as “1998 jet”. Extending the 1998 jet to epoch 2006 (the first epoch of the MOJAVE data used in this paper), we find the component should be at a distance of about 0.525 mas. However, it is not detected in the 2006 VLBA image, suggesting that this component has been considerably dimmed due to the adiabatic radiation loss.

The single-dish light curve data used in this paper were from the archives of the 14-meter Metsähovi radio telescope in Finland and the 40-meter telescope of the Owens Valley Radio Observatory (OVRO) in US.

Metsähovi observations of III Zw 2 were made at 22 and 37 GHz in the period of 1986–2020. The 22- and 37-GHz light curves from 1984 to 1998 have been published in Falcke et al. (1999), and the Metsähovi data from 1986 to 2019 have been published in Chamani et al. (2020). In addition to these data, we also include the new data from 2019 onwards to date. A detailed description of the data reduction and analysis is given in Teraesranta et al. (1998). Under good observing conditions, the detection limit of the Metsähovi telescope at 37 GHz is around 0.2 Jy. Uncertainties in flux densities include the contribution from the root mean square (rms) of the measurements and the errors in the absolute flux density calibration.

The radio monitoring program with the 40 m telescope at the OVRO includes over 1500 sources in the northern sky ($\delta>-20^{\circ}$) (Richards et al., 2011). Each source is observed twice per week. The minimum flux density detected is about 4 mJy. A typical uncertainty of the measurements is $\sim$3%. The high cadence and high sensitivity of the monitoring greatly facilitate the study of the variability properties of radio sources on timescales of months and years. The OVRO monitoring data of III Zw 2 were observed at 15 GHz from 2008 to 2020. The maximum and minimum flux densities are 1825 mJy and 72 mJy, respectively. The integrated flux densities from the 15 GHz VLBA data are superimposed on the OVRO light curve, showing a good match between the two sets of flux densities at adjacent epochs. This suggests that the steep-spectrum extended structures including the kpc-scale jet and lobes seen in the low-frequency images are barely detectable at 15 GHz and above.

The time series analysis is carried out using two methods, the Fourier periodogram (Scargle, 1989; Vaughan, 2005; An et al., 2013; Mohan & Mangalam, 2014) and the wavelet analysis (Torrence & Compo, 1998; An et al., 2013; Mohan et al., 2016b). For both methods, the light curves are first pre-processed by a conversion from a partial uneven sampling with small data gaps to a full even sampling after a linear interpolation and re-sampling. The normalized Fourier periodogram $P(f_{j})$ is evaluated as (Vaughan, 2005; Mohan & Mangalam, 2014; Mohan et al., 2015)

where $\Delta t$ is the sampling time step size, $\overline{x}$ is the mean of the light curve $x(n\Delta t)$ of length $N$, and $F(f_{j})$ is the Fourier transform of $x$ evaluated at the frequencies $f_{j}=j/(N\Delta t)$ with $j=1,.....,(N/2-1)$ (up to the Nyquist frequency). Since astrophysical processes typically produce red noise in the light curve, the periodogram $P(f_{j})$ will be characterized by a power law behavior, especially at low temporal frequencies. This can be captured by parametric model fits to $P(f_{j})$ to estimate the underlying power spectral density (PSD). Here, we use two models, a power law given by

where $A$ is the normalized amplitude, $\alpha$ is the power-law index, $C$ is the ambient noise level, and a bending power law is given by

where $f_{b}$ is a break frequency marking a transition in slopes, the power-law indices are $\alpha$ for $f_{j}>f_{b}$, and $\alpha_{l}$ for $f_{j}<f_{b}$, and $C$ is the ambient noise level.

The fitting of the Fourier periodogram with the above parametric models is carried out using the maximum likelihood estimator. The methodology and estimation of parameters and their errors are discussed in (Mohan & Mangalam, 2014; Mohan et al., 2015). After accounting for model uncertainties, statistical significance testing is carried out to identify outlying peaks in the fit residuals (based on an expected $\chi^{2}$ statistical distribution), which are potential candidates for quasi-periodic oscillations in the light curve.

The statistical significance of all peaks is assessed based on a procedure involving Monte Carlo simulations (Mohan & Mangalam, 2014; Mohan et al., 2015). These are carried out using the algorithm of Timmer & Koenig (1995). The best-fit model and associated parameters are used as trial values to simulate 5000 realizations of the periodogram, oversampled in duration and temporal frequencies, and then re-sampled at the original frequencies. A mean periodogram $\tilde{I}(f_{j})$ is constructed from the simulations and is re-scaled to match the variability properties of the original periodogram; this could be the closest estimate of the underlying PSD. The statistical significance of the periodogram ordinates for any frequency bin is evaluated based on the assumption that the light curve consists of randomly distributed data points (i.e., no periodic behavior). The residuals from the fitting, $P(f_{j})/\tilde{I}(f_{j})$, are then $\chi^{2}_{2}/2$ distributed (Chatfield, 2016), with a conditional probability $p[P(f_{j})|\tilde{I}(f_{j})]=\tilde{I}(f_{j})^{-1}e^{-P(f_{j})/\tilde{I}(f_{j})}$. The cumulative distribution function for the PSD ordinates is then the integral of the $\chi^{2}_{2}$ distribution, i.e., a gamma density function $\Gamma(1,1/2)=\exp{(-x/2)}/2$. Specifying a level of statistical significance $(1-\epsilon)$ in the integral helps to identify outliers (quasi-periodic signals) that may be present in the tail of the distribution. We set a threshold of $\epsilon=0.05$ (95 % level of significance) to identify quasi-periodic signals in the light curves.

The wavelet analysis is employed here in a complementary manner to probe quasi-periodic signals in the light curves and their locations (to infer the total duration and number of cycles). The two-dimensional wavelet power spectrum is a function of the wavelet scale (that can be expressed in units of the sampling wavelength or period) and the time window being sampled, and is evaluated as (Mohan et al., 2016b)

where $W(n,s)$ is the wavelet transform of the evenly sampled light curve $x(n\Delta t)$ at times $(n\Delta t)$ and at scales $s$, $F(2\pi f_{j})$ is the Fourier transform of $x$ evaluated at the circular frequencies $2\pi f_{j}=2\pi j/(N\Delta t)$ with $j=1,.....,N$, and $\Psi^{\ast}(2\pi sf_{j})$ is the complex conjugate of the Fourier transform of the wavelet sampling kernel function. The wavelet transform is the inverse Fourier transform of the convolution product of the above constituents. For a continuous wavelet transform, a commonly used sampling kernel is the Morlet wavelet function (Grossmann et al., 1989) $\psi=\pi^{-1/4}e^{i\omega_{0}t}e^{-t^{2}/2}$ where $\omega_{0}=6$ is a frequency parameter characterizing the wavelet shape and $t$ is the time parameter. In the frequency domain $\Psi(2\pi sf_{j})=\pi^{-1/4}e^{-(2\pi sf_{j}-\omega_{0})^{2}/2}$. The wavelet scales $s\approx 1.03/f_{j}$ for the Morlet function (Torrence & Compo, 1998) is hence in near correspondence with the sampling frequencies. In the current work, the following additional features have been implemented: the use of a cone of influence to help explain the cyclic behavior of the sampling kernel, especially at low temporal frequencies (Torrence & Compo, 1998), the use of a Hann window function to smoothen the noisy features in the power spectrum, and the identification of contiguous features in the power spectrum that may be statistically significant in anticipation of measuring their total duration, number of cycles and the time evolution of the features. These improvements all tend to narrow the search window and consequent computational costs in the identification of statistically significant signals, and considerably improve the contrast between the signal and noise.

The statistical significance of contiguous signals detected in the wavelet analysis employs a two-fold strategy. In the first stage, the algorithm of Timmer & Koenig (1995) and the expected $\chi^{2}_{2}$ statistics of periodogram ordinates (assuming that the light curve ideally consists of random Gaussian noise) are employed to simulate a large number of realizations of the periodogram. The best-fit power-law model of the form $I(f_{m})=Af^{\alpha}_{m}+C$ with associated parameters is used as trial values. The index $\alpha$ is varied in the range from the best-fit value to 0.0 in steps of $-0.2$ and in each case. 1000 realizations of the periodogram are simulated, and they are oversampled in duration and temporal frequencies and then re-sampled at the original frequencies. In each realization, the periodogram is inverse Fourier transformed to obtain a synthetic light curve with similar statistical and variability properties as that of the original light curve; the total number of simulated light curves is typically $\geqslant 20000$. In the second stage, their individual wavelet power spectra are evaluated. For each simulated wavelet power spectrum, the mean of all powers at a given scale $P(n)$ is used to estimate the global wavelet power spectrum GWPS$(s)$ that corresponds to a window function smoothed version of the Fourier periodogram. The candidate signals in the GWPS of the original light curve are compared with their simulated counterparts. The number of times $p$ that a candidate signal in the original GWPS exceeds the values in the simulations (totaling $Q$) at a given wavelet scale is measured in terms of the statistical significance $(1-p/Q)$.

We plot the radio spectrum of III Zw 2 from 72 MHz to 37 GHz in Figure 4. The red-colored squares below 230 MHz are taken from the MWA GLEAM-X survey (Hurley-Walker et al., 2022) and the green-colored diamond is from the GMRT observation (Intema et al., 2017), revealing the extended radio emission characterized by a steep spectrum. The data point in slateblue color is from the ASKAP observation at 888 MHz. The solid circle in magenta color is from the recent uGMRT observation (Silpa et al., 2020) on 2018 November 23. In addition to these low-frequency data, we also include the data points at GHz frequencies observed in the epochs close to the MWA and uGMRT observation: 37 GHz Metsähovi (blue-colored triangle-left), 15 GHz OVRO (navy-colored triangle-right), 8.4 and 2.3 GHz Astrogeo VLBI (purple-colored diamond). From Figure 3 we find that the core dominates the total flux density at GHz frequencies.

The spectrum shows different characteristics above and below 685 MHz: at 685 MHz and above, the core dominates the emission, showing an inverted spectrum; at frequencies below 685 MHz, the flux density from the core is substantially reduced toward the low-frequency end due to increasing synchrotron self-absorption, and the extended jets and lobes become increasingly dominant. The two-component radio spectrum spanning 0.072 GHz – 37 GHz is subjected to a weighted (the measurement errors of flux densities are used as weights) least square fit using the function eqn. D1, in which the first two parts ($A\nu^{\alpha_{\rm low}}+B$) form a power-law spectrum and the last one describes a spectrum of self-absorbed synchrotron radiation (Pacholczyk, 1970; Türler et al., 1999).

where the parameters of the low-frequency power-law spectrum include the amplitude $A$ (mJy), the spectral index $\alpha_{\rm low}$, a baseline flux density $B$ (mJy), and the parameters of the high-frequency section include the amplitude $F_{m}$ (mJy), the frequency of transition from optically thick to thin emission $\nu_{m}$ (GHz), the optically-thick spectral index $\alpha_{\rm thick}$, the optically-thin spectral index $\alpha$, and the optical depth $\tau_{m}$ expressed in terms of $\alpha_{\rm thick}$ and $\alpha$. In the fitting process derived by using the Markov chain Monte Carlo (MCMC) method, $\alpha_{\rm thick}$ is fixed at 2.5, corresponding to the canonical synchrotron self-absorption case, and other parameters are constrained within reasonable ranges $0.0<A\leq 0.8$ Jy, $-1.5\leq\alpha_{\rm low}\leq-0.4$, $0.0<B\leq 0.6$ , $0.0<F_{m}\leq 0.5$ Jy, $7.0\leq\nu_{m}\leq 15.0$ GHz, and $-0.4\leq\alpha\leq 0.0$. This yields $A=27^{+12}_{-9}$ mJy, $\alpha_{\rm low}=-0.97^{+0.21}_{-0.21}$, $B=41^{+8}_{-12}$ mJy, $F_{m}=146^{+11}_{-9}$ mJy, $\nu_{m}=11.24^{+2.01}_{-1.25}$ GHz, and $\alpha=-0.20^{+0.10}_{-0.12}$.

In Figure 4 we also plotted the flux densities of the NE and SW lobes obtained from the VLA images (Brunthaler et al., 2005) and fitted the NE and SW lobes (represented by blue-colored open circles and red-colored open squares, respectively) with power law functions. We then calculated the flux densities of NE and SW lobes at the MWA observation frequencies (i.e., 88, 118, 154, and 200 MHz) by extrapolating the power-law fits. Next, we subtracted the extrapolated flux densities of NE and SW lobes from the measured values to estimate the flux densities from large-scale extended emission contributed mainly by the N and S lobes. Finally, we calculated the spectral index of the extended outer lobes N+S as $\alpha=-1.09\pm 0.12$, which is much steeper than the inner lobes and comparable to the spectral index of the recently observed radio relics at low frequencies, indicating that the outer lobes are more likely to be dying relics, dominated by radiation from aged relativistic electrons.

The synchrotron emission spectrum is assumed to be shaped by a power-law energy distribution of the emitting relativistic electrons $N(E)~{}dE=KE^{-p}~{}dE$ for energy $E$, normalization $K$ and energy index $p\geqslant 2$. The normalization and magnetic field strength $B$ are expressed in terms of the radial distance $r$ from the supermassive black hole as $K=K_{0}(r/r_{0})^{-2}$ and $B=B_{0}(r/r_{0})^{-1}$ with scaling constants $K_{0}$ and $B_{0}$ at a fiducial distance $r_{0}$ (Zdziarski et al., 2015; Mohan et al., 2015; Agarwal et al., 2017). With this, the particle kinetic energy density is

where $\epsilon_{e}$ is the fraction of the total energy density in the particle kinetic energy and $E_{\rm min}$ is the minimum energy required to accelerate electrons (assuming that they are the dominant constituents of the jet) to relativistic energies, here taken to be the rest mass energy $E_{\rm min}=0.51$ MeV. Assuming equipartition of the total energy density between that in the magnetic fields $U_{B}=B^{2}/(8\pi)$ and in the particle kinetic energy $U_{e}$, i.e. $U_{e}/\epsilon_{e}=U_{B}/\epsilon_{B}$, the normalization constant is

where $\epsilon_{B}$ is the fraction of the total energy density in the magnetic fields and $\epsilon=\epsilon_{e}/\epsilon_{B}$. The total energy density is then

The jet luminosity attributable to synchrotron emission in a region of size $\varpi=r\sin\theta_{j}$ (where $\theta_{j}$ is the jet half opening angle) downstream of the jet base is (Ghisellini & Tavecchio, 2010)

where $\beta_{j}$ and $\Gamma_{j}$ are the bulk velocity and Lorentz factor of the jet. The optically thick absorption coefficient (Rybicki & Lightman, 1986; Ghisellini, 2013)

where $e$ is the electron charge, $\nu^{\prime}$ is the emission frequency in the source rest frame, and $f(p)$ is expressed in terms of $p$ and the Gamma function $\Gamma$ dependence as

The optical depth for the emitting region $\varpi$ is $\tau_{\nu^{\prime}}=\alpha_{\nu^{\prime}}\varpi$. From the condition $\tau_{\nu^{\prime}}=1$ that signifies the transition of the emitting region from optically thick to thin, the associated radial distance corresponds to the location of the emitting core. The associated synchrotron self-absorption frequency $\nu^{\prime}=\nu(1+z)/\delta$ where $\nu$ is the frequency in the observer frame, accounting for cosmological redshift through the factor $(1+z)$ for a source at redshift $z$, and relativistic beaming through the Doppler factor $\delta$. Using eqns. E2, E4 and E5 in the above condition,

with the familiar $r\propto\nu^{-1}$ as expected for a self-absorbed radio core (Falcke & Biermann, 1995; Lobanov, 1998).

The mean flux density of the radio core at 15 GHz from the VLBI data is $\left<S_{\nu}(C)\right>=907.15\pm 90.87$ mJy (see Table 2). The associated radio luminosity $L_{R}=4\pi\nu D^{2}_{L}\left<S_{\nu}(C)\right>(1+z)^{-1+\alpha}=2.37\times 10^{42}$ erg s^-1 for $\nu=15$ GHz, $D_{L}=403.1$ Mpc, and $z=0.089$. We use empirical relations to estimate the total jet luminosity (radiative $L_{\rm jet,rad}$ and kinetic $L_{\rm jet,kin}$ components, associated with emission and the acceleration of baryonic constituents of the jet, respectively) from the radio luminosity (Foschini, 2014; An et al., 2020)

The total jet luminosity $L_{\rm jet}=L_{\rm jet,rad}+L_{\rm jet,kin}=1.97\times 10^{44}$ erg s^-1.

Jet properties are estimated from the VLBI 15 GHz data points during the flaring phases ($\delta\geqslant$ 4.3; median Doppler factor). These include an average $\tilde{\delta}=16.9$, an associated $\Gamma=8.5$ based on $\Gamma=(\beta^{2}_{\rm app}+\tilde{\delta}^{2}+1)/(2\tilde{\delta})$. The jet half opening angle $\theta_{j}=1/\Gamma=6\fdg 7$, similar to the lower limit estimated in Chamani et al. (2021). With the estimated $L_{\rm jet}$, $p=2.3$, the assumption of equipartition in the energy density between that in the magnetic fields and particle kinetic energy with $\epsilon=1$ and the above physical properties, the radial distance of the self-absorbed core from eqn. E7, $r_{\rm SSA}\approx 435r_{G}$ (where $r_{G}=GM_{\bullet}/c^{2}=1.47\times 10^{13}$ cm is the gravitational radius for a SMBH of $1.84\times 10^{8}~{}M_{\odot}$ (Grier et al., 2012)). Using this in eqn. E4 ($r_{0}=r_{\rm SSA}$), the associated magnetic field strength is $B_{0}=13.8$ G. For the $B\propto r^{-1}$ scaling relation, this corresponds to $B(r=1~{}{\rm pc})=52.8$ mG. This estimate is consistent with the core-shift based measurement of $\leqslant 60$ mG and of a similar order of magnitude to the SSA-based measurement of $\leqslant 20$ mG (Chamani et al., 2021).

The AGN flares based on generalized shock models can be described by a fast-rising and slowly-declining pattern (Valtaoja et al., 1999), which can be modeled with exponential functions. However, the rising phase of the flux of III Zw 2 does not exhibit a significantly shorter timescale than the declining phase. For example, the rise of the 1998–1999 flare of III Zw 2 is even slower than the decay timescale (Brunthaler et al., 2005). Several lower-amplitude flares were observed prior to the 2009–2010 flare peak, similar to the 1999 and 2004 flares. On the one hand, these smaller flares did not form jets observable in VLBI images, and on the other hand, the opacity of the core and the limited imaging resolution hindered the detection of the corresponding fine structural changes. However, the superposition of these sub-flares leads to a slow increase in flux density before the peak of the major flare. Therefore, it is difficult to obtain a 1:1.3 ratio (as found in many other blazars by Valtaoja et al. 1999) between the rising and declining timescales in III Zw 2. For this reason, we did not fix the decay-to-rise ratio in exponential functions used to model the III Zw 2 flares.

In addition, two $\gamma$-ray flares were found during the 2009-2010 flare, one is associated with the peak of 2009 flare, and the other is in mid of 2010. This motivates us to decompose the light curves with two flare components as an approximate description of the primary 2009 flare and the subsequent 2010 flare (Figure 8), respectively. We need to mention that the model curves shown in Figure 8 are not obtained by least-squares fitting to the observed data, but by selecting model curves with minimum residuals from those obtained with different parameter combinations.

In the main text, we inferred a correlation between $\gamma$-ray and prominent radio bursts, and that the flares are directly connected to the production of jet knots. These can be explained by the shock-in-jet model developed for the blazars (Marscher et al., 2008). Then where did these events occur ?

From Figure 2 we find the jet knots J1 and J2 move along ballistic trajectories. Assuming that the jet speed remains constant, then we can trace back to obtain the position of the jet J2 at the moments of the $\gamma$-ray and radio flares. They are in sequence: 0.041 mas (2009 $\gamma$-ray flare), 0.096 mas (2009 radio flare), 0.163 mas (2010 $\gamma$-ray flare), and 0.205 mas (2010 radio flare). The mass of the SMBH in III Zw 2 is $(1.84\pm 0.27)\times 10^{8}M_{\odot}$ (Grier et al., 2012), therefore its gravitational radius is $R_{g}=9\times 10^{-6}$ pc. The physical size corresponding to the first $\gamma$-ray flare in 2009 is $7.6\times 10^{3}R_{g}$, expressed in units of gravitational radius. This size scale corresponds to the starting section of the jet, and it can be assumed that the collimation of the jet occurs before that distance. A smaller black hole mass of $5\times 10^{7}M_{\odot}$ is derived by Kaastra & de Korte (1988). This implies that jet collimation would occur at a longer distance, measured in units of the black hole’s gravitational radius.

The second $\gamma$-ray flare and subsequent radio flare may happen at a distance of 0.27–0.34 pc, which is consistent with the dimension of the broad line region of III Zw 2  (Kaastra & de Korte, 1988). This suggests that the jet probably collides with the inner boundary of the disk wind within the broad line region, or with the clouds in the torus on a larger scale, resulting in the observed flares (see discussion in Section 3.7).