First GMVA observations with the upgraded NOEMA facility: VLBI imaging of BL Lacertae in a flaring state

We analyze a single-epoch Global mm-VLBI Array (GMVA) observation of the blazar BL Lacertae (BL Lac) at 86 GHz from April 2021. The participation of the upgraded, phased Northern Extended Millimetre Array (NOEMA) adds additional sensitivity to the GMVA, which has facilitated the imaging of BL Lac during an unprecedentedly strong $\gamma$-ray flare. We aim to explore the nature of the inner subparsec jet of BL Lac and the impact of the NOEMA participation in the observation. For the data reduction, we employed two advanced automatic pipelines: rPICARD for the flux density calibration as well as the model-agnostic signal stabilization and GPCAL for the antenna leakage calibration. The conventional hybrid imaging (CLEAN + amplitude and phase self-calibration) was applied to the calibrated visibilities to generate final VLBI images. We performed a ridge-line analysis and Gaussian model-fits on the final jet image to derive the jet parameters. In our data, the presence of NOEMA improves the image sensitivity by a factor of 2.5. The jet shows a clear wiggling structure within 0.4 mas from the core. Our ridge-line analysis suggests the presence of a helical jet structure (i.e., a sinusoidal pattern). Six circular Gaussian components were fitted to the inner jet region. We estimated an apparent brightness temperature of $\sim$3 $\times$ 10$^{12}$ K in the two innermost components. They are likely to be highly boosted by relativistic beaming effect. We find four significant polarized knots in the jet. Interestingly, two of them are located in the core region. Finally, we suggest a number of physical scenarios to interpret our results.


Introduction
At millimeter wavelengths, very long baseline interferometry (VLBI) provides a sub-milliarcsecond scale (sub-mas) resolution, which facilitates the study of the jet-launching in the central region of active galactic nuclei (AGN; Krichbaum et al. 1998).Powerful radio jets launched in the vicinity of supermassive black holes (SMBHs) have been the primary target of VLBI studies owing to their compact, sub-pc scale structures, and bright non-thermal emission across the entire electromagnetic spectrum (Zensus 1997;Blandford et al. 2019).However, our understanding of the jet launching region remains limited due to insufficient angular resolution and jet opacity (i.e., the so-called core shift; Lobanov 1998).The inner jet region is thought to be optically thick at radio due to synchrotron self-absorption.By observing at higher frequencies (e.g., 86 GHz), we can peer deeper into the near-horizon scales of the jet structures.
The Global mm-VLBI Array 1 (GMVA) allows us to explore the inner subpc-scale jets in more detail (e.g., Baczko et al. 2016; Input files for Appendix A are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr(130.79.128.5) or via https:// cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/680/L3 1 https://www3.mpifr-bonn.mpg.de/div/vlbi/globalmm/Boccardi et al. 2016;Hodgson et al. 2017;Kim et al. 2019;Casadio et al. 2021;Paraschos et al. 2022).The GMVA operates at 86 GHz (3.5 mm) and is capable of reaching a very high angular resolution of ∼50 micro-arcsecond (µas).The presence of sensitive and large antennas in the array is the key for better image fidelity and thus high-quality scientific output.The NOrthern Extended Millimetre Array2 (NOEMA) is an interferometer that had 11 antennas at the time of the observations reported here (12 today) that can either be operated in local interferometer mode or phased array mode.The 2SB dual polarization receiver bands of 4 × 7744 MHz processed by the PolyFiX correlator allow phased array operation over the full bandwidth.Based on the measurements in April 2021, NOEMA has a system equivalent flux density (SEFD) of 150-250 Jy with a typical system temperature (T sys ) of ∼80 Kelvin (K) at 86 GHz.This is a major upgrade from the old six-element Plateau de Bure Interferometer (PdBI) that has participated in GMVA observations since 2003.NOEMA has twice the collecting area and four times the spectral bandwidth of the PdBI.In the phased array mode, data rates have passed from maximum 1 GBit/s (PdBI, correlator limited) to 64 GBits/s (NOEMA, recorder limited; the correlator can do up to 128 GBits/s).In April 2021, NOEMA joined regular GMVA observations for the first time.
BL Lacertae (BL Lac; 2200+420) is a prototypical source of the BL Lac class of objects and it is a nearby AGN: its redshift (z) is z ∼ 0.0686 (Vermeulen et al. 1995), corresponding to an image scale of 1.3 pc mas −1 (assuming H 0 = 71 km Mpc s −1 , Ω Λ = 0.73, and Ω m = 0.27).The jet of BL Lac is a strong γ-ray emitter that is variable with time and has been linked to (sub)pc-scale jet activity such as the propagation of moving blobs and shocks (Marscher et al. 2008;Wehrle et al. 2016).The morphology of the jet is known to be very complicated, with prominent bending or wiggling structures.The origin of such behavior is still a matter of debate.Possible interpretations invoke for instance, a helical jet model (Denn et al. 2000), jet precession (Stirling et al. 2003), MHD waves (Cohen et al. 2015), and kink instabilities (Jorstad et al. 2022).It was recently reported (Pushkarev et al. 2023) that the electric vector position angles (EVPAs) of the mas-scale jet are well-aligned with the local jet axis.At higher frequencies ≥43 GHz, however, the appearance of bright polarized emission seems erratic in the upstream region around the radio core (e.g., Marscher et al.Ratio (b) SEFD (c) (mJy beam −1 ) (Jy) All antennas 2.35 -w/o NOEMA (d)  5.92 2.5 ∼190 Notes. (a) Measured at the map center (a region of 4 mas × 4 mas). (b) Ratio of the noise level to the one with all antennas. (c) Overall SEFDs of the station in the data. (d) All antennas but without NOEMA.
2008; Wehrle et al. 2016).These signatures favor the presence of the dominant toroidal magnetic fields in the jet and multiple emitting regions in the unresolved upstream jet region (Hovatta et al. 2012).

GMVA observations in April 2021
In this work, we focus on the analysis of BL Lac that was observed as a calibrator in the project MB018A (PI: B. Boccardi, main target: NGC 315).The observations were performed on 24 April 2021.The duration of the observation was ∼15 h and the actual on-source time on BL Lac was ∼1.6 h; in the case of NOEMA, it was ∼1.1 h.The number of phased antennas for NOEMA was 10 in this experiment.The data consists of 8 IFs over 86.1-86.4GHz with a bandwidth of 64 MHz for each.Overall, 17 radio telescopes participated in the observation: Effelsberg (EF), Onsala (ON), Metsähovi (MH), NOEMA (NN), Pico Veleta (PV; the IRAM 30m telescope), and Greenland (GLT), along with eight VLBA antennas and three KVN antennas.We used the novel rPICARD pipeline (Janssen et al. 2019) for a reproducible data reduction and GPCAL (Park et al. 2021a) for polarization calibration.In Appendix A, we summarize the detailed steps of the data analysis.We also compare the results between the standard data analysis done in AIPS (Greisen 2003) with the new analysis pipeline rPICARD, which relies on CASA (CASA Team et al. 2022).
Figure 1 shows the fully calibrated visibilities of the data.The NOEMA baselines are marked in red color.The shortest baselines include NOEMA-Effelsberg (660 km), the longest baseline is NOEMA-VLBA_MK (Hawaii, 10 670 km).From simple visual inspection of the visibility plots, we can easily find that the groups of the NOEMA visibilities are very compact without large scattering in amplitude/phase, compared to other baselines in the observation.This indicates that they have  smaller visibility errors and higher signal-to-noise ratio (S/N).Indeed, we find that the overall S/Ns of the NOEMA datasets are significantly higher than the others by a factor of at least ≥3 (see Fig. 2).In addition to that, we refer to Table 1 for the impact of NOEMA on imaging sensitivity.

Total intensity structure of the jet
Figure 3 shows a final CLEAN image of the jet obtained from the 2021 GMVA observations.The overall data quality of the VLBA was not the best in the observation and we were only able to image the jet down to ∼0.4 mas from the core due to limited sensitivity.In our maps, the source structure is spatially resolved.We find the common core-jet morphology (e.g., Jorstad et al. 2017), but, interestingly, it displays a clear wiggling (i.e., a multiple-bending feature) structure determined by ridge line analysis (see below).Such behavior has already previously been reported in the source (e.g., Cohen et al. 2015), but at parsec scales on the sky plane.
There is a persistent, bright spot in the jet at a core separation of r = 1.5 mas; Weaver et al. (2022) identified this component as B15.The position angle of B15 (i.e., −166.4 • , CW from 0 • at north to −180 • at south) and the jet flow shown in Fig. 3 (i.e., −172 • on average) are different which suggests a misalignment of the jet axis between subpc and pc scales.
We also fit the data with circular Gaussian components (see Appendix B, for their physical parameters).In both the CLEAN and model-fit images, the overall jet morphology are consistent with each other.From visual inspection, however, they look slightly different.To further examine the differences and measure the jet properties, we draw jet ridge lines by using the two images (see Appendix C, for details of the analysis).Figures 3  and 4 show the results (hereafter, CN ridge for the CLEAN jet and MD ridge for the model-fit jet).As evident from the comparison, the wiggling pattern appears more prominent in the CN ridge than the MD ridge.We consider that the model-fit image describes the jet flow with a number of Gaussian components of certain sizes focused on brightest regions along the jet downstream, whereas the CLEAN process models the jet image with a set of many delta-components (i.e., point sources) that reflects faint emission regions more effectively.Hence, any resolved weak source structures can be better highlighted in the CLEAN image, which result in the discrepancy.This suggests that the CN ridge shows the jet structure in more detail.The reduced Chisquared values of the CLEAN and model-fit images (i.e., ∼1.2 and ∼1.8, respectively) also support this idea.In this context, the MD ridge could be considered as a line that just connects the representative (or brightest) regions of the flow.The jet of BL Lac is known for its complex, wiggling structure that varies over time (e.g., Cohen et al. 2015).Such behavior can be clearly seen in the image with the CN ridge line.The CN ridge shows multiple bending points at ∼0.11 mas, ∼0.23 mas, and ∼0.37 mas from the core.In addition to that, interestingly, the CN ridge line looks rotating around the MD ridge line.The observed features lead us to consider a helical jet structure that theoretical simulations have predicted (e.g., Meier et al. 2001;Mizuno et al. 2011).In Fig. 4, we present several jet parameters estimated by L3, page 3 of 12 the ridge lines.Overall, the distance between the two ridge lines (∆d ridges ) is below the level of angular resolution of the data and the sizes of the Gaussian components.This prevents us from making an argument quantitatively.However, there is an obvious systematic behavior in ∆d ridges that is sinusoidal along the downstream.This might be a hint of the helical jet structure, rather than just a coincidence.We fit a simple power-law to the measured jet width (W): W(r) = ar k where a and k being free parameters.To consider some of the known issues with resolution limit near the core and low S/N at the tail part of the jet, we checked the fitting results in three different areas: (A) the whole jet ridge line, (B) r ≥ 0.07 mas, and (C) between 0.07 mas and 0.28 mas.The best-fit parameters of the fits are on average: a = 0.21 ± 0.02 and k = 0.74 ± 0.05 (with 1σ uncertainty); the three fits resulted in the index k of 0.73-0.75,thus remaining consistent with respect to each other within the errorbars.Previous works measured the collimation profile of the jet at various observing frequencies with different approaches: stacked VLBI images (Kovalev et al. 2020;Casadio et al. 2021) and model-fit components of single-epoch images over a longterm period (Burd et al. 2022).Their results on the inner jet region span a wide range; k spans the interval 0.5-1.2,larger at higher frequencies.Our estimate (i.e., k ∼ 0.74) is in the middle of this range.This tendency suggests that the jet structure at subpc scales is highly complicated.At higher observing frequencies, the source is better resolved and thus the structural complexity (or time-variable local jet axis to left and right, which is probably magnified by projection effects) could better be reflected horizontally on the image plane.Then, the stacking of multi-epoch, highresolution images of such an extreme jet could make its width much wider than a single-epoch image.Lister et al. (2013) considered that in a conical jet model, a single-epoch image may show us a part of the jet rather than the entire jet cross section.However, it is unclear whether the subpc-scale jet of BL Lac is conical (k ∼ 1.0) or parabolic (k ∼ 0.5), considering those different estimates found by the previous works.In our future works, we plan to address this issue in more detail by using VLBI observations of BL Lac with the Event Horizon Telescope (EHT) at 230 GHz.We also note that the overall pattern of our CN ridge line already looks quite similar to the ridge line presented in Casadio et al. (2021) that was obtained from a stacked GMVA image of BL Lac at 86 GHz (i.e., their Fig. 2).
Overall, the apparent jet opening angle (Θ app ) decreases with the distance from the core.Using a single-epoch observation, Hada et al. (2018) and Park et al. (2021b) also found the same trend in the parabolic jet regions of the Narrow-line Seyfert 1 Galaxy 1H0323+342 and the radio galaxy NGC 315, respectively (but see also Ros et al. 2020, for opposite case in a blazar).The apparent opening angle can be varied due to changes in the intrinsic opening angle (Θ int ) and/or the jet viewing angle (θ jet ) which is attributed to the projection effect (e.g., Homan et al. 2002): Θ app = Θ int / sin(θ jet ).We find two local peaks in Θ app at r ∼ 0.12 mas and 0.28 mas.The first one is accompanied by a short dip in W at ∼0.1 mas and a strong emitting component around it.Thus, it might probably be caused by changes in the pressure of the plasma flow.In the case of the second one, we suspect changes in θ jet values around that region.A continuous, increasing pattern in flux density where 0.2 < r < 0.3, coincides with the second peak.This might be caused by smaller θ jet values in this region.mas from the core.The fractional polarization of the four knots are roughly: ∼2% (P1), ∼1% (P2), ∼4% (P3), and ∼17% (P4); measured at their central regions.We refer to Appendix D for details of the polarization detection and parameters of the four emission regions.In the case of the core (or the map peak), it is around 0.5%.The overall EVPAs in the jet are parallel to the local jet direction; i.e., the P3 and P4 knots aptly describe the flow axis.However, the two knots near the core in the upstream region (i.e., P1 and P2) seem to be different.The angular resolution of our data is insufficient to resolve this region and, thus, the detailed jet axis and structure are hidden.Pushkarev et al. (2023) presented a stacked polarization image of the jet of BL Lac at 15 GHz, by using 139 VLBA images of the source spanning over 20 years.In their image, the jet EVPAs are well aligned with the flow down to, for instance, ∼6 mas from the core.The downstream knots in our image (i.e., P3 and P4) confirm that such behavior persists up to about 0.05-0.1 mas from the core.This is also consistent with Rani et al. (2016).Faraday rotation is known to be very weak in the jet on parsec scales (Hovatta et al. 2012); however, it is unclear at the upstream region where the radio core appears.The EVPAs of P3 and P4 suggest the presence of strong toroidal magnetic fields in this subpc-scale region of the jet and probably also further downstream (i.e., on pc scales).The stacked jet image of Pushkarev et al. (2023) supports this; they used 139 VLBA maps of the jet at 15 GHz, which were not Faraday-corrected.This might be an indication of the presence of a helical magnetic field configuration in this inner jet (e.g., Nakamura et al. 2010).

Linearly polarized emission
At lower frequencies, in general, the cores of blazar jets show a single polarized knot that has no significant displacement from it, except when a strong disturbance propagates down the flow (e.g., Kim et al. 2018; see also Issaoun et al. 2022; L3, page 4 of 12 Jorstad et al. 2023, for VLBI images of blazar jets with the extended core polarization structures at 230 GHz).In our image, there are two polarized knots P1 and P2, near the core with a slight offset from the map center (i.e., ∼35 µas for P1 and ∼30 µas for P2).This is an unusual case.However, it seems that these features are also present in VLBI images of some previous works (i.e., Marscher et al. 2008;Rani et al. 2016;Casadio et al. 2017); in particular, both knots can be clearly seen in Casadio et al. (2017).In addition to that, we checked a near-in-time (i.e., 5 April 2021) 86 GHz VLBA image of the jet offered by the Boston University group 3 (BU; see Appendix E).Although this data is resolution-limited, we see hints of the P1 and P2 knots around the core with the similar EVPA patterns in the VLBA image.It is unclear what those features indicate without more detailed information on the core region (e.g., intrinsic EVPAs and opacity).Hovatta et al. (2012) could not measure rotation measure (RM) in the core of BL Lac due to the blending effect of multiple components; they employed the λ 2 -fit to find RM values at cm-wavelengths.Indeed, our image shows such feature (i.e., P1 and P2).At higher frequencies (i.e., 15-43 GHz), Gómez et al. (2016) found RM values of a few thousands of rad m −2 (up to 3000 rad m −2 ) around the core.Thus, the presence of strong Faraday rotation cannot be ruled out in this upstream region of the jet and high resolution multifrequency mm-VLBI observations at, for instance, >43 GHz are required to shed light on this region.It is worthwhile noting that just three days after our observation, there was the strongest γray outburst that is unprecedented in this source over the past ∼15 years (see Appendix F).

On the geometry of the subpc-scale jet
In this section, we suggest a model of the jet geometry in a qualitative way, considering the observed structural features.Our findings favor the hypothesis of helical jet structure (e.g., Meyer 2018).Although it is resolution-limited, the systematic sinusoidal pattern in ∆d ridges looks obvious and this can naturally be explained by the helical jet structure.In this view, the two ridge lines might actually show us the helical jet streamline (i.e., the CN ridge) and the changing axis of the helix (i.e., the MD ridge), respectively.We found three bending points in our analysis.The jet is assumed to have a small viewing angle (e.g., θ jet < 10 • ) and flow downstream helically.Thus, there could be changes in the direction of the actual flow (i.e., θ streamline ) around those bending regions, with respect to the line of sight (i.e., decreasing and increasing).We notice that the flux density increased between 0.2 mas and 0.3 mas in r.An enhanced beaming effect could be responsible for the rise (e.g., Raiteri et al. 2017), thus meaning that the flow (i.e., streamline) was approaching (or better aligned) to the observer in that region.If this is the case, we expect that the jet was initially approaching, and then it turned into receding (by the helical motion) at around r ∼ 0.1 mas.Beyond ∼0.4 mas, the jet disappears quickly and this can be found frequently in previous works (e.g., Casadio et al. 2021).This might be due to a relatively large change in the axis of the jet helix θ jet towards a large angle.
Details of the upstream region (i.e., r < 0.1 mas) are unclear due to the limited angular resolution.We found the observed brightness temperature (T b,obs ) of ∼3 × 10 12 K in this upstream region (C0 and C1).T b,obs can be inferred from the Doppler 3 https://www.bu.edu/blazars/index.htmlfactor (δ) and the intrinsic brightness temperature (T b,int ) as T b,obs = T b,int δ/(1 + z), where z and δ being redshift and δ = [Γ(1 − β cos θ jet )] −1 , respectively; Γ is the Lorentz factor (Γ = 1/ 1 − β 2 ) and β is the jet speed in units of the speed of light.Using the typical δ values of the jet reported by Jorstad et al. (2017) and Weaver et al. (2022): namely, 7.5 and 4.3, respectively, we estimate T b,int of ∼4 × 10 11 K for δ = 7.5 and ∼7 × 10 11 K for δ = 4.3.These estimates are higher than the equipartition limit (T b,eq ) of ∼5 × 10 10 K (Readhead 1994) that assumes equipartition in energy density between the magnetic field and emitting particles; Liodakis et al. ( 2018) recently reported T b,eq to be ∼2.8 × 10 11 K (i.e., the maximum T b,int measured independently of the assumption of equipartition).This assumption has been common practice in the relevant studies (Homan et al. 2021).In this context, we suggest two scenarios to interpret our T b,obs values found in C0 and C1.Firstly, the emission regions might be particle-dominated during the time of our observation (i.e., the flaring state).Secondly, if we follow the assumption of equipartition (i.e., T b,int ∼ T b,eq ), then δ should be ∼64 for Readhead (1994) and ∼11.5 for Liodakis et al. (2018).This indicates the presence of highly enhanced Doppler boosting in these regions.For C1, we consider that opacity is lower than C0 (see e.g., Gómez et al. 2016, for a spectral index map of the jet), thus suggesting lower δ values than the above estimates (e.g., ∼32 and ∼8, respectively).Using a GMVA survey at 86 GHz, Nair et al. (2019) presented a statistics on T b,obs for a large source sample.Interestingly, the T b,obs values we found in the upstream region, are about two times larger than the highest value of Nair et al. (2019), namely: T b,obs ∼ 1.3 × 10 12 K.Thus, the upstream region of the jet is likely more aligned to the line of sight with smaller viewing angle.We also note that if the jet is viewed at a very small viewing angle (i.e., close to ∼0 • ), then the projected EVPAs can be less dependent on the jet orientation (e.g., Jorstad & Marscher 2003), which might be reflected through the P1 and P2 knots in the core region.Lastly, we further add a possibility of the recollimation-shock (RCS) scenario for the core polarization, considering the structured linear polarization of RCS (see e.g., Marscher 2016).

Coincidence with the strongest γ-ray outburst
Luckily, BL Lac was in a peaking stage of its historically huge γray flare and it peaked just three days after our observation.The enhanced jet activity at r < 0.1 mas in terms of polarized emission and T b,obs , is a strong evidence that the origin of the γ-ray event was in this upstream region.Interestingly, the source was also flaring at 230 GHz in April/May 2021, based on the Submillimeter Array (SMA) database4 ; due to limited cadence, it is unclear when the radio flare peaked exactly, but it seems to be late-April or early-May 2021 with a level of ∼6 Jy.Using a spectral energy distribution (SED) modeling, Sahakyan et al. (2023) found that size of the emitting region tends to be smaller during a flaring state, which in our case, favors C1 as the main production site of the γ-ray event.However, the core region (i.e., C0) could also be responsible for the γ-ray emission.MacDonald et al. (2015) proposed the ring of fire model that assumes the presence of a synchrotron-emitting ring of electrons providing soft seed photons for the inverse-Compton process in the jet.The two polarized knots (i.e., P1 and P2) might indicate the ring that is a shocked portion of the jet sheath.Indeed, the transverse emission profiles in MacDonald et al. (2015, i.e., their Figs. 7 and 8) show a similar pattern to the core region of our image with P1 and P2.We assume that the strongest γ-ray outburst of BL Lac was not an orphan event owing to the presence of a near-in-time radio counterpart at 230 GHz (SMA).This could be an indication that P1 and P2 are located downstream of the core (see MacDonald et al. 2015, for an orphan γ-ray flare with the ring located upstream of the core).Then, this scenario suggests that the γ-ray outburst is originated near the C0 region, probably when a moving feature (e.g., Kim et al. 2022) with highly accelerated electrons passes through the ring area.Alternatively, one may also consider relativistic magnetic reconnection as the underlying physical process of the γ-ray outburst.The reconnection dissipates magnetic energy and accelerate particles efficiently in the jets.Zhang et al. (2022) reported that a fast, strong γ-ray flare can be produced by plasmoids (or mini-jets) in the reconnection region.Such plasmoids are thought to show enhanced polarized emission and a large swing in polarization angle caused by plasmoid mergers, which might be reflected in our image as the polarized features in the core region (see also Jorstad et al. 2022, for discussion on a link between magnetic reconnection and kink instability in the jet of BL Lac).

Summary
We have analyzed the first GMVA data, where the upgraded NOEMA interferometer joins as a phased array.From the observations taken in April 2021, we studied the BL Lac data in this work.Fortunately, the observation coincides with the strongest γ-ray flare of this source, which makes an examination of the data all the more interesting.The participation of NOEMA improved the imaging sensitivity largely (i.e., by a factor of ∼2.5 in our data) and the capability of detecting small-scale source features.The advanced calibration pipelines worked on the data nicely.In particular, we found that the visibility phase errors can be reduced significantly with rPICARD, compared to the conventional manual calibration with AIPS.Our GMVA image of the jet of BL Lac gives us hints of the complex nature of the jet at subpc scales from the central engine, such as the wiggling structure, unusual polarization signatures, and enhanced brightness temperatures.We suggest a number of scenarios to explain those interesting features in the presence of helical jet structure.However, we expect that the jet of BL Lac is highly variable with time.Thus, a monitoring of the source with highresolution (<50 µas) mm-VLBI array would be highly useful for shedding more light on the nature of the subpc-scale jet of BL Lac.Naturally weighted image of the jet produced by the conventional manual calibration with AIPS.The map is convolved with a circular restoring beam of radius 0.07 mas (bottom-left).The contours increase by a factor of √ 2 from 1.1 % to 99.6 % of the map peak (i.e., ∼0.67 Jy/beam).The rms noise level is around ∼3.24 mJy/beam.weather information).We carefully inspected the Tsys measurements of all the GMVA antennas and flagged and replaced any erroneous measurements (e.g., 999 K) with reasonable values by using linear interpolation.We found 13968448 and 11634432 visibilities in total from the rPICARD and manual AIPS calibrations, respectively.This means that rPICARD returns ∼20 % of more visibilities; the S/N cutoff in the fringe-fitting was set to a conservative value of 5.5 for rPICARD and 4.5 for AIPS, thus meaning that the threshold for AIPS was even lower than the case of rPICARD. 6For clarity, those visibilities in Figure A.1 are averaged with an interval of 600 seconds.The mean (ph) and standard deviation (1σ ph ) of the phases are: ph ∼ −0.9 • and 1σ ph ∼ 28.5 • for rPICARD and ph ∼ 1.8 • with 1σ ph ∼ 64.1 • for AIPS.These estimates suggest that rPICARD improved the phase stability by a factor of > 2; for comparison, we also provide an image of the jet generated from the same data, but using the manual AIPS calibration (see Figure A.2).
We have identified three features of rPICARD that are responsible for the improved data calibration over the standard AIPS procedures typically employed for the calibration of GMVA data: 1.The "coherence" calibration as part of the signal stabilization prior to solving for instrumental effects substantially reduces the impact of atmospheric phase turbulence and thus leads to more accurate solutions for instrumental calibration steps, such as telescope phase bandpasses.2. By combining all correlation products after the instrumental alignments (including RL phases and delays), higher sensitivities for the fringe search are obtained.3.For the standard AIPS calibration, a single solution interval is used for the fringe-fitting, namely, four minutes in the case of this data.If that interval is too long, atmospheric phase turbulence will not be corrected well.If it is too short, detections may be lost due to poor S/N.Here, rPICARD uses an adaptive solution interval to achieve an optimal calibration.The segmentation time is tuned for each individual scan over a wide search range (starting with five seconds for this data) by picking the shortest interval that yields the maximum number of detections.In some cases (see Janssen et al. 2019, for details), different solution intervals can be used for different antennas in the same scan, for instance, when all the data to single antenna has low S/N.The solutions obtained on different intervals are then combined into a single calibration table via interpolation.
Next to a single Measurement Set that contains the calibrated and raw data of all sources, rPICARD generates UVFITS files of the calibrated data for each source in the data after all the calibration steps are finished.Then, this uv-data was imaged with the software Difmap (Shepherd et al. 1994).We employed the conventional hybrid imaging approach: cycles of the CLEAN algorithm (Högbom 1974) plus amplitude and phase self-calibration, iteratively.To highlight source structures with low-surface brightness and reduce side-lobes, we applied a Gaussian uv-taper to the data at uv radius of 1.2 Gλ and the resultant beam size is 0.184 × 0.058 mas at −6.09 • .During the imaging processes, we used both the uniform and natural weighting schemes to find any missing fluxes of the jet.In this particular dataset, however, we found that the uniform weighting is better than the natural weighting in terms of the rms noise level: 1σ I,rms ∼ 2.35 mJy/beam for natural weighting and 1.29 mJy/beam for uniform weighting.The overall jet shape appears in the final image more smoothly and clearly with the natural weighting without small-scale artifacts that disturb the edges of the jet.Thus, we employ the naturally weighted image in the analysis of the jet structure.
For polarization calibration (i.e., solving for the D-term), we employed the Generalized Polarization CALibration pipeline (GPCAL; Park et al. 2021a).It employs both the Difmap and AIPS.This pipeline was developed recently to overcome a number of limitations of the conventional program LPCAL (Leppänen et al. 1995)  Amplitude -LCP (%) gram PolSolve).In the calibration, BL Lac was used as polarization calibrator.We divided the jet into four separate regions that cover all the CLEAN components in the final image.At such a high frequency (i.e., 86 GHz), polarization structures of the inner subpc-scale jet are likely to be complicated and the similarity assumption (e.g., Leppänen et al. 1995) may not hold.Thus, we employ the instrumental polarization self-calibration technique to reflect those complex polarization structures properly.
We found that some IFs and/or antennas are problematic during the imaging process.We lost the first IF (IF1) in EF and ON.Overall, IF5 was less sensitive than the other IFs.As expected from its high Tsys values (∼450 K on average), MH was very noisy throughout the whole scans and we found mean values of its D-term amplitudes for both RCP and LCP to be ∼30 %.In each antenna, the D-term amplitudes are quite variable with different IFs, particularly FD, NL, OV, and PT (see Figure A.3).However, the overall fitting processes with GPCAL yielded good results with a mean reduced chi-squared value of ∼1.7 for all the IFs.We also find that the D-terms of LCP are slightly higher than the ones of RCP on average.The mean amplitude of the D-terms of RCP over the whole IFs/antennas, is ∼12 % (∼11 % without MH).In the case of LCP, the mean amplitude is ∼13 % (∼12 % without MH).We further refer to Casadio et al. (2017) for previous D-term estimates of the GMVA antennas.For all the antennas except MH, the overall D-terms are below 20 %, with NL & PT having a number of IFs exceeding 20 %; during the data calibration, we noticed that these antennas (including MH) were much worse than the others in this dataset.The best D-terms were found in EF, GL, KP, and NN that range 3-6 % on average.In the polarization imaging, we excluded any IFs or antennas with the D-terms exceeding 20 %.We also note that both with averaged or separated IFs, those polarized features (i.e., P1-4) appeared in the final image.The latter was used in this work due to better rms noise level in polarization.We further compared this image to the one generated without NL and PT completely and there was no significant difference between them.
Although NOEMA joined the observation, quite a number of the other telescopes were not in good condition, affecting the overall data quality.MH and PV suffered from bad weather; in particular, no fringes were found to PV due to bad weather.Unfortunately, most of the VLBA antennas suffered from pointing and/or focus issues in this epoch (plus bad receiver performance in FD and PT).The final image of the jet was also modeled with multiple circular Gaussian components.tainties for their flux and position, respectively (e.g., Lister et al. 2009).The apparent brightness temperatures of the components are calculated following Issaoun et al. (2022).

Appendix B: Gaussian model-fit components
We used the model-fit image to measure the jet width (W).In the case of the CLEAN map, we found that the Gaussian modeling describes the transverse emission profile poorly in some regions due to some extra features.Plus, the profile of the width measured with the CLEAN map was erratic and a power-law could not fit to it.Thus, the model-fit image was used to find the jet width in this work.The width was estimated assuming the ridge line that we found above to be the central region of the jet.For every ridge point, we drew a line orthogonal to the line between two consecutive ridge points (i.e., at r and r −0.01 mas).Then, we fitted the Gaussian form to the profile of the orthogonal line.The FWHM of the Gaussian was found by 2 √ 2 ln2 × σ r , where σ r being the standard deviation (or the fitted Gaussian radius).We followed Pushkarev et al. (2009) to calculate the width and the apparent opening angle (Θ app ) of the jet: Beam and Θ app = 2 arctan(W/2r), with b Beam being the restoring beam size.
Figure 4 shows a number of jet parameters.Overall, we considered the 2σ standard deviation of the fit parameters as their uncertainties.The distance between the two ridge lines (∆d ridges ) was calculated using their x,y-positions on the image plane.To evaluate its uncertainty, we employed the error propagation of their positional errors.In the case of W, the errorbars were too small.However, we found that the position of the fitted Gaussian peak often has a small offset to the actual position of the brightest pixel in each of the transverse slices.We added this offset to the uncertainties of W. The CN ridge was used to find a flux profile along the jet flow.As mentioned in Section 3.1, we consider that the CN ridge line describes the jet structure in more detail.10% of the flux densities were assumed to be their uncertainties.

Fig. 1 .
, visibility amplitude (middle), and phase (bottom) of the fully calibrated data.The NOEMA baselines are indicated by red color.For uniform weighting, a beam size of 0.165 × 0.037 mas at −4.72 • is obtained.

Fig. 2 .
Fig. 2. S/N plot of the calibrated data.The data were scan-and frequency-averaged in this plot.

Fig. 3 .
Fig. 3. Jet of BL Lac at 86 GHz observed by the GMVA on 24 April 2021.The color scale and contours show the total intensity structure of the jet (natural weighting).The contours increase by a factor of 2 from 0.6% to 76.8% of the map peak (∼0.77Jy beam −1 ).The gray circle at bottom-left denotes the restoring circular Gaussian beam of 0.07 mas (=70 µas).The green and blue dots indicate the ridge lines measured with the model-fit and CLEAN maps, respectively.The yellow dashed line points to the direction of the known jet component B15 presented in Weaver et al. (2022).

Fig. 4 .
Fig. 4. Distance between the two ridge lines at every 0.01 mas from the core (∆d ridges ), transverse jet width (W), apparent jet opening angle (Θ app ), and a profile of flux density as a function of the distance from the core (r), shown from top to bottom.A simple power-law form is fitted to the width profile in three different areas divided by r = 0.07 mas and 0.28 mas to take into account any flattening feature resulting from resolution limits and low S/N.The opening angles are considered only when r ≥ 0.07 mas.The outer part of the emission profile (i.e., r > 0.15 mas) is enlarged separately.More details of the measurements can be found in Appendix C.

Figure 5
Figure5shows the distribution of linearly polarized emission of the jet.We found four significant polarized knots within r = 0.25 Fig. A.1.Visibility phases of the calibrated data: the conventional manual calibration with AIPS (A) versus the automatic calibration with rPICARD (B); they come before the imaging processes.The data were averaged into 600 second-bins for clarity.

Fig
Fig. A.2. Naturally weighted image of the jet produced by the conventional manual calibration with AIPS.The map is convolved with a circular restoring beam of radius 0.07 mas (bottom-left).The contours increase by a factor of Fig. A.3.D-terms of the GMVA antennas estimated by GPCAL from the data.From top to bottom: RCP/LCP amplitudes and RCP/LCP phases.The horizontal dashed lines in the amplitude panels indicate 20 %.
Fig. B.1.Gaussian model-fit image of the jet observed by the GMVA at 86 GHz.The Gaussian components are shown in red color.The map is convolved with a circular restoring beam of radius 0.07 mas (bottomright).The contours increase by a factor of √ 2 from 0.8 % to 72.4 %.

Table 1 .
Residual rms noise levels (σ rms ) of the final image.
Martí-Vidal et al. 2021d way most widely used in the community for the polarization calibration (see alsoMartí-Vidal et al. 2021, for another advanced CASA-based pro-

Table B .
1. Properties of the model-fitted Gaussian components.Positions of the components with respect to the core.b From North to South: CW (negative) and CCW (positive).c Observed brightness temperature. a