UvA-DARE (Digital Academic Repository) Unraveling the Innermost Jet Structure of OJ 287 with the First GMVA + ALMA Observations

We present the ﬁ rst 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 2017 April 2. The participation of phased ALMA has not only improved the GMVA north – south resolution by a factor of ∼ 3


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 quasiperiodic 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. 2008Valtonen et al. , 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 ∼ 26 µas; e.g., Valtonen et al. 2008), featuring a very massive primary BH of 1.8 × 10 10 M , and a secondary of 1.5 × 10 8 M (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. 2006Valtonen et al. , 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. 2017Komossa et al. , 2021a,b,c),b,c).
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 Ra-dioAstron space-ground VLBI observations.The images revealed a progressive bending of the inner jet with increasing angular resolution by comparison with multiband 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(Komossa et al. , 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.

OBSERVATIONS AND DATA ANALYSIS
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.

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 Hori-  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.
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 stations 1 .Only the common frequency ranges among all participating stations are used in later processing.

Data Reduction
Data correlation was performed with the DiFX correlator (Deller et al. 2007)  The post-correlation dataset was then processed with the ParselTongue (Kettenis et al. 2006) AIPS (Greisen 2003) interface for fringe-fitting and a priori amplitude calibration.We first performed parallactic angle correction with the AIPS task, clcor 2 , 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 subintervals 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 ∼ 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 ∼ 350 at baselines longer than 1.5 Gλ.
A priori amplitude calibration was performed in AIPS with the task apcal by multiplying the system temperatures (T 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 sys measurements were performed using the hot / cold method and therefore already included the opacity correction.The ALMA T sys values have also taken into account the phasing efficiencies derived during the quality assurance and PolConvert pro-1 The recorded bandwidth for each station is as follows: ALMA 32×62.5MHz,VLBA 2×128MHz, most European stations 1×512MHz 2 We 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.
cesses (e.g., Goddi et al. 2019).The cross-hand phase and delay offsets of the reference station were calibrated using the AIPS procedure, vlbacpol.
After the 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.

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;Chael et al. 2018) and SMILI (Akiyama et al. 2017) libraries, reconstruct images by minimizing an objective function which is a weighted combination of χ 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 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 ∼ 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 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 EV-PAs 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).

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   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 ∼ 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 ∼ 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 be-tween 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 ∼ −15 • , decreases to ∼ −50 • at C1, and starts to increase again near C2.A similar trend can be found also in the SMILI image.

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): where S is the component flux density in Jy, θ obs is the size of the emitting region in mas, and ν 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.1.
We estimate the intrinsic brightness temperature, T int b , by (e.g., Gómez et al. 2016): where δ 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, δ = 8.6 ± 2.8 (Weaver et al. 2022).This gives the intrinsic brightness temperature values T int b,C0 = (3.2± 1.7) × 10 10 K, T int b,C1 = (2.8 ± 1.4) × 10 10 K, T int b,C2a = (3.1 ± 1.4) × 10 10 K, and T int b,C2b = (0.7 ± 0.4) × 10 10 K, for each component, respectively.These values fall below the equipartition value of ∼ 5 × 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.

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 ∼ 8 %, which is in quantitative agreement with the ALMA array results of 8.8 % presented in Goddi et al. (2021).The EVPAs  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 b,max and T b,min , respectively.The orange and purple curves are the rolling mean of the T b,max and T b,min values.
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 ∼ 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 (∼ 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 ∼ 7% and ∼ 13% in the upper and lower sub-component, respectively.These substructures are clearly seen in both the eht-imaging and CLEAN images.

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 paral- lel 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 semidynamical 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. (2011Fromm et al. ( , 2013aFromm et al. ( ,b, 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.Fur- 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 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(Komossa et al. , 2021a) ) including the faint VHE flare detected during February 2017 (Mukherjee & VERITAS Collaboration 2017) 3 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 γ-ray flare in 2015 (Park et al. 2019).

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;Liska et al. 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 ∼ −37 • 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 shortlived, as commented on in Dey et al. (2021).Since the projected separation of the two SMBHs in April 2017 is ∼ 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 × 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. 2017Komossa et al. , 2020Komossa et al. , 2021a,b,c,b,c), we will be able to test or obtain constraints on the phys-ical parameters of the possible jet associated with the secondary SMBH.

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 ∼ 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 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 1 .
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.
Figure3shows 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 Figure3.Component C0, which lies at the southern end of the jet, is compact and shows the highest brightness temperature (Table1).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 (Figure3top-middle & top-

Figure 3 .
Figure3.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×40 µas at a position angle of −86 degrees.For the DIFMAP, eht-imaging, and SMILI images, respectively, the reduced χ 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 Table1.

Figure 4 .
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.
Figure 5.Visibility-based brightness temperature estimates of OJ 287 at 86 GHz using the method described inLobanov (2015).The red and blue dots are the values of T b,max and T b,min , respectively.The orange and purple curves are the rolling mean of the T b,max and T b,min values.

Figure 6 .
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.
(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" Zhao et al.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 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).