Interactions between the Jet and Disk Wind in Nearby Radio-intermediate Quasar III
Zw 2
Ailing Wang1,2 , Tao An1,2 , Shaoguang Guo1,2 , Prashanth Mohan1 , Wara Chamani3,4 , Willem A. Baan5,6 ,
Talvikki Hovatta3,7 , Heino Falcke8 , Tim J. Galvin9,10 , Natasha Hurley-Walker10 , Sumit Jaiswal1 ,
Anne Lahteenmaki3,4 , Baoqiang Lao1,11 , Weijia Lv1, Merja Tornikoski3 , and Yingkang Zhang1
1 Shanghai Astronomical Observatory, CAS, Nandan Road 80, Shanghai, 200030, Peopleʼs Republic of China; antao@shao.ac.cn
2 School of Astronomy and Space Sciences, University of Chinese Academy of Sciences, No. 19A Yuquan Road, Beijing 100049, Peopleʼs Republic of China
3 Aalto University Metsähovi Radio Observatory, Metsähovintie 114, FI-02540 Kylmälä, Finland
4 Aalto University Department of Electronics and Nanoengineering, P.O. Box 15500, FI-00076 Aalto, Finland
5 Xinjiang Astronomical Observatory, Chinese Academy of Sciences, 150 Science 1-Street, 830011 Urumqi, Peopleʼs Republic of China
6 Netherlands Institute for Radio Astronomy ASTRON, NL-7991 PD Dwingeloo, The Netherlands
7 Finnish Centre for Astronomy with ESO, University of Turku, Vesilinnantie 5, FI-20014, Finland
8 Department of Astrophysics, Institute for Mathematics, Astrophysics and Particle Physics, Radboud University, Nijmegen, P.O. Box 9010, 6500 GL, Nijmegen, The
Netherlands
9 CSIRO Space and Astronomy, P.O. Box 1130, Bentley, WA 6102, Australia
10 International Centre for Radio Astronomy Research, Curtin University, Bentley, WA 6102, Australia
11 School of Physics and Astronomy, Yunnan University, Kunming, 650091, Peopleʼs Republic of China
Received 2022 November 17; revised 2022 December 25; accepted 2022 December 27; published 2023 February 23
Abstract
Disk winds and jets are ubiquitous in active galactic nuclei (AGN), and how these two components interact
remains an open question. We study the radio properties of the radio-intermediate quasar III Zw 2. We detect two
jet knots, J1 and J2, on parsec scales that move at a mildly apparent superluminal speed of 1.35c. Two γ-ray flares
were detected in III Zw 2 in 2009–2010, corresponding to the primary radio flare in late 2009 and the secondary
radio flare in early 2010. The primary 2009 flare was found to be associated with the ejection of J2. The secondary
2010 flare occurred at a distance of ∼0.3 pc from the central engine, probably resulting from the collision of the jet
with the accretion disk wind. The variability characteristics of III Zw 2 (periodic radio flares, unstable periodicity,
multiple quasiperiodic signals and the possible harmonic relations between them) can be explained by the global
instabilities of the accretion disk. These instabilities originating from the outer part of the warped disk propagate
inward and can lead to modulation of the accretion rate and consequent jet ejection. At the same time, the wobbling
of the outer disk may also lead to oscillations of the boundary between the disk wind and the jet tunnel, resulting in
changes in the jet–wind collision site. Object III Zw 2 is one of the few cases observed with jet–wind interactions,
and the study in this paper is of general interest for gaining insight into the dynamic processes in the nuclear
regions of AGN.
Unified Astronomy Thesaurus concepts: Active galactic nuclei (16); Galaxies (573); Very long baseline
interferometry (1769); Radio jets (1347)
1. Introduction
Disk winds and jets are ubiquitous in active galactic nuclei
(AGN; Yuan & Narayan 2014; Blandford et al. 2019) and play
an important role in the AGN feedback to their host galaxies
(King & Pounds 2015; Harrison et al. 2018; Shi et al. 2021).
Radio-loud (RL) AGN12 comprise about 10% of the entire
AGN population (e.g., Ivezić et al. 2002); their radio emission
is dominated by relativistic jets (Urry & Padovani 1995;
Blandford et al. 2019). The energy release from radio-quiet
(RQ) AGN, which occupy the majority of the AGN population,
is dominated by thermal emission related to the accretion disk
(Padovani 2016; Panessa et al. 2019). Observational evidence
and magnetohydrodynamic models suggest that low-power jets
and winds may coexist in RQ AGN (e.g., Tombesi et al. 2012;
Fukumura et al. 2014; Giroletti et al. 2017). However, whether
and how the jet and disk wind interact with each other remains
an open question (Panessa et al. 2019). Some studies suggest
that there exists a class of objects with moderate radio
loudness, called radio-intermediate (RI) AGN (Miller et al.
1993). Observing RI AGN is much less difficult than RQ AGN,
and RI AGN have mixed properties of RL and RQ AGN, thus
providing an opportunity to study jet- and wind-driven AGN
feedback and jet–wind interactions.
The unusual AGN III Zw 2 (Zwicky 1967; also named
PG 0007+106 and Mrk 1501) contains many enigmatic
observational characteristics. It is hosted by a spiral galaxy
(Hutchings et al. 1982; Hutchings & Campbell 1983) at redshift
z= 0.0893 (Sargent 1970), with a prominent tidal arm (Surace
et al. 2001) to the N of the nucleus but showing spectroscopic
characteristics of a typical type I Seyfert galaxy (Arp 1968;
Osterbrock 1977). It is further included in the bright quasar
sample (Schmidt & Green 1983) with a bolometric luminosity
up to ≈1045 erg s−1 (Piccinotti et al. 1982; Falcke et al. 1995).
In radio bands, III Zw 2 is identified as an RI AGN
(Falcke et al. 1996b) with a radio loudness of ∼200
The Astrophysical Journal, 944:187 (18pp), 2023 February 20 https://doi.org/10.3847/1538-4357/acaf02
© 2023. The Author(s). Published by the American Astronomical Society.
Original content from this work may be used under the terms
of the Creative Commons Attribution 4.0 licence. Any further
distribution of this work must maintain attribution to the author(s) and the title
of the work, journal citation and DOI.
12 The AGN are divided into RL and RQ AGN by the flux density ratio of the
radio (5 GHz) to optical (4400 Å) continuum, the so-called radio-loudness
parameter (R = S5 GHz/S4400 Å; Kellermann et al. 1989). The RL AGN have
R … 10, and RQ AGN have R < 10.
1
(Kellermann et al. 1989, 1994). The source shows a large
extended structure on kiloparsec scales (Unger et al. 1987;
Brunthaler et al. 2005), but its radio emission on parsec scales
shows blazar-like behavior (Falcke et al. 1999; Liao et al.
2016). The Very Large Array (VLA) images of III Zw 2 show a
triple radio structure extending along the northeast–southwest
(NE–SW) direction with a total extent over 36″ (Brunthaler
et al. 2005). The SW lobe is brighter but shorter than the NE
lobe. The recently published images of III Zw 2 observed by
the upgraded Giant Metrewave Radio Telescope (uGMRT) at
685 MHz (Silpa et al. 2020, 2021) show additional faint
structures beyond the SW and NE lobes previously detected by
the VLA, and these outer lobes extend along the direction
perpendicular to the NE–SW jet.
The most attractive feature of III Zw 2 is its extreme
variability: over twentyfold in radio (Wright et al. 1977;
Schnopper et al. 1978; Aller et al. 1985; Terasranta et al. 1992;
Falcke et al. 1999; Brunthaler et al. 2005) and tenfold in the
X-ray band (Kaastra & de Korte 1988; Salvi et al.
2002a, 2002b). It is also highly variable in the optical
(Lebofsky & Rieke 1980; Sembay et al. 1987), infrared
(Lloyd 1984; Clements et al. 1995), and γ-ray (Liao et al. 2016)
bands. Moreover, the large flares of III Zw 2 in multiple bands
are found to be correlated (Salvi et al. 2002a). The variability
characteristics of III Zw 2 are very rare in RI and RQ AGN and
instead are similar to blazars (Aller et al. 2003; Teräsranta et al.
2004; Richards et al. 2011), making it stand out from the RI
and RQ AGN populations. In addition, the radio flares seem to
exhibit quasiperiodic variations with a period of about 4–5 yr
(Brunthaler et al. 2003; Li et al. 2010), which warrants an in-
depth study of the physical mechanisms underlying the
periodic variability.
Studying the structural changes and variability of the jet
during prominent flare phases can help reveal the mechanism of
jet production and evolution. Imaging the newly born jet
features requires subparsec resolution, which is currently only
possible with very long baseline interferometry (VLBI)
observations. In the past four decades, several large radio
flares have been found in III Zw 2 (Aller et al. 1985; Terasranta
et al. 1992; Brunthaler et al. 2005; Li et al. 2010), with a flux
density approaching or exceeding 3 Jy at 15 GHz. The VLBI
observations during the 1998 flare revealed a compact core–jet
structure within 0.1–0.4 pc (Falcke et al. 1996a, 1999;
Kellermann et al. 1998; Brunthaler et al. 2000). A dramatic
structural change was found from the 43 GHz VLBI data
between 1998 December 12 and 1999 July 15, and an apparent
superluminal speed of 1.25c (Brunthaler et al. 2000) was
inferred, making III Zw 2 the first RI AGN detected with
apparent superluminal jet motion in a spiral galactic nucleus
(Brunthaler et al. 2000, 2005). After the major flare in 2009, the
variability amplitude became smaller (see Figure 4 in
Brunthaler et al. 2005). In the latest VLBI observations of
III Zw 2 in 2017, only one compact core was detected, and the
core was significantly fainter than ∼20 yr ago (Chamani et al.
2021).
In this paper, we investigate the jet kinematics and radio
variability and propose a model to explain the quasiperiodic
variability, the correlation between the prominent flare and jet
production, and the secondary flare created by the jet–wind
collision.
2. Data
The discovery of the southern extended feature from
uGMRT images (Silpa et al. 2020, 2021) motivated us to use
low-frequency interferometric data, including the data from the
Australian Square Kilometre Array Pathfinder (ASKAP; Hotan
et al. 2021) at 888 MHz, GMRT (Swarup et al. 1991) at 150
MHz (the TIFR GMRT Sky Survey; Intema et al. 2017), and
the Murchison Widefield Array (MWA; Tingay et al. 2013)
observations at 72–231 MHz (Appendix A), to study the
complete radio structure, reveal the kiloparsec-scale morph-
ology, and constrain the radio spectrum of the extended
structure. The VLBI data used for the jet kinematics analysis
are obtained from the Monitoring Of Jets in Active galactic
nuclei with VLBA Experiments (MOJAVE; Lister et al. 2018)
program and the archived data from the Astrogeo database.13
Details of the VLBI data reduction are given in Appendix B.
The radio light-curve data are obtained from the single-dish
monitoring programs of the Owens Valley Radio Observatory
(OVRO) at 15 GHz and the Metsähovi Radio Observatory
(MRO) at 37 and 22 GHz (Appendix C).
3. Results and Discussion
3.1. The Jet Structure
The VLA images of III Zw 2 in the literature show a triple
structure along the NE–SW direction (Brunthaler et al. 2005).
The central core C dominates the total flux density at almost all
wavelengths. The SW lobe is located at 15 4 (∼25.5 kpc;
Unger et al. 1987; Kukula et al. 1998; Falcke et al. 1999) and
connected to the central core by a curved jet (Silpa et al. 2020);
the jet initially points to the W and then gradually bends to the
SW. A weaker lobe is detected 21 9 (∼36 kpc) NE of the core
(Brunthaler et al. 2005).
The VLBI image in Figure 1 shows that the parsec-scale jet
extends to the W and is roughly aligned with the SW lobe
(Figure A1); therefore, it is possible that the SW jet is at the
advancing side. If the SW and NE lobes were formed from the
Figure 1. A 15 GHz VLBA image of III Zw 2. The image is created with
natural weighting. The negative red contour is at the −2σ level, and the
positive white contours are in the series of 3σ × (1, 2, 4, 8, 16, 32, 64, 128,
256, 512), where the rms noise is σ = 0.15 mJy beam−1. The color scale shows
the intensity in the logarithmic scale.
13 http://astrogeo.org/
2
The Astrophysical Journal, 944:187 (18pp), 2023 February 20 Wang et al.
AGN activity in the same episode, then the length of the
advancing jet/lobe should be longer than the reverse jet/lobe.
However, the actual situation is the opposite. Although the SW
lobe is brighter than the NE lobe, both the VLA and ASKAP
images (Figure A1) show that the SW lobe is shorter than the
NE lobe. The difference in length between the two lobes seems
to indicate that the ambient interstellar medium (ISM) on two
opposite sides of the nucleus has different properties
(Brunthaler et al. 2005; Silpa et al. 2021); i.e., the ISM at the
SW side has a higher density, and the growth of the SW jet is
more obstructed by the ISM. A companion galaxy (Surace et al.
2001) exists ∼30″ S of III Zw 2 (Figure A2), and gravitational
interactions may have accumulated a large amount of gas in the
intermediate region between the two galaxies. The SW jet
bends southward at the SW lobe, where enhanced polarization
is observed (Silpa et al. 2021), providing observational
evidence for strong interactions between the SW jet and
the ISM.
The VLBI image in Figure 1 shows a compact jet with a
maximum extent of 1.2 mas (corresponding to a projected
distance of ∼2 pc). The integrated flux densities obtained from
the VLBI images match those obtained from single-dish
observations in close epochs, indicating that the contribution
of the optically thin extended jet to the total flux density at
gigahertz frequencies is very small. Another intrinsic reason for
the compact jet is that the radio activity of III Zw 2 is
intermittent, and the VLBI jet is short-lived in nature (see
discussion in Sections 3.2 and 3.4), which also leads to a short
VLBI jet. A similar situation can be seen in another RI AGN,
Mrk 231 (Reynolds et al. 2017; Wang et al. 2021).
Combining images of III Zw 2 observed with different
resolutions at different frequencies, we find that the jet exhibits
an overall S shape. There are many physical mechanisms that
can create S- or Z-shaped jet morphology, including backflow
from hot spots (Leahy & Williams 1984), buoyancy force
bending the outer radio structure in the direction of decreasing
external gas pressure (Worrall et al. 1995), precession of the jet
beam (Ekers et al. 1978), or jet reorientation due to
hydrodynamic processes associated with galaxy mergers
(Gopal-Krishna et al. 2003). Object III Zw 2 lives in a cluster
environment and exhibits ongoing galaxy merger or galaxy–
galaxy interactions (Surace et al. 2001), with observational
evidence of an optical tidal arm N of the III Zw 2 nucleus and a
companion galaxy ∼30″ S of the nucleus (see Figure A2). The
spin flip of the central engine following the galaxy merger can
lead to a reorientation of the jet, which can produce the S- or
X-shaped jet (Gopal-Krishna et al. 2003; Bogdanović et al.
2007). Using a typical advance speed of the terminal hot spot
(e.g., M87, v= 0.11c; Biretta et al. 1995) as a reference, we
obtain a kinematic age of about ∼106 yr for the extended
III Zw 2 lobe. The III Zw 2 jet is not as powerful as the M87 jet
and may have a lower hot-spot advance speed, which would
result in an estimated kinematic age of more than 1Myr but
should not be larger than 107 yr. This kinematic age is within the
timescale of the black hole coalescence (Ebisuzaki et al. 1991)
and the typical lifetime of an AGN, allowing for jet reorientation
as a possible mechanism for the S-shaped morphology.
3.2. Jet Properties at Parsec Scales
The 15 GHz MOJAVE archive data of III Zw 2 covers a time
span of up to ∼18 yr from 1995 July to 2013 June containing
25 epochs. The submilliarcsecond resolution, high sensitivity,
and homogeneous image quality of the MOJAVE VLBI images
make them ideally suited for studying jet kinematics. Although
the jet is very compact, it can be well distinguished from the
core in the Very Long Baseline Array (VLBA) images.
Figure 2(a) shows the variation of the jet position angle with
time, and Figure 2(b) shows the variation of the core–jet
separation with time. It is clear that these jet knots belong to
two separate components (labeled J1 and J2), rather than a
single component. Components J1 and J2 follow their
respective ballistic trajectories at different position angles.
We only used the jet components with the highest
confidence, and the data in the other epochs were discarded
because the data quality was too poor to obtain a reliable model
fit (Appendix B). We obtain proper-motion velocities of
v(J1)= 1.35± 0.13c and v(J2)= 1.21± 0.07c, consistent with
those obtained by the MOJAVE team using all fitted
components (Lister et al. 2019), and also in good agreement
with previous studies based on the 43 GHz VLBI observations
(Brunthaler et al. 2000).
Assuming that the velocity of J1 remains constant from the
time it was created until it disappeared, the ejection time of J1
can be traced back to epoch 2004.74 and is temporally
correlated with the 2004 flare. Similarly, the ejection of J2 is
associated with the 2009 flare. From 1996 onward, three major
radio flares over 3 Jy were observed at 37 GHz, with their
peaks in early 1999 (Brunthaler et al. 2005), late 2004
(Chamani et al. 2020), and 2009 (Figure 3 in the present paper).
Each flare is associated with the creation of a discrete jet knot
(the first discovered superluminal jet knot in Brunthaler et al.
2000 and J1 and J2 reported in this paper), possibly connected
with intermittently enhanced accretion (see more discussion in
Section 3.5).
Some lower-level flares that occurred during the intervals
between these major flares, such as those in 2003 and 2006, did
not produce identifiable long-lived jet knots. Although the lack
of long-term follow-up monitoring of the 1998 jet prevents us
from obtaining information on how it evolved over time, our
observations clearly indicate that the proper-motion speeds of
J1 and J2 are almost the same as that of the 1998 jet, suggesting
that the fundamental condition for the jet production has not
changed significantly during these jet creation events.
The core brightness temperatures (column (7) of Table B2)
calculated from the VLBI-based measurements (core flux
density and size) all exceed the equipartition brightness
temperature limit (Readhead 1994). The Doppler boosting
factor can be estimated as δ= Tb,obs/Tb,eq, and the Lorentz
factor Γ can be estimated by 1 22 app
2( ) ( )d b dG = + + . We
find that large Γ values are related to the major flares. In the
quiescent state, we get a mean value of δ= 2.6. Taking into
account the jet proper-motion speed βapp= 1.35c yields a
viewing angle of approximately 20°, which is a typical value
for a nonblazar radio quasar (Ghisellini et al. 1993). A larger
viewing angle of 35°.4 was derived by Hovatta et al. (2009),
and an upper limit of 41° was derived by Brunthaler et al.
(2000) due to the different Doppler factor and jet speed used.
We calculated the magnetic field strength at 1 pc to be
∼53 mG (Appendix E), similar to that estimated based on the
apparent core shift and synchrotron self-absorption (SSA;
Chamani et al. 2021). This consistency supports the low jet
magnetic flux in III Zw 2 and the inference that the central
engine did not reach the magnetically arrested disk state
(Chamani et al. 2021). Interestingly, our estimate is based on
3
The Astrophysical Journal, 944:187 (18pp), 2023 February 20 Wang et al.
data from the flaring phases, while the earlier estimates by
Chamani et al. are based on observations in a quiescent phase,
suggesting that the jet structure remains relatively unperturbed
by any fresh injection events at the jet base (at a distance of
∼800 gravitational radii).
3.3. Radio Spectrum
We plot the radio spectrum of III Zw 2 from 72 MHz to
37 GHz in Figure 4. We have chosen data points as close as
possible to the time of the first MWA observation epoch (i.e.,
2018 June 15). Due to the lack of low-frequency data in
previous studies, the spectrum below 685 MHz has not been
well constrained. The spectrum components below and above
685 MHz have distinctly different origins, so we fit the whole
spectrum with two radiation components. At ν> 685MHz, the
core dominates the emission, showing an inverted spectrum (or
peaked spectrum); at ν< 685MHz, the flux density of the core
is substantially reduced due to increasing SSA toward low
frequencies, and the extended jets and lobes become gradually
dominant.
We first fit the NE and SW spectra with power-law functions
(i.e., S∝ να) based on the VLA data from the literature
(Brunthaler et al. 2005) and obtained their spectral indices as
αNE= −0.79 and αSW= −0.72, respectively. These are typical
values for optically thin synchrotron emission. We then
extrapolated the flux densities of NE and SW to the MWA
band (central frequencies of 88, 118, 185, and 216 MHz). Next,
the flux densities of the core were fitted with a self-absorbed
synchrotron spectrum model (Pacholczyk 1970; Türler et al.
1999) by fixing the optically thick spectral index αthick to 2.5
(Appendix D). The fit yields a low-frequency power-law
spectrum with the amplitude A 27 9
12= -+ mJy, low-frequency
spectral index 0.97low 0.21
0.21a = - -+ , and baseline flux density
B 41 12
8= -+ mJy and a higher-frequency optically thick
spectrum with the amplitude F 146m 9
11= -+ mJy, turnover
frequency 11.24m 1.25
2.01n = -+ GHz, and optically thin spectral
Figure 2. Properties of VLBI: (a) changes of jet position angle with time and (b) jet proper motion. The blue dashed lines represent the linear regression fit whose
slope gives the jet proper motion; red and green dotted lines mark the time of the γ-ray and radio flare peaks at 15 GHz, respectively. (c) Changes of the VLBI
components’ size with time. (d) Changes of the integrated flux densities of the VLBI components with time. To make the trend of the integrated flux density of the jets
more visible, the integrated flux density of the jets has been multiplied by a factor of 2. The fitted parameters for the VLBI components are given in Table B2. The
linear fitting results are extrapolated to obtainthe jet ejection time, 2004 September and 2009 August for J1 and J2, respectively.
4
The Astrophysical Journal, 944:187 (18pp), 2023 February 20 Wang et al.
index 0.20 0.12
0.10a = - -+ . However, our spectrum shown in
Figure 4, constructed with data points measured when the
source was in a low-activity state (indicated by the red arrow in
Figure 3), is different from that presented in Falcke et al.
(1999), which was in the flaring state. During the 1998 flare,
the core had an inverted spectrum between 1.4 and 666 GHz
peaking around 43 GHz (Falcke et al. 1999), a factor of 3
higher than our fit. In the quiescent state in 2018 depicted by
Figure 4, the emission at gigahertz frequencies is still
dominated by the core C, but the turnover frequency has
shifted to 11.24 1.25
2.01-+ GHz.
We then subtracted the extrapolated flux densities of C, NE,
and SW from the observed MWA flux densities, and the
remaining flux density is mainly from the extended compo-
nents N+S. Finally, we fitted the N+S flux densities with a
power-law function and obtained a spectral index of α=
−1.09± 0.12, which is much steeper than that of the inner
lobes NE and SW but consistent with the spectral indices of
radio relics (Shulevski et al. 2015; Quici et al. 2021).
3.4. Temporal Evolution of the VLBI Components
Figures 2(c) and (d) show the correlation between the core
size and the flux density. The core size gradually increased
when the flux density was in the rising phase before the flare
peak; after the flare peak, the core size gradually decreased.
The value of Spearman’s rank correlation coefficient
(Lehman 2005; Dodge 2008; Myers et al. 2013) between the
core size and flux density is 0.43 with a p-value of 0.14. Their
positive correlation needs to be verified by more data. If this
correlation is further confirmed, this phenomenon may be
naturally explained by the superimposition of a flaring and fast-
moving jet component on a quiescent and stationary core. The
resolution of the VLBI images in this paper is not yet sufficient
to distinguish between these two components in the initial
flaring phase. Only when the jet knot produced after the flare
moves to a certain distance is it sufficient to separate the jet
from the core clearly.
Figure 3. Radio light curves of III Zw 2. The 15 GHz data are observed with the 40 m telescope of OVRO. The 37 GHz data are observed with the 14 m radio
telescope of MRO. The MRO data from 1986 to 2019 have been published in Chamani et al. (2020) and other previous studies; in addition to these data, we include
the new data from 2019 onward to date. The VLBI data are from the MOJAVE program (Lister et al. 2018). The arrow marks the time of the MWA observation
(Hurley-Walker et al. 2022). The red vertical shaded areas mark the time of γ-ray flare occurrence, and the dashed lines correspond to the peaks of the two γ-ray flares
(Liao et al. 2016).
Figure 4. Radio spectrum of III Zw 2. The red filled squares are taken from the
MWA GLEAM-X survey. The green filled circle represents the GMRT 150 MHz
data point (Intema et al. 2017). The magenta data point is the total flux density
from the uGMRT observation (Silpa et al. 2020) on 2018 November 23. The slate
blue data point is the total flux density from the ASKAP observation (McConnell
et al. 2020) on 2020 May 3. The blue, navy, and purple data points are from this
paper, respectively, 37 GHz Metsähovi, 15 GHz OVRO, and 8.4 and 2.3 GHz
Astrogeo VLBI. For these data, we select the epochs close to the MWA and
uGMRT observations in 2018. Core-dominated flux densities show an inverted
spectrum with the best-fit self-absorbed synchrotron spectrum, and more details are
presented in Section 3.3 and Appendix D. The red open squares and blue open
circles are the flux densities of the SW and NE lobes derived from the VLA
images (Brunthaler et al. 2005).
5
The Astrophysical Journal, 944:187 (18pp), 2023 February 20 Wang et al.
As the flaring component (a propagating shock) passed
through the core, it led to a gradual increase in the flux density
of the observed core. As the duration of the radio flare (lasting
several years) was much longer than the cooling time of the
synchrotron radiation of relativistic electrons (∼130 days at 15
GHz), it required a continuous injection of fresh relativistic
electrons into the core to ensure that the core flux density
continued to increase for about 2 yr. The flaring component
continued to move outward and gradually separated from the
stationary core. For a period after the flare (about 1 yr), the
VLBI image resolution was not sufficient to distinguish
between the core and the flaring jet component, but the core
size was observed to increase as the jet moved outward
(Figure 2(c)). With the cessation of the intermittent energy
injection, the core returned to the quiescent state until the next
flaring component passed by. At the same time, the size and
brightness of the core gradually became smaller as the flaring
jet component continued to move outward and separated from
the stationary core. On the other hand, the flaring component
became optically thin, and its surface brightness declined over
time. An alternative model of an inflating balloon has been
proposed to explain the spectrum and structure evolution of
III Zw 2 by Falcke et al. (1999). In their model, the initial phase
of the flare is interpreted as the relativistic jet interacting with
the torus, and the jet–ISM collision causes frustration with the
jet’s advancing motion. Their model is compatible with the one
we propose here to explain the evolution of major flares.
Moreover, their jet–ISM explanation is also consistent with the
jet–wind collision model we propose in Section 3.7 to interpret
the secondary 2010 flare. High-cadence, high-resolution VLBI
monitoring of the flaring jet helps to refine this physical picture.
The flux density variability of the jet in Figure 2(c) shows a
delay relative to the core, although this characteristic is not very
pronounced due to the sparse distribution of the data points. On
the other hand, the temporal evolution of the jet component size
(see Table B2) displays an inverse correlation with the core
size, as is most evident after 2009 (corresponding to J2). We
suggest that after separation from the core region, the moving
discrete jet component underwent adiabatic expansion, leading
to an increase in size, as well as a decrease in surface
brightness. The jet disappeared until the brightness of the shock
fell below the detection threshold.
To summarize the properties of the jet of III Zw 2 described
above, we find that the radio properties of III Zw 2 in the
quiescent state make it look more like a Seyfert 1 galaxy.
However, its strong variability, apparent superluminal jet
motion, and large extended structure are very distinct from
normal RQ quasars, instead resembling RL quasars. Specifi-
cally, its radio properties during the major flares, such as
superluminal jet motion and high brightness temperature, are
consistent with those of a blazar (Falcke et al. 1999). This
hybrid feature, being an RQ quasar at the quiescent state but
behaving like a blazar at the flaring state, has been noted in
previous studies (e.g., Falcke et al. 1999) and might be a
common feature of RI AGN (Section 3.8).
3.5. Periodic Flares and Jet Ejection
Figure 3 shows the 37 GHz light curve with two particularly
major flares in the time period 2003–2020, peaking in late 2004
and late 2009. The 2004 flare has been reported by Li et al.
(2010). The 2009 flare is the highest in magnitude in ≈40 yr of
radio monitoring, with a maximum flux density of 3.12 Jy,
approximately 20 times the baseline in the quiescent state. Such
a large radio variability is rare even in blazars’ light curves
(Richards et al. 2011). After the 2009 flare, the source became
less active, and the flare peak gradually decreased, although
there were a few lower-amplitude flares in 2013, 2015, and
2017. Since 2018.4, the source has entered a quiescent state.
Similar to the 37 GHz light curve, the 15 GHz light curve also
displays a number of flares, with the largest one in late 2009–
early 2010. The maximum flux density increased by a factor of
10.5 compared to the baseline level. After 2016, the source
entered a low-level state.
Brunthaler et al. (2003) found that III Zw 2 has major radio
flares about every 5 yr. Li et al. (2010) found the same flare
period based on the historical light curves of III Zw 2 at 22 and
37 GHz obtained from the MRO database (Teräsranta et al.
2005) covering 18 yr from 1986 to 2004. We continue to
analyze the periodicity of the 37 GHz light curve from 2011.40
and the 15 GHz light curve after 2008 using the Lomb–Scargle
periodogram, yielding a period of P∼ 2.1 yr (Appendix C).
This periodic signal is not sharply distributed in the Lomb–
Scargle periodograms; therefore, it can only be called a
quasiperiodic signal. In the wavelet periodograms, the most
prominent feature is around 1.97–2.02 yr (Figure 5), and a
secondary weaker feature is around 3.97–4.26 yr. The period of
P∼ 4 yr has been steadily maintained for ∼35 yr (Brunthaler
et al. 2005; Li et al. 2010; the present paper) with six to seven
complete cycles, suggesting that the periodicity should not be a
fake signal caused by random red noise but is related to some
intrinsic dynamical process. The 2 yr periodicity has also
appeared in previous periodicity analysis (Li et al. 2010), but
the magnitude is relatively low.
Quasiperiodic variations in AGN light curves are a common
observational phenomenon and often interpreted as being
related to regular perturbations at the jet origin, such as the
precession of the jet nozzle (Valtonen et al. 2008) and rotation
of a hot spot along a helical path (Camenzind &
Krockenberger 1992; Mohan et al. 2016a) or magnetohydro-
dynamic or hydrodynamic instabilities arising from the starting
section of the jet (Hardee 2003; Jorstad et al. 2022) or disk (An
et al. 2013; Wang et al. 2014) and propagating as helical-mode
waves.
In many cases, warping of the accretion disk can occur. In
the case of jet precession, if the spin axis of the primary black
hole is not aligned with the orbital plane of the binary, the
differential precession with the change in radius can cause the
warping of the outer disk (Bardeen & Petterson 1975).
Alternatively, the self-irradiation of a luminous accretion disk
can also lead to the disk warping out of the orbital plane
(Pringle 1996). The dynamic process associated with the
warped disk generally occurs in the outer disk region, while the
jet is launched from the inner region of the disk (Blandford &
Payne 1982) or the vicinity of the black hole (Blandford &
Znajek 1977). There is a big gap between the two physical
processes of jet launching and disk warping. The coupling of
the jet and disk would require that the instabilities originating
from the edge or the outer region of the disk could be
transmitted to the inner region, whereas this inward propaga-
tion is difficult to achieve through viscous processes because
the required timescale is too long. Instead, the year-timescale
periodic variability of III Zw 2 could result from global
acoustic or p-mode oscillations in a thick disk (Rezzolla
et al. 2003a, 2003b). The inferred variability period is on the
6
The Astrophysical Journal, 944:187 (18pp), 2023 February 20 Wang et al.
order of several years for a black hole mass of 108 Me (An
et al. 2013). The trigger of the p-mode oscillations can be a
locally periodic agent at the edge of the disk due to either
gravitational or radiation perturbations. The inward propaga-
tion of global oscillations to the inner region induces
quasiperiodic fluctuations in the accretion flow, which in turn
trigger quasiperiodic injection of plasma into the jet, leading to
the observed harmonic periodic radio variability and the
corresponding periodic jet ejection.
In addition to the primary 4 yr period, there is a secondary
peak at P≈ 2 yr in the periodograms, consistent with four
lower-amplitude flares in 2013, 2015, 2017, and 2019
(Figure 3). The 2 yr periodic component could be a harmonic
of the 4 yr periodicity, with each being dominant in different
activity states. Only the fundamental frequency of the
oscillations (i.e., f= 0.25 yr−1) can produce the largest flare,
as well as the long-lasting jet knot that can be observed in
VLBI images. Lower-amplitude flares associated with harmo-
nic components may also generate jet knots, but they are too
weak to be detected. The presence of multiple harmonic
periodicities on year timescales can be explained by hydro-
dynamic instability, such as the global p-mode oscillations of
the accretion disk. Moreover, the quasiperiodic signals induced
by the accretion disk instability have no fixed frequency and
often vary within a range, which is consistent with the
observations.
3.6. Correlation between γ-Ray and Radio Flares
There are correlations between the prominent variabilities of
III Zw 2 in multiple bands (see Introduction). Liao et al. (2016)
identified two γ-ray flares that peaked in 2009 November and
2010 May, respectively. The properties of the III Zw 2 γ-ray
flares are analogous to those of blazar flares (Liao et al. 2016),
implying that the central engines between III Zw 2 and blazars
are not essentially different. The radio light curves (Figure 3)
also have substructures during the period 2009–2011; in
addition to the strongest flare peaking in late 2009 (the 2009
flare), there is also a lower-amplitude flare on the shoulder of
the declining phase that peaks mid-2010 (the 2010 flare).
The 2009–2011 light curves could be better described with
two flare components than with just one (Appendix F). The
observed light curve matches well with a simulated light
curve consisting of two “exponential rise + exponential
decay” flare components as an approximate description of the
composite behavior of the flares (Figure F1). At both
frequencies, the primary 2009 flare is brighter than the
secondary 2010 flare. The 2009 flare peaks at epoch 2009.96
at 37 GHz, leading the peak epoch 2010.09 at 15 GHz by
∼47 days. The 2010 flare peaks at epoch 2010.56 at 37 GHz
and epoch 2010.58 at 15 GHz, with no significant time delay.
The fact that the 37 GHz radio flare leads the 15 GHz flare
can be naturally explained by the frequency-dependent
Figure 5. Periodogram of III Zw 2. The top panels show the periodograms of 37 (left) and 15 (right) GHz light curves. In the periodogram plots, the horizontal dashed
line is the level of white noise (constant). The curved line represents the mean periodogram (closest to the underlying PSD) from the Monte Carlo simulations using
the best-fit model parameters; this additionally accounts for model uncertainties. The dashed curve is the 95% confidence level (encompasses 95% of the periodogram
ordinates; any outliers are statistically significant quasiperiodic oscillations). A broad peak stands out in both plots, corresponding to a period of ∼2 yr. The bottom
panels show the results from wavelet analysis of the light curves from 2011.6 to 2021: left, 37 GHz; right, 15 GHz. A persistent component corresponding to a
characteristic timescale of 1.97–2.02 yr is clearly seen throughout the time span. A secondary feature appears at a timescale of 3.97–4.26 yr.
7
The Astrophysical Journal, 944:187 (18pp), 2023 February 20 Wang et al.
opacity (Hovatta et al. 2008). The time span between the
2009 and 2010 flares of 0.5–0.6 yr is much longer than the
lifetime of synchrotron electrons (i.e., the cooling time of 85
and 135 days at 37 and 15 GHz) but much shorter than the
time separation between two major flares (∼4 yr). Therefore,
the 2009 and 2010 flares must be associated with the same
episodic nuclear activity but are unlikely to be the same
emission components.
The 2009 flare is reminiscent of the typical core flares
occurring in blazars. In γ-ray AGN, both the γ-ray flare and the
associated delayed millimeter-wavelength flare are created in
the standing shock in the innermost jet (the shock-in-jet model;
Marscher et al. 2008), which is optically thick and shows time
delays in the light curves at different radio frequencies. The
2009 radio flare lagged the γ-ray peak (2009.84) by only ∼40
days (Liao et al. 2016), suggesting that the γ-ray emitting zone
is spatially connected to the radio-emitting zone, and the γ-ray
emitting site is closer to the central engine than the radio flare
site. In contrast, the concurrent 2010 flare at 37 and 15 GHz
frequencies suggests that it should arise from an optically thin
jet component, probably associated with a compressed shock in
the jet propagating downstream. In addition, the time-
integrated flux density is also higher at 37 GHz than at
15 GHz, implying that the energy output of the radio source is
concentrated at higher frequencies during the most active phase
of the flare. However, the peak time and maximum amplitude
of the 2010 flare are similar at both 37 and 15 GHz, suggesting
that the energy dissipation of the 2010 flare is weakly
dependent on frequency, reinforcing the intrinsic difference
between the 2009 and 2010 flares.
3.7. Jet–Wind Collision and the Associated 2010 Flare
The pieces of observational evidence presented in previous
sections, including the temporal evolution of the flux density
and jet component size, optically thin radio spectrum, and
correlated radio/γ-ray flares, lead us to speculate that the 2010
flare could be a compressed shock caused by the downstream
jet–ISM collision. Additional evidence for jet–ISM collision
comes from the VLBA MOJAVE polarimetric observations;
the fractional linear polarization of III Zw 2 increased from
0.2% on 2009 June 3 (prior to the 2019 flare) to 0.7% on 2010
July 12 (corresponding to the peak time of the 2010 flare), and
the polarization angle changed from 6° to 24°. The enhanced
linear polarization and change of polarization angle in the jet
offer direct evidence of a compressed shock created in the jet–
ISM collision, as has been found in other RL quasars (An et al.
2020).
In the warped disk model discussed in Section 3.5, the wind
or outflow driven by the radiation pressure of the disk is
released from the outer region of the disk to form a cylindrical
or conical structure in the broad-line region (BLR). The jet
flushes a tunnel in the axial direction of the wind, in which the
jet moves outward. The axis of the disk wind is bound to the
axis of the outer disk, and the oscillations of the outer disk
would lead to the wind wall oscillating within a certain angle as
well. Thus, the jet axis is not always parallel to the wind axis; at
a certain distance, the jet may hit the inner boundary of the
wind wall, producing an oblique shock. On the jet–wind
interaction interface, the oblique shock is deflected, after which
the shock (jet knot) follows a (new) ballistic trajectory in the
direction of the deflection. The loss of the jet kinetic energy
during this collision leads to the dissipation of radiation,
producing the observed prominent γ-ray and radio flares (e.g.,
the 2010 flare). On the other hand, the oscillations of the wind
boundary would cause the working surface of the jet–wind
collision to vary with time, which seems to coincide with the
observed change in the position angle leading to a change in
the position angle of the (deflected and redirected) VLBI jet, as
seen in Figure 2.
Recent studies (Boccardi et al. 2016, 2021) have shown that
the powerful nuclei of high-excitation galaxies produce disk-
launched winds/outflows that could form the slower jet
sheaths, in addition to highly relativistic jet spines. The slower
sheath may contribute to the collimation of the jet. This is
compatible with the model we propose; in our model, it is the
accretion disk wind rather than the outer-layer jet that
surrounds the jet spine. Collision may occur between the
fast-moving jet spine and the inner boundary of the disk wind
when the axis of the winds is not aligned with the jet axis.
Where did the flares occur? Extrapolating the trajectory of
the jet component J2 back to the peak times of the 2009 and
2010 radio flares, we obtain a distance of J2 of 0.096 mas (a
projected distance of 0.16 pc) and 0.205 mas (0.34 pc) from the
black hole, respectively. Neither of these size scales is
resolvable by the current 15 GHz VLBI observations, and
even the previous 43 GHz VLBI observations (resolution of
0.1 mas; Brunthaler et al. 2005) can only barely resolve this
size. Therefore, all relevant jet–wind collision and flaring
activities are hidden in the 15 GHz VLBI core. Using the same
method, we extrapolate to estimate the distance of the 2009 γ-
ray flare site from the central black hole to be 0.068 pc, a
distance that can be considered as the upper limit of the jet
collimation zone. That is to say, the jet collimation must have
been complete at this distance (7.6× 103Rg in projection).
This distance is comparable with that of the jet collimation
zone derived from other AGN, such as M87 (Asada &
Nakamura 2012; Hada et al. 2013) and other nearby radio
galaxies (Boccardi et al. 2021). The 2010 γ-ray flare site,
0.27 pc from the central engine, favors that the jet–wind
collision occurs farther downstream in the jet.
3.8. Comparison with Other Similar Sources
In this section, we compare the radiation and physical
properties of III Zw 2 with Mrk 231 and explore the generalized
properties of RI AGN. One of the closest known RQ quasars,
Mrk 231 has an extremely strong infrared luminosity and rich
multiphase, multiscale outflows (Wang et al. 2021 and
references therein). Object III Zw 2 shares similar observational
features with Mrk 231: ongoing galaxy mergers, high lumin-
osity, misalignment between the parsec- and kiloparsec-scale
jets, prominent flares, and the associated intermittent jet
ejection. In addition, both sources behave like RQ AGN in
the quiescent state and like a blazar in the flaring state (Wang
et al. 2023). Their flare properties and jet kinematics can be
explained by the classical SSA radiation model. The jet knots
follow a ballistic trajectory on few-parsec scales, but the
position angle varies from one knot to another. These
variability properties and jet structure changes are reflective
of interactions between the jet and the ISM in the BLR or the
jet being reflected by a rotating torus (Wang et al. 2021).
Object III Zw 2 differs from Mrk 231 in that the III Zw 2 jet
extends to a larger distance, while the Mrk 231 jet is confined
within the host galaxy. The ability to develop large-scale jets
depends not only on the initial kinetic properties of the jet but
8
The Astrophysical Journal, 944:187 (18pp), 2023 February 20 Wang et al.
also on the external environment (whether or not it chokes the
jet), and the ability of the central engine to remain active for a
long period has a profound effect on the jet growth (An &
Baan 2012).
4. Summary
In this paper, we analyze in detail the radio structure and
variability properties of III Zw 2, an RI AGN. The main results
are summarized as follows.
1. The overall jet structure from parsec to kiloparsec scales
shows an S-shaped morphology, probably related to the
jet reorientation due to galaxy interaction. The low-
frequency ASKAP and MWA images confirm the
presence of extended emission 27″ to the N and 26″ to
the S from the core. The ultrasteep spectra of these
extended features suggest that they are relics of past AGN
activity.
2. Two jet components, J1 and J2, are detected in the VLBI
images with an apparent superluminal velocity of 1.35c
and an average jet viewing angle of ∼20°.
3. The radio light curves show quasiperiodic flares; before
2008, an ∼4 yr cycle dominates, and after 2008, when the
source is in a low-activity state, a high-frequency
harmonic component of ∼2 yr period becomes dominant.
The variability characteristics (quasiperiodicity, two
periodic signals, and the presence of a harmonic relation
between them) can be explained by the global acoustic
oscillations of the accretion disk. The perturbations
occurring in the outer region of the disk propagate
inward, leading to modulated changes in the accretion
rate and, consequently, the generation of periodic radio
flares and jet ejection.
4. Only major flares associated with the fundamental
frequency oscillation can produce observable jet compo-
nents. The two strongest flares occurring in late 2004 and
late 2009 coincide with the creation of the jet knots J1
and J2 observed in the VLBI images, respectively.
5. The radio flare from late 2009 to mid-2010 can be
decomposed into two subflares, corresponding to two γ-
ray flares, respectively. The 2009 γ-ray flare led the radio
flare, and the high-frequency radio flare led the low-
frequency flare, suggesting that it originated from an
optically thick component, probably in the jet collimation
region. The 2010 radio flare occurred simultaneously at
37 and 15 GHz, suggesting that they occurred in an
optically thin jet zone.
6. The wind or outflow arising from the outer accretion disk
forms a cylinder or cone in the nuclear region. As the axis
of the warped disk is misaligned with the jet, at a certain
distance, the jet flow hits the wind wall, creating an
oblique shock that deflects the jet; at the same time, the
jet–wind collision leads to the production of γ-ray and
radio flares (e.g., those observed in 2010).
Object III Zw 2 has hybrid nature of RQ AGN and blazar; in
the quiescent state, it is a typical RQ AGN, while in the flaring
state, it behaves as a blazar. During the intermittent flares, the
produced jet knots interact with the accretion disk wind in the
BLR, producing γ-ray and radio flares. The characteristics
observed from III Zw 2 may be common to the RI AGN
population, in which both the jets and winds coexist and play
important roles at different spatial scales and timescales.
Detailed studies of typical individual RI AGN will improve our
understanding of the RQ/RL AGN dichotomy and the structure
and dynamics of the AGN nuclear region.
This work was supported by resources provided by the
China SKA Regional Centre prototype (An et al. 2019, 2022)
funded by the Ministry of Science and Technology of China
(MOST; 2018YFA0404603). This research has been supported
by the National SKA Program of China (2022SKA0120102).
S.G. is supported by the CAS Youth Innovation Promotion
Association (2021258). The authors acknowledge the use of the
Astrogeo Center database maintained by L. Petrov. The
National Radio Astronomy Observatory is a facility of the
National Science Foundation operated under cooperative
agreement by Associated Universities, Inc. This publication
makes use of data obtained at the Metsähovi Radio
Observatory, operated by Aalto University in Finland. This
research has made use of data from the MOJAVE database that
is maintained by the MOJAVE team (Lister et al. 2018). This
work makes use of the Murchison Radio-astronomy Observa-
tory, operated by CSIRO. We acknowledge the Wajarri
Yamatji people as the traditional owners of the observatory
site. Support for the operation of the MWA is provided by the
Australian Government (NCRIS) under a contract to Curtin
University, administered by Astronomy Australia Limited. We
acknowledge Paul Hancock, Gemma Anderson, John Morgan,
and Stefan Duchesne for their contributions to the GLEAM-X
pipeline that bring great convenience to us. This research has
made use of data from the OVRO 40 m monitoring program,
which was supported in part by NASA grants NNX08AW31G,
NNX11A043G, and NNX14AQ89G; NSF grants AST-
0808050 and AST-1109911; and private funding from Caltech
and the MPIfR. This research has made use of the NASA/
IPAC Extragalactic Database (NED 2019), which is funded by
the National Aeronautics and Space Administration and
operated by the California Institute of Technology.
Facilities: ASKAP, GMRT, HST, MWA, VLA, VLBA.
Software: AIPS (Greisen 2003), astropy (Astropy Collabora-
tion et al. 2013, 2018), Difmap (Shepherd et al. 1994).
Data Availability
The data sets underlying this paper were derived from the
public domain in the NRAO archive (project codes BU013 and
BA080; https://science.nrao.edu/observing/data-archive) and
Astrogeo archive (project codes BB023, RDV13, BG219D,
UF001B, and UG002U; http://astrogeo.org/). The MOJAVE
data can be found from the MOJAVE website (https://www.
cv.nrao.edu/MOJAVE/sourcepages/0007+106.shtml), MWA
archive (https://asvo.mwatelescope.org), and CSIRO ASKAP
Science Data Archive (CASDA; https://research.csiro.au/
casda/). The MWA GLEAM-X data are not publicly released,
and the calibrated visibility data can be shared on reasonable
request to the corresponding authors. Data from the Owens
Valley Radio Observatory and the Metsähovi Radio Observa-
tory can be requested from the respective data maintainers.
Appendix A
Low-frequency Observations of III Zw 2
We study the low-frequency radio structure of III Zw 2 based
on the recent observations from the MWA (Tingay et al. 2013),
9
The Astrophysical Journal, 944:187 (18pp), 2023 February 20 Wang et al.
GMRT (Swarup et al. 1991), and ASKAP (Johnston et al.
2007; Hotan et al. 2014).
The MWA is the precursor of the Square Kilometre Array
(SKA) low-frequency telescope and located in Western
Australia. We have downloaded the latest unpublished data
of III Zw 2 observed by MWA phase II (Wayth et al. 2018) on
2018 June 15. The extended array of MWA phase II comprises
128 tiles (each consisting of 16 crossed-pair dipole antenna)
distributed over a maximum baseline of ∼6 km, twice as long
as 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 S of decl. +30° 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σ≈ 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
archive14 and processed the data at the China SKA Regional
Centre (An et al. 2019, 2022) using the pipeline developed by
the GLEAM-X team.15 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 dewarping, and
flux density scaling to GLEAM in the CLEAN images. Then the
point-spread function 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 ratio image with reduced noise, allowing for detecting
fainter sources and diffuse structures that are not visible in
individual snapshot images.
Figure A1(a) shows the 216MHz MWA image of III Zw 2,
displaying a compact component. The resolution of the image
is 63 1× 49 1, and the rms noise is 2.1 mJy beam−1 (a factor
of ∼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 δ=+ 11°. Object III Zw 2 is
marginally resolved and shows an elongation along the N–S
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 α= −0.96± 0.09.
Figure A1(b) shows the GMRT image at 150 MHz, showing
a resolved structure along the NE–SW direction with a total
extent of ∼48″. The GMRT image of III Zw 2 is from the TIFR
GMRT Sky Survey Alternative Data Release (TGSS-ADR1;16
Intema et al. 2017). The TGSS-ADR1 covers 90% of the total
sky from δ= −53° to +90° with a noise level of ∼5 mJy
beam−1 and a resolution of 25″× 25″/cos(decl. −19°) at 150
MHz. The total flux density is ∼212 mJy, dominated by the
central core. Besides the core C, there are extensions toward the
NE and 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 S of δ= +41°) pilot survey for
future multiyear surveys with the full-scale ASKAP (McCon-
nell et al. 2020). The observations covering the III Zw 2
position began on 2019 April 21 and lasted 3 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 20× 12 93. The
distance between the NE and SW lobes is ∼33″, and the image
resolution enables one 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 Figure A1(c), which reveals much richer details than
the MWA image; the main jet body is elongated along the 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 Figure A1(c)) is fainter due to their steep-spectrum
nature. The separation between the N and S lobes is ∼55″ (∼98
kpc),17 larger than the size of the previously detected structure
between the NE and SW lobes in VLA images at gigahertz
Figure A1. Radio images of III Zw 2. (a) 215.5 MHz image observed from GLEAM-X on 2018 June 15. The peak contour flux is 126 mJy beam−1, and the contour
levels are 6.42 × (−1, 1, 1.4, 2, 2.8, 4, 5.6, 8, 11.2, 16, 23, 32) mJy beam−1. The rms noise is 2.14 mJy beam−1. The beam size is 63 1 × 49 1 with the major axis
along a position angle of 146°. 8. (b) 150 MHz image acquired by TGSS-ADR1 (Intema et al. 2017) observed on 2016 March 15. The integration time of each pointing
is ∼15 minutes with a beam size of 25″ × 25″. (c) 888 MHz image obtained from RACS (McConnell et al. 2020). The observation was made on 2020 May 3. The
integration time on the source is ∼15 minutes. The peak intensity is 51 mJy beam−1. The lowest contour level is 1 mJy beam−1, and the contours increase in steps of
2 . The rms noise in the image is 0.33 mJy beam−1. The overall size is 64 1 (∼106 kpc). The color scale shows the intensity in the logarithmic scale.
14 https://asvo.mwatelescope.org
15 https://github.com/tjgalvin/GLEAM-X-pipeline maintained by Tim
Galvin.
16 https://vo.astron.nl
17 Adopting the cosmological parameters H0 = 71 km s
−1 Mpc−1, ΩΛ = 0.73,
and Ωm = 0.27 at z = 0.0893 and 1 mas angular separation corresponds to a
1.65 pc linear size in projection on the plane of the sky.
10
The Astrophysical Journal, 944:187 (18pp), 2023 February 20 Wang et al.
frequencies (Brunthaler et al. 2005). In the 685MHz uGMRT
map, the southern extension (labeled ML in Silpa et al. 2021) is
much longer than our Figure A1(c). The total flux density
estimated from the ASKAP images is 92.5 mJy, of which the
main body (C+NE+SW) accounts for ∼86.3 mJy, while the
remaining 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 parts 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)
F
R
D
c p K B
c
4
3 2
, A1
L
p
p
,thin
3
2 5 0
1 2
1
1 2
( ) ( )( )
( )
⎜ ⎟⎛
⎝
⎞
⎠
p n=n +
- -
where R is the size of the extended emission structure, DL is the
luminosity distance, K0 is a normalization factor, B is the
magnetic field strength, ν 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=K0E
− pdE,
and energy equipartition between the magnetic fields and
particle kinetic energy densities, the normalization K0 can be
evaluated similarly to Equation (E2). The synchrotron
frequency corresponds to the minimum injected Lorentz factor
of emitting electrons and can be calculated as
eB
m c2
. A2m
m
e
2
( )n gp=
Using Equation (A2) and the normalization in Equation (A1),
γm= ((p− 2)/(p− 1))(mp/me)òe, in which ò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
B
D F
R c p p E
e
m cc
F R
D
6
2 4
1.6 10 G
mJy 100 kpc
100 Mpc
, A3
L
p
m
e
p
L
2
,thin
3
5 min
2
2
1
1 2 1 3
8 ,thin
1 3 1 3
2 3
( ) ( )
( )
( )
( )
⎜ ⎟
⎜ ⎟ ⎜ ⎟
⎜ ⎟
⎛
⎝⎜
⎛
⎝
⎞
⎠
⎞
⎠⎟
⎛
⎝
⎞
⎠
⎛
⎝
⎞
⎠
⎛
⎝
⎞
⎠
g
p= -
= ´
´
n
n
-
-
-
-

assuming òe= 0.1 and an index p= 1− 2αlow= 2.94 based on
the inferred spectral index α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 μ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
Figure A2. Combined radio (ASKAP) and optical (HST) images. The ASKAP image is obtained from RACS centered at 888 MHz. The HST image is acquired from
the Hubble Legacy Archive at a wavelength of 5446.8 Å.
11
The Astrophysical Journal, 944:187 (18pp), 2023 February 20 Wang et al.
extended region is
E U R
F
R D
4
3
1.82 10 erg
mJy
100 kpc 100 Mpc
, A4
e
L
ext
3 56 ,thin
2 3
7 3 4 3
( )
( )
⎜ ⎟
⎜ ⎟ ⎜ ⎟
⎛
⎝
⎞
⎠
⎛
⎝
⎞
⎠
⎛
⎝
⎞
⎠
p» = ´
´
n
where Ue= òUB= òB
2/(8π), with ò= òe/òB and using B from
Equation (A3). After Equation (A4) relating to the total energy
Eext, from the mean flux density measured in the MWA bands
(88–215 MHz), S ext 0.246( )á ñ =n mJy, and using a spectral
index αlow= −1.15 for the relics and a central frequency
ν= 151.6MHz, the radio luminosity of the extended emission
L D S z4 ext 1 6.18 10LR,ext
2 1 36low( ) ( )p n» á ñ + = ´n a- - erg s−1.
Using the empirical relations in Equations (E8) and (E9), the
associated luminosity of the extended jet is Lext= 5.19×
1039 erg s−1.
Appendix B
High-frequency VLBI Imaging of the Parsec-scale Jet
The source structure of III Zw 2 on parsec scales is revealed
from the VLBI imaging data, which include the archive data
from the MOJAVE (Lister et al. 2018) program and obtained
from the Astrogeo database.18 Details of the VLBI data are
presented in Table B1. 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 2GHz
8GHza for III Zw 2. Since the
source is unresolved in the 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 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, and (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 epochs were removed from the jet
kinematics analysis. The model fitting parameters are given in
Table B2. 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 observa-
tions, the flux density error also includes a certain level of
uncertainty in the flux scale calibration through error propaga-
tion. 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 the squares of the Q and U
components, i.e., p Q U2 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 a short integra-
tion time and sparse (u, v) coverage. Deconvolution and model
fitting are affected by strong side lobes 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 B1. 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
the jet kinematics analysis.
On parsec scales, VLBI images of III Zw 2 reveal either a
single core or a “core + compact jet” structure. The compact
core dominates the total flux density in all images. The jet
extends to the W and is relatively weaker but well
distinguished from the core in seven 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± 0.02 and 0.22± 0.01 mas yr−1,
respectively. These convert to jet transverse speeds of
1.35± 0.13c (J1) and 1.21± 0.07c (J2). The derived jet proper
motions are in good agreement with the previous studies
(Brunthaler et al. 2000) based on the 43 GHz VLBI observa-
tions in 1998. To distinguish it from our jet components, we
name the jet detected by Brunthaler et al. the “1998 jet.”
Extending the 1998 jet to epoch 2006 (the first epoch of the
MOJAVE data used in this paper), we find that 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.
18 http://astrogeo.org/
12
The Astrophysical Journal, 944:187 (18pp), 2023 February 20 Wang et al.
Appendix C
Flux Density Variability and Time Series Analysis
The single-dish light-curve data used in this paper are from
the archives of the 14 m MRO in Finland and the 40 m
telescope of OVRO in the United States.
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 onward to date. A
detailed description of the data reduction and analysis is given
in Terasranta et al. (1998). Under good observing conditions,
the detection limit of the MRO at 37 GHz is around 0.2 Jy.
Uncertainties in flux densities include the contribution from the
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 (δ>
−20°; 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 ∼3%. The high
Table B2
Model Fitting Results of 15 GHz VLBI Data
Time Comp. Sint R PA θFWHM TB δ
(yyyy-mm-dd) (mJy) (mas) (deg) (mas) (1010 K)
(1) (2) (3) (4) (5) (6) (7) (8)
2004-02-11 C 1390 ± 52 0 0 0.098 ± 0.001 79.6 ± 3.00 15.9
2005-03-05 C 1855 ± 186 0 0 0.225 ± 0.003 18.6 ± 2.05 3.7
2005-05-26 C 1890 ± 134 0 0 0.203 ± 0.001 23.9 ± 1.82 4.8
2005-12-22 C 1358 ± 156 0 0 0.380 ± 0.007 4.25 ± 0.61 0.9
2006-06-15 C 642 ± 120 0 0 0.175 ± 0.005 13.5 ± 2.2 2.7
J1 291 ± 9 0.417 ± 0.004 −67.2 ± 0.5 0.263 ± 0.008 0.42 ± 0.07 L
2007-08-09 C 296 ± 58 0 0 0.120 ± 0.004 11.9 ± 2.3 2.4
J1 81 ± 8 0.701 ± 0.006 −70.5 ± 0.5 0.437 ± 0.012 0.15 ± 0.02 L
2008-08-25 C 488 ± 24 0 0 0.072 ± 0.002 52.6 ± 2.6 10.5
J1 20 ± 3 1.035 ± 0.006 −67.8 ± 0.3 0.416 ± 0.012 0.04 ± 0.01 L
2009-06-03 C 1219 ± 44 0 0 0.084 ± 0.001 95.9 ± 3.4 19.2
J1 6.9 ± 1.5 1.123 ± 0.013 −68.1 ± 0.7 0.350 ± 0.026 0.03 ± 0.01 L
2010-07-12 C 1203 ± 74 0 0 0.173 ± 0.001 21.3 ± 1.4 4.3
2011-05-26 C 121 ± 22 0 0 0.147 ± 0.003 3.58 ± 1.38 0.7
J2 46 ± 2 0.390 ± 0.004 −63.2 ± 0.6 0.177 ± 0.008 0.16 ± 0.03 L
2012-04-30 C 141 ± 15 0 0 0.099 ± 0.007 8.11 ± 0.86 1.6
J2 11 ± 1 0.583 ± 0.004 −63.2 ± 0.4 0.328 ± 0.008 0.33 ± 0.01 L
2012-11-02 C 625 ± 18 0 0 0.048 ± 0.001 151.9 ± 4.4 30.4
J2 3.9 ± 0.8 0.714 ± 0.029 −62.6 ± 2.3 0.606 ± 0.058 0.01 ± 0.01 L
2013-06-02 C 565 ± 31 0 0 0.086 ± 0.001 42.3 ± 2.4 8.5
Note. Columns (3)–(7) present the model fitting parameters in sequence: the integrated flux density, radial separation, position angle with respect to the core C
(measured from N through E), component size (FWHM of the fitted Gaussian component), and brightness temperature. Note that the uncertainty of the position (R)
and component size (θFWHM) only takes into account the statistical error, which is inversely proportional to the signal-to-noise ratio of the component.
Table B1
Logs of the 15 GHz VLBI Observations
Date Project Code Bandwidth On-source Time Beam Speak σrms
(yyyy-mm-dd) (MHz) (hr) (Maj., Min., PA) (mJy beam−1) (mJy beam−1)
(1) (2) (3) (4) (5) (6) (7)
2004-02-11 BL111 32.0 1.05 1.2 × 0.6 mas2, −4°. 1 1363.9 0.19
2005-03-05 BL123 32.0 0.93 1.1 × 0.5 mas2, −4°. 3 1680.1 0.29
2005-05-26 BL123 32.0 0.93 1.1 × 0.6 mas2, −6°. 4 1755.0 0.31
2005-12-22 BL123 32.0 0.82 1.2 × 0.6 mas2, −5°. 1 1083.9 0.24
2006-06-15 BL137 32.0 0.93 1.2 × 0.6 mas2, −7°. 3 725.6 0.37
2007-08-09 BL149 32.0 0.93 1.2 × 0.6 mas2, −4°. 4 293.9 0.22
2008-08-25 BL149 32.0 1.05 1.2 × 0.5 mas2, −8°. 1 482.9 0.24
2009-06-03 BL149 64.0 1.05 1.2 × 0.6 mas2, −5°. 2 1204.1 0.16
2010-07-12 BL149CL 64.0 1.05 1.4 × 0.5 mas2, −7°. 7 1138.9 0.18
2011-05-26 BL149DI 64.0 1.17 1.2 × 0.6 mas2, −3°. 5 136.4 0.18
2012-04-30 BL178AJ 64.0 1.28 1.3 × 0.6 mas2, −6°. 6 140.3 0.17
2012-11-02 BL178AR 64.0 0.93 1.2 × 0.6 mas2, −4°. 9 622.7 0.15
2013-06-02 BL178BD 64.0 1.05 1.2 × 0.6 mas2, 0°. 1 557.4 0.19
Note. Column (5) gives the restoring beam in natural weighting; columns (6) and (7) present the peak flux density and rms noise in the images, respectively.
13
The Astrophysical Journal, 944:187 (18pp), 2023 February 20 Wang et al.
cadence and 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 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 kiloparsec-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 preprocessed by a
conversion from a partial uneven sampling with small data gaps
to a full even sampling after a linear interpolation and
resampling. The normalized Fourier periodogram P( fj) is
evaluated as (Vaughan 2005; Mohan & Mangalam 2014;
Mohan et al. 2015)
P f
t
x N
F f
2
, C1j j2
2( ) ∣ ( )∣ ( )= D
where Δt is the sampling time step size, x is the mean of the
light curve x(nΔt) of length N, and F( fj) is the Fourier
transform of x evaluated at the frequencies fj= j/(NΔ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( fj) will be characterized by a power-
law behavior, especially at low temporal frequencies. This can
be captured by parametric model fits to P( fj) to estimate the
underlying power spectral density (PSD). Here we use two
models, a power law given by
I f Af C, C2j j( ) ( )= +a
where A is the normalized amplitude, α is the power-law index,
and C is the ambient noise level, and a bending power law is
given by
I f Af f f C1 , C3j j j b
1l l( ) ( ( ) ) ( )= + +a a a- - -
where fb is a break frequency marking a transition in slopes, the
power-law indices are α for fj> fb and αl for fj< fb, 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 the parameters
and their errors are discussed in Mohan & Mangalam (2014)
and 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 χ2 statistical distribution), which are potential
candidates for quasiperiodic 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 resampled at the
original frequencies. A mean periodogram I fj˜( ) is constructed
from the simulations and rescaled 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 I fj j( ) ˜( ), are then 222c
distributed (Chatfield 2016), with a conditional probability
p P f I f I f ej j j
P f I f1 j j[ ( )∣ ˜( )] ˜ ( ) ( ) ˜( )= - - . The cumulative distribu-
tion function for the PSD ordinates is then the integral of the 2
2c
distribution, i.e., a gamma density function 1, 1 2( )G =
xexp 2 2( )- . Specifying a level of statistical significance
(1− ò) in the integral helps to identify outliers (quasiperiodic
signals) that may be present in the tail of the distribution. We
set a threshold of ò= 0.05 (95% level of significance) to
identify quasiperiodic signals in the light curves.
The wavelet analysis is employed here in a complementary
manner to probe quasiperiodic 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)
P n s W n s
W n s F f sf e
, ,
, 2 2 , C4
j
N
j j
if n t
2
1
2 j
( ) ∣ ( )∣
( ) ( ) ( ) ( )å p p
=
= Y p
=
D*
where W(n, s) is the wavelet transform of the evenly sampled
light curve x(nΔt) at times (nΔt) and scales s, F(2πfj) is the
Fourier transform of x evaluated at the circular frequencies
2πfj= 2πj/(NΔt) with j= 1, ..., N, and Ψ
*
(2πsfj) 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) e ei t t1 4 20 2y p= w- - , where ω0= 6 is a frequency
parameter characterizing the wavelet shape, and t is the
time parameter. In the frequency domain, sf2 j( )pY =
e sf1 4 2 2j 0
2( )p p w- - - . The wavelet scales s≈ 1.03/fj 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 implemen-
ted: 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 smooth 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.
14
The Astrophysical Journal, 944:187 (18pp), 2023 February 20 Wang et al.
The statistical significance of contiguous signals detected
in the wavelet analysis employs a twofold strategy. In the first
stage, the algorithm of Timmer & Koenig (1995) and the
expected 2
2c statistics of the 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 Af Cm m( ) = +a with associated para-
meters is used as a trial value. The index α is varied in the
range from the best-fit value to 0.0 in steps of −0.2 in each
case. A thousand realizations of the periodogram are
simulated, and they are oversampled in duration and temporal
frequencies and then resampled at the original frequencies. In
each realization, the periodogram is inverse Fourier trans-
formed to obtain a synthetic light curve with similar
statistical and variability properties as the original light
curve; the total number of simulated light curves is typically
…20,000. 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).
Appendix D
Integrated Radio Spectrum and Fitting
We plot the radio spectrum of III Zw 2 from 72 MHz to
37 GHz in Figure 4. The red filled squares below 230 MHz
are taken from the MWA GLEAM-X survey (Hurley-Walker
et al. 2022), and the green filled circle is from the GMRT
observation (Intema et al. 2017), revealing the extended radio
emission characterized by a steep spectrum. The slate blue
data point is from the ASKAP observation at 888 MHz. The
magenta filled circle 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
gigahertz frequencies observed in the epochs close to the
MWA and uGMRT observation: 37 GHz Metsähovi (blue
triangle), 15 GHz OVRO (navy triangle), and 8.4 and
2.3 GHz Astrogeo VLBI (purple diamond). From Figure 3,
we find that the core dominates the total flux density at
gigahertz 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 SSA,
and the extended jets and lobes become increasingly dominant.
The two-component radio spectrum spanning 0.072–37 GHz is
subjected to a weighted (the measurement errors of flux
densities are used as weights) least-squares fit using the
function Equation (D1), in which the first two parts
(A Blown +a ) form a power-law spectrum, and the last one
describes a spectrum of self-absorbed synchrotron radiation
(Pacholczyk 1970; Türler et al. 1999),
F A B F
e
e
1
1
3
2
1
8
1 , D1
m
m
m
thick
1 2
m m
m
low
thick
thick
( )
( )
⎜ ⎟
⎜ ⎟
⎜ ⎟
⎛
⎝
⎞
⎠
⎧
⎨⎩
⎫
⎬⎭
⎛
⎝
⎛
⎝
⎞
⎠
⎞
⎠
n nn
t aa
= + +
´ - -
= - -
n a
a
t n n
t
-
-
a a-
where the parameters of the low-frequency power-law
spectrum include the amplitude A (mJy), the spectral index
αlow, and a baseline flux density B (mJy), and the parameters of
the high-frequency section include the amplitude Fm (mJy), the
frequency of transition from optically thick to thin emission νm
(GHz), the optically thick spectral index αthick, the optically
thin spectral index α, and the optical depth τm expressed in
terms of αthick and α. In the fitting process derived using the
Markov Chain Monte Carlo method, αthick is fixed at 2.5,
corresponding to the canonical SSA case, and other parameters
are constrained within reasonable ranges: 0.0 Jy< A„ 0.8 Jy,
−1.5„ αlow„ −0.4, 0.0< B„ 0.6, 0.0 Jy< Fm„ 0.5 Jy, 7.0
GHz„ νm„ 15.0 GHz, and −0.4„ α„ 0.0. This yields A =
27 9
12-+ mJy, 0.97low 0.210.21a = - -+ , B 41 128= -+ mJy, F 146m 911= -+
mJy, 11.24m 1.25
2.01n = -+ GHz, and 0.20 0.120.10a = - -+ .
In Figure 4, we also plot the flux densities of the NE and SW
lobes obtained from the VLA images (Brunthaler et al. 2005)
and fit the NE and SW lobes (represented by blue open circles
and red open squares, respectively) with power-law functions.
We then calculated the flux densities of the 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 the 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 α= −1.09± 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.
Appendix E
Emission Properties of the Parsec–Kiloparsec Scale Jet
The synchrotron emission spectrum is assumed to be shaped
by a power-law energy distribution of the emitting relativistic
electrons N(E)dE= KE− pdE for energy E, normalization K,
and energy index p… 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 r r0 0 2( )= - and
B B r r0 0 1( )= - with scaling constants K0 and B0 at a fiducial
distance r0 (Mohan et al. 2015; Zdziarski et al. 2015; Agarwal
et al. 2017). With this, the particle kinetic energy density is
U N E E dE K
r
r
E
p 2
, E1e e
E
e
p
0
0
2
min
2
min
( )
( )
( )⎜ ⎟⎛
⎝
⎞
⎠ò= = -
¥ - - +
 
where òe is the fraction of the total energy density in the particle
kinetic energy, and Emin is the minimum energy required to
accelerate electrons (assuming that they are the dominant
15
The Astrophysical Journal, 944:187 (18pp), 2023 February 20 Wang et al.
constituents of the jet) to relativistic energies, here taken to be
the rest mass energy E 0.51 MeVmin = . Assuming equiparti-
tion of the total energy density between that in the magnetic
fields UB= B
2/(8π) and in the particle kinetic energy Ue, i.e.,
Ue/òe=UB/òB, the normalization constant is
K
B
p E
8
2 , E2p0
0
2
min
2( ) ( )p= -
-
where òB is the fraction of the total energy density in the
magnetic fields, and ò= òe/òB. The total energy density is then
U U U U
B r
r
1
8
1 . E3e B Btot
0
2
0
2
( ) ( ) ( )⎜ ⎟⎛
⎝
⎞
⎠p
= + = + = +
-
 
The jet luminosity attributable to synchrotron emission in a
region of size r sin jv q= (where θj is the jet half-ope
ning angle) downstream of the jet base is (Ghisellini &
Tavecchio 2010)
L cU
c
B r
8
sin 1 , E4j j
j j
jjet
2 2
tot
2
2
0 0
2( ) ( ) ( )pv b b q= G = G + 
where βj and Γj are the bulk velocity and Lorentz factor of the
jet. The optically thick absorption coefficient is (Rybicki &
Lightman 1986; Ghisellini 2013)
e K
m c
eB
m c
f p
8 2
, E5
e e
p
p
1 2 2 2 2
4 2( ) ( ) ( )
( )
( )⎜ ⎟⎛
⎝
⎞
⎠
a p p n= ¢n¢
+
- +
where e is the electron charge, n¢ is the emission frequency in
the source rest frame, and f (p) is expressed in terms of p and
the gamma function Γ dependence as
f p
p p p
3
3 22
12
3 2
12
6
12
. E6
p
p
1 2
8
12( )( )
( )
( )
⎛
⎝
⎞
⎠
⎛
⎝
⎞
⎠
⎛
⎝
⎞
⎠
=
G
G
+ G + G +
+
+
The optical depth for the emitting region ϖ is t a v=n n¢ ¢ .
From the condition 1t =n ¢ 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 SSA frequency z1( )n n d¢ = + , where ν is
the frequency in the observer frame, accounting for cosmolo-
gical redshift through the factor (1+ z) for a source at redshift z
and relativistic beaming through the Doppler factor δ. Using
Equations (E2), (E4), and (E5) in the above condition,
r
z
e E p f p
m c
e
m c
L
c
1
sin 2
64
2
8
1 sin
, E7
p
j
e
e j j j
1 1 2 2
min
2
jet
2 2
p
p
p
p
p
2
4
2
4
6
2 4
( )
( ) ( )
( )
( )
( ) ( )
( )
( )
( )
( )
⎜ ⎟
⎜ ⎟
⎛
⎝
⎞
⎠
⎛
⎝
⎞
⎠
⎛
⎝⎜
⎞
⎠⎟
n d p q
p
p q b
= +
-
´ + G
- - +
+
+
+
+


with the familiar r∝ ν−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 S C 907.15 90.87( )á ñ = n mJy
(see Table B2). The associated radio luminosity LR =
D S C z4 1 2.37 10L
2 1 42( ) ( )pn á ñ + = ´n a- + erg s−1 for
ν= 15 GHz, DL = 403.1 Mpc, and z = 0.089. We use
empirical relations to estimate the total jet luminosity (radiative
Ljet,rad and kinetic Ljet,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):
L Llog 12.00 0.75 log , E8Rjet,rad ( )= +
L Llog 6.00 0.90 log . E9Rjet,kin ( )= +
The total jet luminosity Ljet= Ljet,rad+ Ljet,kin= 1.97× 10
44
erg s−1.
Jet properties are estimated from the VLBI 15 GHz data
points during the flaring phases (δ… 4.3; median Doppler
factor). These include an average 16.9d˜ = and an associated
Γ= 8.5 based on 1 2app
2 2( ˜ ) ( ˜ )b d dG = + + . The jet half-
opening angle θj= 1/Γ= 6°.7, similar to the lower limit
estimated in Chamani et al. (2021). With the estimated Ljet,
p = 2.3, the assumption of equipartition in the energy density
between that in the magnetic fields and particle kinetic energy
with ò= 1 and the above physical properties, the radial distance
of the self-absorbed core from Equation (E7), rSSA≈ 435rG
(where rG=GM•/c
2= 1.47× 1013 cm is the gravitational
radius for a supermassive black hole, SMBH, of 1.84×
108Me; Grier et al. 2012). Using this in Equation (E4)
(r0= rSSA), the associated magnetic field strength is B0= 13.8
G. For the B∝ r−1 scaling relation, this corresponds to B
(r= 1 pc)= 52.8 mG. This estimate is consistent with the core
shift–based measurement of „60 mG and of a similar order of
magnitude to the SSA-based measurement of „20 mG
(Chamani et al. 2021).
Appendix F
Decomposition of the Major Flaring Phase during
2009–2010 and Origin
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. How-
ever, the superposition of these subflares 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 γ-ray flares were found during the
2009–2010 flare; one is associated with the peak of the 2009
flare, and the other is in the middle 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 F1), respectively. We need to
mention that the model curves shown in Figure F1 are not
obtained by least-squares fitting to the observed data but rather
16
The Astrophysical Journal, 944:187 (18pp), 2023 February 20 Wang et al.
by selecting model curves with minimum residuals from those
obtained with different parameter combinations.
In the main text, we inferred a correlation between γ-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 that the jet knots J1 and J2 move
along ballistic trajectories. Assuming that the jet speed remains
constant, we can trace back to obtain the position of jet J2 at the
moments of the γ-ray and radio flares. They are in sequence:
0.041 (2009 γ-ray flare), 0.096 (2009 radio flare), 0.163 (2010
γ-ray flare), and 0.205 (2010 radio flare) mas. The mass of the
SMBH in III Zw 2 is (1.84± 0.27)× 108 Me (Grier et al.
2012); therefore, its gravitational radius is Rg= 9× 10
−6 pc.
The physical size corresponding to the first γ-ray flare in 2009
is 7.6× 103Rg, 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× 107 Me was 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 γ-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 BLR 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 BLR or the
clouds in the torus on a larger scale, resulting in the observed
flares (see discussion in Section 3.7).
ORCID iDs
Ailing Wang https://orcid.org/0000-0002-7351-5801
Tao An https://orcid.org/0000-0003-4341-0029
Shaoguang Guo https://orcid.org/0000-0003-0181-7656
Prashanth Mohan https://orcid.org/0000-0002-2211-0660
Wara Chamani https://orcid.org/0000-0001-7350-4152
Willem A. Baan https://orcid.org/0000-0003-3389-6838
Talvikki Hovatta https://orcid.org/0000-0002-2024-8199
Heino Falcke https://orcid.org/0000-0002-2526-6724
Tim J. Galvin https://orcid.org/0000-0002-2801-766X
Natasha Hurley-Walker https://orcid.org/0000-0002-
5119-4808
Sumit Jaiswal https://orcid.org/0000-0002-5125-695X
Anne Lahteenmaki https://orcid.org/0000-0002-0393-0647
Baoqiang Lao https://orcid.org/0000-0002-3426-3269
Merja Tornikoski https://orcid.org/0000-0003-1249-6026
Yingkang Zhang https://orcid.org/0000-0001-8256-8887
References
Agarwal, A., Mohan, P., Gupta, A. C., et al. 2017, MNRAS, 469, 813
Aller, H. D., Aller, M. F., Latimer, G. E., & Hodge, P. E. 1985, ApJS, 59, 513
Aller, M. F., Aller, H. D., & Hughes, P. A. 2003, ApJ, 586, 33
An, T., & Baan, W. A. 2012, ApJ, 760, 77
An, T., Baan, W. A., Wang, J.-Y., Wang, Y., & Hong, X.-Y. 2013, MNRAS,
434, 3487
An, T., Mohan, P., Zhang, Y., et al. 2020, NatCo, 11, 143
An, T., Wu, X., Lao, B., et al. 2022, SCPMA, 65, 129501
An, T., Wu, X.-P., & Hong, X. 2019, NatAs, 3, 1030
Arp, H. 1968, ApJ, 152, 1101
Asada, K., & Nakamura, M. 2012, ApJL, 745, L28
Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ,
156, 123
Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A,
558, A33
Bardeen, J. M., & Petterson, J. A. 1975, ApJL, 195, L65
Biretta, J. A., Zhou, F., & Owen, F. N. 1995, ApJ, 447, 582
Blandford, R., & Eichler, D. 1987, PhR, 154, 1
Blandford, R., Meier, D., & Readhead, A. 2019, ARA&A, 57, 467
Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883
Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433
Boccardi, B., Krichbaum, T. P., Bach, U., et al. 2016, A&A, 585, A33
Boccardi, B., Perucho, M., Casadio, C., et al. 2021, A&A, 647, A67
Bogdanović, T., Reynolds, C. S., & Miller, M. C. 2007, ApJL, 661, L147
Brunthaler, A., Falcke, H., Bower, G. C., et al. 2000, A&A, 357, L45
Brunthaler, A., Falcke, H., Bower, G. C., et al. 2003, PASA, 20, 126
Brunthaler, A., Falcke, H., Bower, G. C., et al. 2005, A&A, 435, 497
Callingham, J. R., Ekers, R. D., Gaensler, B. M., et al. 2017, ApJ, 836, 174
Camenzind, M., & Krockenberger, M. 1992, A&A, 255, 59
Chamani, W., Koljonen, K., & Savolainen, T. 2020, A&A, 635, A172
Chamani, W., Savolainen, T., Hada, K., & Xu, M. H. 2021, A&A, 642, A14
Chatfield, C. 2016, The Analysis of Time Series: An Introduction (6th ed.;
Boca Raton, FL: CRC Press)
Clements, S. D., Smith, A. G., Aller, H. D., & Aller, M. F. 1995, AJ, 110, 529
Dodge, Y. 2008, The Concise Encyclopedia of Statistics (Berlin: Springer)
Ebisuzaki, T., Makino, J., & Okumura, S. K. 1991, Natur, 354, 212
Figure F1. Decomposition of the major flare in 2009–2010 with two exponential components. Left: 37 GHz light curve. Right: 15 GHz light curve. The uncertainties
of the fitted parameters only take into account the statistical error.
17
The Astrophysical Journal, 944:187 (18pp), 2023 February 20 Wang et al.
Ekers, R. D., Fanti, R., Lari, C., & Parma, P. 1978, Natur, 276, 588
Falcke, H., & Biermann, P. L. 1995, A&A, 293, 665
Falcke, H., Bower, G. C., Lobanov, A. P., et al. 1999, ApJL, 514, L17
Falcke, H., Malkan, M. A., & Biermann, P. L. 1995, A&A, 298, 375
Falcke, H., Patnaik, A. R., & Sherwood, W. 1996a, ApJL, 473, L13
Falcke, H., Sherwood, W., & Patnaik, A. R. 1996b, ApJ, 471, 106
Fomalont, E. B. 1999, in ASP Conf. Ser. 180, Synthesis Imaging in Radio
Astronomy II, ed. G. B. Taylor, C. L. Carilli, & R. A. Perley (San
Francisco, CA: ASP), 301
Foschini, L. 2014, IJMPS, 28, 1460188
Frank, J., King, A., & Raine, D. J. 2002, Accretion Power in Astrophysics (3rd
ed.; Cambridge: Cambridge Univ. Press)
Fukumura, K., Tombesi, F., Kazanas, D., et al. 2014, ApJ, 780, 120
Ghisellini, G. 2013, Radiative Processes in High Energy Astrophysics, Vol.
873 (Cham: Springer)
Ghisellini, G., Padovani, P., Celotti, A., & Maraschi, L. 1993, ApJ, 407, 65
Ghisellini, G., & Tavecchio, F. 2010, MNRAS, 409, L79
Giroletti, M., Panessa, F., Longinotti, A. L., et al. 2017, A&A, 600, A87
Gopal-Krishna, Biermann, P. L., Wiita, P. J., et al. 2003, ApJL, 594, L103
Greisen, E. W. 2003, in AIPS, the VLA, and the VLBA, ed. A. Heck, Vol. 285
(Dordrecht: Springer), 109
Grier, C. J., Peterson, B. M., Pogge, R. W., et al. 2012, ApJ, 755, 60
Grossmann, A., Kronland-Martinet, R., & Morlet, J. 1989, in Wavelets. Time-
Frequency Methods and Phase Space, ed. J.-M. Combes, A. Grossmann, &
P. Tchamitchian (Berlin: Springer), 2
Hada, K., Kino, M., Doi, A., et al. 2013, ApJ, 775, 70
Hardee, P. E. 2003, ApJ, 597, 798
Harrison, C. M., Costa, T., Tadhunter, C. N., et al. 2018, NatAs, 2, 198
Hotan, A. W., Bunton, J. D., Chippendale, A. P., et al. 2021, PASA, 38, e009
Hotan, A. W., Bunton, J. D., Harvey-Smith, L., et al. 2014, PASA, 31, e041
Hovatta, T., Nieppola, E., Tornikoski, M., et al. 2008, A&A, 485, 51
Hovatta, T., Valtaoja, E., Tornikoski, M., & Lähteenmäki, A. 2009, A&A,
494, 527
Hurley-Walker, N., Callingham, J. R., Hancock, P. J., et al. 2017, MNRAS,
464, 1146
Hurley-Walker, N., Galvin, T. J., Duchesne, S. W., et al. 2022, PASA, 39, e035
Hutchings, J. B., & Campbell, B. 1983, Natur, 303, 584
Hutchings, J. B., Campbell, B., Gower, A. C., Crampton, D., & Morris, S. C.
1982, ApJ, 262, 48
Intema, H. T., Jagannathan, P., Mooley, K. P., & Frail, D. A. 2017, A&A,
598, A78
Ivezić, Ž., Menou, K., Knapp, G. R., et al. 2002, AJ, 124, 2364
Johnston, S., Bailes, M., Bartel, N., et al. 2007, PASA, 24, 174
Jones, F. C., & Ellison, D. C. 1991, SSRv, 58, 259
Jorstad, S. G., Marscher, A. P., Raiteri, C. M., et al. 2022, Natur, 609, 265
Kaastra, J. S., & de Korte, P. A. J. 1988, A&A, 198, 16
Kellermann, K. I., Sramek, R., Schmidt, M., Shaffer, D. B., & Green, R. 1989,
AJ, 98, 1195
Kellermann, K. I., Sramek, R. A., Schmidt, M., Green, R. F., & Shaffer, D. B.
1994, AJ, 108, 1163
Kellermann, K. I., Vermeulen, R. C., Zensus, J. A., & Cohen, M. H. 1998, AJ,
115, 1295
King, A., & Pounds, K. 2015, ARA&A, 53, 115
Kukula, M. J., Dunlop, J. S., Hughes, D. H., & Rawlings, S. 1998, MNRAS,
297, 366
Leahy, J. P., & Williams, A. G. 1984, MNRAS, 210, 929
Lebofsky, M. J., & Rieke, G. H. 1980, Natur, 284, 410
Lehman, A. 2005, JMP for Basic Univariate and Multivariate Statistics: A
Step-by-step Guide (Cary, NC: SAS Institute)
Li, H. Z., Xie, G. Z., Dai, H., et al. 2010, NewA, 15, 254
Liao, N.-H., Xin, Y.-L., Fan, X.-L., et al. 2016, ApJS, 226, 17
Lister, M. L., Aller, M. F., Aller, H. D., et al. 2018, ApJS, 234, 12
Lister, M. L., Homan, D. C., Hovatta, T., et al. 2019, ApJ, 874, 43
Lloyd, C. 1984, MNRAS, 209, 697
Lobanov, A. P. 1998, A&A, 330, 79
Marscher, A. P., Jorstad, S. G., D’Arcangelo, F. D., et al. 2008, Natur, 452, 966
McConnell, D., Hale, C. L., Lenc, E., et al. 2020, PASA, 37, e048
Miller, P., Rawlings, S., & Saunders, R. 1993, MNRAS, 263, 425
Mohan, P., Agarwal, A., Mangalam, A., et al. 2015, MNRAS, 452, 2004
Mohan, P., An, T., Frey, S., et al. 2016a, MNRAS, 463, 1812
Mohan, P., Gupta, A. C., Bachev, R., & Strigachev, A. 2016b, MNRAS,
456, 654
Mohan, P., & Mangalam, A. 2014, ApJ, 791, 74
Müller, A., Pfrommer, C., Ignesti, A., et al. 2021, MNRAS, 508, 5326
Myers, J. L., Well, A. D., & Lorch, R. F., Jr. 2013, Research Design and
Statistical Analysis (New York: Routledge)
NASA/IPAC Extragalactic Database (NED) 2019, NASA/IPAC Extragalactic
Database (NED), IPAC, doi:10.26132/NED1
Osterbrock, D. E. 1977, ApJ, 215, 733
Pacholczyk, A. G. 1970, Radio Astrophysics. Nonthermal Processes in
Galactic and Extragalactic Sources (San Francisco: Freeman)
Padovani, P. 2016, A&ARv, 24, 13
Panessa, F., Baldi, R. D., Laor, A., et al. 2019, NatAs, 3, 387
Piccinotti, G., Mushotzky, R. F., Boldt, E. A., et al. 1982, ApJ, 253, 485
Pringle, J. E. 1996, MNRAS, 281, 357
Quici, B., Hurley-Walker, N., Seymour, N., et al. 2021, PASA, 38, e008
Readhead, A. C. S. 1994, ApJ, 426, 51
Reynolds, C., Punsly, B., Miniutti, G., O’Dea, C. P., & Hurley-Walker, N.
2017, ApJ, 836, 155
Rezzolla, L., Yoshida, S., Maccarone, T. J., & Zanotti, O. 2003a, MNRAS,
344, L37
Rezzolla, L., Yoshida, S., & Zanotti, O. 2003b, MNRAS, 344, 978
Richards, J. L., Max-Moerbeck, W., Pavlidou, V., et al. 2011, ApJS, 194, 29
Rybicki, G. B., & Lightman, A. P. 1986, Radiative Processes in Astrophysics
(Hoboken, NJ: Wiley)
Salvi, N., Page, M. J., Stevens, J. A., Mason, K. O., & Wu, K. 2002a, PASA,
19, 73
Salvi, N. J., Page, M. J., Stevens, J. A., et al. 2002b, MNRAS, 335, 177
Sargent, W. L. W. 1970, ApJ, 160, 405
Scargle, J. D. 1989, ApJ, 343, 874
Schmidt, M., & Green, R. F. 1983, ApJ, 269, 352
Schnopper, H. W., Delvaille, J. P., Epstein, A., et al. 1978, ApJL, 222, L91
Sembay, S., Hanson, C. G., & Coe, M. J. 1987, MNRAS, 226, 137
Shepherd, M. C., Pearson, T. J., & Taylor, G. B. 1994, BAAS, 26, 987
Shepherd, M. C., Pearson, T. J., & Taylor, G. B. 1995, BAAS, 27, 903
Shi, F., Li, Z., Yuan, F., & Zhu, B. 2021, NatAs, 5, 928
Shulevski, A., Morganti, R., Barthel, P. D., et al. 2015, A&A, 583, A89
Silpa, S., Kharb, P., Harrison, C. M., et al. 2021, MNRAS, 507, 991
Silpa, S., Kharb, P., Ho, L. C., et al. 2020, MNRAS, 499, 5826
Surace, J. A., Sanders, D. B., & Evans, A. S. 2001, AJ, 122, 2791
Swarup, G., Ananthakrishnan, S., Kapahi, V. K., et al. 1991, CSci, 60, 95
Teräsranta, H., Achren, J., Hanski, M., et al. 2004, A&A, 427, 769
Terasranta, H., Tornikoski, M., Mujunen, A., et al. 1998, A&AS, 132, 305
Terasranta, H., Tornikoski, M., Valtaoja, E., et al. 1992, A&AS, 94, 121
Teräsranta, H., Wiren, S., Koivisto, P., Saarinen, V., & Hovatta, T. 2005,
A&A, 440, 409
Timmer, J., & Koenig, M. 1995, A&A, 300, 707
Tingay, S. J., Goeke, R., Bowman, J. D., et al. 2013, PASA, 30, e007
Tombesi, F., Sambruna, R. M., Marscher, A. P., et al. 2012, MNRAS,
424, 754
Torrence, C., & Compo, G. P. 1998, BAMS, 79, 61
Türler, M., Courvoisier, T. J. L., & Paltani, S. 1999, A&A, 349, 45
Unger, S. W., Lawrence, A., Wilson, A. S., Elvis, M., & Wright, A. E. 1987,
MNRAS, 228, 521
Urry, C. M., & Padovani, P. 1995, PASP, 107, 803
Valtaoja, E., Lähteenmäki, A., Teräsranta, H., & Lainela, M. 1999, ApJS,
120, 95
Valtonen, M. J., Lehto, H. J., Nilsson, K., et al. 2008, Natur, 452, 851
Vaughan, S. 2005, A&A, 431, 391
Wang, A., An, T., Cheng, X., et al. 2023, MNRAS, 518, 39
Wang, A., An, T., Jaiswal, S., et al. 2021, MNRAS, 504, 3823
Wang, J.-Y., An, T., Baan, W. A., & Lu, X.-L. 2014, MNRAS, 443, 58
Wayth, R. B., Tingay, S. J., Trott, C. M., et al. 2018, PASA, 35, e033
Worrall, D. M., Birkinshaw, M., & Cameron, R. A. 1995, ApJ, 449, 93
Wright, A. E., Allen, D. A., Krug, P. A., Morton, D. C., & Smith, M. G. 1977,
IAUC, 3145, 2
Yuan, F., & Narayan, R. 2014, ARA&A, 52, 529
Zdziarski, A. A., Sikora, M., Pjanka, P., & Tchekhovskoy, A. 2015, MNRAS,
451, 927
Zwicky, F. 1967, AdA&A, 5, 267
18
The Astrophysical Journal, 944:187 (18pp), 2023 February 20 Wang et al.