Multifrequency study of the gamma-ray flaring BL Lacertae object PKS 2233–148 in 2009–2012

ABSTRACT

1 INTRODUCTION

The location of the γ-ray production zone in active galactic nuclei (AGNs) is still an open and actively debated question. The limited angular resolution of γ-ray telescopes means that it is impossible to directly locate the region responsible for the high-energy emission in AGNs. A variety of approaches have been considered to address this problem, and our current understanding is that the regions of γ-ray production may be at different locations in different sources, as is evident from observations. One of the two main competing scenarios is based on the observed rapid γ-ray variability on time-scales of a few hours and suggests that the high-energy emission from blazars is generated on subparsec scales, near the central black hole (e.g. Tavecchio et al. 2010; Yan et al. 2018). Similarly, the observed strength and variability of the absorption of the γ-ray emission in the blazar 3C 454.3 suggests that the γ-ray-emitting zone is within the broad-line region (Bai, Liu & Ma 2009; Poutanen & Stern 2010). The second scenario, in contrast, concludes that the dominant population of γ-ray photons is produced at larger, parsec scales, at distances up to 10–20 pc (Marscher et al. 2010; Agudo et al. 2011; Schinzel et al. 2012; Fuhrmann et al. 2016; Karamanavis et al. 2016), and is based on a joint analysis of data in the γ-ray and radio bands. Kovalev et al. (2009) and Pushkarev, Kovalev & Lister (2010) showed that variability in γ-rays leads that of the 15-GHz radio core on a time-scale of up to a few months. In this paper we are concerned with one particular AGN, the BL Lac object 2233 − 148, which was observed during and after the flare in γ-rays registered in 2010 April by Fermi-LAT.

The structure of the paper is as follows: in Section 2 we describe our and archival observational data and reduction schemes; in Section 3 we discuss our results; and our main conclusions are summarized in Section 4. We use the term ‘core’ to denote the apparent origin of AGN jets, which commonly appears as the brightest feature in very-long-baseline interferometry (VLBI) images of blazars (e.g. Lobanov 1998). The spectral index α is defined as S_ν ∝ ν^α, where S_ν is the observed flux density at frequency ν. All position angles are given in degrees east of north. We adopt a cosmology with Ω_m = 0.27, |$\Omega _\Lambda =0.73$| and H₀ = 71 km s⁻¹ Mpc⁻¹ (Komatsu et al. 2009).

2 OBSERVATIONS AND DATA PROCESSING

2.1 Multi-epoch 4.6–43.2 GHz VLBA observations

For the purposes of our study, we made use of data on the BL Lac object PKS 2233 − 148 as observed (code S2087D) with the very-long-baseline array (VLBA) of the National Radio Astronomy Observatory (NRAO) during four sessions at epochs 2010 May 15, 2010 June 25, 2010 August 1 and 2010 September 9. All 10 VLBA antennas participated in each experiment. The observations were performed in a full polarimetric mode simultaneously in the C, X, U, K and Q frequency bands, which correspond to 6-, 4-, 2-, 1.3- and 0.7-cm wavelengths, respectively (Table 1). Each band was separated into four 8-MHz-wide intermediate frequency channels (IFs), with 16 spectral channels per IF. The signal was recorded with 2-bit sampling at a total recording rate of 256 Mbps with an analogue base-band converter. The data were correlated at the NRAO VLBA Operations Center in Socorro (New Mexico, USA) with an averaging time of 2 s. We split the C and X bands into two subbands (each of 16 MHz width) centred at 4608.5, 5003.5 MHz and 8108.5, 8429.5 MHz, respectively, and in the subsequent analysis the data were processed independently. The U, K and Q bands were not split into subbands, resulting in 32-MHz bandwidths centred at 15365.5, 23804.5 and 43217.5 MHz. The on-source time at each epoch was, in total, about 45 min for the C and X bands, 53 min for the U and K bands, and 83 min the Q band, split into 12 scans distributed over 8 h. The scans were scheduled over a number of different hour angles to maximize the (u, v) plane coverage. The increase of the on-source time with frequency was scheduled with the aim of obtaining comparable image sensitivity at all bands.

The data reduction was performed with the NRAO Astronomical Image Processing System (aips, Greisen 2003) following the standard procedure. The individual IFs for each frequency band were processed separately throughout the data reduction. The antenna gain curves and system temperatures measured during the sessions were used for a priori amplitude calibration. Global gain correction factors for each station for each IF were derived from the results of self-calibration. We applied the significant amplitude scale corrections listed in Table 2 by running the aips task clcor. The phase corrections for station-based residual delays and delay rates were found and applied using the aips task fring in two steps. First, the manual fringe-fitting was run on a short interval on a bright quasar, 3C 454.3 (2251 + 158), in order to determine the relative instrumental phase and residual group delay for each individual IF. Secondly, the global fringe-fitting was run by specifying a point-like source model and a signal-to-noise ratio cutoff of 5 to omit noisy solutions. The fringe-fit solution interval was chosen to be 10, 4, 2, 1.5 and 1 min for the C, X, U, K and Q band, respectively. After fringe-fitting, a complex bandpass calibration was made. The estimated accuracy of the VLBA amplitude calibration in the 5–15 GHz frequency range is about 5 per cent, and at 24–43 GHz is about 10 per cent (see also Kovalev et al. 2005; Sokolovsky et al. 2011).

cleaning (Högbom 1974), phase and amplitude self-calibration (Jennison 1958; Twiss, Carter & Little 1960) and hybrid imaging (Readhead et al. 1980; Schwab 1980; Cornwell & Wilkinson 1981) were performed in the Caltech difmap (Shepherd 1997) package. A point-source model was used as an initial model for the iterative procedure. Final maps were produced by applying a natural weighting of the visibility function. The spanned bandwidth of the IFs in each band is small (<0.2 per cent of the fractional bandwidth in all bands), and thus no spectral correction technique was applied.

In this paper, we present results inferred from the total-intensity images. The polarization calibration and results will be published in a separate paper.

2.2 Multi-epoch 15.4-GHz MOJAVE observations

We also made use of the data at 15.4 GHz from the MOJAVE (Monitoring of Jets in Active Galactic Nuclei With VLBA Experiments) program.¹ The data were obtained at eight more epochs at 15.4 GHz: 2009-12-26, 2010-06-19, 2010-12-24, 2011-09-12, 2012-05-24, 2012-07-12, 2012-12-10 and 2016-09-17. We used the fully calibrated publicly available data. For a more detailed discussion of the data reduction and imaging process schemes, see Lister et al. (2018). The absolute flux density of the observations is accurate within 5 per cent (Lister & Homan 2005; Hovatta et al. 2012).

2.3 15-GHz OVRO observations

We also used public data² for PKS 2233 − 148 observations performed within the Owens Valley Radio Observatory 40-m Telescope monitoring program (Richards et al. 2011). Observations were carried out at 15 GHz in a 3-GHz bandwidth from 2008-10-23 to 2018-02-05 with a typical time sampling of about 4 d. Details of the data reduction and calibration are given in Richards et al. (2011).

2.4 γ-ray Fermi-LAT data

The γ-ray light-curve was generated from data obtained with the LAT (Atwood et al. 2009) onboard the Fermi γ-ray space telescope between 2008-08-09 and 2016-10-17. In the analysis, we used the FermiScienceTools software package³ version v10r0p5 and Pass 8 data. In order to generate the light-curve, we first selected all photons between 100 MeV and 300 GeV within a 10° region of interest (ROI) around the source. In the event selection and analysis we followed the recommendations for Pass 8 data given by the LAT team.⁴

The photon flux over each 7-d bin was calculated using the tool gtlike with instrument response function version P8R2_SOURCE_V6. The source model was generated using the external tool make3FGLxml.py version 01 by selecting all sources within 20° of the target in the 3FGL catalogue (Acero et al. 2015), and also including the Galactic diffuse emission model version gll_iem_v06 and the isotropic diffuse emission model version iso_source_v06. Based on the 3FGL catalogue, the target was modelled with a log-parabola spectrum, defined as |${\rm d}N/{\rm d}E=N_0(E/E_b)^{-(\alpha +\beta \log (E/E_b))}$|⁠. In order to account for the low number of photons in each weekly bin and to reduce the number of free parameters in the fit, we froze the spectral parameters of the target and all other sources in the model to the values reported in 3FGL. For the target, the 3FGL values are α = 2.04, β = 0.09 and E_b = 581.68. In addition, if the source was beyond the 10° ROI or had a test statistic (TS) value (e.g. Mattox et al. 1996) less than 5 in 3FGL, we also froze the flux to the value reported in 3FGL. If the TS of the bin was less than 4 (corresponding to about 2σ) or if the number of predicted photons in that bin was less than 10, we calculated a 95 per cent upper limit of the photon flux (Abdo et al. 2011).

3 RESULTS

3.1 Parsec-scale jet structure

Final naturally weighted VLBA maps of the source brightness distribution at the seven frequencies at each of the four observing epochs are presented in Fig. 1. The source shows a typical parsec-scale AGN morphology of a bright compact core and one-sided jet, which propagates towards the east and is detected up to a distance of about 2 mas at 43 GHz and progressively farther, up to 8 mas, at lower (4.6, 5.0 GHz) frequencies owing to the steep spectrum of the jet emission (see the more detailed discussion in Section 3.6). At the frequency of 8 GHz and higher, the outer jet regions are transversely resolved. The lower-frequency images show a faint emission beyond the core, probably caused by the uncompensated side-lobes owing to the low declination of the source. The images at 8 GHz are the most sensitive, with a typical noise level of about 0.16 mJy beam⁻¹ and a dynamic range of the order of 3000, determined as the ratio of the peak flux density to the rms noise level. The noise level was calculated as the average of rms estimates in three corner quadrants of the image, each of 1/16 of the map size. The fourth quadrant, with a maximum rms, was excluded as being affected by the source structure. The achieved image noise levels are close to the thermal noise estimated for natural data weighting. In Table 3, we summarize the VLBA map parameters.

Figure 1.

Naturally weighted total-intensity contour maps of PKS 2233−148 at four epochs during 2010 at 4.6, 5.0, 8.1, 8.4, 15.4, 23.8 and 43.2 GHz, with a cell size of 0.3, 0.3, 0.2, 0.2, 0.1, 0.06 and 0.03 mas per pixel, respectively. The x- and y-axes are given in mas of relative right ascension and relative declination, respectively. The contours are plotted at increasing powers of |$\sqrt{2}$| starting from the 4 rms level. The full-width at half maximum (FWHM) of the restoring beam is shown as a shaded ellipse in the lower left corner. Note that the scales in the different images are different. The image parameters are listed in Table 3.

Open in new tabDownload slide

Structure modelling of the source brightness distribution was performed with the procedure modelfit in the difmap package by fitting several circular Gaussian components to the calibrated visibility data and minimizing χ² in the spatial frequency plane. We used a minimum number of components (three at lower and four at higher frequencies) that, after being convolved with the restoring beam, adequately reproduce the constructed source morphology. The obtained source models are listed in Table 4 and provide flux densities, positions and sizes of the fitted components. All the positions are given with respect to the core component.

3.2 Radio and γ-ray light-curves

In Fig. 2, we present light-curves of PKS 2233 − 148 based on the Fermi-LAT and OVRO monitoring data, complemented by measurements from the MOJAVE program and our VLBA observations at 15 GHz. The prominent variability at high energies detected during 2010 April and June triggered the four-epoch VLBA multifrequency campaign. The values of the correlated VLBA total flux density are in good agreement with the single-dish OVRO flux density measurements, implying that there is almost no extended emission on kiloparsec scales, as was previously concluded by Drinkwater et al. (1997). We performed a cross-correlation analysis of the light-curves using the z-transformed discrete correlation function (Alexander 1997), specifically developed for sparse, unevenly sampled light-curves. The correlation between the radio and γ-ray light-curves with and without upper limits is insignificant, suggesting that the γ-ray production region in the source might have a complex structure. We discuss this in more detail in Section 3.9.

Figure 2.

(Top) OVRO 15-GHz flux density evolution. Grey rectangles together with red stars indicate the VLBA total flux density at 15.4 GHz, observed within our campaign S2087D and the MOJAVE program. (Bottom)Fermi weekly binned γ-ray light-curve at 0.1–300 GeV. Upper limits are given by blue arrows. Dotted vertical lines indicate the epochs of the multifrequency VLBA observations.

Open in new tabDownload slide

3.3 Core shifts

The VLBI core is believed to represent the apparent jet starting region, located at a distance r_core from the central engine, at which its optical depth reaches τ_ν ≈ 1 at a given frequency. Thus, as a result of nuclear opacity, the absolute position of the radio core is frequency-dependent and varies as |$r_\text{core}\propto \nu ^{-1/k_\text{r}}$| (Blandford & Königl 1979; Konigl 1981). There is growing observational evidence from recent multifrequency studies of the core-shift effect for k_r ≈ 1 (e.g. O’Sullivan & Gabuzda 2009; Fromm et al. 2010; Sokolovsky et al. 2011; Hada et al. 2011; Kravchenko et al. 2016; Lisakov et al. 2017). This is consistent with the Blandford & Königl (1979) model of a synchrotron self-absorbed conical jet in equipartition between energy densities of the magnetic field and the radiating particles. Departures in k_r from unity are also possible (Kutkin et al. 2014, Plavin et al. 2018) and can be caused by pressure and density gradients in the jet or by external absorption from the surrounding medium (Lobanov 1998; Kadler et al. 2004).

We calculated the core-shift vector as |$\boldsymbol {\Delta r}_{\text{core},\nu _1\nu _2} = \boldsymbol {\Delta r}_{12} - (\boldsymbol {r}_1 - \boldsymbol {r}_2)$|⁠, where |$\boldsymbol {\Delta r}_{12}$| is the displacement of the phase centres of the images at different frequencies, and |$\boldsymbol {r}_1$|⁠, |$\boldsymbol {r}_2$| are VLBI core-position offsets from the phase centre. In order to derive the image shift vector |$\boldsymbol {\Delta r}_{12}$|⁠, we used the fast normalized cross-correlation algorithm (Lewis 1995) to align the images to the same position on the sky, selecting the jet regions of optically thin emission and assuming that their positions are achromatic. Every pair of images was restored with the average beam size using a pixel size of 0.03 mas.

In Fig. 3, we present a plot of 65 derived core-shift vectors, where the head of each vector represents the core position at lower frequency, while all core positions at higher frequency of a pair are placed at the origin. The dotted line corresponds the median jet direction of P.A. = 112°. The core-shift effect occurs predominantly along the jet direction. In 68 per cent of cases, the core-shift vectors deviate by less than 10° from the median jet position angle. This good alignment is achievable for a relatively straight jet without substantial curvature in the core region. Assuming that the core shift takes place along the jet and that errors are random in direction, then the standard deviation of the transverse projections of the core-shift vectors on to the jet direction yields a typical error of 45 μas. Thus, in 90 per cent of cases the derived core shifts are significantly (>2σ) different from zero. In Table 5 we list the core-shift measurements: (1) epoch of observations, (2) pair of frequencies, (3) core-shift magnitude, (4) core-shift direction, (5) difference of observing wavelengths.

Figure 3.

Core-shift vectors measured in all frequency pairs. The typical error is 0.045 mas. The shaded grey area encompasses 68 per cent of the vectors deviating by less than 10° from the median jet direction shown by the dotted line.

Open in new tabDownload slide

We studied the frequency dependence of the core shifts (Fig. 4) by fitting a function |$\Delta r_\text{core} = b(\nu ^{-1/k_\text{r}}-\nu _\text{max}^{-1/k_\text{r}})$|⁠, where b and k_r are fitted parameters, and ν_max is fixed to the maximum frequency to which the core shifts were measured (43 GHz for the epochs 2010-05-15, 2010-06-25 and 2010-09-09; 23 GHz for the epoch 2010-08-01, at which we could not reliably measure the core shift with respect to the core position at 43 GHz). The fitted k_r values are smaller then but close to one and not significantly different from it. This can hold even during an outburst. As discussed in Plavin et al. (2018), a flare propagating down the jet disturbs only a limited portion of it, and thus k_r deviates from unity in a limited frequency range, which is significantly narrower than that of our observations. Therefore, for further analysis we use k_r = 1. In this case, r_core ∝ λ. Following the approach proposed in Voitsik et al. (2018), in Fig. 5, we plot the measured core shifts against the difference in observing wavelengths (Table 5) and fit the dependence by a straight line, from which one can also estimate an offset of the apparent jet base at a given wavelength from the true jet origin by setting λ₁ = 0.

Figure 4.

Frequency dependence of core shifts measured relative to the core position at 43 GHz (23 GHz for the epoch 2010-08-01) for all observational multifrequency epochs. Solid lines represent the best power-law fits. Shaded areas show 1σ confidence regions of the fit.

Open in new tabDownload slide

Figure 5.

Core shifts as a function of difference of observing wavelengths. Solid lines represent the best linear fits. Shaded areas show 1σ confidence regions of the fit. Stars denote the expected core shifts from the true jet origin (λ₁ = 0) at the wavelengths of our observations, namely 0.7, 1.3, 2.0, 3.6, 3.7, 6.0 and 6.5 cm.

Open in new tabDownload slide

3.4 Jet shape

The core-shift measurements allow us to perform a jet geometry analysis for the whole set of the fitted components, including the cores. In Fig. 6, we plot the transverse jet widths d as the FWHM of the fitted Gaussian components at all four multifrequency VLBA epochs (Table 4) or the corresponding resolution limits (Kovalev et al. 2005), whichever is larger, as a function of their distance r from the true jet base taking into account the core-shift effect. The respective shifts Δr_core = a(t)λ were added to the fitted core separations, where a(t) is the fitted slope at a corresponding observational epoch (Fig. 5). From this analysis we excluded 16 weak components with SNR < 15 to reduce the influence of low-SNR data points, although the whole set of 96 components yields qualitatively similar result. The BL Lac object PKS 2233 − 148 shows a conical streamline, with d ∝ r^{1.01 ± 0.04} at scales probed by the multifrequency VLBA observations down to 0.1 mas. This is consistent with k_r ≈ 1 derived from the core-shift analysis.

Figure 6.

Jet width versus distance to the jet vertex for 78 components from structure model fits at seven frequencies for four epochs. Cores are marked by larger symbols. The dashed line is the best fit from the least-squares method. The jet shape at scales probed by our multi-frequency VLBA observations is conical.

Open in new tabDownload slide

The apparent jet-opening angle of the source is 38° ± 3°, as reported by Pushkarev et al. (2017), who measured it from a stacked total-intensity image at 15.4 GHz as a result of combining VLBA maps from 11 epochs distributed over a time range of about 3 yr, from 2009 to 2012. The wide opening angle suggests that the jet viewing angle is rather small, of the order of a few degrees.

3.5 Redshift constraint from the jet geometry

The optical spectrum of the source shows no prominent emission lines (Drinkwater et al. 1997) and its redshift is still unknown, but the inferred conical jet shape indicates that this BL Lac object is not too close, and probably located at a redshift exceeding ∼0.1. Otherwise, our VLBA observations would be sensitive enough to reveal a jet geometry transition from parabolic to conical shape, as we detect this transition in a number of nearby (z ≲ 0.1) sources and explain it by a transition from a magnetically dominated to a particle-dominated regime in the outflows (Kovalev et al. 2018, in prep.). More distant AGNs (z ≳ 0.1) typically show close to conical jet streamlines (Pushkarev et al. 2017) because the scales probed by VLBI observations are beyond the shape-transition region.

3.6 Spectral index distribution

The procedure of image registration by means of 2D cross-correlation described in Section 3.3 allows us to align the images at different frequencies and accurately reconstruct a distribution of spectral index α over the source morphology. As an example, in Fig. 7 we present a spectral index map of PKS 2233 − 148 calculated between 4.6 and 23.8 GHz at the epoch 2010-09-09 of our multifrequency VLBA observations. The map shows that the core is partially opaque, with a spectral index of about 0.3, while the outer jet regions are optically thin, with a median value of α_jet = −0.95, which is typical for many other parsec-scale AGN jets (e.g. Pushkarev & Kovalev 2012; Hovatta et al. 2014). The spectral index error map manifests a higher α accuracy in the innermost jet area and progressively larger uncertainties towards regions with lower brightness, where random errors arising from the image noise dominate. Systematic errors (from image alignment) dominate in the core area, especially behind it. The same result was obtained statistically for a large sample of sources in our earlier paper (see Appendix B in Hovatta et al. 2014).

Figure 7.

The distribution of the spectral index in PKS 2233−148 is shown in the top left panelat the epoch 2010 September 9 calculated between 4.6 and 23.8 GHz; it is shown in colour, with the 23.8-GHz total-intensity contours overlaid. The contours are plotted at increasing powers of 2, starting from 0.35 per cent of the peak brightness of 303 mJy beam⁻¹. The white curve denotes the total-intensity ridgeline. The restoring beam is depicted as a shaded ellipse in the lower left corner. The spectral index 1σ error map is shown in the top right panel. The bottom left and right panels show the total-intensity and spectral index profiles along the ridgeline, respectively. The vertical dashed line indicates the edge of the convolved VLBI core along the inner jet direction, while the horizontal dashed line represents the median jet spectral index. The grey area shows 1σ errors on the spectral index.

Open in new tabDownload slide

To analyse how the spectral index changes along the jet, we reconstructed the ridgeline of the outflow in total intensity using a procedure described in Pushkarev et al. (2017). As seen from Fig. 7 (bottom right), the spectral index along the ridgeline slightly flattens in the jet knots, indicating a reacceleration of emitting particles, while between them it is steeper. A similar behaviour was found to be typical in AGN jets (Pushkarev & Kovalev 2012; Hovatta et al. 2014).

The evolution of the spectral index along the ridgeline for all the frequency pairs (except for Q-band data) taken at all four epochs is presented in Fig. 8. Beyond the core region, the spectral index deviates around the value of about −1. The spectral index of optically thin synchrotron radiation parametrizes the energy spectrum of relativistic radiative particles. Assuming a power-law energy distribution N(E) = N₀E^−s, the power index s = 1 − 2α has the mean value of ∼3.0. The evolution of α_jet down the jet shows no effect of spectral aging (steepening downstream) owing to radiative losses of relativistic electrons (Kardashev 1962), which is often seen in AGN jets on parsec scales (Pushkarev & Kovalev 2012). In the case of PKS 2233 − 148, the absence of this effect is probably caused by the dominance of a few quasi-stationary jet features, which could be standing shocks that effectively accelerate the emitting particles.

Figure 8.

Spectral index evolution along the ridgeline for α-maps restored at all frequency pairs except the 43-GHz data. The grey bars show 1σ errors.

Open in new tabDownload slide

3.7 Synchrotron spectrum fitting and magnetic field estimates

Figure 9.

Spectral fits to the core (top) and jet feature J2 (bottom) data. Solid lines represent the spectra derived from the homogeneous synchrotron source model. Dashed lines show a simple power-law model. The best-fit parameters of the models are shown on each plot.

Open in new tabDownload slide

The mean value of the magnetic field inferred from equation (2) for the apparent core at the turnover frequency is (0.03 ± 0.01)δ/(1 + z) G.

The only jet component detected at all the frequencies is J2 (see Fig. 6 and Table 4). It is located at a distance of about 2 mas from the core at 43 GHz. The spectra of this jet feature J2 in addition to a synchrotron fit were also fitted by a simple power law (Fig. 9, bottom). The spectrum of J2 is steep, with a spectral index gradually decreasing from about −1 to −1.6, while the turnover frequency increases slightly from 4.9 ± 0.4  to 5.9 ± 0.3 GHz.

3.8 Evolution of the turnover frequency and source kinematics

The self-absorption turnover frequency derived from the core spectra (Fig. 9, top) gradually decreases from 16.8 ± 5.3 GHz on 2010 May 15 to 6.4 ± 1.1 GHz on 2010 September 9, following an inverse proportionality to time (ν_m ∝ t⁻¹), as predicted by a model of a conical jet with constant plasma speed (Blandford 1990). We interpret these changes as direct observational evidence of a flare propagating downstream. As a result of synchrotron opacity in the nuclear region, the flare developing along the jet becomes detectable at progressively larger distances from the true jet origin, corresponding to the VLBI core locations r_core at longer wavelengths. After the disturbance crosses the apparent core (τ ≈ 1 zone) at a given frequency, its flux density starts to decrease, resulting in a steepening of the core spectrum (Fig. 9, top) owing to energy losses to synchrotron radiation or Compton scattering (Kardashev 1962; Marscher & Gear 1985). The core offset from the jet apex can be calculated as r_core = r_flare = a(t)λ_m(t), where the parameter a(t) was derived from the core-shift analysis for each epoch of the multifrequency VLBA observations (Fig. 5).

In Fig. 10, we show r_flare(λ_m) as a function of time. The slope of the weighted linear fit is 1.17 ± 0.10 mas yr⁻¹. The derived proper motion of the flare propagation is significantly higher than that inferred from a kinematics analysis based on tracing bright jet components. The two jet knots, J2 and J3, of the source studied within the MOJAVE program at 15 GHz are slow pattern features, with an angular speed of −71 ± 35 μas (apparent inward motion) and 45 ± 41 μas (Lister et al. 2016), respectively. These components are quasi-stationary and can be standing recollimation shocks observed in sources with super-magnetosonic jets (e.g. Asada & Nakamura 2012; Cohen et al. 2014) and also obtained in numerical 2D relativistic (magneto)-hydrodynamic simulations (Mizuno et al. 2015; Fromm 2015; Fuentes et al. 2018).

Figure 10.

Flare propagation down the jet. The measurements denote distances from the jet origin to the VLBI core at the frequency of maximum emission (Fig. 9, top) at the four observing VLBA epochs. The solid line represents the best linear fit. Dark and light shaded areas show 1σ and 2σ confidence regions of the fit, respectively. Dot–dashed lines indicate the epoch 2010.31 of the γ-ray flare and the corresponding distance 0.12 mas of the γ-ray emission zone from the central engine. Dashed lines indicate the expected epoch range (2012) when the flare reaches the jet component J2 at 2-mas separation from the true jet base, assuming a constant flare propagation speed.

Open in new tabDownload slide

Because the jet geometry is found to be conical at scales probed by the VLBA observations (Section 3.4), we assume that the regime of constant flow speed holds at least up to 6 mas from the jet apex. Then the expected epoch for the flare to reach the quasi-stationary jet feature J2 at a distance of about 2 mas from the true jet base is ∼ 2012.0 (Fig. 10), at which the turnover frequency is expected to decrease to about 1.5 GHz. This scenario is supported by the flux-density evolution of the component at 15 GHz (Fig. 11). The component shows an increase of the flux density by a factor of about 2 during the period of time from ∼2011.5 to ∼2012.5. The epoch of the peak around 2012.0 was established by fitting the data with Gaussian process regression performed with the PyMC3 python module for Bayesian modelling, for which we used an exponential quadratic covariance function. Note that moving jet components behave in a completely different manner. Typically, their brightness rapidly fades as a result of energy losses and adiabatic expansion (e.g. Pushkarev & Kovalev 2012; Kravchenko et al. 2016; Lister et al. 2016).

Figure 11.

Flux density at 15 GHz of the jet component J2 located at 2 mas from the jet origin. The significant increase of the flux density in 2012 is probably a result of the flare that originated in 2010 and reached this jet region in 2012, as expected from the assumed constant flare propagation speed. Error bars represent 1σ uncertainties of individual measurements. The smooth solid curve shows the Gaussian process fit. Dark and light filled areas correspond to 68 per cent and 95 per cent confidence intervals of the fit, respectively.

Open in new tabDownload slide

It is therefore possible that the flare propagation rate represents the bulk flow speed. Taking into account the lower limit on redshift (z > 0.49) derived from the spectroscopy of the absorption lines formed by the intervening gas (Sbarufatti et al. 2006), the proper motion of the disturbance μ = 1.17 ± 0.10 mas yr⁻¹ corresponds to an apparent speed of |$\beta _\text{app}\gt 34\pm 2\, c$|⁠. This is much faster than the typical apparent speed of |${\approx }4\, c$| derived from a kinematics analysis for a sample of 42 BL Lacs and also significantly higher than the |$\beta _\text{app}^\text{max}=21\, c$| detected in the high-redshift (z = 1.07) BL Lac object 1514 + 197 (Lister et al. 2016). Similarly, by analysing multifrequency time delays of the flares and measuring the core shifts in the blazars 3C 454.3 (Kutkin et al. 2014) and 0235 + 164 (Kutkin et al. 2018) it was found that this approach yields the source jet speed, which is by a factor of a few higher than the estimates based on kinematic analysis.

In Fig. 12, we plot the derived B₁ and B_core for the 43-, 15- and 5-GHz cores as functions of redshift. Assuming that z > 0.5, the magnetic field at a distance of 1 pc from the central engine is of the order of 1 G. We note that the estimates of B_core at 15 GHz derived from the core-shift analysis are comparable (lower by a factor of a few) to those inferred from the synchrotron spectrum fits (Section 3.7), if the source Doppler factor is moderate (δ ≲ 5), as is often observed to be the case in BL Lacertae objects (Hovatta et al. 2009; Liodakis et al. 2017).

Figure 12.

Magnetic field constraints obtained from the core-shift measurements and the flare propagation speed. The dark grey area shows estimates of the magnetic field at 1-pc distance from the jet origin, while the light grey stripes represent the magnetic field at the core at 43.2, 15.4 and 5.0 GHz.

Open in new tabDownload slide

3.9 Location of the γ-ray emission region and the source of seed photons

In order to estimate the location of the γ-ray emission zone in PKS 2233 − 148 we extrapolated the r_flare dependence (Fig. 10) back to the epoch of the γ-ray flare, 2010.31 (Fig. 2). This yields an angular separation of 0.12 ± 0.03 mas from the true jet base, which corresponds to the VLBI core position at 24 GHz (Fig. 6). On a linear scale,this separation is >0.7 ± 0.2 pc in projection, which exceeds a de-projected separation of 8 ± 2 pc if we assume a jet viewing angle of 5°. Considering the second peak of the γ-ray data at epoch 2010.46 (2010 June 17), we inferred a distance of about 0.3 mas from the jet apex. This distancecorresponds to the innermost jet feature J4 detected at 24 and 43 GHz, setting an absolute distance of about 20 pc for another possible location of the γ-ray emission site. These assessments favour a scenario in which the γ-ray production zone is located at large distances from a central energy generator (beyond the broad-line region or torus) and is probably associated with one or more standing shocks in a relativistic outflow of PKS 2233 − 148, as suggested by the complex structure of the major high-energy flares in the source. A similar conclusion regarding the remotness of the γ-ray emission region in blazars on scales of parsecs from the central black hole was also reached from other arguments in a number of recent single-source studies of 1510 − 089 (Marscher et al. 2010), OJ 287 (Agudo et al. 2011), 3C 345 (Schinzel et al. 2012), CTA 102 (Casadio et al. 2015), 1502 + 106 (Karamanavis et al. 2016), BL Lacertae (Wehrle et al. 2016), and also from statistical results from the F-GAMMA project (Fuhrmann et al. 2016). At the same time, the arguments based on short-scale variability and breaks in GeV spectra discussed in the Introduction indicate that the high-energy production site is in the immediate vicinity of the black hole.

Another noticeable feature of the parsec-scale morphology of the source is the presence of a sheath around the jet, as indicated in the 8- and 15-GHz maps (Fig. 1), which provide a combination of a high angular resolution and sensitivity. In order to better visualize the sheath emission we convolved the 8.1-GHz image at the epoch 2010-05-15 with a circular beam, setting its FWHM to that of the minor axis of the original map (Fig. 13). The fact that this sheath is slower than the central spine might mean that it acts as a source of additional seed photons for the γ-ray radiation (e.g. Marscher et al. 2010; Aleksić et al. 2014). Thus, the high-energy emission of the source can be formed through (i) the synchrotron self-Compton mechanism acting in its relativistic outflow and the upscattering of low-energy synchrotron seed photons; and (ii) external Compton scattering resulting from a photon field in the sheath.

Figure 13.

Total-intensity map of PKS 2233−148 at 8.1 GHz at the epoch 2010-05-15 from Fig. 1 but convolved with a circular beam of the size of the minor axis of the original image. Emission at the jet edges indicates the presence of a distinct boundary layer.

Open in new tabDownload slide

While detailed modelling of the spectral energy distribution (SED) is beyond the scope of this paper, based on the publicly available non-simultaneous data in the SSDC SED builder tool⁵ it seems that, when the source is in a high state in γ-rays, the luminosity of the inverse Compton peak is higher than that of the synchrotron peak, indicating that an additional external photon field is indeed needed. However, this should be verified with simultaneous data in all bands, taken at both low and high activity states of the source.

We can also speculate that if instead of the epoch of the γ-ray peak we consider the epoch at which the flare starts rising (a few weeks before the peak), then the γ-ray emission site could be in the immediate vicinity of the central machine. This scenario is vulnerable, however, as a plasma cloud moving fast down the jet leaves the seed photon area rapidly, while the flare is still reaching its maximum.

4 SUMMARY

We performed a radio and γ-ray joint study of the BL Lacertae object PKS 2233 − 148, using multiwavelength data in the period 2009–2012. The VLBA observations at 4.6–43.2 GHz reveal the core-dominated, one-sided and relatively straight jet morphology of the source extending up to 8 mas at a position angle of 112°. By analysing jet widths derived from the structure model fits we established that the outflow has a conical shape. This sets a lower limit of about 0.1 on the source unknown redshift.

We measured the frequency-dependent shift vectors of the apparent core position using a method based on results from (i) structure model fitting and (ii) image alignment achieved by implementing a 2D cross-correlation technique on the optically thin jet regions. The magnitude of the core shifts ranges from 0.04 to 0.7 mas, with a typical uncertainty of 45 μas. The directions of the shift vectors are predominantly aligned with the median jet position angle, deviating from it by ≲ 10° in 68 per cent of cases. The derived core shifts show a frequency dependence |$\propto \nu ^{-1/k_\text{r}}$|⁠, with k_r ≈ 1 indicating that nuclear opacity is dominated by synchrotron self-absorption, and that the physical conditions in the jet on scales probed by the VLBA observations are close to equipartition. We did not find evidence for significant changes in k_r between the observing epochs covering a time-scale of four months, during which a flare was developing down the jet. This suggests that the transverse size of the disturbance area is significantly smaller than the jet part constrained by the magnitude of the core-shift effect within a frequency range of 5–43 GHz. The VLBI core position r_core as a function of wavelength follows an |$r_\text{mas}^\text{core}\approx 0.1\lambda _\text{cm}$| dependence. The magnetic field at a distance of 1 pc from the jet apex derived from the core shift measurements is about 1 G.

We have presented a method of independent assessment of jet kinematics based on core-shift measurements and the evolution of the synchrotron spectrum of the VLBI core. The turnover frequency of the core spectrum linearly shifts towards lower frequencies with time, as the flare originating in 2010 April in γ-rays propagates down the jet. The speed of this propagation is about 1.2 mas yr⁻¹ and probably represents the bulk flow speed. It is much higher than results from traditional kinematics based on tracking bright jet features, namely 0.045 mas yr⁻¹ (Lister et al. 2016).

We have found indications that the γ-ray production zone in the source is located at large distances, 10–20 pc, from a central engine, and can be associated with the stationary radio-emitting jet features observed with VLBI. This favours synchrotron self-Compton scattering as the dominant high-energy radiation mechanism in the relativistic jet of the source. Direct observational evidence for a boundary layer around the jet suggests that the sheath might be an additional source of seed photons for external Compton scattering acting in the source.

ACKNOWLEDGEMENTS

We would like to thank the anonymous referee as well as E. Ros for useful comments and suggestions. The VLBA data processing and core-shift analysis were supported by the Russian Science Foundation grant 16-12-10481. The radio/γ-ray joint analysis was supported by the Academy of Finland projects 296010 and 318431. TH acknowledges support from the Turku Collegium of Science and Medicine. This research has made use of data from the MOJAVE data base, which is maintained by the MOJAVE team (Lister et al. 2018). The MOJAVE project was supported by NASA-Fermi GI grants NNX08AV67G, NNX12A087G and NNX15AU76G. This work made use of the Swinburne University of Technology software correlator (Deller et al. 2011), developed as part of the Australian Major National Research Facilities Programme and operated under licence. This research has made use of data from the OVRO 40-m monitoring program (Richards et al. 2011), which is supported in part by NASA grants NNX08AW31G, NNX11A043G and NNX14AQ89G and NSF grants AST-0808050 and AST-1109911. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

Footnotes

REFERENCES

APPENDIX A: MAGNETIC FIELD FROM SYNCHROTRON SELF-ABSORPTION

The interpretation of a radio spectrum with a low-frequency turnover caused by synchrotron self-absorption, and the determination of physical parameters within this assumption date back to the 1960s (see, for example, one of the pioneering works Slish 1963). In particular, the magnetic field associated with a source of synchrotron emission can be inferred. However, the approximate values of the numerical coefficient in the formula that are a function of a spectral index α (S_ν ∝ ν^α) of the optically thin part of a synchrotron spectrum are tabulated for a limited number of α values, ranging from −0.25 to −1.0 (Marscher 1983). The relation for this coefficient was previously discussed by Gould (1979). In this Appendix, we derive a formula for this coefficient that can be computed precisely.

Figure A1.

Coefficient b(α) and optical depth τ_m(α) as a function of the spectral index α. Dashed lines represent approximate solutions derived using equation (A12).

Open in new tabDownload slide

In Fig. A1, we plot the coefficient b as a function of α = (1 − s)/2. The departure of the approximate solutions from the exact ones exceeds 5 per cent for α ≲ −0.82.