Probing the innermost regions of AGN jets and their magnetic fields with RadioAstron. V. Space and ground millimeter-VLBI imaging of OJ 287

We present the first polarimetric space VLBI observations of OJ 287, observed with RadioAstron at 22 GHz during a perigee session on 2014 April 4 and five near-in-time snapshots, together with contemporaneous ground VLBI observations at 15, 43, and 86 GHz. Ground-space fringes were obtained up to a projected baseline of 3.9 Earth diameters during the perigee session, and at a record 15.1 Earth diameters during the snapshot sessions, allowing us to image the innermost jet at an angular resolution of $\sim50\mu$as, the highest ever achieved at 22 GHz for OJ 287. Comparison with ground-based VLBI observations reveals a progressive jet bending with increasing angular resolution that agrees with predictions from a supermassive binary black hole model, although other models cannot be ruled out. Spectral analyses suggest that the VLBI core is dominated by the internal energy of the emitting particles during the onset of a multi-wavelength flare, while the parsec-scale jet is consistent with being in equipartition between the particles and magnetic field. Estimated minimum brightness temperatures from the visibility amplitudes show a continued rising trend with projected baseline length up to $10^{13}$ K, reconciled with the inverse Compton limit through Doppler boosting for a jet closely oriented to the line of sight. The observed electric vector position angle suggests that the innermost jet has a predominantly toroidal magnetic field, which together with marginal evidence of a gradient in rotation measure across the jet width indicate that the VLBI core is threaded by a helical magnetic field, in agreement with jet formation models.


INTRODUCTION
The BL Lacerta object OJ 287 (0851+202, z=0.306;Stickel et al. 1989) is considered as one of the most remarkable examples of an active galactic nucleus (AGN). The 12-year quasi-periodic outbursts in the optical regime, often accompanied by multiwavelength flaring activity, are attributed to the existence of a supermassive binary black hole (SMBBH) system, hidden in its compact center (Sillanpää et al. 1988;Valtonen et al. 2016, and references thereafter). The well-aligned with our line of sight radio jet of OJ 287 has shown very remarkable behaviour in the past. During the period 2004-2006, Agudo et al. (2012 reported an erratic jetswing by almost 100 • , attributed to asymmetric accretion flow to the central engine. Systematic jet-axis rotation has been observed by Cohen (2017a) and Britzen et al. (2018), describing a 24−30 years rotation cycle of a helical jet, explained by either the Lense-Thirring effect introduced by the wobbling of the accretion disk around the primary black hole (BH), or precession induced by the binary companion (Thirring 1918;Lense & Thirring 1918).
The elucidation of the origin and the phenomenology of relativistic jets emanating from SMBBH belongs to the frontiers of modern astronomical research. Some theoretical models suggest that relativistic AGN jets can be driven by strong magnetic fields anchored on the black hole's accretion disk (Blandford & Payne 1982). Other scenarios suggest that jets are driven by the conversion of the rotational energy of the BH to Poynting flux via the open magnetic field lines, which are attached to the BH ergosphere (Blandford & Znajek 1977). In reality, superposition of both mechanisms is also possible (e.g., Chiaberge et al. 2000). Close to the central engine, the plasma flow is accelerated and collimated under the presence of a coiled magnetic field. Under such extreme conditions, the dynamics and the stability of the relativistic plasma flow can be strongly influenced by disruptions of the accretion flow, pressure mismatches between the jet and the ambient medium, as well as within the jet itself, leading to the formation of moving shocks, standing shocks, and instabilities (Gómez et al. 1997;Komissarov & Falle 1997;Begelman 1998;Aloy et al. 2003).
Very-long-baseline interferometry (VLBI) is a powerful technique that allows us to achieve the angular resolution needed to resolve and study the fine structure of extragalactic jets. By increasing the observed frequency or the baseline length of the array elements, one can probe jet features located in a region of tens of Schwarzschild radii away from the super-massive black hole (SMBH), and ultimately using millimetre wavelengths (Boccardi et al. 2017) achieve the necessary resolution that has enabled the Event Horizon Telescope Collaboration to obtain the first image of a black hole shadow in M87 (EHTCollaboration et al. 2019a,b,c,d,e,f;EHT Collaboration et al. 2021a,b). Besides VLBI observations in total intensity, polarimetric observations are essential (Gómez et al. 2001). They constitute a powerful tool for deriving fundamental constraints on jet physics and magnetic field configuration.
The images we obtain via VLBI observations allow us to illustrate the relativistic plasma flow and trace distinct features known as blobs or knots. These features usually emerge from the bright and typically unresolved end of the radio jet, known as the core. The nature of this feature differs from that of the real jet base and it is not always easily determined. As also described in Marscher (2008), what we see as the VLBI core can correspond to the surface where the opacity (τ v ) reaches unity (Blandford & Königl 1979), or to a standing shock, which is located downstream of the τ v = 1 surface (Daly & Marscher 1988;Gómez et al. 1995;Mizuno et al. 2015). The true nature of OJ 287 VLBI core is still unexplored.
The RadioAstron space VLBI mission (Kardashev et al. 2013), led by the Russian Astro Space Center and the Lavochkin Scientific and Production Association, operated between 2011 and 2019. It featured a 10 m radio telescope on board of the Spektr-R satellite and was equipped with receivers operating 0.32 GHz (Pband), 1.6 GHz (L-band, dual polarization), 4.8 GHz (Cband) and 22 GHz (K-band, dual polarization). With an apogee of ∼350 000 km, space VLBI observations with RadioAstron are capable of imaging blazar jets in total and linearly polarized intensity with an unprecedented resolution of the order of few tens of microarcseconds (µas) when observing at the shorter wavelengths (e.g., . Three Key Science Programmes (KSPs) on AGN imaging have collected data since 2013 to study the launching, collimation, and magnetic field properties of AGN jets (e.g., Zensus AdSPR on strong sources program), while the AGN survey studied the brightness temperature of their cores (e.g., Kovalev et al. 2016;Kovalev et al. 2020). The RadioAstron Polarization KSP has collected data throughout the whole duration of the space VLBI mission, and is aimed to probe the innermost jet regions and their magnetic field in a sample of the most energetic blazars. Results from the Ra-dioAstron Polarization KSP are reported for 0642+449 in Lobanov et al. (2015), BL Lac in , 0716+714 in Kravchenko et al. (2020), 3C 345 in , and 3C 273 in Bruni et al. (2017) and Bruni et al. (2021). Here we present the first RadioAstron observations of OJ 287 performed in 2014 at 22 GHz, in combination with quasi-simultaneous ground VLBI observations at 15, 22, 43, and 86 GHz.
The structure of this article is as follows. In Sec. 2, we present the multi-frequency data set and the data reduction techniques; in Sec. 3 we report on the results from the VLBI study and discuss their implications in Sec. 4. Throughout this paper we have adopted the following cosmological parameters: Ω M = 0.27, Ω Λ = 0.73, H 0 = 71 km s −1 Mpc −1 (Komatsu et al. 2009), which result for OJ 287 in a luminosity distance of D L = 1.577 Gpc and angular scale of 4.48 pc/mas.

OBSERVATIONS AND DATA ANALYSIS
OJ 287 was observed with RadioAstron in 2014 as part of our Polarization KSP. In this section we describe the RadioAstron observations and data analysis, as well as other close in time ground-based VLBI observations at multiple wavelengths. These include global millimeter VLBI array (GMVA) observations, archival 15 GHz VLBA data, and publicly available 43 GHz data from the VLBA-BU-BLAZAR monitoring program. A summary of the RadioAstron and GMVA observations, as well as the used archival data is provided in Tab. 1.

RadioAstron space VLBI observations at 22 GHz
The bright blazar OJ 287 was observed with RadioAstron on 4-5 April 2014 (from 12:00 to 03:45 UT, observing code GA030E) at a frequency of 22.236 GHz combining a 15:45 hours session during the perigee of the spacecraft (total RadioAstron on-source time was 6.3 hours). These imaging observations were supplemented by 5 short (between half an hour and 2 hours long) sessions on 9, 10, 16, 27 March and 18 April, 2014, within the RadioAstron AGN survey program ).
The RadioAstron perigee imaging session of OJ 287 in 2014 April 4-5 was carried out with a ground array of 12 antennas, as listed in Tab. 1. The long-baseline snapshot sessions were conducted in 2014 March 9 (01:00 to 02:00 UT, at a projected baseline of 15.1 Earth diameters, D ⊕ ), March 10 (01:00 to 02:00 UT, 5.8 D ⊕ ), March 16 (2:15 to 03:00 UT, 15.1 D ⊕ ), March 27 (00:00 to 00:25 UT, 4.6 D ⊕ ), and April 18 (00:00 to 01:00 UT, 9.9 D ⊕ ), all during different orbits of the space radio telescope (SRT) than that of the main perigee imaging session. Fringes to the SRT were found only with the Green Bank telescope for all the snapshot sessions except the one on 2014 March 27, for which fringes were found only to Effelsberg. Signal-to-noise ratio (SNR) values of these long-baseline fringes lie in a range from 10 to 20. This is expected, as baseline sensitivity depends on size of the antennas and baseline projection, resulting in no fringe detection for the smaller antennas in our array in the long-baseline snapshots.
The imaging data were recorded in both, left (LCP) and right (RCP) circular polarizations, with a total bandwidth of 32 MHz per polarization, split into two intermediate frequency (IF) bands. The SRT data were recorded by the RadioAstron satellite tracking station in Pushchino, including extended gaps of approximately one-hour duration to allow for cooling of the on board data downlink radio instrumentation of the Spektr-R satellite. The data of the five snapshot survey sessions were recorded in LCP only.
The imaging data were processed using the RadioAstron-dedicated version of the DiFX software correlator , developed at the Max-Planck-Institut für Radiostronomie. After setting the ground stations clocks, we performed fringe-fitting between the largest antennas of the ground array and RadioAstron, separately for each scan involving it. Such a process allows to minimize the effects of the spacecraft acceleration terms along the considered orbit segment, by recentering the signal in the correlation window every few minutes of observations. For scans giving no fringes (on the largest baselines), we estimated clock values by extrapolating them from successful scans on shorter spacebaselines. This provides a first order value for fringes delay and rate, that can be later refined with the postprocessing data reduction software.
The data of the long-baseline snapshot survey sessions were processed by the software correlator developed at the Astro Space Center of Lebedev Physical Institute in Moscow (Likhachev et al. 2017). The initial data reduction of the correlated data was performed with NRAO's AIPS software package (Greisen 1990). This includes a-priori amplitude calibration using the measured system temperatures for the ground antennas and the SRT. Opacity corrections as a function of source elevation were introduced for the ground telescopes, as well as parallactic angle corrections.
Fringe fitting of the imaging data was performed first by solving for the instrumental phase and single band delays of the SRT on a scan near the perigee of the SRT, providing the best fringe solutions for the ground-space baselines. This was followed by a global fringe search to solve for the residual delays and rates of the ground antennas only. Once the ground array was calibrated, fringe fitting of the SRT was performed by performing a baseline stacking of the ground array, and setting an exhaustive baseline search with the most sensitive ground antennas. Both polarizations and IFs were combined to maximize the sensitivity for ground-space fringe detection. Figure 1 shows the residual delay and rate fringe solutions for the SRT corresponding to the perigee imaging session. The strong ground-space fringes present during the perigee passing of the SRT near the end of the observations allowed for solution intervals as short as 10 seconds, capturing the quickly evolving fringe rates due to the SRT acceleration. Progressively larger solution intervals, up to four minutes, were used to increase the sensitivity on the longer projected baseline lengths to the SRT at the beginning of the experiment. Fringes to the SRT were found throughout the whole duration of the perigee imaging session, up to projected baselines of 3.9 Earth diameters in length. The group delay difference between the two polarizations was corrected using AIPS's task RLDLY, and a final complex bandpass function was solved for the receivers before averaging the fringe-fitted data in frequency across each IF and exporting for subsequent imaging.
Fringe fitting, complex bandpass calibration, and apriori amplitude calibration of the long-baseline snapshot survey data were performed using PIMA software (Petrov et al. 2011) as part of the pipeline data processing of the RadioAstron AGN survey .
The resulting coverage of the fringe-fitted data in the Fourier domain for the perigee imaging session and the long-baseline snapshots is shown in Fig. 2.  Imaging was performed using the DIFMAP software (Shepherd 1997) following the standard CLEAN (Högbom 1974) hybrid imaging and self-calibration procedure. Only data from the perigee session (see Fig. 3) were used for the imaging of the RadioAstron observations. Data from the long-baseline snapshots were excluded based on two arguments. First, no visibility closure quantities that could constrain the imaging and self-calibration procedure were obtained during the long-baseline snapshots, for which only single-baseline fringe detection were obtained. Secondly, the VLBI aperture synthesis technique is based on the assumption that the source remains stationary during the whole integration period being considered. This is no longer valid when considering together the perigee imaging session and the long-baseline snapshots, conducted during different orbits of the SRT (weeks apart), and coincident with a high-activity period in the source, as discussed in Sec. 3.2. This is particularly relevant for the extremely small spatial scales probed during the longbaseline snapshots, with fringe detection spacing between 4.6 D ⊕ and 15.1 D ⊕ . Although the long-baseline snapshot data were not used for the imaging of the Ra-dioAstron data, they provide very valuable information regarding the brightness temperature at the smallest spatial scales, which is analyzed in Sec. 3.5.
Perigee-imaging self-calibrated Stokes I visibility amplitudes and phases as a function of Fourier spacing (uv-distance), and CLEAN model fit to these data are shown in Fig. 3, where we also plot the long-baseline complex visibilities for comparison (see Sec. 3.5 for further discussion). Space VLBI fringes to the SRT extend the projected baseline spacing during the perigee imaging session up to ∼ 3.9D ⊕ , increasing accordingly the angular resolution with respect to that provided by ground-based arrays up to ∼ 56 µas for uniform visibility weighting. RadioAstron images of OJ 287 during our April 2014 observations are shown in the right panels of Fig. 4. In this figure we also show for reference the over-resolved image obtained by down-weighting the short baselines and using a Gaussian beam with full width half maximum equal to the nominal resolution of ∼ 12 µas corresponding to the longest projected baseline detection of ∼15.1 D ⊕ obtained during the long-baseline snapshots, although we stress that these data were not used during imaging process.
Calibration of the instrumental polarization and absolute orientation of the polarization vectors is discussed in Appendix A.1.

Complementary ground VLBI observations at 15, 43 and 86 GHz
We complement the analysis of the space VLBI observations of OJ 287 at 22 GHz with data from three ground-based interferometric observations performed at 15, 43, and 8 GHz within less than two months from the RadioAstron observations (see Tab. 1). At 15 GHz, a single epoch of fully calibrated visibilities is provided via the publicly available database of the long-term monitoring MOJAVE program 1 (Lister et al. 2009) (Principal investigator: J. L. Richards,project code S6407B). The dataset at 43 GHz is part of the VLBA-BU-BLAZAR program 2 , which includes monthly observations of 38 radio and γ-ray bright AGNs. The 43 GHz visibilities data available on-line are fully self calibrated, following the analysis described in Jorstad et al. (2017).
The 86 GHz data were obtained from an observation made on 25 May with the GMVA. The duration of each scan was eight minutes. The data were recorded at 2 Gbps rate (512 MHz bandwidth) with 2 bit digitization, apart from for Plateau de Bure observatory, which recorded at 1 Gbps rate mode and Yebes telescope that was equipped with only LCP receiver. During the correlation procedure, a polyphase filter band technology was used, segmented data into 16 IFs of 32 MHz bandwidth (8 physical IFs) per polarisation. The data were correlated with the DiFX correlator (Deller et al. 2007) at the Max-Planck-Institut für Radioastronomie in Bonn, Germany.
Data reduction of the GMVA observations was performed using NRAO's AIPS software. We first applied the parallactic angle correction, followed by the determination of the inter-band phase and delay offsets between the intermediate frequencies, and the phase alignment across the observing band (known also as manual phase calibration). Next, a global fringe fitting was performed, correcting for the residual delays and phases with respect to a chosen reference antenna. The visibility amplitudes were calibrated, considering the contribution of atmospheric opacity effects based on the measured system temperatures and the gain-elevation curves of each telescope. Finally, the a-priori calibrated data were exported to DIFMAP for subsequent imaging following the usual CLEAN and self-calibration procedure. The absolute EVPA calibration, as well as the leakage calibration method, is described in Appendix A.2. Note-a : All images have been obtained using uniform weighting except for this one, for which natural weighting was used. b : Over-resolved image of OJ 287. Columns from right to left: (1) Observing frequency, (2) Peak flux density, (3) Total flux density, (4) Noise level, (5) Peak of linearly polarized flux density, (6) Total linearly polarized flux density, (7) Noise level in the polarization image.

Visibility model fitting
We modeled the observed brightness distribution in the jet of OJ 287 with two-dimensional Gaussian components in the visibility plane using the DIFMAP package and following the procedure described in  RadioAstron 22 GHz polarimetric space VLBI images obtained with natural and uniform weighting in April 2014. The overresolved total intensity RadioAstron image at the nominal resolution corresponding to the maximum projected baseline fringe detection during the long-baseline snapshots is also shown for reference. Total intensity contours, EVPAs, and model-fit components are over-plotted. The lowest contours are ±10 and ±9 times the rms level (see Tab. 2) for the 86 GHz and RadioAstron over-resolved images, respectively, and ±7 for the remaining images, with successive contours in factors of 3 up to 90% of the total flux density peak. Synthesized beams and their parameters are located in the bottom left corner of each map.    Figure 4 shows the ground array (left panels) and Ra-dioAstron space VLBI images (right panels) of OJ 287 taken between April and May 2014, providing a detailed view of the jet at different spatial scales.
The 15 GHz image (left upper panel in Fig. 4) shows a bright VLBI core and a one-sided jet that extends towards the west, followed by a strong jet bending by about 55 • towards the south-west. The core region is modeled by a single Gaussian component labeled as C (upper left frame in Fig. 4), considered as the positional reference for the remaining modelfit components. At higher frequencies, this region is resolved into two subcomponents C1 and C2 of the core and an extra jet component J6 (see Fig. 4).
At 22, 43 and 86 GHz images the modelfit component C1 is considered as the reference point. In the overresolved RadioAstron image, the knot C1 can be further resolved into two sub-components C1a and C1b (bottom right frame in Fig. 4). We note that the sum of the flux densities of C1a and C1b components approximately adds to that of component C1 (see Tab. 3).
The progressively higher angular resolution obtained with the VLBA-BU-BLAZAR 43 GHz image (middle left panel in Fig. 4), and GMVA 86 GHz image (bottom left panel) allows us to map the innermost regions of the VLBI core area up to an angular resolution of about 40 µas. This reveals that the jet bending observed at 15 GHz continues as we probe deeper into the VLBI core, from the west jet direction observed at 15 GHz to the north-west jet direction visible at 86 GHz. The jet curvature at these spatial scales has also been reported by Hodgson et al. (2017), and is clearly imprinted here as well in the position angle of the model fitted components listed in Tab Our space VLBI images of OJ 287 (right panel in Fig. 4) confirm that this progressive jet bending with increasing angular resolution continues up to the smallest scales probed by RadioAstron. At these extremely high angular resolutions the core area can be resolved into two distinct components C1a and C1b. Assuming that the upstream end of the jet characterized by component C1a corresponds to the VLBI core, the innermost jet depicted by the RadioAstron images shows a counter-clock wise rotation from −11.7 • ± 1.0 • for C2 to 13.1 • ± 0.4 • for component C1b, as we probe deeper into the upstream jet.
Downstream from the core region, the jet structure is well represented by up to six Gaussian components, labeled J1 through J6 (see Fig. 4). Except for the outermost component J1, all other jet components are crossidentified at multiple frequencies.
We note here that one can try to decompose the 15 GHz flux density of the component C (3.1 Jy) into contributions from C1 and C2 by calculating their spectral indices between 22 GHz and 43 GHz and using them to estimate their respective flux densities at 15 GHz. This procedure yields the following estimates: S 15GHz,C1 = 0.2 Jy and S 15GHz,C2 ≈ 1.0 Jy. Taken together, they sum up to 1.2 Jy, which is about 1.9 Jy less than the flux density actually measured at 15 GHz for the component C. This discrepancy most likely results from the fact that the 22 GHz RadioAstron observations were made 29 days before the 43 GHz observation and 31 days before (see Fig. 5 and Sec. 3.3), where OJ 287 was in a lower flux density state. As OJ 287 was undergoing the rising stage of a flare during this period (see Sec. 3.2), the discrepancy of 1.9 Jy between the estimated and measured flux density at 15 GHz suggests that the 22 GHz flux density of one or both components (C1 and C2) blended at 15 GHz into the single component C has increased during this period.
The fitted flux densities and sizes of the modelfit components were used to calculate also the brightness temperature of each VLBI knot in the rest frame of the source (e.g., Pushkarev & Kovalev 2012): where S is the component flux density in Jy, θ obs is the size of the emitting region in mas, ν is the observing frequency in GHz, and z is the source redshift. For unresolved modelfit components, we set θ = θ min , where θ min is the resolution limit (Lobanov 2005) obtained for the respective component, and consider the resulting estimate of T b as a lower limit. The brightness temperatures estimated using this procedure are given in column 9 of Tab. 3.

Polarization
Linear polarization is detected in all of our images in Fig. 4, with the exception of the over-resolved RadioAstron image. Tab. 4 summarizes the linear polarization properties of the model fitted components, whereas in the Appendix A.3 we present the polarization error estimation method that we followed. As it is commonly found in VLBI imaging, there is not a one-to-one correspondence between the components seen in total intensity and linear polarization. Hence, while the model fitting for total intensity is performed in the visibility plane, the two-dimensional distribution of the fractional polarization and the electric vector position angle (EVPA) for each knot are calculated in the image domain, based on the amount of polarized emission inside the area of each component. The polarization structure at 15 GHz consists of a low polarized (m ∼ 1.2 %) core with EVPAs in the direction of the innermost jet, and polarized emission in the extended jet between model fitted components J4 and J3. Our characterization of the polarized emission in the jet as components P J4 and P J5 (see Tab. 4) shows an increased degree of polarization with respect to that of the core from 3 % to 7.5%, as expected in the case of lower opacity, and polarization vectors aligned with the local jet direction. As discussed below in Sec. 3.6, the observed EVPAs are not severely affected by Faraday rotation effects, from which we can conclude that the VLBI core and jet in OJ 287 are characterized by a dominant toroidal magnetic field component.
The linearly polarized emission at 43 and 86 GHz shows a progressive increase in degree of polarization with observing frequency in the core, reaching m ∼ 21% at 86 GHz, which is consistent with a transition from an optically thick core at 15 GHz to optically thin at 86 GHz, as discussed in more detail in Sec. 3.3 below, and may be also affected by beam depolarization (Burn 1966). The emission downstream of the VLBI core area at 43 and 86 GHz, characterized as components P J5 and P J6 , shows EVPAs predominantly oriented perpendicular to the local jet direction, and a higher degree of polarization than what is observed further downstream in the jet area corresponding to components P J4 and P J3 at 15 GHz. This suggests that the jet area of P J5 , and maybe the region comprised between components P J5 and P J6 as well, may correspond to a recollimation shock that compresses the magnetic field in the direction of the jet, explaining the different EVPAs. However, we should also note that according to numerical simulations (e.g., Fuentes et al. 2018;Fuentes et al. 2021;Moya-Torregrosa et al. 2021) the observed net polarization in recollimation shocks can change significantly depending on the underlying magnetic field configuration, viewing angle, and other jet parameters. Alternatively, as mentioned in Sasada et al. (2018), the observed polarization may be originated by an oblique shock located where the jet bends to the west.
Alternatively, as suggested by Fuentes et al. (2018);Fuentes et al. (2021), EVPA depends on the helical magnetic field pitch angle. Lastly, is it worth mentioning that the jet presents a bend in that location, so it could correspond to an oblique shock (Sasada et al. 2018).
A higher degree of magnetic field ordering is expected in the shocked plasma, which would also explain the higher degree of polarization in this region (Jorstad et al. 2007;Hovatta et al. 2016).
The space VLBI RadioAstron images from the perigee imaging (top and middle right panels in Fig. 4) show low degrees of linear polarization in the core area corresponding to components C1 and C2, with values of the order of 1% (see Tab. 4), and EVPAs aligned with the local direction of the jet. These low degrees of polarization are expected for optically thick emission, in agreement with what was observed at 15 GHz. The EVPAs are also in concordance with a predominant toroidal magnetic field, or a helical magnetic field with a large pitch angle. The natural weighted image (top right panel in Fig. 4) shows weakly polarized emission in the area associated with component J6, with EVPAs perpendicular to the local jet direction, providing further support for the possibility that this corresponds to a recollimation shock. The measured degree of polarization for component P J6 is smaller than that observed at lower frequencies, which suggests that this component may be also affected by opacity effects at 22 GHz. Polarization information is resolved out in the over-resolved RadioAstron image (bottom right panel in Fig. 4).

Flaring activity of OJ 287 during 2014
OJ 287 showed multi-wavelength flaring activity during April-May 2014, coincident with our VLBI imaging campaign, as shown in the light curves of Fig. 5. The event started at high energies, with γ-ray emission rising on the 25th of February 2014 (shadowed area in Fig. 5). For the complete γ-ray light curve see the VLBA-BU-BLAZAR website 3 . A local optical V band maximum was reported on March 31st (Ganesh et al. 2014), followed few days later by a flux density increase by almost ∼270% from the flux density level recorded between January to February, as measured by singledish observations at 21, 23, 37, and 86 GHz. According to Weaver et al. 2021, (submitted) this flare is associated with the appearance of new features in the OJ 287 jet. The quasi-simultaneous to the VLBI observations light curves were provided by the 40-m telescope of the Owens Valley Radio Observatory (OVRO) at 15 GHz (Richards et al. 2011), the RATAN-600 radio telescope at 5, 8, 11, and 22 GHz (Kovalev et al. 1999(Kovalev et al. , 2002, the F-Gamma multifrequency monitoring program (Angelakis et al. 2019) at 23 GHz, the monitoring project of extragalactic radio sources by Metsähovi Radio Telescope (Teraesranta et al. 1998) at 37 GHz, and the PO-LAMI monitoring program 4 at 86 GHz (Agudo et al. 2018b,a).
We note that the non-simultaneity of VLBI data acquisition under such flaring conditions can introduce a big uncertainty to the spectral index estimation and the Faraday-rotation analysis results, discussed in Sec. 3.3 and Sec. 3.6.

Spectral analysis
The quasi-simultaneous multi-frequency observations in April-May 2014 enable us to analyze the spectral properties of the OJ 287 jet from parsec to sub-parsec scales during its flaring activity. For this analysis, we consider the spectrum, S(ν), of synchrotron selfabsorbed (SSA) emission from a homogeneous spherical region filled with relativistic electron-positron plasma 3 http : //www.bu.edu/blazars/VLBA GLAST/oj287.html 4 www.polami.iaa.es/ with a power-law energy distribution of emitting particles (Türler et al. 2000): where S ν is the observed flux density in Jy, ν m is the turnover frequency in GHz, S m is the turnover flux density in Jy, τ m ∼ 3/2 (1 − (8α/3α t )) 1/2 − 1 is the optical depth at ν m , and α t and α are the spectral indices of the optically thick and thin parts of the spectrum, respectively (using the S ν ∝ ν α definition of spectral index). In our data, crude estimates of the synchrotron turnover point can be made for four components (C1, C2, J6, and J5) by using the Equation 2 and making an assumption about one of the two spectral slopes. To decide on the most plausible distribution of the estimated 1.8 Jy increase of S var,22 , we inspect the component spectra obtained using the actually measured flux densities (see Fig. 6). This inspection suggests that the 22 GHz flux densities of the components C1 is the most likely to be affected by the variability between the epochs t 1 and t 2 . Its flux density at the epoch t 2 should have then increased to ≈ 2.3 Jy. If we use this value to repeat the 15 GHz flux density decomposition of the core component C discussed in Sec. 3.1.1, we obtain S 15GHz,C1 ≈ 1.9 Jy. The respective estimated flux density of the component C then becomes ≈ 3.1 Jy, which is in excellent agreement with the actually measured flux density of this component. This agreement further supports the suggestion that the variability observed at 22 GHz between the epochs t 1 and t 2 can be ascribed to flux density changes of the component C1 in the core region. We therefore apply this assumption for the spectral fitting described below.   For the component C2, the observed spectrum does not warrant estimating the turnover point, and crude limits on S m and ν m can be provided by the component flux density measured at 86 GHz. If we assume that the turnover flux density of C2 is similar to that of C1, then the turnover frequency of C2 should be ≈ 115 GHz. For the component J6, the measured flux densities can be used for estimating S m and ν m only if an assumption is made about the spectral index α of the optically thin part of the spectrum. We assume α = −0.7 for this component. For the component J5, the 15 GHz flux density is likely to be affected by blending with the component J6 upstream. We estimate this blend to contribute ≈ 0.2 Jy to the flux density of J5, and we correct for this blend before fitting the spectrum. The fit is then done with an assumption of α = −0.8 derived from the spectral index measured for J5 between 43 and 86 GHz. The resulting overall constraints for the component spectra and the best fit models are presented in Tab. 5. The fits are also plotted in Fig. 6.
Results of the spectral fitting indicate an unusual spectral evolution, with the turnover frequency first rapidly rising downstream from the component C1 (with ν m > 86 GHz in the component C2) and then falling back. Such a behavior can be explained by relativistic shocks undergoing a transition from Compton-to synchrotrondominated emission regime (Marscher 1990), although reaching viable conclusions on this matter requires taking into account the acceleration of the emitting plasma (Lobanov & Zensus 1999).

Magnetic fields and equipartition Doppler factors
The estimated turnover frequencies and flux densities of the jet components can be used for calculating their respective magnetic field strengths and Doppler boosting factors. The magnetic field strength of a spherical emitting region with an SSA spectrum is given by (Marscher 1983): where b (α) is a coefficient depending on the spectral index (with b(α) = 3.2 for α = −0.5; see Table 1 in Marscher 1983 and Appendix A in Pushkarev et al. 2019), ν m is the spectral turnover frequency in GHz, S m is the spectral turnover flux density in Jy, θ m is the diameter of the emitting region in mas, and δ j is the Doppler boosting factor. The turnover parameters, S m and ν m , are taken from the spectral fits in Tab. 5. The diameters of the emitting regions should be obtained from measurements made at the turnover frequency. In lieu of such measurements, we estimate these diameters by requiring them to reconcile the fitted turnover parameters with the maximum value of brightness temperature listed for a given feature in Tab. 3. For this purpose, we first use Eq. 1 which approximates the brightness distribution by a two-dimensional Gaussian and provides its FWHM. To obtain an estimate of θ m , the FWHM values are further multiplied by a factor of √ 3/ √ 2 ln 2 accounting for the conversion from the Gaussian to the spherical shape. The conversion factor is determined by calculating the diameter of a sphere filled with homogeneous, optically thin plasma such that it provides the same total and peak flux density as those derived for a given two-dimensional Gaussian component of the modelfit. The resulting estimates of θ m and B ssa are listed in the first two columns of Tab. 6. Independent estimates of magnetic field strength can be obtained from the equipartition condition (Pacholczyk 1970;Zdziarski 2014), which yields, with the same units as in Eq. 3: where D L,Gpc is the luminosity distance in Gpc, f is volume filling factor of the emitting plasma (with f =1 adopted hereafter), k u is the ratio of total energy in the emitting region to that carried by the SSA population of electrons (with k u ∼ 1 representing an electronpositron flow), the coefficient c 12 is given in Pacholczyk (1970), and κ ν depends on the low and high frequency cutoffs, ν min and ν max , of the SSA emission. For ν min = 10 7 Hz and ν max = 10 13 Hz assumed in this paper, the ν min contribution to κ ν can be neglected, and κ ν ≈ (ν max /ν m ) 1+α /(1 + α). The resulting estimates of B eq are shown in Tab. 6. Both, B ssa and B eq depend on the Doppler boosting factor, and this dependence can be used for deriving the equipartition Doppler boosting factor, δ eq , by requiring that B ssa = B eq . This condition yields δ eq ∝ (B SSA /B eq ) 7/15 . Using this approach, we calculate δ eq and present them in Tab. 6 for electron-positron (k u = 1) and electron-proton (k u = 100; see Merten et al. 2017) jet.
The observed proper motions of the jet components in OJ 287 at 15 GHz (Lister et al. 2009;Lister et al. 2016) correspond to apparent speeds β app ≤ 15 c, with a median speed of 4.5 c, while at 43 GHz VLBA-BU-BLAZAR monitoring have estimated values between ∼12c and ∼7c (Weaver et al. 2021, submitted). One should therefore expect to find Doppler factors δ j ≥ 1 + β 2 app ≈ 5 − 15. With this estimate, it is feasible to conclude that plasma condition in the features J6 and J5 are likely to be close to the equipartition. The components C1 and C2 in the core region deviate from the equipartition, with a good indication for this deviation to progressively increase at smaller separations from the origin of the jet.

Brightness temperature
The highest values of T b estimated for C1 reach ≈ 5 × 10 12 K (see Tab. 3). This is moderately larger than the inverse-Compton limits of (0.3 − −1.0) ×10 12 K (Kellermann & Pauliny-Toth 1969) and substantially above the equipartition limit of ≈ 5 × 10 10 (Readhead 1994). At the inverse-Compton limit, a jet Doppler boosting factors δ j 12 would be required to explain these brightness temperature estimates. Thus, the kinematic constraints δ j ≈ (5 − 15) discussed Sec. 3.4 can reconcile the highest estimated brightness temperature values (see Tab. 3) with the inverse-Compton limit.
The 22 GHz modelfit estimates of T b can be compared with the estimates of minimum brightness temperature, T b,min made directly from the visibility amplitudes (Lobanov 2015). This comparison is presented in Note-Column designations: θm -estimated effective diameter of the emitting region; Bssa -magnetic field strength for synchrotron self-absorbed spectrum; Beqstrength of the equipartition magnetic field; δeqequipartition Doppler factors calculated for electron-positron (e − e + , ku = 1) and electron-proton (e − p, ku = 100) flow, with δeq ∝ k 2/15 u . For the component C2, the second row lists respective values obtained for the assumed spectral turnover point (see Sec. 3.4). Fig. 7, showing that T b,min reaches ≈ 10 13 K at the longest (u, u)-spacings coming from the auxiliary long-baseline segments of the observations.

Fig. 2 and
In Fig. 7, T b,min shows a continued increasing trend with progressive larger uv-spacing, with no evidence of reaching a plateau. Thus, even larger T b,min could in principle be expected at longer baselines than those observed here.  report T b,min values of the order of 10 13 K at the longest baselines during Ra-dioAstron observations of the jet in BL Lac at 22 GHz. The authors interpret this extremely high T b,min values as resulting from a flaring event that was taking place during the RadioAstron observations, causing the jet flow to depart from equipartition. Similarly, our OJ 287 observations were performed during the onset of a dramatic flaring event in the source (see Sec. 3.2). We should also note that the long baseline snapshots that provide the largest T b,min were obtained during different orbits of the SRT, weeks before and after the perigee imaging session, probing therefore probably different flaring states of the source. This provides a natural explanation for why in Fig. 3 the visibilities obtained during the long-baseline snapshots do not provide a good fit to the CLEAN model obtained using only the perigee imaging session.
The measured T b,min ∼ 10 13 K require Doppler boosting factors, δ j , of the order of 10 -30 to reconcile them with the inverse-Compton limit. This is broadly in agreement with the values estimated from the proper motion of components moving downstream the jet by the MOJAVE (Lister et al. 2016) and VLBA-BU-BLAZAR (Weaver et al., 2021, submitted) monitoring programs of δ j ≈ 5 − 15 and δ j = 8.6 ± 2.8, respectively. At the upper boundary of this range, the respective viewing angle of the jet should be θ j ≈ 3 • − 8 • . Reducing this angle by about 2 • would lead to increasing the Doppler boosting factor up to ∼ 30 even without requiring the bulk Lorentz factor to increase.

Faraday rotation
According to theoretical models, jet launching from super-massive black holes is electromagnetic in nature through the action of helical magnetic fields, either anchored in the accretion disk or the black hole ergosphere (Blandford & Znajek 1977;Blandford & Payne 1982). The observational signature of these fields is imprinted in the polarization information of the synchrotron radiation, and is tightly connected to the phenomenon of Faraday rotation: the rotation of the polarization plane that occurs when a polarized electromagnetic wave passes through a magnetized plasma (the so-called Faraday screen). The rotation of the polarization angle, χ, introduced by the Faraday screen is determined as χ = χ o + RM × λ 2 , where λ is the observed wavelength and χ 0 is the intrinsic electric vector position angle of the emitting region. The rotation measure (RM) is expressed by (e.g., Thompson et al. 2017): where RM is measured in rad/cm 2 , 0 is the permittivity of the vacuum, ν e is the electron density, B is the component of the magnetic field which is parallel to our line of sight, and dl is the path length from the source to the observer through the de-polarizing plasma. The sign of the Faraday rotation coincides with the sign of the line-of-sight magnetic field. In this work, we have combined our highest resolution polarimetric images obtained with the ground array at 43 and 86 GHz, with that of RadioAstron at 22 GHz (natural weighted image) to produce the RM map of the innermost jet regions in OJ 287. For this we first convolved the images at the three different frequencies with a common restoring beam of 0.2 × 0.1 mas at position angle of 0 • , which slightly over-resolves the 7 mm image while still preserving a fraction of the higher resolution achieved at 86 GHz and 22 GHz. Alignment of the images at the three different frequencies was obtained by performing a cross-correlation of the total intensity images following the approach described in , and references therein. The estimated shifts were all of the order of the image pixel size of 6µas, and therefore no corrections were applied.
The obtained RM map, with overlaid Faradaycorrected EVPAs is shown in Fig. 8. It should be noted that the images used to compute the RM were obtained at different epochs (see Tab. 1) during a flaring state of the source (see Sec. 3.2), and therefore the reliability of the obtained RM map relies on the assumption of negligible polarization structural changes in the time span covering the considered observations. Figure 8 reveals a region with enhanced RM in the VLBI core area of OJ 287, with a median RM of −1000±300 rad/m 2 . This is broadly consistent (in magnitude) with a previous core RM estimation of −367±71 rad/m 2 based on simultaneous VLBI observations between 8 GHz and 15 GHz taken in 2006 April 28 .
The size of the RM region is of the order of our resolution beam, and therefore we lack the necessary angular resolution to accurately measure gradients in the RM (e.g., Taylor et al. 2010). With this limitation in mind, we have plotted in the right panel of Fig. 8 several cuts of the RM perpendicular and parallel to the local jet direction, which are indicative for the presence of gradients across and along the jet. Gradients in RM along the jet direction are expected to be associated with a progressive decrease in the magnetic field strength and electron energy density with distance along the jet (Jorstad et al. 2007). On the other hand, gradients across the jet width are indicative for the presence of a toroidal magnetic field component. This is also consistent with the measured Faraday-corrected EVPAs, which are aligned with the local jet direction, suggesting also the presence of a predominant toroidal component in the magnetic field. The negative gradient in RM from north-east to south-west direction suggests a toroidal magnetic field oriented clockwise as seen in the direction of flow motion.
The indicative RM gradient across the jet width and EVPAs are both in agreement with the presence of a helical magnetic field threading the innermost regions of the jet in OJ 287, as predicted by jet formation models (Blandford 1993;Blandford & Payne 1982) and observed previously in a number of sources (e.g., Gabuzda et al. 2004;O'Sullivan & Gabuzda 2009;Gabuzda 2017).

DISCUSSION AND CONCLUSIONS
OJ 287 is considered as one of the best candidates for harboring a binary black hole (BBH) system in its center. The updated BBH central engine description for OJ 287 allows us to track the changes in the orientation of the accretion disk and the primary BH spin (Valtonen & Pihajoki 2013). However, additional assumptions are needed to explain its decades long radio jet observations. Interestingly, the assumption that the jet is launched perpendicular to the innermost disk axis leads to BBH orbital time scale observational implications (Valtonen et al. 2012;Dey et al. 2021). Further, the time evolution of the innermost disk axis shows up as a bending of the radio jet. This is because the changes at the launch angle are propagated to the more distant parts of the jet with a time delay. The rather small variations of the launch angle are further amplified by projection ef-fects since the viewing angle of the jet is small. This is in agreement with our observations that confirm a progressive jet bending with increasing angular resolution up to the smallest spatial scales probed by RadioAstron.
Employing VLBI data sets at 15, 43 and 86 GHz, a consistent description for the temporal evolution of OJ 287 radio jet was provided in Dey et al. (2021) making use of a helicity parameter that allows for outward jet motion that is not exactly in a straight line, as noted in Valtonen & Pihajoki (2013). In addition, one may use the information of the time evolution of optical polarization and determine the jet orientation close to the jet launch site (Valtonen & Wiik 2012;Sasada et al. 2018).
In the jet distance range from 0.2 mas to 1 mas the component position angles (PA) listed in Tab. 3 are in agreement with Dey et al. (2021) within the errors. For example, the PA of component J6 is −55 ± 10 • in the model while here we find −36 ± 13 • (86 GHz). The optical polarization data arise from a knot at unknown distance, but judging from the agreement of the PA of our innermost knot C2 of the 86 GHz map and the PA of the model (Valtonen & Wiik 2012), it is possible that the optical emission region is not far from knot C2. Beyond 1 mas our PA values agree with the earlier model of Valtonen & Wiik (2012).
From our BBH central engine prescription, we expect the inner components C2 and J6 to rotate counter clockwise (i.e., increasing values of PA) by about 15 • degrees between 2014 and 2017 while the changes further out in the jet would be barely noticeable.
Even though the observed periodicity in the light curve of OJ 287 and innermost position angle changes can be explained by a BBH model, alternative physical models can also lead to a similar phenomenology. Britzen et al. (2018), based on 22 years observations at 15 GHz found that the jet of OJ 287 is rotating with a period of ∼30 yr. Modeling of OJ 287 radio data showed that this rotation can be explained by a combination of jet precession and nutation. The physical cause of the precession can be driven by a binary system in OJ 287 center, as well as by the mechanism of the Lense-Thirring precession by the tilted accretion disk of a single BH (e.g., Liska et al. 2018), which in the case of OJ 287 provides realistic parameters. Agudo et al. (2012) propose that variable asymmetric injection of the jet flow, perhaps related to turbulence in the accretion disc, coupled with hydrodynamic instabilities leads to the non-ballistic dynamics that causes the observed non-periodic changes in the direction of the inner jet. Cohen (2017b) presented evidence that OJ 287 is behaving as a rotating helix based on the study of MOJAVE 15 GHz VLBA images from 1995 to 2015. The results of the ridge line analysis of the data showed that the jet is rotating with a period of possibly ∼30 yr. The inner jet apparently seems to have moved to a new direction after the rotation, indicating that jet nozzle has been re-oriented. A model of a helical jet, observed from a small and varying viewing angle, had been proposed earlier by Valtonen & Pihajoki (2013). Another suggested scenario by Hodgson et al. (2017) proposes that changes in the position angle of the jet are due to opacity shifts of the observed core in a bent jet.
On the other hand, the 12-year periodicity of OJ 287 optical light-curves, that usually is attributed to a binary black hole system in the center of OJ 287 can be also explained by the existence of a non-radially moving feature along a helical jet. Butuzova & Pushkarev (2020), showed that the differences between the peak flux values of the periodic optical flares, as well as the time-lag between optical and radio repeated variability, can be caused by the development of helical mode of the Kelvin-Helmholtz instability, inside a, well-aligned with our line of sight, helical jet.
Our polarization images are consistent with the innermost jet in OJ 287 being threaded by a predominantly toroidal magnetic field, while our Faraday rotation map is indicative of the presence of a gradient in rotation measure across and along the jet. Both pieces of information suggest that the innermost jet region in OJ 287 is threaded by a helical magnetic field, as expected from jet formation models, and in agreement with previous studies. In particular Myserlis et al. (2018) report a clockwise EVPA rotation by ∼ 340 • during a multifrequency single-dish campaign taken with the 100 m Effelsberg radio telescope in 2016, which is interpreted as produced by a polarized component propagating on bent jet threaded by a helical magnetic field.
The spectral analysis combining our RadioAstron observations and ground-based VLBI observations is in agreement with the parsec-scale jet being in equipartition between the particles and magnetic field. However, there are some clear evidence for the jet being dominated by the internal energy of the emitting particles as we probe progressively closer to the central engine in the VLBI core area, in agreement with the onset of a large multi-wavelength flare that peaked few months after our VLBI observations. Brightness temperatures have been estimated from the model fitted components, as well as from the visibility amplitudes. The maximum observed brightness temperature of 5.2×10 12 K for the VLBI core can be reconciled with the inverse-Compton limit assuming a moderate Doppler boosting factor of the order of 5-15, in agreement with those estimated from the proper motion of superluminal components. Ground-space fringes have been detected up to a record projected baseline distance of 15.1 Earth diameters in length (one of the longest ever obtained with RadioAstron at 22 GHz), from which we have estimated a minimum brightness temperature of T b,min ∼ 10 13 K. The rising T b,min trend with projected baseline length with no indication of reaching a plateau suggests that even larger brightness temperatures could be measured with longest baselines, although such extremely high brightness temperatures could be explained by larger Doppler boosting factors than those measured at parsec-scales if the innermost jet is pointing close to our line of sight. Alternatively, they may be an indication for the presence of other physical phenomena, as demonstrated by the RadioAstron observations of the quasar 3C 273 from Kovalev et al. (2016).
Further RadioAstron observations of OJ 287 have been performed in 2016, 2017, 2018, which are also in combination with quasi-simultaneous ground-based millimeter VLBI observations taken with the Event Horizon Telescope (in 2017 and 2018) at 230 GHz, as well as with the Global Millimeter VLBI Array (GMVA), including phased-ALMA, at 86 GHz. These observations, together with accompanying multi-wavelength coverage (Komossa et al. 2020(Komossa et al. , 2021 have the potential to either spatially resolve the hypothetical binary black hole system in OJ 287, or to put strong constrains on binary black hole and alternative models that predict the changing innermost jet position angle and overall periodic flaring activity that characterize this enigmatic source.

ACKNOWLEDGEMENTS
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 for the Instituto de Astrofísica de Andalucía (SEV-2017-0709). YYK, PAV, ABP acknowledge support from the Russian Science Foundation grant 21-12-00241. The RadioAstron project is led by the Astro Space Center of the Lebedev Physical Institute of the Russian Academy of Sciences and the Lavochkin Scientific and Production Association under a contract with the Russian Federal Space Agency, in collaboration with partner organizations in Russia and other countries. The European VLBI Network is a joint facility of of independent European, African, Asian, and North American radio astronomy institutes. Scientific results from data presented in this publication are derived from the EVN project code GA030E. This research is partly based on observations with the 100 m telescope of the MPIfR at Effelsberg. The VLBA is an instrument of the National Radio Astronomy Observatory, a facility of the National Science AIPS's task LPCAL (Leppanen et al. 1995) was used to estimate the instrumental polarization leakage, also known as D-terms, independently for each IF. The target source OJ 287 was used for computing the D-terms since it was the only source observed simultaneously with the SRT while providing the best parallactic angle coverage for the ground antennas. Fig. 9 shows the obtained D-terms for all the participating antennas during the perigee imaging session. We find consistent values across both IFs, confirming the reliability of the estimated D-terms. For the SRT we obtained (−0.88−4.55j)±(0.30+0.15j) and (−4.99 + 2.23) ± (0.63 + 0.55j), for RCP and LCP respectively. The dispersion across the two IFs has been used to estimate the errors.
Absolute calibration of the EVPA was obtained through comparison with VLA observations of OJ 287 at 22.295 GHz performed in 2014 May 1 under a different observing program (Marscher et al., private communication). Given the time span between our RadioAstron observations and the VLA ones we estimate an error for the absolute EVPA calibration, ∆χ, of the order of 10 • .

A.2. GMVA polarization calibration
Similarly to the RadioAstron observations, the calibration of the instrumental polarization leakage for the 86 GHz data was performed by employing AIPS' task LPCAL.
LPCAL assumes that the source structure can be described as a collection of polarization components each one with a constant fractional polarization, known as the self-similarity assumption. A good parallactic angle coverage is also required for the LPCAL fitting algorithm. Hence, ideally LPCAL should be used on sources with a simple polarization structure (preferably with a low fractional polarization), and a good parallactic angel coverage, requirements that are rarely met in millimeter VLBI observations. Alternative methods should then be considered when these requirements are not fully met. For this reason we have tested three different approaches to compute the D-terms for the 86 GHz GMVA observations.
First we computed the D-terms using OJ 287, which provided consistent values across all IFs, with a small dispersion of the order of ∼2% in amplitude and 13 • in phase. Following Casadio et al. (2019), we also computed the D-terms for each one of the bright and compact sources that were observed in the same session together with OJ 287 and had a large parallactic angle coverage (≥80 • ), taking the median values between the different intermediate frequency channels for each antenna as the representative D-terms. Fig. 10 shows the obtained D-terms for each source, including our target, and the median values. Finally we also tested which one of the D-terms obtained for each individual source provided the highest polarization image dynamic range across all observed sources. Out of these tests we found that the D-terms provided by OJ 287 yielded the polarization image with the highest dynamic range.
Calibration of the absolute EVPA for the GMVA 86 GHz observations was carried out through comparison with single dish IRAM 30 m telescope of OJ 287 as part of the POLAMI monitoring program, with an estimated uncertainty of the order of 5 • (Agudo et al. 2018b,a).

A.3. Estimations and uncertainties of the polarimetric parameters
The polarized flux density in this work is computed for all data sets by the standard formula P = Q 2 + U 2 , whereas, the uncertainty of P , σ P , is estimated by taking into account a calibration uncertainty of about 10% of the polarized flux density and a statistical error provided by the map thermal noise (Lico et al. 2014) as σ P = (0.1 × P ) 2 + rms 2 P , based on the polarization flux estimations that were obtained from the image domain of each image.
The fractional polarization percentage is determined as m = 100 (P/S), and the error on m is given by: The term σ D−term represents the systematic polarization calibration error, which is defined as (Roberts et al. 1994;: where σ amp is the standard deviation of the D-term amplitudes, N ant is the number of antennas, N IF is the number of the IFs, and N scan is the number of independent scans with different parallactic angles. For the GMVA image we have σ amp ∼ 2% , N ant = 7, N scan = 5, N IF = 16, S peak = 3.47 Jy, and S total = 5.07 Jy, which results in σ D−term = 4.62 mJy/beam. Similarly, for the RadioAstron data we have σ amp ∼ 1%, N ant = 13, N scan = 6, N IF = 2, S peak = 1.18 Jy, and S total = 2.38 Jy, which results in σ D−term = 1.92 mJy/beam. For the 15 GHz data, σ D−term ≈ 2 × rms P and the error for the absolute EVPA calibration is of the order of 5 • , whereas for the 43 GHz data, ∆m = 1% and ∆χ = 5 • % (Jorstad et al. 2005).