© The Authors 2025

1. Introduction

Among the various classes of active galactic nuclei (AGNs), BL Lacertae (BL Lac) objects stand out for their rapid, large-amplitude variability and significant polarization across multiple wavebands, which is attributed to relativistic jets aligned closely with our line of sight (e.g., Padovani 2017). A prominent member of this subclass is OJ 287, situated at a redshift of z = 0.306 (Stickel et al. 1989). Optical observations of OJ 287 extend back to the 1880s (Sillanpää et al. 1988), resulting in an exceptionally long light curve spanning nearly 150 years. This extensive dataset reveals quasiperiodic brightness fluctuations, including a prominent ∼60-year cycle (Valtonen et al. 2006) and recurrent, doubly peaked high-luminosity flares approximately every 12 years (Valtonen et al. 2006; Dey et al. 2018). These periodic variations are explained well by a supermassive black hole binary (SMBHB) model in which a secondary supermassive black hole follows a precessing, eccentric orbit around a more massive primary. Flares are generated each time a smaller component crosses the primary’s accretion disk (Lehto & Valtonen 1996; Sundelius et al. 1997; Dey et al. 2019). Within this framework, the well-known “impact outbursts” of OJ 287 (e.g., Valtonen et al. 2008) are the direct manifestations of these disk crossings.

In recent years, OJ 287 has gained significant attention due to its potential connection to the long-anticipated direct gravitational wave (GW) detections by the LIGO-VIRGO collaboration in 2015 (Abbott et al. 2016). The following reports of a GW background from pulsar timing arrays such as NANOGrav (Agazie et al. 2023) have reinforced the idea that SMBHBs are primary contributors to the GW background in the range of frequencies probed with pulsar timing. Also, OJ 287 proved to be a “testing ground” for multi-messenger studies based on electromagnetic and GW observations even before the first GW detections (see, e.g., Valtonen et al. 2008). (Sub)milliarcsecond-scale very long baseline interferometry (VLBI) images of OJ 287 provide crucial information for determining the stage of SMBHB evolution along the track defined by Begelman et al. (1980).

The parsec-scale jet of OJ 287 shows interesting behavior too; for example, it undergoes dramatic position angle (PA) variations of up to ∼130° (Agudo et al. 2012; Cohen 2017), which frequently coincide with major flaring episodes. These reorientations likely arise from a combination of geometric effects (such as jet precession) and dynamical processes, including the propagation of powerful shocks. Several mechanisms could explain these PA variations, including plasma instabilities and frame-dragging from a jet–disk misalignment, similar to patterns seen in other AGN jets (Abraham 2000; Caproni & Abraham 2004; Tateyama & Kingham 2004; Britzen et al. 2018; Qian 2018). Past VLBI observations have connected these PA rotations with structural changes in the jet, including component ejections, core flux variations, and transient features in the innermost regions (e.g., Agudo et al. 2012; Hodgson et al. 2017). However, attempts to model the high-frequency radio behavior using purely geometric or instability-driven approaches have yielded inconsistent results, as these models often fail to account for the complex, multi-zone nature of relativistic jets. For instance, Agudo et al. (2012) demonstrated that single-zone models cannot reproduce the observed rapid polarization variability, suggesting the necessity for multi-zone models to explain the simultaneous occurrence of polarization angle swings and flux density flares. Similarly, Hodgson et al. (2017) find that purely geometric models, such as those involving helical magnetic fields, are insufficient to explain the observed polarization variability, indicating that additional physical processes, like shock propagation or magnetic reconnection, must be considered to accurately model the polarization behavior in AGN jets. Additionally, the failure to detect the predicted October 2022 outburst (Komossa et al. 2023b) posed a significant challenge to certain SMBHB models. In response, refinements to these frameworks have incorporated additional factors, such as the geometry of the accretion disk, to improve predictive accuracy (Valtonen et al. 2023).

Space VLBI offers unprecedented angular resolution on the order of tens of microarcseconds, enabling detailed studies of AGN jets in total and linearly polarized intensity (e.g., Fuentes et al. 2023). The 10-meter space radio telescope (SRT) RadioAstron (Kardashev et al. 2013) operated from 2011 to 2019. It served as the primary instrument aboard the Spektr-R spacecraft, which functioned as its orbital platform. RadioAstron conducted observations at frequencies ranging from 1.6 to 22 GHz, and its orbit extended to ∼350 000 km, providing exceptional angular resolution at its shortest wavelengths. Results from the RadioAstron polarization key science program have significantly advanced our understanding of jet physics. The program encompassed diverse source classes, including BL Lac objects (Gómez et al. 2016), 0716+714 (Kravchenko et al. 2020), 0642+449 (Lobanov et al. 2015), 3C 345 (Pötzl et al. 2021), 3C 273 (Bruni et al. 2017), and more.

Dedicated observations of OJ 287 with RadioAstron at 22 GHz from April 4–5, 2014, achieved an unprecedented angular resolution of 12 μas, revealing a progressively bending jet extending down to ∼0.3 mas, providing key evidence of a predominantly poloidal magnetic field in the core, and suggesting that the parsec-scale jet maintains equipartition between particles and magnetic field energies (Gómez et al. 2022). Close in time, 1.68 GHz RadioAstron observations probed larger jet scales, resolving structures of up to ∼10 pc in size (Cho et al. 2024) and confirming a similar energetic equilibrium. Meanwhile, Global Millimeter VLBI Array (GMVA) observations at 3.5 mm (86 GHz) in April 2017, combined with new imaging techniques such as regularized maximum likelihood (RML), revealed a twisted, complex inner jet structure and polarization signatures indicative of recollimation shocks and a helical magnetic field configuration (Zhao et al. 2022). In the same year, the Event Horizon Telescope observed the source for the first time at 1 mm (Gómez et al. 2025). Nevertheless, for all these observations, at ∼20 μas resolution the source showed a jet extending only down to 0.2–0.3 mas.

In this paper we extend previous studies by presenting RadioAstron observations of OJ 287 at 22 GHz from April 26, 2016. For the first time, we image a continuous and sharply bent jet morphology down to 1 mas, enabled by the combination of high sensitivity, high angular resolution, and a favorable jet orientation. Although previous high-resolution VLBI studies (e.g.,Agudo et al. 2012; Hodgson et al. 2017) have resolved features within this scale, the multiband, ribbon-like structure revealed here had remained undetected, probably due to a limited dynamic range or a less optimal viewing geometry.

By juxtaposing high-resolution space VLBI images with quasi-simultaneous multiwavelength data, we investigate how the jet morphology and polarization structure have evolved.

2. Observations and data analysis

The RadioAstron mission observed OJ 287 at 22.24 GHz (λ = 1.35 cm) on April 25, 2016, from 16:00 to April 26 at 06:50 UT, spanning a total of approximately 15 hours, including about 12 hours on source (experiment code GG079C). The observations included calibrators 0716+714, 3C 345, 0059+581, 4C +38.41, and 4C +45.51 for ground-based radio telescopes. A total of 27 ground-based telescopes and the SRT participated in the experiment, providing extensive u-v coverage (see Table 1 and Fig. 1). The data were recorded in both left (LCP) and right (RCP) circular polarizations, with a total bandwidth of 64 MHz per polarization (512 Mbps with 2-bit sampling), split into four intermediate frequency (IF) bands for the ground antennas. For the SRT, the bandwidth was 32 MHz per polarization, split into two IFs and sampled at 1-bit. Due to operational issues, data from the Badary (BD), Ulsan (KU), and Yonsei (KY) stations were excluded from the analysis.

Fig. 1. (u, v) coverage of our observation at 22 GHz. Left: ground-only array. Right: full observations with ground-based arrays plus the SRT. Different colors show the ground stations providing the baseline with the SRT (see Table 1 for station codes). The maximum baseline lengths of the space baseline and the ground baseline are shown with dashed circles, which correspond to angular resolutions of 47 μas and 220 μas, respectively. Each point has been averaged with a 2-minute interval in both panels.

The SRT data were downlinked in real time to the RadioAstron tracking stations in Pushchino and Green Bank, with three-hour gaps to cool the onboard data downlink system of the Spektr-R satellite. The data were processed using the RadioAstron-dedicated version of the DiFX software correlator (Bruni et al. 2016), developed at the Max-Planck-Institut für Radioastronomie. The SRT provided baselines longer than the Earth’s diameter (D_⊕), resulting in a u − v extension up to 4.3 Gλ (4.6 D_⊕, Fig. 1), achieving an angular resolution of ∼47 μas. The spacecraft’s highly elongated orbit in the east–westdirection yields higher resolution along the east-southeast to west-northwest direction.

The calibration of the RadioAstron data was performed using NRAO’s AIPS software package (Greisen et al. 2003), following similar procedures to those described in Gómez et al. (2022) and Cho et al. (2024). The process involved solving for residual single- and multi-band delays, phases, and phase rates through gradual fringe fitting of the data. In the first iteration, RadioAstron was excluded, and a global fringe search was performed on the ground array with a solution interval of 60 s, with a signal-to-noise ratio (S/N) threshold of 5. In subsequent fringe fitting that included the SRT, a lower S/N threshold of 4 was used for space–ground baselines to account for their typically lower fringe amplitudes. Effelsberg and the Green Bank Telescope were used as reference antennas throughout the process.

Manual phase calibration was done iteratively, starting with the ground array and then progressively incorporating the RadioAstron SRT. Fringe fitting was performed between the ground telescopes and the SRT in each segment separately. The ground array was then coherently combined through baseline stacking to increase the S/N of possible fringe detections to the SRT (e.g., Bruni et al. 2017).

To account for the spacecraft’s residual acceleration near perigee, we adopted different solution intervals ranging from 10 to 180 seconds during fringe fitting. While the acceleration term was checked during correlation, it was not applied, as it proved unnecessary in practice. On shorter baselines, the higher S/N permitted shorter solution intervals, making acceleration corrections straightforward in post-correlation analysis. On longer baselines, the effect of acceleration was negligible and did not require correction. Additionally, since the longest projected baselines were less sensitive to acceleration, we optimized the total data bandwidth by combining IFs. Polarization alignment was performed using RLDLY task on a bright scan of the source. Amplitude calibration was performed by loading and applying the system temperature and gain curve information for each antenna, and deriving amplitude corrections using APCAL, while accounting for issues at specific antennas. The amplitude solutions were smoothed using SNSMO. Due to the low S/N of many scans, instead of performing a bandpass calibration, the outer six spectral channels on each side of the IFs were flagged to eliminate signal attenuation and leakage.

To complement our RadioAstron data, we also incorporated archival 43 GHz Very Long Baseline Array (VLBA) observations to study the jet morphology over time. The 43 GHz data used in this work (presented in Fig. 3) were obtained with the VLBA, spanning the period from 2014 to 2017 with annual cadence, within the VLBA-BU-BLAZAR program¹. The program conducts regular monthly observations of a sample of γ-ray bright AGNs. A detailed description of the observations and data reduction can be found in Jorstad et al. (2017) and Weaver et al. (2022). For this study, we reimaged a total of four images.

2.1. Imaging and model-fitting: Total intensity

Image reconstruction was performed using the RML method, implemented in the eht-imaging software package ehtim (Chael et al. 2016, 2018). Unlike traditional techniques, the RML approach avoids using the inverse Fourier transform of the visibilities, V, during the imaging process. Instead, it integrates forward modeling with regularization to derive an image, I, that minimizes the objective function:

$$\begin{array}{c}\hfill J(I)={\displaystyle \sum _{\mathrm{data}\mathrm{terms}}}{\alpha }_D{\chi }_D^2(I,V)-{\displaystyle \sum _{\mathrm{reg}.\mathrm{terms}}}{\beta }_R{S}_R(I).\end{array}$$(1)

The first term of this equation promotes consistency with the observed data, while the second term enforces regularization to ensure a physically meaningful image. The relative contributions of these terms are controlled by the hyperparameters α_D and β_R, which balance data fidelity against regularization. Using this method, we achieve higher angular resolution and image fidelity compared to traditional techniques like CLEAN (e.g., Event Horizon Telescope Collaboration 2019, 2022).

The imaging 22 GHz RadioAstron data obtained on April 25, 2016, was carried out in several steps, beginning with the reconstruction of the ground array. We initialized the process with a Gaussian prior image (full width at half maximum = 150 μas), refining the solution iteratively. This was done across a 2 mas field of view, discretized onto a 300 × 300 pixel grid. We incorporated closure phases and logarithmic closure amplitudes, accounting for non-closing systematic errors of 1% of visibility amplitude. After achieving convergence, we performed self-calibration to derive gain solutions and improve image fidelity. In subsequent iterations, we expanded our approach. We first included visibility phases and amplitudes, followed by a comprehensive integration of both closure quantities and complex visibilities. The final stage involved applying this strategy to include the SRT, using the ground array image as a prior.

To determine optimal relative weights for the regularization terms, we systematically varied their values across a search grid and selected the combination yielding the highest image fidelity and best data fit. This involved testing various values for relative entropy, total variation (tv), total squared variation (tv2), total variation ℓ2 with logarithmic regularizer (tv2log), and the ℓ1 norm. The tv2log regularizer combines the tv norm with an ℓ2 norm and incorporates a logarithmic transformation, balancing between promoting sparsity and smoothness while handling large dynamic ranges in image intensities (see also Savolainen et al. 2023). The final combination of values were: relative entropy: 0.1, tv: 0, tv2: 0.01, tv2log: 0.05, and ℓ1: 1, resulting χ_phas² : 1.04, and χ_logamp²: 1.98. Details of the parameter selection process are provided in Appendix A. We applied the same imaging procedure to the 43 GHz data, with the resulting images shown in Figs. 2 and 3.

Fig. 2. 22 GHz image of OJ 287 on April 25, 2016. Left: total intensity image shown in brightness temperature (Tb) units, calculated using the relation Tb = Iνc2/(2kν2), where Iν is the specific intensity per pixel. The pixel size used is 5 μas, derived from the total field of view and the image resolution (320 pixels across). The peak Tb reaches ∼6 × 1012 K in the core region. Right: total intensity contours (dashed white lines) spaced logarithmically between 0.01% and 0.1% of the peak intensity. The color map shows the magnitude of the linearly polarized flux. Regions where polarization is not reliably detected downstream have been masked. White ticks indicate the EVPA orientation. Yellow circles mark the centroids of the two-dimensional Gaussian-fitted components. Polarized features P1, P2, and P3 are shown in purple. The central dashed green curve traces the jet ridgeline, as derived from the smoothed image (see Sect. 3).

Fig. 3. Multi-epoch 43 GHz VLBI images of OJ 287 observed on May 3, 2014, April 11, 2015, April 22, 2016, and April 16, 2017. Each panel shows the source morphology in total intensity as reconstructed by ehtim. The color scale is the same for all days and is in units of brightness temperature; an angular scale indicator is included in the bottom-right corner of each image (200 μas) for reference. The χ2 statistics used as diagnostics of the model fit quality – the closure phase reduced χ2 (χcphase2) and the logarithmic closure amplitude reduced χ2 (χlogamp2) – are displayed in the bottom-left corner of each image. The ridgelines, delineating the jet’s axis for each epoch, are overlaid as dashed lime lines, with a black outline for enhanced visibility. In this image sequence we can clearly see the jet’s remarkable PA evolution over the four years.

We imported the fully calibrated RadioAstron data into Difmap (Shepherd 1997), an interactive program for synthesis imaging. Using the MODELFIT algorithm, we parameterized the jet brightness distribution with two-dimensional circularGaussian components. The uncertainties for each component parameter were formally assessed based on the local S/N in the surrounding image (Fomalont 1999; Lobanov 2005; Schinzel et al. 2012). For the flux density measurements, we adopted a more conservative 20% uncertainty, to account for the increased calibration uncertainties inherent to space VLBI observations at 22 GHz. This choice is more appropriate than the 10% typically used in ground-based VLBI at lower frequencies (Lister et al. 2009), given the limited system temperature and gain information available for some antennas and the lower S/N on the space baselines. All parameters for the fittedGaussian components are listed in Table 2.

2.2. Imaging and model-fitting: Polarization

Following total intensity imaging and polarimetric calibration, we reconstructed polarized images using ehtim. We minimized the objective function in Eq. (1), incorporating both polarimetric visibilities ($\stackrel{\sim }{P}=\stackrel{\sim }{Q}+i\stackrel{\sim }{U}$) and the visibility-domain polarization ratio ($\stackrel{\sim }{m}=\stackrel{\sim }{P}/\stackrel{\sim }{I}$), the latter being immune to residual station gain errors from Stokes I imaging. For image reconstruction, we employed two regularizers: (1) the Holdaway-Wardle regularizer (hw; Dayton et al. 1990) to constrain pixel polarization fraction below 0.75 (the theoretical maximum for synchrotron radiation), and (2) the polarimetric total variation (ptv) regularizer (Rudin et al. 1992) to ensure smooth polarization transitions. We reconstructed total intensity and linearly polarized intensity images independently, solving for the fractional polarization (m) and the electric vector position angle (EVPA; χ) in each pixel. The iterative gradient descent process used the previous iteration’s output blurred by a 20 μas Gaussian kernel, with final values of hw:1 and ptv:1 chosen for optimal data fit and image fidelity. The pipeline alternated between polarimetric objective function minimization and D-term calibration, maximizing consistency between self-calibrated data and corrupted image reconstructions. After D-term solutions were obtained (see Appendix A), we blurred the polarimetric image and repeated the imaging-calibration cycle untilconvergence.

In addition to correcting for instrumental polarization, VLBI polarimetric analyses require accurate absolute EVPA calibration. To this end, we applied a correction based on a close-in-time VLBA observation at 43 GHz from the BEAM-ME project (successor to the VLBA-BU-BLAZAR program), which reported an EVPA of −12° ±2° for OJ,287 on April 22, 2016. This value was used to calibrate the RadioAstron EVPA measurements. As an independent reference, we also considered a single-dish Effelsberg observation of OJ,287 at 10.45 GHz, conducted on April 12, 2016, which measured EVPA = 3° .36 ± 0.75 (Myserlis et al. 2018), obtained within the framework of the MOMO monitoring program (Komossa et al. 2015, 2023a). Finally, we computed the net polarization, |m|_net, which reflects the degree of alignment in polarization direction across the source, and the average linear polarization fraction, ⟨|m|⟩, which quantifies the typical polarization strength, independent of the directional coherence following Event Horizon Telescope Collaboration (2021):

$$\begin{array}{c}\hfill {|m|}_{\mathrm{net}}=\sqrt{{\left({\displaystyle \sum _i}{Q}_i\right)}^2+{\left({\displaystyle \sum _i}{U}_i\right)}^2}/{\displaystyle \sum _i}{I}_i\end{array}$$(2)

$$\begin{array}{c}\hfill ⟨|m|⟩={\displaystyle \sum _i}\sqrt{Q_i^2+U_i^2}/{\displaystyle \sum _i}{I}_i.\end{array}$$(3)

These equations yield 0.45% and 5.5%, respectively. Although several bright calibrator sources (0716+714, 3C 345, 0059+581, 4C +38.41, and 4C +45.51) were observed as part of the experiment to serve as fringe finders and potential D-term calibrators, we ultimately derived the instrumental polarization calibration using OJ 287 itself, which provided more reliable solutions due to its strong, compact polarized structure and longer integration time.

3. Results

3.1. Jet morphology

The left panel of Fig. 2 shows the total intensity image of OJ 287 from our RadioAstron observations. The core region exhibits a remarkable and atypical structure, deviating from the standard core-dominated morphology typically observed in blazars. The jet emanating from the core initially propagates northwest before undergoing a sequence of dramatic directional changes. The first significant bend occurs at ∼100 μas from the core, where the jet trajectory sharply changes direction. Subsequently, two additional bends are observed at ∼300 μas and ∼550 μas to the west, forming a ribbon-like morphology. This complex trajectory may reflect a possible interaction with the surrounding medium, precession of the jet axis, inherent instabilities within the jet flow, or all mechanisms working in combination.

Our model-fitting analysis reveals four distinct Gaussian components. The core region is resolved into: C1, with a flux density of 1.5 Jy, located upstream and likely representing the VLBI core, and C2, appearing as the brightest knot with a flux density of 1.8 Jy. We note that compared to our 2014 observations where C1 and C2 showed flux densities of 0.5 Jy and 1.17 Jy, respectively (Gómez et al. 2022), the current enhanced brightness of C1 suggests a recent ejection event where a new component is becoming optically thin at 22 GHz. Further downstream, components B1 and B2 are located at the bending points, marking the transition from the compact core region to the broader jet structure. We computed the brightness temperatures for each component using the relation T_b = 1.22 × 10¹²(S_ν/(θ_obs²ν²))(1+z) (K) (e.g., Pushkarev & Kovalev 2012), where S_ν is the component flux density in Jy, θ_obs is the component full width at half maximum in mas, and ν is the observing frequency in GHz, and z is the redshift. The values decrease systematically with distance from the core down the jet: T_b, C1 = (1.3 ± 0.28)×10¹³ K, T_b, C2 = (1.2 ± 0.31)×10¹² K, T_b, B1 = (4.4 ± 1.5)×10¹¹ K, and T_b, B2 = (8.6 ± 3.0)×10⁹ K.

The exceptionally high T_b of C1 exceeds both the inverse Compton limit (∼10¹² K; Kellermann & Pauliny-Toth 1969) and the equipartition brightness temperature (∼10¹¹ K) limits (Readhead 1994), requiring Doppler factors δ_j ∼ 10–30 for reconciliation. This agrees with estimates from VLBA-BU-BLAZAR monitoring at 43 GHz (δ_j = 8.6 ± 2.8; Weaver et al. 2022), suggesting a jet viewing angle θ_j ≈ 3° −8°. Such extreme brightness temperatures could indicate strong Doppler boosting, departure from equipartition conditions, or potentially more exotic emission mechanisms such as relativistic proton emission (Kardashev 2000) or coherent emission processes (Benford & Lesch 1998). See also the results and discussion in Kovalev et al. (2016, 2020). Meanwhile, the decreasing T_b trend along B1 and B2 reflects the expected jet expansion and energy dissipation with distance from the core (e.g., Röder et al. 2025).

Our polarimetric analysis revealed a prominent structure in the inner jet (extending within the first ∼100 μas), as shown in the right panel of Fig. 2. In this region, the EVPAs are perpendicular to the jet axis, highlighting the dominance of the poloidal magnetic field component, favoring the development of kink instabilities (e.g., Nakamura & Meier 2004). We describe the polarization structure in terms of three distinct features, namely P1, P2, and P3. P1, the innermost component, exhibits a polarized intensity of P₁ = 0.07 ± 0.01 Jy, a fractional polarization of m₁ = 3.0%±0.3%, and an EVPA of χ₁ = 42° ±5.0°. P2, located farther downstream, exhibits P₂ = 0.07 ± 0.01 Jy, m₂ = 18.5%±2.0%, and χ₂ = 48° ±5.0°. The distinction between P1 and P2 was made based on their proximity, morphology, and the presence of a gradient in both polarized intensity and total intensity. While there is a continuity in emission between P1 and P2, P1 is identified as the brightest feature in the southernmost region, whereas P2 corresponds to a distinct intensity peak farther downstream. P3, the outermost feature, displays the lowest polarized intensity (P₃ = 0.05 ± 0.01 Jy), but maintains significant fractional polarization (m₃ = 20.5%±3.3%), with an EVPA of χ₃ = 38.8° ±5.0°. The elevated fractional polarization values in P2 and P3 indicate that, despite their lower polarized intensities, the magnetic field retains a significant degree of alignment, which could arise from the effects of shocks or compression in the jet flow. All uncertainties in P, m, and χ were calculated following Gómez et al. (2022).

3.2. Jet ridgeline analysis

To quantify the jet structure and PA of OJ 287, we performed a ridgeline analysis as described in Zhao et al. (2022). First, we generated a smoothed version of the reconstructed image by applying a Gaussian blur with a width equivalent to twice the image resolution (47 μas). The smoothed images were transformed into polar coordinates, centered on the jet origin. To construct the ridgeline, the polar-transformed images were sliced radially, and transverse intensity profiles were extracted along the jet axis. For each slice, we fitted a Gaussian function to the intensity profile to identify the peak flux density position, which was subsequently transformed back into Cartesian coordinates. Finally, the resulting set of peak positions was interpolated using a cubic spline to obtain a continuous ridgeline tracing the jet path. All uncertainties were estimated using the covariance matrix from Gaussian fits to derive the uncertainty of the peak position (σ_x0). The radial uncertainty (σ_r) was set equal to σ_x0, while the angular uncertainty (σ_θ) was approximated as 1/r, accounting for the decreasing angular precision with increasing radial distance. Failed Gaussian fits were excluded from the analysis. Using standard error propagation techniques, the ridgeline positions and their associated uncertainties were transformed back into Cartesian coordinates. In addition to the ridgelines derived from the data analyzed in this work, we also included ridgeline data for OJ 287 from four selected epochs at 43 GHz, obtained from the publicly available VLBA-BU-BLAZAR database. These epochs span April 2014 to 2017 (May 3, 2014, April 11, 2015, April 22, 2016, and April 16, 2017). Their comparison allowed us to investigate annual changes in jet orientation and check their agreement with theoretical predictions (see the discussion in Sect. 4.2).

As can be seen in Fig. 3, the jet demonstrates a clear temporal evolution. In 2014, OJ 287 appears smooth and extended, displaying modest curvature and reaching ∼600 μas apparent length from the core. The 2015 epoch reveals a shorter jet configuration, with emission confined to ∼300 μas. By 2016, the jet regains its extended structure to ∼600 μas, showing also strong morphological agreement with 22 GHz RadioAstron image, reinforcing the jet structure across frequencies. The 2017 epoch exhibits again a compact, bright inner region with minimal downstream emission, similar to the confined structure observed in 2015. These variations in apparent jet extension and morphology can result from multiple factors, including changes in Doppler boosting or jet precession due to instabilities, Lense-Thirring effects, or the presence of an SMBHB system. Additionally, differences in data quality, such as variations in dynamic range or noise level, can affect the detectability of extended jet features and may also contribute to the observed morphological differences.

4. Discussion

4.1. OJ 287 multi-scale jet bending

The revelation of a ribbon-like morphology with three distinct bends within 650 μas in the jet of OJ 287 revealed by RadioAstron observations is a unique occurrence. Unlike previous high-resolution VLBI studies of OJ 287 (e.g., Murai et al. 2017; Lico et al. 2022; Gómez et al. 2022; Zhao et al. 2022), which could only resolve the first bend, our RadioAstron observations from April 25, 2016, have provided the first direct imaging of the full multi-scale jet reorientation within theacceleration and collimation zone. According to Marscher et al. (2008), Homan et al. (2015), the acceleration and collimation zone in blazars typically extends to ∼10⁴–10⁶ Schwarzschild radii (R_S) from the central black hole. For OJ 287, using the latest mass estimate for the central supermassive black hole of ∼10⁸ M_⊙ (Komossa et al. 2023b,a) and a luminosity distance of 1.6 Gpc, this corresponds to an angular scale of approximately 6.1–610 μas on the sky. The resolution of ∼47 μas achieved in this work (improved to ∼15 μas after ehtim imaging), offered us the ability to capture these multiple bends and reveal that the jet undergoes a substantially more complex evolution than previously assumed.

4.2. Jet position angle variations and SMBHB precession predictions

Position angle variations in the OJ 287’s jet can provide crucial insights into properties of the binary black hole in this system. Dey et al. (2018) conducted a Bayesian analysis to explore how the PA variations observed in the source at multiple frequencies (15, 43, and 86 GHz) could be explained by precession mechanisms linked to the SMBHB system dynamics. Their study explored two potential origins for the precession: the evolution of the primary black hole’s spin, which attributes the variations to Lense-Thirring precession caused by the primary black hole’s spin (spin model), and the precession of the angular momentum of the inner accretion disk, which links to the hydrodynamical evolution of the inner disk’s angular momentum under the gravitational influence of the secondary black hole (disk model). While both models describe the data well, the accretion disk precession model provided a more consistent description of PA variations observed at higher frequencies (e.g., 86 GHz in Dey et al. 2018), making it particularly relevant for comparison with our data here.

To test these theoretical predictions, we analyzed jet PAs using both our RadioAstron observations and the VLBA-BU-BLAZAR data at 43 GHz. We derived the mean PA for each epoch by averaging ridgeline points within the first bright jet component from the core, consistent with the model resolution limit and the region of the jet consistently visible in all epochs. To establish robust uncertainty estimates, we implemented a statistical bootstrapping method. This involved randomly resampling with replacement from our original PA measurements 10 000 times, creating 10 000 simulated datasets. For each of these simulated datasets, we calculated the circular standard deviation, producing a distribution of possible standard deviations. This approach provides a reliable error estimation that accounts for the circular nature of angle measurements (where 0° and 360° are identical), while reducing potential biases from outliers or epochs with limited sampling. For the 22 GHz 2016 data, the presence of RadioAstron ensured angular resolutions comparable to the 86 GHz observations.

Our analysis yielded an innermost PA (≤100 μas translate to ∼0.45 pc projected distance)of −25° ±5° for the RadioAstron 2016 observations. For the VLBA-BU-BLAZAR 43 GHz data, we obtained values of −12° ±6° on May 3, 2014, −38° ±6° on April 11, 2015, −6° ±5° on April 22, 2016, and −32° ±6° on April 16, 2017. These measurements reveal a clear evolution of the jet PA within ∼100 μas, with a change of ∼25° in the 43 GHz measurements across the four epochs. This change represents a significant reorientation of the jet direction. In the context of the disk model, the PA variations observed in our VLBI datasets align well with the predictions only for the first 2 epochs, with the 2016 data (at both 43 and 86 GHz) deviating from the models, while the 2017 data follow the 86 GHz prediction as shown in Fig. 4.

Fig. 4. Fits to the variations in the PA of the radio jet at 43 GHz and 86 GHz based on the expected jet precession in the disk model by Dey et al. (2018). The red and black lines correspond to the model predictions at 43 GHz and 86 GHz, respectively. The innermost PA value from RadioAstron is marked with black diamond, and VLBA-BU-BLAZAR data points are marked with red circles.

4.3. Relation with multiwavelength flares and jet dynamics

A major optical outburst occurred in OJ 287 on December 5, 2015 (MJD 57363, corresponding to 2015.93), with the source reaching its brightest optical level in three decades. This event is also known as the “Centenary Flare,” as it was predicted with astonishing accuracy by the SMBHB model (Valtonen et al. 2006). The model suggested that the outbursts are attributed to the secondary black hole periodically impacting the accretion disk of the primary supermassive black hole, producing thermal bremsstrahlung radiation visible as bright optical flares (Lehto & Valtonen 1996; Prochaska & Wolfe 1998). It is expected that such interactions in binary systems can enhance jet activity through multiple mechanisms: disk perturbations that increase accretion rates (Lehto & Valtonen 1996; Sundelius et al. 1997; Ricarte et al. 2016), magnetic field intensification associated with the orbital motion of the secondary (e.g., Palenzuela et al. 2010; Gold et al. 2014), and direct gravitational torques capable of influencing the jet-launching region (e.g., Gorjian et al. 2008; D’Orazio et al. 2013).

Subsequently, on March 15, 2016 (2016.21), OJ 287 exhibited another strong flare detected in both the optical and radio bands, originating from the jet. The optical flare, comparable in strength to the December 2015 outburst, displayed polarization characteristics indicative of synchrotron emission, suggesting a jet-driven origin rather than a direct disk impact (Gupta et al. 2017; Chang & Murray 2018). Simultaneously, a significant radio flare was observed, likely triggered by a propagating shock wave within the relativistic jet (Gupta et al. 2017). Kinematic analysis of VLBA-BU-BLAZAR data (Weaver et al. 2022) confirmed the ejection of component B15 at 2016.23 ± 0.10, characterized by a proper motion of 0.28 ± 0.02 mas/year (5.38 ± 0.32 c). One month before this ejection, the core experienced a prominent flare where the flux density increased from 1.7 Jy to 3.9 Jy (130%), consistent with the typical behavior observed when a new component traverses the radio core (Marscher et al. 2008).

An independent study by Lico et al. (2022) at 86 GHz detected the same component, designating it as K, with a measured propagation speed of ∼0.32 mas/year (∼4.8 c), yielding an ejection time of ∼2016.3. Their analysis suggests that the TeV flaring activity observed February 1–4, 2017, was triggered when this new jet feature passed through a recollimation shock S1 at ∼0.1 mas from the radio core, a stationary jet feature mentioned in multiple previous studies (Hodgson et al. 2017; Gómez et al. 2022). Figure 5 illustrates the emergence and propagation of this new component, tracking both its motion away from the core and its flux density evolution over time, while also showing the temporal relationship between optical, radio and TeV flaring events.

Fig. 5. Flux density and radial core separation evolution of jet components within 0.2 mas of the core in OJ 287 at 22, 43, and 86 GHz from September 2015 to April 2017. Top panel: Flux density evolution over time, showing the core (red circles), a knot identified as B15 at 43 GHz and K at 86 GHz (green square and yellow triangle, respectively), and components C1 (magenta diamond) and C2 (purple cross) identified in this work at 22 GHz. Error bars indicate measurement uncertainties. Bottom panel: Radial core separation projected distances. A TeV flare (dashed cyan vertical line) and the centenary optical flare (dashed gray vertical line) are marked. The dashed red vertical line and the gray shaded region show the estimated ejection time and its uncertainty for knot B15 or K (Weaver et al. 2022).

4.4. Comparison with single-dish radio polarization monitoring

Myserlis et al. (2018) demonstrated that OJ 287 underwent a prolonged clockwise rotation of its radio EVPA during 2016, based on single-dish Effelsberg observations carried out within the framework of the MOMO monitoring program (Komossa et al. 2015, 2023a). Their analysis revealed that this rotation occurred within the jet core as observed by the VLBA-BU-BLAZAR monitoring program at 43 GHz. This region corresponds to a distance of approximately ∼10⁴–10_S^R from the central black hole (assuming a mass of 10⁸–10¹⁰ M_⊙), placing it in the outer acceleration and collimation zone. Our RadioAstron data now allow us to examine this region with enhanced resolution.

The radio EVPA evolution in 2016 reported in Myserlis et al. (2018) was attributed to the possible bending of the inner jet within the 43 GHz core, which has a projected size of about 0.15–0.2 mas. Our RadioAstron imaging results of the same year (Fig. 2) show that there is indeed strong bending of the inner jet at those spacial scales.

A noteworthy detail emerges when comparing the EVPA rotation rate in 2016 reported in Fig. 4 of Myserlis et al. (2018) with the inner jet morphology of the same year as revealed by RadioAstron. The EVPA rotation rate decreased from 4°/day to 1°/day shortly after the onset of the rotation event, as seen only in the 10.45 GHz data of that figure². This deceleration can be attributed to a more tightly wound helical structure near the jet origin, precisely what our 2016 RadioAstron image reveals.

To quantitatively assess this relationship, we applied ridgeline analysis to our 2016 RadioAstron image. Assuming the ridgeline represents a projected helical trajectory, its resemblance with a sine wave of increasing period along the (mean) jet direction suggests that the step of the helical trajectory also increases gradually along that direction. We estimated the rotation rate of the helix along the (mean) jet direction by splitting the ridgeline in seven segments, each one corresponding to a quarter of the (increasing) period, which corresponds to a 90° change along the helix for each segment. The resultant helical trajectory rotation rate is shown as the black line in Fig. 6.

Fig. 6. Comparison of the rotation rate of the polarization angle rotation event at 10.45 GHz shown in Myserlis et al. (2018) and the rotation rate of the helical jet shape along the jet direction as seen in the 2016 RadioAstron observations of OJ 287. In both cases we see a tighter winding of the jet helical, with it bending closer to the jet base.

For comparison, the red line in Fig. 6 represents the absolute 10 GHz EVPA rotation rate reported in Fig. 4 of Myserlis et al. (2018), converted from °/day versus days to °/μas versus μas. To change from temporal (days) to spacial (μas) scales, we assumed that the polarized component travels at a relativistic velocity and also moved from the observed to the emission reference frame using the mean Doppler factor and the redshift of the source as done also in Eq. 1 of Myserlis et al. (2018). Finally, for theconversion of linear to angular distance we adopted 4.48 pc/mas, as expected at the redshift of OJ 287.

The striking similarity between the two rotation rates in Fig. 6 indicates that the polarized component responsible for EVPA rotation indeed traversed a more tightly wound jet bend initially, consistent with the increased ridgeline rotation rate near the jet base. Since the exact region of the jet responsible for the single-dish EVPA rotation might be different from the one depicted in the 2016 RadioAstron image, we arbitrarily shifted the EVPA rotation rate in Fig. 6 by 43 μas to align the two peaks and enable a direct comparison.

In summary, our high-resolution RadioAstron images from 2016 provide a cohesive picture that aligns with the high-cadence single-dish polarization monitoring data and the EVPA variability at the inner jet within the 43 GHz core, as reported by Myserlis et al. (2018).

5. Summary and conclusions

We have presented new, high-resolution RadioAstron observations of OJ 287 at 22 GHz obtained on April 25, 2016. These observations, conducted at an angular resolution of ∼47 μas, reveal complex jet dynamics on extreme spatial scales. Our main findings are as follows:

1. Ribbon-like jet structure with multiple sharp bends. Our space VLBI images reveal an exceptionally complex inner jet morphology in OJ 287 featuring three distinct bends within ∼650 μas of the VLBI core; this is the first time such a morphology has been observed. This observation provides direct evidence that the full multi-scale reorientation of the jet occurs within the acceleration and collimation zone. The observed structure indicates a more complex jet evolution than previously assumed, possibly influenced by jet precession, magnetohydrodynamic instabilities, or interactions with the surrounding medium.

2. Extreme brightness temperatures exceeding 10¹³K. Our model-fitting analysis identifies four key components: C1 (the VLBI core), C2 (further downstream), and B1 and B2 (located at bending points). The core exhibits an exceptionally high brightness temperature of T_b ∼ 1.3 × 10¹³ K, exceeding the inverse Compton and equipartition limits. This implies strong Doppler boosting (δ ∼ 10 − 30), which is consistent with prior estimates of OJ 287’s jet viewing angle (θ_j ≈ 3° −8°).

3. Ordered magnetic field revealed by polarization structure. Our polarimetric analysis shows EVPAs predominantly transverse to the jet axis in the core region, indicating a dominant poloidal magnetic field component.

4. Jet PA variations partially support precession models. We tracked the evolution of the innermost jet PA using VLBA-BU-BLAZAR observations at 43 GHz from 2014 to 2017. The jet shows a PA evolution of ∼30° over four years. These variations are partially consistent with precession models, particularly the disk precession model of Dey et al. (2018), which aligns with some of our 43 GHz data points. However, we note that not all epochs fit the predictions, andfurther observations would be necessary to draw firm conclusions. Therefore, while our data do not directly confirm the SMBHB scenario, they offer valuable structural and temporal constraints that can inform and refine such models.

5. Connection between the jet’s new components and high-energy flares. Existing kinematic data at 43 and 86 GHz combined with our 22 GHz Gaussian model-fitting confirm the emergence of a new jet component (B15 or K) around March 2016, shortly after the December 2015 optical“Centenary Flare.” This component later interacted with a stationary recollimation shock (S1), and this interaction coincided with a TeV flare in early 2017. These findings further support a physical connection between disk impacts, jet activity, and high-energy emission in OJ 287.

6. Consistency between jet bending and EVPA rotation. A comparison with the high-cadence single-dish polarization monitoring dataset presented in Myserlis et al. (2018) shows a remarkable agreement between the observed EVPA rotation and the jet ridgeline rotation rate. This suggests that the EVPA swings observed in 2016 trace an underlying helical bending of the jet.

Acknowledgments

We thank Etienne Bonnassieux for his useful comments and fruitful discussions on the manuscript. Authors E. Traianou, and J. L. Gomez acknowledge financial support from the Severo Ochoa grant CEX2021-001131-S funded by MCIN/AEI/10.13039/501100011033. The work at the IAA-CSIC was supported in part by the Spanish Ministerio de Economía y Competitividad (grant number PID2022-140888NB-C21). I.C. is supported by the KASI-Yonsei Postdoctoral Fellowship program. This study makes use of VLBA data from the VLBA-BU Blazar Monitoring Program (BEAM-ME and VLBA-BU-BLAZAR; http://www.bu.edu/blazars/BEAM-ME.html), funded by NASA through the Fermi Guest Investigator grants, the latest is 80NSSC23K1508. 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 State Space Corporation ROSCOSMOS, in collaboration with partner organizations in Russia and other countries. The VLBA is an instrument of the National Radio Astronomy Observatory. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated by Associated Universities, Inc. This research has made also use of data obtained with the Global Millimeter VLBI Array (GMVA), which consists of telescopes operated by the MPIfR, IRAM, Onsala, Metsähovi, Yebes, the Korean VLBI Network, the Green Bank Observatory and the Very Long Baseline Array (VLBA). The VLBA is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The data were correlated at the special version (Bruni et al. 2016) of the RadioAstron DiFX correlator (Deller et al. 2011) of the MPIfR in Bonn, Germany. This research has made use of NASA’s Astrophysics Data System. MW is supported by a Ramón y Cajal grant RYC2023-042988-I from the Spanish Ministry of Science and Innovation. LIG gratefully acknowledges support by the Chinese Academy of Sciences PIFI program, grant No. 2024PVA0008. ABP is supported in the framework of the State project “Science” by the Ministry of Science and Higher Education of the Russian Federation under the contract 075-15-2024-541. MML, YYK, and APL weresupported by the M2FINDERS project that had received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 101018682). This work is partly based on observations carried out with the IRAM 30m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). The research at Boston University was supported in part by NASA Fermi GI grant 80NSSC20K1567.

References

Appendix A: Image fidelity and data issues

During the data processing, we flagged multiple stations and IFs due to poor data quality. Specifically, stations KT, KU, BD, and RO were excluded because of the absence of fringes, while the first and fourth IFs of station ON were flagged for corrupted signals. We used GB as the reference antenna throughout the analysis. Station LA displayed erratic amplitude fluctuations despite having stable phases. Therefore, we adopted a special strategy: LA was included in closure-phase imaging but excluded from amplitude-based imaging, and it was used only during self-calibration steps, ensuring that LA amplitudes followed the more robustly constrained model derived from the other antennas.

Fig. A.1. Selected closure amplitudes and phases from coherently averaged visibilities on triangles as a function of time.

Fig. A.2. Self-calibrated visibility amplitudes and phases as a function of uv distance from the RadioAstron observations of OJ 287 from April 24–25, 2016, at 22 GHz. The magenta points represent the fit to the data using the ehtim model derived from the imaging process. Reliable space-ground fringe detections were achieved up to a projected baseline length of 4.6 Earth diameters. The dot-dashed black line separates the ground-only baselines from those involving the SRT.

Fig. A.3. Multiplicative gain correction factors at each station.

Appendix B: Instrumental polarization

We note that station JB was flagged due to an excessively large gain value of $2.85\pm 0.4$. For polarization imaging, we flagged additional stations HH, KT, MK, NT, PT, EF, and SR because of insufficient parallactic angle coverage or unstable signal fluctuations (particularly noted at PT). The D-term calibration for RadioAstron yielded $(-0.0174+0.0111i)$ for RCP and $(0.0198-0.0001i)$ for LCP.

Fig. B.1. Time evolution of selected baseline polarization amplitudes and phases derived from coherently averaged visibility measurements.

Fig. B.2. D-term solutions for each station and each polarization. Labels mark the different radio telescopes.

All Tables

All Figures

Fig. 1. (u, v) coverage of our observation at 22 GHz. Left: ground-only array. Right: full observations with ground-based arrays plus the SRT. Different colors show the ground stations providing the baseline with the SRT (see Table 1 for station codes). The maximum baseline lengths of the space baseline and the ground baseline are shown with dashed circles, which correspond to angular resolutions of 47 μas and 220 μas, respectively. Each point has been averaged with a 2-minute interval in both panels.

In the text

Fig. 2. 22 GHz image of OJ 287 on April 25, 2016. Left: total intensity image shown in brightness temperature (Tb) units, calculated using the relation Tb = Iνc2/(2kν2), where Iν is the specific intensity per pixel. The pixel size used is 5 μas, derived from the total field of view and the image resolution (320 pixels across). The peak Tb reaches ∼6 × 1012 K in the core region. Right: total intensity contours (dashed white lines) spaced logarithmically between 0.01% and 0.1% of the peak intensity. The color map shows the magnitude of the linearly polarized flux. Regions where polarization is not reliably detected downstream have been masked. White ticks indicate the EVPA orientation. Yellow circles mark the centroids of the two-dimensional Gaussian-fitted components. Polarized features P1, P2, and P3 are shown in purple. The central dashed green curve traces the jet ridgeline, as derived from the smoothed image (see Sect. 3).

In the text

Fig. 3. Multi-epoch 43 GHz VLBI images of OJ 287 observed on May 3, 2014, April 11, 2015, April 22, 2016, and April 16, 2017. Each panel shows the source morphology in total intensity as reconstructed by ehtim. The color scale is the same for all days and is in units of brightness temperature; an angular scale indicator is included in the bottom-right corner of each image (200 μas) for reference. The χ2 statistics used as diagnostics of the model fit quality – the closure phase reduced χ2 (χcphase2) and the logarithmic closure amplitude reduced χ2 (χlogamp2) – are displayed in the bottom-left corner of each image. The ridgelines, delineating the jet’s axis for each epoch, are overlaid as dashed lime lines, with a black outline for enhanced visibility. In this image sequence we can clearly see the jet’s remarkable PA evolution over the four years.

In the text

	Fig. 4. Fits to the variations in the PA of the radio jet at 43 GHz and 86 GHz based on the expected jet precession in the disk model by Dey et al. (2018). The red and black lines correspond to the model predictions at 43 GHz and 86 GHz, respectively. The innermost PA value from RadioAstron is marked with black diamond, and VLBA-BU-BLAZAR data points are marked with red circles.
In the text

Fig. 5. Flux density and radial core separation evolution of jet components within 0.2 mas of the core in OJ 287 at 22, 43, and 86 GHz from September 2015 to April 2017. Top panel: Flux density evolution over time, showing the core (red circles), a knot identified as B15 at 43 GHz and K at 86 GHz (green square and yellow triangle, respectively), and components C1 (magenta diamond) and C2 (purple cross) identified in this work at 22 GHz. Error bars indicate measurement uncertainties. Bottom panel: Radial core separation projected distances. A TeV flare (dashed cyan vertical line) and the centenary optical flare (dashed gray vertical line) are marked. The dashed red vertical line and the gray shaded region show the estimated ejection time and its uncertainty for knot B15 or K (Weaver et al. 2022).

In the text

	Fig. 6. Comparison of the rotation rate of the polarization angle rotation event at 10.45 GHz shown in Myserlis et al. (2018) and the rotation rate of the helical jet shape along the jet direction as seen in the 2016 RadioAstron observations of OJ 287. In both cases we see a tighter winding of the jet helical, with it bending closer to the jet base.
In the text

	Fig. A.1. Selected closure amplitudes and phases from coherently averaged visibilities on triangles as a function of time.
In the text

Fig. A.2. Self-calibrated visibility amplitudes and phases as a function of uv distance from the RadioAstron observations of OJ 287 from April 24–25, 2016, at 22 GHz. The magenta points represent the fit to the data using the ehtim model derived from the imaging process. Reliable space-ground fringe detections were achieved up to a projected baseline length of 4.6 Earth diameters. The dot-dashed black line separates the ground-only baselines from those involving the SRT.

In the text

	Fig. A.3. Multiplicative gain correction factors at each station.
In the text

	Fig. B.1. Time evolution of selected baseline polarization amplitudes and phases derived from coherently averaged visibility measurements.
In the text

	Fig. B.2. D-term solutions for each station and each polarization. Labels mark the different radio telescopes.
In the text