EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 695 (March 2025) A&A, 695 (2025) A103 Full HTML
Open Access
Issue
A&A
Volume 695, March 2025
Article Number A103
Number of page(s) 9
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/202452530
Published online 11 March 2025

© The Authors 2025

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model.

Open Access funding provided by Max Planck Society.

1. Introduction

The first active galactic nucleus (AGN) that could be linked to the emission of an extremely high-energy (EHE) neutrino was TXS 0506+056 (IceCube Collaboration 2018a) in 2017. The identification was possible thanks to the multiwavelength flaring of TXS 0506+056 (from radio to TeV), observable by many ground- and space-based telescopes (e.g., Padovani et al. 2019). In addition, an analysis of archival IceCube data revealed evidence of enhanced neutrino activity between September 2014 and March 2015 (IceCube Collaboration 2018b, 2024). Most recently, further neutrino emission was observed from the direction of TXS 0506+056 on April 18, 2021 (2021.30) by the Baikal Gigaton Volume Detector (Baikal-GVD; (Allakhverdyan et al. 2024)) and on September 18, 2022 (2022.72) by IceCube (e.g., Becker Tjus et al. 2022).

TXS 0506+056 is a BL Lac object at a redshift of z = 0.3365 ± 0.0010 (Paiano et al. 2018). In Britzen et al. (2019a,b) we proposed that the enhanced neutrino activity in 2014–2015 and the EHE neutrino were generated in a collision that took place within the precessing jet. We discussed two scenarios: a collision with another jet in a supermassive binary black hole system and a collision between newly emitted and older, strongly curved jet material within one jet. Further very long-baseline interferometry (VLBI) studies of the parsec-scale jet structure were performed by, for example, Li et al. (2022), Ros et al. (2020), and Kun et al. (2019).

Interestingly, the radio light curve recently showed an increasing flux trend that is not seen in other wavelengths (e.g., Acciari et al. 2022). With more Very Long Baseline Array (VLBA) data now available, we re-sorted the nuclear radio structure to reconsider possible explanations for the atypical and still enigmatic parsec-scale kinematics and their possible connection to neutrino emission.

Throughout the paper we adopt the following parameters: a luminosity distance DL = 1.8277 Gpc at the source redshift of z = 0.3365 with cosmological parameters corresponding to a Λ cold dark matter Universe with Ωm = 0.308, Ωλ = 0.692, and H0 = 67.8 km s−1 Mpc−1 (Planck Collaboration XIII 2016). Thus, a proper motion of 1 mas yr−1 corresponds to an apparent superluminal speed of 21.64c, while 1 mas = 4.961 pc.

2. Observations and data analysis

2.1. VLBA data analysis

We remodeled and reanalyzed 12 VLBA observations (15 GHz, MOJAVE1) performed between 2018.96 and 2022.85. In addition, 12 VLBA observations (8 GHz; web site of the Astrogeo Center2), observed between 2010.81 and 2018.84, were remodeled and reanalyzed. The 8 GHz and 15 GHz maps with model-fit components superimposed are available online in Figs.B.1–B.6. In the following, we discuss the results obtained by combining these new datasets with those analyzed and published in Britzen et al. (2019a).

Gaussian circular components were fitted to the data to obtain the optimum set of parameters within the difmap-modelfit program (Shepherd et al. 1997). The average beam full width at half maximum at 8/15 GHz is 1.94 × 0.80 or 1.06 × 0.48 mas, respectively. According to experience, the typical angular resolution is one fifth of the beam, but it increases with signal-to-noise ratio (Lee et al. 2008). Thus, the 8/15 GHz observations have a typical angular resolution of 0.39 × 0.16 mas or 0.21 × 0.10 mas, respectively. The median image sensitivity at 8/15 GHz is 0.15 and 0.54 mJy/beam, respectively.

Every epoch was fitted independently from all the other epochs. The model-fitting procedure was performed blindly so as not to impose any specific outcome. Since the unique neutrino detection from this source, we expected TXS 0506+056 to reveal special properties. Therefore, and similar to our approach for the first publication on this source (Britzen et al. 2019a), without knowing what these peculiar properties would be, we also remodeled very faint components. With this approach, any unusual morphology and morphological changes should appear in the model-fitting process. The excellent data quality of the Monitoring Of Jets in Active galactic nuclei with VLBA Experiments (MOJAVE) program allows such a deep analysis. Special care was taken to correctly identify the core component in every individual dataset.

In Fig. 1 we show, as an example, one of the maps with the Gaussian components superimposed. The brighter features are clearly found in a central north-south-directed stripe.

thumbnail Fig. 1.

Image obtained from data observed on June 11, 2022, shortly before the recording of the fourth neutrino event. Gaussian model-fit components are superimposed.

To obtain the conservative uncertainty estimates of the model component parameters, we followed Britzen et al. (2023), who employed the uncertainty estimates of Lee et al. (2008), but we calibrated them using several sources from the MOJAVE sample that showed a tight power-law dependence of the brightness temperature (Tb) on the component size (R), accounting for the residual uncertainty in the amplitude scale at each epoch. The resulting uncertainty estimates are conservative (i.e., upper limits) because some dispersion in the Tb(R) dependence could have an intrinsic origin. This is similar to the method of Lister et al. (2009) for estimating the positional uncertainty of the model components from the fit of the multi-epoch kinematic data. The obtained uncertainties are typically larger that those estimated in Britzen et al. (2019a) following the bootstrap approach. To account for the uncertainty in amplitude scale left after a self-calibration we added a 5% error to our flux errors in quadrature (Kovalev et al. 2005; Lister et al. 2018).

2.2. Optical data analysis

The optical observations were conducted in the R band using the Apogee Ap6E CCD camera (1024 × 1024, 24 × 24 μm, quantum efficiencies are 0.40 and 0.72 at 400 and 560 nm, respectively) attached to the prime focus of the 70 cm meniscus telescope (f/3) at Abastumani Observatory. The readout, digitizing, downloading time is 1 s. We only used the central portion of 350 × 350 sq. pixels (15 × 15 sq. arcmin), while the entire field of view is 40 × 40 sq. arcmin. The exposure times are from 30 to 180 s, depending on the activity state of the source. The images were reduced using Daophot II and calibrated using three comparison stars in the field of the source.

2.3. Fermi-LAT data analysis

We obtained the Pass-8 photon data for TXS 0506+056 from the Fermi Science Support Center for the time interval MJD 54683-60487 (August 5, 2009, until June 26, 2024). We created a model file using the script make4FGLxml. In the most recent Fermi observations with the Large Area Telescope (LAT) source catalog (Ballet et al. 2023), TXS 0506+056 has the spectral shape of a log-parabola,

d N d E = N 0 ( E E b ) − ( α + β ln ( E / E b ) ) , $$ \begin{aligned} \frac{\mathrm{d}N}{\mathrm{d}E} = N_0\left(\frac{E}{E_{\rm b}}\right)^{-(\alpha + \beta \ln (E/E_{\rm b}))}, \end{aligned} $$(1)

with normalization N0, spectral parameters α and β, and scale parameter Eb. To generate a light curve, we carried out an unbinned likelihood analysis to time bins of width 45 days, using the FermiTools. For this fit, the normalization for all sources within a radius of 5° around our source of interest was left as a free parameter, while all other parameters were fixed to their catalog values. We chose the energy range 0.1–300 GeV. The resulting light curve is shown in Fig. 7 b, where we only show data with a test statistic TS ≥ 25, corresponding to a detection threshold of 5σ. This cut led to the rejection of three data points.

3. Results

The analysis of the new VLBA data again reveals an unexpected pc-scale jet morphology of TXS 0506+056. The distribution of jet components is puzzling and does not correspond to typical blazar jet morphologies. In the following, we discuss our findings in more detail.

3.1. Inner jet components

In Britzen et al. (2019b) we identified jet components across the epochs. For these inner jet components, we identified a1, 1a1, 1a2. In Fig. 4a we again plot these features (in greenish colors).

In addition, we show the inner jet components that could be identified based on the more recent 15 GHz data. 1a3, 1a4, and 1a5 are newly identified. For the older components, we observe that their paths are displaced with time (from right to left) and components move away from the core. The “new” jet components (in magenta, red, and orange) show a significantly different behavior with time. 1a3 and 1a5 move on the same, rather straight path, which is perpendicular to the general jet direction (Fig. 4a). 1a5 starts when 1a3 is no longer seen. The direction of motion for 1a3 and 1a5 is the same: from right to left (indicated by an arrow in magenta and an arrow in orange, respectively). However, 1a5 starts farther to the right, as seen in Fig. 4a. 1a3 is seen for the last time (2020.45) when 1a5 appears.

Interestingly, the last five components (between 2011.16 and 2019.96) of the earlier observed component 1a trace a similar path as 1a3 and 1a5. We speculate that neither the path of 1a3 nor the path of 1a5, nor the last part of the path of 1a, are component paths in a traditional sense. All three might be related to lensing.

It is remarkable that we can perform such detailed studies of a small region that is within 1 mas of the core along the y-axis, and within 0.3 mas along the x-axis. This confirms our hypothesis that this region might be gravitationally lensed and thus appears enlarged.

The apparent speeds, based on the linear regression for those newly identified components shown in Fig. 4b are 2.43 ± 0.65c (0.15 ± 0.04 mas/yr) for 1a3, 4.86 ± 0.97c (0.30 ± 0.06 mas/yr) for 1a4, and 0.81 ± 0.17c (0.05 ± 0.01 mas/yr) for 1a5. We note that these apparent speeds refer to the radial motion only.

The flux-density evolution of these jet features is shown in Fig. 4c. There is a rise in flux for 1a3. Interestingly, when 1a5 becomes visible on the right, the flux is the same as when 1a3 stopped to be visible. 1a4 reveals the highest flux of all jet features and reveals a similar drop in flux as 1a5.

3.2. Ordering of jet components into a ring structure after 2016

The distribution of xy-positions of the jet components changes drastically with time. In Fig. 3a we show the component positions before 2016 and in Fig. 3b we show the positions after 2016. Figure 3a reveals an extended jet morphology with a diagonal structure that Britzen et al. (2019a) interpreted as possibly due to a collision with a second jet.

thumbnail Fig. 2.

All radio component positions in xy-coordinates. Light blue symbols denote the 15 GHz data from Britzen et al. (2019a), and dark blue and orange the new 15 GHz and 8 GHz results obtained in this study. In panel a the symbols are of equal size, while in panel b they are scaled according to the flux density (log). A zoomed-in view of the central region is shown, and all components below a flux limit of 1 Jy are displayed. For better visibility, the core flux has been removed. In this and all subsequent figures, the core position is marked with a black cross. The data of the two plots are available in TablesC.1 and C.2.

thumbnail Fig. 3.

xy-positions before 2016 (a) and after 2016 (b). The two plots are shown with the same scale to facilitate their comparison. An ordering of component positions seems to take place after 2016 and is visible as a ring structure around the radio core in panel b. The data are available in Tables D.1, D.2, and E.1.

thumbnail Fig. 4.

Evolution of the inner jet component properties. (a) Inner jet component paths in xy-coordinates. 1a, 1a1, and 1a2 (green shades) were identified in Britzen et al. (2019b). 1a3, 1a4, and 1a5 are identified in the 15 GHz data from this paper. The most recent data point (cyan) is from epoch 2022.85 and lies close to the path that seems to also be followed by 1a3 and 1a5. 1a3 starts where 1a5 ends; this position is marked by an open black square. The directions of motion for 1a3 and 1a5 are marked by arrows in magenta and orange, respectively. The radio core is indicated by a black filled circle (not to scale). (b) Evolution of the distance from the core for the features shown above (with linear regression as a solid line). (c) Flux densities of those inner jet components that could be traced across the epochs.

thumbnail Fig. 5.

Jet structures as recovered by model fitting. Panel a shows five epochs, and the potentially unaltered jet of TXS 0506+056 is visible. Panels b and c show selected epochs with peculiar, but similar, arc-like structures. Superimposed lines do not represent fits to the data but accentuate the detected radio features.

The jet from 2016 onward (Fig. 3b) is – in contrast – dominated by a ring-like structure (with a radius of about 1.5–2 mas) around the core region. Radio emission is also detected on the counter-jet side (at distances to almost 2 mas from the core). We thus observe a drastic change of jet component positions with time, which cannot be reconciled with the typically observed outward directed motion of jet components in blazars.

3.3. Arc-like structures in the jet

While the map in Fig. 1 shows a radio structure at 15 GHz that resembles typical AGN jets, a more detailed look at individual epochs reveals a so-called classical jet in only a subset of the analyzed epochs (Fig. 5a). Most of the epochs show structures that resemble arcs (see Figs. 5b, and c). We find that the jet features tend to rearrange rapidly to new configurations with time.

3.4. Jet structures at times of neutrino detections

To explore whether the rapidly varying jet morphology is related with the neutrino detections, we show in Figs. 6a–d the jet features at those epochs close to the reported neutrino detections. Obviously, the number and arrangement of the jet features changes significantly between the detections.

thumbnail Fig. 6.

Positions of the radio features model-fitted in the jet of TXS 0506+056 at the epochs closest in time to the neutrino detections.

The broadest distribution with most jet features is detected during the time of the neutrino flare (Fig. 6a). At this time, the jet appears like two jets crossed. It seems that the jet emission is enhanced and widely distributed. At the time around IceCube-170922A (Fig. 6b), the distribution of jet features appears focused around a central point at 0.1, −2.5 mas in xy-coordinates.

GVD210418CA (Fig. 6c) and IceCube-220918 (Fig. 6d) were at times detected with similar jet structures. In addition to the probably pure jet (in the north-south direction), a few other features indicate a ring structure. At these times, the jet and core are less affected by the observed rearrangements that took place at the times of the other neutrino emissions.

3.5. Atypical long-term radio flaring and dip in the optical data

In Fig. 7a we show the decomposed total flux density obtained in the VLBA observations. A major flare of the total flux density, as well as of the core flux density, started in 2016 and continued beyond 2022. A double peak at the maximum of the total radio flux-density flaring is observed. The core flux reveals a fast beating with three maxima at the peak. Atypical is that this major flare is not accompanied by γ-ray flaring, which can be seen when looking at the γ-ray light curve in Fig. 7b. This finding is confirmed by, for example, Allakhverdyan et al. (2024).

thumbnail Fig. 7.

Radio and Fermi-LAT γ-ray fluxes. (a) Total radio flux density (light green), the core flux density (black), and the jet flux density (light blue). The red symbols denote the flux of component 1a4. Vertical solid brown lines indicate the dates of the neutrino detections. (b) Fermi-LAT γ-ray fluxes (red points). Epochs of neutrino detections are marked by the three vertical black lines. The time of increased neutrino flux is highlighted by the gray-shaded area. The tabulated data are available online in Tables C.1 and C.2.

Figure 7a also shows that the 15 GHz jet flux density (light blue) appeared rather constant until 2020. However, a sharp increase in flux density (by a factor of about five) was observed at the time of the first minimum at the peak of the core flux. Jet feature 1a4 is seen in the jet at the same time (at a core distance of about ∼0.3 mas) with a high flux density.

The optical light curve, shown in Fig. 8, reveals a dip at 2019.8 – 2020.1, shortly before the peak in the radio core flux density. The dip in the optical data is accompanied by short flares at the beginning and the end of the dip. The optical data reveal a flare within the time of the beating of the radio core (October 18, 2020).

thumbnail Fig. 8.

Optical light curve as measured by the Abastumani Observatory.

4. Discussion

The origin of astrophysical very high-energy neutrinos detected by IceCube still remains unsolved, although several potential source classes have been proposed (for recent reviews see, e.g., Halzen et al. 2023; Murase et al. 2023). The jets of radio-loud AGNs are among the best candidates with evidence of a significant spatial correlation between IceCube neutrino arrival directions and radio-loud and γ-ray bright AGN monitoring (e.g., Plavin et al. 2020, 2021, 2023; Buson et al. 2022).

TXS 0506+056, based on its now four detected neutrino events, seems to be one of the best sources to study this phenomenon in more detail. The high resolution of VLBI observations seems well suited for a thorough investigation of structural peculiarities of the jet, which may be associated with neutrino production in blazars.

We reported evidence of jet precession and a potential interaction: either between the jet of TXS 0506+056 and another jet in a binary black hole system, or a collision between newly emitted and older, bent jet material within the TXS 0506+056 jet at the time before the EHE neutrino event IceCube-170922A (Britzen et al. 2019a). Precession of the jet in TXS 0506+056 has also been modeled by Li et al. (2022) and de Bruijn et al. (2020).

In the meantime, the number of VLBA observations of TXS 0506+056 has increased considerably and allows us to perform a more thorough study of the nuclear jet morphology and its evolution with time. We again find signatures of an atypical jet structure and kinematics in TXS 0506+056 that cannot easily be reconciled with typical AGN jet properties.

4.1. Possible evidence of neutrino production via jet interaction in or with the broad line region

The dominant structure seen in the jet component positions (see Figs. 2a and b for the composite plot of older and more recent data) is the arc-structure at about 1.5–2 mas core distance. As a possible physical origin of the arc-structure we discussed a gravitational lensing scenario (Sect. 4.2). Another possible explanation might be that observed jet features are decelerated or stopped by some kind of interaction at a radial distance of about 1.5–2 mas from the core. The result is that the jet appears shorter at certain epochs and jet components appear at large position angle differences along a bow-like structure. If the arc indeed represents an interaction site, then it is likely that the jet components in the inner jet (within 1 mas core distance) have a different physical nature compared to the jet features that make up the emission at the interaction site (1.5–2 mas distance from the core).

According to Padovani et al. (2019), TXS 0506+056 is not a true BL Lac object, but intrinsically a flat-spectrum radio quasar with hidden broad lines and a standard accretion disk. The existence of a dusty torus is also implied by this new classification as quasar. The authors estimate the radius of the broad line region (BLR) as RBLR ≈ 7 × 1016 cm, using the scaling relation of Ghisellini & Tavecchio (2008), and assume that the BLR is a spherical emitting shell extending from 0.8 to 1 RBLR. Based on variability and causality arguments and assuming a conical jet structure, they constrain the distance of the γ-ray emission region from the central black hole to be up to ∼1018 cm, depending on the Doppler factor.

We next considered a scenario in which the interaction of the jet with the BLR material leads to the acceleration of ultra-relativistic particles, followed by the emission of electromagnetic radiation and neutrinos. We assumed that this acceleration is sufficiently efficient to accelerate protons to the required energies of Ep ≫ 1 PeV to produce neutrinos in the TeV to PeV regime. We postulate that photo-pion production commences at the onset of the interaction, while the emission region is still inside the BLR, and the BLR radiation field acts as the most intense target photon field for pγ interactions.

Padovani et al. (2019) estimate the luminosity of the BLR as LBLR ∼ 8 × 1043 erg s−1. With RBLR given in the previous subsection, this yields a quasi-isotropic energy density of BLR photons of uBLR ∼ 0.043 erg cm−3. The BLR radiation field is expected to be strongly dominated by the Lyα emission line at ELyα = 10.2 eV, corresponding to a normalized energy ϵLyα ≡ ELyα/(mec2)∼2 × 10−5. A simple δ-function approximation for the γγ cross section yields a rough estimate for the γγ optical depth to photons around a characteristic energy of Eγ ∼ 25 (1 + z)−1 GeV of τγγ ∼ (σT/3) (uBLR/ELyα) RBLR ∼ 40. Thus, we do not expect high-energy γ-rays produced co-spatially with the neutrino emission to be able to escape from within the BLR. This is consistent with the estimates by Padovani et al. (2019) of the location of the γ-ray emission region well outside the BLR.

Photo-pion production on the Lyα photons from the BLR (Doppler boosted by ∼Γ ≡ 10 Γ1) in the Δ+ resonance at 1232 MeV requires protons of energies E′p ∼ 15 Γ1−1 PeV in the emission-region rest frame, leading to neutrinos (each neutrino attaining ∼5% of the original proton energy) of observed energy Eν ∼ 7.5 δ1 Γ1−1 PeV, assuming that neutrinos, like photons isotropically produced in the emission region, are Doppler boosted by a Doppler factor δ ≡ 10 δ1. The large variety of propagation directions of individual jet components that we have observed seems to suggest that at least some of the components may be on misaligned trajectories, with Γ ≫ δ, plausibly leading to neutrinos in the ∼100 TeV regime.

Reimer et al. (2019, 2020) derived an estimate of the required jet power to produce a neutrino flux corresponding to the neutrino flare from TXS 0506+056 in 2014 – 2015. Using their Eq. (6) (corrected in the erratum), with u′t ≈ Γ2 uBLR ∼ 4.3 Γ12 erg cm−3 and ϵ′t ∼ 2 × 10−4 Γ1 yields a normalization of a power-law proton spectrum N p ( γ p ) = N 0 γ p − 2 $ N_{\mathrm{p}} (\gamma_{\mathrm{p}}) = N_0 \, \gamma_{\mathrm{p}}^{-2} $ of N0 ∼ 6.0 × 1052 δ1−4 Γ1−1. Assuming that the same power-law proton spectrum extends down to γp, min = 1, we find a required jet power of L p ∼ 2.0 × 10 46 R 16 − 1 δ 1 − 4 Γ 1 $ L_{\mathrm{p}} \sim 2.0 \times 10^{46} \, R_{16}^{-1} \, \delta_1^{-4} \, \Gamma_1 $ erg s−1, where R16 is the radius of the assumed spherical emission region in units of 1016 cm. The black-hole mass of TXS 0506+056 was given in Padovani et al. (2019) as MBH ∼ 3 × 108 M⊙, corresponding to an Eddington luminosity of LEdd ≈ 3.8 × 1046 erg s−1. Thus, the proposed scenario of pγ-induced neutrino production on the BLR radiation field seems plausible, with reasonable power requirements.

4.2. Gravitational lensing of the jet

The drastic morphological changes of the jet described earlier might correspond to different phases of gravitational lensing. We summarize and categorize the major transformations observed.

  • 2009–2011: The nuclear jet appears as two crossed jets (Fig. 3a).

  • 2011–2016: The jet appears segmented into arc-like strings. Instead of jet knots moving along a jet ridgeline, we detect arc-like structures (Fig. 5). Between 2014 and 2015, the first segment of the ring starts to form in the southwest direction (Fig. 6a).

  • 2016–2019: Emission on the counter-jet side is detected as well, which gives the impression of a full ring around the radio core when looking at all the positions of radio features model-fitted (Fig. 3b). The long rise of the total VLBA flux density is accompanied by a similar rise of the core flux density, which oscillates at the peak (Fig. 7a). The jet flux density abruptly rises at the time of the peak. Shortly before the peak in the radio, the optical data reveal a dip (Fig. 8).

  • 2019–2022.85: Either the original jet is seen or arc-like strings (Figs. 5a and c; Figs. 6c and d).

We discuss the possibility of gravitational lensing and its implications in the following. Strong gravitational lensing of a distant object results in the appearance of multiple source images and/or elongated arc-like structures (Narayan & Bartelmann 1996).

Today, ∼102 gravitationally lensed quasars have been found (e.g., Barnacka 2018), and their number is expected to grow significantly with future surveys. Although just two gravitationally lensed blazars, namely B 0218+357 and PKS 1830-211, are known to date (see O’Dea et al. 1992; Jauncey et al. 1991; Larchenkova et al. 2011a,b; Barnacka 2018, among others), they represent a unique class of objects that offer an opportunity to reveal the jet structure at high energies (Barnacka et al. 2014). The relation between the γ-ray emission and the radio/submillimeter emission in these two lensed blazars has been studied by, for example, Cheung et al. (2014) and Spingola et al. (2016).

Time delay measurements between images in lensed blazars can be used to provide additional constraints to a lens model and to determine the Hubble parameter, H0 (e.g., Biggs 2018; Muller et al. 2020). Strong lensing has also been used to study variability in blazars. By analyzing Atacama Large Millimeter/submillimeter Array (ALMA) full polarization observations of the lensed blazar PKS 1830−211 Martí-Vidal et al. (2015) and Marti-Vidal et al. (2020) confirm models of magnetic turbulence to explain the observed polarization variability.

Usually, a galaxy cluster or an individual galaxy acts as a lens. In the optical data, there is no indication of the presence of a galaxy close to the line of sight toward TXS 0506+056. Hence, a lens should be “invisible,” with a black hole being the most “natural” candidate. We checked whether gravitational lensing on a black hole can qualitatively reproduce the observed geometry of TXS 0506+056. The lensing black hole is assumed to be quiescent, acting only gravitationally. In the lensing scenario, the arc-shaped structure clearly seen in Fig. 3b could be actually an image of the jet lensed at a black hole. The jet is directed almost toward the observer and the lens lies close to the line of sight.

The characteristic lensing scale – the Einstein angle – is defined in angular units (i.e., in radians) as

Θ E = 4 G M c 2 D ls D l D s , $$ \begin{aligned} \Theta _{\rm E}=\sqrt{\frac{ 4GM }{c^2}\frac{D_{\rm ls}}{D_{\rm l} D_{\rm s}}} ,\end{aligned} $$(2)

where Dls, Dl, and Ds are the angular diameter distances between the lens and the source, the observer and the source, and the observer and the lens, respectively, M is the mass of the lens, and G is the gravitational constant. For a point mass lens, if the lens and the source are exactly aligned, the image will consist of a ring with radius ΘE. If the angular separation between the source and the lens is less than ΘE, then two images of the source are produced whose relative brightness depends on the misalignment between the source and the lens (see Fig. 9 for an illustration).

thumbnail Fig. 9.

Schematic diagram illustrating gravitational lensing of a distant source. Light from a source is bent as it travels toward an observer. As a result, for a point mass lens, an observer sees two images of the source. In the sky plane (right panel), images appear distorted. If the source, the lens, and the observer are positioned in a straight line, then the source will appear as a ring around the lens in the sky plane. In the case of any misalignment, the observer will see two elongated images instead. The left and right panels show the same lensing system but seen at different viewing angles. The drawing is not to scale and is for illustrative purposes only.

At the source redshift z = 0.3365, the angular diameter distance to the source is Ds = DL/(1 + z)2 = 1023 Mpc, where DL = 1.8277 Gpc is the corresponding luminosity distance. If we assume that the observed arc at a distance of 1.5 mas from the core (see Fig. 3b) is an image of the jet, then by solving the equation ΘE = 1.5 mas, we get a relation between the lensing black hole mass and its distance from the observer. At intermediate distances, the required black hole mass is on the order of 105–106 solar masses. As mentioned above, there are no intervening galaxies in the optical data, so the lens should be an isolated non-accreting supermassive black hole (SMBH) without a luminous host galaxy. Even though such objects could in principle exist (van Dokkum et al. 2023), there have been no reliable detections so far.

We next considered another possibility, that the lens is close to TXS 0506+056, that is, TXS 0506+056 is actually a binary or dual system with one active SMBH and a second non-emitting black hole acting as the lens. Then, the distance between the lens and the source Dls ≪ Dl, Dls ≪ Ds. To obtain an arc at a distance of 1.5 mas from the core, we need a lens with a mass

M [ M ⊙ ] ≃ 2.9 × 10 11 D ls [ kpc ] , $$ \begin{aligned} M\,[M_{\odot }] \simeq \frac{2.9\times 10^{11}}{D_{\rm ls}\,[\mathrm{kpc}]}, \end{aligned} $$(3)

that is, M ∼ 1010 M⊙ for a line-of-sight separation of Dls ∼ 10 kpc between the source and the lens is a possible solution. We stress here that the lens mass estimate should not be confused with the BH mass in TXS 0506+056. Gravitational lensing allows one to reconstruct a mass distribution in a lens but not in a source. For “realistic” binary systems with a parsec-scale (or even sub-parsec, De Rosa et al. 2019) projected separation, the Einstein angle is on the order of 1 μas or less (i.e., well below the resolution limit). Dual black hole systems with kiloparsec-scale separation are possible and detected (e.g., De Rosa et al. 2023). Catalogs of dual SMBHs (e.g., Comerford & Greene 2014; Stemo et al. 2021) also include systems that consist of an AGN and a quiescent SMBH. In cosmological simulations, a population of “wandering” SMBHs is naturally produced by the hierarchical merger-driven galaxy assembly (Tremmel et al. 2018; Ricarte et al. 2021b, 2021b). Ricarte et al. (2021b) also demonstrated that wandering SMBHs can be detected even in galaxies without obvious signatures of an ongoing merger. Thus, we conclude that lensing of a jet on a non-emitting SMBH might be responsible for the observed geometry of TXS 0506+056.

Gravitational lensing causes not only distortions in an appearance of a background source, it also amplifies an observed flux from the source. If the ring-like structure observed in TXS 0506+056 after 2016 is indeed due to gravitational lensing, then one should observe an increase in flux from the source. The exact amplification (and also duration in time) depends strongly on the relative position between the source and the lens, their line-of-sight distance, a mass and a velocity of the lens. By measuring the distance to the observed ring-like structure and equating it to the Einstein angle, we get only an approximate relation between the mass and the distance to the lens. A flux magnification and a characteristic timescale of the change in flux remain unconstrained. Moreover, TXS 0506+056 demonstrates different behavior at different wavelengths.

From Fig. 7, we clearly see a long rise (from 2016 to 2020) in the radio emission from the source accompanied by a short outburst in 2017 in γ-rays, while the optical light curve does not show any increase in brightness in the considered period of time. An absence of an expected flux magnification at all wavelengths may be interpreted as an argument against the gravitational lensing scenario, since the lensing is frequency-independent. Although, it should be noted that light curves of TXS 0506+056 at different wavelengths are collected within regions of extremely different areas. High spatial resolution in radio allows one to measure flux density of individual jet components, while the optical light curve comes from the arcminute2-scale area (see Sect. 2.2), and the γ-ray photons are typically gathered within a circle with radius of several degrees. If γ-rays are produced at a distance exceeding Θe = 1.5 mas (in projection) from the core, then γ-rays are not expected to experience strong lensing. Moreover, lensing is achromatic only when emitting regions at different wavelengths coincide in space and have equal spatial extent (the statement stays true for neutrino and gravitational waves as well). Examples of chromatic (micro)lensing in lensed quasars and supernovae can be found in, for instance, Wambsganss (2006), Goldstein et al. (2018), and Vernardos et al. (2024). Detailed modeling of the lens system reproducing the exact relative positions of images, their flux ratio, and additional data on the time delay between images are needed to provide further constraints to the lens system.

The above discussed lensing scenario implies that a lens and the source are roughly co-aligned, and two elongated images appear. One image is the arc-like structure 1.5 mas away from the core and the other image is the emission on the counter-jet side (see Fig. 3b, green dots with y > 0).

We also considered another possibility, that inner jet components (see Fig. A.1) and the arc-like structure are actually both the images resulting from gravitational lensing. The arc-like structure covers a much larger area in the sky plane than the inner jet; thus, it is expected to be much more magnified than the inner jet. But according to Fig. 2b, the situation is just the opposite. Thus, with gravitational lensing we cannot explain the observed “unusual” direction of motion of inner jet components.

5. Conclusions

We have explored the available VLBA observations at 15 GHz (and 8 GHz) to study whether the neutrino-emitting blazar TXS 0506+056 reveals peculiarities that might help us identify differences with non-neutrino-emitting blazars. The jet structure reveals drastic changes, and the jet kinematics is atypical and not comparable to any other blazar jets studied so far.

The data analyzed here are consistent with a neutrino origin in an interaction of the jet with the BLR material. However, we also find several pieces of evidence that, when put together, indicate a strong gravitational lensing of TXS 0506+056 by another SMBH.

Further support for the strong gravitational lensing of IceCube AGN neutrino candidate emitters comes from a study of IceCube AGN neutrino candidate PKS 1717+177. In this blazar, the peculiar meandering of the jet might be due to gravitational lensing by an intervening SMBH (Britzen et al. 2024).

Data availability

Figure A.1 as well as all model-fit files of 8 GHz and 15 GHz data (Figs. B.1–B.3 and B.4–B.6) are available online on Zenodo at https://doi.org/10.5281/zenodo.14724643. All the data displayed in Figs. 2, 7a, and 7b are tabulated (Tables C.1, C.2, D.1, D.2, and E.1) and also available on Zenodo.

Acknowledgments

We thank the anonymous referee for valuable comments which helped to improve the manuscript. We are thankful for the helpful and supportive comments by the internal reviewer S. v. Fellenberg. The authors thank F. Halzen, A. Witzel, J. Becker Tjus, and X. Rodrigues for very inspiring discussions. E.K. thanks the Alexander von Humboldt Foundation for its Postdoctoral Fellowship. E. Kun acknowledges support from the German Science Foundation DFG, via the Collaborative Research Center SFB1491: Cosmic Interacting Matters – from Source to Signal. This work has made use of public Fermi data obtained from the High Energy Astrophysics Science Archive Research Center (HEASARC), provided by NASA Goddard Space Flight Center. This research has made use of data from the MOJAVE database, which is maintained by the MOJAVE team (Lister et al. 2018). The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. We acknowledge The Astrogeo VLBI FITS image database. The use of Very Long Baseline Array (VLBA) data is acknowledged. VLBA is operated by the National Radio Astronomy Observatory, which is a facility of the National Science Foundation, and operated under cooperative agreement by Associated Universities, Inc. The use of the VLBA under the US Naval Observatory’s time allocation is acknowledged. This work supports USNO’s ongoing research into the celestial reference frame and geodesy. These data were retrieved from the VLBA public archive. FJ acknowledges funding by the Austrian Science Fund (FWF) [P35920]. The computational results presented have been achieved in part using the Vienna Scientific Cluster (VSC). MZ was supported in part by the Czech Science Foundation Junior Star grant no. GM24-10599M. LIC, AC and FCP acknowledges funding by ESA PRODEX DPLISAM and by the Romanian Ministry of Research, Innovation and Digitalization under the Romanian National Core Program LAPLAS VII – contract no. 30N/2023.

References

  1. Acciari, V. A., Aniello, T., Ansoldi, S., et al. 2022, ApJ, 927, 197 [NASA ADS] [CrossRef] [Google Scholar]
  2. Allakhverdyan, V. A., Avrorin, A. D., Avrorin, A. V., et al. 2024, MNRAS, 527, 8784 [Google Scholar]
  3. Ballet, J., Bruel, P., Burnett, T. H., Lott, B., & The Fermi-LAT Collaboration 2023, ArXiv e-prints [arXiv:2307.12546] [Google Scholar]
  4. Barnacka, A. 2018, Phys. Rep., 778, 1 [CrossRef] [Google Scholar]
  5. Barnacka, A., Geller, M. J., Dell’Antonio, I. P., & Benbow, W. 2014, ApJ, 788, 139 [NASA ADS] [CrossRef] [Google Scholar]
  6. Becker Tjus, J., Jaroschewski, I., Ghorbanietemad, A., et al. 2022, ApJ, 941, L25 [NASA ADS] [Google Scholar]
  7. Biggs, A. D. 2018, MNRAS, 481, 1000 [NASA ADS] [CrossRef] [Google Scholar]
  8. Britzen, S., Fendt, C., Böttcher, M., et al. 2019a, A&A, 630, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Britzen, S., Fendt, C., Böttcher, M., et al. 2019b, A&A, 632, C3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Britzen, S., Krishna, G., Kun, E., et al. 2023, Universe, 9 [Google Scholar]
  11. Britzen, S., Kovačević, A. B., Zajaček, M., et al. 2024, MNRAS, 535, 2742 [NASA ADS] [CrossRef] [Google Scholar]
  12. Buson, S., Tramacere, A., Pfeiffer, L., et al. 2022, ApJ, 933, L43 [NASA ADS] [CrossRef] [Google Scholar]
  13. Cheung, C. C., Larsson, S., Scargle, J. D., et al. 2014, ApJ, 782, L14 [NASA ADS] [CrossRef] [Google Scholar]
  14. Comerford, J. M., & Greene, J. E. 2014, ApJ, 789, 112 [NASA ADS] [CrossRef] [Google Scholar]
  15. de Bruijn, O., Bartos, I., Biermann, P. L., & Tjus, J. B. 2020, ApJ, 905, L13 [NASA ADS] [CrossRef] [Google Scholar]
  16. De Rosa, A., Vignali, C., Bogdanović, T., et al. 2019, New Astron. Rev., 86, 101525 [Google Scholar]
  17. De Rosa, A., Vignali, C., Severgnini, P., et al. 2023, MNRAS, 519, 5149 [NASA ADS] [CrossRef] [Google Scholar]
  18. Ghisellini, G., & Tavecchio, F. 2008, MNRAS, 387, 1669 [NASA ADS] [CrossRef] [Google Scholar]
  19. Goldstein, D. A., Nugent, P. E., Kasen, D. N., & Collett, T. E. 2018, ApJ, 855, 22 [NASA ADS] [CrossRef] [Google Scholar]
  20. Halzen, F., & Kheirandish, A. 2023, in The Encyclopedia of Cosmology. Set 2: Frontiers in Cosmology. Volume 2: Neutrino Physics and Astrophysics, ed. F. W. Stecker, 107 [NASA ADS] [CrossRef] [Google Scholar]
  21. IceCube Collaboration (Abbasi, R., et al.) 2024, 38th International Cosmic Ray Conference, 1465 [Google Scholar]
  22. IceCube Collaboration (Aartsen, M. G., et al.) 2018a, Science, 361, eaat1378 [NASA ADS] [Google Scholar]
  23. IceCube Collaboration (Aartsen, M. G., et al.) 2018b, Science, 361, 147 [NASA ADS] [Google Scholar]
  24. Jauncey, D. L., Reynolds, J. E., Tzioumis, A. K., et al. 1991, Nature, 352, 132 [Google Scholar]
  25. Kovalev, Y. Y., Kellermann, K. I., Lister, M. L., et al. 2005, AJ, 130, 2473 [Google Scholar]
  26. Kun, E., Biermann, P. L., & Gergely, L. Á. 2019, MNRAS, 483, L42 [NASA ADS] [CrossRef] [Google Scholar]
  27. Larchenkova, T. I., Lutovinov, A. A., & Lyskova, N. S. 2011a, Astronomy Letters, 37, 233 [NASA ADS] [CrossRef] [Google Scholar]
  28. Larchenkova, T. I., Lyskova, N. S., & Lutovinov, A. A. 2011b, Astronomy Letters, 37, 441 [NASA ADS] [CrossRef] [Google Scholar]
  29. Lee, S.-S., Lobanov, A. P., Krichbaum, T. P., et al. 2008, AJ, 136, 159 [Google Scholar]
  30. Li, Z., Qiao, J., Zhao, W., & Er, X. 2022, Journal of Cosmology and Astroparticle Physics, 2022, 095 [CrossRef] [Google Scholar]
  31. Lister, M. L., Cohen, M. H., Homan, D. C., et al. 2009, AJ, 138, 1874 [NASA ADS] [CrossRef] [Google Scholar]
  32. Lister, M. L., Aller, M. F., Aller, H. D., et al. 2018, ApJS, 234, 12 [CrossRef] [Google Scholar]
  33. Martí-Vidal, I., Muller, S., Vlemmings, W., Horellou, C., & Aalto, S. 2015, Science, 348, 311 [Google Scholar]
  34. Marti-Vidal, I., Muller, S., Mus, A., et al. 2020, A&A, 638, L13 [CrossRef] [EDP Sciences] [Google Scholar]
  35. Muller, S., Jaswanth, S., Horellou, C., & Martí-Vidal, I. 2020, A&A, 641, L2 [EDP Sciences] [Google Scholar]
  36. Murase, K., & Stecker, F. W. 2023, in The Encyclopedia of Cosmology. Set 2: Frontiers in Cosmology. Volume 2: Neutrino Physics and Astrophysics, ed. F. W. Stecker, 483 [NASA ADS] [CrossRef] [Google Scholar]
  37. Narayan, R., & Bartelmann, M. 1996, ArXiv e-prints [arXiv:astro-ph/9606001] [Google Scholar]
  38. O’Dea, C. P., Baum, S. A., Stanghellini, C., et al. 1992, AJ, 104, 1320 [CrossRef] [Google Scholar]
  39. Padovani, P., Oikonomou, F., Petropoulou, M., Giommi, P., & Resconi, E. 2019, MNRAS, 484, L104 [NASA ADS] [CrossRef] [Google Scholar]
  40. Paiano, S., Falomo, R., Treves, A., & Scarpa, R. 2018, ApJ, 854, L32 [NASA ADS] [CrossRef] [Google Scholar]
  41. Planck Collaboration XIII 2016, A&A, 594, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. Plavin, A., Kovalev, Y. Y., Kovalev, Y. A., & Troitsky, S. 2020, ApJ, 894, 101 [Google Scholar]
  43. Plavin, A. V., Kovalev, Y. Y., Kovalev, Y. A., & Troitsky, S. V. 2021, ApJ, 908, 157 [Google Scholar]
  44. Plavin, A. V., Kovalev, Y. Y., Kovalev, Y. A., & Troitsky, S. V. 2023, MNRAS, 523, 1799 [NASA ADS] [CrossRef] [Google Scholar]
  45. Reimer, A., Böttcher, M., & Buson, S. 2019, ApJ, 881, 46 [Google Scholar]
  46. Reimer, A., Böttcher, M., & Buson, S. 2020, ApJ, 899, 168 [NASA ADS] [CrossRef] [Google Scholar]
  47. Ricarte, A., Tremmel, M., Natarajan, P., & Quinn, T. 2021b, ApJ, 916, L18 [NASA ADS] [CrossRef] [Google Scholar]
  48. Ricarte, A., Tremmel, M., Natarajan, P., Zimmer, C., & Quinn, T. 2021b, MNRAS, 503, 6098 [CrossRef] [Google Scholar]
  49. Ros, E., Kadler, M., Perucho, M., et al. 2020, A&A, 633, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  50. Shepherd, M. C. 1997, in Astronomical Data Analysis Software and Systems VI, eds. G. Hunt, & H. Payne, ASP Conf. Ser., 125, 77 [Google Scholar]
  51. Spingola, C., Dallacasa, D., Orienti, M., et al. 2016, MNRAS, 457, 2263 [NASA ADS] [CrossRef] [Google Scholar]
  52. Stemo, A., Comerford, J. M., Barrows, R. S., et al. 2021, ApJ, 923, 36 [NASA ADS] [CrossRef] [Google Scholar]
  53. Tremmel, M., Governato, F., Volonteri, M., Pontzen, A., & Quinn, T. R. 2018, ApJ, 857, L22 [Google Scholar]
  54. van Dokkum, P., Pasha, I., Buzzo, M. L., et al. 2023, ApJ, 946, L50 [NASA ADS] [CrossRef] [Google Scholar]
  55. Vernardos, G., Sluse, D., Pooley, D., et al. 2024, Space Sci. Rev., 220, 14 [NASA ADS] [CrossRef] [Google Scholar]
  56. Wambsganss, J. 2006, in Saas-Fee Advanced Course 33: Gravitational Lensing: Strong, Weak and Micro, eds. G. Meylan, P. Jetzer, P. North, et al., 453 [Google Scholar]

All Figures

thumbnail Fig. 1.

Image obtained from data observed on June 11, 2022, shortly before the recording of the fourth neutrino event. Gaussian model-fit components are superimposed.
In the text
thumbnail Fig. 2.

All radio component positions in xy-coordinates. Light blue symbols denote the 15 GHz data from Britzen et al. (2019a), and dark blue and orange the new 15 GHz and 8 GHz results obtained in this study. In panel a the symbols are of equal size, while in panel b they are scaled according to the flux density (log). A zoomed-in view of the central region is shown, and all components below a flux limit of 1 Jy are displayed. For better visibility, the core flux has been removed. In this and all subsequent figures, the core position is marked with a black cross. The data of the two plots are available in TablesC.1 and C.2.
In the text
thumbnail Fig. 3.

xy-positions before 2016 (a) and after 2016 (b). The two plots are shown with the same scale to facilitate their comparison. An ordering of component positions seems to take place after 2016 and is visible as a ring structure around the radio core in panel b. The data are available in Tables D.1, D.2, and E.1.
In the text
thumbnail Fig. 4.

Evolution of the inner jet component properties. (a) Inner jet component paths in xy-coordinates. 1a, 1a1, and 1a2 (green shades) were identified in Britzen et al. (2019b). 1a3, 1a4, and 1a5 are identified in the 15 GHz data from this paper. The most recent data point (cyan) is from epoch 2022.85 and lies close to the path that seems to also be followed by 1a3 and 1a5. 1a3 starts where 1a5 ends; this position is marked by an open black square. The directions of motion for 1a3 and 1a5 are marked by arrows in magenta and orange, respectively. The radio core is indicated by a black filled circle (not to scale). (b) Evolution of the distance from the core for the features shown above (with linear regression as a solid line). (c) Flux densities of those inner jet components that could be traced across the epochs.
In the text
thumbnail Fig. 5.

Jet structures as recovered by model fitting. Panel a shows five epochs, and the potentially unaltered jet of TXS 0506+056 is visible. Panels b and c show selected epochs with peculiar, but similar, arc-like structures. Superimposed lines do not represent fits to the data but accentuate the detected radio features.
In the text
thumbnail Fig. 6.

Positions of the radio features model-fitted in the jet of TXS 0506+056 at the epochs closest in time to the neutrino detections.
In the text
thumbnail Fig. 7.

Radio and Fermi-LAT γ-ray fluxes. (a) Total radio flux density (light green), the core flux density (black), and the jet flux density (light blue). The red symbols denote the flux of component 1a4. Vertical solid brown lines indicate the dates of the neutrino detections. (b) Fermi-LAT γ-ray fluxes (red points). Epochs of neutrino detections are marked by the three vertical black lines. The time of increased neutrino flux is highlighted by the gray-shaded area. The tabulated data are available online in Tables C.1 and C.2.
In the text
thumbnail Fig. 8.

Optical light curve as measured by the Abastumani Observatory.
In the text
thumbnail Fig. 9.

Schematic diagram illustrating gravitational lensing of a distant source. Light from a source is bent as it travels toward an observer. As a result, for a point mass lens, an observer sees two images of the source. In the sky plane (right panel), images appear distorted. If the source, the lens, and the observer are positioned in a straight line, then the source will appear as a ring around the lens in the sky plane. In the case of any misalignment, the observer will see two elongated images instead. The left and right panels show the same lensing system but seen at different viewing angles. The drawing is not to scale and is for illustrative purposes only.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.