The X-shaped radio galaxy J0725+5835 is associated with an AGN pair

X-shaped radio galaxies (XRGs) are those that exhibit two pairs of unaligned radio lobes (main radio lobes and"wings"), one of the promising models for the peculiar morphology is jet re-orientation. To clarify it, we conducted the European VLBI Network (EVN) 5 GHz observation of an XRG J0725+5835, which resembles the archetypal binary AGNs 0402+379 in radio morphology but it is larger in angular size. In our observation, two milliarcsec (mas) scale radio components with non-thermal radio emission are detected, each of them coincides with an optical counterpart with similar photometric redshift and (optical and infrared) magnitude, corresponding to dual active nuclei. Furthermore, with the improved VLA images, we find a bridge between the two radio cores and a jet bending in the region surrounding the companion galaxy, which further supports the interplay between the main and companion galaxies. In addition, we also report the discovery of an arcsec-scale jet in the companion. Given the projected separation of $\sim100$ kpc between the main and companion galaxies, XRG J0725+5835 is likely associated with a dual jetted-AGN system. In both EVN and VLA observations, we find signatures that the jet is changing its direction, which is likely responsible for the X-shaped morphology. On the origin of jet re-orientation, several scenarios are discussed.


INTRODUCTION
Astronomical observations show that supermassive black holes (SMBHs) are located at the center of most massive galaxies (e.g., Kormendy & Ho 2013). According to the hierarchical model of structure formation, galaxy mergers occur frequently in the early Universe (Volonteri et al. 2003), which naturally results in two SMBHs that sink to the center of a newly formed galaxy due to dynamical friction and form a binary system. Moreover, galaxy mergers can efficiently drive gas inflow to fuel the SMBHs, thus stimulating the active galactic nuclei (AGNs) and nuclear starburst (Hopkins et al. 2006;Di Matteo et al. 2007; Montuori et al. 2010). To establish a complete view of galaxy merger evolution, one needs to further explore, e.g. how the strength of star formation, AGN accretion, and feedback evolve with the merger sequence. Even though dual and binary AGNs candidates with separations ranging from 100 kpc down to few parsec scales have been detected (Komossa et al. 2003;Bianchi et al. 2008;Fu et al. 2011;Koss et al. 2011;Comerford et al. 2011;Woo et al. 2014;Deane et al. 2014;Fu et al. 2015;Husemann et al. 2018;Yang et al. 2017;Liu et al. 2018;Goulding et al. 2019;Silverman et al. 2020), the number of confirmed dual and binary AGNs is still smaller than predicted , and one reason for this deficiency is the lack of effective candidates. Here we define dual AGNs as gravitational interacting systems with a separation smaller than ≈ 100 kpc.
Under suitable conditions, an AGN can launch a pair of relativistic jets of non-thermal radio emission with dimensions that can evolve up to megaparsecs (e.g., Heckman & Best 2014; Dabhade et al. 2017). Ekers et al. (1978) noted that a symmetric jet distortion in galaxy NGC 326 can be caused by jet precession, while further studies indicate that a more rapid change of jet direction is required to form a nearly orthogonal distortion in NGC 326 (e.g. Murgia et al. 2001) and other morphologically similar sources (i.e. X-shaped radio galaxies, XRGs, Dennett- Thorpe et al. 2002). There are two scenarios concerning the jet re-orientation: jet precession due to binary black hole evolution (Zhang et al. 2007;Bogdanović et al. 2007;Krause et al. 2019;Horton et al. 2020) and the sudden spin-flip after binary black hole coalescence (Merritt & Ekers 2002). On the other hand, backflowing plasma from the jet may also be warped as the "wing-shape" in XRGs (Rossi et al. 2017;Cotton et al. 2020).
This work focuses on the XRGs and their origin is still debated (see Yang et al. 2019, and references therein). The jet precession model requires that binary black holes are located in parsec-scale (Krause et al. 2019), while it lacks observational support. Interestingly, we find observational evidence that several XRGs are in kilo-parsecscale galaxy pairs (e.g. Battistini et al. 1980;Leahy & Williams 1984;Murgia et al. 2001;Madrid et al. 2006). A study of ∼100 radio sources associated with dumbbell galaxies, i.e., two nearly equally bright elliptical galaxies within a common envelope, has revealed remarkably distorted radio structures in about a dozen of them (Wirth et al. 1982). Currently, there is enough observational evidence to show that galaxy pairs are associated with jet distortion (e.g. Parma et al. 1991;Borne & Colina 1993;Hardcastle et al. 2019), therefore, it occasionally forms X-shaped structures. J0725+5835 was selected as an X-shaped radio galaxy (Cheung 2007) because the directions of its primary and secondary jets are almost orthogonal. What caught our attention is its similar morphology with the parsec-scale binary AGN 0402+379 (Rodriguez et al. 2006). There is a radio core at the symmetric center of J0725+5835 driving the large-scale radio jet, while a compact (at arcsec-scales) radio component to its north-east direction is likely a massive neighboring galaxy that could affect the jet shape of the main radio galaxy. In the following, we refer to the radio core of the primary jet as core B and the compact radio component in the northeast direction as core A (see Figure 1); both core A and core B have Panoramic Survey Telescope and Rapid Response System (Pan-STARRS) counterparts (see Pan-STARRS1 i-band image in Figure 1) and show similar photometric redshift z ∼ 0.42 (see Section 4.1). In this paper, we present the results of observations of the XRG J0725+5835 by the European Very Long Baseline Interferometry (VLBI) Network (EVN) and the Jansky Very Large Array (VLA). This paper is organized as follows: we describe the EVN and VLA observations and the data reduction in Section 2 and present the imaging results of J0725+5835 in Section 3. We discuss the interplay between cores A and B, jet re-orientation, and simultaneous scenarios in Section 4. Throughout the paper, a standard ΛCDM cosmology is used with H 0 = 71 km s −1 Mpc −1 , Ω m = 0.27, and Ω Λ = 0.73.

VLBI observation and data reduction
We observed J0725+5835 on 2017 November 14 with ten antennas of the EVN at C-band (5 GHz, the project code RSY06, PI: Xiaolong Yang). The total observation time is 2 h with a data recording rate of 2 Gbps. The observation was performed in phase-referencing mode, using 0724+571 (R.A. = 07 h 28 m 49 s .6317, DEC. = 57 • 01 24 .374) as the phase reference calibrator.
We calibrated the EVN data in the Astronomical Image Processing System (AIPS), a software package developed by the National Radio Astronomy Observatory (NRAO) of U.S. (Greisen 2003), following the standard procedure. Apriori amplitude calibration was performed using the system temperatures and antenna gain curves provided by each station. The Earth orientation parameters were obtained and calibrated using the measurements from the U.S. Naval Observatory database, and the ionospheric dispersive delays were corrected from a map of the total electron content provided by the Crustal Dynamics Data Information System (CD-DIS) of NASA 1 . The opacity and parallactic angles were also corrected using the auxiliary files attached to the data. The delay in the visibility phase was solved using the phase reference calibrator 0724+571. A global fringe-fitting on the phase-referencing calibrator 0724+571 was performed by taking the calibrator's model to solve miscellaneous phase delays of the target. The target source's data were exported to DIFMAP (Shepherd 1997) for self-calibration and model fitting. Especially, in our VLBI data reduction, we take the radio coordinate from the VLA B-array S-band (see Section 2.2) as a reference for core A, which proved to be sufficiently accurate (see Table 1). No self-calibration was applied to the target source since it is too weak (∼ 7σ and 13σ for cores A and B, respectively). The final image was created using natural weighting; see Figure 2.

VLA observation and data reduction
This source is a part of the VLA L-(1.4 GHz) and S-band (3 GHz) surveys of 89 XRGs . In this work, we have retrieved two VLA Barray datasets at L-and S-band (project code 16A-220, PI: Lakshmi Saripalli) and one VLA A-array dataset at L-band (project code 16B-023, PI: Lakshmi Roberts et al. (2018) and , in order to ensure uniformity in the analysis across all the parameters we have re-reduced the data using the Common Astronomy Software Application (CASA v5.1.1, McMullin et al. 2007).
Our data analysis followed the standard routines described in the CASA Cookbook 3 . We used 3C 138 and 3C 147 as the primary flux calibrators for VLA B-array and A-array datasets, respectively, and adopted the flux density standards 'Perley-Butler 2017' (Perley & Butler 2017) to set the overall flux density scale; after that, we bootstrapped to the secondary flux density calibrator (J0713+4349) and the target. Antenna delay and bandpass corrections were also determined by fringe-fitting the visibilities. Furthermore, J0713+4349 was also used as a phase calibrator in the observations. We determined the complex gains from the phase calibrator and applied it to the target. We also performed an ionospheric correction using the data obtained from the CDDIS archive. Deconvolution, self-calibration, and model fitting were performed in the DIFMAP software package. For the good uv-coverage of the VLA observations and the high signal-to-noise ratio (SNR>9), the VLA data allowed self-calibration. This was initially performed only on phase, and subsequently on both phase and amplitude when we achieved a good model. Natural weighting was used to create the final image.  Table 1). Where the arrows are defined as R0: the direction of fossil radio emission, R1: the wings, R2: the direction of lobes, R3: the direction from core to the tip of hotspots, and the green arrow and R4: the arrows are obtained from left and right panel, respectively, indicate the direction from VLBI to VLA position.  Note-Columns give (1) source name/facilities, (2) frequency, (3-6) J2000 positions and errors, (7) total flux density, (8) source angular size, (9) brightness temperature, and (10) peak brightness.
In Figure 3, we show the VLA B-array images obtained at L-and S-band and the VLA A-array image at L-band. In the B-array L-band observation, we obtained a synthesized beam size of 6.24 × 4.83 at a position angle of 89.6 • , while in the A-array L-band observation, we obtained a synthesized beam size of 1.83 × 1.21 at a position angle of 78.1 • . In the VLA B-array S-band observation, we obtained a synthesized beam size of 2.83 × 2.13 at a position angle of 61.5 • . We perform a two-dimensional Gaussian model-fitting for the VLA Aarray L-band and B-array S-band data, and the results are listed in Table 1, where the uncertainties of the inte-grated flux density and peak flux density are estimated using the method described by Yang et al. (2020). In our EVN observation, both core A and core B were detected at 4.92 GHz (see Figure 2). The peak flux densities of A and B are 0.084 ± 0.024 and 0.280 ± 0.038 mJy beam −1 , yielding an SNR of 6.8 and 13.3, respectively. For the EVN data, two-dimensional Gaussian models were used to fit the visibility data for each target to obtain the integrated flux density and size of the Gaussian components (full width at half maximum, FWHM). The model-fitting results are listed in Table 1. For the VLA data, we fit two-dimensional Gaussian models to cores A and B only to obtain their integrated flux densities of 3.865 ± 0.178 and 0.973 ± 0.174 mJy at 1.54 GHz (VLA A-array), and 2.822 ± 0.036 and 1.196 ± 0.026 mJy at 3 GHz (VLA B-array), respectively. By combining the radio flux densities at different frequencies, we can obtain the spectral indices (S ∝ ν +α ) of core A and core B as −0.47 ± 0.07 and +0.30 ± 0.27, respectively.
By fitting two-dimensional Gaussian models to the EVN data, we can estimate the brightness temperatures using the formula (e.g. Ulvestad et al. 2005) are 10 6.37 and 10 7.28 K, respectively, consistent with a non-thermal origin for the mas-scale radio emission.
In the phase-referenced EVN observation, the astrometric accuracy can be measured by considering the po-sitional uncertainty of phase calibrator 0724+571 (σ p = 0.14 mas), the errors in phase-transferring from calibrator to target (σ pr,α = 0.035 mas and σ pr,δ = 0.117 mas) (Pradel et al. 2006) and the rms error of target. The rms error of the target can be estimated as θ FWHM /2SNR, where θ FWHM is the FWHM of the beam and SNR is the signal-to-noise ratio. While for VLA observations, only the rms astrometric error is considered.

Interplay between A and B
The EVN detection of compact radio cores with high brightness temperatures (> 10 6 K) and flat radio spectra (α > −0.5, spectral index of core A marginally satisfies the criteria, but see Figure 5) supports radio-frequency active AGNs located at both galactic nuclei in A and B. The Dark Energy Camera Legacy Survey (DECaLS, Dey et al. 2019) measured a photometric redshift of 0.399±0.057 and 0.442±0.028 for cores A and B, respectively. The similar redshifts between A and B support the possibility that they are in a close pair system.
Since we have no spectroscopic redshift measurements on this source, we checked radio properties to further explore the features supporting an interplay between A and B. Clearly, the radio bridge between cores A and B, as well as the deflection of the eastern radio outflow (jet relics) in the core A region (see the blue arrow in the middle panel of Figure 2), both hint the interplay between cores A and B. Again, the radio structure around core A is unlikely from itself alone, because an arcsec-scale jet is found in both VLA A-array L-band and B-array S-band images (see Figure 4), and shows a different extension. Additional evidence of the radio bridge between A and B is from the radio spectral index map ( Figure 5), which shows the continuous distribution between A and B. Furthermore, the interaction between the radio outflow and core A may also be inferred from the enhanced polarization in the circum-A region (supplementary data of Roberts et al. 2018). In the case here, the explanation for the deflection of the radio outflow is straightforward, which could be a combined effect of galactic ram pressure (Fiedler & Henriksen 1984), density gradients (Smith 1984) or oblique magnetic field (Koide 1997;Chibueze et al. 2021) from core A. Since core A is radio-frequency active, it is reasonable to suppose that multiple mechanisms are at play. According to our EVN 5 GHz observation, the angular distance between cores A and B is 19.5 arcsec. If we take the averaged photometric redshift between cores A and B (z = 0.42, the scale is 5.504 kpc arcsec −1 ) as a reference, the projected distance between the two AGNs is about 107 kpc, a marginally dual AGN system (Liu et al. 2011). This also results in a projected distance of ∼ 450 kpc between the northern and southern hotspots. Cores A and B have nearly equal apparent r-band magnitudes of 20.54 and 20.70 in the Pan-STARRS sky survey (Chambers et al. 2016), and there appears to be a tidal-tail structure in the image of galaxy A (see Figure  1), which further hints the interplay between two galax- Figure 5. Radio spectral index map between 1.5 and 3 GHz. The spectral index map was created from the VLA B-array S-band and B-array L-band data. The corresponding spectral index scale is linked with a color map in the color bar, where the yellow regions are with the flat spectrum, while the green regions are with the steep spectrum. Here only the regions with radio flux density above 7σ in both frequencies are presented. The cyan contours are from VLA B-array L-band data and same with the left panel of Figure 3, and its beam shows at the lower-right corner with a cyan ellipse.
ies. The association of an X-shaped jet morphology in interacting galaxies is quite common (e.g. Wirth et al. 1982;Murgia et al. 2001;Machacek et al. 2007).

Jet is changing its direction?
With the EVN and VLA observations, we can measure the coordinates of cores A and B (see Table 1). The coordinates measured by the VLA are consistent between A-array L-band and B-array S-band, while the B-array S-band has higher accuracy, so we will use it in the discussion below. At the position of core B, our EVN 5 GHz measurement has a deviation from the VLA. Interestingly, this deviation is not caused by an error but is intrinsic, i.e. if we take the EVN 5 GHz location of core B as the reference, then the location of core B measured by the VLA B-array at S-band is at R.A. offset= −67 ± 9 mas and DEC offset= −74 ± 9 mas (red cross in the right panel of Figure 2), beyond the 7σ astrometric error of VLA. Furthermore, it's unlikely to be a systematic error, compared to the VLA coordinates of core A which are close to the EVN measured position (see left panel of Figure 2).
A reasonable explanation for the position offset is that there is an inner < 10 kpc-scale jet in the direction −124 • ± 26 • . Since the inner kpc-scale jet is not resolved by VLA observations, the VLA position of core A is likely an average between the radio core (i.e., the EVN position) and the jet. Such a structure-based peakshift along the jet is common in AGNs (e.g. Boccardi et al. 2017;EHT MWL Science Working Group et al. 2021). Furthermore, a frequency-based core shift may contribute to a few mas-scale positional shifts (Shabala et al. 2012). The frequency-based core shift is along with the jet, so it will not affect the hypothesis of an inner jet. Though other scenarios can not be ruled out, e.g. an off-nuclear AGN or a radio bright star-forming region in the circumnuclear region of galaxy B, it seems like the jet re-orientation is self-consistent with other signatures below. We marked the direction from the EVN to VLA positions of core B as R 4 in the middle and left panels of Figure 2. It seems the bumps near R 4 in the VLA S-band image (see arrows in the middle panel of Figure  3) give some hints about the inner jet direction. Indeed, the arcsecond-scale jet direction of core A (Figure 4) is precisely consistent with the direction from the EVN to VLA coordinates (red cross in the left panel of Figure  2). The inner jet direction of core B is clearly different from the major jet axis direction on arcsec scales (see Figure 2), thus hinting at an ongoing jet re-orientation at the nuclei.
In the VLA S-band image (Figure 3), there are signatures supporting jet re-orientation in our target (see Krause et al. 2019, for a summary of signatures that supports jet re-orientation), here we marked it follow Krause et al. (2019), see Figure 3: (R) ring-like structure and wide hotspot on the northern side, and the southern hotspot is not at the edge of its lobe, (S) the northern hotspot on the opposite side of lobe/jet as the southern hotspot. The signatures of jet re-orientation were seen in the jetted kiloparsec-scale (NGC 326, Murgia et al. 2001) and parsec-scale (0402+379, Rodriguez et al. 2006) dual/binary AGN systems. More interestingly, the VLA structure of J0725+5835 is similar to 0402+379 (Rodriguez et al. 2006) and appears as a scaled-up version. Again, lateral flow signatures were found in the spectral index map of the southern hotspot (i.e. the spectral index gradient in Figure 5) which supports the ongoing jet re-orientation. Further evidence that supports jet precession is from 150 MHz image (Figure 6), where the northern jet extension is different from the VLA images.
There is no evidence, e.g. the FRI type main jets, double boomerang morphology, and bright apex in each boomerang (see Cotton et al. 2020), support the backflow is the main mechanism that is responsible for the X-shaped radio morphology. Furthermore, the major jet direction of core B is nearly aligned with the minor axis of the host (see Figure 7 in Appendix), which also disfavors the backflow model (Gillone et al. 2016). According to the above evidence of jet distortion, we favor an explanation in terms of a fast realignment of the jet axis for the X-shaped jet in J0725+5835. In addition, the backflowing plasma from the major jets maybe at work in assisting to form the 'wings'. For example, Chon et al. (2012) propose that after a re-orientation the new jet backflow is deflected into the fossil cavities created in the previous active phase.
Though the association between kpc scale binary (or dumbbell) galaxies and distorted jet structure was well addressed (Battistini et al. 1980;Wirth et al. 1982;Parma et al. 1991;Borne & Colina 1993;Murgia et al. 2001), it's unlikely in XRGs that the evolution of kiloparsec-scale galaxy pairs can drive dramatic jet realignment from wings to major jet (see Krause et al. 2019). Alternative models include a sudden spin-flip in galaxy merger (Merritt & Ekers 2002) and interaction of pcscale binary SMBHs (Liu 2004). Therefore, it is unlikely that galaxies A and B, as the pair in J0725+5835, are responsible for the X-shaped jet. In the spin-flip model, the jet axis will only change once, it can explain the quick jet realignment from R 1 (the wings) to R 2 (major jet direction). While the deficit of the spin-flip model is that the later jet re-orientation cannot be explained. Again, we will see below 4.3 that an on-going merger signature appears in galaxy B, which still disfavors the spin-flip model. On the other hand, the jet precession model can explain the ongoing jet re-orientation. However the jet precession model may require pc-scale binary SMBHs, which is not seen in our VLBI observation, a Here the TGSS contours are in steps of (0.008, 0.032, 0.128) Jy beam −1 , and the VLA contours are in steps of (1, 2, 4, ...) × 0.309 mJy beam −1 . The pseudo-color image illustrate TGSS radio flux density distribution.
possible explanation here could be its binary companion is not radio-frequency active.

Host galaxies
Taking an averaged redshift z ∼ 0.42, we can estimate black hole mass by using (a) the correlation between (5 GHz) radio power P 5 and black hole mass M BH (Lacy et al. 2001). While the L 5 − M BH correlation couples with Eddington ratio (Ho 2002), here we roughly assume an Eddington ratio of 0.1 in our case (e.g. Lacy et al. 2001). The 5 GHz radio flux density for both core A and core B is measured by taking VLA A-array Lband and B-array S-band flux density. The 5 GHz radio powers for cores A and B are P 5,A = 10 24.0 W Hz −1 and P 5,B = 10 23.7 W Hz −1 , respectively, corresponding to a black hole mass M BH ∼ 10 9 M (Lacy et al. 2001) for both galaxies A and B with an intrinsic dispersion of 1 dex due to Eddington ratio uncertainties; (b) the correlation between M BH and B-band absolute magnitude of bulge (despite the host classification, we assume B-band magnitude of our target is bulge-dominated, Graham & Scott 2013). Here the Pan-STARRS1 rband corresponding to a rest-frame B-band in our case (the effective wavelength λ ef f,B = 4400Å). Therefore, the B-band absolute magnitude (with Galactic absorption correction based on Schlafly & Finkbeiner 2011) for galaxies A and B are M r,A = −21.41 ± 0.03 and M r,B = −21.25 ± 0.05, respectively, which yield a black hole mass M BH = 10 9 − 10 10 M for both galaxies A and B; (c) the correlation between M BH and K-band absolute magnitude of bulge (assuming K-band magnitude of our target is bulge-dominated, Kormendy & Ho 2013). Here the WISE W 1-band corresponding to a rest-frame 2MASS K s -band in our case (the effective wavelength λ ef f,Ks = 2.159 µm). Therefore, the K-band absolute magnitude (with Galactic absorption corrections based on Schlafly & Finkbeiner 2011) for galaxies A and B are M K,A = −26.54 ± 0.03 and M K,B = −26.31 ± 0.04, respectively, which again yield a black hole mass M BH = 10 9 − 10 10 M for both galaxies A and B, here the M K − M BH correlation is tighter than (a) and (b) (see Kormendy & Ho 2013). The three estimates of black hole mass are consistent with each other, which induce a total galaxy stellar mass M stellar = 10 11 − 10 12 M (Reines & Volonteri 2015;Greene et al. 2020). The stellar mass is consistent with 100 kpc scale galaxy pairs (Shah et al. 2020).
The mid-infrared magnitudes for the two cores have been collected from AllWISE Data Release (Cutri et al. 2021) as m W 1,A = 15.296 ± 0.038, m W 2,A = 15.018 ± 0.075, m W 3,A = 12.200 for core A, and m W 1,B = 15.520±0.041, m W 2,B = 15.444±0.097, m W 3,B = 12.089 for core B, which yield the mid-infrared colors m W 1,A − m W 2,A = 0.27 ± 0.08 and m W 2,A − m W 3,A = 2.81 ± 0.07 for galaxy A, and m W 1,B − m W 2,B = 0.07 ± 0.10 and m W 2,B − m W 3,B = 3.35 ± 0.09 for galaxy B. Galaxy A is located in the transition region from early-type spirals (or Sa/SBa-type in Hubble's classification scheme, with semiquiescent star formation, see descriptions in Jarrett et al. 2019) to active star-forming galaxies, while galaxy B is located in the region of active star-forming galaxies (Jarrett et al. 2019). It is shown that 100 kpc scale galaxy pairs already have tidal interaction (Liu et al. 2011), while it has no help on the AGN activity (Shah et al. 2020). Therefore, the star-forming activity in galaxy B is more likely self-induced, e.g. an ongoing merger (Hopkins et al. 2006;Di Matteo et al. 2007).

SUMMARY
We perform a high-resolution EVN 5 GHz observation of the XRG J0725+5835, two compact and non-thermal radio cores are detected. We have also re-processed archival VLA 1.5 and 3 GHz data, and the radio spectral index map and astrometric coordinate have been measured in this work. Comparing with the VLA and EVN observations: (1) We find XRG J0725+5835 is associated with a 100 kpc-scale dual jetted-AGN system, where one of them is Faranoff-Riley type I (FR I), the other one is FR II; (2) We find evidence that the jet in core B is changing its direction, which naturally forms an X-shaped structure in J0725+5835. Two models are discussed on the origin of jet re-orientation, including the sudden flip of black hole spin and jet precession due to the interaction of parsec-scale binary SMBHs. Future low-frequency radio observations with high sensitivity may help to reveal the fossil radio emission; furthermore, obtaining a spectroscopic redshift and high signal-to-noise ratio image in optical is crucial to identify the proposed mechanisms in this work. This work is supported by the National Science Foundation of China (12103076) and the National Key R&D Programme of China (2016YFA0400702, 2018YFA0404602, 2018YFA0404603). XLY is supported by the Shanghai Sailing Program (21YF1455300) and China Postdoctoral Science Foundation (2021M693267). LCH was supported by the National Science Foundation of China (11721303, 11991052). XLY and TA thank the financial support of the Bureau of International Cooperation, Chinese Academy of Sciences (114231KYSB20170003). Scientific results from data presented in this publication are derived from the EVN project RSY07 and VLA project 16A-220 and 16B-023. The European VLBI Network (EVN) is a joint facility of independent European, African, Asian, and North American radio astronomy institutes. The VLBI data processing in this work made use of the compute resource of the China SKA Regional Centre prototype, funded by the Ministry of Science and Technology of China and the Chinese Academy of Sciences. This work made use of Karl G. Jansky Very Large Array (VLA) data. The National Radio Astronomy Observatory is a facility of the U.S. National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This work made use of TGSS data. We thank the staff of the GMRT that made these observations possible. GMRT is run by the National Centre for Radio Astrophysics of the Tata Institute of Fundamental Research. This research has made use of the NASA/IPAC Infrared Science Archive, which is funded by the National Aeronautics and Space Administration and operated by the California Institute of Technology. The Photometric Redshifts for the Legacy Surveys (PRLS) catalog used in this paper was produced thanks to funding from the U.S. Department of Energy Office of Science, Office of High Energy Physics via grant DE-SC0007914. This work made use of data from the Pan-STARRS1 Surveys (PS1).    Figure 7). The histograms on the diagonal show the marginalized posterior densities for each parameter, i.e. position angles (θ, in radians) and standard deviation (in arcsec) of the major (σY ) and minor (σX ). We take 50% of the distributions as the best-fit values, which are marked in red cross-hairs and red lines. The uncertainties are computed as the 16% and 84% of the distributions, thus representing 1σ confidence ranges shown with black dashed lines. the standard deviations (∼ F W HM 2.355 ) of the major and minor axes. We take 16% and 84% of the distributions as the lower and upper limits, respectively.