Probing the Innermost Regions of AGN Jets and Their Magnetic Fields with RadioAstron. V. Space and Ground Millimeter-VLBI Imaging of OJ 287

The bright blazar OJ 287 was observed with RadioAstron on 2014 April 4–5 (from 12:00 to 03:45 UT, observing code GA030E) at a frequency of 22.236 GHz combining a 15:45 hr session during the perigee of the spacecraft (total RadioAstron on-source time was 6.3 hr). These imaging observations were supplemented by five short (between half an hour and 2 hr long) sessions on 2014 March 9, 10, 16, 27, and April 18 within the RadioAstron AGN survey program (Kovalev et al.2020).

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 Table 1. The long-baseline snapshot sessions were conducted on 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 of the snapshot sessions except the one on 2014 March 27, for which fringes were found only to Effelsberg. Signal-to-noise ratio (S/N) 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 1 hr duration to allow for cooling of the onboard 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 (Bruni et al. 2016), 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 us to minimize the effects of the spacecraft acceleration terms along the considered orbit segment, by re-centering 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 space baselines. This provides a first-order value for fringes delay and rate, which can be later refined with the post-processing 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).

2.1.1. Initial Data Processing and Imaging

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 s, 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.

Zoom In Zoom Out Reset image size

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 the AIPS 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 a priori 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 (Kovalev et al. 2020).

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 Figure 2.

Zoom In Zoom Out Reset image size

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 Figure 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. Second, 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 Section 3.2. This is particularly relevant for the extremely small spatial scales probed during the long-baseline 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 RadioAstron data, they provide very valuable information regarding the brightness temperature at the smallest spatial scales, which is analyzed in Section 3.5.

Zoom In Zoom Out Reset image size

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 Figure 3, where we also plot the long-baseline complex visibilities for comparison (see Section 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 2014 April observations are shown in the right panels of Figure 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 FWHM 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 the imaging process.

Zoom In Zoom Out Reset image size

Calibration of the instrumental polarization and absolute orientation of the polarization vectors is discussed in Appendix A.1.

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 2 months from the RadioAstron observations (see Table 1 and and Table 2). At 15 GHz, a single epoch of fully calibrated visibilities is provided via the publicly available database of the long-term monitoring MOJAVE program ²¹ (Lister et al. 2009; Principal investigator: J. L. Richards, project code S6407B). The data set at 43 GHz is part of the VLBA-BU-BLAZAR program ²² , which includes monthly observations of 38 radio and γ-ray bright AGNs. The 43 GHz visibilities data available online are fully self-calibrated, following the analysis described in Jorstad et al. (2017).

The 86 GHz data were obtained from an observation made on May 25 with the GMVA. The duration of each scan was 8 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 (eight physical IFs) per polarization. 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.

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 Traianou et al. (2020). The imaging and modeling of the multifrequency VLBI data allowed us to characterize the central compact region of the jet down to parsec-scale. The cross-identification of the individual components across epochs and frequencies was based on consistency of their positions, flux densities, and sizes. The uncertainties of the component parameters are estimated using the measured rms noise in the respective images and the S/N of detection of individual components (see Fomalont 1999; Schinzel et al. 2012; Pötzl et al. 2021, for details). All of the modelfit components are described in Table 3 and displayed in Figure 4.

3.1.1. Total Intensity

Figure 4 shows the ground array (left panels) and RadioAstron space VLBI images (right panels) of OJ 287 taken between 2014 April and May, providing a detailed view of the jet at different spatial scales.

The 15 GHz image (upper-left panel in Figure 4) shows a bright VLBI core and a one-sided jet that extends toward the west, followed by a strong jet bending by about 55° toward the southwest. The core region is modeled by a single Gaussian component labeled as C (upper-left frame in Figure 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 Figure 4).

At 22, 43, and 86 GHz images, the modelfit component C1 is considered to be the reference point. In the over-resolved RadioAstron image, the knot C1 can be further resolved into two subcomponents, C1a and C1b (lower-right panel in Figure 4). We note that the sum of the flux densities of C1a and C1b components approximately adds to that of component C1 (see Table 3).

The progressively higher angular resolution obtained with the VLBA-BU-BLAZAR 43 GHz image (middle-left panel in Figure 4), and GMVA 86 GHz image (lower-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 northwest 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 Table 3. The position angle of the innermost component rotates from − 584 ± 04 for component J5 at 15 GHz, to − 184 ±14 for C2 at 43 GHz, and finally − 81 ± 11 for component C2 at 86 GHz.

Our space VLBI images of OJ 287 (right panel in Figure 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 counterclockwise rotation from − 117 ± 10 for C2 to 131 ± 04 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 Figure 4). Except for the outermost component J1, all other jet components are cross-identified 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 Figure 5 and Section 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 Section 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.

Zoom In Zoom Out Reset image size

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):

$$\begin{eqnarray}&&{T}_{b}=1.22\times {10}^{12}\displaystyle \frac{S}{{\theta }_{\mathrm{obs}}^{2}{\nu }^{2}}\left(1+z\right)\,\,\left({\rm{K}}\right),\end{eqnarray} \tag{ 1 }$$

where S is the component flux density in Jy, θ_obs is the size of the emitting region in milliarcseconds, ν is the observing frequency in GHz, and z is the source redshift. For unresolved modelfit components, we set $\theta ={\theta }_{\min }$, where ${\theta }_{\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 Table 3.

3.1.2. Polarization

OJ 287 showed multiwavelength flaring activity during 2014 April–May, coincident with our VLBI imaging campaign, as shown in the light curves of Figure 5. The event started at high energies, with γ-ray emission rising on 2014 February 25 (shadowed area in Figure 5). For the complete γ-ray light curve, see the VLBA-BU-BLAZAR website ²³ . A local optical V-band maximum was reported on March 31 (Ganesh et al. 2014), followed a few days later by a flux density increase by almost ∼270% from the flux density level recorded between January and February, as measured by single-dish 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, 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 POLAMI monitoring program ²⁴ at 86 GHz (Agudo et al. 2018b, 2018a).

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 Section 3.3 and Section 3.6.

The quasi-simultaneous multifrequency observations in 2014 April–May 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 self-absorbed (SSA) emission from a homogeneous spherical region filled with relativistic electron-positron plasma with a power-law energy distribution of emitting particles (Türler et al. 2000):

$$\begin{eqnarray}&&{S}_{\nu }={S}_{{\rm{m}}}{\left(\displaystyle \frac{\nu }{{\nu }_{{\rm{m}}}}\right)}^{{\alpha }_{{\rm{t}}}}\displaystyle \frac{1-\exp [-{\tau }_{{\rm{m}}}{\left(\nu /{\nu }_{{\rm{m}}}\right)}^{\alpha -{\alpha }_{{\rm{t}}}}]}{1-\exp (-{\tau }_{{\rm{m}}})}\,\,\,\,[\mathrm{Jy}],\end{eqnarray} \tag{ 2 }$$

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, ${\tau }_{{\rm{m}}}\sim 3/2\left[{\left(1-\left(8\alpha /3{\alpha }_{{\rm{t}}}\right)\right)}^{1/2}-1\right]$ 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. As the 22 GHz observation with RadioAstron (epoch t₁ = 2014.258, in further discussion) preceded the observations at other frequencies by up to 50 days, we attempt to account for the potential effect of source variability by estimating the 22 GHz flux density at the epoch of May 14 (epoch t₂ = 2014.368, hereafter) situated in the middle of the time period between the VLBI observations at 43 GHz and 86 GHz.

At the epoch t₁ of the RadioAstron observation, a total flux density S_t1 = 2.98 ± 0.04 Jy was measured at the Effelsberg 100-m telescope. For the epoch t₂, we use the Gaussian process predictor with the 21–23 GHz measurements shown in Figure 5 to estimate S_t2 = 4.84 ± 0.07 Jy. Hence the 22 GHz emission of OJ 287 should have increased by about 1.8 Jy between the two epochs.

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 Figure 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₁ and t₂. Its flux density at the epoch t₂ 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 Section 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₁ and t₂ 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.

Zoom In Zoom Out Reset image size

A clear peak in the spectrum is observed only in the component C1, while the component J6 has a rising spectrum, and the component J5 has a falling spectrum. 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 Table 5. The fits are also plotted in Figure 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 synchrotron-dominated emission regime (Marscher 1990), although reaching viable conclusions on this matter requires taking into account the acceleration of the emitting plasma (Lobanov & Zensus 1999).

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):

$$\begin{eqnarray}&&{B}_{\mathrm{ssa}}[{\rm{G}}]={10}^{-5}b(\alpha )\,{\theta }_{{\rm{m}}}^{4}{\nu }_{{\rm{m}}}^{5}{S}_{{\rm{m}}}^{-2}\displaystyle \frac{{\delta }_{{\rm{j}}}}{1+z},\end{eqnarray} \tag{ 3 }$$

where $b\left(\alpha \right)$ 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 milliarcseconds, and δ_j is the Doppler boosting factor.

The turnover parameters, S_m and ν_m, are taken from the spectral fits in Table 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 Table 3. For this purpose, we first use Equation (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 $\sqrt{3}/\sqrt{2\mathrm{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 Table 6.

Table 6. Estimates of Magnetic Field and Doppler Factors in the Jet

	θm	Bssa	Beq	δeq	δeq
	(mas)	(G)	(G)	(e− e+)	(e− p)
C1	0.04	0.002 δj	53 ${\delta }_{{\rm{j}}}^{-8/7}$	103	173
C2	0.02	>0.011 δj	>81 ${\delta }_{{\rm{j}}}^{-8/7}$	<63	<105
		0.027 δj	92 ${\delta }_{{\rm{j}}}^{-8/7}$	44	75
J6	0.03	0.24 δj	39 ${\delta }_{{\rm{j}}}^{-8/7}$	11	18
J5	0.20	0.13 δj	4.1 ${\delta }_{{\rm{j}}}^{-8/7}$	5.0	8.1

Note. Column designations: θ_m—estimated effective diameter of the emitting region; B_ssa—magnetic field strength for synchrotron self-absorbed spectrum; B_eq—strength of the equipartition magnetic field; δ_eq—equipartition Doppler factors calculated for electron-positron (e⁻ e⁺, k_u = 1) and electron-proton (e⁻ p, k_u = 100) flow, with ${\delta }_{\mathrm{eq}}\propto {k}_{{\rm{u}}}^{2/15}$. For the component C2, the second row lists respective values obtained for the assumed spectral turnover point (see Section 3.4).

Download table as:  ASCIITypeset image

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 Equation (3):

$$\begin{eqnarray}&&{B}_{\mathrm{eq}}={10}^{-4}{\left[\displaystyle \frac{(1+{k}_{{\rm{u}}}){c}_{12}\,{\kappa }_{\nu }}{f}\displaystyle \frac{{S}_{{\rm{m}}}\,{\nu }_{{\rm{m}}}}{{\theta }_{{\rm{m}}}^{3}\,{D}_{{\rm{L}},\mathrm{Gpc}}}\displaystyle \frac{{\left(1+z\right)}^{10}}{{\delta }_{{\rm{j}}}^{4}}\right]}^{2/7},\end{eqnarray} \tag{ 4 }$$

where D_L,Gpc is the luminosity distance in Gpc, f is the 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 electron-positron flow), the coefficient c₁₂ is given in Pacholczyk (1970), and κ_ν depends on the low- and high-frequency cutoffs, ${\nu }_{\min }$ and ${\nu }_{\max }$, of the SSA emission. For ${\nu }_{\min }={10}^{7}\,\mathrm{Hz}$ and ${\nu }_{\max }={10}^{13}$ Hz assumed in this paper, the ${\nu }_{\min }$ contribution to κ_ν can be neglected, and ${\kappa }_{\nu }\approx {({\nu }_{\max }/{\nu }_{{\rm{m}}})}^{1+\alpha }/(1+\alpha )$. The resulting estimates of B_eq are shown in Table 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 ${\delta }_{\mathrm{eq}}\propto {({B}_{\mathrm{SSA}}/{B}_{\mathrm{eq}})}^{7/15}$. Using this approach, we calculate δ_eq and present them in Table 6 for electron-positron (k_u = 1) and electron-proton (k_u = 100; see Merten et al. 2017) jets.

The observed proper motions of the jet components in OJ 287 at 15 GHz (Lister et al. 2009, 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 ${\delta }_{{\rm{j}}}\geqslant \sqrt{1+{\beta }_{\mathrm{app}}^{2}}\approx 5\mbox{--}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.

The highest values of T_b estimated for C1 reach ≈ 5 × 10¹² K (see Table 3). This is moderately larger than the inverse-Compton limits of $\left(0.3\mbox{--}1.0\right)\times {10}^{12}$ K (Kellermann & Pauliny-Toth 1969) and substantially above the equipartition limit of ≈ 5 × 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 ${\delta }_{{\rm{j}}}\approx \left(5\mbox{--}15\right)$ discussed Section 3.4 can reconcile the highest estimated brightness temperature values (see Table 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}_{{\rm{b}},\min }$, made directly from the visibility amplitudes (Lobanov 2015). This comparison is presented in Figures 2 and 7, showing that ${T}_{{\rm{b}},\min }$ reaches ≈ 10¹³ K at the longest $\left(u,u\right)$-spacings coming from the auxiliary long-baseline segments of the observations.

Zoom In Zoom Out Reset image size

In Figure 7, ${T}_{{\rm{b}},\min }$ shows a continued increasing trend with progressive larger uv-spacing, with no evidence of reaching a plateau. Thus, even larger ${T}_{{\rm{b}},\min }$ could in principle be expected at longer baselines than those observed here. Gómez et al. (2016) report ${T}_{{\rm{b}},\min }$ values of the order of 10¹³ K at the longest baselines during RadioAstron observations of the jet in BL Lac at 22 GHz. The authors interpret these extremely high ${T}_{{\rm{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 Section 3.2). We should also note that the long-baseline snapshots that provide the largest ${T}_{{\rm{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 Figure 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}_{{\rm{b}},\min }\sim {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.

According to theoretical models, jet launching from SMBHs is electromagnetic in nature through the action of helical magnetic fields, either anchored in the accretion disk or the BH 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 × λ², where λ is the observed wavelength and χ₀ is the intrinsic electric vector position angle of the emitting region. The rotation measure (RM) is expressed by (e.g., Thompson et al. 2017):

$$\begin{eqnarray}&&\mathrm{RM}=\displaystyle \frac{{e}^{3}}{8{\pi }^{2}{\epsilon }_{0}{m}_{e}^{2}{c}^{3}}\int {\nu }_{e}{{\boldsymbol{B}}}_{\parallel }d{\bf{l}},\end{eqnarray} \tag{ 5 }$$

where RM is measured in rad/ cm², ₀ is the permittivity of the vacuum, ν_e is the electron density, B _∥ is the component of the magnetic field that is parallel to our line of sight, and d l 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 a 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 Gómez et al. (2016), 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 Faraday-corrected EVPAs is shown in Figure 8. It should be noted that the images used to compute the RM were obtained at different epochs (see Table 1) during a flaring state of the source (see Section 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.

Zoom In Zoom Out Reset image size

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⁻². This is broadly consistent (in magnitude) with a previous core RM estimation of −367 ± 71 rad m⁻² based on simultaneous VLBI observations between 8 GHz and 15 GHz taken in 2006 April 28 (Hovatta et al. 2012).

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 Figure 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 the northeast to southwest 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 & Payne 1982) and observed previously in a number of sources (e.g., Gabuzda et al. 2004; O'Sullivan & Gabuzda 2009; Hovatta et al. 2012; Gabuzda 2017).

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. Figure 9 shows the obtained D-terms for all of 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.

Zoom In Zoom Out Reset image size

Absolute calibration of the EVPA was obtained through comparison with Very Large Array (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°.

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 IF channels for each antenna as the representative D-terms. Figure 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.

Zoom In Zoom Out Reset image size

Calibration of the absolute EVPA for the GMVA 86 GHz observations was carried out through comparison with the single-dish Institut de Radioastronomie Millimétrique (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, 2018a).

The polarized flux density in this work is computed for all data sets by the standard formula $P=\sqrt{{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 ${\sigma }_{P}=\sqrt{{\left(0.1\times P\right)}^{2}+{{rms}}_{{\rm{P}}}^{2}}$, 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\left(P/S\right)$, and the error on m is given by:

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{m}=\displaystyle \frac{1}{S}\sqrt{{\sigma }_{P}^{2}+{\left(m\times {\sigma }_{S}\right)}^{2}+{\sigma }_{{\rm{D}}-\mathrm{term}}^{2}}.\end{eqnarray} \tag{ A1 }$$

The term σ_D−term represents the systematic polarization calibration error, which is defined as (Roberts et al. 1994; Hovatta et al. 2012):

$$\begin{eqnarray}&&{\sigma }_{{\rm{D}}-\mathrm{term}}={\sigma }_{\mathrm{amp}}\displaystyle \frac{\sqrt{{S}_{\mathrm{total}}^{2}+{\left(0.3\times {S}_{\mathrm{peak}}\right)}^{2}}}{\sqrt{{N}_{\mathrm{ant}}{N}_{\mathrm{IF}}{N}_{\mathrm{scan}}}}\end{eqnarray} \tag{ A2 }$$

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° (Hovatta et al. 2012), whereas for the 43 GHz data, Δm = 1% and Δχ = 5°% (Jorstad et al. 2005).