Collimation of the Relativistic Jet in the Quasar 3C 273

Collimation of the Relativistic Jet in the Quasar 3C 273 Hiroki Okino 1,2 , Kazunori Akiyama 3,4,5 , Keiichi Asada 6 , José L. Gómez 7 , Kazuhiro Hada 5,8 , Mareki Honma 5,8,1 , Thomas P. Krichbaum 9 , Motoki Kino 10,5 , Hiroshi Nagai 2,8 , Uwe Bach 9 , Lindy Blackburn 11,4 , Katherine L. Bouman 12 , Andrew Chael 13,30 , Geoffrey B. Crew 3 , Sheperd S. Doeleman 11,4 , Vincent L. Fish 3 , Ciriaco Goddi 16,14,15,17 , Sara Issaoun 11,30 , Michael D. Johnson 11,4 , Svetlana Jorstad 18 , Shoko Koyama 19,6 , Colin J. Lonsdale 3 , Ru-Sen Lu 20,21,9 , Ivan Martí-Vidal 22,23 , Lynn D. Matthews 3 , Yosuke Mizuno 29 , Kotaro Moriyama 24,5 , Masanori Nakamura 25,6 , Hung-Yi Pu 26 , Eduardo Ros 9 , Tuomas Savolainen 27,28,9 , Fumie Tazaki 5 , Jan Wagner 9 , Maciek Wielgus 9 , and Anton Zensus 9 1 Department of Astronomy, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan; [email protected] 2 National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan 3 Massachusetts Institute of Technology Haystack Observatory, 99 Millstone Road, Westford, MA 01886, USA 4 Black Hole Initiative, Harvard University, 20 Garden Street, Cambridge, MA 02138, USA 5 Mizusawa VLBI Observatory, National Astronomical Observatory of Japan, 2-12 Hoshigaoka, Mizusawa, Oshu, Iwate 023-0861, Japan 6 Institute of Astronomy and Astrophysics, Academia Sinica, 11F of Astronomy–Mathematics Building, AS/NTU No. 1, Sec. 4, Roosevelt Rd, Taipei 10617, Taiwan, R.O.C. 7 Instituto de Astrofísica de Andalucía-CSIC, Glorieta de la Astronomía s/n, E-18008 Granada, Spain 8 Department of Statistical Science, The Graduate University for Advanced Studies (SOKENDAI), 10-3 Midori-cho, Tachikawa, Tokyo 190-8562, Japan 9 Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany 10 Kogakuin University of Technology & Engineering, Academic Support Center, 2665-1 Nakano, Hachioji, Tokyo 192-0015, Japan 11 Center for Astrophysics |Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA 12 California Institute of Technology, 1200 East California Boulevard, Pasadena, CA 91125, USA 13 Princeton Gravity Intiative, Jadwin Hall, Princeton University, Princeton, NJ 08544, USA 14 Dipartimento di Fisica, Universit degli Studi di Cagliari, SP Monserrato-Sestu km 0.7, I-09042 Monserrato, Italy 15 INAF—Osservatorio Astronomico di Cagliari, via della Scienza 5, I-09047 Selargius (CA), Italy 16 Universidade de São Paulo, Instituto de Astronomia, Geofísica e Ciências Atmosféricas, Departamento de Astronomia, São Paulo, SP 05508-090, Brazil 17 INFN, Sezione di Cagliari, Cittadella Univ., I-09042 Monserrato (CA), Italy 18 Institute for Astrophysical Research, Boston University, 725 Commonwealth Avenue, Boston, MA 02215, USA 19 Niigata University, 8050 Ikarashi 2-no-cho, Nishi-ku, Niigata 950-2181, Japan 20 Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan Road, Shanghai 200030, People’s Republic of China 21 Key Laboratory of Radio Astronomy, Chinese Academy of Sciences, Nanjing 210008, People’s Republic of China 22 Departament d’Astronomia i Astrofísica, Universitat de València, C. Dr.Moliner 50, E-46100 Burjassot, Spain 23 Observatori Astronòmic, Universitat de València, C. Catedrático Beltrán 2, E-46980 Paterna, Spain 24 Institut für Theoretische Physik, Goethe-Universität Frankfurt, Max-von-Laue-Straße 1, D-60438 Frankfurt am Main, Germany 25 National Institute of Technology, Hachinohe College, 16-1 Uwanotai, Tamonoki, Hachinohe City, Aomori 039-1192, Japan 26 Department of Physics, National Taiwan Normal University, No. 88, Sec. 4, Tingzhou Rd., Taipei 116, Taiwan, R.O.C. 27 Aalto University Department of Electronics and Nanoengineering, PL 15500, FI-00076 Aalto, Finland 28 Aalto University Metsähovi Radio Observatory, Metsähovintie 114, FI-02540 Kylmälä, Finland 29 Tsung-Dao Lee Institute and School of Physics and Astronomy, Shanghai Jiao Tong University, 800 Dongchuan Road, Shanghai, 200240, People’s Republic of China Received 2021 December 22; revised 2022 October 4; accepted 2022 October 4; published 2022 November 22 Abstract The collimation of relativistic jets launched from the vicinity of supermassive black holes (SMBHs)at the centers of active galactic nuclei (AGNs)is one of the key questions to understand the nature of AGN jets. However, little is known about the detailed jet structure for AGN like quasars since very high angular resolutions are required to resolve these objects. We present very long baseline interferometry (VLBI)observations of the archetypical quasar 3C 273 at 86 GHz, performed with the Global Millimeter VLBI Array, for the first time including the Atacama Large Millimeter/submillimeter Array. Our observations achieve a high angular resolution down to ∼60 μas, resolving the innermost part of the jet ever on scales of ∼10 5 Schwarzschild radii. Our observations, including close-in-time High Sensitivity Array observations of 3C 273 at 15, 22, and 43 GHz, suggest that the inner jet collimates parabolically, while the outer jet expands conically, similar to jets from other nearby low-luminosity AGNs. We discovered the jet collimation break around 10 7 Schwarzschild radii, providing the first compelling evidence for structural transition in a quasar jet. The location of the collimation break for 3C 273 is farther downstream from the sphere of gravitational influence (SGI)from the central SMBH. With the results for other AGN jets, our results show that the end of the collimation zone in AGN jets is governed not only by the SGI of the SMBH but also by the more diverse properties of the central nuclei. Unified Astronomy Thesaurus concepts: Active galactic nuclei (2035);Blazars (164);Quasars (1319);Relativistic jets (1390);Very long baseline interferometry (1769) 1. Introduction Relativistic jets ejected from active galactic nuclei (AGNs) are tightly collimated plasma outflows from galactic centers, known as the most energetic persistent phenomena in the universe. Their formation, acceleration, and collimation mechanisms constitute one of the most important questions in The Astrophysical Journal, 940:65 (17pp), 2022 November 20 https://doi.org/10.3847/1538-4357/ac97e5 © 2022. The Author(s). Published by the American Astronomical Society. 30 NASA Hubble Fellowship Program, Einstein Fellow. 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. 1 AGN jet physics (e.g., reviewed in Blandford et al. 2019).In particular, jet collimation processes have been studied by focusing on jet structures since they provide various information not only about the jet itself but also about their surrounding environment contributing to jet confinement (Begelman & Li 1994; Komissarov et al. 2007,2009; Fromm et al. 2018). The detailed nature of the jet collimation can be addressed by high-angular-resolution observations using very long baseline interferometry (VLBI; e.g., Boccardi et al. 2017; Hada 2019), allowing direct comparisons with special or general relativistic (SR or GR)magnetohydrodynamic simulations (Chael et al. 2018; Nakamura et al. 2018; Fromm et al. 2019). In particular, intensive studies of the M87 jet in the last decade have provided new insights about jet collimation since that jet is the best source for investigating the global jet structures from the immediate vicinity of the central black hole (BH)to beyond the host galaxy (Asada & Nakamura 2012; Hada et al. 2013,2016; Nakamura & Asada 2013; Kim et al. 2018; Walker et al. 2018; Event Horizon Telescope Collaboration et al. 2019a; Park et al. 2019a). Following the aforementioned pioneering works on M87, jet shapes have been measured mainly from several nearby sources: NGC 6251 (Tseng et al. 2016), NGC 4261 (Nakahara et al. 2018), 1H 0323+342 (Hada et al. 2018), Cygnus A (Boccardi et al. 2016; Nakahara et al. 2019), NGC 1052 (Nakahara et al. 2020; Baczko et al. 2022), NGC 315 (Boccardi et al. 2021; Park et al. 2021),3C84(Nagai et al. 2014; Giovannini et al. 2018), and for a large number of jet samples collected by the long-term monitoring MOJAVE program 31 (Pushkarev et al. 2017). More recently, Kovalev et al. (2020)investigated jet widths from over 300 sources and found the transitions of jet shapes from 10 nearby sources. They also found that the range of the transition locations is ∼10 5–6 gravitational radii from the core, which may indicate a common property of nearby AGN jets. The collimation properties of distant quasars are still unclear due to the difficulties in resolving the transverse jet shape. Indeed, Algaba et al. (2017)systematically measured core sizes for distant radio-loud AGNs to investigate the upstream jet structures. However, the core size is often unresolved even at an extremely high angular resolution with a space radio dish (e.g., Gómez et al. 2016), and is often measured with high systematic uncertainties resulting from the limited angular resolution. Although the jet structure of the flat-spectrum radio quasar (FSRQ)4C 38.41 was measured by Algaba et al. (2019),they could not detect significant collimation breaks because of the limited range of the measured jet scales. In addition, the physical parameters for this source, such as the BH mass and the jet viewing angle, have large uncertainties. Therefore, to understand the jet collimation for high-powered AGN like quasars, observations with a high angular resolution for the well-studied sources at other wavelengths are required to transversely resolve the detailed structures down to the bases of the jets. In this paper, we study the quasar 3C 273 (1226+023),well known as one of the brightest extragalactic objects and the closest quasars (Schmidt 1963)with prominent jets (Davis et al. 1985; Conway et al. 1993; Bahcall et al. 1995; Jester et al. 2005; Perley & Meisenheimer 2017). The 3C 273 jets have been observed many times since their discovery because the elongated jet is a unique target for resolving down to the central (sub)parsec scale with VLBI observations (Krichbaum et al. 1990; Lobanov & Zensus 2001; Savolainen et al. 2006; Kovalev et al. 2016; Bruni et al. 2017,2021; Jorstad et al. 2017; Lister et al. 2019,2021). Furthermore, recent infrared interferometric observations with the GRAVITY instrument on the Very Large Telescope Interferometer (Gravity Collaboration et al. 2017) precisely estimated the BH mass and viewing angle for 3C 273 as M BH =(2.6 ±1.1)×10 8 M e and θ=12°±2°(Gravity Collaboration et al. 2018). 32 The estimated BH mass yields a linear angular relation of 1 mas ∼2.7 pc ∼1.2 ×10 5 R s , which makes 3C 273 one of the best-resolved quasars. For these reasons, 3C 273 is the ideal target for investigating the global structure of a quasar jet. For 3C 273, a preliminary detection of the transition from the parabolic to conical shape has been reported in an early work of a subset of the authors in a conference proceeding (Akiyama et al. 2018), based on a marginal detection of a parabolic jet from single-band data sets at 43 GHz only covering a narrow range of the spatial area. It also lacks a careful consideration of frequency-dependent locations of radio cores at different frequencies due to the effect of synchrotron self-absorption relevant for the innermost jet probed at millimeter wavelengths. The presence of the transition has remained inconclusive in the earlier work. In this paper, we report new multifrequency observations of the jet of the quasar 3C 273 observed with several global VLBI networks. In particular, the Global Millimeter VLBI Array (GMVA)at 3.5 mm/86 GHz including the phased Atacama Large Millimeter/submillimeter Array (ALMA; Matthews et al. 2018)provided a strong advantage, remarkably increasing the north–south resolution and sensitivity (see also Issaoun et al. 2019 for the first published image with GMVA+ALMA on Sgr A * ). These observations provided us with high-fidelity imaging of the finest jet structure of 3C 273 from the innermost subparsec scale region with a maximum resolution of several tens of microarcseconds (μas). The paper is organized as follows. In Section 2, we describe the observations and data reduction. Images were reconstructed with the state-of-the-art regularized maximum likelihood methods as described in Section 3, and further analyzed in Section 4. We show the global jet structure of 3C 273, including the core shift measurements and the jet collimation profile in Section 5. The physical implications of our observations are discussed in Section 6. Finally, we summarize our results in Section 7.Inthis paper, we assume a flat ΛCDM cosmology (e.g., Planck Collaboration et al. 2014)with H 0 =70 km s −1 Mpc −1 ,Ω M =0.3, and Ω Λ =0.7, and we adopt values of M BH =2.6 ×10 8 M e and a viewing angle of θ=9° 33 for the 3C 273 jet. 2. Observations and Data Analysis 2.1. GMVA 86 GHz We observed 3C 273 with the GMVA at 86 GHz (λ∼3.5 mm)on 2017 April 3 (project code MA008), as one of the first 31 https://www.physics.purdue.edu/MOJAVE/index.html 32 Li et al. (2022)reported a new estimation of BH mass of ∼10 9 M e and viewing angle of ∼5°for 3C 273 during the review process of our paper. These new estimates are based on a joint analysis of spectroastrometric data in Gravity Collaboration et al. (2018)and new data from reverberation mapping with a more generalized model. We confirm that these estimates do not affect our main results and conclusions. 33 Previous studies have reported various values for the viewing angle of the 3C 273 jet. Jorstad et al. (2017)reported θ∼6°from VLBI monitoring at 43 GHz. Meyer et al. (2016)showed the possible range of 3°.8–7°.2. However, other studies have reported larger values of θ10°(Savolainen et al. 2006; Gravity Collaboration et al. 2018). Therefore, we assume θ=9°, which is broadly consistent with previously reported values of jet viewing angles. 2 The Astrophysical Journal, 940:65 (17pp), 2022 November 20 Okino et al. VLBI observations with ALMA (see Goddi et al. 2019, for a detailed description). The phased ALMA, eight VLBA stations, and four European stations participated in our observation, as summarized in Table 1. 3C 273 and calibrators (3C 279)were observed for a track of ∼16 hr, almost a full track for the given observing array. ALMA participated for ∼5.2 hr in the middle of the track, providing baselines to both VLBA and European stations. Data were recorded at a total bandwidth of 256 MHz per polarization, which was further subdivided into four 58 MHz intermediate frequencies (IFs)of 116 channels each. Data were correlated using the DiFX correlator (Deller et al. 2011)at the Max Planck Institute for Radio Astronomy in Bonn, Germany. We note that Maunakea (MK)did not provide robust fringe detections on any of the observed sources because of the bad weather at MK throughout the GMVA campaign during spring 2017 (see also Issaoun et al. 2019). We also note that the visibility phases of PV baselines at each subband IF had an instrumental offset of 180°in subsequent channels of 32 MHz widths owing to a misconfiguration of the subband alignment in the correlation stage. Initial data calibration was performed using the Astronomical Image Processing System (AIPS; Greisen 2003). Preceding standard calibration, the inner-IF phase offsets in the PV baselines were corrected using a manually created bandpass table. Visibility amplitudes were a priori calibrated in the standard manner; visibilities at each baseline were first normalized with the available autocorrelation spectra and then scaled with the system-equivalent flux densities of the corresponding stations. Phases were calibrated in a standard manner with multiple fringe-fitting runs after parallactic angle correction. First, the instrumental phases stable across the track, such as phase bandpass, inter-IF phase, and delay offsets, were corrected with a scan of 3C 273, providing strong detection to all stations except MK. This enabled coherent integration across the entire bandwidth and more sensitive fringe-fitting by combining all IFs. Second, the delay and rate offsets of each scan were calibrated at the solution interval of the scan length, which is typically several minutes. Then, short-timescale phase rotations were fringe-fitted at solution intervals of 10 s. The new addition of ALMA to GMVA has significantly improved the overall array performance. Figure 1shows the uv coverage for our GMVA observation of 3C 273. The long baselines beyond ∼1.5 Gλin north–south (N–S)direction correspond to ALMA. The phased ALMA has improved the angular resolution of the GMVA observations in the N–S direction by more than a factor of two (see Section 5.1). Furthermore, ALMA provides detection at very high signal-tonoise ratio on long baselines to both VLBA and European stations at ∼1.5–2.0 Gλ(Figure 2). The sensitive ALMA secured the detection of fringes to most of the stations in the array while it was participating in the observation. 2.2. HSA 15/22/43 GHz We performed multifrequency observations of 3C 273 at 15, 22, and 43 GHz with the High Sensitivity Array (HSA)on 2017 March 26, eight days before our GMVA observations (project code: BA122). The observing array consisted of the 10 VLBA stations and the 100 m Effelsberg (EB)telescope in Germany, as summarized in Table 1. Full-track observations were performed for ∼11 hr, of which EB participated in the first ∼3.4 hr. Data were recorded at dual circular polarizations with four 64 MHz IFs subdivided into 256 channels, providing Table 1 Summary of Observations of 3C 273 Array Project Code Stations Frequency Obs. Date Beam Size Geometric Mean (GHz)(yyyy/mm/dd)(mas ×mas, deg)(mas) (1)(2)(3)(4)(5) VLBA BH151 VLBA(10)1.667 2008/02/02 10.312 ×4.460, −4.3 6.778 HSA BA122 EB+VLBA(10)15.368 2017/03/26 1.025 ×0.408, −6.3 0.647 HSA BA122 EB+VLBA(10)23.768 2017/03/26 0.577 ×0.273, −3.1 0.397 HSA BA122 EB+VLBA(10)43.168 2017/03/26 0.325 ×0.128, −6.1 0.204 GMVA MA008 AA+EB+ON+PV+YS+VLBA(8)86.268 2017/04/03 0.061 ×0.052, −42.7 0.057 Note. (1)The stations that participated in the observations. The two-letter codes for each station are as follows: ALMA (AA), Effelsberg (EB), Onsala (ON), Pico Veleta (PV), Yebes (YS), Brewster (BR), Fort Davis (FD), Hancock (HN), Kitt Peak (KP), Los Alamos (LA), Maunakea (MK), North Liberty (NL), Owens Valley (OV), Pie Town (PT), and St. Croix (SC). Ten VLBA stations (BR, FD, HN, KP, LA, MK, NL, OV, PT, and SC)are included in VLBA(10), and eight stations, excluding HN and SC, in VLBA(8).(2)The observed central frequency. (3)The observing date. (4)Sizes of the major/minor axes and position angles of the synthesized beam with uniform weighting. (5)The beam size of an equivalent circular Gaussian with the same beam solid angle derived from the geometric mean of the majorand minor-axis beam sizes. Figure 1. uv coverage of the GMVA observation at 86 GHz. Baselines to ALMA are shown in red. Data are averaged at 60 s. 3 The Astrophysical Journal, 940:65 (17pp), 2022 November 20 Okino et al. an aggregate of 256 MHz per polarization. Data were correlated using the VLBA DiFX correlator. The initial data calibration for the HSA data was performed using AIPS. We note that the EB baselines in several IFs showed lower visibility amplitudes and larger scattering in visibility phases due to the bandpass issue, which cannot be calibrated in this stage. These data were removed prior to the data calibration. Data were calibrated in a standard manner similar to the GMVA data sets described in Section 2.1; amplitudes were a priori calibrated, and then phases were calibrated with fringe-fitting. Fringes were detected on all stations, and similar baseline coverages were obtained at all three observing frequencies. 2.3. VLBA 1.7 GHz To complement our data on the larger-scale jet structure, we reduced an archival VLBA 1.7 GHz data set of 3C 273 observed in 2008 February (project code: BH151). The data set has one of the best hour-angle coverage among existing archival data sets at 1.7 GHz. 34 Data were reduced in a standard manner in AIPS. Fringes were detected on all 10 VLBA stations, as summarized in Table 1. 3. Imaging We reconstructed multifrequency images of 3C 273 from the GMVA, HSA, and VLBA data sets presented in Section 2with regularized maximum likelihood (RML)methods implemented in SMILI (Akiyama et al. 2017a,2017b). Our GMVA and HSA arrays are heterogeneous and were expected to have more residual calibration errors at higher frequencies, such as 43 and 86 GHz. Furthermore, our GMVA observations had sparse uv coverages lacking intermediate baselines as shown in Figure 1. The traditional iterative hybrid mapping using CLEAN (e.g., Högbom 1974)with self-calibration implemented in popular packages such as AIPS and Difmap (Shepherd 1997)is generally more challenging than for the lower-frequency data sets and does not offer the flexibility. RML methods provide a more flexible imaging framework, by directly using various types of data sets including robust closure quantities (e.g., Thompson et al. 2017)and a wider range of assumptions for the source images. We employed an imaging approach inspired by recent Event Horizon Telescope (EHT)imaging of M87 (Event Horizon Telescope Collaboration et al. 2019b)that explored a wide range of parameters, implying assumptions of source images, and assessed uncertainties in images and their deliverables more conservatively. Using a scripted RML imaging pipeline, we explored the distribution of images giving a satisfactory fit to the data, which were all further used in the image analysis in Section 4. We briefly introduce the RML method in Section 3.1. Then, we describe the imaging pipeline used in our imaging of 3C 273 in Section 3.2, and the details of the imaging parameter survey in Section 3.3. 3.1. RML Methods RML methods are a new class of imaging techniques that were conceived to overcome many of the technical challenges of millimeter VLBI imaging, particularly with the EHT. RML imaging takes a forward-modeling approach inspired by Bayesian statistics, directly solving for an image without using a dirty beam or dirty map. RML methods have been demonstrated to improve the overall quality of image reconstruction not only for synthetic observations (e.g., Honma et al. 2014; Chael et al. 2016,2018; Akiyama et al. 2017a,2017b; Kuramochi et al. 2018)but also for actual interferometric measurements with VLBI arrays (e.g., Event Horizon Telescope Collaboration et al. 2019b; Issaoun et al. 2019; Kim et al. 2020; Janssen et al. 2021)and connected interferometers (e.g., Matthews et al. 2018; Yamaguchi et al. 2020). RML methods derive a reasonable or conservative image from an infinite number of images consistent with given interferometric measurements by solving for an image that minimizes the sum of data consistency metrics, such as χ 2 terms, and regularization functions that mathematically describe prior assumptions for the source morphology. This framework allows a flexible choice of input data, for instance, the direct use of closure quantities free from antenna-based calibration errors (e.g., Chael et al. 2016,2018; Akiyama et al. 2017b; Blackburn et al. 2020)and further inclusion of various observing effects such as systematic non-closing errors (Event Horizon Telescope Collaboration et al. 2019b). Popular regularization functions include ℓ 1 -norm, total variation (TV), and total squared variation (TSV), enforcing sparsity in some basis of the image (Honma et al. 2014; Ikeda et al. 2016; Akiyama et al. 2017a,2017b; Kuramochi et al. 2018)and the information entropy of the image (e.g., Chael et al. 2016). By combining various regularization functions, RML methods can explore a wide range of images consistent with the data, often leading to reconstruction with more reasonable assumptions of the target source (e.g., Akiyama et al. 2017a,2017b; Event Horizon Telescope Collaboration et al. 2019b). Further technical and mathematical details of the RML approaches adopted in this study are described in Event Horizon Telescope Collaboration et al. (2019b). 3.2. Imaging Pipeline We used a scripted pipeline in Python for our RML imaging with SMILI. In the pipeline, 3C 273 images were reconstructed from AIPS-calibrated data (see Section 2)by utilizing Figure 2. The signal-to-noise ratios of visibilities on 3C 273 as a function of projected baseline length. Data are averaged at 60 s. Baselines to ALMA are shown in red. 34 We cross-checked archival data of a more close-in-time VLBA observation at 1.6 GHz in 2014 with the experiment code of BG216H. We obtained a jet width profile consistent with the data sets used in this paper, confirming that the width of the 1.6 GHz jet does not change significantly over the timescale of 10 yr. 4 The Astrophysical Journal, 940:65 (17pp), 2022 November 20 Okino et al. weighted-ℓ 1 , TV, TSV, and maximum entropy method (MEM) regularizers (see Appendix A of Event Horizon Telescope Collaboration et al. (2019b)for mathematical definitions). The pixel size of the image is typically set to one-tenth of the mean FWHM size of the uniform-weighted synthesized beam summarized in Table 1. Prior to imaging, the uv coordinates of the data were rotated by 45°counterclockwise, resulting in the same rotation of the image-domain axis, and the horizontal axis approximately aligned with the jet direction. This allows us to use a narrower rectangular field of view of the image and minimizes the computational cost of imaging. Then, the visibilities were coherently averaged for 60 s, after manually flagging outliers. If specified, systematic errors were added in quadrature to the thermal noise of the time-averaged complex visibilities to account for non-closing errors. We explored 0 or 1% of non-closing errors within the range expected for polarization leakages (e.g., Zhao et al. 2022). Our imaging procedure is iterative, with four stages of imaging and self-calibration. In the first stage, the initial image is set to a circular Gaussian with the size corresponding to the geometric mean of the major and minor-axis sizes of the synthesized beam at each observed frequency (see Table 1). The subsequent stages of imaging started with the final image of the previous stage convolved with the above circular Gaussian. The first two stages begin with the initial image and use visibility amplitudes, log closure amplitudes, and closure phases for imaging. At each stage, the visibility amplitudes of all baselines added a 5% error in quadrature to the thermal noise to account for the amplitude calibration uncertainties. In the final two stages, the imaging uses complex visibilities and closure quantities. 3.3. Imaging Parameter Survey Using the pipeline described in Section 3.2, we explored a wide range of parameters, as summarized in Table 2. At each frequency, the parameter survey examined several tens to a few hundred parameter sets with four to five parameters: potential systematic non-closing errors (denoted as Sys. Error)and the regularization parameters for the ℓ 1 w , TV, TSV, and MEM regularizers. Among the sets of imaging parameters explored, we selected “top sets”of the parameters that provide reconstructed images with reasonable fits to data. The distributions of the corresponding top-set images allow us to assess and identify the morphologies commonly seen and insensitive to the imaging choices as well as their uncertainties among the reconstructed images consistent with our interferometric measurements. We started imaging with the HSA data at 15 GHz. We adopted the published VLBA image at 15 GHz and Stokes Ion 2017 May 25 (project code BL229AH)from the MOJAVE program (Lister et al. 2018), which was the closest epoch to our HSA observations, as the soft mask of the imaging region with ℓ 1 w regularization. This allows noise suppression outside of the area, which has historically had no significant emission. 35 Then, the best-fit image among the top-set 15 GHz images is adopted as the soft mask for subsequent 22 GHz image processing to maximize the consistency between adjacent frequencies. We processed data at higher frequencies (43 and 86 GHz)in the same manner using best-fit images at the adjacent lower frequencies. For L-band data, we used the Lband VLBA image in Akiyama et al. (2018)as its soft mask. We selected the top-set images with good fits to the data using χ 2 statistics. For the selection of top-set parameters, we adopted the minimum threshold of 1.5 for the reduced χ 2 values of full complex visibilities, amplitudes, and closure quantities of self-calibrated images. Consequently, 15 to 45 images were selected at each frequency, as summarized in Table 2, and used in post-imaging analysis. 3.4. Total Flux Scaling for GMVA 86 GHz The overall scaling of visibility amplitudes depends strongly on the accuracy of the a priori calibration of visibility amplitudes based on measurements of the system-equivalent flux density at each station, which often have large systematic errors in high-frequency VLBI observations above ∼86 GHz (see, e.g., Koyama et al. 2016; Kim et al. 2018, for previous GMVA observations). We assumed Effelsberg (EB)and Pico Veleta (PV), often considered as reliable stations for a priori calibrations, as the reference stations for overall gain scaling (see, e.g., Angelakis et al. 2015; Fuhrmann et al. 2016; Agudo et al. 2018). We scaled the total flux density of our GMVA 86 GHz images such that the median gain amplitude of EB and PV is unity. The derived scaling factor is ∼1.46, and we obtained the scaled total flux (median in the top set)as ∼3.1 Jy. The scaling of amplitudes and the total flux density described here will not affect our main results based on the collimation profile measured from the source morphology (Section 4.2). We note that our 86 GHz images show reasonable spectral indices to images at lower frequencies (see Section 5.3). 4. Image Analysis With all top-set images from the imaging parameter survey (Section 3.3), we investigate the jet profile as a function of distance from the central black hole. First, we identified the position of the central black hole using the core shift effect as described in Section 4.1. Then, the jet radii at each frequency were measured, as described in Section 4.2. Table 2 List of Imaging Parameters Freq. Sys. Error ℓ 1 w TV TSV MEM Total No. of No. of Top-set (GHz)(%)Parameters Parameters 1.6 0, 1 1 1 0, 10 −2 ,10 −1 ,1,10 0,10 −2 ,10 −1 , 1 40 15 15 0, 1 10 −1 , 1, 10 1 0, 10 −2 ,10 −1 ,1,10 0,10 −2 ,10 −1 , 1, 10 150 45 22 0, 1 10 −1 , 1, 10 1 0, 10 −2 ,10 −1 ,1,10 0,10 −2 ,10 −1 , 1, 10 150 30 43 0, 1 10 −1 , 1, 10 1 0, 10 −2 ,10 −1 ,1,10 0,10 −2 ,10 −1 , 1, 10 150 35 86 0, 1 10 −1 ,110,10 −2 ,10 −1 , 1 0, 10 −2 ,10 −1 , 1 64 24 35 The typical jet proper motion of 3C 273 is a few mas yr −1 , which does not change the overall jet emission structure in a few months. 5 The Astrophysical Journal, 940:65 (17pp), 2022 November 20 Okino et al. 4.1. Core Shift Measurements Core shift is a positional shift of the radio core between two different observing frequencies due to the frequency dependence of the optical depth of the synchrotron self-absorption (Blandford & Königl 1979). In our study, we employed the widely used self-referencing method (e.g., Lobanov 1998; O’Sullivan & Gabuzda 2009; Pushkarev et al. 2012; Fromm et al. 2013; Hada et al. 2018)to measure the core shift of 3C 273, based on a well-vetted assumption that the emitting regions from the extended jet are optically thin features, and their positions do not change at different frequencies. From this assumption, one can obtain the core offset values between two images at different frequencies after aligning the corresponding optically thin emission regions. We only used the images at 15, 22, 43, and 86 GHz observed within a week (see Table 1). The offsets were measured for four combinations (15/22, 22/43, 43/86, 15/43 GHz). 36 For each frequency pair, we derived the offsets using all combinations of the top-set images and adopted their mean values as the measured offsets. First, we needed to identify the core position for each jet image to compare them at different frequencies. We decided on the locations of some components in the upstream regions by several circular Gaussian model fittings using Difmap.We applied this analysis to all the top-set images and identified four to seven components for the 15–86 GHz images. We set the core at the most upstream component for each image. To obtain the core offset between two frequencies, we used a two-dimensional cross-correlation of the optically thin emission regions in the jet images following the widely used method described in Croke & Gabuzda (2008). First, the pixel size of both images was set to a twentieth of the restoring beam size at the lower frequency, which was much smaller than the angular resolution at both frequencies. Next, we convolved both images with the restoring beam at the lower frequency. After masking the optically thick core regions of both images, we computed the cross-correlation coefficient r xy defined by (¯)( ¯) (¯)( ¯) ()=SS - - SS - SS - nnnn nn nn == == == rIIII II II 1 xy i n j nij ij i n j nij i n j nij 11 1, 1 2, 2 11 1, 1 211 2, 2 2 where I ν,ij is the intensity at pixel (i,j)at frequency ν, and ¯ I is the averaged intensity in the calculated region. The cross-correlation coefficient r xy is a function of the relative positional offset (Δx,Δy). The positional shift between two images at frequencies ν 1 and ν 2 is given by the location of the global maximum r xy . The core offset was defined as the remaining positional offset of the peaks after aligning the corresponding optically thin emission regions. 4.2. Jet Profile Measurements We investigated the jet collimation in 3C 273 using multifrequency images. Prior to analysis, we cut off the brightness distribution for all the images below the rms values (see Section 5.1)at each frequency to remove the noise outside the jet regions. Then, the jet radius was measured as a function of distance from the core at each frequency in two steps. First, we measured the position angle (PA)of the jet at each distance from the core. Following Pushkarev et al. (2017),we took a circular slice of the image centered at the core with a radius corresponding to each distance and adopted the intensity-weighted centroid as the PA of the jet at the corresponding distance. After deriving the PA, we measured the jet radius at each distance. Because the jet bends on milliarcsecond scales, we examined two ways of slicing the image to measure the jet radius. The first one is to take a slice perpendicular to the PA at each distance, and the other is to slice the image perpendicular to the local tangent line of the PA profile, following Pushkarev et al. (2017). The effective jet radius was measured by taking the second moment σ w of the cross section along with each slice at each distance. Then, the FWHM of the jet was obtained as s= F WHM 2 2ln2 w. We derived the deconvolved FWHM (w)defined as q=- w FWHM2beam 2, correcting the blurring effects with the restoring beam at the FWHM of θ beam 37 (see Table 1). Finally, the jet radius rat each distance was obtained from half the size of the deconvolved FWHM (i.e., r=w/2). The jet radius profile was measured with an interval of half of the restoring beam size (see Table 1).Weonlyacceptmeasurements of jet radius that satisfy the following three criteria: (i) the FWHM is larger than θ beam so that the vertical jet emission is well resolved, (ii)the distance from the radio core is greater than 2×θ beam to conservatively remove the effect of the unresolved core, and (iii)the peak intensity of the slice is three times larger than the mean of the residual rms noise derived from all top-set images (see Section 5)at each frequency except for 22 GHz. At 22 GHz, instead, the threshold was set to be nine times larger than the mean residual rms noise to make reliable measurements only at bright locations in the sparse intensity distribution (see the image in Figure 5). The jet radius profile at each frequency was measured for all top-set images with the above two slicing methods. Then, the mean and standard deviation of the jet radius at each distance and frequency were adopted as the corresponding jet radius and its uncertainty, respectively. 5. Results 5.1. GMVA Image In Figure 3, we show the mean total intensity image of the 3C 273 jet with GMVA at 86 GHz derived by averaging all top-set reconstructions (24 images)obtained from the imaging parameter survey in Section 3.3. Here we show the mean image restored at two characteristic resolutions of our GMVA observations equivalent to uniform weighting. 38 Our new image provides the sharpest view of the inner jet at 86 GHz thanks to the critical addition of ALMA. Without ALMA, the corresponding beam sizes will be 265 ×50 μas in uniform weighting, providing significantly worse resolutions along the N–S direction by a factor greater than 2.5. The spatial resolution of our observations is also higher than that of previous 36 The other two (i.e., 15/86 and 22/86 GHz)pairs were omitted due to the lack of extended emission commonly seen at both frequencies. 37 Although interferometric imaging with RML, in principle, does not require a post-imaging beam convolution to obtain piecewise smooth images, the raw image may contain some over-resolved features that are not strongly supported by data. To reduce potential bias from such features, in this work, we adopt the conventional approach that derives deconvolved widths from ones measured from beam-convolved images widely used with CLEAN imaging. 38 The restoring beam is arbitrary for RML methods, given that the nature of methods is such that the synthesized beam is not used for deconvolution. Data weighting is equivalent to the natural weighting of CLEAN imaging since data are only weighted by thermal noise and not by uv density, while the reconstructed images are known to have a high fidelity at the resolution comparable to, or higher than, uniform weighting. 6 The Astrophysical Journal, 940:65 (17pp), 2022 November 20 Okino et al. observations, such as the high-sensitivity VLBA +GBT observations in Hada et al. (2016)with a resolution of 340 ×110 μas at PA =−10°. The high resolution of our new observation allows us to resolve the vertical structure of the inner jet within a few milliarcseconds from the core at 86 GHz for the first time. To assess the image uncertainty quantitatively, we obtained residual images using the top-set reconstructions (not blurred) as the image models in the Difmap software. We calculated the image rms values for each residual in the top set, and derived the mean value of ∼3 mJy beam −1 as the representative image uncertainty at 86 GHz with a natural weighting scheme. We then determined the lowest contour level for the uniformweighting beam to be the same dynamic range (∼100)as the one for natural weighting. In addition, the final top-set reconstructions have some variations in the image domain depending on the imaging parameters. Figure 4shows the standard deviation image of GMVA 86 GHz across the top-set reconstructions. The brightest core shows the largest emission (∼200 mJy beam −1 ), while the faint and extended emissions show smaller variation. These results suggest that the variability in the image domain, especially in the bright region, is dominated by the choice of imaging parameter rather than the rms of the residual images. The median value of the scaled total flux densities among the top sets of 86 GHz images is ∼3.1 Jy (Section 3.4), corresponding to ∼30% of the arcsecond-scale core flux density measured with ALMA in Band 3 during this GMVA campaign (Goddi et al. 2019). The extended jet emission at 3 mm is broadly consistent with the 7 mm jet, implying optically thin spectra (see Section 5.3)and suggesting that the rest of the ALMA core flux density has been resolved out at the shortest VLBI baselines with fringe spacings equivalent to a few milliarcseconds. 5.2. HSA and VLBA Images In Figure 5, we show lower-frequency images of 3C 273 observed with the HSA and VLBA, which were all averaged over the top-set reconstructions similar to the GMVA images in Figure 3. The overall jet morphology is consistent with previous ground-based (MOJAVE: Lister et al. 2018; VLBA-BUBlazar 39 : Jorstad et al. 2017)and space-based VLBI observations (Radio Astron: Bruni et al. 2017,2021). Regardless of the observing frequency, the 3C 273 jet extends in the southwest direction with multiple bright knots whose relative locations are consistent between different frequencies. The HSA 43 GHz image shows the inner jet structure within ∼2 mas, and two extended components around 4 and 7 mas from the core. These overall morphology is consistent with a Figure 3. The inner-jet image of 3C 273 obtained with the GMVA observation including ALMA at 86 GHz. Here we show the mean of the top-set reconstructions restored at the resolution with a uniform weighting scheme. The image is restored with a circular beam of 57 μas (indicated in the bottom left corner)corresponding to the geometric mean of the uniform-weighting beam of 61 ×52 μas. The peak intensity is 0.69 Jy beam −1 . The contours are multiplies by factors of 2, from the lowest contour level of 6.5 mJy beam −1 . 39 http://www.bu.edu/blazars/VLBAproject.html 7 The Astrophysical Journal, 940:65 (17pp), 2022 November 20 Okino et al. close-in-time image in the database of the VLBA-BU-Blazar program observed a week previously (2017 March 19). The HSA 15 and 22 GHz images show the same regions with a bright core and extended jet emissions up to ∼25 mas from the core. One can see that the direction of the jet changes at ∼10 mas from the core, with the PA changing from ∼−140° (upstream)to ∼−120°(downstream). In the jet, several components are located at the same position at both frequencies. However, the faint jet component is not detected at 22 GHz unlike at 15 GHz, due to the optically thin synchrotron emission from the jet (see Section 5.3)and also the limited uv coverage for short baselines at 22 GHz compared with 15 GHz. In our VLBA 1.7 GHz image, the jet emissions downstream of ∼20 mas are reconstructed with our imaging. The image is similar to the 3C 273 jet shown in Kovalev et al. (2016)at the same frequency band. One can see the bright component at ∼40 mas and the extended jet emission at ∼100 mas from the core. 5.3. Spectral Index Maps We measured the spectral index (defined αas I ν ∝ν α )of the jet emission between two adjacent frequency bands. To derive the spectral index map for each pair of images, we first aligned the two images by adopting the offset determined by the twodimensional cross-correlation (see Section 4.1). Then, the spectral index maps were derived for all possible pairs of topset images at two frequencies. Then, we averaged over all spectral index maps for each frequency pair to create the mean spectral index map. In Figure 6, we show the mean spectral index map of the 3C 273 jet for each pair of neighboring frequency bands. In all frequency pairs, the 3C 273 jet consistently shows spectra that are optically thick upstream of the core, become flat near the core, and then become optically thin downstream. These features are commonly observed in radio jet sources (e.g., Hovatta et al. 2014), therefore suggesting that the core offsets derived in Section 4.1 are plausible. 5.4. Core Shift In Table 3, we show the measured core offsets for each pair of frequencies. The measured offsets are in the x-direction, parallel to the large-scale jet axis. The standard deviations of the measured offsets are the same or less than the pixel sizes of the regridded images used for these measurements. Therefore, the measured offsets do not depend on the different combinations of imaging parameter sets used to reconstruct the top-set images. We fit a power-law model of the core shift caused by the synchrotron self-absorption of the jet emission (e.g., Lobanov 1998). We assume that the core location at the frequency ν i relative to the 86 GHz core is given by () ()nnD= - -- xA ,2 ix i kk 86 () ()nnD= - -- yA ,3 iyi kk 86 for each direction, where ν 86 =86.3 GHz. The power-law index kdepends on the optically thin spectral index of the synchrotron emission, magnetic field, and particle density distributions. Assuming a conical jet with a constant velocity and in the equipartition state between magnetic and particle energy densities, k≈1 has been shown from the synchrotron self-absorption model proposed by Blandford & Königl (1979), which was confirmed in many sources by previous core shift measurements (see also Pushkarev et al. 2012, and references therein). For 3C 273, Lisakov et al. (2017)has performed multiepoch core shift measurements and reported values of k close to 1 for each epoch. Therefore, given the small sets of frequency pairs available in the presented observations, we assume k=1 and only derive coefficients (A x ,A y )in our analysis. In this model, the positional offsets (μ x ,μ y )of the radio cores between each frequency pair (ν i ,ν j )are described as () ()mnn=D -D = - -- xxA ,4 xijx i k j k () ()mnn=D -D = - -- yyA .5 yij yi k j k We derive the best-fit model parameters by comparing the model offsets above with measured ones through the leastsquares fitting based on χ 2 defined by ⎧ ⎨ ⎩⎫ ⎬ ⎭ ⎧ ⎨ ⎩⎫ ⎬ ⎭ () () () () () cnn m s nn m s =å -- +å -- -- -- A A.6 ij xi k j k x x ij yi k j k y y 2 , 2 , 2 We found the best-fit parameters of A x =−5.635 ±2.424 and A y =0.278 ±2.424, which minimize χ 2 . The errors here indicate a 68.3% confidence interval derived from the χ 2 surface. With these values, Equations (2)and (3)can be used to derive the location of the upstream end of the jet, often considered the black hole position (e.g., Hada et al. 2011),by Figure 4. Standard deviation of the top-set reconstructions of GMVA 86 GHz restored with the circular Gaussian beam of 57 μas shown by the white circle in the bottom left corner. The contours represent the mean image shown in Figure 3. 8 The Astrophysical Journal, 940:65 (17pp), 2022 November 20 Okino et al. taking ν i →∞. From the best-fit parameters, we obtained the location of the jet apex at x 0 =65 μas and y 0 =−3μas upstream 40 of the 86 GHz core for 3C 273. Figure 7shows the resultant core positions along the large-scale jet axis (i.e., xdirection). The lower-frequency core is located at a greater distance from the central black hole, as anticipated for the core shift caused by synchrotron self-absorption (Blandford & Königl 1979; Lobanov 1998). 5.5. Jet Collimation Profile of 3C 273 Figure 8shows the radius and position angle profile of the 3C 273 jet as a function of the deprojected distance from the central black hole. Our VLBI measurements cover deprojected distances ranging from ∼10 5 R s to ∼10 8 R s . The radius profile at each frequency smoothly connects. Near z∼10 7 R s , the jet radius profile shows a sharp increase by a factor of ∼2, where the PA of the jet sharply changes by ∼20°to the north. The overall PA profile meanders around the mean PA of −138°, coinciding with the large-scale jet (Conway et al. 1993). The jet radius profile in the outer region beyond ∼10 7 R s is steeper than that in the inner region. Thus, we fit the power-law functions to the jet radius profile to quantify the jet shape in each region. The details of our model selection are described in the Appendix. Here, we excluded the region from 4 ×10 6 R s to 10 7 R s (the gray shaded region in Figure 8)from our fitting analysis because of its peculiar and complex morphologies not seen in the rest of the jet. In Figure 9, we show the streamline of the jet image of HSA 15 GHz to illustrate the region that we excluded from our fitting analysis. This region shows a rapid swing in the streamline from the global position angle, and simultaneously the vertical profiles of the jet show highly asymmetric ridge-like structures not seen in the rest of the global jet. These peculiar morphologies make it difficult to define and measure the width of the jet in a consistent way with the rest of the global jet. First, we model the jet radius r(z)with a couple of single power-law functions (r(z)∝z a ), separately fitting to each of the upstream (z<4×10 6 R s )and downstream (z>10 7 R s ) regions. The best-fit values and associated uncertainties of the parameters are estimated using the percentile bootstrap method, where in each trial the parameters are derived by applying an unweighted least-squares method to the data sets randomly resampled from the original data, allowing duplication. After 10,000 bootstrap trials, we adopt the median values as the bestfit parameters, and one-third of the 99.7% percentile confidence intervals as their 1σuncertainties. The power-law index is found to be =- + a 1.28 0.05 0.08, indicating a conical/hyperbolic shape in the downstream region, whereas =- + a 0.66 0.05 0.04 is indicative of parabolic collimation found in the upstream region. Furthermore, we also examine a joint fitting of both regions with a broken power-law function defined by ⎜⎟ ⎡ ⎣ ⎢⎛ ⎝⎞ ⎠⎤ ⎦ ⎥ () () () () / / =+ - - rz R z z 21 7 aa aa 02 b 2 ud ud following Nakahara et al. (2018).R 0 is the scale factor, z b is the location of the break, and a u and a d are the power-law indices of the upstream and downstream regions, respectively. Here sis the parameter controlling the sharpness of the curvature at the transition point, which is fixed to be s=10. The best-fit parameters and their errors are derived using the bootstrapping method in the same manner. The power-law indices in the downstream and upstream regions are found to be =- + a 1.31 d0.06 0.10 and =- + a 0.67 u0.04 0.02, respectively, which are both consistent with the results of the above double power-law fitting. The jet shape transition is found to occur at =´ - + z R1.1 10 b0.2 0.1 7 s . These results show that the 3C 273 jet has a transition in its shape from semiparabolic to conical/ hyperbolic at the deprojected distance of ∼10 7 R s from the central black hole. 6. Discussion 6.1. Nature of the Jet Collimation in 3C 273 Our new quasi-simultaneous observations clearly show the transition of the jet shape from semiparabolic to conical/ hyperbolic at ∼10 7 R s from the central black hole, providing the first compelling example of a collimation break in a quasar jet. Interestingly, the transitional location roughly coincides with the region at ∼4×10 6 to 10 7 R s where the jet width (2r) Figure 5. Multifrequency images of 3C 273 jets with the HSA at 43, 22, and 15 GHz and the VLBA at 1.7 GHz. Each panel shows the mean top-set image at each frequency restored by the circular Gaussian beam with the beam solid angle of the uniform-weighting beam size in Table 1. The lowest contour levels are 6.5 mJy beam −1 at 43 GHz, 8.2 mJy beam −1 at 22 GHz, 7.1 mJy beam −1 at 15 GHz, and 11.6 mJy beam −1 at 1.7 GHz. The contours for all images are multiplied by a factor of 2. The lowest contour level of each image is estimated from the mean residual image rms of the top-set reconstructions. 40 The location of the jet apex was derived for k=1, and may vary by at most a factor of two in the range k∼0.8–1.2 constrained by Lisakov et al. (2017). These uncertainties do not affect the main conclusions of this paper. 9 The Astrophysical Journal, 940:65 (17pp), 2022 November 20 Okino et al. ORCID iDs Hiroki Okino https://orcid.org/0000-0003-3779-2016 Kazunori Akiyama https://orcid.org/0000-0002-9475-4254 Keiichi Asada https://orcid.org/0000-0001-6988-8763 José L. Gómez https://orcid.org/0000-0003-4190-7613 Kazuhiro Hada https://orcid.org/0000-0001-6906-772X Mareki Honma https://orcid.org/0000-0003-4058-9000 Thomas P. Krichbaum https://orcid.org/0000-00024892-9586 Motoki Kino https://orcid.org/0000-0002-2709-7338 Hiroshi Nagai https://orcid.org/0000-0003-0292-3645 Uwe Bach https://orcid.org/0000-0002-7722-8412 Lindy Blackburn https://orcid.org/0000-0002-9030-642X Katherine L. Bouman https://orcid.org/0000-0003-0077-4367 Andrew Chael https://orcid.org/0000-0003-2966-6220 Geoffrey B. Crew https://orcid.org/0000-0002-2079-3189 Sheperd S. Doeleman https://orcid.org/0000-00029031-0904 Vincent L. Fish https://orcid.org/0000-0002-7128-9345 Ciriaco Goddi https://orcid.org/0000-0002-2542-7743 Sara Issaoun https://orcid.org/0000-0002-5297-921X Michael D. Johnson https://orcid.org/0000-0002-4120-3029 Svetlana Jorstad https://orcid.org/0000-0001-6158-1708 Shoko Koyama https://orcid.org/0000-0002-3723-3372 Colin J. Lonsdale https://orcid.org/0000-0003-4062-4654 Ru-Sen Lu https://orcid.org/0000-0002-7692-7967 Ivan Martí-Vidal https://orcid.org/0000-0003-3708-9611 Lynn D. Matthews https://orcid.org/0000-0002-3728-8082 Yosuke Mizuno https://orcid.org/0000-0002-8131-6730 Kotaro Moriyama https://orcid.org/0000-0003-1364-3761 Masanori Nakamura https://orcid.org/0000-00016081-2420 Hung-Yi Pu https://orcid.org/0000-0001-9270-8812 Eduardo Ros https://orcid.org/0000-0001-9503-4892 Tuomas Savolainen https://orcid.org/0000-0001-6214-1085 Fumie Tazaki https://orcid.org/0000-0003-0236-0600 Figure 13. Upper panel: radial profile of the jet in 3C 273 with different models of power-law functions. Measurement data are the same as in Figure 8. The gray dashed line shows the best-fit single power-law function (r∝z a )with a=0.80 obtained from all the data except for the gray shaded region. The red solid line shows the best-fit broken power-law function as presented in Figure 8. Lower panel: fractional deviations of the jet radii from the obtained models. The red and gray squares correspond to the deviations from the lines in the upper panel. 16 The Astrophysical Journal, 940:65 (17pp), 2022 November 20 Okino et al. Jan Wagner https://orcid.org/0000-0003-1105-6109 Maciek Wielgus https://orcid.org/0000-0002-8635-4242 Anton Zensus https://orcid.org/0000-0001-7470-3321 References Agudo, I., Thum, C., Molina, S. N., et al. 2018, MNRAS,474, 1427 Akiyama, K., Asada, K., Fish, V., et al. 2018, Galax,6, 15 Akiyama, K., Ikeda, S., Pleau, M., et al. 2017a, AJ,153, 159 Akiyama, K., Kuramochi, K., Ikeda, S., et al. 2017b, ApJ,838, 1 Algaba, J. C., Nakamura, M., Asada, K., & Lee, S. S. 2017, ApJ,834, 65 Algaba, J. C., Rani, B., Lee, S. S., et al. 2019, ApJ,886, 85 Angelakis, E., Fuhrmann, L., Marchili, N., et al. 2015, A&A,575, A55 Asada, K., & Nakamura, M. 2012, ApJL,745, L28 Asada, K., Nakamura, M., Doi, A., Nagai, H., & Inoue, M. 2014, ApJL,781, L2 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 Baczko, A. K., Ros, E., Kadler, M., et al. 2022, A&A,658, A119 Bahcall, J. N., Kirhakos, S., Schneider, D. P., et al. 1995, ApJL,452, L91 Begelman, M. C., & Li, Z.-Y. 1994, ApJ,426, 269 Blackburn, L., Pesce, D. W., Johnson, M. D., et al. 2020, ApJ,894, 31 Blandford, R., Meier, D., & Readhead, A. 2019, ARA&A,57, 467 Blandford, R. D., & Königl, A. 1979, ApJ,232, 34 Boccardi, B., Krichbaum, T. P., Bach, U., Bremer, M., & Zensus, J. A. 2016, A&A,588, L9 Boccardi, B., Krichbaum, T. P., Ros, E., & Zensus, J. A. 2017, A&ARv,25, 4 Boccardi, B., Migliori, G., Grandi, P., et al. 2019, A&A,627, A89 Boccardi, B., Perucho, M., Casadio, C., et al. 2021, A&A,647, A67 Bromberg, O., Nakar, E., Piran, T., & Sari, R. 2011, ApJ,740, 100 Bruni, G., Gómez, J. L., Casadio, C., et al. 2017, A&A,604, A111 Bruni, G., Gómez, J. L., Vega-García, L., et al. 2021, A&A,654, A27 Cappi, M., Matsuoka, M., Otani, C., & Leighly, K. M. 1998, PASJ,50, 213 Chael, A., Rowan, M., Narayan, R., Johnson, M., & Sironi, L. 2018, MNRAS, 478, 5209 Chael, A. A., Johnson, M. D., Narayan, R., et al. 2016, ApJ,829, 11 Conway, R. G., Garrington, S. T., Perley, R. A., & Biretta, J. A. 1993, A&A, 267, 347 Croke, S. M., & Gabuzda, D. C. 2008, MNRAS,386, 619 Davis, R. J., Muxlow, T. W. B., & Conway, R. G. 1985, Natur,318, 343 Deller, A. T., Brisken, W. F., Phillips, C. J., et al. 2011, PASP,123, 275 Deller, A. T., Tingay, S. J., Bailes, M., & West, C. 2007, PASP,119, 318 Doi, A., Hada, K., Kino, M., Wajima, K., & Nakahara, S. 2018, ApJL,857, L6 Event Horizon Telescope Collaboration, Akiyama, K., Alberdi, A., et al. 2019a, ApJL,875, L1 Event Horizon Telescope Collaboration, Akiyama, K., Alberdi, A., et al. 2019b, ApJL,875, L4 Fromm, C. M., Perucho, M., Porth, O., et al. 2018, A&A,609, A80 Fromm, C. M., Ros, E., Perucho, M., et al. 2013, A&A,557, A105 Fromm, C. M., Younsi, Z., Baczko, A., et al. 2019, A&A,629, A4 Fuhrmann, L., Angelakis, E., Zensus, J. A., et al. 2016, A&A,596, A45 Fujita, Y., Kawakatu, N., & Shlosman, I. 2016, PASJ,68, 26 Giovannini, G., Savolainen, T., Orienti, M., et al. 2018, NatAs,2, 472 Goddi, C., Martí-Vidal, I., Messias, H., et al. 2019, PASP,131, 075003 Gómez, J. L., Lobanov, A. P., Bruni, G., et al. 2016, ApJ,817, 96 Gravity Collaboration, Sturm, E., Dexter, J., et al. 2018, Natur,563, 657 Gravity Collaboration, Dexter, J., Shangguan, J., et al. 2020, A&A,635, A92 Gravity Collaboration, Abuter, R., Accardo, M., et al. 2017, A&A,602, A94 Greisen, E. W. 2003, in AIPS, the VLA, and the VLBA, Astrophysics and Space Science Library, ed. A. Heck, Vol. 285 (Dordrecht: Kluwer),109 Hada, K. 2019, Galax,8, 1 Hada, K., Doi, A., Kino, M., et al. 2011, Natur,477, 185 Hada, K., Doi, A., Wajima, K., et al. 2018, ApJ,860, 141 Hada, K., Kino, M., Doi, A., et al. 2013, ApJ,775, 70 Hada, K., Kino, M., Doi, A., et al. 2016, ApJ,817, 131 Haga, T., Doi, A., Murata, Y., et al. 2015, ApJ,807, 15 Högbom, J. A. 1974, A&AS, 15, 417 Homan, D. C., Lister, M. L., Kovalev, Y. Y., et al. 2015, ApJ,798, 134 Honma, M., Akiyama, K., Uemura, M., & Ikeda, S. 2014, PASJ,66, 95 Hovatta, T., Aller, M. F., Aller, H. D., et al. 2014, AJ,147, 143 Hunter, J. D. 2007, CSE,9, 90 Husemann, B., Bennert, V. N., Jahnke, K., et al. 2019, ApJ,879, 75 Ikeda, S., Tazaki, F., Akiyama, K., Hada, K., & Honma, M. 2016, PASJ,68, 45 Issaoun, S., Johnson, M. D., Blackburn, L., et al. 2019, ApJ,871, 30 Janssen, M., Falcke, H., Kadler, M., et al. 2021, NatAs,5, 1017 Jester, S., Röser, H.-J., Meisenheimer, K., & Perley, R. 2005, A&A,431, 477 Jorstad, S. G., Marscher, A. P., Lister, M. L., et al. 2005, AJ,130, 1418 Jorstad, S. G., Marscher, A. P., Morozova, D. A., et al. 2017, ApJ,846, 98 Kameno, S., Sawada-Satoh, S., Inoue, M., Shen, Z.-Q., & Wajima, K. 2001, PASJ,53, 169 Kim, J.-Y., Krichbaum, T. P. & Event Horizon Telescope Collaboration 2020, A&A,640, A69 Kim, J.-Y., Krichbaum, T. P., Lu, R.-S., et al. 2018, A&A,616, A188 Komissarov, S. S., Barkov, M. V., Vlahakis, N., & Königl, A. 2007, MNRAS, 380, 51 Komissarov, S. S., Vlahakis, N., Königl, A., & Barkov, M. V. 2009, MNRAS, 394, 1182 Kormendy, J., & Ho, L. C. 2013, ARA&A,51, 511 Kovalev, Y. Y., Kardashev, N. S., Kellermann, K. I., et al. 2016, ApJL, 820, L9 Kovalev, Y. Y., Pushkarev, A. B., Nokhrina, E. E., et al. 2020, MNRAS, 495, 3576 Koyama, S., Kino, M., Giroletti, M., et al. 2016, A&A,586, A113 Krichbaum, T. P., Witzel, A., Booth, R. S., et al. 1990, A&A, 237, 3 Kuramochi, K., Akiyama, K., Ikeda, S., et al. 2018, ApJ,858, 56 Li, Y.-R., Wang, J.-M., Songsheng, Y.-Y., et al. 2022, ApJ,927, 58 Lisakov, M. M., Kovalev, Y. Y., Savolainen, T., Hovatta, T., & Kutkin, A. M. 2017, MNRAS,468, 4478 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 Lister, M. L., Homan, D. C., Kellermann, K. I., et al. 2021, ApJ,923, 30 Lobanov, A. P. 1998, A&A, 330, 79 Lobanov, A. P., & Zensus, J. A. 2001, Sci,294, 128 Matthews, L. D., Crew, G. B., Doeleman, S. S., et al. 2018, PASP,130, 015002 McKinney, W. 2010, in Proc. of the 9th Python in Science Conf., ed. S. van der Walt & J. Millman (Austin: SciPy),56 Mertens, F., & Lobanov, A. 2015, A&A,574, A67 Meyer, E. T., Sparks, W. B., Georganopoulos, M., et al. 2016, ApJ,818, 195 Mizuno, Y., Gómez, J. L., Nishikawa, K.-I., et al. 2015, ApJ,809, 38 Nagai, H., Haga, T., Giovannini, G., et al. 2014, ApJ,785, 53 Nakahara, S., Doi, A., Murata, Y., et al. 2018, ApJ,854, 148 Nakahara, S., Doi, A., Murata, Y., et al. 2019, ApJ,878, 61 Nakahara, S., Doi, A., Murata, Y., et al. 2020, AJ,159, 14 Nakamura, M., & Asada, K. 2013, ApJ,775, 118 Nakamura, M., Asada, K., Hada, K., et al. 2018, ApJ,868, 146 Nokhrina, E. E., Beskin, V. S., Kovalev, Y. Y., & Zheltoukhov, A. A. 2015, MNRAS,447, 2726 Nokhrina, E. E., Kovalev, Y. Y., & Pushkarev, A. B. 2020, MNRAS, 498, 2532 O’Sullivan, S. P., & Gabuzda, D. C. 2009, MNRAS,400, 26 Park, J., Hada, K., Kino, M., et al. 2019a, ApJ,871, 257 Park, J., Hada, K., Kino, M., et al. 2019b, ApJ,887, 147 Park, J., Hada, K., Nakamura, M., et al. 2021, ApJ,909, 76 Perley, R. A., & Meisenheimer, K. 2017, A&A,601, A35 Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2014, A&A,571, A16 Pushkarev, A. B., Hovatta, T., Kovalev, Y. Y., et al. 2012, A&A,545, A113 Pushkarev, A. B., Kovalev, Y. Y., Lister, M. L., & Savolainen, T. 2017, MNRAS,468, 4992 Russell, H. R., Fabian, A. C., McNamara, B. R., et al. 2018, MNRAS, 477, 3583 Russell, H. R., Fabian, A. C., McNamara, B. R., & Broderick, A. E. 2015, MNRAS,451, 588 Savolainen, T., Wiik, K., Valtaoja, E., & Tornikoski, M. 2006, A&A,446, 71 Schmidt, M. 1963, Natur,197, 1040 Shepherd, M. C. 1997, in ASP Conf. Ser. 125, Astronomical Data Analysis Software and Systems VI, ed. G. Hunt & H. Payne (San Francisco, CA: ASP),77 Sun, M. 2009, ApJ,704, 1586 Thompson, A. R., Moran, J. M., & Swenson, G. W., Jr. 2017, Interferometry and Synthesis in Radio Astronomy (3rd ed.; Berlin: Springer) Tseng, C.-Y., Asada, K., Nakamura, M., et al. 2016, ApJ,833, 288 van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, CSE,13, 22 Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe,17, 261 Walker, R. C., Hardee, P. E., Davies, F. B., Ly, C., & Junor, W. 2018, ApJ, 855, 128 Yamaguchi, M., Akiyama, K., Tsukagoshi, T., et al. 2020, ApJ,895, 84 Zhao, G.-Y., Gómez, J. L., Fuentes, A., et al. 2022, ApJ,932, 72 17 The Astrophysical Journal, 940:65 (17pp), 2022 November 20 Okino et al.