Signatures of a spinning supermassive black hole binary on the mas-scale jet of the quasar S5 1928+738 based on 25 yr of VLBI data

ABSTRACT

In a previous work, we have identified the spin of the dominant black hole of a binary from its jet properties. Analysing Very Long Baseline Array (VLBA) observations of the quasar S5 1928+738, taken at 15-GHz during 43 epochs between 1995.96 and 2013.06, we showed that the inclination angle variation of the inner (<2 mas) jet symmetry axis naturally decomposes into a periodic and a monotonic contribution. The former emerges due to the Keplerian orbital evolution, while the latter is interpreted as the signature of the spin-orbit precession of the jet emitting black hole. In this paper, we revisit the analysis of the quasar S5 1928+738 by including new 15-GHz VLBA observations extending over 29 additional epochs, between 2013.34 and 2020.89. The extended data set confirms our previous findings which are further supported by the flux density variation of the jet. By applying an enhanced jet precession model that can handle arbitrary spin orientations κ with respect to the orbital angular momentum of a binary supermassive black hole system, we estimate the binary mass ratio as ν = 0.21 ± 0.04 for κ = 0 (i.e. when the spin direction is perpendicular to the orbital plane) and as ν = 0.32 ± 0.07 for κ = π/2 (i.e. when the spin lies in the orbital plane). We estimate more precisely the spin precession velocity, halving its uncertainty from |$(-0.05\pm 0.02)$| to |$(-0.04\pm 0.01)^{\circ }\, \mathrm{yr}^{-1}$|⁠.

techniques: interferometric, galaxies: active, galaxies: jets, quasars: individual: S5 1928+738, quasars: supermassive black holes

1 INTRODUCTION

When supermassive black holes (SMBHs), located in the centre of large galaxies merge as their host galaxies undergo a major merger (e.g. Begelman, Blandford & Rees 1980), first the dynamical friction, then the gravitational radiation governs the dynamical evolution of their orbit (e.g. Binney & Tremaine 1987; Merritt & Milosavljević 2005). Once the dissipation of energy and angular momentum becomes dominated by gravitational radiation, the already sub-pc separated binary goes through three main phases: inspiral, plunge, and ringdown (e.g. Kidder 1995). Post-Newtonian (PN) techniques are involved in the description of the binary dynamics in the inspiral phase for which the PN parameter is constrained as 0.001 < ε < 0.1 (e.g. Gergely & Biermann 2009). In the plunge stage, numerical calculations are required to characterize the final SMBH (e.g. Buonanno & Damour 2000). The ringdown is characterized in terms of quasi-normal modes of black holes (Berti, Cardoso & Starinets 2009).

A large number of archival studies identified indirect signals of sub-pc separated SMBH binaries in radio-loud active galactic nuclei (AGNs), including the recent works (e.g. Eracleous et al. 2012; Graham et al. 2015; Charisi et al. 2016; Liu et al. 2016; Chen et al. 2020). SMBH binaries are also emerging in the multimessenger discussion of AGNs, mostly of blazars (e.g. Britzen et al. 2019, 2023, 2021; Kun, Biermann & Gergely 2019; de Bruijn et al. 2020; ecker Tjus et al. 2022; Broschewski, Becker Tjus & Biermann 2023).

The source we study here, S5 1928+738 (redshift z_spect = 0.302; Lawrence et al. 1986) is a low-synchrotron-peaked, low polarization quasar (⁠|$\mathrm{RA}=19^\mathrm{h} 27^\mathrm{m} 48{_{.}^{\rm s}}495$|⁠, Dec. = +73°58′1|${_{.}^{\prime\prime}}$|570), with an elongated jet structure, extending over ∼20 mas at 15 GHz. Three decades ago, the presence of an SMBH binary was already suspected at the centre of this quasar (Roos, Kaastra & Hummel 1993). Based on the long-term data of the jet of the quasar S5 1928+738 from the Monitoring Of Jets in Active galactic nuclei with VLBA Experiments (MOJAVE) program, obtained with the Very Long Baseline Array (VLBA), we further strengthened the idea that this quasar harbours an SMBH binary (Kun et al. 2014). In that work, we have shown the spin-precession of the dominant black hole. We monitored the direction of the symmetry axis of the inner 2 milliarcseconds (mas) of the jet of S5 1928+738 through its inclination and position angles. By fitting the shape of a projected helical jet to the position of the VLBA-detected jet components, we derived the value of the inclination and position angles. Then we decomposed the observed time variation of the inclination of the symmetry axis into a periodic term with an amplitude of ∼0.°89 and a linear decreasing trend with a rate of |$\sim 0.05^{\circ }\, \mathrm{yr}^{-1}$|⁠.

We explained the periodic component in the direction of the inner jet as generated by the SMBH binary with total mass m ≈ 8.13 × 10⁸ M_⊙ (m = m₁ + m₂, m₂ < m₁). This periodic change in the direction of the symmetry axis of the jet is induced by the orbital motion of the jet emitter dominant black hole. We derived the binary parameters consistent with the VLBA measurements spanning through 18 yr: orbital period T = (4.78 ± 0.14) yr, separation r = (0.0128 ± 0.0003) pc and PN parameter ε ≈ 0.003, lower limit of the mass ratio ν ≳ 1/5 (ν = m₂/m₁). The upper limit of ν could not be estimated from the observational data. Therefore, we adopted the upper limit of the mass ratios typical for SMBH mergers, ν ≲ 1/3 (Gergely & Biermann 2009).

Next, we identified the slow, secular reorientation of the jet (superposed to the periodic component) as due to the spin-orbit precession of the jet emitter SMBH. From the observations, we derived an upper limit for the precession period of the more massive black hole (T_SO < 5500 yr) and the gravitational time-scale of the merger (T_GW < 1.64 × 10⁶ yr). These results apply for zero spin angle (κ = 0), i.e. assuming the spin of the jet emitter SMBH perpendicular to the orbital plane.

Periodic jet structures are attributed either to the orbital motion of supermassive black hole binaries (Hummel et al. 1992; Roos et al. 1993; Britzen et al. 2001; Deane et al. 2014), or instabilities in the plasma. Our model presented in Kun et al. (2014) introduces a scenario in which the orbital motion of the jet emitter black hole shakes up an intrinsically helical jet (for visualization, see fig. 9. in Kun et al. 2018).

In this paper, we extend the time range of the 15-GHz VLBI observations of the jet of S5 1928+738, by adding new VLBA observations taken in 29 epochs between 2013.34 and 2020.89. In Section 2, we describe the analysis of these MOJAVE data. In Section 3, we investigate the overall jet morphology and kinematics between 1995.96 and 2020.89. We study the structure of the inner 2 mas of the jet of S5 1928+738 to derive the variation of the inclination angle of the inner jet symmetry axis. In Section 4, we develop a new binary SMBH jet model enabling arbitrary spin angles. In Section 5, we apply this model to the inclination angle variation of the inner jet axis of S5 1928+738 and constrain the mass ratio as function of the spin angle κ, generalizing our earlier results. In this section, we also present a sanity check with the 15-GHz flux densities. Finally, in Section 6, we discuss and summarize our results. We adopt the cosmological parameters as follows: H₀ = 67.8 km s⁻¹ Mpc⁻¹, Ω_m = 0.308, and |$\Omega _\Lambda =0.692$|⁠.

2 VLBA DATA ANALYSIS

We analysed calibrated visibilities in 29 epochs between 2013.34 and 2020.89 taken with VLBA at 15 GHz observing frequency, provided by the MOJAVE program (Lister et al. 2009). All ten VLBA stations participated in the observations at 19 epochs, while single antennas were not available in 10 epochs (see Table 1). MOJAVE publishes amplitude- and phase-calibrated visibility data and we performed the imaging and model fitting by employing the difmap software (Shepherd 1997). All non-flagged data were used. The summary of the image parameters of the observing epochs between 2013.34 and 2020.89 for S5 1928+738 is given in Table 1.

Table 1.

Summary of the image parameters of the observing epochs between 2013.34 and 2019.65 for S5 1928+738. (1) Date of observation, (2) VLBA project code, (3) participating VLBA stations BR – Brewster, FD – Fort Davis, HN – Hancock, KP – Kitt Peak, LA – Los Alamos, MK – Mauna Kea, NL  – North Liberty, OV – Owens Valley, PT – Pie Town, SC – St. Croix (‘All’ means all telescopes are participating, and e.g. ‘–BR’ means ‘All but Brewster’), (4) peak brightness, (5) total flux density, (6) root mean square (rms) of the brightness in the residual image, (7) major axis of the elliptical Gaussian restoring beam (FWHM), (8) minor axis of the restoring beam (FWHM), (9) position angle of the major axis of the restoring beam, (10) number of jet components, (11) reduced χ² of the model fit. The full table is available in electronic format online.

Date(1)	Code(2)	Array(3)	|$S_\mathrm{p}^{(4)}$|	|$S_\mathrm{tot}^{(5)}$|	rms(6)	|$b_\mathrm{maj}^{(7)}$|	|$b_\mathrm{min}^{(8)}$|	|$b_\mathrm{pa}^{(9)}$|	N(10)	|${\chi ^2_\mathrm{red}}^{(11)}$|
			(Jy beam−1)	(Jy)	(mJy beam−1)	(mas)	(mas)	(°)
2013-05-05	BL178BC	All	2.115	3.405	0.17	0.751	0.715	−83.4	12	1.09
2013-08-12	BL178BH	−FD	2.300	3.626	0.20	0.647	0.588	29.0	11	1.07
2013-12-15	BL193AA	All	2.848	4.086	0.09	0.835	0.704	43.4	11	1.27
2014-09-01	BL193AN	All	3.675	4.878	0.12	0.653	0.621	−35.1	12	1.36
2015-06-16	BL193AT	All	5.199	6.350	0.10	0.702	0.628	−22.5	10	1.92

To minimize the degrees of freedom of the model fitting, we used circular Gaussian components only to describe the surface brightness distribution of the jet. The fitted parameters were the total flux density, the position, and the full width at half maximum (FWHM) size of the components. The errors were calculated as in Kun et al. (2014) and Kun et al. (2015). The components were identified by requiring their appearance in at least 5 subsequent epochs and smooth changes of their fitted parameters (especially the positions and flux densities). As in Kun et al. (2014), the jet components are labelled with letter ‘C’ and a number, such that the larger number denotes a component closer to the VLBI core. This way it was straightforward to label new jet components, that were ejected after the last observing epoch considered in our earlier work (Kun et al. 2014). The results of the model fits for the epochs between 2013.34 and 2020.89 are listed in Table 2.

Table 2.

Circular Gaussian model fit results for S5 1928+738. (1) epoch of observation, (2) flux density, (3) and (4) relative right ascension and declination of the component centre with respect to the core, (5) component diameter (FWHM), (6) jet component identifier. The full table is available in electronic format online.

Epoch(1)	S(2)	x(3)	y(4)	d(5)	ID(6)
(yr)	(Jy)	(mas)	(mas)	(mas)
2013.34	0.241 ± 0.009	0.000 ± 0.019	0.000 ± 0.018	0.023 ± 0.009	CS
	1.855 ± 0.025	0.080 ± 0.020	−0.240 ± 0.038	0.128 ± 0.001	Cg
	0.505 ± 0.013	0.306 ± 0.021	−0.774 ± 0.020	0.169 ± 0.004	C16
	0.346 ± 0.011	0.363 ± 0.021	−1.057 ± 0.020	0.182 ± 0.006	C15
	0.145 ± 0.007	0.646 ± 0.024	−1.709 ± 0.023	0.288 ± 0.017	C14
	0.054 ± 0.005	1.074 ± 0.039	−2.217 ± 0.039	0.694 ± 0.069	C13
	0.100 ± 0.006	1.584 ± 0.032	−3.149 ± 0.032	0.529 ± 0.030	C11
	0.080 ± 0.006	1.375 ± 0.051	−3.800 ± 0.050	0.938 ± 0.066	C9
	0.037 ± 0.004	1.517 ± 0.059	−4.946 ± 0.059	1.127 ± 0.187	C8
	0.013 ± 0.003	1.442 ± 0.065	−6.747 ± 0.065	1.240 ± 0.326	C7
	0.013 ± 0.003	2.469 ± 0.082	−9.089 ± 0.082	1.595 ± 0.855	C4
	0.016 ± 0.006	2.393 ± 0.141	−17.060 ± 0.141	2.789 ± 2.973	C1

In its core region at 15 GHz (e.g. Kun et al. 2014) and 43 GHz observing frequencies (e.g. Lister & Smith 2000), S5 1928+738 includes two components ‘CS’ and ‘Cg’. As in our earlier work (Kun et al. 2014), the component CS is assumed to be the VLBI core.

3 MORPHOLOGY, KINEMATICS, AND JET GEOMETRY

3.1 Morphology and proper motions

The jet structure of S5 1928+738 is very rich (see fig. 1 of Kun et al. 2014) and the jet itself extends up to 18 mas projected radial distance at 15 GHz observing frequency, measured from the VLBI core. In Fig. 1 we plot the core separation of the jet components as a function of time. Additionally to the two quasi-stationary components CS and Cg, we identified 20 outward-moving components.

Figure 1.

Core separation of the jet components as a function of time with error bars and their fitted proper motion represented by continuous lines. Accelerating components are marked by ⋆. The separation of the jet components is relative to the VLBI core ‘CS’ at (0,0).

We calculated the proper motion of the components first by fitting linear proper motions as

$$\begin{eqnarray} r_\mathrm{L}=\mu (t-t_\mathrm{ej}), \end{eqnarray}$$

(1)

and then by fitting accelerating proper motions as (e.g. Homan et al. 2001; Lister et al. 2009)

$$\begin{eqnarray} r_\mathrm{A}=r_\mathrm{L} + \frac{\dot{\mu }}{2} \left(t-t_\mathrm{mid}\right)^2, \end{eqnarray}$$

(2)

where μ is the linear proper motion measured in mas yr⁻¹, t_ej the ejection time, |$\dot{\mu }$| the angular acceleration, and t_mid the half of the sum of the maximum and the minimum epochs when the respective components are detected. The proper motion of the components can be converted to apparent speed in the units of the speed of light, such that

$$\begin{eqnarray} \beta _\mathrm{app}=0.0158 \frac{\mu D_L}{(1+z_\mathrm{spect})}, \end{eqnarray}$$

(3)

where μ is the proper motion, D_L the luminosity distance in Mpc and z_spect the redshift. We summarize the linear and accelerated proper motions in Table 3.

Table 3.

Component speed fits. (1) jet component identifier. From linear fit (L): (2) linear proper motion, (3) ejection time, (4) χ². From accelerated fit (A): (5) non-linear proper motion, (6) acceleration rate, (7) ejection time, (8) χ². (9) Difference between the AICs (Δ = AIC_lin − AIC_acc), (10) type of the proper motion based on the AIC (L – linear, A – accelerated), see details in the text, (11) apparent speed.

ID(1)	|$\mu _\mathrm{L}^{(2)}$|	|$t_\mathrm{ej,L}^{(3)}$|	χ2(4)	|$\mu _\mathrm{A}^{(5)}$|	|${\dot{\mu }_\mathrm{A}}^{(6)}$|	|$t_\mathrm{ej,A}^{(7)}$|	χ2(8)	Δ(9)	PM(10)	βapp,lin(11)
	(mas yr−1)	(yr)		(mas yr−1)	(mas yr−2)	(yr)				(c)
C1	0.305 ± 0.010	1959.5 ± 1.4	270.2	0.312 ± 0.009	0.011 ± 0.003	1961.4 ± 1.3	261.0	7.2	L	5.97 ± 0.20
C2	0.373 ± 0.032	1972.3 ± 2.3	102.1	0.375 ± 0.034	−0.011 ± 0.037	1972.4 ± 2.3	101.6	−1.5	L	7.30 ± 0.63
C3	0.438 ± 0.026	1985.5 ± 0.8	192.1	0.407 ± 0.024	−0.066 ± 0.021	1984.2 ± 0.9	187.3	2.8	L	8.57 ± 0.51
C4	0.315 ± 0.005	1985.2 ± 0.3	593.4	0.325 ± 0.006	0.006 ± 0.002	1986.4 ± 0.5	512.8	78.6	A	6.16 ± 0.10
C5	0.392 ± 0.019	1990.2 ± 0.4	146.2	0.407 ± 0.020	0.058 ± 0.030	1990.8 ± 0.4	119.6	24.6	A	7.67 ± 0.37
C6	0.328 ± 0.013	1990.0 ± 0.3	162.3	0.343 ± 0.012	0.038 ± 0.010	1990.7 ± 0.3	104.0	56.3	A	6.42 ± 0.25
C7	0.287 ± 0.004	1989.9 ± 0.2	144.1	0.292 ± 0.004	0.005 ± 0.001	1990.3 ± 0.2	105.6	36.5	A	5.62 ± 0.08
C8	0.263 ± 0.002	1993.5 ± 0.1	760.2	0.263 ± 0.002	−0.001 ± 0.001	1993.4 ± 0.2	758.1	0.1	L	5.15 ± 0.04
C9	0.224 ± 0.003	1994.3 ± 0.2	328.1	0.225 ± 0.003	0.003 ± 0.001	1994.8 ± 0.3	278.2	47.9	A	4.38 ± 0.06
C10	0.202 ± 0.006	1996.5 ± 0.2	22.8	0.198 ± 0.006	−0.013 ± 0.005	1996.3 ± 0.2	16.8	4.0	L	3.96 ± 0.12
C11	0.341 ± 0.005	2002.8 ± 0.2	523.9	0.346 ± 0.003	0.010 ± 0.002	2003.3 ± 0.1	247.7	274.2	A	6.71 ± 0.08
C12	0.180 ± 0.009	2000.7 ± 0.3	3.2	0.179 ± 0.009	0.013 ± 0.016	2000.8 ± 0.3	2.7	−1.5	L	3.52 ± 0.18
C13	0.367 ± 0.010	2005.9 ± 0.3	1303.1	0.375 ± 0.009	0.018 ± 0.005	2006.5 ± 0.3	995.3	305.8	A	7.18 ± 0.18
C14	0.266 ± 0.009	2006.4 ± 0.3	989.3	0.267 ± 0.009	−0.003 ± 0.005	2006.4 ± 0.3	985.7	1.6	L	5.22 ± 0.16
C15	0.267 ± 0.010	2009.3 ± 0.2	705.4	0.268 ± 0.007	0.023 ± 0.005	2009.6 ± 0.2	329.9	373.6	A	5.32 ± 0.20
C16	0.209 ± 0.009	2009.8 ± 0.3	619.3	0.208 ± 0.009	0.011 ± 0.007	2009.9 ± 0.3	591.6	25.7	A	4.05 ± 0.16
C17	0.140 ± 0.017	2010.7 ± 0.9	484.6	0.127 ± 0.015	−0.071 ± 0.023	2009.5 ± 1.1	350.3	132.3	A	2.60 ± 0.27
C18	0.125 ± 0.006	2011.8 ± 0.3	120.5	0.127 ± 0.006	−0.009 ± 0.008	2011.8 ± 0.3	96.3	22.2	L	2.31 ± 0.10
C19	0.067 ± 0.012	2011.2 ± 1.5	18.0	0.068 ± 0.012	−0.020 ± 0.030	2011.2 ± 1.5	15.8	0.1	L	1.70 ± 0.16
C20	0.097 ± 0.013	1994.1 ± 1.1	0.2	0.096 ± 0.006	0.076 ± 0.013	1994.3 ± 0.5	0.0	−1.8	L	2.49 ± 0.33

In order to decide whether linear or accelerated motion gives the better fit to the proper motion of jet components, we applied the Akaike Information Criterion (AIC) method (Akaike 1974). AIC is useful when model comparison by exact statistical tests gives inconclusive results, or for nested models. The linear proper motion (equation 1) is nested into the accelerated proper motion (equation 2). First we calculated AIC = χ² + 2N (where N is the number of fitted parameters in a model) for the linear (AIC_lin) and acceleration model (AIC_acc) from component to component, then compared them (Δ = AIC_lin − AIC_acc). The following thresholds are usually employed in the literature (see e.g. in Burnham & Anderson 2003): Δ ≤ 2 indicates approximately comparable performances, 4 ≤ Δ ≤ 7 indicates a measurable difference in the fits, while Δ > 10 clearly favours one fit over the other. From these definitions, we only use Δ > 10 to distinguish between the models. After the calculation of Δ, ten components emerged as linearly moving components (C1, 2, 3, 8, 10, 12, 14, 18, 19, 20), while ten other components (C4, 5, 6, 7, 9, 11, 13, 15, 16, 17) showed signs of acceleration.

Applying different core separation upper cuts (r ≤ 2 mas, r ≤ 3 mas, r ≤ 4 mas, r ≤ 5 mas), we recalculated the Δ values between the performances of the linear and accelerating proper motion fits for the accelerating jet components, that start their observed trajectory in the innermost jet region. We show the results in Table 4. It seems that, for most of the components, increasing the cut-off core separation makes the linear proper motion fits worse compared to the accelerating proper motion fits. In the core separation region r ≤ 2 mas, for four components out of six, the linear proper motion fit gives better performance (Δ ≤ 10). This might be a sign of an acceleration zone somewhere between 2 and 3 mas, therefore we adopt the core separation cut of 2 mas from our earlier work (Kun et al. 2014) to define the inner jet in S5 1928+738, where jet precession is supposed to be mostly unaltered by external forces.

Table 4.

Differences in the AIC_lin and AIC_acc for the components (Δ = AIC_lin − AIC_acc), that show accelerating motion, in different regions of the jet.

ID	Δr ≤ 2 mas	Δr ≤ 3 mas	Δr ≤ 4 mas	Δr ≤ 5 mas
C9	−1.5	−2.0	29.2	14.2
C11	−1.1	28.3	98.4	239.0
C13	2.8	99.9	698.9	424.5
C15	2.7	151.4	373.6	373.6
C16	42.8	25.7	25.7	25.7
C17	132.3	132.3	132.3	132.3

The apparent speed of the jet components is related to the bulk jet speed β and their inclination angle ι, as (Urry & Padovani 1995)

$$\begin{eqnarray} \beta _\mathrm{app}=\frac{\beta \sin \iota }{1-\beta \cos \iota }, \end{eqnarray}$$

(4)

where the jet speed β is related to the Lorentz factor γ as

$$\begin{eqnarray} \beta =\sqrt{1-\frac{1}{\gamma ^2}}. \end{eqnarray}$$

(5)

The maximal β_app occurs for the critical inclination cos ι_c = β, sin ι_c = γ⁻¹, having the value β_{app, max} = (γ² − 1)^1/2.

In practice, when analysing the jet components, it is unlikely that any of them will have exactly the direction cos ι_c = β, sin ι_c = γ⁻¹. Hence, the jet component with maximal speed |$\beta _\mathrm{app,max} ^\mathrm{obs}\le \beta _\mathrm{app,max}$| leads to a lower limit

$$\begin{eqnarray} \gamma _\mathrm{low}=\sqrt{1+(\beta _\mathrm{app,max}^\mathrm{obs})^2} \end{eqnarray}$$

(6)

of the Lorentz factor of the jet.

The highest linear proper motion for the time period between 2013.34 and 2020.89, μ_max = (0.44 ± 0.03) mas yr⁻¹, is exhibited by component C3. This is in good agreement with our earlier results for a time period between 1995.96 and 2013.06, when we found C3 as the fastest component with the linear proper motion μ_max = (0.43 ± 0.02) mas yr⁻¹. With redshift z_spect = 0.302, luminosity distance D_L = 1612.6 Mpc, this proper motion translates to |$\beta _\mathrm{app,max} ^\mathrm{obs}=8.60\pm 0.51$|⁠. For the minimum Lorentz factor given by equation (6), we obtain γ_low ≈ 8.66 and the minimum intrinsic jet speed is β_low ≈ 0.992 in the units of the speed of the light. Then the related critical inclination angle is |$\iota _\mathrm{c}=\arcsin (\gamma _\mathrm{low}^{-1})\approx 7^{\circ }$|⁠.

3.2 Evolution of the orientation of the inner jet symmetry axis between 1995.96 and 2020.89

In Kun et al. (2014), we described the shape of the inner jet (core separation <2 mas) of the 15-GHz VLBI jet of S5 1928+738 as a helical structure, with the following geometry, given in a system |$\mathcal {K}_\mathrm{j}$| with the z_j coordinate along the jet symmetry axis:

$$\begin{eqnarray} &&{x_\mathrm{j}(u) = \frac{b}{2\pi }u\cos u,} \\ &&{y_\mathrm{j}(u) = \frac{b}{2\pi }u\sin u,} \\ &&{z_\mathrm{j}(u) = \frac{a}{2\pi }u,} \end{eqnarray}$$

(7)

where u is the parameter of the helical structure, with a and b ≪ a being the axial and radial growth rates of the jet per turn, respectively. This shape appears in projection in the plane of the sky. An adapted coordinate system |$\mathcal {K}$| has its x and y axes pointing towards north and east, respectively, with its z axis along the line of sight (LOS). The symmetry axis of the helical jet has polar angle ι₀ (the inclination angle measured from the LOS, z-axis) and azimuthal angle λ₀ (measured from the north, x-axis). For a jet pointing towards us, ι₀ ∈ [0, π/2].

The coordinates in |$\mathcal {K}_\mathrm{j}$| are related to those in |$\mathcal {K}$| as follows:

$$\begin{eqnarray} {\begin{pmatrix}x \\ y \\ z \\ \end{pmatrix}} = R_z(\lambda _0) R_y(\iota _0) {\begin{pmatrix}x_\mathrm{j} \\ y_\mathrm{j} \\ z_\mathrm{j} \\ \end{pmatrix}}, \end{eqnarray}$$

(8)

where the rotation matrices are

$$\begin{eqnarray} R_y(\iota {_0}) \!=\! {\begin{pmatrix}\cos \iota {_0} &\quad 0 &\quad \sin \iota {_0} \\ 0 &\quad 1 & 0 \\ -\sin \iota {_0} &\quad 0 & \quad\cos \iota {_0} \end{pmatrix}} \end{eqnarray}$$

(9)

and

$$\begin{eqnarray} R_z (\lambda {_0}) = {\begin{pmatrix}\cos \lambda {_0} &\quad -\sin \lambda {_0} &\quad 0 \\ \sin \lambda {_0} &\quad \cos \lambda {_0} &\quad 0 \\ 0 &\quad 0 &\quad 1 \end{pmatrix}}. \end{eqnarray}$$

(10)

Hence, the coordinates of the helical structure in the plane of sky, employing equation (7) are

$$\begin{eqnarray} {x}(u)=F(u)\cos \lambda {_0}-G(u)\sin (\lambda {_0}), \end{eqnarray}$$

(11)

$$\begin{eqnarray} {y}(u)=F(u)\sin \lambda {_0}+G(u)\cos (\lambda {_0}), \end{eqnarray}$$

(12)

where

$$\begin{eqnarray} F(u)=u \left(\frac{b}{2 \pi } \cos (u)\cos (\iota {_0})+\frac{a}{2\pi }\sin (\iota {_0}) \right), \end{eqnarray}$$

(13)

$$\begin{eqnarray} G(u)=\frac{b}{2\pi } u\sin (u). \end{eqnarray}$$

(14)

Then the inclination angle and the position angle of the helical jet in the source frame arise as

$$\begin{eqnarray} \iota (u)=\arcsin {\sqrt{x(u)^2+y(u)^2}}, \end{eqnarray}$$

(15)

$$\begin{eqnarray} \lambda (u)=\arctan \left({y(u)}/{x(u)}\right). \end{eqnarray}$$

(16)

The intrinsic half-opening angle of the jet cone (ψ₀) can be given by the radial and axial increments as

$$\begin{eqnarray} \tan \psi _0=\frac{\sqrt{{x_\mathrm{j}}^2+{y_\mathrm{j}}^2}}{z_\mathrm{j}}=\frac{b}{a}. \end{eqnarray}$$

(17)

The apparent half-opening angle of the cone ψ^app is obtained through the following reasoning. The generating lines of the cone

$$\begin{eqnarray} \hat{G}_{\mathrm{j}\pm }= {\begin{pmatrix}\sin \psi _0 \cos (\lambda _0 \pm \frac{\pi }{2}) \\ \sin \psi _0 \sin (\lambda _0 \pm \frac{\pi }{2}) \\ \cos \psi _0 \\ \end{pmatrix}}, \end{eqnarray}$$

(18)

which lie in the plane perpendicular to the (x_j, y_j) projection of the LOS are transformed into the system |$\mathcal {K}$| system through |$\hat{G}_\pm = R_y(\iota {_0}) \hat{G}_{\mathrm{j}\pm }$|⁠, obtaining

$$\begin{eqnarray} \hat{G}^{x}_\pm = \sin \iota _0 \cos \psi _0 +\cos \left(\lambda _0\pm \frac{\pi }{2}\right)\cos \iota _0 \sin \psi _0, \end{eqnarray}$$

(19)

$$\begin{eqnarray} \hat{G}^{y}_\pm = \sin \left(\lambda _0 \pm \frac{\pi }{2} \right) \sin \psi _0, \end{eqnarray}$$

(20)

$$\begin{eqnarray} \hat{G}^{z}_\pm =\cos \iota _0 \cos \psi _0-\cos \left(\lambda _0\pm \frac{\pi }{2}\right)\sin \iota _0 \sin \psi _0. \end{eqnarray}$$

(21)

Note that the scalar product

$$\begin{eqnarray} \hat{G}_{\mathrm{j}+} \cdot \hat{G}_{\mathrm{j}-}=\cos 2\psi _0 \end{eqnarray}$$

(22)

gives the intrinsic opening angle, while the projections of |$\hat{G}_\pm$| to the (x, y) plane span the apparent opening angle 2ψ^app according to

$$\begin{eqnarray} \hat{G}^x_+ \hat{G}^x_- + \hat{G}^y_+ \hat{G}^y_-=\cos 2\psi ^\mathrm{app}. \end{eqnarray}$$

(23)

The apparent half-opening angle of the jet cone ψ^app is related to its intrinsic value ψ₀ as (0° < ι₀ < 90°, ψ₀ < ι₀)

$$\begin{eqnarray} \tan \psi ^\mathrm{app}=\frac{ \tan \psi _0}{\sin \iota _0} \left(1-\frac{\tan ^2 \psi _0}{\tan ^2 \iota _0}\right)^{-1/2}. \end{eqnarray}$$

(24)

For small values of ψ₀ the above equation simplifies (in radians) as

$$\begin{eqnarray} \psi ^\mathrm{obs} \approx \frac{\psi _0}{\sin \iota _0}. \end{eqnarray}$$

(25)

In Kun et al. (2014), we derived the parameters a = (10.6 ± 0.4) mas, b = (0.18 ± 0.06) mas and the initial values ι₀ = 7.0° and λ₀ = (160 ± 2)° (χ²_red = 0.31, D = 11). Here and throughout the paper, we used the Levenberg–Marquardt algorithm for fitting purposes, such that the χ² was minimized during the fitting. Then keeping a and b as fixed parameters, we fitted the time variation of ι₀ and λ₀ to the 15-GHz VLBI jet components of the inner jet of S5 1928+738 in 43 epochs between 1995.96 and 2013.06. The best-fitting results showed that the variation of the inclination of the inner 2 mas of the jet of S5 1928+738 is best described with a periodic term with an amplitude of ∼0.89° and a linear decreasing trend with a rate of |$\sim 0.05^{\circ }\, \mathrm{yr}^{-1}$| (see in Table 5).

We repeated the fit of the projected helical jet to the (x, y) positions of the inner components of the 15-GHz VLBI jet of S5 1928+738 employing new 15-GHz VLBA observations in 29 additional epochs between 2013.34 and 2020.89. We present the best-fitting results for ι₀ in Table 6. We show the phenomenological fit to the inclination angle change of the symmetry axis of the 15-GHz inner jet of S5 1928+738 with a periodic and a monotonic term (see in Kun et al. 2014) described as

$$\begin{eqnarray} \iota _0 (t)=A_0\!+\!A_1\sin \left(\frac{2\pi }{T}t\!-\!\phi \right)\!+\!A_2 t \end{eqnarray}$$

(26)

in Fig. 2, where χ² = 37.8, and the degrees of freedom D = 67. We show the previous phenomenological fit that we got based on the inclination angle variation until 2013.06 (χ² = 20.9, D = 38) and the prediction based on it, together with the extended time series and the fit on it. In Table 5, we present the best-fitting parameters of the phenomenological fit that we acquired in Kun et al. (2014) and the best-fitting values of the same parameters that we got by fitting the extended time-series of the inclination angle variation of the inner jet symmetry axis of S5 1928+738 (a total of 72 epochs between 1995.96 and 2020.89). The new fitted parameters agree within the errors of the previous parameters obtained in Kun et al. (2014). Employing the extended inclination angle curve, we are able to constrain the errors on the fitted parameters more tightly. Most importantly, the relative error on the linear trend decreased from 45 to 22 per cent, which makes it possible to constrain the velocity of the spin precession more precisely.

Figure 2.

Phenomenological fit of the inclination angle variation of the inner jet symmetry axis of S5 1928+738 at 15 GHz. Data points already published in Kun et al. (2014; squares with error bars), the phenomenological fit to them (continuous line between 1995.96 and 2013.34, χ ²_red = 0.55, D = 38, where D is the degrees of freedom), the prediction based on it (dotted line), the new data points (circles with error bars) and the phenomenological fit to the whole data set (continuous line across the whole plot, χ ²_red = 0.56, D = 67) are shown.

Table 5.

The parameters of the analytical fits to the time series of the inclination angles from Kun et al. (2014) and for the extended data set. These parameters are the average (A₀), the amplitude of the variation (A₁), monotonic term (A₂), period of the variation (T) and initial phase (ϕ), everything in the framework of the observer.

ι0(t), 1995.96–2013.06 (Kun et al. 2014), χ2red = 0.55
A0 [°]	A1 [°]	A2 [° yr−1]	T [yr]	ϕ [°]
8.67 ± 0.19	0.89 ± 0.17	−0.05 ± 0.02	6.22 ± 0.19	11.95 ± 15.01
ι0(t), 1995.96–2020.89, extended, χ2red = 0.56
A0 [°]	A1 [°]	A2 [° yr−1]	T [yr]	ϕ [°]
8.64 ± 0.16	0.72 ± 0.11	−0.04 ± 0.01	5.98 ± 0.10	24.63 ± 16.1

Table 6.

The flux density of the 15 GHz jet (F′) and inclination of the jet axis (ι₀) obtained from the parametric fit on the 15-GHz data between 2013.34 and 2020.89. The full table is available in electronic format online.

Epoch	F′(Jy)	ι0(°)
2013.34	3.405 ± 0.035	9.81 ± 1.04
2013.61	3.626 ± 0.042	9.50 ± 1.11
2013.95	4.037 ± 0.029	9.30 ± 1.03
2014.66	4.878 ± 0.036	7.65 ± 0.84
2015.45	6.350 ± 0.037	7.70 ± 0.85

4 THE GENERAL JET EMITTER BINARY BLACK HOLE MODEL

4.1 Newtonian order: the effect of orbital motion on the jet

In our model, the jet velocity is the vectorial sum of the intrinsic jet velocity and the orbital velocity of the jet emitter black hole in the instant of the jet emission. In the following we define these vectors and derive the jet velocity.

To describe geometry of our jet model, we define an orthonormal basis |$\mathcal {K}_\mathrm{LN}$|⁠, such that the direction of the Laplace–Runge–Lenz vector |$\mathbf {\hat{A}_N}$| and |$\mathbf {\hat{Q}_N}=\mathbf {\hat{L}_N}\times \mathbf {\hat{A}_N}$| form the orbital plane together with |$\mathbf {\hat{L}_N}$| (where |$\mathbf {\hat{L}_N}$| is the direction of the Newtonian orbital angular momentum, perpendicular to the orbital plane). We show the geometric configuration of the black hole binary system fixed on its barycentre in Fig. 3.

Figure 3.

Geometric configuration of the black hole binary system centred on its barycentre. The black dot marks the position of the jet-emitting SMBH along its elliptical orbit. LOS indicates the line of sight, L_N is the Newtonian orbital angular momentum, A_N is the Laplace–Runge–Lenz vector, and |$\mathbf {\hat{Q}_N}=\mathbf {\hat{L}_N}\times \mathbf {\hat{A}_N}$|⁠. The true anomaly χ is measured from |$\mathbf {\hat{A}_N}$|⁠, the argument of the periapsis is ω, the orbital inclination is I, the spin angle of the dominant SMBH (which is the jet emitter) with respect to the orbital normal is κ₁, the angle between the projection of the dominant spin on to the orbital plane and the periapsis line is ψ₁ and the inclination angle of the spin with respect to the LOS is ϕ₁. The position angle of the spin projected onto the plane of the sky (τ₁) is measured from north through east. Furthermore, v₁ is the orbital velocity vector of the dominant SMBH at the instant of the jet component launching, v_s is the original jet velocity vector (that is parallel to the spin), and v_j is the vectorial sum of these two. Finally, ζ is the instantaneous half-opening angle of the jet.

The radial vector of the reduced mass in the orbital plane is given as

$$\begin{eqnarray} \mathbf {r} = {\begin{pmatrix}r\cos \chi \\ r\sin \chi \\ 0 \end{pmatrix}}. \end{eqnarray}$$

(27)

Then orbital velocity of the reduced mass in a black hole binary is derived as

$$\begin{eqnarray} \mathbf {v} = {\begin{pmatrix}\dot{r} \cos {\chi }-r \dot{\chi } \sin \chi \\ \dot{r} \sin {\chi }+r \dot{\chi } \cos \chi \\ 0 \end{pmatrix}}, \end{eqnarray}$$

(28)

with

$$\begin{eqnarray} \dot{\chi }=\frac{L_\mathrm{N}}{\mu r^2}, \end{eqnarray}$$

(29)

$$\begin{eqnarray} \dot{r}=\frac{A_\mathrm{N}}{L_\mathrm{N}} \sin \chi , \end{eqnarray}$$

(30)

$$\begin{eqnarray} r=\frac{p}{1+e\cos \chi }, \end{eqnarray}$$

(31)

$$\begin{eqnarray} p=\frac{L_\mathrm{N}^2}{Gm\mu ^2}, \end{eqnarray}$$

(32)

$$\begin{eqnarray} e=\frac{A_\mathrm{N}}{Gm \mu }, \end{eqnarray}$$

(33)

here L_N is the length of the Newtonian orbital momentum vector L_N, A_N is the length of the Laplace–Runge–Lenz vector A_N, χ is the true anomaly (measured from the direction |$\mathbf {\hat{A}_\mathrm{N}}$| in the plane of motion), e is the eccentricity, G is the gravitational constant, m = m₁ + m₂ (m₁ > m₂) is the total mass, and μ = m₁m₂/m is the reduced mass.

Then the velocity of the reduced mass in |$\mathcal {K}_\mathrm{LN}$| becomes

$$\begin{eqnarray} \mathbf {v} = \frac{G m \mu }{L_\mathrm{N}} {\begin{pmatrix}-\sin \chi \\ e+\cos \chi \\ 0 \end{pmatrix}}, \end{eqnarray}$$

(34)

where E is the eccentric anomaly that is related to the true anomaly χ as

$$\begin{eqnarray} \chi =2\arctan \left[\sqrt{\frac{1+e}{1-e}}\tan \frac{E}{2}\right]. \end{eqnarray}$$

(35)

The Kepler equation gives the time dependence of the eccentric anomaly as

$$\begin{eqnarray} E-e\sin E=n t, \end{eqnarray}$$

(36)

with the mean anomaly n = 2π/T.

The two orbiting bodies form a gravitationally bound system, in which the bodies orbit under the forces in 0PN order (Newton’s second law)

$$\begin{eqnarray} m_1 \ddot{\mathbf {r_1}}=\mathbf {F_{12}}, \end{eqnarray}$$

(37)

$$\begin{eqnarray} m_2 \ddot{\mathbf {r_2}}=\mathbf {F_{21}}. \end{eqnarray}$$

(38)

Since F₁₂ and F₂₁ are equal in length and opposite in direction (F₁₂ + F₂₁ = 0, r₁ − r₂ = r), one can derive the position of the orbiting bodies as

$$\begin{eqnarray} \mathbf {r_1}=\frac{m_2}{m} \mathbf {r}, \end{eqnarray}$$

(39)

$$\begin{eqnarray} \mathbf {r_2}=-\frac{m_1}{m} \mathbf {r}, \end{eqnarray}$$

(40)

where r₁ and r₂ point to the locations of m₁ and m₂ from the position of the centre of mass that travels with constant velocity (r₁ = R₁ − R_CM and r₂ = R₂ − R_CM, where R_CM is the location of the centre of mass, while R₁ and R₂ are vectors describing the motion of m₁ and m₂). From this, the velocity of the two bodies are

$$\begin{eqnarray} \mathbf {v_1} = \frac{m_2}{m} \mathbf {v}, \end{eqnarray}$$

(41)

$$\begin{eqnarray} \mathbf {v_2} = -\frac{m_1}{m} \mathbf {v}, \end{eqnarray}$$

(42)

where v is the relative velocity between v₁ and v₂. Then, from equation (34)

$$\begin{eqnarray} \mathbf {v_1} = \frac{G m \mu }{L_\mathrm{N}}\frac{\nu }{(1+\nu)} {\begin{pmatrix}-\sin \chi \\ e+\cos \chi \\ 0 \end{pmatrix}}, \end{eqnarray}$$

(43)

$$\begin{eqnarray} \mathbf {v_2} = - \frac{G m \mu }{L_\mathrm{N}}\frac{1}{(1+\nu)} {\begin{pmatrix}-\sin \chi \\ e+\cos \chi \\ 0 \end{pmatrix}}, \end{eqnarray}$$

(44)

As mentioned above, the jet velocity vector v_j in |$\mathcal {K}_\mathrm{LN}$| is the vector sum of v_s and the orbital velocity of the jet emitter black hole v₁ (both in |$\mathcal {K}_\mathrm{LN}$|⁠), i.e. v_j = v_s + v₁. Now we derive v_j.

The jet power (P_j) is predicted to depend on the spin of the black hole as P_j ≈ (SB)² (Blandford & Znajek 1977), where S = aM is the spin, a is the dimensionless spin parameter, M is the mass), and B is the magnetic field strength at the horizon of the black hole. Indeed, observations of radio-loud AGN show that the accretion disc alone is not enough to power up the relativistic jets of AGNs (Ghisellini et al. 2014). It is plausible to assume SMBHs in a tight binary system are co-evolving in a similar environment acquiring similar dimensionless spin parameters. In case of similar rotation velocities the ratio of the black hole spins, therefore the ratio of the jet power in a double-jetted binary system depends mostly on the mass ratio of the black holes. Since we do not see the signs of another jet in S5 1928+738 that would indicate the presence of a comparable powered pair of jets, we assume the dominant SMBH is the jet emitter black hole, similarly to our earlier work (Kun et al. 2014). In the |$\mathcal {K}$| coordinate system (attached to the LOS, see in Section 3.2), the inclination and position angle of the dominant spin S₁ are ϕ₁ and τ₁ (see Fig. 3). In 0PN order, the direction of the dominant spin does not change with time |$\dot{\phi }_1^{0PN}=\dot{\tau }_1^{0PN}=0$|⁠. In |$\mathcal {K}_\mathrm{LN}$|⁠, the spin of the jet emitter BH S₁ points to the direction

$$\begin{eqnarray} \mathbf {S_\mathrm{1}}= {\begin{pmatrix}v_\mathrm{s} \sin \kappa _1 \cos \psi _1\\ v_\mathrm{s} \sin \kappa _1 \sin \psi _1\\ v_\mathrm{s} \cos \kappa _1 \end{pmatrix}}, \end{eqnarray}$$

(45)

with the polar angle |$\kappa _1= \arccos (\mathbf {\hat{S}_1} \cdot \mathbf {\hat{L}_N})$| measured between S₁ and L_N and azimuthal angle ψ₁ measured between the projection of the dominant spin S₁ onto the orbital plane and the periapsis line. Assuming that the dominant SMBH emits the jet via the Blandford–Znajek mechanism (Blandford & Znajek 1977), the unperturbed jet velocity is then |$\mathbf {v_\mathrm{s}}=v_\mathrm{s} \mathbf {\hat{S}_1}$|⁠.

The jet velocity vector v_j in |$\mathcal {K}_\mathrm{LN}$| is the vector sum of v_s and the orbital velocity v₁ (both in |$\mathcal {K}_\mathrm{LN}$|⁠), such that

$$\begin{eqnarray} \mathbf {v_\mathrm{j}}(\mathcal {K}_\mathrm{LN})= {\begin{pmatrix}v_\mathrm{s} \sin \kappa _1 \cos \psi _1-v_{1} \sin \chi \\ v_\mathrm{s} \sin \kappa _1 \sin \psi _1+v_{1} (e+\cos \chi)\\ v_\mathrm{s} \cos \kappa _1 \end{pmatrix}}, \end{eqnarray}$$

(46)

where

$$\begin{eqnarray} v_1=\frac{G m \mu }{L_\mathrm{N}}\left(\frac{\nu }{1+\nu }\right). \end{eqnarray}$$

(47)

Then the half-opening angle of the precession cone induced by the orbital motion in the source frame ζ can be given as

$$\begin{eqnarray} \cos \zeta &=& \left[ \frac{\mathbf {v_{j}^2+\mathbf {v_\mathrm{s}}^2}-\mathbf {v_1}^2}{2 {v_\mathrm{s}} {v_\mathrm{j}}} \right] \\ &=&\frac{1+X f(\kappa ,e,\chi ,\psi _1)}{((1+e^2)X^2+1+2 e X^2 \cos \chi +2X f(\kappa ,e,\chi ,\psi _1))^{1/2}},\\ \end{eqnarray}$$

(48)

with

$$\begin{eqnarray} X=\frac{v_1}{v_\mathrm{s}} \end{eqnarray}$$

(49)

and

$$\begin{eqnarray} f(\kappa _1, e,\chi ,\psi _1)=\sin \kappa _1 \left[e \sin \psi _1-\sin {(\chi -\psi _1)}\right]. \end{eqnarray}$$

(50)

Since ζ and χ correspond to the polar and azimuthal angles of the vector v_j, respectively, in a spherical polar system, we can write down the jet velocity given in |$\mathcal {K}_\mathrm{j}$| as

$$\begin{eqnarray} \mathbf {v_\mathrm{j}}(\mathcal {K}_\mathrm{j})=v_\mathrm{j} {\begin{pmatrix}-\sin \zeta \sin \chi \\ \sin \zeta (e+\cos \chi)\\ \cos \zeta \end{pmatrix}}, \end{eqnarray}$$

(51)

where ζ is the instantaneous half-angle of the precession cone and χ is the instantaneous true anomaly.

The jet velocity in |$\mathcal {K}_\mathrm{j}$| is related to that in |$\mathcal {K}$| as follows:

$$\begin{eqnarray} {\begin{pmatrix}v_{\mathrm{j},x}^{\prime \prime } \\ v_{\mathrm{j},y}^{\prime \prime } \\ v_{\mathrm{j},z}^{\prime \prime } \\ \end{pmatrix}} = R_z(\tau _1) R_y(\phi _1) {\begin{pmatrix}v_{\mathrm{j},x} \\ v_{\mathrm{j},y} \\ v_{\mathrm{j},z} \\ \end{pmatrix}}, \end{eqnarray}$$

(52)

where the forms of R_y(ϕ₁) and R_z(τ₁) are given by equation (9) and equation (10), respectively. Then the jet velocity projected in the plane of the sky becomes

$$\begin{eqnarray} v_{\mathrm{j},x}^{\prime \prime }={v_{\mathrm{j},x}^{\prime }}\cos \tau _1^\mathrm{0PN}- {v_{\mathrm{j},y}^{\prime }}\sin \tau _1^\mathrm{0PN}, \end{eqnarray}$$

(53)

$$\begin{eqnarray} v_{\mathrm{j},y}^{\prime \prime }={v_{\mathrm{j},x}^{\prime }}\sin \tau _1^\mathrm{0PN}+ {v_{\mathrm{j},y}^{\prime }}\cos \tau _1^\mathrm{0PN}, \end{eqnarray}$$

(54)

 where

$$\begin{eqnarray} v_{\mathrm{j},x}^{\prime }=-v_\mathrm{j}\sin \zeta \sin \chi \cos \phi _1^\mathrm{0PN}+v_\mathrm{j}\cos \zeta \sin \phi _1^\mathrm{0PN}, \end{eqnarray}$$

(55)

$$\begin{eqnarray} v_{\mathrm{j},y}^{\prime }=v_\mathrm{j} \sin \zeta (e+\cos \chi). \end{eqnarray}$$

(56)

Then the inclination angle and the position angle of the jet ejection in 0PN order:

$$\begin{eqnarray} \iota _\mathrm{j}^\mathrm{0PN}=\arcsin {\sqrt{{v_\mathrm{jx}^{\prime \prime }}^2+{v_\mathrm{jy}^{\prime \prime }}^2}}, \end{eqnarray}$$

(57)

$$\begin{eqnarray} \lambda _\mathrm{j}^\mathrm{0PN}=\arctan ({v_{\mathrm{j},y}^{\prime \prime }}/{v_{\mathrm{j},x}^{\prime \prime }}). \end{eqnarray}$$

(58)

Making a series expansion in ζ, keeping only the leading order (since the orbital velocity of the jet emitter black hole is much smaller then the jet velocity, and ζ should be small), we get the inclination angle of the jet ejection from equation (57) (for |$\phi _1^{0PN}\gt 0$|⁠) in the source frame as

$$\begin{eqnarray} \sin \iota _\mathrm{j}^{0PN}=\sin \phi _1^\mathrm{0PN}-\zeta \cos \phi _1^\mathrm{0PN} \sin \chi , \end{eqnarray}$$

(59)

where ζ comes from equation (48).

4.2 The 1.5 PN order: inclusion of the spin precession

Up to the second PN order, the binary dynamics is conservative, the total energy E and the total momentum J = S₁ + S₂ + L_N being the constants of motion. If the S₁ and S₂ spins (m₁ > m₂) are not parallel with the Newtonian orbital momentum L_N, the spins begin to precess (Barker & O’Connell 1975, 1979)

$$\begin{eqnarray} \mathbf {\dot{S}}_i=\Omega _i \times \mathbf {S}_i, \end{eqnarray}$$

(60)

where Ω_i is the angular momentum of the i spin and it contains spin-orbit (1.5PN), spin–spin (2PN) and quadrupole–monopole contributions (2PN) up to 2PN order. In higher PN orders magnetic dipole, self-spin and back-radiation effects appear in the binary dynamics. Due to the gravitational radiation, the eccentricity of the orbits disappears (Peters 1964), therefore from here we take only circular orbits into account.

Gergely & Biermann (2009) have shown that the typical mass ratio range of the merging supermassive black holes is ν ∈ [1/30: 1/3], where the spin-flip of the dominant spin happens sometime during the binary evolution. In the one-spin case, when the second spin S₂ is neglected due to ν ∈ [1/30: 1/3], the Newtonian orbital momentum L_N = μr × v (where v is the velocity of the reduced mass μ), the dominant spin S₁ and the total momentum J lay in the same plane, such that

$$\begin{eqnarray} \alpha _1=\arccos (\mathbf {\hat{L}_N} \cdot \mathbf {\hat{J}}), \end{eqnarray}$$

(61)

$$\begin{eqnarray} \beta _1=\arccos (\mathbf {\hat{S}_1} \cdot \mathbf {\hat{J}}), \end{eqnarray}$$

(62)

$$\begin{eqnarray} \kappa _1=\arccos (\mathbf {\hat{S}}_{1}\cdot \mathbf {\hat{L}_N}), \end{eqnarray}$$

(63)

where κ₁ = α₁ + β₁ is the spin angle. Then the angular velocity vector Ω₁ of the dominant spin S₁ with only spin–orbit (SO) contribution in 1.5PN order (Gergely 2010) is:

$$\begin{eqnarray} \mathbf {\Omega _1}^\mathrm{SO}=\frac{G(4+3\nu)}{2c^2 r^3}L_\mathrm{N}\mathbf {\hat{L}_N}, \end{eqnarray}$$

(64)

where r = |r| is the separation along the unit vector |$\mathbf {\hat{r}}=\mathbf {r}/|\mathbf {r}|$| and e = 0, from which

$$\begin{eqnarray} \Omega _1^\mathrm{SO}=\frac{c^{3}}{2Gm}\varepsilon ^{5/2}\eta \left(4+3\nu \right). \end{eqnarray}$$

(65)

Gergely & Biermann (2009) have shown that in the inspiral phase of the merger, when L_N > S₁ holds, the radiation time-scale of the orbital momentum due to the gravitational waves and the angular velocity of the spin-orbit precession are:

$$\begin{eqnarray} \frac{1}{T_\mathrm{GW}} =-\frac{\dot{L}}{L}=\frac{32c^{3}}{5Gm}\varepsilon ^{4}\eta , \end{eqnarray}$$

(66)

$$\begin{eqnarray} \Omega _\mathrm{p} =\frac{2c^{3}}{Gm}\varepsilon ^{5/2}\eta , \end{eqnarray}$$

(67)

respectively, where η = ν(1 + ν)⁻² ∈ [0: 0.25] is the asymmetrical mass ratio.

In Section 4.1, we derived the inclination under which the jet ejects the blobs |$\iota _\mathrm{j}^\mathrm{0PN}$| (⁠|$\dot{\iota }_\mathrm{j}^\mathrm{0PN}=0$|⁠), such that a periodicity in this direction arose due to the orbital motion. This angle, described by equation (57), is valid only in the 0PN order while the dominant spin does not move and |$\dot{\phi }_1^{0PN}=0$|⁠. This changes in 1.5 PN order because the spin precesses due to the spin-orbit precession and |$\dot{\phi _1}^\mathrm{1.5PN}\ne 0$|⁠. Considering that the angle between the spin and the total angular momentum |$\beta _1= \arccos (\mathbf {\hat{S}_1} \cdot \mathbf {\hat{J}})$| does not change in the conservative dynamics (⁠|$\dot{\beta }_1=0$|⁠) and that the jet emission alone (Blandford & Znajek 1977) can carry away the spin only typically over the Hubble time (Armitage & Natarajan 1999; Mangalam, Gopal-Krishna & Wiita 2009), we can approximate the inclination angle of the dominant spin such that (similarly to equation 59)

$$\begin{eqnarray} \phi _1^\mathrm{1.5PN}(t)=j_0+\beta _1 \cos j_0\sin \Omega _\mathrm{p} t, \end{eqnarray}$$

(68)

where j₀ is the inclination of J. Inserting equation (68) to equation (59), the inclination angle of the jet symmetry axis (or the direction of the jet ejection in case of an intrinsically straight jet ridge line) becomes

$$\begin{eqnarray} \sin \iota _\mathrm{j}^\mathrm{1.5PN}(t)&=&\sin \phi _1^\mathrm{1.5PN}(t)-\zeta (t)\cos \phi _1^\mathrm{1.5PN}(t) \sin (\chi -\psi _1), \\ \iota _\mathrm{j}^\mathrm{1.5PN}(t) &\approx& j_0+\beta _1 \cos j_0\sin \Omega _\mathrm{p} t - \zeta (t) \sin (\chi -\psi _1), \end{eqnarray}$$

(69)

in radians, where we used the first-order paraxial approximation of ϕ_1.5PN(t) since ι₀ < 10° all the time, see Fig. 2. For circular orbits, the true anomaly equals the mean anomaly χ = nt, where n = 2π/T with T being the orbital period. Therefore for circular orbit equation (48) simplifies to

$$\begin{eqnarray} \cos \zeta (t)=\frac{1-X f_\mathrm{c}}{(X^2-2X f_\mathrm{c})^{1/2}}, \end{eqnarray}$$

(70)

where f_c = sin κ₁sin (χ − ψ₁) (setting e = 0 in equation 50) and

$$\begin{eqnarray} v_{1}=\sqrt{\frac{Gm}{r}}\left(\frac{\nu }{1+\nu }\right) \end{eqnarray}$$

(71)

is the circular velocity of the jet emitter (dominant) black hole. Then, after correcting for the cosmological redshift, equation (69) can be fitted to the inclination angle of the symmetry axis of the (intrinsically) helical jet.

5 TEST AND RESULTS

5.1 When the dominant spin is perpendicular to the orbital plane (κ₁ = 0)

In Section 3, we derived the observed period as T′ = (5.98 ± 0.10) yr, from which, considering the redshift of S5 1928+738 as z_spect = 0.302, we get the orbital period in the source frame as T ≈ 4.6 yr. Employing the total mass m = 8.13 × 10⁸ M_⊙ as in Kun et al. (2014), we get the orbital separation in the source frame as r ≈ 0.013 pc and then the PN parameter becomes ε ≈ 0.003. From jet kinematics, we identify v_s = 0.992. We define

$$\begin{eqnarray} \frac{v_1}{v_\mathrm{s}}=\frac{\varepsilon ^{1/2}}{\beta _\mathrm{s}}\left(\frac{\nu }{1+\nu }\right)=K\left(\frac{\nu }{1+\nu }\right)=0.055 \left(\frac{\nu }{1+\nu }\right), \end{eqnarray}$$

(72)

and then equation (69) becomes

$$\begin{eqnarray} \iota _\mathrm{j}^\mathrm{1.5PN}(t)&=& j_0+\beta _1 \sin \Omega _\mathrm{p} t-(1+z_\mathrm{spect}) \\ &&\times \, \arccos \left[\frac{1-K \left(\frac{\nu }{1+\nu }\right) \sin \kappa \sin (n^{\prime }t-\psi _1)}{\left(K^2 \left(\frac{\nu }{1+\nu }\right)^2+1- 2 K \left(\frac{\nu }{1+\nu }\right)\sin \kappa \sin (n^{\prime }t-\psi _1)\right)^{1/2}}\right] \\ && \times \, \sin (n^{\prime }t-\psi _1). \end{eqnarray}$$

(73)

From our fit in Section 3 (see Table 5), we identified the constant term as A₀ = j₀ = (8.64 ± 0.16)°, the linear trend as A₂ = β₁cos j₀sin Ω_p = (− 0.044 ± 0.01)° yr⁻¹, the observed period as T′ = (5.98 ± 0.10) yr, and the spin angle as ϕ = ψ₁ = (24.9 ± 16.1)°. Then we fitted equation (73) to the variation of the inclination angle of the symmetry axis of the 15-GHz VLBI jet of S5 1928+738 (see in Fig. 4.), keeping κ₁ zero (i.e. the dominant spin is perpendicular to the orbital plane, S₁||L_N) and K = 0.055 (see equation 72). Then the best-fitting mass ratio emerged as ν(κ = 0) = 0.209 ± 0.041.

Figure 4.

Upper panel: Inclination variation of the symmetry axis of the 15-GHz jet of S5 1928+738 (dots with error bars), and the dynamical fit with κ ₁ = 0° (continuous line) and with κ ₁ = 75° (dashed line). Lower panel: Observed total flux density of the 15-GHz VLBI jet of S5 1928+738 (dots with errors, note the error bars are usually smaller than the symbols), overlaid with the prediction based on the dynamical fit of the inclination variation of the symmetry axis of the inner jet (continuous line).

5.2 When the dominant spin direction is arbitrary (0 ≤ κ₁ ≤ π): constrains on the mass ratio

We repeated the test described in Section 5.1, but now fitting the mass ratio for fixed values of κ₁ between 0 and π. We show the resulted ν (κ₁) values and their ±1σ errors in Fig. 5. The reduced χ² of the fits changes between |$\chi ^2_\mathrm{red}=0.53$| and |$\chi ^2_\mathrm{red}=0.55$|⁠, indicating that we cannot choose a best model solely based on the χ² values. We can constrain the mass ratio from ν = 0.21 ± 0.04 (κ₁ = 0) to ν = 0.32 ± 0.07 (κ₁ = π/2).

Figure 5.

The best-fitting mass ratio ν and its ±1σ error for a given spin angle κ₁.

5.3 Sanity check with the total flux density

We show the total flux density of S5 1928+738 observed at 15 GHz in the lower panel of Fig. 4. This radiation comes dominantly from the inner jet. The apparent flux density F′ depends on the Doppler factor δ = γ⁻¹(1 − βcos ι)⁻¹, such that

$$\begin{eqnarray} F^{\prime }(\nu)=F (\nu)\delta ^{n-\alpha } , \end{eqnarray}$$

(74)

where F ∼ ν^α is the intrinsic flux density, α is the spectral index and the value of n in the exponent depends on the shape of the jet structure (n = 2 for a linear jet and n = 3 for a spherical outflow, see in, e.g. Urry & Padovani 1995). Since S5 1928+738 has a long jet, in our model n = 2. We estimated the apparent flux density variation of the inner jet of S5 1928+738 by substituting the fitted inclination angle variation of the jet symmetry axis ι(t) into the Doppler factor as δ(t) = γ⁻¹(1 − βcos ι(t))⁻¹. We already derived the characteristic jet velocity as β ≈ 0.992 (from which γ = (1 − β²)^−1/2 ≈ 7.92) and a flat spectrum (α ≈ −0.054 for S5 1928+738, see in Kun et al. 2014). Since we have no information on the intrinsic flux-density, this is the only free parameter of the fit. We show the resulted F′(t) curve in the lower panel of Fig. 4, where the best-fitting flux density in the source frame is F = 0.083 ± 0.002 Jy. Though the fit is only marginal in terms of the statistical goodness (χ²_red ≈ 219), we still got back the main features of the observed flux density variations solely by assuming the best-fitting parameters from our inclination fit ι(t) and freeing the intrinsic flux density.

There might be several reasons to explain why the prediction of the flux density variation based on the binary model is somewhat poor. Once, our jet model does not contain synchrotron emissivity structure (e.g. a non-uniformly distributed electron population, for example a spine–sheath structure), it is a kinematic model of the time variation of the symmetry axis of the inner jet. In our model, the jet components are ejected either under smaller or larger viewing angles compared to the symmetry axis because of the inner helicity. Therefore deviations from the average behaviour, ups and downs on the average flux density curve, are expected.

6 DISCUSSION AND CONCLUSIONS

In Kun et al. (2014), we argued that the periodic and monotonic terms in the change of the inclination of the 15 GHz inner jet symmetry axis of S5 1928+738 can be best explained by spin-orbit precession. Jet precession arising from other phenomena, e.g. Lense–Thirring (LT) precession, magnetic instabilities, water-hose models, etc. could explain only either the periodic, or the monotonic term, but not both simultaneously.

The jet of a rotating black hole can precess due to the LT precession of a massive accretion disc that leads to the Bardeen–Petterson effect (Bardeen & Petterson 1975). In this case, the viscosity of a massive accretion disc of a Kerr black hole, combined with the LT effect, forces the innermost region of the accretion disc to be parallel to the rotation axis of the black hole (see e.g. Liska et al. 2019), while the outermost regions of the accretion disc keep the original direction of the momenta. This combined effect leads to the precession of the accretion disc and to the precession of the jet (Blandford & Znajek 1977; Liska et al. 2018). In this case, the direction of the symmetry axis of the jet helix changes either periodically, or in a monotonic way, depending on the time-scale of the precession (see an application of the LT model in, e.g. fig. 3 of Caproni, Abraham & Mosquera Cuesta 2006). This is contrary to our results, because what wee see is a periodic and monotonic change, acting together on the jet structure.

Magnetohydrodynamic (MHD) instabilities, most prominently Kelvin–Helmholtz (KH) instabilities, can also lead to wiggled jets (for application, see e.g. Perucho et al. 2006). In this case, what we see is not a pattern motion, but jet components moving along curved paths. This might be the cause of the internal helicity of the jet in S5 1928+738. What we see in the case of S5 1928+738 is that the direction of the symmetry axis of this jet helix changes periodically, on top of a much smaller, seemingly monotonic trend. Models relying on MHD (or specifically KH) instabilities, can not explain the observed complexity in S5 1928+738: A periodically changing inclination angle of the symmetry axis on top of a monotonic decrease.

Assuming m = 8.13 × 10⁸ M_⊙ as the total mass of the binary (ε ≈ 0.003, r ≈ 0.0129), we derived the orbital period in the source frame as 4.78 yr and the lower limit on the mass ratio as ν_lo = 0.21 from the jet kinematics and the periodic term of the inclination variation of the symmetry axis of the inner jet of S5 1938+738. For ν > 0.21, the spin-orbit period is T_SO < 5500 yr and the merger time is T_GW < 1.64 × 10⁶ yr. Assuming the upper limit of typical mass ratio range of merging SMBHs ν < 1/3, we constrained the spin-orbit period with lower limits as T_SO > 4207 yr and the merger time as T_GW > 1.25 × 10⁶ yr.

We summarize our results regarding the 1.5 PN dynamics of the hypothetical binary of S5 1928+738 in Table 7. In the recent analysis, we derived the orbital period as ≈4.6 yr. Then, assuming the total mass as m = 8.13 × 10⁸ M_⊙, the orbital separation (for circular orbit) arises as r ≈ 0.0125 pc. From the jet kinematics, we got ν = 0.21 ± 0.04 for κ₁ = 0 (i.e. when the spin direction is perpendicular to the orbital plane) and ν = 0.32 ± 0.07 for κ₁ = π/2 (i.e. when the spin lies in the orbital plane). For the mass ratio ν = 0.21 ± 0.04, the spin-orbit precession is |$T_\mathrm{SO}=5171^{+802}_{-536}$| yr and the merger time is |$T_\mathrm{GW}=(1.48^{+0.23}_{-0.15})\times 10^6$| yr. For the mass ratio ν = 0.32 ± 0.07, the spin-orbit precession is |$T_\mathrm{SO}=4039^{+597}_{-364}$| yr and the merger time is |$T_\mathrm{GW}=(1.16^{+0.17}_{-0.10})\times 10^6$| yr.

Table 7.

Estimation of the mass ratio, the spin-orbit precession period and the merger time for S5 1928+738, assuming the total mass as m = 8.13 × 10⁸ M_⊙ and the orbital period as T ≈ 4.6 yr in the source frame.

	ν	TSO	TGW
κ1 = 0	0.21 ± 0.04	|$5171^{+802}_{-536}$| yr	|$(1.48^{+0.23}_{-0.15})\times 10^6$| yr
κ1 = π/2	0.32 ± 0.07	|$4039^{+597}_{-364}$| yr	|$(1.16^{+0.17}_{-0.10})\times 10^6$| yr

The dynamical evolution of an SMBH binary rotates the spin to be parallel with the Newtonian orbital momentum L_N, which means the spin is perpendicular to the orbital plane. For κ₁ = 0°, it is an aligned, for κ₁ = π, it is an anti-aligned configuration. Once κ₁ becomes 0 or π, the configuration itself becomes stable. In this case, the spin-orbit precession disappears. Since we see the slow, secular reorientation of the jet symmetry axis (on top of a much faster periodicity, due to the orbital motion), our model is consistent with the observations and phenomenological model when we exclude the possibility of κ₁ = 0 and κ₁ = π.

Based on the periodic modulation and the monotonic trend in the inclination angle of the symmetry axis of the inner jet, we further supported the existence of the SMBH binary in the heart of the quasar S5 1928+738. We made a more precise estimation of the velocity of the spin precession, from (−0.05 ± 0.02) to (−0.04 ± 0.01)° yr⁻¹. Using this rate, we estimate the symmetry axis of the jet goes through our line of sight sometime between 2180 and 2200, which is far from the present. Finally, we conclude that the precession of the 15-GHz VLBI inner jet of S5 1928+738 is still ongoing, showing a similar behaviour over 25 yr between 1995.96 and 2020.89. Since the jet is approaching our line of sight and the Doppler boosting strengthens non-linearly with the decreasing inclination angle, we advise to monitor the quasar S5 1928+738 across all regimes of the electromagnetic spectrum. Depending on the geometry, either the jet will go through our line of sight, or turn back toward larger inclination angles. The merging SMBH system in S5 1928+738 is in the beginning of its inspiral phase, emitting gravitational waves at a frequency of 10⁻⁸ Hz and with a characteristic strain of the order of 10⁻¹⁴ (based on the calculations in Sesana 2016). Further multimessenger observations might prove further or challenge the existence of the SMBH binary harboured in S5 1928+738.

ACKNOWLEDGEMENT

We thank Zoltán Haiman for valuable comments on our earlier related work. EK thanks the Alexander von Humboldt Foundation for its Fellowship. We acknowledge support from the Deutsche Forschungsgemeinschaft DFG, within the Collaborative Research Center SFB1491 ‘Cosmic Interacting Matters – From Source to Signal’ (project No. 445052434). Support by the Hungarian National Research, Development and Innovation Office (NKFIH) is acknowledged (grant numbers OTKA K134213 and K123996). This research has made use of data from the MOJAVE database that is maintained by the MOJAVE team (Lister et al. 2018). The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

DATA AVAILABILITY

The data underlying this article are available at https://www.cv.nrao.edu/MOJAVE/sourcepages/1928 + 738.shtml as part of the MOJAVE Survey.