A cosmic collider: Was the IceCube neutrino generated in a precessing jet-jet interaction in TXS 0506+056?

Context. The neutrino event IceCube−170922A appears to originate from the BL Lac object TXS 0506+056. To understand the neutrino creation process and to localize the emission site, we studied the radio images of the jet at 15 GHz. Aims. Other BL Lac objects show properties similar to those of TXS 0506+056, such as multiwavelength variability or a curved jet. However, to date only TXS 0506+056 has been identified as neutrino emitter. The aim of this paper is to determine what makes the parsec-scale jet of TXS 0506+056 specific in this respect. Methods. We reanalyzed and remodeled 16 VLBA 15 GHz observations between 2009 and 2018. We thoroughly examined the jet kinematics and flux-density evolution of individual jet components during the time of enhanced neutrino activity between September 2014 and March 2015, and in particular before and after the neutrino event. Results. Our results suggest that the jet is very strongly curved and most likely observable under a special viewing angle of close to zero. We may observe the interaction between jet features that cross each other’s paths. We find subsequent flux-density flaring of six components passing the likely collision site. In addition, we find a strong indication for precession of the inner jet, and model a precession period of about 10 yr via the Lense-Thirring effect. We discuss an alternative scenario, which is the interpretation of observing the signature of two jets within TXS 0506+056, again hinting toward a collision of jetted material. We essentially suggest that the neutrino emission may result from the interaction of jetted material in combination with a special viewing angle and jet precession. Conclusions. We propose that the enhanced neutrino activity during the neutrino flare in 2014–2015 and the single EHE neutrino IceCube-170922A could have been generated by a cosmic collision within TXS 0506+056. Our findings seem capable of explaining the neutrino generation at the time of a low gamma-ray flux and also indicate that TXS 0506+056 might be an atypical blazar. It seems to be the first time that a potential collision of two jets on parsec scales has been reported and that the detection of a cosmic neutrino might be traced back to a cosmic jet-collision.


Introduction
It has been speculated for many years that cosmic neutrino production can be assigned to the most powerful accelerators in space, the jets of active galactic nuclei (AGN; e.g., Ginzburg & Syrovatskii 1963;Rachen & Biermann 1993). Only recently, however, could the origin of the first neutrino be linked to a BL Lac object: TXS 0506+056 (IceCube Collaboration et al. 2018a). The identification was possible thanks to the multiwavelength flaring (from radio to TeV) of TXS 0506+056, observable by many ground-and space-based telescopes (e.g., Padovani et al. 2018). The neutrino can be associated with the blazar TXS 0506+056 with chance coincidence being rejected at the ∼3σ level (e.g., Ansoldi et al. 2018). An analysis of archival IceCube data revealed evidence of an enhanced neutrino activity between September 2014 and March 2015 (IceCube Collaboration et al. 2018b).
TXS 0506+056 is a BL Lac object at a redshift of z = 0.3365 ± 0.0010 (Paiano et al. 2018). While the identification of the origin of the detected neutrino with this BL Lac object looks solid, many other BL Lac objects show properties similar to those of TXS 0506+056, but have not been identified as neutrino emitters. Multiwavelength flaring, observed at different frequencies, is a regular phenomenon of BL Lac objects or blazars. To our knowledge only one other flaring object, a flatspectrum radio quasar (FSRQ), has been discussed as a candidate for neutrino emission in combination with multiwavelength flaring (PKS B1424−418, Kadler et al. 2016). With this paper we  (Nov. 14, 2010) marked. Shown is an image with the difmap modelfit components superimposed (circles) and the components belonging to the inner jet (light blue arrow) and those features that seem to belong to the spike (red arrow). The arrows indicate the dominant direction of motion. The filled ellipse in the bottom left corner gives the beam size (point spread function). Panel b: image from a later epoch. The two potential jet scenarios are indicated (see text for details). The physical nature of the spike is unclear. The spike in the two-jet scenario might be the jet of a potential second black hole. The second core (core II) is therefore marked with a question mark. study the question of why and how the neutrinos were produced if they originate in TXS 0506+056. We later propose that the enhanced neutrino activity in 2014-15 as well as the extremely high-energy (EHE) neutrino were generated in a collision within the jet.
Throughout the paper we adopt the following parameters: a luminosity distance D L = 1.8277 Gpc at the source redshift of z = 0.3365 with cosmological parameters corresponding to a ΛCDM Universe with Ω m = 0.308, Ω λ = 0.692, and H 0 = 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 16.18c, while 1 mas = 4.961 pc.

Observations and data reduction
We remodeled and reanalyzed 16 VLBA observations (15 GHz, MOJAVE 1 ) performed between January 2009 and March 2018. The typical beam (point spread function) of these observations is 0.984 × 0.444 mas. Observations on September 6, 2015, had the smallest beam with 0.942 × 0.406 mas. The beam indicates the resolution, which is usually assumed to be one-fifth of the beam. The beam of the observations depends on the projected length and orientation of the vector between the antennas as viewed by the celestial source. This changes as the Earth rotates.
Gaussian circular components were fitted to the data to obtain the optimum set of parameters within the difmap-modelfit program Shepherd (1997). Every epoch was fitted independently from all the other epochs. The following parameters were fitted to the data: the flux-density of the component, the radial distance of the component center from the center of the map, the position angle of the center of the component (degrees north -> east) with respect to an imaginary line drawn vertically through the map center, the full width at half maximum (FWHM) major axis of the circular component.
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, without knowing what these peculiar properties would be, we remodeled all the available MOJAVE data to also allow very faint components.
With this new approach, any unusual morphology and morphological changes should appear in the model-fitting process. To allow for unusual morphologies, we also applied a new technique in the identification of the jet features. We identified patterns in the jet motion by searching for a smooth evolution in x-and y-coordinates with time for individual jet features. In particular, faint components do reappear in later epochs at similar places, which makes us confident that these jet features are in fact real. The excellent data quality of the MOJAVE program allows such a deep analysis.
In Fig. 1a we show one of the images with difmap model components superimposed. The point spread function (beam) is shown as well.
Special care was taken to correctly identify the core component in every individual data set. We used the brightest jet feature as reference point for all the jet components. Uncertainties Panel a: core distances of the jet features that could be traced through the epochs in the one-jet scenario. The proper motions were determined via linear regression (dashed lines), for feature 1h a quadratic fit is shown. For jet components 1c and 1h some components are crossed by an orange line; these lines mark those epochs where the component motion seems to be along a curve. The colored rings refer to the component's maximum flux-density (see panel d). Panel b: component paths in x and y for those features shown in panel a. Again, the orange lines mark the epochs where motion along a curve (or strongly bent structure) seems to be observed. The black arrows indicate the direction of motion in the case of the one-jet scenario. Panel c: paths of the three innermost components that could be traced in panels a and b. Panel d: flux-density evolution of the jet components 1b-1h as a function of time. The colored circles indicate the times of their maxima. The same circles and colors are used in panels a and b.
of the modelfit component parameters were determined using bootstrap (Efron 1979). For bootstrap applications to estimate uncertainties of very-long-baseline interferometry (VLBI) results and details of our implementation, see Pashchenko (2019). We bootstrapped the adjusted residuals between selfcalibrated interferometric visibilities and the best-fit difmap model. The residuals were first filtered from outliers and then centered, and the distributions of the adjusted residuals were fitted using kernel density estimates (KDEs). The kernel width was determined from the data using cross-validation (Hastie et al. 2001). This was done independently for each baseline, correlation, and frequency sub-band. We then added samples of the residuals from the fitted KDE to the model visibilities obtained from our best difmap model. The resulting bootstrapped visibility data sets were fitted in difmap using the original bestfit model as the initial guess. It was done 300 times for each epoch. Thus, we obtained a distribution for each parameter. The standard deviation was used to estimate the corresponding uncertainty.
The obtained errors are random errors due to thermal noise (Pashchenko 2019). They are larger for weak components and data with lower signal-to-noise ratio as expected (Fomalont 1999). This can be seen from Fig. 2c where the position uncertainty is lower for recent epochs with high-sensitivity data. Although our error estimates capture the dependence between difmap model parameters, they do not account for a systematic error uncertainty due to self-calibration and core shuttle effect (Lisakov et al. 2017). The former results from the true brightness distribution that is non-Gaussian. For high-sensitivity data sets the systematic errors dominate the total uncertainty. 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). The uncertainty of the component size has been rounded up. Table A.1 lists all the parameters of the model-fitting procedure and their uncertainties. The parameters are the epoch of the observation, the flux-density of the component, the position in x and y, the size of the component, and its identification. The core component (reference position) is labeled with a "c." Components with no identification could not be reliably traced between epochs. The jet features that could be traced reliably through at least three epochs are labeled 1a-1a2, 1b-1h.

Results of the VLBA data analysis
We tested several identification scenarios based on the available data and present the kinematic scenarios which seemed to be the most reliable and robust. In the following, we propose two scenarios: one strongly curved jet-structure or two jets in collision.

Core identification
In general, the absolute position information is lost in the VLBI data reduction. We usually assume that the brightest component in a VLBI map corresponds to the core position. Any motion of jet features is calculated with respect to this core position as a reference point. For TXS 0506+056 the brightest component of the jet can easily be identified and traced across the epochs. We thus use this feature as the reference position and calculate all component motions and apparent speeds with regard to this core position. The morphology of the jet of TXS 0506+056 is more complex compared to other AGN jets. As a comparison, for the quasar 1308+326 (Britzen et al. 2017) a single jet with outward directed motion is found. The kinematics is not so clear in TXS 0506+056. In the following, we discuss a scenario, where a second jet might be involved. This second jet would require a different origin (core), not coinciding with the brightest feature of this jet. The position of this second core (as indicated in Fig. 1b) cannot be determined unambiguously from the currently available data. We thus use the one brightest component as reference position for all the jet components and their apparent motions.
3.2. Two scenarios: one strongly curved jet-structure versus a binary jet At first sight the jet of TXS 0506+056 reveals a core-jet structure typical of BL Lac objects. However, a striking feature in the radio emission of the parsec-scale jet is a straight "spike" that seems to enter from the right in the xy-plane. In Fig. 1a we indicate this atypical spike structure with a red arrow. This spike appears to be disconnected from the inner jet structure visible close to the core (blue arrow in Fig. 1a). In the following we analyze whether the spike is part of one strongly curved jet (one-jet scenario) or whether the spike belongs to a second jet (two-jet scenario). We analyze the kinematics of the individual features in both parts of the jet, the inner jet and the spike.
To visualize the jet features and their positions in TXS 0506+056, we show in Fig. 1a one example of a map and the identified jet features. This is 1 of 16 epochs of data that have been reanalyzed, and only those components that could be model-fitted in this one epoch are labeled. The number of jet components that can be model-fitted varies with time and data set. None of the epochs of data shows the whole set of jet components fitted to the data over the nine years of data reanalyzed. The components labeled 1a, 1a1, 1a2 belong to the inner jet; 1b, 1d, 1f, 1g are part of the spike-structure. Components 1c and 1h might be part of a curved structure that connects the inner jet with the spike-structure. The evolution of the jet features in position as a function of time is shown in Figs. 2a-c.

First scenario: a strongly curved jet-structure
We analyzed the jet kinematics assuming that the jet is one jet only. We show the results in the plots in Fig. 2. Figure 2a displays the core distance as function of time for those jet components which could be reliably traced. Figure 2b shows the paths of those jet components in the xy-plane. In Fig. 2c we zoom into the paths of the inner three jet components (1a, 1a1, and 1a2). Their paths are clearly displaced with time (from right to left). The most recently appearing jet feature moves on a straight path in contrast to the other features. While the first two jet components (1a, 1a1) reveal strongly curved paths, the third jet feature (1a2) moves along a rather straight path. These three jet features closest to the core move at slower apparent velocities compared to the rest of the jet components that could be traced. Table 1 lists the calculated proper motions and apparent speeds.
Tracing the jet features in the spike is difficult; however, we think that we managed to find a solid identification of components 1b-1h. All of them move almost perpendicular to the direction of the inner three jet features (Fig. 2b). For at least two outer jet features (1c, 1h) we find evidence for a change in the direction of motion in a pronounced curve. These two components moving in a curve are crossed in Figs. 2a and b by orange lines. A detailed investigation of the spike reveals that some of the components seem to move on the same path (1g, 1h). If the jet of TXS 0506+056 consists of only one jet, then it seems that we are observing very strong curvature effects in the jet of TXS 0506+056. We discuss implications for the angle to the line of sight in Sect. 8.

Second scenario: two jets on a collision course
We are not aware of any similarly strongly curved jet with components moving on such curved paths. We thus consider another scenario in which the spike belongs to a second jet. The parsecscale jet of TXS 0506+056 would thus be comprised of two jets. We sketch this scenario in Fig. 1b. Jet I would then be the dominant jet with the core of the parsec-scale jet and the inner jet components. Jet II would be the spike that seems to enter from the right in the xy-plane. Jet I shows apparent speeds away from the core while jet II reveals apparent motion in the direction from A103, page 4 of 17 S. Britzen et al.: A cosmic collider: Was the IceCube neutrino generated in a precessing jet-jet interaction in TXS 0506+056? right to the left in the image plane. We may speculate that this second jet could originate from a second black hole (likely direction to second black hole indicated in Fig. 1b with a question mark). The paths of the two jets seem to intersect.
Jet II seems to be closer to our line of sight and more strongly Doppler beamed. The apparent speeds of the jet features are higher (Table 1). The neutrino(s) might have been generated in an interaction of two jets in a possible active supermassive binary black hole system. We discuss this in more detail in Sect. 8.
Tracing the flux-density of each individual component (1b-1h) reveals that each of these jet components goes through a maximum in the flux-density evolution (see Fig. 2d, marked by differently colored rings). We checked where in the xy-plane this maximum of flux-density occurs (Fig. 2b), and for all the components 1b to 1h this flux peak is at the projected crossing site of the inner jet features (1a-1a2) with the outer jet. This brightening in flux-density is of small amplitude since all these components appear to be fainter than the inner jet components. This brightening and passing of the collision site occurs for 1h roughly at the time of the neutrino event (Fig. 2d). Since all the jet components "lighten up" in the same region, this is further support for the collision scenario we proposed earlier.
We can estimate the interaction angle following the typical interrelation between the viewing angle of the jet components and the apparent superluminal speed. We sketch the model geometry of the proposed jet interaction in the two-jet scenario in Fig. 3. For jet II with apparent velocities of about β II 10, a viewing angle between α II = 1.5 • and α II = 3 • is needed in combination with a proper Lorentz factor γ II = 10 − 16. In contrast, for jet I with an apparent velocity of about β I 2, a viewing angle between α I = 5 • and α I = 20 • is likely in combination with a proper Lorentz factor γ I = 2 − 4. We can therefore estimate the collision angle between Θ = α I + α II and Θ = α I − α II , which also depends on the 3D alignment of the jet sources. In numbers, this gives a Θ 15 • and a large difference in jet speed, or a Θ 2 • and a smaller difference in jet speed. In any case the difference in the Lorentz factor is ∆γ 10, which may indeed lead to a powerful collision between the two jets.

Apparent speeds of jet features
The kinematic parameters derived for the components that could be traced are listed in Table 1. As mentioned before, we used the core feature as reference point for all the components, which could be reliably identified. We refer to Fig. 1a for the labeling of the components. The estimated time of ejection can only be listed for those features that are part of the inner jet (1a, 1a1, 1a2) as we assume their origin is the core. The position angles for only these components are listed in Col. 6. For the remaining components it is not clear when and from which core they were ejected.
The proper motions listed have all been determined by applying linear regression fits to the data (see Fig. 2a). If the one-jet scenario is the proper model describing our results, all components (1a-1h) will belong to jet I. If the two-jet scenario is the proper description, components 1a, 1a1, 1a2 will belong to jet I, and components 1b-1h to jet II. Components 1c and 1h move along a curved path in the xy-plane. In the case that the twojet scenario is correct, the first data point for 1c might not be correctly identified and for 1h the first three data points might be misidentified. Only the component identifications for 1c and 1h need to be modified for switching between the one-jet and two-jet scenario. The proper motions of components 1c and 1h in Table 1 have been calculated for the two-jet scenario, without those components that move on the curved path.

Modeling the precessing inner jet
As listed in Table 1 and indicated in Fig. 2c, the jet components of jet I, 1a, 1a1, and 1a2, start under different position angles with time. The flux-density and apparent velocity of 1a2 is higher compared to those of 1a and 1a1. Interpreted in geometrical terms, the jet in more recent times is pointing closer to the line of sight (i.e., moving toward us). This might indicate that this jet is precessing.
We tested this hypothesis by constructing a simple precession model, which was previously used to explain periodic radio variability and component positional and apparent velocity changes in blazars (e.g., for OJ287 by Abraham 2000;Britzen et al. 2018) and a radio galaxy (for 3C120 by Caproni 2004). Within the precession model the jet is envisaged to undergo a bulk precession motion, which leads to the periodic changes of the Doppler boosting factor due to the changes of the viewing angle of jet components. As a consequence, the underlying emission of the jet is affected with respect to the observer since the flux density in the observer's frame depends on the Doppler factor as where S jet is the underlying non-thermal synchrotron jet emission in the source frame, S jet ∝ ν −α , with α the spectral index, and δ(γ, φ) is the Doppler-beaming factor, which depends on the Lorentz factor γ, the jet component velocity (1), the power-law index ξ is equal to the sum of the geometrical factor p and the spectral index α, ξ = p + α, where p = 2 for a continuous cylindrical jet and p = 3 for the discrete jet composed of spherical components. In the following, we focus on the components of jet I (1a, 1a1, 1a2) and their variable apparent velocities β app , position angles η, and the variable flux density of the core I. The model is based on the precession equations that are discussed in detail in Abraham (2000), Caproni (2004), Britzen et al. (2018), among others, which we introduce here for clarity. The apparent velocity β app of jet components is given by where the component velocity β is related to the Lorentz factor γ L by γ L = 1/ 1 − β 2 . The viewing angle at any time Φ(t) is the angle between the direction of the motion of the component and the line of sight, while the position angle η(t) is the projection of the component position on the sky plane. The temporal evolution of the two angles can be expressed simply as where the Cartesian coordinates in the observer frame are given by the relations where η 0 is the position angle of the precession cone on the sky plane. The time-dependent coefficients A and B can be expressed as functions of the half-opening angle Ω of the precession cone, the viewing angle of the precession-cone axis φ 0 , and of the constant angular velocity of the bulk precession of the jet ω prec = 2π/P obs prec , where P obs prec is the precession period in the observer frame. The temporal evolution of A and B coefficients with respect to the reference epoch t 0 is then given as follows: In addition, this set of equations is complemented by the observed flux density S obs modulated by the Doppler-boosting factor δ (see Eq. (1)). The illustration of the precession model with the geometry description is in Fig. 4. Using least-squares fitting, we fit the temporal evolution of flux density, apparent velocity, and position angle simultaneously to obtain the eight best-fit parameters that are listed in Table 2. The reduced χ 2 ν of the global least-squares fitting is χ 2 ν = χ 2 /ν = 17.78/(21 − 8) = 1.37. The corresponding plots for the flux density, the apparent velocity, and the position angle are in Fig. 5 (top row), which also includes fit residuals in the bottom part of each plot. All fitted parameters have large errors comparable to the best-fit values. This is mainly due to the small number of measured points (21) with respect to the number of parameters (8). In the next step, we fixed the value of the Dopplerboosting exponent ξ. We selected two extreme cases with the realistic range: -We fixed ξ to 1.5, which corresponds to the spectral-index value of α = −0.5 (cylindrical geometry) and α = −1.5 (spherical discrete geometry). Both cases represent optically thick synchrotron emission, S ∝ ν 0.5−1.5 , typical of radio cores; -We increased the value of ξ to 3.0, which corresponds to α = 1 for the cylindrical geometry (optically thin synchtrotron emission, S ∝ ν −1 , typical of older radio lobes) or α = 0 for discrete emission components (flat spectrum). With these two values of the Doppler boosting exponent, we essentially cover the whole range of the radio-spectral indices that were found within optical diagnostic diagrams of galaxies (see Zajaček et al. 2019, for details). The blazars that are usually associated with LINERs having lower line-ratio log ([O iii]/Hβ), typically have flat to inverted radio spectra with an offset in the radio-spectral index distribution with respect to Seyfert galaxies (Zajaček et al. 2019).
Both fits with the fixed Doppler-boosting exponent are shown in the middle and bottom rows of Fig. 5. The remaining seven parameters are generally better constrained (see Table 2) being consistent within uncertainties with the general fit. The precession period is ∼10 years and the Lorentz factor is ∼2 for both fits with the fixed Doppler-boosting exponent. The precession angles remain largely unconstrained, mainly because we fit the apparent velocities and the position angles for only three components of jet I. However, the half-opening angle of the precession cone and its viewing angle converge toward smaller values (Ω ∼ 20 • , φ 0 ∼ 10 • ), which is consistent with the blazar jet being viewed closer to the line of sight.
In the case of TXS 0506+056, the jet precession may be driven by the torques of the secondary massive black hole, which could be associated with core II, but is not necessarily. In the following analysis we consider the estimate of the black hole mass based on the bulge mass of the elliptical host (Padovani et al. 2019), M • ∼ 3 × 10 8 M . It is difficult to infer the real supermassive black hole binary (SMBHB) distance from the VLBI image at 15 GHz, as the core distance is comparable to the beam size of ∼1 mas = 4.961 pc. If the projected angular distance corresponded to the real black hole separation, then the precession period would have to be generally longer than which is simply estimated from the third Kepler law approximately applicable to the bound binary system of two massive black holes. This value is almost four orders of magnitude larger than that inferred from the precession model fit. However, the angular distance between core I and core II does not simply   6. Panel a: observed period of the SMBHB system in years as a function of the primary mass fraction and of the decadic logarithm of the orbital separation of the black holes in milliparsecs. The considered total black hole mass is M tot = 3 × 10 8 M according to Padovani et al. (2019). The plot represents the complete parameter space in terms of (x p , r ps , P obs ps ) that would yield the observed disk-jet precession period of 10 yr. Panel b: color-coded merger timescale in years as a function of the primary black hole mass ratio x p and of the decadic logarithm of the mutual component distance r ps in milliparsecs. The whole parameter space (x p , r ps ) is shown that yields merger timescales longer than the source-frame SMBHB orbital period, τ merge > P ps .
correspond to the physical distance between two putative black holes. Core II could be an optically thick component of jet II, not really associated with the secondary black hole.
To estimate the actual separation between two components of the SMBHB, we consider the model of a rigid disk precession around the primary black hole, which can be driven by the torques from the secondary black hole. The rigid precession works under the assumption that the secondary black hole is located outside the outer disk radius, r out d < r ps . This scenario is expected in a merging system, where a misalignment between the accretion disk around a primary black hole and the orbital plane of a secondary component can be natural. The precession driven by a SMBHB was already applied to explain optical flares of OJ287 with the period of ∼12 years (Katz 1997). The brightening occurs when the precessing jet approaches the light of sight. The precessing disk-jet in the SMBHB systems has a typical period of several years (Caproni 2004).
We consider the component masses M 1 and M 2 with the total mass of M tot = M 1 + M 2 and the orbital period of P ps in the source frame. The separation between the primary (core I) and the secondary (core II) component is then given by the third Kepler law, where the period in the source frame is related to the observed period simply as P ps = P obs ps 1 + z · Precession will proceed when the disk plane around the primary black hole (core I) and the orbital plane of a secondary black hole (associated with core II) are misaligned by an angle θ, which is equal to the half-opening angle of the precession cone. In our case, we consider θ ≈ Ω ∼ 20 • as the rounded value obtained from the fits with the fixed Doppler-boosting exponent (see Cols. 4 and 5 of Table 2). The equality Ω ∼ θ applies under the assumption of an efficient disk-jet coupling. For the observed precession period, we take the averaged best-fit value of P obs prec = 10.0 yr, which is supposed to be equal to the sourceframe disk and jet precession, P disk ≈ P jet = 7.48 yr. Then the outer radius of the precessing disk around the primary can be expressed as a function of the primary mass ratio, x p = M p /M tot (Papaloizou & Terquem 1995), where n is the polytropic index of the gas in the accretion disk, with n = 3/2 and n = 3 for non-relativistic and relativistic gas, respectively, which we set equal to 3/2 in our calculations, considering the non-relativistic case. In Fig. 6, we plot the colorcoded observed orbital period of the binary as a function of the primary mass ratio (x p = 0.5 − 1) and of the logarithm of the binary separation r ps between the two black holes in milliparsecs. For the rigid precession of the primary disk to take place, the minimal orbital separation of the two black holes is r pc = 0.01 mpc (the mass ratio in the full range x p = 0.5−1.0), which corresponds to the observed orbital period of P obs ps = 2.3×10 −4 yr. The maximum separation is r pc = 8.5 mpc (with the mass ratio x p ∼ 0.5), which corresponds to the observed orbital period of P obs ps = 5.7 yr. The separation between the two black holes can be further constrained by the merger timescale τ merge , which is the timescale on which the two black holes at the initial distance of r ps will merge. The merger timescale can be expressed as (see Shapiro & Teukolsky 1983, for a general expression) where x p and x s are the primary and the secondary black hole mass ratios, respectively. Under the assumption that the binary will merge on a timescale longer than the SMBHB source-frame orbital period t merge > P ps , it limits the binary distance to the milliparsec scale, specifically to r ps = 0.0832 − 8.511 mpc (see Fig. 6, right panel). r out d = 100 r ms s=0 s=-0.5 s=-1 s=-2 Fig. 7. Lense-Thirring precession period compared to the potential co-moving precession period of the inner jet of TXS 0506+056 (green solid line; P LT = P disk = 7.48 yr). The bottom axis represents the black hole spin value from the maximally retrograde (a = −1) to maximally prograde (a = 1). Different lines represent the Lense-Thirring precession period for different slopes of the surface density profile of the inner accretion disk, using the notation Σ ∝ r s . Left panel: outer radius is r out d = 10 r ms , where r ms is the marginally stable radius, which depends on the spin value. Middle panel: r out d = 50 r ms . Right panel: r out d = 100 r ms .
Another potential source of disk-jet precession is the Lense-Thirring precession due to the misalignment of the accretion disk angular momentum and the spin of the corresponding primary black hole, presumably within core I. For the TXS 0506+056 source this situation may occur for the case when the two black holes are be located at a relatively large projected distance of 1.5 mas = 7.5 pc. The secondary black hole (within core II)surrounded by an accretion disk as well, and possibly by a bound stellar system -can be the source of outflows, both a disk wind or stellar winds. The primary black hole (core I) may then capture the outflow material associated with core II with an initial angular momentum given by the orbital speed of core II. The captured material has an angular momentum that is in general misaligned with respect to the spin of the primary black hole, which therefore implies the non-zero Lense-Thirring torque on concentric rings within the accretion disk. The misalignment may of course occur independently of core II. However, the two-core scenario provides a possibility of the disk precession around the first black hole triggered by the second black hole system, which is distant and thus cannot induce direct torques on the primary disk on timescales of 10 years.
We note that the scenario discussed above has similarities to X-ray binaries hosting a component launching an outflow (star) and a compact accreting object. For some sources, a jet precession has been reported, namely for SS 433 (Monceau- Baroux et al. 2015) or Cyg V404 (Miller-Jones et al. 2019). The jet precession is therefore not a rare phenomenon, but is generally expected to occur during the lifetime of supermassive black holes and X-ray binaries as they evolve along the hardnessluminosity diagram (Fender et al. 2004).
As was demonstrated for the blazar OJ287 with the observed precession period of the same order of magnitude (20 years, Britzen et al. 2017), the precessing inner part of an accretion disk must be spatially compact to yield a small enough precession period, with the outer disk radius of the order of 100 r ms , where r ms is the radius of marginally stable orbit. Blazars exhibit very low quiescent optical emission; therefore, the inner parts of their accretion flows are expected to be hot and optically thin, geometrically thick, advection-dominated accretion flows (ADAFs, or hot flows; Yuan & Narayan 2014), with the surface density slope of s ≈ −0.5, where the surface density of the accretion flow is Σ ∝ r s .
In general, the amplitude of the precession frequency ω LT decreases with the cube of the distance from the rotating black hole. Considering the precession of an accretion disk that behaves like a rigid body, the period is determined by the black hole angular momentum J • = GM 2 • |a|/c 2 , where a is the dimensionless spin parameter 2 . Another crucial parameter is the radial extent of the precessing part of the accretion flow, with the inner radius R in = ξ in R g and the outer radius R out = ξ out R g , where R g = GM • /c 2 ∼ 4.4 × 10 13 (M • /3 × 10 8 M ) cm is the gravitational radius. The inner radius is set equal to the marginally stable orbit, which depends on the spin of the Kerr black hole, • +A 2 1 ) 1/2 . For the maximum prograde spin, we obtain ξ ms = 1, and for the maximum retrograde spin we have ξ ms = 9. As already mentioned, the precession period also depends on the surface density of the accretion flow, Σ = Σ 0 ξ s . The Lense-Thirring precession period in the source frame may then be expressed as ) where Υ(ξ) = ξ 3/2 + a • and Ψ(ξ) = 1 − (1 − 4a • ξ −3/2 + 3a 2 • ξ −2 ) 1/2 . The outer radius of the precessing disk ξ out is set to 10, 50, 100 times ξ ms , for which the precession period partially overlaps the source-frame disk-jet precession of 7.48 years. Figure 7 shows the results of our calculations for TXS 0506+056, formally for a range of surface density slopes, s = {0, −0.5, −1, −2}, but focusing mainly on the ADAF slope of s = −0.5, for which the case of a solid-body precession applies. This stems from the consideration that in the large scale-height accretion flows, warps induced by a misaligned black hole propagate as bending waves that approximately spread on a soundcrossing timescale, which is shorter than the precession period. This is a fundamental difference with respect to the classical thin disk where warps propagate on a viscous timescale and the stable Bardeen-Petterson configuration can develop (with the spinalignment up to the Bardeen-Petterson radius, see Ingram & Done 2012, and references therein). Therefore, thick hot flows essentially behave as solid bodies, which was also confirmed in global numerical simulations (Fragile et al. 2007). We conclude that the Lense-Thirring precession with a period in the co-moving frame of P LT = P disk = 7.48 yr is possible for an outer (precessing) disk radius in the range of r out d = 10−100 r ms , with a spin ranging from small or intermediate values (both prograde and retrograde) to large values (only prograde) of the primary black hole of M • = 3 × 10 8 M , for an increasing radial extent of the precessing part of the accretion flow (see Table 3 for an overview). Figure 8a shows the single-dish flux-density obtained in 15 GHz monitoring observations with the Owens Valley Radio Observatory (OVRO). In addition, we show the VLBI core flux density obtained from the VLBA data analysis. TXS 0506+056 has been observed by OVRO more frequently than via VLBA observations. Thus, the single-dish light curve is significantly better sampled compared (Fig. 8a) to the VLBI flux-density curves (see Fig. 8b), and the VLBA observations do not cover all the peaks in the single-dish light curve.

Correlation between the radio and the γ-ray light curve
In general, the radio light curve obtained by the OVRO program is dominated by the flux of the VLBI core. The flux-density of the core is the main contributor to the overall single-dish flux. The major flare in flux-density between 2008 and 2010, and the steep flare starting in 2016 all seem to come from the core region. However, between mid-2010 and 2016, the core flux is relatively stable, but the single-dish flux reveals some flaring. The smaller flares around 2014 and 2015.5 seem to originate in the jet features. In particular, the source reveals a rising radio flux at 15 GHz at the time of the enhanced neutrino activity. The major contributor to the flaring in addition to the core flux at the time of the enhanced neutrino activity is jet feature 1a1, which is brighter than the other jet components at that time. In the following, we study a possible correlation between the radio and the γ-ray light curve.
For the GeV regime we used Pass 8 photon data from the Fermi-LAT (Atwood et al. 2009) of the time span August 4, 2008 until August 3, 2018 (MJD 54682−58333), downloaded from the Fermi-LAT data server 3 . We used version v10r0p5 of the Fermi ScienceTools 4 . We computed a light curve from the photon data by fitting a power law of the form with a prefactor N 0 and a photon index γ to be determined by fitting the data and the scale E 0 fixed to its catalog value. Including data within 15 • around the position of TXS 0506+056 and the energy range E = 0.1−300 GeV, this fit is performed for 3 https://fermi.gsfc.nasa.gov/cgi-bin/ssc/LAT/ LATDataQuery.cgi 4 https://fermi.gsfc.nasa.gov/ssc/data/analysis/ software/ OVRO flux (15 GHz) VLBI core flux (15 GHz) jet feature 1a1 core + inner jet features (1a, 1a1, 1a2) core + VLBI jet features flux (15 GHz) Fig. 8. Flux-density of jet feature 1a1, VLBI core, core plus inner jet features, core plus all jet features, and single-dish flux-density observed with the OVRO telescope. All measurements were obtained at the same frequency of 15 GHz. The major increase in the radio flux-density at 15 GHz from 2016 onward seems to be produced by the core region.
every time bin of width 45 days. All parameters of all sources within 5 • around TXS 0506+056 were left free for the fit as well. Sources between 5 • and 25 • were included in the likelihood analysis with all parameters fixed to their catalog values. The Galactic model file gll_iem_v06.fits is used to model the diffuse emission, and we use the isotropic spectral template iso_P8R2_SOURCE_V6_v06.txt.
In the upper panel of Fig. 9 the GeV light curve (blue) is plotted together with the radio light curve (red) resulting from longterm blazar monitoring by OVRO (Richards et al. 2011 The single neutrino IceCube-170922A is clearly coincident with a GeV flare visible in our light curve and in agreement with the results by Padovani et al. (2018), among others. This single flare is superimposed onto a general increasing trend of the GeV flux starting at 2015 and declining since 2018, soon after the neutrino event. The radio light curve presents a similar increasing trend, but lagging behind the GeV fluxes by about one year. In addition, the radio light curve is superimposed by three subflares, the second of which coincides with the GeV flare and the neutrino event. The radio observations stop here with a flaring in mid-2018 similar to the GeV flare in 2017.
The neutrino detection prior to the single event IceCube-170922A, indicated by the gray-shaded area in Fig. 9, does not correspond to a distinct GeV flare. However, it is worth noting that these neutrinos were detected during the rising part of a radio flare, similar to the single neutrino event. A striking feature of the radio-GeV relation during these neutrino detections is that they are anti-correlated throughout the gray-shaded area, whereas in the other parts of the light curves they seem to be more correlated. In the following we attempt to quantify this visual impression.
In order to investigate the local correlation between the radio and the GeV light curves, we compute the Pearson correlation coefficient for adjacent parts of the data (i.e., for time windows of 100 days). Since the radio and the GeV data are not simultaneous we linearly interpolate between data points with A103, page 10 of 17 S. Britzen et al.: A cosmic collider: Was the IceCube neutrino generated in a precessing jet-jet interaction in TXS 0506+056? a time resolution of 10 days. The result of this procedure is shown in the middle panel of Fig. 9. The visual impression of anti-correlation during the gray-shaded part is indeed confirmed quantitatively by negative values of the correlation coefficient during that time interval. In addition, the single neutrino event is accompanied by a negative correlation between the radio and the GeV data. However, all the data are characterized by an ongoing change between positive and negative correlation coefficients, and is not restricted to the epoch of neutrino detection alone.

High-energy emission from blazars
TXS 0506+056 belongs to the class of blazars comprising BL Lac objects and γ-ray loud FSRQs. These blazars are the most extreme active galaxies known. They emit over a broad range of frequencies, from the radio to very high-energy (VHE) γ rays. Variability is detected at all wavelengths on different timescales. The generally accepted model is that the low-frequency component of the blazar spectral energy distribution (SED) might be synchrotron radiation from non-thermal, ultrarelativistic electrons. Two different approaches are being studied with regard to the high-energy emission.
The leptonic origin of gamma rays can be due to inverse Compton interaction (relativistic electrons locally accelerated collide with low-energy photons). The photons can be produced within the jet (SSC process: e.g., Marscher & Gear 1985) or external to the jet. External seed photons can be produced in the accretion disk (Dermer & Schlickeiser 1993) or can be reprocessed in the circumnuclear material, for example in the broadline region of quasars (e.g., Sikora et al. 1994;Dermer et al. 1997).
A hadronic origin can be due to p-γ interactions (relativistic protons collide with low-energy photons) or pp interactions (relativistic protons collide with low-energy protons). If a significant fraction of the kinetic power in the jet is converted into proton acceleration, these protons can reach the threshold for pγ pion production. Then synchrotron-supported pair cascades will develop (e.g., Mannheim & Biermann 1992). In order to reach the interaction threshold for pγ interactions on soft target photons and/or to contribute significantly to the radiative output via proton-synchrotron radiation, the protons need to be accelerated to ultra-relativistic energies (E p 10 17 eV). To accelerate PeV protons, large turbulent magnetic fields are required. During pion-and muon-decay, neutrinos can be generated that cannot be produced within purely leptonic models. For a review on the different blazar models, see Böttcher (2007).

Photo-hadronic neutrino production in interacting jets
In this section we consider the feasibility of the production of 100 TeV neutrinos detected by IceCube from the direction of TXS 0506+056 by photo-hadronic interactions in the jets of this blazar. The underlying assumptions are that relativistic protons in the ∼100 TeV-10 PeV range are accelerated in jet II and that the synchrotron photon field of jet I acts as the target photon field for pγ pion production. As elaborated in Reimer et al. (2019), the production of ∼ 100 TeV neutrinos requires protons of energy E p 200E 14 /δ II,1 TeV (i.e., γ p 2 × 10 5 E 14 /δ II,1 ) in the co-moving frame of the electromagnetic and neutrino emission region, where E ν ≡ 100E 14 TeV is the neutrino energy and δ II ≡ 10δ II,1 is the Doppler factor with which electromagnetic and neutrino emission from the jet source is boosted into the observer's frame. IceCube neutrino production will be most efficient with target photons of co-moving energy E t ≥ 1.4δ II,1 /E 14 keV, in the co-moving frame, corresponding to p-γ interactions at the ∆ + resonance energy. As argued in Reimer et al. (2019), protons of the required energies are easily confined, and thus plausibly accelerated to ∼PeV energies in the typically assumed ∼10-100 G magnetic fields typical of (lepto-)hadronic blazar models.
The superluminal motions measured from the two interacting jets indicate characteristic Lorentz factors of their bulk motion of Γ II ∼ Γ I ∼ 10. The relative motion of two jet components with respect to each other is then characterized by the Lorentz factor Γ rel = Γ II Γ I (1 − β II · β I ). Due to projection effects, the actual angle between the velocities of the two jets is difficult to constrain, but we can conservatively estimate the relative Doppler factor as Γ rel = 10 Γ rel,1 .
We further assume that the distance between the neutrinoproducing component of the inner jet and the synchrotron X-ray emitting component of the spike are separated by an angular distance of d ab,ang = 1 d mas mas. The angular diameter distance to TXS 0506+056 is d A = d L /(1 + z) 2 ∼ 1.0 Gpc ∼ 3.0 × 10 27 cm, so that the physical (projected) distance between the two components is d ab = 4.8 d mas pc. Thus, the neutrino-producing component of jet II will receive a synchrotron flux from jet I, which corresponds to a co-moving energy density of u t,a ∼ 16 Γ 4 rel L X,I 4 π d 2 ab c ∼ 2.2 × 10 6 Γ 4 rel,1 d −2 mas erg cm −3 .
As detailed in Reimer et al. (2019), such a dense target photon field will ensure that photo-pion induced pair cascades will be strongly in the Compton-supported regime, in which case the intrinsic γγ opacity is so high that any emerging γ-ray flux is far below the observed Fermi-LAT flux from TXS 0506+056, even in quiescence (as during the 2014-2015 neutrino flare). Thus, the production of neutrinos via this mechanism is in agreement with the absence of γ-ray activity during the 2014 -2015 neutrino flare.
Assuming that the proton spectrum continues as an unbroken power law with index α p = 2.0 from γ p,min = 1 to γ p,max > 10 6 and the proton population is constrained in an emission region of size R b = 10 16 R 16 cm, this proton spectrum normalization corresponds to a total kinetic power in relativistic protons of which is a moderate energy requirement compared to the Eddington luminosity, L Edd ∼ 1.3 × 10 46 M 8 erg s −1 of a M BH = 10 8 M 8 M black hole.
Since the neutrino event in 2017 was only a single neutrino event, we can only calculate an upper limit on the neutrino flux. This leaves too much freedom for a meaningful model interpretation. However, a scenario with lower γγ-opacity (and hence, lower photo-pion efficiency) would be in agreement with an enhanced gamma-ray flux, along with the significantly lower neutrino flux compared to the 2014-2015 flare.
We thus conclude that the interaction of the two jets in TXS 0506+056 can provide a suitable target photon field for photo-hadronic production of IceCube neutrinos as observed from this source.

Discussion and conclusions
In this paper we presented a detailed analysis of the kinematics of the parsec-scale jet of TXS 0506+056 based on archival VLBA data. In the following we discuss and summarize the implications of our findings.

Atypical jet kinematics in TXS 0506+056
Several groups have analyzed the kinematics of the parsec-scale jet of TXS 0506+056 based on MOJAVE observations before. Kun et al. (2019) find that the radio jet-components maintain peculiar quasi-stationary core distances. They find support for a helical jet structure and a small inclination angle, less than 8.2 • derived from the average brightness temperature of the core. Lister et al. (2013) report proper motions for four moving features in the parsec-scale jet of TXS 0506+056 based on the MOJAVE data: 332 ± 82, 49 ± 45, 100 ± 14, and 52 ± 16 µas yr −1 . They mention that one weak component at the western edge of the jet reveals statistically significant inward motion of 100 ± 14 µas yr −1 .
The comparison with the results of Kun et al. (2019) is complicated due to the different numbers of jet components that have been modeled to the data (four for Kun et al. 2019 and ten in our approach). Kun et al. (2019) do not detect any significant motion due to a large spread in the core separations (larger than 1.2 mas for their C4) of the individual components in their Fig. 2. In contrast, we find systematic trends in the apparent speeds of the components: an increase in apparent velocities for the inner jet components (1a-1a2) and remarkable stability in the apparent jet speeds for five of the outer jet components (1d-1h). This consistency in jet speeds supports our identification schemes. In comparison to Kun et al. (2019) we find that the curvature of the jet in the one-jet scenario is atypically large for blazars.
The jet kinematics of TXS 0506+056 is complex and its analysis complicated in particular due to specific and a priori unknown source properties that may be related to the neutrino production. Overall we find evidence for strong curvature most likely due to an extreme viewing angle which might be close to zero. As a consequence, the jet is likely to be subject to significant beaming. Jetted material within the source obviously moves in different directions so that interactions between the jet material seems unavoidable.
We note that AGN with extreme curvature in the jet kinematics have been observed before. For example, Savolainen et al. (2006) investigated PKS 2136+141 and characterized its jet as extremely curved, bending by a significant angle (210 • ) in the sky plane. Our interpretation of TXS 0506+056 is somewhat different as we think that we are looking directly into the jet cone.
Having presented evidence for a strongly curved jet in TXS 0506+056, we discussed two model scenarios, both of which might explain the observed jet kinematics. The first is the scenario of one strongly curved jet, the other the scenario of a structure made up of two jets. We find evidence for different apparent speeds of the inner part of the jet (inner jet) and the speeds of those components that either belong to the outer and older part of the same jet or to a second jet.

Precessing inner jet
We detected changes in the paths in the inner jet, which we suggest are related to the precession of the jet source. In Sect. 4 we model the precession of the inner jet components (1a-1a2). The precession period is on the order of ∼10 yr. Precession is to be expected in a binary system of supermassive black holes, and A103, page 12 of 17 has been observed in several AGN jets (e.g., Britzen et al. 2018 for OJ 287). OJ 287 is certainly the best AGN for precession studies as this source has been observed in the radio regime for several decades (for details, see Britzen et al. 2018). Precession naturally causes a change in the jet direction and speed (see, e.g., Jorstad et al. 2004 for 3C 279). For 3C 279, a precession model has been applied to explain the observed changes with regard to jet direction and component speeds (Abraham & Carrara 1998).
As our precession analysis for TXS 0506+056 shows, the precession observed for the components of the inner jet in TXS 0506+056 cannot be explained by the second core (of the two-jet scenario). The two cores would be at a projected separation that is too large to explain the observed precession period. The Lense-Thirring effect, however, might explain the observed precession. Independent support for the precession hypothesis might come from EGRET observations. TXS 0506+056 was identified by EGRET as the highest energy emitting blazar in observations in 1994 (Dingus & Bertsch 2002). We find that the increase in the radio flux-density around 2017 coincides with a major flare in the GeV data. We can only hypothesize that the observed 1994 gamma-ray brightness may have been related to another bright precession phase in TXS 0506+056.

Collision of jetted material
Tracing the individual components in this source is complicated due to an atypical source morphology and source evolution. In Table A.1 we labeled the components that could be reliably traced. As can be seen in Table A.1, jet components were modeled that could not be reliably traced in several epochs, and thus are not labeled. We find some evidence for unlabeled jet emission around the labeled jet components. Compared to the analysis of the kinematics of other AGN (e.g., OJ287 by Britzen et al. 2018), the situation for TXS 0506+056 is much more difficult. On the other hand, the additional jet emission around the identified jet features and atypical "messy" morphology supports our collision or interaction model. We do not expect to trace the moment of collision between the individual jet features and it is not clear to us whether this would be detectable at all. However, considering all the information we have (projected direction of motion of individual jet features from the inner jet and the spike, precessing motion, temporal evolution, and flux-density flaring of subsequent individual jet features while crossing the same projected position) we can fairly conclude that these features are all consistent with a jet collision scenario. To our knowledge, this is the first time that such a scenario in an AGN jet can be suggested (and therefore no comparison with other collision scenarios is possible).
We note that solid evidence for partly colliding jets has been recently claimed for a protostellar binary outflow system based on Atacama Large Millimeter/submillimeter Array (ALMA) observations (Zapata et al. 2018). This finding, derived for a non-relativistic stellar binary system, shows remarkable similarities to the scenario that we are proposing in this paper.
Our considerations are also sustained by the reasoning presented in the next subsection: the interaction of jetted material is in principle capable of generating the neutrinos and explaining the low gamma-ray flux at the time of the enhanced neutrino emission.

Neutrino production by photo-hadronic interaction
Essentially, for both of the scenarios (one jet or two jets) that we consider as possible explanations for the kinematics in the substructure in TXS 0506+056, the paths of the jet components seem to collide. These collisions may then lead to the generation of high-energy emission and in particular to the generation of a high-energy neutrino.
Standing shocks have been discussed as a consequence of collisions (see, e.g., Rani et al. 2015). However, for the case of TXS 0506+056, we do not think that a standing shock provides the proper physical scenario, since all the components we trace appear to be moving. Overall, we do not find any indication of a standing shock.
Hydrodynamical simulations of colliding jets may actually support a collision scenario. Molnar et al. (2017), modelling the jet(s) of 3C 75, find that colliding jets with a rather large velocity difference actually stick together after the collision. Jets with similar velocity bounce after the collision. The velocities we derived and that may be assigned to the two jets, γ II = 10 − 16 and γ I = 2 − 4 indeed support a scenario of a collision of two jets that after colliding and shocking stick together.
Based on the presented kinematic analysis and discussion we conclude that an interaction of jetted material, which to our knowledge has not been observed in this form before, provides a proper scenario to explain the observations and neutrino emission.
Our analysis revealed a likely site of the interaction with the lightning up of six jet components that pass this interaction site. At the time of the EHE neutrino, the inner jet part seems to be directed more toward us: the apparent speed is higher, and the jet path is straight. Around the time of the neutrino generation, jet II seems to be closer to jet I. It seems possible that the spike and in particular jet components 1g and 1h interact with the components from the inner jet at the time of the detected enhanced neutrino activity and single neutrino detection.
We describe in Sect. 7 how an interaction of two jets in TXS 0506+056 or in a strongly curved one-jet scenario can provide a target photon field for photo-hadronic production of Ice-Cube neutrinos.
8.5. Implications of the discovery of a binary AGN-jet on parsec scales We speculate above that our observations can be interpreted as a two-jet scenario. If the proposed two-jet scenario indeed provides the proper physical explanation of our observational results, this may imply important consequences. To date, AGN with jets from two active cores have been discovered with larger (kiloparsec) spatial separations (e.g., NGC 6240 based on Chandra observations by Komossa et al. 2002). AGN with double cores at small separation (parsec scales) have been pursued for many years. However, they seem to be rare and difficult to identify. A candidate for a sub-pc binary black hole has recently been discovered: the Seyfert galaxy NGC 7674 with two VLBI cores separated by projected 0.35 pc and also showing a 0.7 kpc radio jet (Kharb et al. 2017). It would be a breakthrough if our analysis presented in this manuscript had provided another candidate source for a binary black hole jet source. We conclude that TXS 0505+056 is an atypical AGN and might not be representative of the blazar class of AGN. However, this source might provide the proper setup for an interaction of jetted material and the generation of neutrinos under a special viewing angle. This could provide further high-energy events in the future.