Collimation of the Relativistic Jet in the Quasar 3C 273

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 ∼105 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 107 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.

1. INTRODUCTION Relativistic jets ejected from active galactic nuclei (AGN) 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 AGN jet physics (e.g., reviewed in Blandford et al. 2019).Particularly, 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. 2007Komissarov et al. , 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 (MHD) 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 the M87 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;Nakamura & Asada 2013;Hada et al. 2013Hada et al. , 2016;;Kim et al. 2018;Walker et al. 2018;Park et al. 2019a;Event Horizon Telescope Collaboration 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), 3C 84 (Nagai et al. 2014;Giovannini et al. 2018), and for a large number of jet samples collected by the long-term monitoring MOJAVE program 1 (Pushkarev et al. 2017).More recently, Kovalev et al. (2020) investigated jet widths from over 300 sources and found the 1 https://www.physics.purdue.edu/MOJAVE/index.htmltransitions of jet shapes from ten 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 AGN 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 wellstudied 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 (pc) scale with VLBI observations (Krichbaum et al. 1990;Lobanov & Zensus 2001;Savolainen et al. 2006;Kovalev et al. 2016;Bruni et al. 2017Bruni et al. , 2021;;Jorstad et al. 2017;Lister et al. 2019Lister et al. , 2021)).Furthermore, recent infrared interferometric observations with the GRAVITY instrument on the Very Large Telescope Interferometer (VLTI) (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 and θ = 12 • ± 2 • (Gravity Collaboration et al. 2018) 2 .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 frequencydependent locations of radio cores at different frequencies due to the effect of synchrotron self-absorption relevant for the inner-most 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 phased Atacama Large Millimeter/submillimeter Array (ALMA) (Matthews et al. 2018) provided the 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 sub-pc 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. In this paper, we assume a flat ΛCDM cosmology (e.g., Planck Collaboration 2014) with H 0 = 70 km s −1 Mpc −1 , Ω M = 0.3 and Ω Λ = 0.7, and we adopt 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 sub-band IF had an instrumental offset of 180 degrees in subsequent channels of 32 MHz widths owing to a misconfiguration of the sub-band 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 (BP) table.Visibility amplitudes were a-priori calibrated in the standard manner; visibilities at each baseline were first normalized with the available auto-correlation spectra and then scaled with the system equivalent flux densities (SEFDs) 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 the 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.
3 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   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 1 shows 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, the ALMA provides very high SNR detection 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 the observation.

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 March 26, 2017, eight days before our GMVA observations (project code: BA122).The observing array consisted of the ten VLBA stations and the 100-m Effelsberg (EB) telescope in Germany, as summarized in Table 1.Full-track observations were performed for ∼11 h, of which EB participated in the first ∼3.4h.Data were recorded at dual circular polarizations with four 64 MHz IFs subdivided into 256 channels, providing 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 apriori 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.

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 February 2008 (project code: BH151).The data set has one of the best hour-angle coverage among existing archival data sets at 1.7 GHz 4 .Data were reduced in a standard manner in AIPS.
Fringes were detected on all ten VLBA stations, as summarized in Table 1.

IMAGING
We reconstructed multifrequency images of 3C 273 from the GMVA, HSA, and VLBA data sets presented in Section 2 with regularized maximum likelihood (RML) methods implemented in SMILI (Akiyama et al. 2017a,b).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 uvcoverages lacking intermediate baselines as shown in Figure 1.The traditional iterative hybrid mapping using CLEAN (e.g., Högbom 1974) and 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 EHT imaging of M 87 (Event Horizon Telescope Collab-4 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 jet width of the 1.6 GHz jet does not change significantly over the time scale of 10 years. oration 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 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.

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. 2016Chael et al. , 2018;;Akiyama et al. 2017a,b;Kuramochi et al. 2018) but also for actual interferometric measurements with VLBI arrays (e.g., Issaoun et al. 2019;Event Horizon Telescope Collaboration 2019b;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 mathematically describing 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. 2016Chael et al. , 2018;;Akiyama et al. 2017a;Blackburn et al. 2020), and further inclusion of various observing effects such as systematic non-closing errors (Event Horizon Telescope Collaboration 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,b;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,b;Event Horizon Telescope Collaboration 2019b).Further technical and mathematical details of the RML approaches adopted in this study are described in Event Horizon Telescope Collaboration 2019b.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 , 1  1  0, 10 −2 , 10 −1 , 1  0, 10 −2 , 10 −1 , 1  64  24 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 weighted-1 , TV, TSV and MEM regularizers (see Appendix A of Event Horizon Telescope Collaboration 2019b, for mathematical definitions).The pixel size of the image is typically set to one-tenth of the mean full width at half maximum (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 to use a narrower rectangular image field of view and minimizes the computational cost of imaging.Then, the visibilities were coherently averaged for 60 sec, after manual 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).

Imaging Pipeline
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.

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 hundreds, of parameter sets with four to five parameters: potential systematic non-closing errors (denoted as Sys err) and the regularization parameters for the w 1 , TV, TSV, and MEM regularizers.Among the sets of the 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 I on May 25, 2017 (project code BL229AH) from the MOJAVE program (Lister et al. 2018), which was the closest epoch to our HSA observations, as the softmask of the imaging region with w 1 regularization.This allows noise suppression outside of the area, which has historically had no significant emission5 .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 L-band 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, fifteen to forty-five images were selected at each frequency, as summarized in Table 2, and used in post-imaging analysis.

Total flux scaling for GMVA 86 GHz
The overall scaling of visibility amplitudes strongly depends on the accuracy of the a-priori calibration of visibility amplitudes based on measurements of the systemequivalent 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) is ∼ 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).

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.

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 selfabsorption (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)6 .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 (1) where I ν,ij is the intensity at pixel (i, j) at frequency ν, and Ī 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.

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 taking a slice perpendicular to the PA at each distance, and the other is slicing 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 width of the jet was obtained by FWHM = 2 √ 2 ln 2σ w .We derived the deconvolved FWHM width (w) defined as w = FWHM 2 − θ 2 beam , correcting the blurring effects with the restoring beam at the FWHM of θ beam7 (see Table 1).Finally, the jet radius r at each distance was obtained from half size of the deconvolved FWHM width (i.e., r = w/2).
The jet radius profile was measured with an interval of half of the restoring beam size (see Table 1).We only accept the jet radius measurements satisfying 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.

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 weighting8 .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 NS direction by a factor greater than 2.5.The spatial resolution of our observations is also higher than previous 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 mas 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 as the representative image uncertainty at 86 GHz with natural weighting scheme.We then determined the lowest contour level for the uniform weighing 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 4 shows the standard deviation image of GMVA 86 GHz across the top-set reconstructions.The brightest core shows the largest (∼ 200 mJy/beam), 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 the ALMA in Band 3 during this GMVA campaign (Goddi et al. 2019).The extended jet emission at 3mm is broadly consistent with the 7mm jet implying optically thin spectra (see Section 5.3), 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.

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

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 two-dimensional crosscorrelation (see Section 4.1).Then, the spectral index maps were derived for all possible pairs of top-set 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 optically thick spectra upstream of the core, turns into flat spectra near the core, and then becomes optically thin downstream.These features are commonly observed in radio jet sources (e.g., Hovatta et al. 2014), therefore supporting that the core offsets derived in Section 4.1 are plausible.

Core shift
In Table 3, we show the measured core offsets for each combination of frequency pairs.The measured offsets are along 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 re-gridded 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 ∆x for each direction, where ν 86 = 86.3GHz.The powerlaw index k depends on the optically thin spectral index  In this model, the positional offsets (µ x , µ y ) of the radio cores between each frequency pair (ν i , ν j ) are described as We derive the best-fit model parameters by comparing the model offsets above with measured ones through the least-squares fitting based on χ 2 defined by 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 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 upstream10 of the 86 GHz core for 3C 273. Figure 7 shows the resultant core positions along the large-scale jet axis (i.e.x-direction).The lower-frequency core is located at a farther distance from the central black hole, as anticipated for the core shift caused by synchrotron selfabsorption (Blandford & Königl 1979;Lobanov 1998).

Jet collimation profile of 3C 273
Figure 8 shows 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 the de-projected 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 • 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 at the outer region greater than ∼ 10 7 R s is steeper than that of 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 Appendix A. 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 in 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 (> 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 unweighted leastsquares method to the data sets randomly re-sampled from the original data allowing duplication.After 10,000 times bootstrap trials, we adopt the median values as the best-fit 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.08 −0.05 , in-dicating a conical/hyperbolic shape in the downstream region, whereas a = 0.66 +0.04 −0.05 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 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 s is 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 d = 1.31 +0.10 −0.06 and a u = 0.67 +0.02 −0.04 , respectively, which are both consistent with the results of the above double power-law fitting.The jet shape transition is found to be z b = 1.1 +0.1 −0.2 × 10 7 R 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.Our new quasi-simultaneous observations clearly show the transition of the jet shape from semi-parabolic 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 − 10 7 R s where the jet width (2r) rapidly expands by a factor of ∼ 2 along the significant bending, which was not used in the fitting because of its peculiarity (see Figure 9).Recent space VLBI observations with Radio Astron at 1.6 GHz revealed the complicated limbbrightened jet emission in this area (Bruni et al. 2021).Mertens & Lobanov (2015) shows non-ballistic motions of jet components along the northern and southern limbs over 14 years, which is longer than the jet crossing time of this region considering the typical proper motion of the jet (∼ 1 mas/year, Lister et al. 2019).After passing through the location of the collimation break, the northern and southern limbs smoothly merge into a single limb, and its position angle farther downstream is timely and spatially stable over 26 years of 15 GHz observations of the MOJAVE program11 .The stability of the broadened, double-ridge structure, apparently around the collimation break, suggests a stationary magnetohydrodynamic feature, possibly triggered by some changes in the  , which is broadly consistent with the axis of the large-scale jet (e.g., Davis et al. 1985;Conway et al. 1993).The gray shade across both panels shows the region excluded from the fitting analysis because of a large bending seen in the jet (see the details in Section 5.5).circum-jet environment causing pressure mismatch with the jet (e.g., Mizuno et al. 2015), as the HST-1 knot in M87 near the collimation break is indicated to be (Asada & Nakamura 2012).
The jet velocity field is another important physical quantity that is closely related to jet collimation.This prediction is supported by various VLBI observations of M87 (Nakamura & Asada 2013;Asada et al. 2014;Park et al. 2019bPark et al. , 2021)), 1H 0323+342 (Hada et al. 2018), and statistical studies of large samples of AGN jets (Homan et al. 2015).In Figure 10, we show the profile of the apparent jet velocities of 3C 273 from long-term monitoring observations at 15 GHz by the MOJAVE (Lister et al. 2019) and 43 GHz by the BU-Blazar (Jorstad et al. 2017) programs, compared to the collimation profile.Except for the stationary knot features at 1 mas identified in the literature (Savolainen et al. 2006;Jorstad et al. 2005Jorstad et al. , 2017;;Lisakov et al. 2017;Bruni et al. 2017), Figure 10 indicates that the jet has already been accelerated to relativistic speed with an apparent velocity of at least ∼ 5c before the collimation break.This trend in the 3C 273 jet differs from that of M87 showing the gradual accelerations in the parabolic regions and velocity saturation at the collimation break (Asada et al. 2014;Park et al. 2019a).Further observations of the jet kinematics in both the inner (< 1 mas) and outer regions (> 10 mas) would be useful for investigating the detailed velocity field in the collimation zone of 3C 273, providing useful comparisons to radio galaxies.
For the inner-most jet at the inner sub-pc scale, one can compare our results with recent sub-milliarcsecondscale observations with near-infrared optical interferometry.3C 273 has been intensively studied with the Very Large Telescope Interferometer (VLTI), which recently resolved the broad line regions (BLRs) with a size of 46 µas ∼ 0.1 pc as an ionized gas disk (Gravity Collaboration et al. 2018), and the dust torus emission with a size of ∼ 150 µas ∼ 0.4 pc (Gravity Collaboration et al. 2020).The innermost collimation profile of the 3C 273 jets, expected by extrapolating to a farther inner region, is smaller than the measured geometry of the resolved BLRs and hot dust torus distributions.Such a region, where the jet is surrounded by the BLRs and dust torus, would be reachable with future higher frequency observations, such as with the EHT at 230 or 345 GHz.Future higher-frequency VLBI observations will be critical for understanding the role of such an external medium common in quasars in jet collimation.
Downstream of the collimation break, the jet structure up to the ∼kpc scale shows a conical/hyperbolic expansion with a power-law index of a = 1.28.These unconfined structures downstream of the jets have been found in other AGN jets (Pushkarev et al. 2017;Kovalev et al. 2020).To investigate the jet shape beyond the VLBI scales (>∼ 1 kpc), we show additional measurements based on CLEAN images from MERLIN and VLA observations at 1.6 GHz and 5 GHz in Akiyama et al. (2018) and Perley & Meisenheimer (2017), respectively, in Figure 11.These additional measurements are obtained from the kpc-scale radio lobe at the terminal end of the 3C 273 jet.Although previous observations detected the intermediate-scale (∼ 1 − 10 arcsec) jet emission between the central core and the kpc-scale lobe (e.g., Perley & Meisenheimer 2017), the transverse structure is not well resolved, and the radius cannot be measured reliably.This causes a lack of measurements in 10 8 − 10 10 R s .The inner semi-parabolic line, if extrapolated, significantly underestimates the jet radius at the kpc-scale.On the other hand, the jet radius at the kpc-scale is smaller than the extrapolated conical/hyperbolic line, and these kpc-scale jets can be more conically connected to 1.7 GHz measurements.Future observations of the intermediate scale would be useful for constraining the propagation of the jet into the kpc-scale radio lobe, and investigating detailed jet features interacting with the ambient ionized/molecular gas (Husemann et al. 2019) and hot gas cocoons (Bromberg et al. 2011).

Comparison with Other AGN sources
The jet structure has been investigated for many objects in the last decades, and the transitions of jet collimation profiles have been discovered in several jets of radio galaxies and BL Lac objects.The presented re-sults on 3C 273 have added a case as a quasar jet.To compare our 3C 273 results with other AGN jets, in Figure 12, we show the relation between the black hole mass and the deprojected distance of the jet collimation break for both 3C 273 and other AGN sources.The break locations are widely distributed from ∼ 10 4 R s to ∼ 10 8 R s , where 3C 273 is located on the further side.
A popular interpretation for the origin of the jet collimation break is that the transition of the external pressure profile triggers the end of the jet collimation, where the external pressure does not provide enough support to confine the jet (e.g., Asada & Nakamura 2012;Tseng et al. 2016).For the pioneering case of M87, Asada & Nakamura (2012) discussed that the transition of the jet shape may result from the transitional change in the external pressure profile of the circum-jet medium near the Bondi radius and the end of the sphere of gravitational influence (SGI).This interpretation is motivated by X-ray observations, which indicate that the different pressure profiles of the ambient gas are realized inside and outside the Bondi/SGI radius in M87 (Russell et al. 2015(Russell et al. , 2018)).
We further investigated the hot gas properties for other sources to check whether similar conditions to M87 could be established.However, the number of sources is limited whose proximity to the central core can be spatially resolved by X-ray observations, making it difficult to give accurate estimations of the Bondi radii for those sources.Therefore, we consider the SGI radius instead of the Bondi radius for the gravitational spheres of the central black holes.Note that we checked the temperature of the bright cores of several sources hosted by elliptical galaxies in our collected samples and confirmed that the Bondi radii are approximately of the same order of magnitude as the SGI radii12 (e.g., Sun 2009;Fujita et al. 2016).
Here, the SGI radius is estimated as , where σ d is the stellar velocity dispersion.To investigate the relationship between the location of the SGI boundary and the jet collimation break among our samples with various BH masses, we introduced an emprical M BH − σ d relation derived in Kormendy & Ho (2013).This relation is particularly strong for the classical bulge and elliptical galaxies, and our sample is supposed to follow this relation (e.g., see Table 2 in Kormendy & Ho 2013, and references therein).Indeed, 3C 273 follows the relation from the stellar velocity dispersion of σ d = 210 km/s (Husemann et al. 2019) and SMBH mass of M BH = 2.6 × 10 8 M , which gives us an estimation of the SGI radius as R SGI ∼ 10 6 R s (see the solid line in Figure 12).The red star shows our results from the 3C 273 jet, and the cyan circles show results from other sources reported in the literature.The red horizontal bar indicates the region of a persistent feature seen near the break distance with a complex bending double-ridge-like morphology potentially associated with the transition in the jet shape (see Figure 9 and Section 6.1).To compare with the location of the sphere of the gravitational influence (SGI), we show the SGI locations of the classical bulges and elliptical galaxies in gray squares taken from Kormendy & Ho (2013).The black solid line and gray shaded area are the best-fit model and the 1-σ error range of the relation between the black hole mass and SGI location, transformed from the MBH − σ d relation in Kormendy & Ho (2013).The dotted lines show the actual distances in parsecs.The results of other sources are obtained from the following references; M87: Asada & Nakamura 2012, NGC 6251: Tseng et al. 2016, NGC 4261: Nakahara et al. 2018, Cygnus A: Nakahara et al. 2019, NGC 1052: Nakahara et al. 2020, NGC 315: Park et al. 2021, 4C38.41: Algaba et al. 2019, 3C 264: Boccardi et al. 2019, 1H 0323+342: Hada et al. 2018, and other sources: Kovalev et al. 2020.Note that we show two cases of different black hole masses for 1H 0323+342 (see Section 6.2).
The most important finding here is that the measured locations of the collimation breaks, including 3C 273, are more widely distributed than the intrinsic scatter of the SGI locations seen in the samples of Kormendy & Ho (2013).In particular, NGC 4261 (Nakahara et al. 2018), NGC 1052 (Nakahara et al. 2020) andNGC 315 (Boccardi et al. 2021;Park et al. 2021) imply a transition of the jet collimation in a substantially inner area than the SGI boundary.For NGC 315, Park et al. (2021) discussed that the non-relativistic outflow, such as wind from the accretion disk, may play a critical role in jet collimation -for instance, the inner break of the jet collimation seen in NGC 315 may occur if the wind does not reach the Bondi radius or the end of the SGI.
The other two sources have a common property at the jet base.They are surrounded by a dense obscuring disk/torus causing free-free absorption (FFA).For NGC 1052, the physical size of the FFA torus is ∼ 0.5 pc (∼ 3.4 × 10 4 R s ) in radius (Kameno et al. 2001), which is close to the de-projected break location of ∼ 0.15 pc reported from Nakahara et al. (2020) 13 For NGC 4261, Haga et al. (2015) estimated the transition radius in the accretion disk of ∼ 2 × 10 3 R s from FFA measurements, which is close to the jet shape break of ∼ 10 4 R s for NGC 4261 (Nakahara et al. 2018).
The break location of 1H 0323+342 is not well constrained because of the significantly uncertain black hole mass.It shows a distant collimation break for the lower black hole mass (case B in Figure 12), while the jet collimation break occurs near the SGI for the larger black hole mass (case A).Hada et al. (2018) discussed that radiation driven outflow/wind associated with the narrow-line region (NLR) may confine the jet of 1H 0323+342 (case B), while 1H 0323+342 and M87 may share the common jet collimation mechanism (case A), which is supported by 1H 0323+342 having a stationary recollimation shock at the transitional region like the HST-1 knot in M87 (see Doi et al. 2018).
Our results and Figure 12 indicate that the transition distance of the jet shape is not necessarily determined by the SGI boundary, but rather by its diverse environment, such as the presence of the disk, torus, or disk wind and their spatial extent.The hot gas cocoon surrounding the jet may also be an influence (e.g., Bromberg et al. 2011).
Recently, Boccardi et al. (2021) discussed the jet collimation properties of various radio galaxies with respect to the activity of the central accretion disk by categorizing samples into low-and high-excitation radio galaxies (LERG and HERG).They found that HERGs had larger radii and longer shape transition distances for the collimation properties.For 3C 273, it is unlikely that a component near the central engine, such as the torus, is the critical factor causing the break because the collimation break in the 3C 273 jet is outside the range of ∼ 100 pc.Then 3C 273, if viewed from a typical viewing angle of radio galaxies, would be categorized as a HERG given the high ratio of the X-ray luminosity L X (2 − 10 keV) to the Eddington luminosity L Edd implied by L X > 1.4 × 10 46 erg/s (Cappi et al. 1998) and L Edd = 3.4 × 10 46 erg/s given for M BH = 2.6 × 10 8 M (Gravity Collaboration et al. 2018).Therefore, our results of 3C 273, which show relatively large jet widths (2r >∼ 10 4 R s in z >∼ 10 5 R s ) and distant location of the jet collimation break (z b > 10 6 R s ), are consistent with the statistical trend reported in Boccardi et al. (2021).
Finally, we briefly note a scenario in which the jet collimation break is governed by the internal jet physics proposed by Kovalev et al. (2020).They proposed that under a single ambient pressure profile, such as p ∝ z −2 , the collimation break occurs in the region where the magnetic and particle energies are equivalent.In their model, the break location mainly depends on the initial magnetization parameters (σ M ) at the jet base, as well as the black hole mass and spin (see Equation 7in Nokhrina et al. 2020).The large σ M for 3C 273 estimated in Nokhrina et al. (2015) may be qualitatively consistent with the distant jet collimation break found in this study.However, it is still hard to tightly constrain the value of σ M , and it is still an open problem.Therefore, it remains unclear whether this scenario can explain the diverse distributions of the break locations ranging from 10 4 R s to 10 8 R s shown in Figure 12.If the value of σ M can be estimated accurately for more sources, the physical connection between σ M the break location could be discussed in more detail.

SUMMARY
In this paper, we have investigated the global jet structure of the archetypical quasar 3C 273 with VLBA, HSA, and GMVA observations.In particular, we have reported on new GMVA observations at 86 GHz conducted in the first session involving ALMA which significantly enhances the sensitivity and angular resolution of the images, providing the detailed morphology of the innermost jet.With the robust imaging analysis with the state-of-the-art RML imaging techniques, we obtained the jet collimation profile over a wide range of 10 5 − 10 8 R s in the deprojected distance.We summarize our results as follows.
• The quasar 3C 273 jet is found to have the structural transition from semi-parabolic to conical/hyperbolic shape at ∼ 10 7 R s , providing the first clear example of a quasar jet.Our results suggest the existence of the jet collimation break for sources exhibiting widely different accretion rates.
• The collimation break is located in the area at ∼ 4 × 10 6 − 10 7 R s where the jet shows a bending streamline for ∼ 10 mas persistent more than the typical jet crossing time.The area is also known as the broadened vertical structure with brightening two limbs resolved by recent Space VLBI observations.This peculiar and persistent feature, colocated with the jet break, may be explained by a stable magneto-hydrodynamic feature, as HST-1 may be for M87.
• The extrapolation of the collimation profile is consistent with the spatial distribution of BLRs and dusty winds resolved with recent GRAVITY observations.Our observations suggest that future higher-frequency observations, for instance, with the Event Horizon Telescope at 230 and 345 GHz can probe further inner regions where the jet might interact with such hot dusty ionized gas.
• Comparison to the collimation profile of the 3C 273 jet with the velocity field obtained by long-term monitoring observations at 15 and 43 GHz shows that the 3C 273 jet is already accelerated to its apparent maximum velocity at the upstream of the collimation break, which is apparently different from the M87 jet where the ends of the collimation and acceleration zones are co-located.
• The collimation break of 3C 273 is located nearly several-to-ten times further than the estimated SGI location.Furthermore, we compared all available locations of the jet break from the literature with the SGI locations derived from M BH − σ d relation, clearly showing that the transition of the jet collimation is governed not only by the SGI of the central black hole.
• Relatively further location of the jet collimation break discovered for 3C 273 seems broadly consistent with the statistical trend of radio galaxies that higher excitation AGN sources tend to end their jet collimation at more distant locations from the central black hole.
Our results demonstrate that ultra-high-resolution observations with new millimeter facilities can provide many new insights for high-powered AGN which were not previously suitable to study jet collimation.Further study using full polarimetric data and Faraday rotation synthesis can reveal the distribution of the three-dimensional magnetic field and the circumnuclear medium in the active collimation region (e.g., Park et al. 2019a).These explorations with polarimetric analysis using presented multifrequency data including GMVA+ALMA will be presented in a forthcoming paper.

Figure 1 .
Figure 1.uv-coverage of the GMVA observation at 86 GHz.Baselines to ALMA are shown in red.Data are averaged at 60 s.

Figure 2 .
Figure 2. The SNRs 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.
. The overall jet morphology is consistent with previous ground-based (MOJAVE: Lister et al. 2018; VLBA-BU-Blazar 9 : Jorstad et al. 2017) and Space VLBI observations (Radio Astron: Bruniet 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 extended two components around 4 and 7 mas from the core.These overall morphology is consistent with a close-in-time image in the database of the VLBA-BU-Blazar program observed a week ago(March 19, 2017).

Figure 3 .
Figure3.The inner-jet images 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 resolutions in uniform).The image restored with a circular beam of 57 µas (labeled 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.The contours are multiplies by the factor of 2, from the lowest contour level of 6.5 mJy/beam.

Figure 4 .
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 the left panel of Figure 3.
Frequency pairs; (2) Total number of combinations; (3) Pixel size of re-gridded images used to measure offsets; (4) Mean offset in x-direction; (5) Mean offset in ydirection; (6) Standard deviation of the offsets in x-direction; (7) Standard deviation of the offset in y-direction.ofthe synchrotron emission, magnetic field, and particle density distributions.Assuming a conical jet having 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 byBlandford & Königl (1979), which was confirmed in many sources by previous core shift measurements (see alsoPushkarev et al. 2012, and references therein).For 3C 273,Lisakov et al. (2017) has performed multi-epoch 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.

Figure 5 .Figure 6 .
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 Table1.The lowest contour levels are 6.5 mJy/beam at 43 GHz, 8.2 mJy/beam at 22 GHz, 7.1 mJy/beam at 15 GHz, and 11.6 mJy/beam 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.
Figure 7.The core position along the large-scale jet axis as a function of the frequency based on the fitted power-law model.The core positions, denoted by the black dots, are relative to the 86 GHz core.The red line shows the bestfit power-law model, and the black line shows the best-fit location of the central black hole.The error bars of the bestfit core positions and the shaded area of the best-fit model indicate the 68% confidence interval estimated from the χ 2 surface of the fitted parameters (Section 5.4).
Nature of the Jet Collimation in 3C 273

Figure 8 .
Figure 8.The radius and PA of the 3C 273 jet as a function of the de-projected distance from the central black hole.The de-projected distance from the central black hole is derived using the best-fit core shift model at a viewing angle of i = 9 • (see Section 1) for the black hole mass of MBH = 2.6 × 10 8 M (Gravity Collaboration et al. 2018).(The upper panel ): the jet radius profile.The colored circles and associated error bars indicate the mean values and standard deviations of the measurements, respectively.Different colors show measurements from different observations, as shown in the legend.The solid and dashed lines show single power-law fits (r ∝ z a ) to the inner side, semi-parabolic (a = 0.66), and to the outer side, conical/hyperbolic (a = 1.28), respectively.The surrounding shaded areas represent the uncertainty of the single power-law fits (1, 2, 3 × σ).The green line represents the best-fit broken power-law function, giving the best-fit distance of the jet shape break of z b = 1.1×10 7 Rs shown on the vertical dotted line.(The lower panel ): the jet PA profile.The color convention is the same as in the upper panel.The horizontal dashed line shows PA = −138 • , which is broadly consistent with the axis of the large-scale jet (e.g.,Davis et al. 1985;Conway et al. 1993).The gray shade across both panels shows the region excluded from the fitting analysis because of a large bending seen in the jet (see the details in Section 5.5).

Figure 9 .Figure 10 .
Figure9.One of the top-set images from HSA 15 GHz observation (gray contours) overlaid with the jet streamline of the 3C 273 jet (filled and open circles).The excluded regions from our fitting are indicated by gray-filled circles, showing a stable complex and bending double-ridge-like morphology not seen in the rest of the global jet (see Section 6.1).A black dashed line indicates the global position angle of the jet estimated by the large scale jet image (e.g.,Davis et al. 1985).

Figure 11 .
Figure 11.The jet radius and PA profiles of the 3C 273 jet, extended from Figure 8 to include measurements of the large-scale jet beyond the deprojected distance of ∼200 kpc traced with VLA 5 GHz and MERLIN 1.6 GHz.
Figure12.The relation between the black hole mass and the de-projected distance of the jet collimation break from the central engine.The red star shows our results from the 3C 273 jet, and the cyan circles show results from other sources reported in the literature.The red horizontal bar indicates the region of a persistent feature seen near the break distance with a complex bending double-ridge-like morphology potentially associated with the transition in the jet shape (see Figure9and Section 6.1).To compare with the location of the sphere of the gravitational influence (SGI), we show the SGI locations of the classical bulges and elliptical galaxies in gray squares taken fromKormendy & Ho (2013).The black solid line and gray shaded area are the best-fit model and the 1-σ error range of the relation between the black hole mass and SGI location, transformed from the MBH − σ d relation inKormendy & Ho (2013).The dotted lines show the actual distances in parsecs.The results of other sources are obtained from the following references; M87: Asada & Nakamura 2012, NGC 6251:Tseng et al. 2016, NGC 4261:  Nakahara et al. 2018, Cygnus A: Nakahara et al. 2019, NGC 1052: Nakahara et al. 2020, NGC 315: Park et al. 2021, 4C38.41:Algaba et al. 2019, 3C 264: Boccardi et al. 2019, 1H 0323+342: Hada et al. 2018, and other sources: Kovalev et al. 2020.Note that we show two cases of different black hole masses for 1H 0323+342 (see Section 6.2).

Figure 13 .
Figure 13.Upper panel: Radial profile of jet in 3C 273 with different models of power-law functions.Measurements data are same as Figure 8. 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-colored squares correspond to the deviations from the lines in the upper panel.

2
Li et al. (2022)reported a new estimation of BH mass of ∼ 10 9 M 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 GravityCollaboration 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.values of M BH = 2.6 × 10 8 M and a viewing angle of θ = 9 •3 for the 3C 273 jet.

Table 1 .
3C 273 and calibrators (3C 279) were observed for a track of ∼16 hours, almost a full track for the given observing array.ALMA participated for ∼5.2 hours 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.

Table 1 .
Summary of observations of 3C 273

Table 2 .
List of Imaging ParametersFreq.Sys err

Table 3 .
The positional offsets of the core between different frequencies