A self-lensing supermassive binary black hole at radio frequencies: the story of Spikey continues

The quasar J1918+4937 was recently suggested to harbour a milliparsec-separation binary supermassive black hole (SMBH), based upon modeling the narrow spike in its high-cadence Kepler optical light curve. Known binary SMBHs are extremely rare, and the tight constraints on the physical and geometric parameters of this object are unique. The high-resolution radio images of J1918+4937 obtained with very long baseline interferometry (VLBI) indicate a rich one-sided jet structure extending to 80 milliarcseconds. Here we analyse simultaneously-made sensitive 1.7- and 5-GHz archive VLBI images as well as snapshot 8.4/8.7-GHz VLBI images of J1918+4937, and show that the appearance of the wiggled jet is consistent with the binary scenario. We develop a jet structural model that handles eccentric orbits. By applying this model to the measured VLBI component positions, we constrain the inclination of the radio jet, as well as the spin angle of the jet emitter SMBH. We find the jet morphological model is consistent with the optical and radio data, and that the secondary SMBH is most likely the jetted one in the system. Furthermore, the decade-long 15-GHz radio flux density monitoring data available for J1918+4937 are compatible with a gradual overall decrease in the the total flux density caused by a slow secular change of the jet inclination due to the spin-orbit precession. J1918+4937 could be an efficient high-energy neutrino source if the horizon of the secondary SMBH is rapidly rotating.


INTRODUCTION
Recently Hu et al. (2020) interpreted a narrow spike in the densely-sampled Kepler optical light curve of the quasar J1918+4937 (also known as KIC 11606854, dubbed as Spikey by Hu et al. 2020) as a result of gravitational selflensing in a supermassive black hole binary (SMBHB) system. The quasar has a spectroscopic redshift of z sp = 0.926 (Healey et al. 2008). In this scenario, the orbital plane of the binary lies sufficiently close to the line of sight so that when one of the companions -the black hole with the larger mass -passes in front of the other, the optical emission of the latter active galactic nucleus (AGN) is significantly enhanced. Taking two relativistic effects, the binary self-lensing and the orbital Doppler boosting into account, Hu et al. (2020) modeled the Kepler light curve containing the spike. They found that the system is composed of two black holes (BHs), ⋆ Email: kun.emma@csfk.mta.hu with masses of 2.5 × 10 7 M ⊙ and 5.0 × 10 6 M ⊙ . The eccentric orbit (e ≈ 0.52) has a period of T = 418 d in the rest frame of the object. From our point of view, the orbital plane is seen almost edge-on, within an angle of ∼ 8 • .
Studying binary AGNs is an active field of both observational and theoretical astrophysics, due to its connection to cosmological structure formation, galaxy evolution, and most recently gravitational waves. Observations of such objects are very challenging (for a review see e.g. Komossa & Zensus 2016), and securely confirmed cases are extremely rare (De Rosa et al. 2019). Spikey stands out from the very few SMBHB candidates because the binary selflensing model (Hu et al. 2020) constrains the orbital parameters, the geometry, and the masses of the companions very accurately. The model also provides a testable prediction that the next flaring will occur in 2020.
Apart from being a moderately bright X-ray AGN (Hu et al. 2020), J1918+4937 is also a prominent radio-loud quasar. Variations in its ∼ 100-mJy level flux density at 15 GHz are being monitored at the Owens Valley Radio Observatory (OVRO, Richards et al. 2011). The source is known to have a compact radio jet structure at milliarcsec (mas) angular scales, as revealed by very long baseline interferometry (VLBI) imaging observations (e.g. Kovalev et al. 2007). Detecting binary AGNs separated by a small fraction of a pc is practically impossible with direct imaging observations, even with the high resolution offered by VLBI and in the most nearby universe (e.g. An et al. 2018). But radio interferometric observations could help in another way, by detecting a discernible effect of a binary companion on the appearance of the relativistic jet produced by the other AGN in the system, because the orbital motion of the jetted AGN may result in a helical shape of the jet (for a recent review, see De Rosa et al. 2019, and references therein).
A growing number of studies propose radio-loud AGNs as strong candidates for efficient high-energy (HE) neutrino emitter objects, especially the blazars (e.g. Kadler et al. 2016;IceCube Collaboration et al. 2018a;Garrappa et al. 2019;Giommi et al. 2020), whose jets point close to the line of sight of the observer. The underlying physical mechanisms involve light-matter and/or matter-matter interactions in a relativistically moving plasma. Kun et al. (2017Kun et al. ( , 2019 proposed a model in which the radio and neutrino observations were put into a common physical picture involving the spinflip of a SMBH in a merging binary. Recently the γ-ray flaring blazar TXS 0506+056, an efficient particle accelerator, turned out to be the source of several IceCube neutrinos (IceCube Collaboration et al. 2018a;IceCube Collaboration 2018b). Studies indicate that neutrino emission might be due to a recent merger activity (Britzen et al. 2019;Kun et al. 2019). Although there is no indication yet of an observed neutrino event near its position, Spikey is a VLBI source directing its jet close to our line of sight, and also a SMBH binary candidate, making this AGN an object of great interest as a potential neutrino source.
In this paper, we investigate whether the available radio data are consistent with the behaviour of the object, the gravitational self-lensing model, and in particular the parameters derived for Spikey by Hu et al. (2020). Based on archival data from 2008, we present sensitive and detailed VLBI images of J1918+4937 obtained at 1.7 and 5 GHz for the first time, and show the OVRO flux density curve (Sect. 2). By modeling the source brightness distribution at ∼ 1 − 10 mas scale, we derive parameters describing the relativistic jet, estimate the apparent speed of the jet based on snapshot VLBI observations conducted at the 8.4/8.7 GHz frequency band, and put forward a scenario where a jet is launched from one of the accreting BH components of the system in Sect. 3. Here we also investigate the case whether the OVRO flux density curve is compatible with the binary model. We discuss our findings based on the OVRO singledish and VLBI radio observations in Sect. 4. We also discuss whether Spikey could be an efficient high-energy neutrino emitter in the near future based on the behaviour of its VLBI jet and the proposed SMBH merger scenario. Finally we conclude the paper with a summary in Sect. 5.
We assume a flat ΛCDM cosmological model with H 0 =70 km s −1 Mpc −1 , Ω m = 0.3, and Ω Λ = 0.7 in this paper. In this model, an object at z sp = 0.926 has a luminosity distance of D L ≈ 6 Gpc, and 1 mas angular size corresponds to 7.853 pc projected linear size (Wright 2006 Kharb et al. (2010) scheduled the observations in 5min switching cycles, with 2 min spent on the calibrator J1918+4937 and 3 min on the weak target NGC 6764, including antenna slewing times. As a valuable byproduct of this phase-referencing experiment, ∼1.6 and 1.8 h of VLBA data accumulated on our object of interest, J1918+4937, at 1.7 and 5 GHz, respectively.
We downloaded the raw VLBA data of the BK154 experiment from the public archive 1 of the U.S. National Radio Astronomy Observatory (NRAO). For the data calibration, we used the NRAO Astronomical Image Processing System (AIPS, Greisen 2003) in a standard way (e.g. Diamond 1995). We started with calibrating the ionospheric delays based on total electron content measurements, and corrected for the measured Earth orientation parameters. We applied digital sampler corrections, then used a short 1-min scan on the bright fringe-finder source 1758+388 to solve for instrumental phases and delays. Bandpass correction was performed using the same scan. We used the gain curve and system temperature information from the participating VLBA stations for a-priori amplitude corrections. Finally fringefitting was done and the solutions were applied to the data.
The calibrated visibility data of J1918+4937 were exported to the Difmap software (Shepherd et al. 1994). After the standard hybrid mapping procedure involving several iterations of clean decomposition, phase-only selfcalibration, and finally phase and amplitude self-calibration, we obtained the naturally-weighted 1.7-and 5-GHz VLBI images of J1918+4937 shown in Fig. 1. As a finishing step of data reduction, we fitted circular Gaussian brightness distribution model components directly to the self-calibrated visibility data in Difmap. This allows us to describe the radio structure with a limited set of parameters that are listed in Table 1, where the errors were calculated as in Kun et al. (2014). These model components will be used for determining the shape of the jet. This way we can also gain information on the Doppler boosting of the relativistic jet.
Further VLBI imaging observations of J1918+4937 made at the 8.4/8.7-GHz frequency band are available in Figure 1. VLBA images of J1918+4937. Left: at 1.7 GHz. The peak intensity is 128 mJy beam −1 . The lowest contours are drawn at ±0.3 mJy beam −1 . The elliptical Gaussian restoring beam is 8.9 mas × 6.4 mas (FWHM) at a position angle −28 • . Right: at 5 GHz. The peak intensity is 119 mJy beam −1 . The lowest contours are drawn at ±0.28 mJy beam −1 . The restoring beam is 2.9 mas × 2.1 mas (FWHM) at a position angle −38 • . In both images, the positive contour levels increase by a factor of 2, and the restoring beam size is indicated in the bottom-left corner.  Table 2 and plotted as a function of time in Fig. 2. We also included here the 5-GHz data point (Table 1) because it was obtained at a close frequency and it helps filling the gap in the time coverage of the 8.4/8.7-GHz measurements. Although the data points have a scatter beyond the formal uncertainties due to the complications with imaging a complex structure from limited observations, the component separation clearly increases with 2 http://astrogeo.org/cgi-bin/imdb get source.csh?source=J1918%2B4937 time. The line in Fig. 2 indicates the linear fit for estimating the apparent angular proper motion, 0.10 ± 0.01 mas yr −1 .

Total flux density monitoring
The quasar J1918+4937 is included in the sample of extragalactic sources regularly monitored with the 40-m OVRO radio telescope at 15 GHz frequency (Richards et al. 2011). The total flux density variability of J1918+4937 from early 2008 to date can be seen in Fig. 3 which we constructed from the monitoring data available at the OVRO website 3 . A few flux density data points (∼ 9% of the total number) with excessively large error bars were discarded as unreliable.
The data after 2015 November 28 (2015.9) have error bars typically a factor of ∼ 2 smaller than before. This is likely due to the new receiver installed at OVRO in 2014 June, and a new data processing pipeline used. Keeping in mind that the time sampling of the flux density curve is more or less uniform, we re-scaled the error bars, in order to associate comparable weights to the older and the more recent data for a subsequent model fitting. The procedure was as follows. For the n-th data point, we calculated the standard deviation of the flux densities in the range [(n − k) . . . (n + k)] with k = 10, and assigned it to the n-th data point as its new error bar. This way, the first and last k data points had to be dropped from the light curve. The 15-GHz flux densities with the smoothed error bars are shown in Fig. 3, overlaid on the original light curve. Figure 1 shows an asymmetric radio structure with a compact core and a one-sided extension. It is typical for bright radio-loud quasars where the emission from one of the intrinsically symmetric jets that is pointing close to the observer's line of sight is enhanced by relativistic beaming (for a recent review, see Blandford et al. 2019). In the case of J1918+4937, the approaching jet is pointing towards the Northwest as projected on the sky. The radio emission can be traced out to about 80 mas at the lower observing frequency, 1.7 GHz, then it becomes diffuse and resolved out on the long interferometer baselines. This angular extent corresponds to a projected linear size of 630 pc.

Jet parameters
The bright VLBI core at the southeastern end of the nearly straight structure (Fig. 1) is in fact the base of jet where it becomes optically thick at the given observing frequency. The fitted Gaussian model parameters of the core (Table 1) can be used to calculate the apparent brightness temperature, where F is the flux density measured in Jy, θ the diameter of the circular Gaussian component in mas (full width at halfmaximum, FWHM), and ν the observing frequency in GHz.
The ratio between the apparent and the intrinsic brightness temperatures gives the Doppler-boosting factor, δ = T b /T b,int . If we follow the usual practice and assume the equipartition brightness temperature (Readhead 1994) as T b,int ≈ 5 × 10 10 K, then the Doppler factor is δ ≈ 5. On the other hand, based on measurements of a sample of pc-scale jets, Homan et al. (2006) arrived at a somewhat lower typical intrinsic brightness temperature value, T b,int ≈ 3× 10 10 K. Considering this, the Doppler factor of the jet in J1918+4937 would become δ ≈ 9.
The amount of Doppler boosting depends on two fundamental jet parameters, the bulk Lorentz factor (Γ) of the plasma flow (i.e. the intrinsic jet speed) and the jet inclination with respect to the line of sight (ι 0 ). If the apparent proper motion of the jet components can be measured based on VLBI imaging observations conducted at multiple epochs, it is possible to estimate values of Γ and ι 0 as well (e.g. Urry & Padovani 1995). Even though sensitive imaging data are found in the archives for J1918+4937 at a single epoch only at the above frequencies (1.7 and 5 GHz), from the available multi-epoch snapshot 8-GHz VLBI observations we were able to track the motion of one of the inner jet components. Assuming a linear outward motion (Fig. 2), we estimate its apparent speed in the units of the speed of light (c) as β app = 5.33 ± 0.65. If we consider β app as a representative estimate of the apparent jet speed in J1918+4937, and take the possible values of the Doppler factors derived above, we can obtain (see e.g. Urry & Padovani 1995) the bulk Lorentz factor Γ = β 2 app + δ 2 + 1 2δ (2) and the jet inclination angle For δ = 5, we get Γ ≈ 5.4 and ι 0 ≈ 11. • 5, and for δ = 9, we get Γ ≈ 6.1 and ι 0 ≈ 5. • 6.

Jet structural model utilizing eccentric SMBH orbit
While the jet shape in Fig. 1 seems remarkably straight on scales of several tens of mas, some wiggling is also apparent, especially at 5 GHz where the angular resolution is higher. Here we build up a structural (morphological) model of the jet as seen projected onto the plane of the sky, based on the fitted circular Gaussian model component positions (Table 1). We assume that these compact radio components were launched by a jetted supermassive black hole (SMBH) moving along an eccentric orbit in the binary system, and the jet launching is affected by the periodically changing orbital velocity of the jet emitter SMBH. This idea was applied earlier in several studies (Roos et al. 1993;Kun et al. 2014Kun et al. , 2015 but for circular orbits. Here we further develop the model, to allow for eccentric binary orbits with arbitrary spin angles. Note that in the jet model below, the jet components themselves move along ballistic trajectories and not along helical paths. Rather we see a helical pattern on the sky formed by the subsequently emitted components, as the angle of the jet launching changes periodically. We assume that this pattern motion preserves the jet launching angle at least up to tens of mas from the central engine. Meanwhile, the physical distances between the components are growing as the time passes. Let us assume an orthogonal coordinate system K in which the z axis is parallel to the orbital angular momentum L N = L NLN (z||L N ) (hereL N denotes the unit vector pointing to the direction of the orbital angular momentum), and the x axis is directed towards the pericentre of the orbit. The orbital configuration is depicted in Fig. 4. The instantaneous orbital velocity vector of the i-th BH in the orbital plane as a function of the eccentric anomaly E is where v 0,i is the circular orbital speed of the jet emitter SMBH (i = 1 for the dominant, and i = 2 for the secondarymass BH), is its true anomaly, is the semi-major axis of the orbit, G is the gravitational constant, T is the orbital period, m = m 1 + m 2 is the total mass, and e is the orbital eccentricity. If the dominant BH is the jet emitter, then its velocity should be considered in Eq. 4, which is and if the secondary BH is the jetted one, its velocity is The direction of the jetted BH spin S i in K is the unit vector where κ i = arccos(Ŝ i ·L N ) is the angle between S i and the orbital angular momentum L N , and ψ i is the angle between the projection of the spin onto the orbital (x, y) plane and the periapsis line. We assume that one of the two BHs emits the jet via the Blandford-Znajek mechanism (Blandford & Znajek 1977). In this case, the jet symmetry axis is directed along the BH spin S i , consequently the unperturbed jet velocity vector becomes v s = v sŜi in K, and its components are The jet velocity vector v jet is the vectorial sum of the unperturbed jet velocity vector v s and the orbital velocity v i , such that Let ζ be the angle between v jet and v s , which is calculated as vs (||S) Plane of the sky O r b i t a l p l a n e I Figure 4. Geometric configuration of the Spikey system centred on its barycentre. The black dot marks the position of the jetemitting SMBH along its elliptical orbit. LOS indicates the line of sight, and L N is the Newtonian orbital angular momentum. The true anomaly is χ, the argument of the periapsis is ω, the orbital inclination is I , the BH spin angle with respect to the orbital normal is κ 1 , the angle between the projection of the spin onto the orbital plane and the periapsis line is ψ, and the inclination angle of the spin with respect to the LOS is ι 0 . The position angle of the spin projected onto the plane of the sky (λ 0 ) is measured from North through East. Furthermore, v 1 is the orbital velocity vector of the jet-emitting SMBH at the instant of the jet component launching (if the secondary BH emits the jet, for its argument of periapsis ω 2 = ω 1 + π holds in radians), v s is the original jet velocity vector (that is parallel to the spin) and v jet is the vectorial sum of the above two. Finally, ζ is the instantaneous half-opening angle of the jet. For the sake of clarity, we shifted the jet velocity vector to the barycentre. In reality, the jet launches from the immediate vicinity of the emitting SMBH. where with C 1 = 1 + e 2 + 2e cos χ, with For the orbital velocities in Spikey, even at this sub-pc separation, v 0,i ≪ v s , and then the series expansion of their ratio (Eq. 12) in leading order gives Now let us define a new orthogonal coordinate system K ′ , such that its z ′ axis is parallel to the spin of the jetted BH. In this system, the jet morphological model turns to where u is the polar angle (Kun et al. 2014), φ is the initial phase of u, and B is the jet growth in mas perpendicular to its symmetry axis while u changes by 2π over the time period T u . This latter quantity is measured in the observer's frame as where v 0,i cos κ i is the orbital velocity perpendicular to S i , s is the scale factor that relates projected linear size to the measured angular size (in pc mas −1 ). Another parameter, A is the jet growth in mas parallel to its symmetry axis while u changes by 2π over the time period T u . The quantity A is measured in the observer's frame as Then we define a coordinate system K ′′ , such that the x ′′ and y ′′ axes point to East and North in the plane of the sky, respectively, and the z ′′ axis coincides with the direction of the line of sight (LOS), as shown in Fig. 4. The inclination angle between the LOS and the spin of the jet emitter BH is ι 0 (which we call spin inclination angle), and λ 0 is its position angle measured from North (y ′′ axis) through East (x ′′ axis). Employing the same rotational matrices as in Kun et al. (2014), In this model, the helical jet shape is in fact the pattern drawn by the perturbed jet components ejected at different epochs from the central engine. In other words, the individual components do not move along a helix, rather the pattern they collectively form grows both in the direction of the spin and perpendicular to it, as described by the parameters A and B, respectively.

Application of the jet model to the VLBI data
After setting up the model to describe the jet structure, we now take into account the measured VLBI component positions at 1.7 and 5 GHz (Table 1) Figure 5. The modeled jet shape fitted to the 1.7-and 5-GHz VLBI component positions, marked by blue triangles and red squares, respectively. Left: the right ascension and declination of the components relative to the VLBI core and the best-fit jet shape assuming Doppler factor δ = 5, once with the dominant-mass (m 1 ) SMBH being the jet emitter (dotted black curve) and then with the secondarymass (m 2 ) SMBH (continuous black curve). Middle: the right ascension of the components relative to the VLBI core as a function of u. Right: the declination of the components relative to the VLBI core as a function of u. The jet structure is rotated by 90 • towards East. Table 3. Grid parameters of the best-fit models, such as projected pitch along the spin (A ′′ ), Lorentz factor (Γ), spin angle (κ), and those derived from them, the spin inclination angle (ι 0 ) and the jet speed β, assuming that either the dominant-mass SMBH (top) or the secondary SMBH (bottom) launches the jet. We show the lowest χ-square values (χ 2 min ). We also list the parameter ranges in which the models are indistinguishable from each other (i.e. ∆AIC ≤ 2) and give the averages and standard deviations of the parameters in those ranges.
The jet is emitted by the larger mass SMBH (m 1 ) δ = 5 (χ 2 min = 34.54) δ = 9 (χ 2 min = 33.84) but at different frequencies, and the position of the optically thick core components (i.e., the base of the jet used as a reference for the relative position of other components further along the jet) is known to depend on the observing frequency, an effect called core shift (Blandford & Königl 1979;Lobanov 1998). Sokolovsky et al. (2011) conducted a dedicated survey with the VLBA at nine frequencies in the 1.4 − 15.4 GHz range to quantify the core-shift effect in 20 AGN jets. The average (and median) core shift between 1.7 and 5 GHz was found to be approximately 0.9 mas. This is comparable to the uncertainties of our component positions (Table 1). Therefore we used model components fitted at both 1.7 and 5 GHz together in the further analy-sis. Note that the angular resolution of the interferometer is about 3 times better at the higher frequency. Thus the inner section of the jet is characterised by more components at 5 GHz, while the outer section is only seen at 1.7 GHz, where the array is more sensitive to the weaker, extended, steep-spectrum features.
We are in a unique situation because some of the key parameters of the Spikey SMBHB system are accurately known from Hu et al. (2020). Therefore we adopt the (restframe) orbital period T = 1.14 yr, the orbital inclination I = 1.43 rad, the mass of the primary SMBH m 1 = 2.5 × 10 7 M ⊙ , the mass of the secondary SMBH m 2 = 5.0 × 10 6 M ⊙ , and the orbital eccentricity e = 0.52. These numbers imply a = 5.1 × 10 13 m = 0.002 pc and the circular velocities v 0,1 ≈ 0.006c and v 0,2 ≈ 0.03c. The bulk jet speed (expressed in the units of c) and the spin inclination angle with respect to the LOS are and respectively (e.g. Urry & Padovani 1995), where β s = v s /c. We apply non-linear least squares curve (parametric) fitting with σ −2 weights by employing the Levenberg-Marquardt algorithm to get the best-fit jet model, such that the χ 2 was minimized during the process. As a next step, we characterise the reliability of our best-fit model and investigate whether there are other solutions that cannot be discriminated from the above one, solely based on their χ 2 value. The Akaike information criterion (AIC, Akaike 1974) estimates the quality of each model relative to each of the others, i.e. it is a tool for model selection, either for nested or not nested models. The lower the AIC, the better the performance of the given model. Models in which the difference in AIC relative to AIC min is ≤ 2 perform approximately equally (Burnham & Anderson 2002), therefore the selection of any of them might lead to inconclusive statements. Here, as the number of parameters is the same, we select the models of approximately equal quality solely based on their χ 2 values.
If we apply the Doppler factor δ = 5 (see Sect. 3.1), then a lower limit for the Lorentz factor is Γ min = 2.6. This corresponds to the case when the jet is seen exactly poleon (i.e. ι 0 = 0). For a numerical parameter estimation, we set up a grid where the projected jet growth along the spin direction, A ′′ = A ′ sin ι 0 changes from 10 to 30 mas (in steps of ∆A ′′ = 0.1 mas), Γ changes from 3 to 20 (in steps of ∆Γ = 0.5), and κ i changes from 1 • to 90 • (in steps of ∆κ i = 1 • ). The bulk jet speed varies from 0.9428 c to 0.9987 c on the grid as Γ changes between 3 and 20. The only parameter we have to solve for is φ i , while A ′′ , Γ, and κ i are changing along the grid as described above. Since v 0,i ≪ v jet , we neglect the term corresponding to v 0,i in A ′ (Eq. 20), and the jet grows along the spin direction solely as a result of the non-zero jet velocity v jet .
By fitting the jet model described by Eqs. 21-22, the following best-fit parameters emerged if we assume the dominant-mass BH as the jetted one: A ′′ = 17.9 mas, Γ = 20.0, κ 1 = 0 • (with the lowest χ 2 = 34.54, reduced χ 2 R = 1.38). The modeled jet shape corresponding to these values and the measured VLBI component positions are plotted in Fig. 5. After considering the best-quality models leading to AIC difference from AIC min as ∆AIC ≤ 2, we calculate the average value and the standard deviation of the grid parameters. We repeated the process with Doppler factor δ = 9 (corresponding to T b,int = 3 × 10 10 K; see Sect. 3.1). Selecting the best-quality models, we get again the average value and the standard deviation of their grid parameters. The bestfit grid parameters, as well as parameters of models giving the same performance are summarized in Table 3 for δ = 5 and δ = 9. We also show here the spin inclination angles (ι 0 ) and jet speeds (β) derived from the corresponding grid parameters. It seems that the model is not very sensitive to the Lorentz factor, which is not surprising because the same projected jet opening angle can be generated with a variety of parameter pairs if we allow to simultaneously change the jet growth in the direction to the spin and perpendicular to it. Note that the best-fit jet structure model ( χ 2 = 34.54) is achieved with κ 1 = 0 • , and the parameter range of κ 1 giving models with comparable quality emerged as [0 • :30 • ]. The orbital velocity of the more massive SMBH is relatively small compared to the jet velocity along the spin because of the small BH mass ratio in Spikey. The fitting process tries to balance it with increasing the cos κ i term in Eq. 19 in order to model the observed jet growth perpendicular to the spin as closely as possible.
We repeated the jet-shape-fitting process, now assuming the m 2 mass SMBH as the jetted one. The grid parameters of the best-fit jet model, the parameter ranges in which the models lead to ∆AIC ≤ 2, as well as the averages and standard deviations of the parameters in these ranges are summarized in Table 3. The modeled jet shape corresponding to these values is plotted in Fig. 5. If the secondary-mass SMBH is assumed as the jet emitter, then the best-fit model gives κ 2 = 73 • if δ = 5, with a slightly lower χ 2 = 31.42 compared to the value found for the primary SMBH. The case is similar for δ = 9, and the corresponding parameter ranges leading to ∆AIC ≤ 2 are much more tightly constrained, without containing the limiting κ i = 0. This is because the velocity of the secondary SMBH is much larger compared to the more massive one, and the observed jet growth perpendicular to the spin can be modeled without maximizing the v 0,i cos κ i term in Eq. 19.

Total flux density variations
The observed period in the optical light curve of Spikey is T obs = (1 + z sp )T = 805 d, which was recognised as the observed orbital period of the SMBH binary (Hu et al. 2020). The 15-GHz radio flux density curve (Fig. 3, Sect. 2.2) measured at OVRO (Richards et al. 2011) indicates a decreasing trend on a long term, together with some flaring activity, and possibly a longer flare started at around 2015 November. If we interpret the flux density changes as quasi-periodic, a signal with 500 − 600-d period might seem superimposed on the linear trend. This period is ∼ 200−300 d shorter than the one in the optical light curve and therefore we could not reliably fit a periodic component by employing the proposed binary parameters of Spikey (Hu et al. 2020), where the periodically strengthened Doppler boosting would readily explain the radio flux density variations. The expected periodic effect is most likely masked by the episodic activity of the jetted AGN in the system. Instead we fitted a simple linear function to the smoothed data (see Sect. 2.2), resulting in a slope of (−4.67 ± 0.13) mJy yr −1 . This trend is also shown in Fig. 3. Also, we cannot exclude the possibility that some level of radio emission is associated with the second SMBH component.
The decreasing trend in the flux density curve might indicate that the average inclination angle of the jet becomes larger with time. In the framework of the SMBHB model, this can be interpreted as the jet direction gradually moving away from the line of sight, therefore decreasing the Doppler boosting effect on the observed radio emission. Below we in-vestigate whether this scenario is consistent with the known binary parameters (Hu et al. 2020).
The orbital period in the order of years and the subpc separation in Spikey indicate that the binary has already progressed into the inspiral evolutionary phase of the merger, i.e., the third and final stage where the gravitational radiation becomes the dominant dissipative effect over dynamical friction and gravitational slingshot interactions (e.g. Merritt & Milosavljević 2005). In the inspiral phase, the dynamical evolution of the binary can be treated analytically by expanding the equation of motion in terms of the socalled post-Newtonian (PN) parameter as follows (Kidder 1995): (25) where r is the binary separation being in the coordinate system K. Here ε = Gmc −2 r −1 is the PN parameter with r = a(1 − e 2 )(1 + e cos χ) −1 , and O(ε n ) represents the n-th PN order. For the eccentric orbit in Spikey, we average the PN parameter for one orbit: which value suggests Spikey recently entered into its inspiral phase, where 0.001 ε 0.1 (Gergely & Biermann 2009;Levin et al. 2011).
Up to 2PN orders, the merger dynamics is conservative, the constants of motion being the total energy and the total angular momentum vector J = S 1 + S 2 + L, where L is the total orbital angular momentum. The SMBH spins obey precessional motion (Barker & O'Connell 1975: where the i index refers to the first or second component of the binary. The angular velocity Ω i of the i-th spin S i contains up to 2PN order spin-orbit (1.5PN), spin-spin (2PN), and quadrupole momentum contributions (2PN). For the typical mass ratios ν ∈ [1/30 . . . 1/3], only the dominant spin counts (Gergely & Biermann 2009). The mass ratio in Spikey is ν ≈ 1/5, so it falls into the above range implying the second spin might be neglected in the binary dynamics. In 1.5PN, the spin-orbit precession of the spins S 1 and S 2 occurs with angular velocities respectively, where L N = µr × v is the Newtonian orbital angular momentum, µ = m 1 m 2 /m is the reduced mass which moves with velocity v. Employing the formulae of the instantaneous separation given in Eq. 26 and the orbital velocity vector given in Eq. 4 (both expressed in K) The time dependence of r, χ, and E can be given by solving the Kepler equation E(t) − e sin E(t) = 2π/T(t − τ), where τ is the time of pericentre passage. Substituting Eq. 31 into Eqs. 29-30, and averaging the spin-orbit precession period T SO = 2πΩ over one orbit, we get a value for the dominant-mass SMBH as T SO,1 (t) (1 + z sp ) ≈ 15, 700 yr and for the secondary SMBH as T SO,2 (t) (1 + z sp ) ≈ 3, 800 yr in the observer's frame.
Assuming that the bulk Lorentz factor in the jet remains constant with time, and the long-term decreasing trend in the OVRO flux density curve (Fig. 3) is solely due to the secular change of the spin inclination angle, we calculate the possible jet inclination angles at two epochs of the OVRO flux density monitoring period by employing the flux density ratio below: where the indices 1, 2 mark the flux density and spin inclination angle at two arbitrary epochs. We assumed a flat radio spectrum. In Fig. 3, we marked three different epochs, t A , t B , and t C , which are the starting epoch of the smoothed OVRO flux density curve, the epoch of the 1.7-and 5-GHz VLBA observations, and the last epoch of the smoothed OVRO flux density curve, respectively. The mean 15-GHz flux densities at these three epochs are F(t A ) = 132 mJy, F(t B ) = 130 mJy, and F(t C ) = 81 mJy, respectively, based on the (−4.67 ± 0.13) mJy yr −1 slope of the linear function fitted to the flux density data. Employing the minimum and maximum spin inclination angles allowed by the VLBI measurements at epoch t B in the framework of the present binary model, ι 0,min = 3. • 7 (with Γ = 5, δ = 9) and ι 0,max = 11. • 5 (with Γ = 5, δ = 5), and the flux density ratio in Eq. 32, we calculate minimum and maximum spin inclination angles at the starting and finishing OVRO epochs. The resulting possible spin inclination angles are summarized in Table 4. According to our results, the spin inclination angle could have changed by 2. • 5-2. • 6 over 11 yr in the framework of the present model. By expanding the equation of motion in terms of the PN parameter, as we have seen, the dynamical evolution of the binary can be treated analytically while it progresses through the inspiral phase where 0.001 ε 0.1 (Gergely & Biermann 2009;Levin et al. 2011). The time scale of the spin-flip is proportional to ε −9/2 , while the time scale of the spin-orbit precession is proportional to ε −5/2 , when the spin is comparable to the orbital angular momentum (S 1 ≈ L). For Spikey, ε(t) ≈ 0.001 means that if the flip occurs, it happens on a time scale more than 10 6 times longer than the time scale of the precession. We can safely state that if the slow decrease in the total flux density of Spikey is indeed due to the increase of the spin inclination angle, then the underlying mechanism should be the spinorbit precession, not the spin-flip.

No spike in the radio light curve
The long-term 15-GHz OVRO monitoring (Richards et al. 2011) covers the time of the Kepler spike (Smith et al. 2018) occured in 2011 June. Since the optical flare lasted only for about 15 days, it was poorly sampled in the radio. However, Table 4. Possible jet inclination angles at the beginning (2008.22) and at the end (2019.22) of the OVRO flux density curve, based on the estimated minimum and maximum jet inclination values at the epoch of VLBI observations (2008.87). If we assume the minimum (maximum) inclination angle at 2008.87, the inclination angle changes from 3. • 6 (10. • 9) to 6. • 2 (13. • 4) in 11 years of the OVRO observations. These angles are marked by boldface (italic) in the table, respectively. there are 3 measurement points available in the OVRO data set for J1918+4937 in this time range, roughly at the beginning, middle, and end of the optical spike. From these data, there is no evidence for any radio brightening around Julian Date 2455724. On the contrary, the 15-GHz flux density stays constant within the measurement errors. Why is the radio emission unaffected in the SMBHB self-lensing scenario that Hu et al. (2020) proposed for the optical spike? There are two possible reasons. First of all, if only one of the BHs powers a radio jet, and this one is the lensing object in the foreground, then a radio magnification is obviously not expected. But even if the lensed object in the background is a radio-loud AGN, an optical spike is not necessarily expected to be coupled with a radio brightening. The optical emission of AGNs is known to originate mainly from the accretion disk on the scale of ∼ 10 −5 pc (e.g. Koratkar & Blaes 1999). On the other hand, most of the 15-GHz radio emission comes from an ultracompact region downstream the jet, on ∼ 0.1 − 1 pc projected scale (e.g. Lobanov 1998). However, according to the model of Hu et al. (2020), the SMBHB separation in the Spikey system is at least two orders of magnitude smaller. The entire binary system is therefore located well inside the region where the 15-GHz radio emission originates from. There is nothing to be gravitationally lensed in the Spikey system in radio, and even the superior angular resolution of VLBI is insufficient to directly resolve the companions.

Jet modeling with accurate binary parameters
Modeling the observed high-resolution structure and kinematics of VLBI jets in quasars is usually applied to infer parameters of suspected SMBH binaries (e.g. Britzen et al. 2001;Lobanov & Roland 2005;Britzen et al. 2012;Valtonen & Wiik 2012;Caproni et al. 2013;Kun et al. 2014Kun et al. , 2015Kun et al. , 2018. In some of these cases, there is independent indication for the existence of the binary, e.g. from periodic optical variability. However, in the case of J1918+4937 (Spikey), the analysis of the Kepler light curve by Hu et al. (2020) offers more than simply an indication. The measured optical spike is a unique phenomenon requiring special circumstances, therefore its successful modeling with gravitational self-lensing and orbital Doppler boosting provides accurately determined BH masses, orbital parameters and geometric constraints for the system (Hu et al. 2020). Unlike the usual practice, these parameters could therefore be fed directly into the VLBI jet model presented here (Sect. 3).
It was necessary to refine this model to allow for highly eccentric binary orbits. In all earlier modeling, circular orbits were assumed for simplicity, as no reliable information about the binary orbital parameters were available.
We used VLBI imaging data taken at 1.7 and 5 GHz frequencies for Spikey, and also investigated the long-term OVRO flux density monitoring measurements at 15 GHz in the context of the SMBHB model proposed by Hu et al. (2020). The shape of the VLBI jet represented by the individual component positions is remarkably consistent with the Spikey binary parameters. The constraints we obtained on Γ based on the single-epoch deep VLBI imaging of J1918+4937 at these two frequencies are not particularly strong (see Table 3). Indeed, qualitatively, a jet with a given Doppler boosting factor can be produced either by relatively slowly-moving plasma blobs directed very close to the line of sight, or a fast jet with comparably larger inclination. Plausible values of Γ and the mean jet inclination angle with respect to the LOS (ι 0 ) could be provided only with multiepoch VLBI jet kinematic studies (e.g. Lister et al. 2019).
However, utilizing also the available multi-epoch snapshot VLBI imaging observations of J1918+4937 at the 8.4/8.7-GHz frequency band, we were able to estimate the apparent speed (β app ) in the jet from the measured linear proper motion of an inner jet component. The values of Γ and ι 0 derived from β app for the two possible values of the Doppler factor (δ = 5 and 9) fall within the parameter ranges set by our VLBI jet stucture model based on the parameters of the orbital motion of a SMBHB along eccentric orbit (Hu et al. 2020). Moreover, the values estimated from jet kinematics, Γ ≈ 5 − 6 and ι 0 ≈ 6 − 12 • , appear more consistent with the solutions in Table 3, where the jet emitter is the smaller BH with mass m 2 .
As the optical emission likely arises from the gas bounded to the individual SMBHs in the binary system, the luminosity of the brighter minidisk (e.g. Ryan & MacFadyen al. 2017) would be Doppler boosted and this minidisk is likely the one associated with the fastestmoving secondary SMBH (D'Orazio et al. 2016;Hu et al. 2020). The spike in the Kepler optical light curve of Spikey can be explained with the gravitational self-lensing if the larger-mass SMBH passes between the smaller-mass SMBH and the observer (Hu et al. 2020), magnifying the optical emission of the minidisk around the smaller SMBH. Also, Hu et al. (2020) successfully explained the long-term variability in the light curve of Spikey by variable Doppler boosting due to the motion of the secondary SMBH. This means that at least the smaller SMBH has an accretion disk what we see in optical. Our VLBI jet model that utilizes the binary model of Hu et al. (2020) is indeed more consistent with the jet parameters derived from VLBI monitoring of Spikey if we assume the secondary SMBH is the jetted one in the system. Notably, the χ 2 values also indicate slightly better fits for those solutions, and if the secondary-mass SMBH is the jetted one, the parameter κ i is much better constrained, without reaching the limiting value 0 • .
As it is often seen in radio-loud AGNs, the OVRO monitoring light curve of Spikey (Fig. 3) is rather complex. Variations with characteristic time scales of ∼ 1 yr and shorter are superimposed on a generally decreasing trend in flux density. We attempted to relate this long-term trend seen during the entire monitoring period of more than 11 yr to the SMBHB model in which one of the companions launches the relativistic jet and is therefore responsible for the synchrotron radio emission. Spin-orbit precession in a close binary SMBH system that is already in its inspiral phase can cause a change in the orientation of the BH spin and the jet orientation. Considering the Spikey parameters, we found that this change (about 2. • 5-2. • 6 during ∼ 11 yr) should have a noticeable effect during the OVRO monitoring period by driving away the jet from the line of sight and thus decreasing the Doppler boosting, effectively causing the observed gradual dimming of the radio source.
Based on our study, we can confidently say that the Spikey jet and the radio light curve are fully consistent with the binary SMBH model of Hu et al. (2020). Both the jet shape and the long-term decreasing flux density trend can be reconciled with the proposed binary parameters and standard jet physics. However, alternative explanations cannot be excluded for the observed VLBI jet pattern, as well as for the radio light curve. Precessing jets can also be produced by tilted accretion discs around rapidly spinning BHs (Liska et al. 2018), without invoking the presence of a binary companion. Any periodic or quasi-periodic effect related to the jet itself, its surrounding medium, or the jet feeding mechanism can in principle affect its observed structure. For example, plasma instabilities along the jet (e.g. Nakamura & Meier 2004) and quasi-periodic instabilities in the accretion flow (e.g. Honma et al. 1992) can also cause wiggled jet structures. Similarly, total flux density variations can be produced by a multitude of physical effects, not only the change in the jet inclination angle. In particular, a longterm change in the bulk jet Lorentz factor could result in a similar trend seen in Fig. 3. The main point of why the SMBHB scenario is the most favourable one to explain the GHz VLBI jet structure of Spikey is that we already have an indication that Spikey hides a SMBHB based on the gravitational self-lensing model of Hu et al. (2020) and the spike seen in the optical light curve of the object.

Could Spikey become a neutrino emitter
AGN? Kun et al. (2017Kun et al. ( , 2019 proposed a scenario of binary SMBH evolution which naturally explains the observed high-energy (HE) neutrino emission, and leads to the emission of gravitational waves (GWs) through a sequence induced by the merger. For the typical mass ratio of merging SMBH binaries (ν ∈ [1/3 . . . 1/30]), L > S 1 is always transformed into L < S 1 (Gergely & Biermann 2009). It means the spin of the dominant BH usually flips, while spin-orbit precessing. Three main phases of the emission of HE particles are expected in this scenario. The first one is the process of spinflip, when the jets sweep through a large cone. The second one is after the spin-flip, when a new jet is boring into the environment, leading to the injection of more seed particles to create HE nuclei, γ-photons and neutrinos. The third one is probably in the instant of the coalescence of SMBHs, when a giant shock wave may be generated by low-frequency GWs to accelerate particles to high energies, leading to a final burst of HE nuclei, γ-rays and neutrinos.
To speculate if Spikey could be a neutrino emitter based on the available data, it is vital to conclude which spin the jet is connected to. The jet power (P jet ) is proportional to the mass of the central object (m) and the square of its dimensionless spin parameter a * (e.g. Narayan & McClintock 2012;Steiner et al. 2013). We have seen that the optical light curve and the VLBI observations of Spikey together are slightly more consistent with the secondary SMBH being the jetted one in the system. The ratio of the spin magnitudes in Spikey is which means the spin of the secondary SMBH might not be neglected in the binary dynamics only if its horizon rotates much faster compared to the horizon of the dominant one, i.e. if a * 2 ≫ a * 1 holds for the dimensionless spin parameters. If it is the case, then the jet power would be much larger if the secondary SMBH emits the jet, because P jet ∝ (m, a * 2 ). So the physical picture in Spikey becomes self-consistent if the horizon of the dominant-mass SMBH rotates much slower compared to the secondary SMBH. In this case, S 2 could be in the order of S 1 , and eventually flip in the inspiral phase.

SUMMARY
J1918+4937 (Spikey) is so far a unique extragalactic object hosting a closely-separated (∼ 0.002 pc) SMBHB system where the masses of the companions, as well as the orbital and geometric parameters could be accurately determined from a narrow spike in its Kepler optical light curve, using a combined gravitational self-lensing and orbital Doppler boosting model (Hu et al. 2020). At least one of the SMBH companions is a radio-loud AGN with a prominent relativistic plasma jet. Archival high-resolution radio interferometric imaging observations made with the VLBA at 1.7 and 5 GHz (Kharb et al. 2010) allowed us to study its structure. We estimated the Doppler boosting caused by the small inclination angle of the jet to the line of sight. We then set up a model describing a jetted SMBH in a binary system with eccentric orbit, and investigated whether the apparently helical jet shape is consistent with the binary parameters derived for Spikey (Hu et al. 2020). By successfully applying our structural model to Spikey, we could derive the jet Lorentz factor and viewing angle, albeit with loose constraints. A comparison with the jet parameters inferred from multi-epoch VLBI monitoring data at 8.4/8.7 GHz, together with the somewhat better fits provided by the jet structural model suggest that the smaller-mass (m 2 ) component of the binary might be the jet-emitting BH.
We also studied the long-term single-dish 15-GHz flux density curve (Richards et al. 2011). While spikes similar to the optical one are not expected in the radio, the long-term behaviour of light curve may bear the imprint of a close binary companion to the radio-loud AGN. Indeed, the gradually decreasing trend is consistent with the expected spinorbit precession which slowly increases the viewing angle of the jet.
Recent developments in extragalactic neutrino astronomy suggest that AGN with jets inclined close to our line of sight might be strong sources of the high-energy neutrinos reconstructed in the IceCube Neutrino Detector. Based on the properties of its VLBI jet, the binary parameters proposed by Hu et al. (2020), and the merger-induced neutrino emission scenario proposed by Kun et al. (2017Kun et al. ( , 2019, we found that Spikey could become an efficient high-energy neutrino source if the horizon of the secondary SMBH is rapidly rotating. While the observed VLBI jet structure and the longterm trend in the flux density monitoring could possibly be explained with other effects, the consistency of both types of measurements with the Spikey binary parameters is remarkable, and can be considered as a support for the model of Hu et al. (2020). The jet parameters could be determined with higher confidence and our values confirmed in the future with further frequent sensitive multi-epoch VLBI imaging observations. Our jet structural model involving eccentric orbit can later be applied for similar binary candidate AGNs with a jetted companion.

DATA AVAILABILITY
The datasets underlying this article were derived from sources in the public domain as given in the respective footnotes.