The Parsec-scale Jet of the Neutrino-emitting Blazar TXS 0506+056

Figure 1 shows the radio jet morphology of TXS 0506+056 on four representative epochs at 15 GHz. A core and four jet components (named as J1 to J4 from the outermost inwards) are fitted. This nomenclature is different from the previous VLBI studies (Kun et al. 2019; Ros et al. 2020) and ensures a uniformity in the description of emitting jet components in this study. The largest detectable jet extent is related to the quality of individual images. The rms noise levels in most images are in the range of 0.1–0.3 mJy beam⁻¹, except the last two images, which show slightly higher noise levels of 0.4 and 0.5 mJy beam⁻¹ (Table A1). In order for a proper comparison, we fixed the lowest contours in Figure 2 to a uniform value of 1 mJy beam⁻¹, which corresponds to two to five times the rms noise (Table A1). The maximal jet length at 15 GHz was observed on 2018 May 31, which we estimate to be about 4.2 mas in projection on the plane of the sky.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The jet shows a prominent bending morphology, as seen in Figure 1. The jet body extends first to the south, then turns to southeast at about 1.5–2 mas away from the core. Two jet components, J2 and J3, are detected at the section of jet bending. The 8 and 5 GHz data indicate a consistent curved jet structure as inferred from the 15 GHz images. The PAs of the innermost jet component J4 show a distinct gradual change from ∼197° (or −163° in Table A2) in 2009 to 178° in 2018. Similar changes of PAs are also seen in J2 and J3. From Table A1, it can be inferred that the restoring beam of each individual image is similar due to the almost identical (u, v) coverage of these observations. The distinct change in the jet structure may not simply originate from the image quality systematics or deconvolution beams but may be genuinely related to the jet dynamical structure.

The sizes of the fitted jet components increase with the radial distance from the center. The outermost component J1 is the most extended with the largest size, and the innermost J4 has the smallest size. Therefore a conical geometry could be a good description of the TXS 0506+056 jet. Since the jet body is curved, the opening angle of the jet may not be appropriately represented by the envelope containing all jet components. Figure 2 shows the jet width (w_r) versus radial distance of the jet component from the core (r) for all epoch data at 15 GHz. The jet width is estimated based on the jet component size from the model fitting of the visibility data with circular Gaussian models (refer to Table A2) rather than from the fitting of the brightness cross section perpendicular to the jet ridgeline in the image. This is as the former can yield reduced errors due to sparse sampling in the u − v space and an increase in the weights assigned to bright structures. The scatter in the data points for each component, mainly along the jet width axis owes to errors arising from the model fitting of their corresponding brightness distributions, and random errors. The jet width at increasing radial distances is well fitted with a linear function, supporting the inference of a conical jet body. The projected half opening angle of the jet beam $\tan 2\phi \,={\rm{\Delta }}{w}_{r}/{\rm{\Delta }}r$ from which ϕ = 110 ± 03.

Figure 3 shows the radial distance of the jet component from the center of the core versus observing time. A linear regression gives jet component proper motion speeds of 0.071 ± 0.018 mas yr⁻¹ (1.49 ± 0.34c) for J1, and 0.019 ± 0.003 mas yr⁻¹ (0.41 ± 0.06c) for J4 (Table 1). It indicates that the apparent jet speed is higher in the outer region (J1: ∼90 pc) than the inner portion (J4: ∼13 pc). Accounting for the conical jet geometry (Section 3.1), this may not suggest genuine jet acceleration at de-projected distances of tens of pc. The fitting of J2 and J3 gives negative values, −0.013 ± 0.007 mas yr⁻¹ (−0.28 ± 0.15c) and −0.017 ± 0.007 mas yr⁻¹ (−0.35 ± 0.15c), respectively. Since J2 and J3 are located at the jet bending location, their motion is largely affected by the helical trajectory. In the following discussion, we will use the speed estimated for J1 as a representative of the jet speed. Compared with the typical γ-ray-emitting blazars (Lister et al. 2016), the jet proper motion speed of TXS 0506+056, even referring to the fastest jet component J1, is considerably low.

Zoom In Zoom Out Reset image size

Table 1.  Vector Motion Fit Properties of Jet Components

Comp.	N	$\left\langle S\right\rangle$	$\left\langle r\right\rangle$	$\left\langle \vartheta \right\rangle$	μapp	βapp
		(mJy)	(mas)	(°)	(mas yr−1)	(c)
(1)	(2)	(3)	(4)	(5)	(6)	(7)
J1	17	20.8	3.669	172.1	0.071 ± 0.016	1.49 ± 0.34
J2	17	46.8	1.865	188.0	0.013 ± 0.007	0.28 ± 0.15
J3	15	11.4	1.265	218.8	0.017 ± 0.007	0.35 ± 0.15
J4	17	98.2	0.524	186.6	0.019 ± 0.003	0.41 ± 0.06

Note. (1) Component name, (2) number of fitted epochs, (3) mean flux density at 15 GHz, (4) mean distance from core, (5) mean PA with respect to the core feature, (6) proper motion, (7) apparent speed in units of the speed of light.

Download table as:  ASCIITypeset image

The east and west boundaries of the jet structure are marked by the components J1 and J3, respectively. A mean PA of the jet PA_mean = 1954 is estimated by averaging the PAs of J1 and J3. For a scenario involving a precessing motion, such as along a helical trajectory, the jet component PA swings around PA_mean. Figure 4 shows the relative jet PAs dPA (defined as $\mathrm{dPA}=\mathrm{PA}-{\mathrm{PA}}_{\mathrm{mean}}$) versus the observing time. Components J2 and J3 exhibit this wiggling behavior, indicated by the oscillations of dPA(J2) and dPA(J3). A sinusoidal fit to dPA versus time gives periodicities of about 5.15 yr for J2 (χ² = 0.5), and 5.23 yr for J3 (χ² = 0.4). The relative PA dPA_J1 ranges between −51 and 24. When fitting a sinusoidal function with J1, it indicates a period of about 17.61 yr (χ² = 2.3), but only covers half a cycle. The increasing cycle period with radial distance in the jet components indicates that they trace helical paths. The projection of a conical helical line with a constant helix angle parallel to the axis of the cone (here, PA = 1954) onto the plane of the sky is a logarithmic spiral showing increasing pitch with radial distance. This is consistent with the increasing trend of the apparent jet speed with radial distance.

Zoom In Zoom Out Reset image size

The component J4 also shows a significant variation, but a periodic pattern is not easily discernible. This is as the model fitting is affected by the brighter core component with J4 in close proximity. Whereas, the dPA range of J4, ±10°, is substantially larger than those of J1, J2, and J3 which are within ∼6°, implying that the innermost cone within 0.5 mas (converted to a deprojected physical scale of ∼14 pc) may have a larger opening angle. Such a geometry with a wider initial opening angle at the jet base has been noticed before (Ros et al. 2020) and has been observed in other AGN (e.g., Giovannini 2018), and can resemble an inverted funnel structure as expected in kinematic models of the innermost jet (e.g., Mohan & Mangalam 2015).

The brightness temperature (T_b) of a jet component can be employed as an indicator of Doppler boosting. The study of Kellermann & Pauliny-Toth (1969) first suggested that a maximal brightness temperature of 10¹¹–10¹² K for a compact self-absorbed radio source may be set by a catastrophic inverse Compton cooling when exceeding 10¹² K. For an incoherent synchrotron radio source, the condition of equipartition of energy between the radiating particles and magnetic fields provides a constraining ${T}_{{\rm{b}},\mathrm{eq}}=5\times {10}^{10}$ K (Readhead 1994). Brightness temperature can be calculated from the VLBI observables (see Equation (4) in Kovalev et al. 2005). The observed brightness temperature ${T}_{{\rm{b}},\mathrm{obs}}$ is amplified compared to the intrinsic ${T}_{{\rm{b}},\mathrm{int}}$ of the synchrotron emission by the Doppler factor δ when the relativistic jet from the AGN is directed very close to the observer's line of sight. Taking ${T}_{{\rm{b}},\mathrm{eq}}$ as an upper limit of ${T}_{{\rm{b}},\mathrm{int}}$, the Doppler factor is $\delta ={T}_{{\rm{b}},\mathrm{obs}}/{T}_{{\rm{b}},\mathrm{int}}$. The ${T}_{{\rm{b}},\mathrm{obs}}$ and δ for the TXS 0506+056 core component are listed in columns 8–9 in Table A2; δ ranges between 1.6 and 5.8 during 2009 January—2016 June (relatively quiescent period) with an average of 2.65.

The bulk Lorentz factor Γ and jet inclination θ can be calculated from the measurement of jet apparent speed and Doppler boosting factor. The jet inclination angle is θ = 20° ± 2° assuming the above β_app and δ = 2.65 and the intrinsic jet half opening angle is 38 ± 05. The bulk Lorentz factor during the quiescent period is Γ_q = 1.9 ± 0.2. After 2016 November, δ rises to a maximum of 13.6 on 2017 June 17 followed by a decrease to 9.0 on 2018 December 16. The averaged Doppler factor during the flaring period is $\left\langle {\delta }_{{\rm{f}}}\right\rangle =9.8$ with a corresponding Γ_f = 5.1 ± 0.1. The relatively lower Doppler factor and slower jet speed in TXS 0506+056 compared to the γ-ray-emitting blazar population (Lister et al. 2016) suggests only a moderately relativistic jet in TXS 0506+056.

The frequency-dependent apparent core position shift in blazars is attributable to synchrotron self-absorption in the radio core (Lobanov 1998; Hirotani 2005); it can be deduced through a frequency-dependent time lag in radio flares in single-dish measurements (if dominated by the core, e.g., Mohan et al. 2015; Agarwal et al. 2017) or by an observed shift in the core position from near-simultaneous multi-band VLBI images (e.g., Pushkarev et al. 2012; Li et al. 2018) and can be used to estimate the distance of the core from the central engine, the magnetic field strength in the parsec-scale jet, jet luminosity, and power (e.g., Lobanov 1998; Zdziarski et al. 2015). Here, we select the multiple frequency VLBI data observed quasi-simultaneously, i.e., the separation of observation time is less than two months to minimize effects of variability, to estimate the core shift and magnetic field (Table 2).

The projected core shift is

$$\begin{eqnarray}&&{\rm{\Delta }}{r}_{\mathrm{proj}}=\displaystyle \frac{{\rm{\Delta }}d}{2\tan {\phi }_{\mathrm{proj}}}=\displaystyle \frac{{\rm{\Delta }}d\sin \theta }{2\tan \phi },\end{eqnarray} \tag{ 1 }$$

where ${\rm{\Delta }}{r}_{\mathrm{proj}}$ is core shift in mas, ${\rm{\Delta }}d={d}_{1}-{d}_{2}$ is core size difference between frequency ν₁ and ν₂ in mas, and ϕ_proj and ϕ are the projected and real jet half opening angle (in degrees) respectively. The following jet parameters are estimated using the formulation prescribed in Lobanov (1998). The core-shift measure

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{r\nu }=4.85\times {10}^{-9}\displaystyle \frac{{D}_{L}{\rm{\Delta }}{r}_{\mathrm{proj}}}{{(1+z)}^{2}({\nu }_{1}^{-1}-{\nu }_{2}^{-1})}\,\mathrm{pc}\,\mathrm{GHz},\end{eqnarray} \tag{ 2 }$$

where z is redshift, D_L is luminosity distance in pc, and ν₁ and ν₂ are observation frequency in GHz. This corresponds to the core distance from the central engine of

$$\begin{eqnarray}&&{r}_{\mathrm{core}}=\displaystyle \frac{{{\rm{\Omega }}}_{r\nu }}{\nu \sin \theta }.\end{eqnarray} \tag{ 3 }$$

The magnetic field strength at a fiducial 1 pc is

$$\begin{eqnarray}&&{B}_{1}\cong {\left(\displaystyle \frac{{{\rm{\Omega }}}_{r\nu }^{3}{(1+z)}^{2}}{{\delta }^{2}\phi {\sin }^{2}\theta }\right)}^{1/4},\end{eqnarray} \tag{ 4 }$$

with a corresponding scaled magnetic field strength in the core ${B}_{\mathrm{core}}={B}_{1}{({r}_{\mathrm{core}}/(1\mathrm{pc}))}^{-1}$.

We use the redshift z = 0.3365, luminosity distance D_L = 1762 Mpc, apparent jet speed β_app = 1.49, Doppler factor in the quiescent period δ_q = 2.65 and during the flare δ_f = 9.82, inclination angle θ = 20°, and jet half opening angles ϕ_proj = 11° and ϕ = 38. The core-shift measure Ω_rν ranges between 12.1 ± 0.1 pc GHz and 36.1 ± 0.3 pc GHz (Table 2). The magnetic field at 1 pc, B₁ ranges between (0.4 ± 0.1) G—(0.9 ± 0.2) G, and that in the core B_core ranges between (0.13 ± 0.04) G—(0.47 ± 0.15) G. The distance of the radio core r_core ranges between 0.8 ± 0.1 pc and 6.9 ± 0.7 pc.

The magnetic field B₁ during the quiescent and flaring periods are 0.69 ± 0.03 G and 0.57 ± 0.06 G, respectively, indicating that a portion of the magnetic energy density may have been converted to the radiative energy output during the flare, confined to within the core. The jet parameters and the estimates from the above equations are tabulated in Table 3.

Table 3.  Magnetic Field Strength Before and During 2017 Flare

Parameter	Values
z	0.3365
DL	1762 Mpc
βapp	1.49 ± 0.34
$\left\langle \delta \right\rangle$	2.65 ± 0.01
$\left\langle {{\rm{\Gamma }}}_{q}\right\rangle$	1.9 ± 0.2
$\left\langle \theta \right\rangle$	20 ± 2°
ϕproj	11.0 ± 03
ϕint	3.8 ± 05
${B}_{1\,\mathrm{pc}}^{\mathrm{before}}$	0.69 ± 0.03 G
${B}_{1\,\mathrm{pc}}^{\mathrm{during}}$	0.57 ± 0.06 G

Note. z: redshift; D_L: luminosity distance; ${\beta }_{{app}}$: the apparent speed in units of the speed of light; $\left\langle \delta \right\rangle$ and $\left\langle {{\rm{\Gamma }}}_{q}\right\rangle$: the mean Doppler factor and bulk Lorentz factor during the quiescent period, respectively; $\left\langle \theta \right\rangle$: the jet inclination angle; ${\phi }_{{proj}}$: the apparent half opening angle; ϕ_int: the intrinsic jet half opening angle; ${B}_{1\,\mathrm{pc}}^{\mathrm{before}}$ and ${B}_{1\,\mathrm{pc}}^{\mathrm{during}}$: magnetic field at 1 pc of actual distance from the jet vertex of quiescent and flaring periods, respectively.

Download table as:  ASCIITypeset image

Figure 5 shows the 15 GHz single-dish radio light curve of TXS 0506+056 based on observations made with the OVRO in the top panel. The compiled data points are between 2008 January and 2019 May (about 11.5 yr). The source shows a relatively low-level activity from 2008 to early 2017, with the total flux density ranging between 0.29 Jy and 0.74 Jy and a mean of 0.48 ± 0.1 Jy. After 2017 September, the source entered a flaring state with the flux density continuing to increase. The flux density is currently larger by a factor of >3 compared to the mean value during 2008–2017. This is the largest flare in the observation period of ∼11.5 yr. The available data does not allow the identification of the exact peak as it continues to increase.

Zoom In Zoom Out Reset image size

The integrated VLBI flux densities (red-colored diamond) at 15 GHz are overlaid onto the OVRO light curve (black-colored dots) in the top panel of Figure 5. It shows that the compact VLBI emission dominates the total flux density and follows the same variation trend as the single-dish light curve, and is consistent with the study of Britzen et al. (2019) where this VLBI core dominance was noted. Quantitative comparison between the VLBI flux density and the closest-epoch single-dish one suggests that only 1%–18% is attributable to extended emission, which is resolved out in VLBI images. The VLBI core accounts for 53%–78% of the integrated VLBI flux density in each individual epoch (see Table A2 and the middle panel of Figure 5) and also shows a significantly increasing trend since 2017. A moderate level of variability is found in jet components J1, J2, and J3. The innermost jet component J4 shows notable variation, with a minimum of 45 mJy on 2013 February 28 and a maximum of 154 mJy on 2016 November 18. This may largely be due to rapid local changes in the inclination angle and the corresponding Doppler beaming factor, such as in a turbulent environment or along a helical trajectory.

As the radio flux density is mainly dominated by the core, the aperiodic variability may be attributed to episodes of injection of energetic particles at the jet base. This can result in increased flux density from relativistically beamed propagating shock fronts (Blandford & Königl 1979). Equating the variability timescale t_var to the synchrotron cooling timescale,

$$\begin{eqnarray}&&\displaystyle \frac{\delta {t}_{\mathrm{var}}}{1+z}=\displaystyle \frac{3{m}_{e}c}{4{\sigma }_{T}{\beta }^{2}\gamma {U}_{B}}\end{eqnarray} \tag{ 5 }$$

where m_e = 9.11 × 10⁻²⁸ g is the electron mass, σ_T = 6.65 × 10⁻²⁵ cm² is the Thomson cross section, β ≈ 1 is the velocity, and γ ≈ 100 is the Lorentz factor of the injected electrons, and ${U}_{B}={B}^{2}/(8\pi )$ is the magnetic field energy density. Using t_var = 240 days obtained as the average flare duration from Table 1 of Kun et al. (2019), the size of the emission region ${\rm{\Delta }}r\leqslant \delta {t}_{\mathrm{var}}/(1+z)=0.3$ pc assuming that δ = δ_q = 2.65 is representative of the true size (during the quiescent period). This corresponds to a distance $r={\rm{\Delta }}r/\sin \phi =4.5$ pc, which is an upper limit but similar to the typical r_core. The magnetic field strength is $B=0.7\,{\rm{G}}\,{\delta }^{-1/2}$, which is 0.43 G during the quiescent period and 0.23 G during the flaring period marking a more stark contrast. These values are nearly similar to the core-shift estimates, and suggest a conversion of magnetic energy density to particle kinetic energy which then helps accelerate the particles upon injection at the jet base.

The helical geometry at larger scales is suggestive of an origin from instabilities in the jet. This has been suggested as a likely mechanism in the blazars 1156+295 (Zhao et al. 2011), 3C 48 (An et al. 2010), and 3C 273 (Perucho et al. 2006) where the Kelvin–Helmholtz instability operational from the subparsec scales causes the helically structured jet. This scenario may likely be occurring in TXS 0506+056 as well, since the jet launching region does not exhibit periodic signatures, though this requires further high-resolution VLBI monitoring observations to confirm.

The neutrino and associated γ-ray flares were discovered in the early rising phase of the flare (IceCube Collaboration et al. 2018a), marked with a vertical blue line in Figure 5. Another period, in which a 3.5σ neutrino event was detected, is indicated with a red-colored shade zone. It is puzzling that the 2017 neutrino flare occurred in the still-rising phase of the radio flux density; after the neutrino event, the flux density keeps rising to a maximum of 1.55 Jy on 2019 April 25. It is possible that the neutrino event and γ-ray flares were associated with the onset of the injection of energetic particles into the optically thick jet base.