# Unravelling the Innermost Jet Structure of OJ 287 with the First GMVA+ALMA
Observations

Guang-Yao Zhao Instituto de Astrofísica de Andalucía-CSIC, Glorieta de la
Astronomía s/n, E-18008 Granada, Spain  José L. Gómez Instituto de Astrofísica
de Andalucía-CSIC, Glorieta de la Astronomía s/n, E-18008 Granada, Spain
Antonio Fuentes Instituto de Astrofísica de Andalucía-CSIC, Glorieta de la
Astronomía s/n, E-18008 Granada, Spain  Thomas P. Krichbaum Max-Planck-
Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany
Efthalia Traianou Instituto de Astrofísica de Andalucía-CSIC, Glorieta de la
Astronomía s/n, E-18008 Granada, Spain  Rocco Lico Instituto de Astrofísica de
Andalucía-CSIC, Glorieta de la Astronomía s/n, E-18008 Granada, Spain  INAF-
Istituto di Radioastronomia, Via P. Gobetti 101, I-40129 Bologna, Italy  Ilje
Cho Instituto de Astrofísica de Andalucía-CSIC, Glorieta de la Astronomía s/n,
E-18008 Granada, Spain  Eduardo Ros Max-Planck-Institut für Radioastronomie,
Auf dem Hügel 69, D-53121 Bonn, Germany  S. Komossa  Max-Planck-Institut für
Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany  Kazunori Akiyama
Massachusetts Institute of Technology Haystack Observatory, 99 Millstone Road,
Westford, MA 01886, USA  National Astronomical Observatory of Japan, 2-21-1
Osawa, Mitaka, Tokyo 181-8588, Japan  Black Hole Initiative at Harvard
University, 20 Garden Street, Cambridge, MA 02138, USA  Keiichi Asada
Institute of Astronomy and Astrophysics, Academia Sinica, 11F of Astronomy-
Mathematics Building, AS/NTU No. 1, Sec. 4, Roosevelt Rd, Taipei 10617,
Taiwan, R.O.C.  Lindy Blackburn Black Hole Initiative at Harvard University,
20 Garden Street, Cambridge, MA 02138, USA  Center for Astrophysics $|$
Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  Silke
Britzen Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121
Bonn, Germany  Gabriele Bruni INAF – Istituto di Astrofisica e Planetologia
Spaziali, via Fosso del Cavaliere 100, 00133 Roma, Italy  Geoffrey B. Crew
Massachusetts Institute of Technology Haystack Observatory, 99 Millstone Road,
Westford, MA 01886, USA  Rohan Dahale Instituto de Astrofísica de Andalucía-
CSIC, Glorieta de la Astronomía s/n, E-18008 Granada, Spain  Indian Institute
of Science Education and Research Kolkata, Mohanpur, Nadia, West Bengal
741246, India  Lankeswar Dey  Department of Astronomy and Astrophysics, Tata
Institute of Fundamental Research, Mumbai 400005, India  Roman Gold
CP3-Origins, University of Southern Denmark, Campusvej 55, DK-5230 Odense M,
Denmark  Achamveedu Gopakumar  Department of Astronomy and Astrophysics, Tata
Institute of Fundamental Research, Mumbai 400005, India  Sara Issaoun Center
for Astrophysics $|$ Harvard & Smithsonian, 60 Garden Street, Cambridge, MA
02138, USA  Department of Astrophysics, Institute for Mathematics,
Astrophysics and Particle Physics (IMAPP), Radboud University, P.O. Box 9010,
6500 GL Nijmegen, The Netherlands  NASA Hubble Fellowship Program, Einstein
Fellow  Michael Janssen Max-Planck-Institut für Radioastronomie, Auf dem Hügel
69, D-53121 Bonn, Germany  Svetlana Jorstad Institute for Astrophysical
Research, Boston University, 725 Commonwealth Ave., Boston, MA 02215, USA
Astronomical Institute, St. Petersburg University, Universitetskij pr., 28,
Petrodvorets,198504 St.Petersburg, Russia  Jae-Young Kim Department of
Astronomy and Atmospheric Sciences, Kyungpook National University, Daegu
702-701, Republic of Korea  Korea Astronomy and Space Science Institute,
Daedeok-daero 776, Yuseong-gu, Daejeon 34055, Republic of Korea  Max-Planck-
Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany  Jun Yi
Koay Institute of Astronomy and Astrophysics, Academia Sinica, 11F of
Astronomy-Mathematics Building, AS/NTU No. 1, Sec. 4, Roosevelt Rd, Taipei
10617, Taiwan, R.O.C.  Yuri Y. Kovalev Lebedev Physical Institute of the
Russian Academy of Sciences, Leninsky prospekt 53, 119991 Moscow, Russia
Moscow Institute of Physics and Technology, Institutsky per. 9, Dolgoprudny,
Moscow region, 141700, Russia  Max-Planck-Institut für Radioastronomie, Auf
dem Hügel 69, D-53121 Bonn, Germany  Shoko Koyama Niigata University, 8050
Ikarashi-nino-cho, Nishi-ku, Niigata 950-2181, Japan  Institute of Astronomy
and Astrophysics, Academia Sinica, 11F of Astronomy-Mathematics Building,
AS/NTU No. 1, Sec. 4, Roosevelt Rd, Taipei 10617, Taiwan, R.O.C.  Andrei P.
Lobanov Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121
Bonn, Germany  Laurent Loinard Instituto de Radioastronomía y Astrofísica,
Universidad Nacional Autónoma de México, Morelia 58089, México  Instituto de
Astronomía, Universidad Nacional Autónoma de México, CdMx 04510, México  Ru-
Sen Lu Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80
Nandan Road, Shanghai 200030, People’s Republic of China  Key Laboratory of
Radio Astronomy, Chinese Academy of Sciences, Nanjing 210008, People’s
Republic of China  Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69,
D-53121 Bonn, Germany  Sera Markoff Anton Pannekoek Institute for Astronomy,
University of Amsterdam, Science Park 904, 1098 XH, Amsterdam, The Netherlands
Gravitation and Astroparticle Physics Amsterdam (GRAPPA) Institute, University
of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands  Alan P.
Marscher Institute for Astrophysical Research, Boston University, 725
Commonwealth Ave., Boston, MA 02215, USA  Iván Martí-Vidal Departament
d’Astronomia i Astrofísica, Universitat de València, C. Dr. Moliner 50,
E-46100 Burjassot, València, Spain  Observatori Astronòmic, Universitat de
València, C. Catedrático José Beltrán 2, E-46980 Paterna, València, Spain
Yosuke Mizuno Tsung-Dao Lee Institute, Shanghai Jiao Tong University,
Shengrong Road 520, Shanghai, 201210, People’s Republic of China  School of
Physics and Astronomy, Shanghai Jiao Tong University, 800 Dongchuan Road,
Shanghai, 200240, People’s Republic of China  Institut für Theoretische
Physik, Goethe-Universität Frankfurt, Max-von-Laue-Straße 1, D-60438 Frankfurt
am Main, Germany  Jongho Park Institute of Astronomy and Astrophysics,
Academia Sinica, 11F of Astronomy-Mathematics Building, AS/NTU No. 1, Sec. 4,
Roosevelt Rd, Taipei 10617, Taiwan, R.O.C.  Tuomas Savolainen Aalto University
Department of Electronics and Nanoengineering, PL 15500, FI-00076 Aalto,
Finland  Aalto University Metsähovi Radio Observatory, Metsähovintie 114,
FI-02540 Kylmälä, Finland  Max-Planck-Institut für Radioastronomie, Auf dem
Hügel 69, D-53121 Bonn, Germany  Teresa Toscano Instituto de Astrofísica de
Andalucía-CSIC, Glorieta de la Astronomía s/n, E-18008 Granada, Spain

(Received 21 March, 2022; Accepted 27 April, 2022)

###### Abstract

We present the first very-long-baseline interferometric (VLBI) observations of
the blazar OJ 287 carried out jointly with the Global Millimeter VLBI Array
(GMVA) and the phased Atacama Large Millimeter/submillimeter Array (ALMA) at
3.5 mm on April 2, 2017. Participation of phased-ALMA not only has improved
the GMVA north-south resolution by a factor of $\sim$,3, but also has enabled
fringe detection with signal-to-noise ratios up to 300 at baselines longer
than 2 G$\lambda$. The high sensitivity has motivated us to image the data
with the newly developed regularized maximum likelihood imaging methods,
revealing the innermost jet structure with unprecedentedly high angular
resolution. Our images reveal a compact and twisted jet extending along the
northwest direction with two bends within the inner 200 µas that resembles a
precessing jet in projection. The component at the southeastern end shows a
compact morphology and high brightness temperature, and is identified as the
VLBI core. An extended jet feature that lies at $\sim$ 200 µas northwest of
the core shows a conical shape in both total and linearly polarized intensity,
and a bimodal distribution of the linear polarization electric vector position
angle. We discuss the nature of this feature by comparing our observations
with models and simulations of oblique and recollimation shocks with various
magnetic field configurations. Our high-fidelity images also enabled us to
search for possible jet features from the secondary supermassive black hole
(SMBH) and test the SMBH binary hypothesis proposed for this source.

galaxies: active – galaxies: individual (OJ 287) – galaxies: jets –
polarization – radio continuum: galaxies

††journal: ApJ††software: $\mathcal{AIPS}$ (Greisen, 2003), DiFX (Deller et
al., 2007), DIFMAP (Shepherd et al., 1995), eht-imaging (Chael et al., 2018),
ParselTongue (Kettenis et al., 2006), SMILI (Akiyama et al., 2017), Altair
(VanderPlas et al., 2018), Vega-Lite (Satyanarayan et al., 2017).

##  1 Introduction

The BL Lac type object OJ 287 ($z=0.306$; Stickel et al., 1989) is a well-
studied low synchrotron peaked BL Lac object (LBL) that has attracted great
interest as it shows quasi-periodic optical outbursts with a cycle of about 12
years. These outbursts appear to come in pairs with separations of one to two
years and have been suggested to originate due to the presence of a
supermassive binary black hole (SMBBH) system at its center (e.g., Sillanpää
et al., 1988; Lehto & Valtonen, 1996). According to this model, the observed
quasi-periodic double-peaked optical outbursts are triggered when the
secondary supermassive black hole (SMBH) impacts the accretion disk of the
primary in its orbit. Further advances of this model have accounted for
general relativistic effects and the parameters are also further constrained
with follow-up observations (e.g., Valtonen et al., 2008, 2011; Dey et al.,
2018). The model requires a compact binary with a major axis of the orbit of
0.112 pc (corresponding to an angular scale of $\sim$ 26 µas; e.g., Valtonen
et al., 2008), featuring a very massive primary BH of $1.8\times
10^{10}M_{\odot}$, and a secondary of $1.5\times 10^{8}M_{\odot}$ (e.g.,
Valtonen et al., 2012; Dey et al., 2018). This model is not only successful in
reproducing the observed light curves of OJ 287, but also in predicting impact
outbursts that were later confirmed by observations (e.g., Valtonen et al.,
2006, 2016; Laine et al., 2020; Komossa et al., 2020). Independent of the
binary model of OJ 287, dedicated multi-wavelength observation and modeling of
the OJ 287 (MOMO) project has led to the discovery of several bright flare
events and long-lasting deep fades, and monitoring spectroscopy of the last
two decades has established OJ 287 as one of the most spectrally variable
blazars in the soft X-ray band (e.g., Komossa et al., 2017, 2021a, 2021b,
2021c).

Another observational signature of OJ 287 is that the position angle (PA) of
the parsec-scale jet was found to be “wobbling” by previous very long baseline
interferometric (VLBI) observations (e.g., Tateyama & Kingham, 2004; Agudo et
al., 2012; Cohen, 2017; Britzen et al., 2018). Such changes of the inner jet
PA could also be explained by the SMBBH model (e.g., Dey et al., 2021), but
alternative models could not be fully ruled out. For instance, Agudo et al.
(2012) suggest instabilities coupled to the accretion disk as likely origin
for the non-periodic changes in the inner jet orientation. Britzen et al.
(2018) suggest the flux variation could be explained by viewing angle changes
and Doppler beaming effects of a precessing jet. The precession could be
driven by either the binary motion (e.g., Dey et al., 2021) or the Lense-
Thirring effect due to the misalignment between the BH spin and the accretion
disc (e.g., Chatterjee et al., 2020; Liska et al., 2018).

The massive central black hole (BH), the relatively low redshift, and the
bright close to line-of-sight relativistic jet also make OJ 287 one of the
nearest high-luminosity AGN in which the magnetic launching and acceleration
of jets can be studied through high-resolution VLBI observations. Two
competing scenarios have been proposed for the formation of relativistic jets.
The main difference between them is whether the magnetic fields are twisted by
the rotational energy of the BH (BZ model; Blandford & Znajek, 1977) or its
accretion disk (BP model; Blandford & Payne, 1982). It is also possible that
both mechanisms are at work (e.g., Chiaberge et al., 2000). In the innermost
region of the jet, the plasma flow is accelerated and collimated in the
presence of a spiral magnetic field, while the jet expands in width and
propagates downstream into the interstellar space. The disruption of the
accretion flow and the interaction with the ambient medium often result in the
formation of moving and standing shocks. The detailed process of jet
formation, acceleration, and collimation remains unclear as it requires
extremely high angular resolution to probe into the innermost region in the
vicinity of the central black hole.

High-resolution VLBI observations are ideal for probing the compact structure
near the central engine. Previous VLBI observations of OJ 287 have provided
key information on the parsec-scale structure and dynamics of the jet (e.g.
Hodgson et al., 2017; Cohen, 2017; Britzen et al., 2018). In particular, Gómez
et al. (2022) recently presented 22 GHz images of OJ 287 with unprecedented
angular resolution for the source obtained with the RadioAstron space-ground
VLBI observations. The images revealed a progressive bending of the inner jet
with increasing angular resolution by comparison with multi-band ground-based
VLBI images. The inner jet components show high brightness temperatures that
exceed the inverse Compton limit, indicating strong Doppler boosting in the
jet. The polarized images show electric vector position angles (EVPAs) aligned
with the jet axis, which indicates the jet has a predominantly toroidal
magnetic field. Multi-frequency analysis shows hints for a rotation measure
gradient across the jet, which suggests the VLBI core is threaded by a helical
magnetic field.

VLBI observations at wavelengths shorter than 7 mm hold the potential of
probing areas closer to the central engine that are optically thick at lower
frequencies (see e.g., Boccardi et al., 2017). Previous VLBI observations at
3.5 mm with the Global Millimeter VLBI Array (GMVA) show the existence of
quasi-stationary components and changes in the morphology and PA in the
innermost jet region (e.g., Hodgson et al., 2017). However, most of the
previous GMVA observations are limited in sensitivity due to typically shorter
atmospheric coherence times, lower antenna efficiencies, and thus higher
system equivalent flux densities (SEFDs) compared to longer wavelengths.
Participation of large sensitive stations in mm-VLBI observations are
desirable alongside with further developments of the instruments and
calibration methods (e.g., Rioja & Dodson, 2011; Rioja et al., 2017; Zhao et
al., 2018).

In this paper, we present the first VLBI observations of OJ 287 with the GMVA
and phased Atacama Large Millimeter/submillimeter Array (ALMA) on April 2,
2017. These observations are accompanied by a multi-wavelength campaign
including the first 1.3 mm observation of the source with the Event Horizon
Telescope (EHT; Event Horizon Telescope Collaboration et al., 2019a), the
results of which will be presented in a forthcoming paper. The campaign was
carried out during a major outburst event of OJ 287 in 2016-17 with the
largest X-ray outburst recorded so far (Komossa et al., 2017, 2021c) and the
first very high energy (VHE) flare detection (Mukherjee & VERITAS
Collaboration, 2017).

We summarize the details of the GMVA+ALMA observations and the methods we use
to calibrate, image, and analyze the data in section 2; we present our
observational results including total intensity and linear polarization images
in section 3; in section 4, we discuss the nature of the components in the jet
and possible constraints on the theoretical models, followed by a summary in
section 5.

##  2 Observations and data analysis

Figure 1:  $(u,v)$-coverage of the fringe-fitted interferometric visibilities
of OJ 287, observed with GMVA+ALMA on April 2, 2017 at 86 GHz. The baselines
to ALMA are plotted in red color and the other GMVA baselines are plotted in
blue.  Figure 2:  Self-calibrated visibility amplitudes (top) and phases
(middle) as a function of $(u,v)$-distance of the GMVA+ALMA observation of OJ
287 on April 2, 2017 at 86 GHz. The data were averaged every 15 seconds and
all channels in each IF are averaged. Over-plotted in orange are the fit to
the data of the reconstructed image obtained with SMILI. The bottom panel
shows the fringe signal-to-noise ratio as a function of $(u,v)$-distance, with
the data on ALMA baselines plotted in red and the other baselines in blue.

In this section, we describe the details of our 3.5 mm observations of OJ 287
with GMVA + ALMA, the data calibration procedure, and the methods used to
obtain subparsec-scale images of OJ 287.

###  2.1 Observations

We carried out high-resolution VLBI observations towards OJ 287 at 3.5 mm with
GMVA on April 2, 2017. These observations mark the first VLBI observations
with the phased-ALMA which consists of 37 ALMA antennas and is equivalent to a
70-meter dish (Event Horizon Telescope Collaboration et al., 2019a). The
participating stations also include 8 Very Long Baseline Array (VLBA) stations
and 5 European stations (Effelsberg, IRAM-30m, Metsähovi, Onsala, Yebes-40m).
The on-source time was around 375 minutes between UT 17 and UT 7 the next day
(April 3).

Most stations had good or typical weather conditions during the observation
except for the VLBA Mauna Kea (MK) and Pie Town (PT) stations, which resulted
in few fringe detections with limited signal-to-noise ratios (S/N) on
baselines to these two stations. No fringes were found on baselines to
Metsähovi due to a faulty backend setup. All data were recorded in full
polarization mode, with most stations recorded on a circular polarization
basis, while the ALMA data were converted from a mixed linear-circular basis
to circular polarization mode using PolConvert (Martí-Vidal et al., 2016).
Yebes-40m telescope recorded only left hand circular polarization (LCP). The
bandwidth and frequency range recorded are not the same at all stations111The
recorded bandwidth for each station is as follows: ALMA 32$\times$62.5MHz,
VLBA 2$\times$128MHz, most European stations 1$\times$512MHz. Only the common
frequency ranges among all participating stations are used in later
processing.

###  2.2 Data Reduction

Data correlation was performed with the DiFX correlator (Deller et al., 2007)
at the Max-Planck-Institut für Radioastronomie in Bonn, Germany. The final
correlated data have a total bandwidth of 232 MHz which were further divided
into four 58 MHz intermediate frequency (IF) bands.

The post-correlation dataset was then processed with the ParselTongue
(Kettenis et al., 2006) $\mathcal{AIPS}$ (Greisen, 2003) interface for fringe-
fitting and a priori amplitude calibration. We first performed parallactic
angle correction with the $\mathcal{AIPS}$ task, clcor 222We note the mount
types for IRAM-30m (Nasmyth-Left) and Yebes-40m (Nasmyth-Right) are different
from the rest of antennas in the array (altitude-azimuth). The Yebes-40m data
were not used for polarimetric analysis as they were only recorded in LCP.,
and manual phase calibration using short segments of data to remove
instrumental phase offset between different IFs. We then perform a global
fringe-fitting of the data using the task fring with a solution interval of 10
seconds and sub-intervals down to 2 seconds and by integrating over the whole
232 MHz bandwidth and averaging parallel-hand polarizations (RR & LL).

The $(u,v)$-coverage towards OJ 287 for all baselines with fringe detections
is shown in Figure 1. We note that the participation of ALMA has provided an
increase in the north-south resolution by a factor of $\sim$ 3 for
observations of OJ 287. ALMA has also significantly improved fringe detection
due to its high sensitivity (see Figure 2) with the maximum fringe S/N
reaching $\sim$ 350 at baselines longer than 1.5 G$\lambda$.

A priori amplitude calibration was performed in $\mathcal{AIPS}$ with the task
apcal by multiplying the system temperatures ($T_{\text{sys}}$) and gain
curves of each antenna. Opacity corrections were applied to stations that
measure the system temperatures with the noise diode method (VLBA &
Effelsberg). For ALMA, IRAM-30m, and Yebes-40m, the $T_{\text{sys}}$
measurements were performed using the hot / cold method and therefore already
included the opacity correction. The ALMA $T_{\text{sys}}$ values have also
taken into account the phasing efficiencies derived during the quality
assurance and PolConvert processes (e.g., Goddi et al., 2019). The cross-hand
phase and delay offsets of the reference station were calibrated using the
$\mathcal{AIPS}$ procedure, vlbacpol.

After the $\mathcal{AIPS}$ calibration, the data were averaged in time (with
an interval of 15 s) and frequency (with all channels within each IF averaged)
for further processing.

###  2.3 Imaging $\&$ Model-fitting

We performed imaging and self-calibration of the data independently with three
different imaging softwares: DIFMAP, eht-imaging, and SMILI. DIFMAP is the
software commonly used for the conventional CLEAN method for interferometric
imaging (Shepherd et al., 1995). It interactively establishes a collection of
point source models from the inverse Fourier transform of the visibilities,
i.e., the dirty map. CLEAN windows, which define the regions to search for
CLEAN components, are used during our imaging process. Phase-only self-
calibration is performed after each step of cleaning. Amplitude and phase
self-calibration is performed once a good fit to the visibilities is
established through the multiple steps of cleaning and phase self-calibration.
We repeat the clean and self-calibration loops several times during our
imaging process by gradually decreasing the solution interval of the amplitude
and phase self-calibration. On the other hand, the regularized maximum
likelihood (RML) methods, employed by the eht-imaging (Chael et al., 2016,
2018) and SMILI (Akiyama et al., 2017) libraries, reconstruct images by
minimizing an objective function which is a weighted combination of $\chi^{2}$
of the data and various regularizer terms. The data terms may include the
closure quantities (closure phases and amplitudes; e.g., Thompson et al.,
2017), visibility amplitudes, and complex visibilities. Common regularizers
include the maximum entropy (e.g., Chael et al., 2018), the $\ell_{1}$-norm
(e.g., Honma et al., 2014; Akiyama et al., 2017), the total variation (TV) and
the total squared variation (TSV) of the brightness (e.g., Kuramochi et al.,
2018). With RML methods, it is possible to achieve an angular resolution a few
times finer than the nominal interferometric beam (e.g., Akiyama et al., 2017;
Event Horizon Telescope Collaboration et al., 2019b). During our imaging
process with eht-imaging and SMILI, we started with a Gaussian prior image and
reconstruct images with only the closure quantities, or a combination of
closure quantities and low-weighted visibility amplitudes. After a few
iterations of imaging and self-calibrating, we include full complex
visibilities into the optimization process, further constraining the
reconstructed images. To determine the best set of regularizer combinations,
we survey a range of different weights of each regularizer, in total $\sim$
128 combinations, and select the one that results in the best fit to the
closure quantities.

After imaging of the total intensity, we estimate the instrumental
polarimetric leakage (known as D-terms) for each station using the self-
calibrated data of OJ 287. This process was carried out independently with two
pipelines: the $\mathcal{AIPS}$ task lpcal, and the eht-imaging library, each
based on a particular set of self-calibrated dataset generated during the
total intensity imaging process, i.e., DIFMAP and eht-imaging, respectively.
Both approaches provide consistent values of D-terms. Details of leakage
calibration are described in Appendix A. Polarization imaging of the lpcal
processed data was carried out with DIFMAP. With eht-imaging, the imaging were
performed iteratively with the D-term calculation. Calibration of the absolute
orientation of the EVPAs was performed through comparison with the ALMA array
data (Goddi et al., 2021).

We also carried out non-imaging analysis of the data to measure the properties
of the jet. We perform circular Gaussian model-fitting to the SMILI self-
calibrated visibility data with DIFMAP. The results indicate that the jet
structure can be represented by four Gaussian components. We label the
components following the convention described in Gómez et al. (2022). The
total flux, size, and position offset with respect to the core (the component
at the southeastern end of the jet; see section 3 below) of all components are
listed in Table 1. The uncertainties of the fitted parameters are derived
following the equations outlined in Nair et al. (2019).

Table 1: Model-fitting parameters of OJ 287 with GMVA+ALMA on April 2, 2017.

Comp | S | r | PA | FWHM | $T^{\mathrm{obs}}_{b}$
---|---|---|---|---|---
Name | (Jy) | (µas) | (∘) | (µas) | ($10^{10}$ K)
(1) | (2) | (3) | (4) | (5) | (6)
C0 | 1.25 $\pm$ 0.13 | 0.0 | 0.0 | 31.1 $\pm$ 2.4 | 21.2 $\pm$ 4.1
C1 | 0.96 $\pm$ 0.10 | 63.9 $\pm$ 2.4 | -41.5 $\pm$ 2.2 | 28.9 $\pm$ 2.4 | 18.8 $\pm$ 3.9
C2a | 2.49 $\pm$ 0.25 | 172.9 $\pm$ 2.4 | -40.0 $\pm$ 0.8 | 44.4 $\pm$ 2.5 | 20.6 $\pm$ 3.1
C2b | 0.51 $\pm$ 0.05 | 212.7 $\pm$ 3.5 | -47.7 $\pm$ 0.9 | 42.5 $\pm$ 4.6 | 4.7 $\pm$ 1.2

Note. — Columns from left to right: (1) Component ID, (2) Flux density, (3)
Radial distance from the core component (C0), (4) Position angle, (5)Component
full width at half maximum (FWHM), (6) Observed brightness temperature

##  3 Results

Figure 3:  From left to right: total intensity maps of OJ 287 at 3.5 mm
obtained with GMVA+ALMA observation on April 2, 2017 reconstructed with
DIFMAP, eht-imaging, and SMILI, respectively. The x and y-axis in each image
represent the right ascension and declination axis on the sky, respectively.
The DIFMAP image is convolved with the natural-weighted beamsize of the array,
which is 64 µas$\times$40 µas at a position angle of $-86$ degrees. For the
DIFMAP, eht-imaging, and SMILI images, respectively, the reduced $\chi^{2}$ of
closure phases is: 1.21, 1.22, 1.19; and that of log closure amplitudes is:
1.18, 1.22, 1.08. The second row shows the same images but convolved with a
circular beam of 40 µas. The bottom-right panel shows the model-fitted
circular Gaussian components overlaid on the convolved SMILI total intensity
map. The flux, location, and size of each component are listed in Table 1.

###  3.1 Jet Morphology

Figure 3 shows the total intensity maps of OJ 287 obtained with our GMVA+ALMA
observations, achieving the highest angular resolution to date of the source
at the wavelength of 3.5 mm. The imaging results are consistent across
different imaging methods (CLEAN & RML). Under the nominal resolution, the jet
appears to consist of three major features, extending along the southeast to
northwest direction. We denote the three features as components C0, C1, and
C2, as shown in the bottom right panel of Figure 3.

Component C0, which lies at the southern end of the jet, is compact and shows
the highest brightness temperature (Table 1). This feature is more likely to
be the VLBI core at 3.5 mm. The component C2 has the highest flux density
among the three components. This feature shows complex substructures under the
fine resolution of the RML images (Figure 3 top-middle & top-right). We see
hints of the jet bending and extending towards the western direction
downstream of C2. This bend is more obvious in the lower frequency maps which
are more sensitive to the extended lower brightness regions despite the lower
angular resolutions (e.g., Cohen, 2017; Jorstad et al., 2017). The downstream
jet is largely resolved out and not well-constrained in our high-resolution
images because of their steep spectra and extended structure. Our higher-
resolution images reveal for the first time the twisted morphology of the
innermost, ultra-compact jet region. The first bending occurs between C0 and
C1, with the jet axis gradually changing from north to northwest (clockwise).
We see also hints for a subsequent bending happening downstream of C1 where
the jet axis turns towards the counter-clockwise direction.

The three-component structure is also consistent with the recent 22 GHz
RadioAstron space-ground VLBI observations of OJ 287 made at a similar
resolution (Gómez et al., 2022). However, a position angle difference of
$\sim$ 50∘ of the inner jet can be found when comparing with the RadioAstron
image obtained in 2014. Such a difference could be attributed to the variation
in the position angle in $\sim$ three years. A detailed analysis of the inner
jet position angle variation on a yearly scale and the comparison with
theoretical predictions will be presented in a forthcoming paper (Zhao et al.
in preparation).

In order to quantify the position angle evolution along the jet, we fit the
jet ridge-line on the eht-imaging map. First, we transform the image to polar
coordinates centered at the jet origin and slice it transversely. For each
slice, we store the flux density peak position and then transform them back to
Cartesian coordinates. Thus, we obtain a collection of positions tracing the
jet axis between C0 and C2. The results are presented in Figure 4, where we
also show a sketch to trace the conical structure of C2 in the figure. The jet
axis near C0 extends along a position angle of $\sim$ $-15^{\circ}$, decreases
to $\sim$ $-50^{\circ}$ at C1, and starts to increase again near C2. A similar
trend can be found also in the SMILI image.

Figure 4:  The continuous blue line traces the ridge line of the inner jet of
OJ 287 overplotted on the eht-imaging reconstructed super resolution image.
The dashed blue lines represent the conical structure of the C2 component.

###  3.2 Brightness Temperatures

We investigate the brightness temperature of the OJ 287 jet using two
independent approaches: 1) we calculate the observed brightness temperature of
each Gaussian component from the model-fitting results using the following
equation (e.g., Tingay et al., 2002):

$T_{\mathrm{b}}^{\mathrm{obs}}=1.22\times
10^{12}\frac{S}{\theta_{obs}^{2}\nu^{2}}.$ (1)

where $S$ is the component flux density in Jy, $\theta_{\mathrm{obs}}$ is the
size of the emitting region in mas, and $\nu$ is the observing frequency in
GHz. 2) we calculate the minimum and maximum brightness temperature directly
from the visibilities using the method described in Lobanov (2015). The model
fitting results, which are listed in Table 1, show the observed brightness
temperature of the jet components at 86 GHz ranges from $10^{10}$ to $10^{11}$
K. This is in agreement with the values calculated from the visibility
amplitudes as shown in Figure 5. The brightness temperature values agree
quantitatively with the typical values at the same frequency band (e.g., Lee
et al., 2008; Nair et al., 2019). The 86 GHz brightness temperatures are about
one order of magnitude lower compared to those at 22 GHz obtained from the
RadioAstron results (Gómez et al., 2022). This can be attributed to
differences in intrinsic brightness and opacity between the two frequencies.

We estimate the intrinsic brightness temperature,
$T_{\mathrm{b}}^{\mathrm{int}}$, by (e.g., Gómez et al., 2016):

$T_{\mathrm{b}}^{\mathrm{int}}=\left(1+z\right)\delta^{-1}T_{\mathrm{b}}^{\mathrm{obs}}$
(2)

where $\delta$ stands for the Doppler factor. We adopt the value of the latest
estimates based on the proper motion of moving components by the VLBA-BU-
BLAZAR monitoring program, $\delta=8.6\pm 2.8$ (Weaver et al., 2022). This
gives the intrinsic brightness temperature values
$T_{\mathrm{b,C0}}^{\mathrm{int}}=(3.2\pm 1.7)\times 10^{10}\text{ K,
}T_{\mathrm{b,C1}}^{\mathrm{int}}=(2.8\pm 1.4)\times 10^{10}\text{ K,
}T_{\mathrm{b,C2a}}^{\mathrm{int}}=(3.1\pm 1.4)\times 10^{10}\text{ K, and
}T_{\mathrm{b,C2b}}^{\mathrm{int}}=(0.7\pm 0.4)\times 10^{10}\text{ K}$, for
each component, respectively. These values fall below the equipartition value
of $\sim 5\times 10^{10}$ (Readhead, 1994), indicating possible magnetic
dominance in the innermost jet. However, this is quite uncertain as the errors
in the Doppler factor and brightness temperature values are large.

Figure 5:  Visibility-based brightness temperature estimates of OJ 287 at 86
GHz using the method described in Lobanov (2015). The red and blue dots are
the values of $T_{\mathrm{b,max}}$ and $T_{\mathrm{b,min}}$, respectively. The
orange and purple curves are the rolling mean of the $T_{\mathrm{b,max}}$ and
$T_{\mathrm{b,min}}$ values.

###  3.3 Polarization

We perform polarimetric imaging of the instrumental polarization calibrated
data with CLEAN and eht-imaging independently. The corresponding images are
shown in Figure 6. Our images show that the overall degree of polarization of
OJ 287 is $\sim$ 8 %, which is in quantitative agreement with the ALMA array
results of 8.8 % presented in Goddi et al. (2021). The EVPAs extend mostly
along the mean jet axis, which suggests that the magnetic field in the jet has
a predominant toroidal component. Again, the image reconstructed by eht-
imaging shows fine structure because of the super-resolution that is naturally
achieved by the forward modeling method. However, even in the CLEAN image,
which is convolved with the nominal beam, we see a remarkable polarimetric
structure in the inner jet. The overall structure is consistent between the
two images reconstructed independently with different approaches. We notice
that the apparent difference in the fractional polarization between the two
maps is due to the fact that CLEAN images are convolved with the nominal beam.
The overall degree of polarization and the EVPA distributions agree well
between the two images.

Among the several jet components, C0 shows the lowest fractional polarization
of $\sim$ 5% as measured from the eht-imaging map. This further supports this
component as the jet core, which is usually depolarized (e.g., Lister & Homan,
2005). C1 exhibits a high level of polarization ($\sim$ 16%) which indicate
the magnetic field is more ordered in this region. C2 shows a conspicuous
polarimetric structure which can be further divided into two subcomponents
with the EVPAs lying perpendicularly to each other. The EVPA in the upper
subcomponent also lies perpendicularly to the direction along which the
brightness extends, while in the bottom subcomponent they lie nearly in
parallel. The degree of polarization is $\sim$ 7% and $\sim$ 13% in the upper
and lower sub-component, respectively. These substructures are clearly seen in
both the eht-imaging and CLEAN images.

Figure 6: Polarized images of OJ 287 produced by lpcal +clean method (left)
and the RML imaging method eht-imaging (right). The total intensity image is
shown in a grayscale. The contours represents the linearly polarized flux
density. The ticks show the orientation of the EVPAs where the lengths
indicate the polarization intensity magnitude, and the color represents the
fractional polarization. Only the lpcal+clean image is convolved with the
beam, shown in the bottom-left.

##  4 Discussion

###  4.1 Nature of the C2 Region

Our GMVA+ALMA observations have revealed a remarkable structure of the inner
jet of OJ 287 because of the improved $(u,v)$-coverage and high sensitivity.
In particular, component C2 shows a complex conical structure in both total
and linearly polarized intensity, and a bimodal distribution in the EVPAs.
Previous multi-epoch observations show that this component is nearly
stationary (e.g., Jorstad et al., 2017; Hodgson et al., 2017; Lico et al.,
2022). In the following, we discuss the possible nature of this component.

Oblique shocks could result from the jet striking a cloud of interstellar
media. Under the precessing jet model, this would naturally happen for some
period as the jet sweeps through the ambient material. Since the location of
C2 coincides with where the jet bends, the northeastern section of C2 could be
interpreted as an oblique shock on one side of the jet. The oblique shock is
in a plane making a small angle to the jet boundary on the north-east side.
The flow is then bent by the shock toward the west. The magnetic field could
get compressed to strengthen the component nearly parallel to the jet.
Therefore, the EVPA on the north-east side is roughly perpendicular to the
jet. The southwestern section of C2 could then just be the main jet after the
bend, with the magnetic field transverse to the jet direction at that point,
as usual for a BL Lac object.

Conical shock waves can be formed when there is a pressure imbalance between
the jet plasma and the ambient medium. The properties of shocks in
relativistic jets have been explored by numerical and semi-dynamical
simulations. Gomez et al. (1995) carried out relativistic hydrodynamics (RHD)
simulations of a parsec-scale jet surrounded by ambient medium with constant
or decreasing pressure. The simulations confirmed the existence of stationary
components associated with recollimation shocks. Gómez et al. (1997) simulated
the interaction of standing shocks and relativistically moving perturbations
propagating down the stable jet and found that the shock could enhance the
emission of the moving feature and the stationary component could be
temporarily “dragged” downstream. Further simulations of the interaction
between recollimation shocks and traveling shocks are presented in Fromm et
al. (2016), based on the observations presented in Fromm et al. (2011, 2013a,
2013b, 2015), for the particular case of CTA 102. Various configurations of
the upstream magnetic field components are also included in subsequent
numerical simulations (e.g., Broderick & McKinney, 2010; Porth et al., 2011;
Fuentes et al., 2018). In particular, Mizuno et al. (2015) studied the
kinematically-dominated jets with different magnetic field configurations
including axial, toroidal, and helical based on a relativistic
magnetohydrodynamics (RMHD) simulation code. Fuentes et al. (2018)
characterized the properties of recollimation shocks in RMHD simulations of
jets at the parsec scale as a function of the dominant type of energy:
internal, kinetic, or magnetic. By solving the radiative transfer equations
for synchrothron radiation using as input these simulations, they analyzed the
total intensity and linear polarization signatures imprinted in the stationary
components associated with these shocks. Fuentes et al. (2021) extended the
analysis to RMHD jet models threaded by helical magnetic fields with larger
magnetic pitch angles, and explored as well the effect of different non-
thermal particle populations on the polarimetric properties of stationary
features and the overall observed synchrotron radiation.

On the other hand, Cawthorne & Cobb (1990) established a semi-dynamical model
assuming only the shock front is emitting and found that conical shock waves
could result in polarization angles either parallel or perpendicular to the
jet axis. This model also considered only random magnetic fields in the
upstream jet. In Cawthorne (2006), a poloidal magnetic field component was
added to the model, and the results can explain well the observed polarization
of the knot K1 in 3C 380. Furthermore, Cawthorne et al. (2013) extended this
model to include a paired collimating and decollimating shock and the
predicted EVPA could successfully describe the observational results of the BL
Lac object 1803+784.

Comparing our observational results of C2 with the numerical and semi-
dynamical studies, we find that the conical shape of the emitting region is
quite consistent between our observation and the simulation works. Numerical
simulations predict a series of stationary shocks along the jet that can be
triggered by a pressure imbalance between the jet and the external medium. The
reason we find only one conical-shaped component is most likely the adiabatic
expansion of the jet. As also shown in Gomez et al. (1995), with decreasing
pressure downstream of the jet, the intensity of the stationary components
gradually decreases and the separation between components increases, so the
downstream shocks may be too faint and become undetectable at our observing
frequency. Regarding the polarized emission, the semi-dynamic simulations show
different EVPA distributions across the cone. However, the EVPA pattern is
more symmetric with respect to the cone axis. Numerical simulations also show
that the EVPA pattern will depend on the upstream magnetic field configuration
and the viewing angle (e.g., Mizuno et al., 2015; Gómez et al., 2016; Fuentes
et al., 2021). Fuentes et al. (2021) pointed out that jets with a large
magnetic pitch angle, i.e., threaded by a helical magnetic field dominated by
its toroidal component, can exhibit a bimodal EVPA distribution around
recollimation shocks for small viewing angles. This EVPA configuration could
imply a sign flip of the Stokes $\mathcal{Q}$ parameter that leads to a EVPA
flip, which then results in a dip in the linearly polarized emission, as we
observe in the C2 component from the reconstructed polarimetric images.

Alternative to the standing shock scenario, the observed properties of the C2
component could be a result of geometric effects due to the bending of the jet
axis towards the line of sight. With a decreasing viewing angle, the enhanced
Doppler boosting could amplify the emission in this region and make C2 the
brightest component in the inner jet. If the viewing angle becomes smaller
than the jet opening angle, the bimodal distribution of the EVPAs could be
produced by the existence of helical magnetic fields in the jet as the
direction of the projected magnetic field is different across the component
(Fuentes et al., 2021). This scenario is supported by previous observations
which revealed the existence of a bending around C2 (e.g., Jorstad et al.,
2017; Hodgson et al., 2017; Gómez et al., 2022). However, it is difficult to
explain the conical shape of the emission region with this assumption.

Moreover, by means of multi-epoch GMVA observations, Lico et al. (2022)
identified a new jet feature in the region of C2, in a quasi-concurrent GMVA
observing epoch. The authors argue that the passage of this new jet component
through the stationary feature at 0.1 mas core-separation (i.e., C1) triggered
the high energy outburst during 2016-2017 (Komossa et al., 2017, 2021a)
including the faint VHE flare detected during February 2017 (Mukherjee &
VERITAS Collaboration, 2017) 333In fact, the high X-ray flux detected during
the Swift MOMO program of OJ 287 triggered the VHE observations, which led to
the first VHE detection. and moved down to the C2 jet region at the time of
these observations. In this scenario, the component C2 in our observations
could correspond to the blending of the new feature and the standing shock.
The observed bimodal distribution of the EVPAs could be due to different
polarimetric properties of the two components. A similar case was found in the
core region of PKS 1510-089 during a $\gamma$-ray flare in 2015 (Park et al.,
2019).

###  4.2 Testing the SMBBH model

OJ 287 is one of the most promising candidates to harbor a SMBBH system at the
center. In fact, OJ 287 is among the candidates for hosting a nano-Hz
gravitational wave emitting SMBBH system (Valtonen et al., 2021). The binary
model has been successful in explaining the periodic light curves and
predicting upcoming impact flares, which were confirmed by observations within
a few hours (e.g., Laine et al., 2020). The direction of the jet axis was also
found to be varying with time and this could be also related to the orbital
motion of the BHs (Dey et al., 2021). Models that do not require a secondary
BH to explain the observed variability have also been proposed. For instance,
the flux variation could be explained by viewing angle changes and Doppler
beaming effects of a precessing jet. The precession could be driven by the
Lense-Thirring effect due to the misalignment between the BH spin and the
accretion disc (e.g., Chatterjee et al., 2020; Liska et al., 2018, 2021;
Britzen et al., 2018). MHD instabilities (current-driven or Kelvin-Helmholtz)
would be also possible to produce helical distorted jet structure (e.g.,
Mizuno et al., 2012; Perucho et al., 2012; Vega-García et al., 2019).

Dey et al. (2021) established a model to explain the parsec-scale jet
direction variations at different frequencies in which the jet precession is
powered by the SMBBH with parameters constrained by optical observations. This
model predicts the 86 GHz jet axis should be $\sim$ $-37^{\circ}$ around
April, 2017 assuming a disc model. The position angles of the inner jet
components (e.g., C1, & C2a) measured in our GMVA+ALMA observation agree well
with this prediction (see Table 1). However, we note that this agreement is
partially due to the observing epoch being not far apart from the 86 GHz GMVA
data used to constrain the model. Furthermore, this agreement will not rule
out other possible scenarios. For example, the tilted accretion could also
result in precession of the inner jet. Britzen et al. (2018) argue that the PA
change observed at 15 GHz can be modeled by a jet precession combined with a
nutation of the axis. The precession could be a result of Lense-Thirring
effects and a secondary BH is not always required. Furthermore, our RML images
also revealed a twisted pattern of the innermost jet that resembles a
precessing jet in projection.

Future kinematic studies with multi-epoch GMVA and EHT observations will
hopefully provide further insights to distinguish among different theoretical
models for the underlying nature of the source.

Dey et al. (2021) also explored the possibility of the existence of a jet from
the secondary SMBH based on the SMBBH model. With the high sensitivity and
improved north-south resolution because of the participation of ALMA, we found
no evidence for a secondary jet, even in the eht-imaging and SMILI images with
super resolution. There could be several possible reasons for such a non-
detection. First, the jet is likely to be short-lived, as commented on in Dey
et al. (2021). Since the projected separation of the two SMBHs in April 2017
is $\sim$ 10 µas (Dey et al., 2018), the current image resolution is not
sufficient to spatially resolve the binary system if there is no extended jet
emission from the secondary SMBH. The same would apply if the secondary jet
extends in a similar direction as the primary jet. If the secondary jet is
present and points in a different direction, the non-detection implies that
the brightness temperature of the jet must be lower than $4\times 10^{9}$ K,
which corresponds to three times the r.m.s. level of the eht-imaging map. We
note the dynamic range of our image reconstruction is much higher than the
mass ratio of the two BHs.

We further note that the GMVA+ALMA observations presented in this work are
part of a multiwavelength observing campaign of OJ 287. Close in time
observations with the EHT at 230 GHz (on April 4, & 9, 2017) and with the
RadioAstron space-VLBI mission at 22 GHz (on March 7, 2017) could provide even
higher angular resolutions and probe slightly different regions of the inner
jet. Together with the observations at X-ray and optical bands (e.g., Komossa
et al., 2017, 2020, 2021a, 2021b, 2021c), we will be able to test or obtain
constraints on the physical parameters of the possible jet associated with the
secondary SMBH.

##  5 Summary

We have carried out GMVA+ALMA observations of OJ 287 on April 2, 2017, which
is the first VLBI observation with the phased-ALMA. The improved north-south
resolution and array sensitivity together with the newly developed RML methods
have enabled us to obtain high fidelity, super-resolved images of the OJ 287
jet with unprecedentedly high angular resolution. The convolved RML images
also agree with the CLEAN reconstruction. The images have revealed a twisted
structure in the innermost region of the jet. Our result suggests that the C0
component lying at the southeastern end of the jet is more likely the VLBI
core as it is bright, compact, and relatively depolarized. The component C2
located at $\sim$ 200 µas northwest of the core shows a conical morphology and
complex substructures in polarization. We argue that this component could be
an oblique or recollimation shock, or related to a traveling component passing
through a stationary feature in the jet. We have also carried out the first
attempt to search for a jet from the secondary black hole as proposed by Dey
et al. (2021) based on the SMBBH model. The non-detection could be due to the
small projected separation, the short lifetime, or the difference in the
physical conditions of the secondary jet. The EHT and RadioAstron observations
carried out in 2017 and later could provide further tests of the SMBBH model.

## acknowledgments

The work at the IAA-CSIC is supported in part by the Spanish Ministerio de
Economía y Competitividad (grants AYA2016-80889-P, PID2019-108995GB-C21), the
Consejería de Economía, Conocimiento, Empresas y Universidad of the Junta de
Andalucía (grant P18-FR-1769), the Consejo Superior de Investigaciones
Científicas (grant 2019AEP112), and the State Agency for Research of the
Spanish MCIU through the “Center of Excellence Severo Ochoa” award to the
Instituto de Astrofísica de Andalucía (SEV-2017-0709). This publication
acknowledges the project M2FINDERS that is funded by the European Research
Council (ERC) under the European Union’s Horizon 2020 research and innovation
programme (grant agreement No 101018682). L.L. acknowledges the support of the
DGAPA/PAPIIT grants IN112417 and IN112820, the CONACyT-AEM grant 275201 and
the CONACyT-CF grant 263356. TS was supported by the Academy of Finland
projects 274477, 284495, 312496, and 315721. Y.K. was supported in the
framework of the State project “Science” by the Ministry of Science and Higher
Education of the Russian Federation under the contract 075-15-2020-778. R.-S.
L is supported by the Max Planck Partner Group of the MPG and the CAS, the Key
Program of the National Natural Science Foundation of China (grant No.
11933007), the Key Research Program of Frontier Sciences, CAS (grant No. ZDBS-
LY-SLH011), and the Shanghai Pilot Program for Basic Research – Chinese
Academy of Science, Shanghai Branch (JCYJ-SHFY-2022-013).

This paper makes use of the following ALMA data: ADS/JAO.ALMA#2016.1.01116.V.
ALMA is a partnership of ESO (representing its member states), NSF (USA), and
NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI
(Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA
Observatory is operated by ESO, AUI/NRAO, and NAOJ. The ALMA data required
non-standard processing by the VLBI QA2 team (C. Goddi, I. Martí-Vidal, G. B.
Crew, H. Rottmann, and H. Messias). This work is based partially on
observations with the 100-m telescope of the MPIfR (Max-Planck-Institut für
Radioastronomie) at Effelsberg. This research has made use of data obtained
with the Global Millimeter VLBI Array (GMVA), which consists of telescopes
operated by the MPIfR, IRAM, Onsala, Metsähovi, Yebes, the Korean VLBI
Network, the Greenland Telescope, the Green Bank Observatory (GBT), and the
Very Long Baseline Array (VLBA). The VLBA and the GBT are facilities of the
National Science Foundation operated under a cooperative agreement by
Associated Universities, Inc. The data were correlated at the correlator of
the MPIfR in Bonn, Germany.

ALMA, GMVA, VLBA, GBT, EVN, Effelsberg, Metsähovi Radio Observatory, Onsala
Space Observatory, IRAM-30m RT, Yebes-40m RT

## References

  * Agudo et al. (2012) Agudo, I., Marscher, A. P., Jorstad, S. G., et al. 2012, ApJ, 747, 63
  * Akiyama et al. (2017) Akiyama, K., Ikeda, S., Pleau, M., et al. 2017, AJ, 153, 159
  * Blandford & Payne (1982) Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883
  * Blandford & Znajek (1977) Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433
  * Boccardi et al. (2017) Boccardi, B., Krichbaum, T. P., Ros, E., & Zensus, J. A. 2017, Astron Astrophys Rev, 25, 4
  * Britzen et al. (2018) Britzen, S., Fendt, C., Witzel, G., et al. 2018, MNRAS, 478, 3199
  * Broderick & McKinney (2010) Broderick, A. E., & McKinney, J. C. 2010, ApJ, 725, 750
  * Cawthorne (2006) Cawthorne, T. V. 2006, MNRAS, 367, 851
  * Cawthorne & Cobb (1990) Cawthorne, T. V., & Cobb, W. K. 1990, ApJ, 350, 536
  * Cawthorne et al. (2013) Cawthorne, T. V., Jorstad, S. G., & Marscher, A. P. 2013, ApJ, 772, 14
  * Chael et al. (2018) Chael, A. A., Johnson, M. D., Bouman, K. L., et al. 2018, ApJ, 857, 23
  * Chael et al. (2016) Chael, A. A., Johnson, M. D., Narayan, R., et al. 2016, ApJ, 829, 11
  * Chatterjee et al. (2020) Chatterjee, K., Younsi, Z., Liska, M., et al. 2020, MNRAS, 499, 362
  * Chiaberge et al. (2000) Chiaberge, M., Celotti, A., Capetti, A., & Ghisellini, G. 2000, A&A, 358, 104
  * Cohen (2017) Cohen, M. H. 2017, Galaxies, 5, doi:10.3390/galaxies5010012
  * Deller et al. (2007) Deller, A. T., Tingay, S. J., Bailes, M., & West, C. 2007, Publications of the Astronomical Society of the Pacific, 119, 318
  * Dey et al. (2021) Dey, L., Valtonen, M. J., Gopakumar, A., et al. 2021, MNRAS, 503, 4400
  * Dey et al. (2018) Dey, L., Valtonen, M. J., Gopakumar, A., et al. 2018, ApJ, 866, 11
  * Event Horizon Telescope Collaboration et al. (2019a) Event Horizon Telescope Collaboration, Akiyama, K., Alberdi, A., et al. 2019a, The Astrophysical Journal Letters, 875, L2
  * Event Horizon Telescope Collaboration et al. (2019b) —. 2019b, The Astrophysical Journal Letters, 875, L4
  * Fromm et al. (2016) Fromm, C. M., Perucho, M., Mimica, P., & Ros, E. 2016, A&A, 588, A101
  * Fromm et al. (2015) Fromm, C. M., Perucho, M., Ros, E., Savolainen, T., & Zensus, J. A. 2015, A&A, 576, A43
  * Fromm et al. (2013a) Fromm, C. M., Ros, E., Perucho, M., et al. 2013a, A&A, 557, A105
  * Fromm et al. (2011) Fromm, C. M., Perucho, M., Ros, E., et al. 2011, A&A, 531, A95
  * Fromm et al. (2013b) Fromm, C. M., Ros, E., Perucho, M., et al. 2013b, A&A, 551, A32
  * Fuentes et al. (2018) Fuentes, A., Gómez, J. L., Marti, J. M., & Perucho, M. 2018, ApJ, 860, 121
  * Fuentes et al. (2021) Fuentes, A., Torregrosa, I., Martí, J. M., Gómez, J. L., & Perucho, M. 2021, A&A, 650, A61
  * Goddi et al. (2019) Goddi, C., Martí-Vidal, I., Messias, H., et al. 2019, PASP, 131, 075003
  * Goddi et al. (2021) —. 2021, ApJ, 910, L14
  * Gómez et al. (1997) Gómez, J. L., Marti, J. M., Marscher, A. P., & Ibanez, J. M. 1997, Vistas in Astronomy, 41, 79
  * Gomez et al. (1995) Gomez, J. L., Marti, J. M. A., Marscher, A. P., Ibanez, J. M. A., & Marcaide, J. M. 1995, ApJ, 449, L19
  * Gómez et al. (2016) Gómez, J. L., Lobanov, A. P., Bruni, G., et al. 2016, ApJ, 817, 96
  * Gómez et al. (2022) Gómez, J. L., Traianou, E., Krichbaum, T. P., et al. 2022, ApJ, 924, 122
  * Greisen (2003) Greisen, E. W. 2003, in Astrophysics and Space Science Library, Vol. 285, Information Handling in Astronomy - Historical Vistas, ed. A. Heck, 109
  * Hodgson et al. (2017) Hodgson, J. A., Krichbaum, T. P., Marscher, A. P., et al. 2017, A&A, 597, A80
  * Honma et al. (2014) Honma, M., Akiyama, K., Uemura, M., & Ikeda, S. 2014, PASJ, 66, 82
  * Jorstad et al. (2017) Jorstad, S. G., Marscher, A. P., Morozova, D. A., et al. 2017, The Astrophysical Journal, 846, 98
  * Kettenis et al. (2006) Kettenis, M., van Langevelde, H. J., Reynolds, C., & Cotton, B. 2006, in Astronomical Society of the Pacific Conference Series, Vol. 351, Astronomical Data Analysis Software and Systems XV, ed. C. Gabriel, C. Arviset, D. Ponz, & S. Enrique, 497
  * Komossa et al. (2021a) Komossa, S., Grupe, D., Gallo, L. C., et al. 2021a, ApJ, 923, 51
  * Komossa et al. (2020) Komossa, S., Grupe, D., Parker, M. L., et al. 2020, MNRAS, 498, L35
  * Komossa et al. (2017) Komossa, S., Grupe, D., Schartel, N., et al. 2017, in New Frontiers in Black Hole Astrophysics, ed. A. Gomboc, Vol. 324, 168–171
  * Komossa et al. (2021b) Komossa, S., Grupe, D., Kraus, A., et al. 2021b, Universe, 7, 261
  * Komossa et al. (2021c) Komossa, S., Grupe, D., Parker, M. L., et al. 2021c, MNRAS, 504, 5575
  * Kuramochi et al. (2018) Kuramochi, K., Akiyama, K., Ikeda, S., et al. 2018, ApJ, 858, 56
  * Laine et al. (2020) Laine, S., Dey, L., Valtonen, M., et al. 2020, ApJ, 894, L1
  * Lee et al. (2008) Lee, S.-S., Lobanov, A. P., Krichbaum, T. P., et al. 2008, AJ, 136, 159
  * Lehto & Valtonen (1996) Lehto, H. J., & Valtonen, M. J. 1996, ApJ, 460, 207
  * Lico et al. (2022) Lico, R., Casadio, C., Jorstad, S. G., et al. 2022, A&A, 658, L10
  * Liska et al. (2018) Liska, M., Hesp, C., Tchekhovskoy, A., et al. 2018, MNRAS, 474, L81
  * Liska et al. (2021) Liska, M., Hesp, C., Tchekhovskoy, A., et al. 2021, MNRAS, 507, 983
  * Lister & Homan (2005) Lister, M. L., & Homan, D. C. 2005, AJ, 130, 1389
  * Lobanov (2015) Lobanov, A. 2015, A&A, 574, A84
  * Martí-Vidal et al. (2016) Martí-Vidal, I., Roy, A., Conway, J., & Zensus, A. J. 2016, A&A, 587, A143
  * Mizuno et al. (2015) Mizuno, Y., Gómez, J. L., Nishikawa, K.-I., et al. 2015, ApJ, 809, 38
  * Mizuno et al. (2012) Mizuno, Y., Lyubarsky, Y., Nishikawa, K.-I., & Hardee, P. E. 2012, ApJ, 757, 16
  * Mukherjee & VERITAS Collaboration (2017) Mukherjee, R., & VERITAS Collaboration. 2017, The Astronomer’s Telegram, 10051, 1
  * Nair et al. (2019) Nair, D. G., Lobanov, A. P., Krichbaum, T. P., et al. 2019, A&A, 622, A92
  * Park et al. (2019) Park, J., Lee, S.-S., Kim, J.-Y., et al. 2019, ApJ, 877, 106
  * Perucho et al. (2012) Perucho, M., Kovalev, Y. Y., Lobanov, A. P., Hardee, P. E., & Agudo, I. 2012, ApJ, 749, 55
  * Porth et al. (2011) Porth, O., Fendt, C., Meliani, Z., & Vaidya, B. 2011, ApJ, 737, 42
  * Readhead (1994) Readhead, A. C. S. 1994, ApJ, 426, 51
  * Rioja & Dodson (2011) Rioja, M., & Dodson, R. 2011, AJ, 141, 114
  * Rioja et al. (2017) Rioja, M., Dodson, R., Gómez, J., et al. 2017, Galaxies, 5, 9
  * Satyanarayan et al. (2017) Satyanarayan, A., Moritz, D., Wongsuphasawat, K., & Heer, J. 2017, IEEE transactions on visualization and computer graphics, 23, 341
  * Shepherd et al. (1995) Shepherd, M. C., Pearson, T. J., & Taylor, G. B. 1995, Bulletin of the Astronomical Society, 27, 903
  * Sillanpää et al. (1988) Sillanpää, A., Haarala, S., Valtonen, M. J., Sundelius, B., & Byrd, G. G. 1988, Astrophysical Journal, 325, 628
  * Stickel et al. (1989) Stickel, M., Fried, J. W., & Kuehr, H. 1989, A&AS, 80, 103
  * Tateyama & Kingham (2004) Tateyama, C. E., & Kingham, K. A. 2004, ApJ, 608, 149
  * Thompson et al. (2017) Thompson, A. R., Moran, J. M., & Swenson, George W., J. 2017, Interferometry and Synthesis in Radio Astronomy, 3rd Edition, doi:10.1007/978-3-319-44431-4
  * Tingay et al. (2002) Tingay, S. J., Reynolds, J. E., Tzioumis, A. K., et al. 2002, ApJS, 141, 311
  * Valtonen et al. (2012) Valtonen, M. J., Ciprini, S., & Lehto, H. J. 2012, MNRAS, 427, 77
  * Valtonen et al. (2011) Valtonen, M. J., Lehto, H. J., Takalo, L. O., & Sillanpää, A. 2011, ApJ, 729, 33
  * Valtonen et al. (2006) Valtonen, M. J., Nilsson, K., Sillanpää, A., et al. 2006, ApJ, 643, L9
  * Valtonen et al. (2008) Valtonen, M. J., Lehto, H. J., Nilsson, K., et al. 2008, Nature, 452, 851
  * Valtonen et al. (2016) Valtonen, M. J., Zola, S., Ciprini, S., et al. 2016, ApJ, 819, L37
  * Valtonen et al. (2021) Valtonen, M. J., Dey, L., Gopakumar, A., et al. 2021, Galaxies, 10, 1
  * VanderPlas et al. (2018) VanderPlas, J., Granger, B., Heer, J., et al. 2018, Journal of Open Source Software, 3, 1057
  * Vega-García et al. (2019) Vega-García, L., Perucho, M., & Lobanov, A. P. 2019, A&A, 627, A79
  * Weaver et al. (2022) Weaver, Z. R., Jorstad, S. G., Marscher, A. P., et al. 2022, arXiv e-prints, arXiv:2202.12290
  * Zhao et al. (2018) Zhao, G.-Y., Algaba, J. C., Lee, S. S., et al. 2018, AJ, 155, 26

##  Appendix A Calibration of the Instrumental Polarization

Calibration of the instrumental polarization leakage (also known as the
D-terms) is required to obtain reliable polarimetric maps of the target. Each
of the two pipelines that we used to perform polarimetric imaging (see section
2) has independently implemented this calibration step.

The lpcal pipeline loads the self-calibrated visibility data and the CLEAN
Stokes I image of OJ 287 produced by DIFMAP and runs the $\mathcal{AIPS}$ task
lpcal to solve for the D-terms. lpcal assumes that the source can be divided
into a few sub-components, each with a constant fractional polarization. lpcal
solves the D-terms for each IF independently; the results are shown in the top
left panel of Figure 7. We have flagged the stations that only have data for
one circular polarization (Yebes-40m) and stations that show low S/N on cross-
hands (RL & LR) polarization data (VLBA NL, & ON).

On the other hand, the eht-imaging pipeline performs the instrumental
polarization calibration in parallel with the imaging of the polarimetric data
products. The pipeline computes the leakage terms by minimizing the difference
between the self-calibrated data and the sampled data from the corrupted
reconstructions. For details of the polarimetric imaging with eht-imaging,
refer to Chael et al. (2016). The eht-imaging software by default averages the
data at different IFs, so we have flagged the stations that show large
differences in the D-terms across IFs (VLBA BR, & OV) in our polarimetric
analysis. The eht-imaging results are shown in the top-right panel of Figure
7.

Despite the different approaches for solving the instrumental leakages, the
two pipelines provide very consistent results of D-term estimation, as shown
in the bottom panels of Figure 7 which validates our polarization calibration.
The absolute calibration of the EVPA was obtained by comparison with the ALMA
observations of OJ 287 at the same frequency performed during the same
observation campaign (Goddi et al., 2021).

Figure 7: Top: D-term solutions for each station obtained with lpcal (left)
and eht-imaging (right). lpcal solves for each of the four IFs independently,
while in eht-imaging the data for different IFs are averaged. Bottom:
comparison of the Real (left) and Imaginary (right) components of the D-term
solutions between eht-imaging ($x$-axis) and lpcal ($y$-axis). In all panels,
the circular and square symbols represent the data for RCP and LCP,
respectively. Solutions for different stations are plotted in different
colors. The station name each abbreviation stands for is as follows: AA: ALMA;
EB: Effelsberg (RDBE); EF: Effelsberg (DBBC2); FD: VLBA Fort Davis; KP: VLBA
Kitt Peak; LA: VLBA Los Alamos; PT: VLBA Pie Town; PV: IRAM-30m.