A&A 374, 784-799 (2001)
DOI: 10.1051/0004-6361:20010685

On the origin of compact radio sources

The binary black hole model applied to the gamma-bright quasar PKS 0420-014

S. Britzen1 - J. Roland2 - J. Laskar3 - K. Kokkotas4 - R. M. Campbell5 - A. Witzel6


1 - Netherlands Foundation for Research in Astronomy, Oude Hoogeveensedijk 4, 7991 PD Dwingeloo,
The Netherlands
2 - Institut d'Astrophysique, 98bis boulevard Arago, 75014 Paris, France
3 - Bureau des Longitudes, 77 avenue Denfert Rochereau, 75014 Paris, France
4 - Department of Physics, Aristotle University of Thessaloniki, Thessaloniki 54006, Greece
5 - Joint Institute for VLBI in Europe, Oude Hoogeveensedijk 4, 7991 PD Dwingeloo, The Netherlands
6 - MPIfR, Auf dem Hügel 69, 53121 Bonn, Germany

Received 11 July 2000 / Accepted 10 May 2001

Abstract
VLBI monitoring of PKS 0420-014 at 3.6 cm performed during the period 1989.32-1992.48 reveals superluminal jet components in a strongly curved jet. We present a binary black hole model capable of simultaneously explaining the flux-density variations and the motions of jet components on helical trajectories in active galactic nuclei. From VLBI observations we derive the characteristics of the parsec-scale jet of PKS 0420-014. Moreover, we use the optical light curve of the quasar PKS 0420-014 over the time span 1989.67-1993.93 to determine the characteristics of the binary black hole system. We propose that galactic and extragalactic radio sources with jets are due to the existence of binary systems. Extragalactic binary systems are the results of mergers and provide a natural explanation why only a small fraction of quasars are radio sources and why extragalactic radio sources are associated with elliptical galaxies. These binary systems will be one of the primary targets for the Laser Interferometer Space Antenna (LISA).

Key words: black holes physics - gravitational waves - quasars: individual: PKS 04020-014 - radio continuum: galaxies


1 Introduction

Extragalactic compact radio sources associated with nuclei of galaxies often show flux-density variations across the wavelength spectrum on different timescales. VLBI-component ejection in these sources is often associated with an optical and/or a $\gamma$-ray outburst. Typical sources of this kind are 3C 273 and 3C 345 (see Zensus 1997 and references therein). Moreover, VLBI observations show that jet components of several sources move on helical trajectories. The radio source 3C 273 shows cycles of quasi-periodic activity and ejections of VLBI components on helical paths close to the nucleus. It has been shown by Roland et al. (1994) that the variability of 3C 273 can be explained by a precessing accretion disk. In the case of the BL Lac Object OJ 287, the optical variability recorded over a time span of one hundred years can be successfully explained by a binary black hole (BBH) model, e.g., Sillanpää et al. (1988) and Lehto & Valtonen (1996). According to Lehto & Valtonen (1996) the smaller black hole of a binary pair crosses the accretion disk of the larger black hole during the orbit, and causes the substructure inside the major outbursts. Vicente et al. (1996) find evidence from 3.6 cm geodetic VLBI monitoring of OJ 287 for a helical morphology of the jet and for component ejection during optical outbursts. Their findings are consistent with the supermassive BBH model proposed to explain the periodic optical outbursts. Other sources where BBH models in the nuclei have been discussed are Mrk 501 (Conway & Wrobel 1995) and 1928+738 (Roos et al. 1993). PKS 0420-014 is known for its pronounced flux-density variability across several wavelength bands. Wagner et al. (1995) predict jet component ejection and component motion on a curved path within the framework of the lighthouse model (Camenzind & Krockenberger 1992) based on the observed optical flaring. Recently, Britzen et al. (2000a) published the results of nine epochs of VLBI monitoring of PKS 0420-014 taken over a period of three years (1989.32-1992.48). Wagner et al. (1995) used four of these epochs to demonstrate for the first time superluminal motion on a curved path in this source. Britzen et al. (2000a) studied the structural evolution of PKS 0420-014 and found superluminal motion with $\beta_{{\rm app\/}}$ ranging from 2 to 14 c for five jet components. In addition they discuss evidence for a correlation between component ejection and flux-density outbursts.

Helical trajectories and the flaring behavior seen in the optical light curve (1989-1993) of PKS 0420-014 have been explained in the framework of the lighthouse model (Camenzind & Krockenberger 1992; Wagner et al. 1995). However, in this scenario the mass of the central object is found to be very large: i.e., for an orbital radius of the knot of 0.2 light years it is $M\simeq1.05\times10^{10}~M_{\odot}$ and for an orbital radius of 0.75 light years it is $M\simeq3.85\times10^{10}~M_{\odot}$ (Wagner et al. 1995). The observational data base necessary to test and restrict theoretical models for PKS 0420-014 has recently been significantly improved due to a more detailed analysis of the VLBI epochs (Britzen et al. 2000a). This enables a new attempt in modelling this source. To explain the origin of the precession of the accretion disk and the origin of the optical light curve, we propose that PKS 0420-014 contains a binary system of massive black holes. As we will see, the BBH model requires significantly fewer massive black holes compared to the lighthouse model.

We first develop the mathematical formalism of the binary system model and the precession of the accretion disk. We then briefly discuss the overall flux-density variability observed for PKS 0420-014 and the kinematic information deduced from the VLBI observations of the jet components. We apply the BBH model to show that the kinematics of the jet components can be explained by the precession of the accretion disk and that the light curve of the optical outbursts can be explained by the motion of the jet-emitting black hole around the center of gravity of the system. Using VLBI and optical observations, we deduce the characteristics of the binary system. Finally, we discuss the hypothesis that binary systems underlie the existence of compact radio sources, and we investigate the possibility of detecting low frequency gravitational waves from such systems using the LISA experiment.

2 The binary system model

2.1 The two-fluid model

  In this section, we develop the mathematical formalism to describe a system consisting of two black holes and the influence of this system on the jet and its emission. We assume that the nuclei of extragalactic radio sources eject two fluids (the two-fluid model: Sol et al. 1989; Pelletier & Roland 1989, 1990; Pelletier & Sol 1992). These are:
(i) an e $^{-}-{\rm p}$ plasma called the jet. It is mildly relativistic, ejected by the nucleus at speeds $v_{j}\leq0.5\;c$;
(ii) an e$^{\pm}$ plasma called the beam. It is relativistically ejected by the nucleus at speeds $v_{\rm b}\simeq
c,$ i.e., $3\leq\gamma_{\rm b}\leq20,$ where $\gamma_{\rm b}$ is the bulk Lorentz factor.

The e $^{-}-{\rm p}$ jet carries most of the mass and the kinetic energy ejected by the nucleus, and is responsible for the formation of kpc-jets, hot spots, and extended lobes. The relativistic e$^{\pm}$ beam moves in a channel through the mildly relativistic jet and is responsible for the formation of superluminal sources and their $\gamma$-ray emission (Roland et al. 1994). The relativistic beam can propagate if the magnetic field B is parallel to the flow in the beam and in the mixing layer between the beam and the jet and if it is greater than a critical value (Pelletier et al. 1988; Achatz & Schlickeiser 1993). The magnetic field in the e $^{-}-{\rm p}$ jet rapidly becomes toroidal (Pelletier & Roland 1990). The observational evidence for the two-fluid model has been discussed by Roland & Hetem (1996) for instance. Recent observational evidence for relativistic ejection of e$^{\pm}$ beams comes from the $\gamma$-ray observations of MeV sources (Roland & Hermsen 1995; Skibo et al. 1997) and from VLBI polarization observations (Attridge et al. 1999). The formation of the X-ray and $\gamma$-ray spectra, assuming relativistic ejections of e$^{\pm}$ beams, has been investigated by Marcowith et al. (1995, 1998) in the case of Centaurus A. The possible existence of VLBI components with two different speeds has been recently pointed out in the case of the radio galaxies Centaurus A (Tingay et al. 1998) and M 87 (Biretta et al. 1999). In addition to observations of one-sided superluminal sources, the existence of two-sided subluminal VLBI sources has been well established in the case of Cygnus A (Krichbaum et al. 1998), 3C 338 (Feretti et al. 1993), and Centaurus A (Tingay et al. 1998).

2.2 The precession of the accretion disk

We assume that the relativistic beam moves along the magnetic field lines, which are perturbed by the precession of the accretion disk (Roland et al. 1994). We suppose that the accretion disk makes an angle $\Omega_{\circ}$ with the plane (xOy), such that the perturbed beam is on the surface of a cone with symmetry axis $(z^{\prime}Oz)$ and opening angle $\Omega_{\circ}$. Moreover, we suppose that the ejected component moves with a constant speed characterized by its bulk Lorentz factor $\gamma_{\rm c}$ and follows the magnetic field lines. The perturbation of the magnetic field propagates with a speed given by the Alfvén speed $V_{\rm A}$. The coordinates of the component moving in the perturbed beam are given by

   
$\displaystyle %
x_{\rm c} =R(z_{\rm c}(t))\cos(\omega_{\rm p}t-k_{\rm p}z_{\rm c}(t))$     (1)
$\displaystyle y_{\rm c} =R(z_{\rm c}(t))\sin(\omega_{\rm p}t-k_{\rm p}z_{\rm c}(t))$     (2)
$\displaystyle z_{\rm c} =z_{\rm c}(t)$     (3)

where $\omega_{\rm p}=2\pi/T_{\rm p}$, $T_{\rm p}$ is the precession period, and $k_{\rm p}$ is defined by
$\displaystyle %
k_{\rm p}=\frac{2\pi}{T_{\rm p}V_{\rm A}}\cdot$     (4)

We take the form of the amplitude $R(z_{\rm c}(t))$ to be

\begin{eqnarray*}R(z_{\rm c}(t))=\frac{R_{0}z_{\rm c}(t)}{a+z_{\rm c}(t)},
\end{eqnarray*}


where a is $a=R_{0}/(2\tan\Omega_{\circ}).$ The differential equation governing the evolution of $z_{\rm c}(t)$ can be obtained via the relation for the speed of the component
 
$\displaystyle %
v_{\rm c}^{2}=\left( \frac{{\rm d}x_{\rm c}(t)}{{\rm d}t}\right...
...{{\rm d}t}\right) ^{2}+\left( \frac{{\rm d}z_{\rm c}(t)}{{\rm d}t}\right) ^{2},$     (5)

where $v_{\rm c}$ is related to the bulk Lorentz factor by $v_{\rm c}/c=\sqrt
{(1-1/\gamma_{\rm c}^{2}).}$ Using (1) and (2), we find from (5) that ${\rm d}z_{\rm c}(t)/{\rm d}t$ is the solution of the equation
 
$\displaystyle %
\left(1+R^{2}(z_{\rm c})k_{\rm p}^{2}+({\rm d}R/{\rm d}z_{\rm c...
...\rm d}z_{\rm c}}{{\rm d}t}
+R^{2}(z_{\rm c})\omega_{\rm p}^{2}-v_{\rm c}^{2}=0,$     (6)

which admits two solutions corresponding to a jet and a counter-jet. We assume that the line of sight is in the plane (xOz) and makes an angle i with the z axis. Following Camenzind & Krockenberger (1992), if we call the angle between the velocity of the component and the line of sight $\theta$, we have

\begin{eqnarray*}\cos(\theta(t))=\left( \frac{{\rm d}x_{\rm c}}{{\rm d}t}\sin i+\frac{{\rm d}z_{\rm c}}{{\rm d}t}\cos
i\right) /v_{\rm c}.
\end{eqnarray*}


The Doppler beaming factor $\delta$, characterizing the anis- otropic emission of the moving component, is

\begin{eqnarray*}\delta(t)=\frac{1}{\gamma_{\rm c}[1-\beta_{\rm c}\cos(\theta(t))]}\;,
\end{eqnarray*}


where $\beta_{\rm c}=v_{\rm c}/c.$ The observed flux density is

\begin{eqnarray*}S_{\rm c}=\frac{1}{D^{2}}\delta(t)^{2+\alpha_{\rm r}}(1+z)^{1-\alpha_{\rm r}}\int_{c}%
j_{\rm c}{\rm d}V,
\end{eqnarray*}


where D is the distance to the source, z is its redshift, $j_{\rm c}$ is the emissivity of the component, and $\alpha_{\rm r}$ is the synchrotron spectral index. As the component is moving relativistically toward the observer, the observed time is shortened and is given by

\begin{eqnarray*}t_{{\rm obs\/}}=\int_{0}^{t}[1-\beta_{\rm c}\cos(\theta(t^{\prime
}))](1+z){\rm d}t^{\prime}.
\end{eqnarray*}


If we consider $i\simeq0$, the observed quasi-period of a quasi-periodic flaring source is shortened approximately in the following way

\begin{eqnarray*}P_{{\rm obs\/}}\simeq\frac{P(t=0)}{2\gamma_{\rm c}^{2}}(1+z).
\end{eqnarray*}


2.3 The binary system

In order to explain the origin of the precession, we propose that the nucleus is part of a binary system (see Appendix I). In the following, index 1 denotes the black hole ejecting the VLBI components, index d denotes the accretion disk and index 2 denotes the second black hole (the correspondence with indices in Appendix I is 0, 1, 2 respectively). The two black holes orbit in the plane (xOy), and the origin is centered at the mass center of the system. The elliptical orbit is given by

\begin{eqnarray*}r=\frac{p}{1+e\cos(\varphi)},
\end{eqnarray*}


where e and p are respectively the eccentricity and the parameter of the orbit. The coordinates of black hole 1 are given by
x1(t)=$\displaystyle \frac{m_{2}}{m_{1}+m_{2}}\frac{p}{1+e\cos(\varphi(t))}%
\cos(\varphi(t))\;$ (7)
y1(t)=$\displaystyle \frac{m_{2}}{m_{1}+m_{2}}\frac{p}{1+e\cos(\varphi(t))}%
\sin(\varphi(t))\;.$ (8)

Then the system of Eqs. (1)-(3) becomes
   
$\displaystyle %
x_{\rm c}$=$\displaystyle x_{1}(t)\cos(\omega_{\rm b}t{-}k_{\rm b}z_{\rm c}(t)){+}R(z_{\rm c}(t))\cos(\omega
_{\rm p}t-k_{\rm p}z_{\rm c}(t))\;$  
      (9)
$\displaystyle y_{\rm c}$=$\displaystyle y_{1}(t)\sin(\omega_{\rm b}t{-}k_{\rm b}z_{\rm c}(t)){+}R(z_{\rm c}(t))\sin(\omega
_{\rm p}t-k_{\rm p}z_{\rm c}(t))\;$  
      (10)
$\displaystyle z_{\rm c}$=$\displaystyle z_{\rm c}(t),$     (11)

where $\omega_{\rm b}=2\pi/T_{\rm b}$, $T_{\rm b}$ is the binary system period, and $k_{\rm b}$ is defined by

\begin{eqnarray*}k_{\rm b}=\frac{2\pi}{T_{\rm b}V_{\rm A}}\cdot
\end{eqnarray*}


As in Eq. (6) ${\rm d}z_{\rm c}(t)/{\rm d}t$ is the solution of an equation of the following form

\begin{eqnarray*}A\left( \frac{{\rm d}z_{\rm c}}{{\rm d}t}\right) ^{2}+B\left( \frac{{\rm d}z_{\rm c}}{{\rm d}t}\right) +C=0,
\end{eqnarray*}


where the coefficients A, B, and C are derived in Appendix II.

3 PKS 0420-014

3.1 The variability of PKS 0420-014

Two sorts of variability have been detected in the source PKS 0420-014. The source is known for its pronounced flux-density variability in several parts of the electromagnetic spectrum (Wagner et al. 1995; Britzen et al. 2000a and references therein). PKS 0420-014 did undergo a simultaneous optical/gamma-ray flaring at the time of a radio outburst, providing an opportunity to investigate the relation between the structural evolution and the flaring activity seen in the radio, optical, and gamma-ray regimes (Wagner et al. 1995; Britzen et al. 2000a).

A detailed description and analysis of five years of optical monitoring within the Hamburg Quasar Monitoring Program and at the Landessternwarte in Heidelberg is presented in Wagner et al. (1995). PKS 0420-014 reveals two modes of variability at optical wavebands (see Fig. 1, this paper). On top of an almost regular major flaring cycle having a period of $\sim$13 months (peaks in early 1990, 1991, and 1992) in which each flare lasts about one month, variations on shorter timescales are detected. The EGRET instrument detected gamma-radiation above 100 MeV during two observation periods in early 1992 (Hartman et al. 1999). The epoch of the strong gamma-ray flaring in 1992 coincided with the strongest outburst of optical emission observed in 5 yrs of monitoring since 1989 (Wagner et al. 1995).

VLBI monitoring performed at 3.6 cm revealed structural variability in the jet of PKS 0420-014 (Britzen et al. 2000a) possibly correlated with brightness variations in several bands. The BBH model introduced in the previous section is capable of explaining the radio morphology and its changes, as well as the flux-density flaring observed in the optical waveband. In the following section, we apply this model to the quasar PKS 0420-014. A summary of the flux-density variability and a detailed description of the VLBI monitoring can be found in Britzen et al. (2000a); first results of the BBH modelling have already been presented in Britzen et al. (1998b).

3.2 VLBI morphology

Britzen et al. (2000a) present nine VLBI observations of the OVV quasar PKS 0420-014 performed at $\lambda=3.6$ cm between 1989.32 and 1992.48 and one VLBI observation at 43 GHz. A subsample of four of the nine VLBI observations have already been presented in Wagner et al. (1995). In these four epochs the motion of the two brightest jet components with respect to the core was traced and for the first time superluminal motion as well as motion on a curved path could be proven for this source. The large number of VLBI observations available now enables a more detailed analysis of the jet structure out to core separations of $\sim$3 mas, and out to $\sim$6 mas with higher uncertainties in parameter estimates. The jet can be well described by up to eight distinct jet components. In this analysis, we concentrate on the five innermost components (A-E), labelled according to their increasing separation from the core, see Figs. 3-5. A linear regression analysis yields superluminal motion of value $13.8 \pm 8.7~c$ for component A, $7.4 \pm 0.5~c$ (B), $5.0 \pm 0.7~c$ (C), $3.8 \pm 0.7~c$ (D), and $2.1 \pm 0.6~c$ (E). The components seem to separate faster in the vicinity of the core; the "older'' jet components, at larger core separations, move more slowly.

In this article, we assume:

Another possibility, not investigated in this article, is to assume that all the VLBI components move on the same path and that their bulk Lorentz factors are the same at the beginning of the ejection but they decrease as a function of time; the outer components are then observed with lower apparent speeds (let us note that this change does not affect the results found for the BBH system).

The global VLBI observations performed at 43 GHz indicate a similar path for the jet components, but one in which the bending continues towards the innermost region and becomes even more pronounced with a ${\sim}90^{\circ}$ turn within 0.5 mas of the core. The total bending within the inner 2 mas must be at least ${\sim}160^{\circ}$. Britzen et al. (2000a) estimate epochs of component ejection from the core using the derived component velocities and trajectories. They discuss possible correlations between epochs of component ejection and multiwavelength flux-density activity. Although the sampling of the flux-density data and the large error bars of the calculated component ejection epochs do not allow a definite correlation of such events, the onset of ejection of the relativistic particles responsible for the radio emission of components A and B appears to occur before the ejection of the relativistic particles responsible for their simultaneous optical/gamma-ray flaring. The observational data shown in this paper result from elliptical Gaussian modelfitting and have been discussed in detail in Britzen et al. (2000a). The uncertainties for each component assumed there are: for the flux S, $\Delta S$ ${\sim} 10{-}20\%$; for the distance r from the core, $\Delta r \sim 0.1$ mas.

4 Binary Black hole model for PKS 0420-014

4.1 General constraints for the radio and the optical fit

We propose to explain the phenomena seen in PKS 0420-014 in the radio and the optical regime using only the BBH model. The precession of the accretion disk and the motion of the black holes around the center of gravity of the BBH system produce two different perturbations of the magnetic field lines which confine the e$^{\pm}$ beam. The perturbations of the magnetic field lines propagate along the beam with a speed given by the Alfvén speed. Within this model, the VLBI observations can be explained by the precession of the accretion disk, and the optical light curve by the motion of the black hole ejecting the radio jets around the center of gravity of the BBH system. The problem that complicates the modelling of both wavelength regimes is the different variability timescales and amplitudes of the observed phenomena. The optical light curve of PKS 0420-014 shows variations with peaks of about 2 magnitudes with a typical time scale of about 1 year. The radio flux-density of the VLBI component B increases over several years with brightness variations somewhat smaller than a factor of 10. As the e$^{\pm}$ responsible for the optical light curve and the e$^{\pm}$ responsible for the radio emission follow the same trajectories, the different behaviors observed in radio and optical ranges imply that the relativistic e$^{\pm}$ responsible for the optical light curve are emitted during a much shorter time than the relativistic e$^{\pm}$ responsible for the radio emission.

Hence, one can explain the variability of PKS 0420-014 as arising from a long burst of relativistic e$^{\pm}$ responsible for the synchrotron radio emission and a shorter burst of very energetic relativistic e$^{\pm}$ responsible for the synchrotron optical emission, where the emission of the e$^{\pm}$ responsible for the synchrotron radio emission precedes the short burst of more energetic e$^{\pm}$ responsible for the optical synchrotron emission. This behavior is typical of compact variable sources (Tornikoski et al. 1994). Thus, the optical emission is modeled as the emission of a point-source component and the radio emission is modeled as the integrated emission of an elongated component along the beam. Let us note that the difference between the optical and the radio light curve cannot be explained by a possible deceleration of the VLBI component. We aim to obtain one set of parameters that allows us to fit the optical and VLBI data simultaneously. As indicated previously we assume the geometry of the system is the same for the different VLBI components and the bulk Lorentz factors and the duration of the emission of the e$^{\pm}$ responsible for the synchrotron radio emission of the different VLBI components are different.

4.2 Characteristic of the precession

In this section we describe how we derived the model parameters that give the best representation of the VLBI observations. In particular, the results of the VLBI monitoring include the flux-density evolution with time (Fig. 2), the linear separation of the components from the core as functions of time (Figs. 3-5), and the apparent speed of the components as a function of time (Fig. 6). We want to determine the geometry of the perturbed magnetic field lines in which the e$^{\pm}$ plasma (the beam) propagates. The geometry of the system is characterized by three parameters: the inclination angle i, the opening angle of the ejection cone $\Omega_{\circ}$, and the radius of the helical trajectory $R_{\circ}$. In addition, the model depends on the precession period $T_{\rm p}$, the propagation speed of the perturbations along the magnetic field lines (i.e., the Alfvén speed) $V_{\rm A}$, and finally the bulk Lorentz factor $\gamma_{\rm c}$ that characterizes the speed of the ejected component. We assume a priori that the geometry is the same for the various components. We then model the results for component B under the assumption that the bulk Lorentz factors of the various components are different. The observed brightness and kinematic evolution of jet component B which constrain our model parameters are:

The bulk Lorentz factor of component B, $\gamma_{\rm B}$, and the duration of the emission of the relativistic e$^{\pm}$ responsible for the synchrotron radio emission of component B have to be chosen to fit simultaneously the radio and optical observations. The values found to be consistent with the observations are, respectively, $\gamma_{\rm B}\simeq8$ derived from the observed motion (Britzen et al. 2000a), and $\tau_{\rm B}\simeq25$ yrs of proper time derived from the fit of the flux density variations and of the core separation. When the characteristics of the binary system are known, it is possible to check the values chosen for the bulk Lorentz factor and the duration of the emission of the relativistic e$^{\pm}$ as discussed in Sect. 4.4. The numerical values obtained from the fit provide a coherent solution; however, the solution presented is not the only possible one. To determine the parameters, we plotted the results of the model with the observed data and choose the values of the parameters which provided the best fit with the data, we did not perform a $\chi^{2}$ analysis and therefore we cannot directly estimate the errors of the deduced parameters. For a given bulk Lorentz factor, there may be 2 possible inclination angles which can explain an observed apparent speed. However, from the shape of the flux variations and the distance of the component after $\tau\simeq3.3$ yrs, it is possible to exclude the large inclination angles. Indeed, it is only for small inclination angles that the flux density increases monotonically during the first 4 yrs. We investigated the following values for the inclination angle: $i=3^{\circ}$, $4^{\circ}$, and $5^{\circ}$ (smaller inclination angles are possible but are not investigated here). Even if the e$^{\pm}$ beam is a pure e$^{\pm}$ plasma, the mixing layer between the e$^{\pm}$ beam and the e--p jet contains protons and the Alfvén speed given by $V_{\rm A}=B/\sqrt{4\pi\rho}$ can be very small compared to the speed of light. We investigated the four possible values for the Alfvén speed: $V_{\rm A}=0.0125\;c$, $0.025\;c$, $0.05\;c$ and $0.1\;c$. From the VLBI constraints for component B, the value of the bulk Lorentz factor, $\gamma_{{\rm B\/}}\simeq8,$ and the value of the duration of the emission of the relativistic e$^{\pm}$ responsible for the synchrotron radio emission, $\tau\simeq25$ yrs of proper time, we obtained the best fit for the following set of 5 parameters i, $\Omega_{\circ}$, $R_{\rm o}$, $V_{\rm A}$ and $T_{\rm p}$ (see Tables 1 and 2). The geometrical parameters corresponding to the 3 inclination angles are given in Table 1. The values $T_{\rm p}$ corresponding to the various $V_{\rm A}$ are given in Table 2; they are valid for the 3 values of the inclination angle.


 

 
Table 1: The geometrical parameters.
  $i=3^{\circ}$ $i=4^{\circ}$ $i=5^{\circ}$
$\Omega_{\circ}(^{\circ})$ 11 12 13
$R_{o}({\rm pc\/})$ 10 10 10



 

 
Table 2: Values of $V_{\rm A}$ and $T_{\rm p}$ corresponding to the geometrical parameters.
$V_{\rm A}/c$ 0.1 0.05 0.025 0.0125
$T_{\rm p}({\rm yr})$ 4400 $10\,000$ $20\,000$ $40\,000$


4.3 Characteristic of the binary system

The beginning of the variability pattern ( ${\sim}1990.0{-}1992.0$) visible in the optical light curve (Wagner et al. 1995) is caused by the birth of component B. As can be seen in Fig. 1, the ejection of component B can not account for all the variability features and we assume that the ejection of component A contributes as well to the optical flaring later. We therefore fit the 1989.67-1993.93 optical light curve as a blend of emissions of components A and B. Wagner et al. (1995) supposed that an optical light curve with half the period indicated in Fig. 1 is possible, in that case the period of the BBH system will be reduced by a factor 2. Assuming that the nucleus of PKS 0420-014 is a binary system, we then have the following free parameters: the masses of the two black holes m1 (mass of the black hole ejecting the radio plasma) and m2 (mass of the second black hole), the eccentricity of the orbit e, and the rotation period $T_{\rm b}$. The main constraints from the optical light curve are the amplitudes of the optical peaks and the time separation between the peaks. In Sect. 4.2 we underlined the necessity to investigate 4 different Alfvén speeds; additionally we investigated the range $10^{6}~M_{\odot}\leq m\leq
10^{9}~M_{\odot}$ for the masses of the 2 black holes and $0\leq e<1$ for the orbital eccentricity. Our calculations were done for the 3 possible inclination angles $i=3^{\circ},\;4^{\circ}$ and $5^{\circ}$ with the bulk Lorentz factors $\gamma_{{\rm B}}=8$ and $\gamma_{{\rm A}}=10$ for the VLBI components B and A respectively. A fit of the optical light curve is possible for the 4 values of the Alfvén speed; however the variations of the VLBI flux density of component B constrain the Alfvén speed reliably to ${\simeq}0.05~c$ value. We find that the parameter combination given in Table 3 allows us to fit the data best. A fit of the optical light curve is possible for values of the eccentricity $0\leq
e\leq0.3$, and in the following we will adopt the value e=0 for the eccentricity, remarking that the characteristics of the binary system do not change significantly for $0\leq
e\leq0.3$.

 

 
Table 3: Characteristics of the binary system (e=0, $V_{\rm A}=0.05~c$ and $T_{\rm p}({\rm yr\/})=10^{4}$ yrs).
i $3^{\circ}$ $4^{\circ}$ $5^{\circ}$
$m_{1}(M_{\odot})$ $5.4\times10^{7}$ $7\times10^{7}$ $9\times10^{7}%
$
$m_{2}(M_{\odot})$ $1.6\times10^{8}$ $2.1\times10^{8}$ $2.7\times10^{8}%
$
$T_{\rm b}({\rm yr})$ 160 150 140


We show the fit of the optical light curve in Fig. 1 obtained for $i=3^{\rm o}$, $m_{1}=5.4\times10^{7}~M_{\odot}$, $m_{2}=1.6\times10^{8}~M_{\odot}$ and $T_{\rm b}=160~{\rm yrs}$.

  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{10101fig1.eps}\end{figure} Figure 1: Fit of the optical light curve of PKS 0420-014. The optical light curve (Wagner et al. 1995) is shown with dots, superimposed are the emission of component A indicated by the dashed line and the emission of component B indicated by the continuous line.
Open with DEXTER

We determine the characteristic sizes of the binary system we found. As the 3 inclination angles are possible, to simplify the discussion and to determine some mean values for the binary system, we will assume in the following that the inclination angle is $i\simeq 4^{\circ }$. As the mass of the accretion disk is much smaller than the masses of the two black holes, the precession period $T_{\rm p}$ of the accretion disk is given by (see Appendix I)[*]

$\displaystyle %
T_{\rm p}=\frac{2\pi}{s}=\frac{4}{3}\frac{m_{1}+m_{2}}{m_{2}}\frac{T_{\rm b}^{2}%
}{T_{\rm d}},$     (12)

where m1 is the mass of the black hole ejecting the material and m2is the mass of the other black hole. We can deduce the mean rotation period $T_{\rm d}$ of the disk, the corresponding mean size $R_{\rm d}$ of the disk, and the distance D12 between the two black holes. The mean rotation period of the disk is

\begin{displaymath}%
T_{\rm d}=\frac{4}{3}\frac{m_{1}+m_{2}}{m_{2}}\frac{T_{\rm b}^{2}}{T_{\rm p}}\cdot
\end{displaymath} (13)

The characteristic sizes of the binary system are given in Table 4. We note that the accretion disk radius is $R_{\rm d}\simeq750\times R_{\rm g}$ where $R_{\rm g}$ is the Schwarzschild radius.

The amplitude of the motion of the black hole ejecting radio jets is about 30 $\mu$as in the case of PKS 0420-014.


 

 
Table 4: Characteristic sizes of the binary system ( $i\simeq 4^{\circ }$, $e\simeq 0$ and $V_{\rm A}=0.05~c$).
$m_{1}(M_{\odot})$ $7\times10^{7}$
$m_{2}(M_{\odot})$ $2.1\times10^{8}$
$T_{\rm p}({\rm yr})$ $10\,000$
$T_{\rm b}({\rm yr})$ 150
$T_{\rm d}({\rm yr})$ 4.0
$R_{\rm d}({\rm pc})$ 0.0050
$D_{12}({\rm pc})$ 0.090


4.4 Characteristics of the different VLBI components

  An important constraint concerning the modelling of the characteristics of VLBI component-emission comes from the difference between the flux-density variations observed in the optical and the radio. Indeed, the optical light curve shows flaring with an amplitude of 2 magnitudes with a typical time-scale of about 1 year and the radio flux density of the VLBI component B increases over years with brightness variations somewhat smaller than a factor of 10. As mentioned previously this important difference can be explained if the durations of the emission of the e$^{\pm}$ responsible for the optical synchrotron emission and for the radio emission are different, namely the relativistic e$^{\pm}$ responsible for the optical light curve are emitted during a much shorter time than the relativistic e$^{\pm}$ responsible for the radio emission. This difference can be understood if we consider that the optical emission is due to a point source component and that the radio emission is due to an extended component whose observed flux is integrated along the beam. The typical time scale of the duration of the emission of the relativistic particles for component B is found to be about 25 yrs of proper time which corresponds to about 1.2 yrs of observed time (see Fig. 7). The VLBI flux-density variations of component B are shown in Fig. 2 (during the first year, the observed flux-density is reduced by synchrotron self absorption).

  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{10101fig2.eps}\end{figure} Figure 2: The flux-density of component B; VLBI observations are shown with dots and the calculated flux-density is indicated by the continuous line. The derived duration of the emission of the relativistic particles responsible for the radio emission is $\tau =26$ yrs.
Open with DEXTER

If the typical time scale of the duration of the emission of the relativistic particles for components A and B is found to be about 25 yrs, for the other components that have smaller bulk Lorentz factors it would be about 10 yrs (see Table 5). For the case of PKS 0420-014, we determined the bulk Lorentz factors of components C, D and E using their mean apparent speeds obtained from VLBI observations. We show the core separations of components B, C, D and E in Figs. 3-5. A newly ejected VLBI component first looks like a point component and then becomes elongated along the beam. So we expect to see a trajectory corresponding to a point source at the beginning and then a mean trajectory corresponding to an elongated component later. In Figs. 3-5 we plot the core separations as a function of the time corresponding both to a point component (dashed line) and to a component having an extended emission duration as given in Table 5, i.e., about 25 yrs for components A and B and about 10yrs for the other components (continuous line).

 

 
Table 5: Characteristics of the VLBI components of PKS 0420-014.
Component $\gamma_{\rm b}$ Duration $\Delta t$ Birth
A 10 ${\simeq}25$ yr 1.22 1991.13
B 8 ${\simeq}25$ yr 3.14 1989.21
C 5.2 ${\simeq}10$ yr 4.54 1987.81
D 3.9 ${\simeq}10$ yr ${\simeq}8.1$ ${\simeq}1984.3$
E 2.5 ${\simeq}10$ yr ${\simeq}18.8$ ${\simeq}1973.5$



  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{10101fig3.eps}\end{figure} Figure 3: Core separation of component B. The dashed line corresponds to a point component and the continuous line corresponds to a component which duration of the relativistic particles is given in Table 5.
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{10101fig4.eps}\end{figure} Figure 4: Core separation of component C.
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{10101fig5.eps}\end{figure} Figure 5: Core separations of components D and E.
Open with DEXTER

The variations of the apparent speeds of the VLBI components A, B and C are shown in Fig. 6. We mention an important consequence of the geometry determined previously. If we increase the bulk Lorentz factor to $\gamma_{\max}=20$, the maximum apparent speed is limited to $v_{\rm app}\simeq10\;c$. So a change of the bulk Lorentz factor between 10 and 20 does not change the apparent speed predicted by the model significantly. Should the apparent speed of a new component increase towards values of $10\;c$, this could be an indication for a changing geometry of the source; e.g., the angle $\Omega_{\circ}$ between the accretion disk and the plane of the two black holes can become smaller, producing a smaller opening angle of the ejection cone of the newly ejected components.

  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{10101fig6.eps} %
\end{figure} Figure 6: Apparent speeds of the VLBI components A, B and C. The asymtotic apparent speed corresponding to $\gamma_{\max}=20$ is also shown.
Open with DEXTER

We determined the birth dates of the components and calculated the ages corresponding to epoch 1992.35. We list them together with the bulk Lorentz factors and the duration of the emission of the relativistic particles in Table 5. We obtained estimates on the age of component A from the age of component B determined by the VLBI observations of the core separation and applying the delay of 1.92 yrs between peaks corresponding to each component, as deduced from the optical light curve (see Fig. 1). The times given in Table 5 correspond to observed times. The proper times of the VLBI components A, B, C and D have been plotted in Fig. 7. From the observed age of the components we can deduce the proper times corresponding to the beginning of the emission of relativistic particles for each VLBI component. The proper times corresponding to the beginning of the emission of the VLBI components are about 20 yrs for component A, 50 yrs for component B, 60 yrs for component C, and 75 yrs for component D. The determination of the proper times of the VLBI components serves to check the consistency of the bulk Lorentz factors with the durations of the emission of relativistic particles found for the components. Indeed, the fit using an arbitrary value of $\gamma_{\rm B}$, different from $\gamma_{{\rm B}}=8$ as discussed in Sect. 4.2 yields bulk Lorentz factors and durations for the different VLBI components that lead to incompatible proper times for the different components. Let us remark that for PKS 0420-014, the time delay between the ejection of new VLBI components is not related a priori to the rotation period of the binary system, which is about 150 yrs. In the case of component B, the beginning of the emission of the e$^{\pm}$ responsible for the radio emission is 1989.21 (see Table 5) and the beginning of the emission of the e$^{\pm}$ responsible for the optical emission is about 1989.8. This difference in the observed times between the ejection of the e$^{\pm}$ responsible for the radio and the optical emissions corresponds to about 10 years of proper time (see Fig. 7), and has to be compared to the duration of the emission of the e$^{\pm}$responsible for the radio emission, which is about 25 yrs of proper time. As indicated previously, this difference between these two times indicates that the variability of PKS 0420-014 can be explained by a long burst of relativistic e$^{\pm}$ responsible for the synchrotron radio emission and a shorter burst of very energetic relativistic e$^{\pm}$ responsible for the synchrotron optical emission, where the emission of the e$^{\pm}$ responsible for the synchrotron radio emission precedes the short burst of more energetic e$^{\pm}$ responsible for the optical synchrotron emission.

  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{10101fig7.eps}\end{figure} Figure 7: The proper times of VLBI components A, B, C and D. Crosses indicates the begining of the emission of the relativistic particles for each VLBI component.
Open with DEXTER

5 Discussion

5.1 Extragalactic radio sources

Helical trajectories and the flaring behavior seen in the optical light curve of PKS 0420-014 have been explained in the framework of the lighthouse model (Camenzind & Krockenberger 1992; Wagner et al. 1995). However, in this scenario the mass of the central object is found to be very large, i.e. $M\geq10^{10}~M_\odot$.

The BBH model we propose requires less massive black holes. Due to the relativistic speed of the components, the observed time is shortened, so that the binary period and the precession period can produce effects observable on time scales ranging from 1 to 10 yrs. The optical light curve can show peaks separated by about one year (Fig. 1). The observed separation of the light-curve peaks in time is related to the rotation period of the binary system, but is not simply the rotation period of the BBH system shortened by the relativistic speed (one has to take into account the effect of the propagation speed of the beam perturbations, i.e., the Alfvén speed). However, the observed delay between the emission of VLBI components and the optical brightness variations is not necessarily related a priori to the rotation period of the binary system. Further, the amplitude of the peaks of the optical light curve depends strongly of the ratio m1/m2. In the case of PKS 0420-014, we have $m_{1}/m_{2}\simeq0.3$ and the amplitude of the peaks is about 2 magnitudes. However, radio sources containing a binary system whose mass ratio is $m_{1}/m_{2} \gtrsim10$ will show helical VLBI trajectories but optical light curves with peaks of very small amplitude after a new VLBI component is ejected. To fit the VLBI observations, we assumed that the bulk Lorentz factors of the various components were different. This is possible if the environment of the clouds emitting the photons responsible for the Compton drag changes with time. As indicated previously, within the geometry we determined the ejected components cannot have apparent speeds larger than $v_{\rm app}\simeq10\;c$. If VLBI observations of new components of PKS 0420-014 show apparent speeds exceeding $10\;c$, this can indicate a change of the geometry, for instance the angle $\Omega_{\circ}$ of the accretion disk with the rotation plane of the binary system could be decreasing.

An important consequence that we can deduce from the difference between the optical light curve and the flux-density variations of the ejected VLBI components is the following. The optical light curve shows peaks of about 2 magnitudes and the radio flux density evolution does not show similar peaks during the five first years. This can be understood if the relativistic e$^{\pm}$ responsible for the radio emission are ejected over a much longer time than those responsible for the optical emission. This points out an important characteristic of the variability of nuclei of extragalactic radio sources: the emission of the e$^{\pm}$ (more precisely: the density and the energy of the emitted e$^{\pm}$) is strongly variable with time. As an important consequence, not all radio sources necessarily show relativistic ejections of e$^{\pm}$, and subluminal radio sources are observable if the mildly relativistic e $^{-}-{\rm p}$ jet contains relativistic e- accelerated by magnetic field reconnection close to the accretion disk. Another important consequence found from the fit of the optical light curve and the radio data is that the burst of the most energetic e$^{\pm}$ responsible for the optical emission occurs after the beginning of the emission of the e$^{\pm}$responsible for the radio emission. VLBI observations of compact radio sources at mm-wavelengths show either helically curved paths within one milliarcsecond separation of the core (see e.g., Krichbaum et al. 1993) or jets with small changes in PA due to intrinsic bends in the first milliarcsecond (Rantakyrö et al. 1998). These observations can be explained if the opening angle of the ejection cone, which is equal to the angle between the accretion disk and the plane of the orbit of the two black holes, is large. The angle found for PKS 0420-014 is particularly large. Britzen et al. (2000b) found a tendency for the gamma-bright sources to have slightly higher curvature of their VLBI trajectories than for the gamma-weak sources.

The observations of PKS 0420-014 and some other compact sources suggest that their nuclei contain a BBH system. More generally, we will assume that extragalactic radio sources are associated with galaxies formed after the merging of galaxies and that the formation of extragalactic radio sources is related to the presence of binary black hole systems in their nuclei. This hypothesis could then explain why

Generally, the binary system contains a single radio source. However, there exists at least one radio source that contains a binary system with two compact sources: it is 3C 75, whose VLA map shows 4 jets (Owen et al. 1985).
Finally, let us mention that the precession of the accretion disk could be explained by the Lense-Thirring effect due to the rotation of a massive central black hole (Landau & Lifshitz 1989). The precession of the disk is given by
$\displaystyle %
\Omega\simeq\frac{4GM^{\prime}}{c^{2}R_{\rm d}^{3}},$     (14)

where $M^{\prime}$ is the angular momentum of the black hole. To get a precession period $T_{\rm p}\simeq10\,000$ yrs with a mean disk radius $R_{\rm d}
\simeq0.0050$ pc, in the case of PKS 0420-014 the angular momentum of the black hole should be $M^{\prime}\simeq2.4\times10^{65}$ (g cm s-1). As the angular momentum $M^{\prime}$ of a black hole of mass m is bounded by
$\displaystyle %
M^{\prime}\leq\frac{1}{2}cm\,R_{\rm g},$     (15)

where $R_{\rm g}$ is the Scharschild radius, the minimum mass to explain the precession of the disk of PKS 0420-014 is
$\displaystyle %
m\geq\sqrt{\frac{cM^{\prime}}{G}}\simeq1.6\times10^{8}\;M_{\odot}.$     (16)

This value is in principle possible and Wagner et al. (1995) explained the optical light curve of PKS 0420-014 and the correlation with elliptical galaxies using a massive spinning BH. However, as indicated previously, the BBH system requires less massive black holes and can explain the difference between the shapes of the optical light curve and the radio light curve, the association of radio sources with elliptical galaxies and the small fraction of radio quasars.

5.2 Galactic radio sources

Galactic radio sources that contain radio jets are especially suited to test the "binary hypothesis'' for the formation of extragalactic radio sources and the two-fluid model of radio sources. Ten galactic radio sources have shown evidence of relativistic jets (Mirabel & Rodríguez 1999). Some like SS 433 (Romney et al. 1987) and GRO J1655-40 (Hjellming & Rupen 1995) show helical jet paths or "wiggles'' on the VLBI scale. Moreover, we can observe:

-
subluminal sources, i.e., SS 433 (Vermeulen et al. 1993);

-
superluminal sources, i.e., GRS 1915+105 (Mirabel & Rodríguez 1994);

-
sources with two different speeds, i.e., Cygnus X3.
It is remarkable that galactic radio sources show the same characteristics as extragalactic radio sources, i.e., subluminal and superluminal speeds and helical VLBI paths. In addition, all of them belong to binary systems.

5.3 Emission of gravitational waves

Let us finish by discussing an important consequence of the fact that all radio sources with jets might be associated with binary systems. These systems will emit gravitational waves and may be detectable by the LISA experiment. Assuming that the two black holes are moving circularly, then the energy losses of the system due to the emission of gravitational waves are given by (Landau & Lifshitz 1989)

$\displaystyle %
\frac{{\rm d}\mathcal{E}}{{\rm d}t}=-\frac{32}{5}\frac{G^{4}m_{1}^{2}m_{2}^{2}(m_{1}%
+m_{2})}{c^{5}r^{5}},$     (17)

where
$\displaystyle %
\mathcal{E}=-\frac{Gm_{1}m_{2}}{2r}\;,$     (18)

is the Newtonian energy of the system. The system continuously loses orbital energy and the coalescing time is given by the formula
$\displaystyle %
T_{{\rm coal}}=\frac{5}{256}\frac{c^{5}r_{\circ}^{4}}{G^{3}m_{1}m_{2}%
(m_{1}+m_{2})},$     (19)

where $r_{\rm o}$ is the present separation between the two black holes. For the given system using the values of Table 4, the coalescence time is $T_{{\rm coal}}\simeq9\times10^{9}$ yrs. We should point out, that in the formula above we have not included the cosmological corrections, this means a multiplication of the above results by (1+z) and possible spin effects which will not significantly affect the duration of the evolution. Practically, the evolution of the system can be described by the standard Newtonian and post-Newtonian relations (inspiral phase) using the weak-field approximation. Only during the very end of its life (merging phase) will strong field effects become important, but this will be only for the last few cycles. This phase is not yet well understood but as we discuss later, the signal will be quite strong and the a priori knowledge of the waveform will not be needed for the detection. Finally, right after the merging there will be an observable quasi-normal mode ringing of the final BH, which is well understood from perturbation theory.

We can estimate the typical timespan in which this system will be observable by LISA. If we assume an optimistic lower limit in the bandwidth of LISA of $5\times10^{-6}$ Hz, then the corresponding period of the system will be $1.3\times10^{7}$ s and the separation between the 2 black holes will be 0.00037 pc. The collision of such a system occurs after about 2.6 yrs, so due to its huge mass the binary system will spend a very short period of its life in the expected bandwidth. Still during the final phase the signal will be detectable. Even if the signal of the inspiral phase is outside the bandwidth of LISA (Larson et al. 1999), the signal from the merging phase will be detectable. Finally, when the two BHs coalesce, the resulting BH will have a total mass ${\simeq}2.8\times10^{8}~M_{\odot}$ and the quasinormal ringing will be at a frequency of $4.3\times10^{-5}$ Hz (Kokkotas & Schmidt 1999).

Assuming that the lifetime of a active radio galaxy is about 107 yrs, there will be about 103 times more galaxies containing a binary black hole system than active radio galaxies. As the space density of active radio galaxies is about $10^{-3}/{\rm Mpc}^{3}$ (Condon 1996), within a sphere of 104 Mpc (which corresponds to a redshift $z\mathrel{\mathchoice {\vcenter{\offinterlineskip\halign{\hfil
$\displaystyle ..., assuming $H_{0}\simeq70$ km s-1 Mpc-1 and q0=0.5) we will find about $4\times10^{9}$ galaxies containing a binary system of massive black holes. Thus, at any moment we can expect to have about one galaxy that contains a binary system of massive black holes in a phase that should be observable with LISA. Let us mention an important characteristic of the evolution of extragalactic radio sources (Roland et al. 1985; Roland & Rhee 1989; Roland et al. 1990). When the nucleus stops its activity, the evolution of extended lobes of radio sources associated with clusters of galaxies is dominated by synchrotron and inverse Compton losses rather than diffusion losses and then one can observe that the radio spectrum evolves following the two Kardashev laws (when expansion loses are stopped by the pressure of the intergalactic medium, an equilibrium between the injection of new material from the nucleus and radiation losses is reached and a first break happens in the radio spectrum i.e. the radio spectrum is characterized by the two spectral indices $\alpha_{\rm o}=0.8$ and $\alpha_{1}=\alpha_{\rm o}+0.5=1.3$, latter when the nucleus stops to inject new material in the extended lobe a second break happens in the radio spectrum and the high frequency spectrum is characterized by the spectral index $\alpha_{2}=(4/3)\times\alpha_{1}+1=2.7$). The age of the extended lobes without input of new material from the nucleus of the galaxy is of order of 108 yrs. These observations exclude a scenario where the nucleus restarts its activity every 107 yrs. Moreover, the energy contained in the extended lobes can reach 1061 ergs, i.e. corresponds to about $6\times10^{6}~M_{\odot}\times c^{2}$, so the nucleus cannot have a recurrent activity every say 107 yrs and the radio source cannot have 103 phases of activity.

Although these estimates are quite uncertain, massive BBH systems will be one of the primary targets for LISA (Larson et al. 1999). Moreover, a simple estimation shows that such systems, due to their large masses, will be observable by LISA regardless of their distance. Indeed, the amplitude h of the gravitational wave on earth is

$\displaystyle %
h=0.237\times\frac{\mu^{1/2}M^{1/3}}{Df^{1/6}},$     (20)

where $\mu$ is the reduced mass, M the total mass and D the distance of the system and f is the observation frequency.

For massive BBH systems we obtain:

$\displaystyle h =1.87\times10^{-16}\left( \frac{\displaystyle \mu}{\displaystyl...
...ft( \frac{\displaystyle 10^{-4}\;{\rm Hz}}
{\displaystyle f}\right) ^{1/6}\cdot$     (21)

Thus assuming a sensitivity of the order of 10-20 for LISA, these sources will be detectable up to 106 Gpc, i.e. for any source with a known redshift.

6 Conclusion

The structural evolution of the quasar PKS 0420-014 can be modelled by a precession of the accretion disk. We explain the precession of the accretion disk by the presence of a binary system of black holes in the nucleus. We show that the light curve of the optical outburst associated with the birth of a new VLBI component can be explained by the motion of the beam-emitting black hole around the center of gravity of the BBH system. The amplitude of the optical peaks is very sensitive to the ratio of the masses of the binary system.

More generally, we propose that the formation of galactic and extragalactic compact radio sources containing jets can be related to the existence of binary systems. The existence of binary systems due to mergers in extragalactic radio sources can provide a natural explanation to understand why only about 5-10% of the quasars are radio sources and why extragalactic radio sources are associated with elliptical galaxies. The hypothesized existence of binary systems in all radio sources provides possible targets for the research of gravitational waves using the LISA experiment.

Acknowledgements
Part of this work was supported by the European Commission, TMR Programme, Research Network Contract ERBFMRXCT97-0034 CERES.

7 Appendix I: Precession in the three-body problem

7.1 Expansion in Jacobi coordinates

Let P0,P1,P2 be three objects of mass m0,m1,m2, and ${\vec u}_{i}=\overrightarrow{OP_{i}}$, where O is at the origin of an inertial reference frame. Note that this treatment places no assumptions or specific assignments on the masses 0, 1, 2; In the body of the text, the black hole that emits the VLBI components corresponds to mass 0, the accretion disk to mass 1, and the other black hole to mass 2.

In the canonical variables $({\vec u}_{i},{\vec{\tilde{u}}}_{i}=m_{i}\dot{{\vec u}}
_{i})_{i=0,2}$, their motion is described by the Hamiltonian

$\displaystyle %
\mathcal{H}=T+U=\frac{1}{2}
\sum_{i=0}^{2}\frac{\left\Vert {\ve...
...}_{i}\right\Vert ^{2}}{m_{i}}-G\sum_{0\leq i<j}
\frac{m_{i}m_{j}}{\Delta_{ij}},$     (22)

where $\Delta_{ij}=\left\Vert {\vec u}_{i}-{\vec u}_{j}\right\Vert $, and Gis the gravitational constant.

The Jacobi coordinates are obtained through the linear transformation A:

$\displaystyle \left[ \begin{array}{c} {\vec r}_0 \\  {\vec r}_1 \\  {\vec r}_2\...
...[
\begin{array}{c}
{\vec u}_0 \\  {\vec u}_1 \\  {\vec u}_2 \end{array}\right],$     (23)

where M01=m0+m1 and M02=m0+m1+m2. This linear transformation can be extended in a canonical transformation by applying tA-1 to the associated momentum, which gives
$\displaystyle %
\left[ \begin{array}{c} \vec{\tilde r}_0 \\  \vec{\tilde r}_1 \...
...\vec{\tilde u}_0 \\  \vec{\tilde u}_1 \\  \vec{\tilde u}_2 \end{array}\right] .$     (24)

If we use a reference frame with origin at the center of mass G of the system, we obtain ${\vec r}_{0}=0$, ${\vec{\tilde{r}}}_{0}=0$, and thus
$\displaystyle %
{\vec u}_0$=$\displaystyle -{m_1}/{ M_01}{\vec r}_1 -{m_2}/{ M_02}
{\vec r}_2$  
$\displaystyle {\vec u}_1$=$\displaystyle {m_0}/{ M_01}{\vec r}_1 -{m_2}/{ M_02} {\vec r}_2$  
$\displaystyle {\vec u}_2$=$\displaystyle { M_01}/{ M_02}{\vec r}_2$ (25)


$\displaystyle %
\vec{\tilde u}_0 = -\vec{\tilde r}_1 -{m_0}/{ M_01}\vec{\tilde r}_2$      
$\displaystyle \vec{\tilde u}_1 =
\vec{\tilde r}_1 -{m_1}/{ M_01} \vec{\tilde r}_2$      
$\displaystyle \vec{\tilde u}_2 = \vec{\tilde r}_2.$     (26)

7.2 The Hamiltonian in Jacobi coordinates

Using these formulas, we obtain for the kinetic energy

 
$\displaystyle %
T = \frac{1}{2}\frac{M_{01}}{m_{0} m_{1}}\left\Vert
\mathbf{\ve...
...{02}}}{{ M_{01} m_{2}}}}\left\Vert \mathbf{\vec{\tilde r}}_{2}\right\Vert ^{2}.$     (27)

The computation of the potential energy U is a little more complex. The goal is to isolate in its expression the part that can be associated with the kinetic energy (Eq. (27)) in order to form "disjoints-Kepler's'' problems, and to treat the remainder as a perturbation. The natural problem in Jacobi coordinates consists of a mass mi, attracted by a body of mass $m_{0}+m_{1}+\cdots+m_{i-1}$, located at the center of mass Gi-1 of $m_{0},\dots,m_{i-1}$, (Fig. 8), which corresponds to the potential (see Laskar 1990)
$\displaystyle %
U_{0}=-G\frac{m_{0}m_{1}}{r_{1}}
-G\frac{M_{01}m_{2}}{r_{2}},$     (28)

where $r_{k}=\left\Vert {\vec r}_{k}\right\Vert $. The complete potential can then be written as U=U0+U1 with, as $\Delta_{01}=v_{1}$,
 
$\displaystyle U_1 = G \frac{ M_01m_2}{r_2} -G \frac{m_0m_2}
{\Delta_02}-G \frac{m_1m_2}{\Delta_12}$ $\textstyle = Gm_0m_2 \left[ \frac{1}
{r_2}- \frac{1}{\Delta_02}\right] + Gm_1m_2 \left[ \frac{1}{r_2} - \frac
{1}{\Delta_12}\right] ,$   (29)

with
$\displaystyle %
\Delta_02 = \left\Vert{\vec u}_2-{\vec u}_0\right\Vert = \left\...
...}_1\right\Vert = \left\Vert{\vec r}_2
-\frac{m_0}{ M_01}{\vec r}_1\right\Vert .$     (30)

Evaluating Eq. (25) yields
$\displaystyle %
\Delta_02 = r_2 \left(1+2\frac{m_1}{ M_01}\ \rho\cos S +
\frac{m_1^2}{ M_01^2}\ \rho^2\right)^1/2$     (31)
$\displaystyle \Delta_12 = r_2 \left(1-2\frac{m_0}{
M_01}\ \rho\cos S + \frac{m_0^2}{ M_01^2}\ \rho^2\right)^1/2,$    

where S denotes the angle $({\vec v}_{1},{\vec v}_{2})$, and $\rho
=r_{1}/r_{2}$. As $\rho$ is small, we can expand (Eq. (29)) in powers of $\rho$. The two first terms of this series cancel and only the terms starting with the contribution of the quadrupolar momentum remain, i.e.,
$\displaystyle %
U_{1}=-\frac{Gm_{0}m_{1}m_{2}}{r_{2}M_{01}}\rho^{2}\ \frac{1}{2}(3\cos^{2}S-1)\ +O(\rho^{3}).$     (32)


  \begin{figure}
\includegraphics[width=8cm,clip]{10101fig8.eps}\end{figure} Figure 8: Jacobi coordinates.
Open with DEXTER

7.3 Poincaré variables and equations of motion


  \begin{figure}
\par\includegraphics[width=5.5cm,clip]{10101fig9.eps}\end{figure} Figure 9: Elliptical elements. The longitude of te ascending node, $\Omega $, and yhe inclination, i, define the position of the plane of the orbit; the semimajor axis, a and the eccentricity, e, the shape of the ellipse; the argument of perihelion $\omega $, the orientation of the orbit in its plane, while the mean anomaly, M, which is proportional to the area swept from the perihelion direction fixes the position of the body P on its orbit.
Open with DEXTER

The Hamiltonian $\mathcal{H}_{0}=T+U_{0}$ corresponds to a disjoint sum of two Keplerian problems: motion of m1 around m0, and motion of m2around the center of mass of m0 and m1. These two-body problems are expressed naturally in elliptical variables (ak: semi-major axis, ek: eccentricity, ik: inclination, Mk: mean anomaly, $\omega_{k}$: argument of perihelion, $\Omega_{k}$: longitude of the node, for k=1,2). In order to have canonical variables, we can use the Delaunay elliptical variables, but here we will prefer the rectangular complex variables of Poincaré defined from the classical elliptical elements as (Laskar 1990; Laskar & Robutel 1995):

$\displaystyle %
x_k = \sqrt{\Lambda_k} \sqrt{1-\sqrt{1-e_k^2}} \exp
({i\varpi_k})$      
$\displaystyle y_k = \sqrt{\Lambda_k} \sqrt{\sqrt{1-e_k^2}(1-\cos i_k)}
\exp({i\Omega_k}),$     (33)

where
$\displaystyle %
\Lambda_{k}=\beta_{k}\sqrt{\mu_{k}a_{k}}\ ;\qquad\mu_{k}=GM_{0k}$     (34)

and
$\displaystyle %
\beta_{1}=\frac{m_{0}m_{1}}{M_{01}};\qquad\beta_{2}=\frac{M_{01}m_{2}}{M_{02}}\cdot$     (35)

The angle $\varpi_{k}=\omega_{k}+\Omega_{k}$ is the longitude of perihelion. The canonical conjugate variables are then $(\Lambda_{k},$ $\lambda_{k},$xk, $-i\bar{x}_{k},$ yk, $-i\bar{y}_{k})$, where $\lambda_{k}
=\varpi_{k}+M_{k}$ is the mean longitude, and the associated Hamilton equations of motion become:
 
$\displaystyle %
\frac{{\rm d}\l _k}{{\rm d}t}= \frac{\partial H}{\partial \Lamb...
...} ;\qquad\frac{{\rm d}\Lambda_k}{{\rm d}t} = -\frac{\partial H}{\partial \l _k}$      
$\displaystyle \frac{{\rm d}{\bar x}_k}{{\rm d}t}= i\frac{\partial H}{\partial x_k} ;\qquad\frac{{\rm d}x_k}{{\rm d}t} = -i
\frac{\partial H}{\partial {\bar x}_k}$      
$\displaystyle \frac{{\rm d}{\bar y}_k}{{\rm d}t}= i\frac{\partial H}{\partial y...
... \qquad\frac{{\rm d}y_k}{{\rm d}t} = -i\frac{\partial H}{\partial {\bar y}_k}
,$     (36)

where the fully perturbed Hamiltonian is $\mathcal{H}=\mathcal{H}_{0}+U_{1}$, with
$\displaystyle %
\mathcal{H}_{0}=-\frac{\mu_{1}^{2}\beta_{1}^{3}}{2\Lambda_{1}^{2}}- \frac{\mu_{2}^{2}\beta_{2}^{3}}{2\Lambda_{2}^{2}}$     (37)

and
$\displaystyle %
U_{1}=-\frac{Gm_{0}m_{1}m_{2}}{M_{01}}\frac{a_{1}^{2}}{a_{2}^{3...
...}\right) ^{2}\left( \frac{a_{2}}{r_{2}}\right) ^{3}\ \frac{1}{2}(3\cos^{2}S-1).$       (38)

7.4 Expansion of the secular Hamiltonian

Now we use the classical expansions of r/a, a/r, and $\cos S$ in elliptical variables in order to expand $\mathcal{H}$ in elliptical variables. We retain only up to second order terms in this expansion, and then only the secular part of the Hamiltonian, i.e., the terms that do not depend of the longitudes. This system will thus reflect the long term precession motion of the orbits, and not the short term oscillations of small amplitude. With $\gamma_{k}=\sin i_{k}/2$, we obtain

 \begin{displaymath}%
U_1s = -{\cal K}\bigg(\frac{1}{4}+\frac{3}{8}e_1^2+\frac
{3...
...}\gamma_2^2 +3\gamma_1\gamma
_2\cos(\Omega_1-\Omega_2) \bigg),
\end{displaymath} (39)

where
$\displaystyle %
\mathcal{K}=\frac{Gm_{0}m_{1}m_{2}}{M_{01}}\frac{a_{1}^{2}}{a_{2}^{3}}\cdot$     (40)

In rectangular Poincaré variables this becomes

\begin{displaymath}%
U_1s = -{\cal K}\bigg( + \frac{1}{4} + \frac{3}{4}\frac{x_1...
...}{4} \frac{y_1 {\bar y}_2}{\sqrt{\L _1\Lambda_2}} \bigg)\cdot
\end{displaymath} (41)

7.5 Motion of the node and perihelion

As the secular Hamiltonian does not depend on the longitudes, from e.g., (37), we have immediately

$\displaystyle %
\frac{{\rm d} \Lambda_{k}}{{\rm d}t} =
-\frac{\partial H}{\partial\lambda_{k}} = 0,$     (42)

and thus the semi-major axes a1, a2 are constant. For the variations of the other variables ( $e_{k},i_{k},\varpi_{k}, \Omega_{k}$), we use the expressions in term of the rectangular variables ( xk,yk) and we obtain the linear equations

 \begin{displaymath}%
\frac{{\rm d} }{{\rm d}t} \left[
\begin{array}{c}x_1 \\ x_2...
...}\right] \left[ \begin{array}{c}x_1 \\ x_2\end{array}\right] \
\end{displaymath} (43)


 \begin{displaymath}%
\frac{\rm d}{{\rm d}t} \left[
\begin{array}{c}y_1 \\ y_2\en...
...}\right] \left[ \begin{array}{c}y_1 \\ y_2\end{array}\right] .
\end{displaymath} (44)

We also have
$\displaystyle %
\Lambda_1$ $\textstyle = \beta_1\sqrt{\mu_1 a_1} = \frac{m_0m_1}{
M_01}n_1a_1^2$    
$\displaystyle \Lambda_2$ $\textstyle = \beta_2\sqrt{\mu_2 a_2} = \frac{ M_01m_2}{
M_02}n_2a_2^2 ,$   (45)

where nk are the mean motion of the Keplerian problems, defined by $n_{k}^{2}a_{k}^{3}=\mu_{k}$ (or $n_{k}=2\pi/T_{k}$, where Tk are the periods of the Keplerian orbits). Thus
$\displaystyle %
\frac{{\cal K}}{\Lambda_1}$ $\textstyle =\frac{m_2}{M_02}\left( \frac
{n_2}{n_1}\right) n_2$    
$\displaystyle \frac{{\cal K}}{\Lambda_2}$ $\textstyle =\frac{m_0m_1}{M_01^2}\left(
\frac{a_1}{a_2}\right) ^2n_2 .$   (46)

The solutions for the eccentricities (x1,x2) are easily obtained as the system of formulas (44) is diagonal in this approximation. We have
$\displaystyle %
x_{1}=x_{10}\exp({ig_{1}t});\quad x_{2}=x_{20}\exp({ig_{2}t}),$     (47)

with
 
$\displaystyle %
g_{1}=\frac{3}{4}\frac{\mathcal{K}}{\Lambda_{1}}\ ;\qquad
g_{2}=\frac{3}{4}\frac{\mathcal{K}}{\Lambda_{2}},$     (48)

that is, in the original elliptical variables

\begin{displaymath}%
\left\{ \begin{array}{llll}
a_1 & = a_10={\rm cte}; & a_2 ...
...\quad & \varpi_2 & = \varpi_20 + g_2\, t \,.\end{array}\right.
\end{displaymath} (49)

One should thus notice that in this approximation, not only are the semi-major axes constant, but so are the eccentricities. We have a precession of the perihelions of the orbits in the plane, without any change in their shapes. In order to compute the motion of the inclination and nodes, we need to diagonalize the linear system Sy. Since the Sy matrix is symmetric, this can be done with orthogonal transformations, and since it is singular, s0=0 is one of the eigenvalues, and the other one, s, is equal to its trace
$\displaystyle %
s_{0} = 0\ ; \qquad s = -\frac{3}{4}\left( \frac{\mathcal{K}}{\Lambda_{1}} +\frac{\mathcal{K}}{\Lambda_{2}}\right) \cdot$     (50)

We use the orthogonal matrix S of the eigenvectors for the change of variables

\begin{displaymath}%
S =\frac{1}{\sqrt{ 1+\delta}}
\left[
\begin{array}{cc} 1 ...
..., \qquad
\hbox{with }
\delta=\frac{\Lambda_{1}}{\Lambda_{2}},
\end{displaymath} (51)

and we introduce the proper variables (w1,w2) defined as [y] = S[w], that is

\begin{displaymath}%
\left[\begin{array}{c}w_1 \\ w_2\end{array}\right] =S^-1\le...
...ray}\right]\left[\begin{array}{c}y_1 \\ y_2\end{array}\right].
\end{displaymath} (52)

With these new variables, the system Sy becomes

\begin{displaymath}%
\frac{\rm d}{{\rm d}t}\left[
\begin{array}{c}w_1 \\ w_2\end...
...}\right] \left[ \begin{array}{c}w_1 \\ w_2\end{array}\right] ,
\end{displaymath} (53)

and we easily obtain the solutions
$\displaystyle %
w_{2}(t)={\rm cte}=w_{20}\ ;\qquad w_{1}(t)=w_{10}\exp({i\,st}).$     (54)

The solution in the original variables ( y1,y2) is then

\begin{displaymath}%
\left[ \begin{array}{c} y_1 \\ y_2\end{array}\right] =\frac...
...\right] \left[
\begin{array}{c}
w_1 \\ w_2\end{array}\right] .
\end{displaymath} (55)

7.6 Simplifications when $\Lambda _{1} << \Lambda _{2}$

When $\Lambda _{1} << \Lambda _{2}$, we have $\delta<< 1$, and the solutions becomes $y_{k} \approx w_{k}$, which gives

\begin{displaymath}%
\left\{ \begin{array}{llll} i_1 & \approx i_10={\rm cte} ; ...
...& \Omega_2 &
\approx\Omega_20 = {\rm cte} .
\end{array}\right.
\end{displaymath} (56)

8 Appendix II: The $\mathsf{d{\it z}/d{\it t}}$ equation

In the case of the binary system, we can obtain the ${\rm d}z/{\rm d}t$ equation similar to Eq. (6) in the following way. Calling

\begin{eqnarray*}\eta(t) & =&\omega_{\rm b}t-k_{\rm b}z_{\rm c}(t)\;\\
\psi(t) & =&\omega_{\rm p}t-k_{\rm p}z_{\rm c}(t)\;,
\end{eqnarray*}


Eqs. (9) and (10) become

\begin{eqnarray*}x_{\rm c}&=&x_{1}(t)\frac{\exp(i\eta(t))+\exp(-i\eta(t))}{2}\\
&&+R(z_{\rm c}(t))\frac{\exp(i\psi(t))+\exp(-i\psi(t))}{2}\;
\end{eqnarray*}



\begin{eqnarray*}y_{\rm c}&=&y_{1}(t)\frac{\exp(i\eta(t))-\exp(-i\eta(t))}{2i}\\
&&+R(z_{\rm c}(t))\frac{\exp(i\psi(t))-\exp(-i\psi(t))}{2i}\cdot
\end{eqnarray*}


Then we calculate $({\rm d}x_{\rm c}/{\rm d}t)^{2}+({\rm d}y_{\rm c}/{\rm d}t)^{2}$ as

\begin{eqnarray*}\left( \frac{{\rm d}x_{\rm c}}{{\rm d}t}\right) ^{2}+\left( \fr...
...x_{\rm c}}{{\rm d}t}+i\frac{{\rm d}y_{\rm c}}{{\rm d}t}\right)},
\end{eqnarray*}


and Eq. (5) can be written in the form

\begin{eqnarray*}A\left( \frac{{\rm d}z_{\rm c}}{{\rm d}t}\right) ^{2}+B\frac{{\rm d}z_{\rm c}}{{\rm d}t}+C=0,
\end{eqnarray*}


where
A = $\displaystyle 1+\frac{\omega_{\rm p}^{2}R^{2}(z_{\rm c})}{V_{\rm A}^{2}}+\left( \frac{{\rm d}R}{{\rm d}z_{\rm c}
}\right) ^{2}$  
    $\displaystyle {+}\frac{\omega_{\rm b}^{2}\,x_{1}^{2}(t)}{2V_{\rm A}^{2}}\,[1{-}...
...}\frac
{\,\omega_{\rm b}^{2}\,y_{1}^{2}(t)}{2V_{\rm A}^{2}}[1{+}\cos(2\eta(t))]$  
    $\displaystyle +\frac{\omega_{\rm b}}{V_{\rm A}}\,\left( x_{1}(t)-\,\,y_{1}(t)\right) \frac
{{\rm d}R}{{\rm d}z_{\rm c}}\sin(\eta(t)+\psi(t))$  
    $\displaystyle +\frac{\omega_{\rm b}}{V_{\rm A}}\left( x_{1}(t)+\,\,y_{1}(t)\right) \,\frac
{{\rm d}R}{{\rm d}z_{\rm c}}\sin(\eta(t)-\psi(t))$  
    $\displaystyle +\frac{\omega_{\rm b}\omega_{\rm p}}{V_{\rm A}^{2}}R(z_{\rm c})\left[ x_{1}(t)+\,\,y_{1}
(t)\right] \cos(\eta(t)-\psi(t))$  
    $\displaystyle +\frac{\omega_{\rm b}\omega_{\rm p}}{V_{\rm A}^{2}}R(z_{\rm c})\left[ x_{1}(t)-\,\,y_{1}
(t)\right] \cos(\eta(t)+\psi(t)),$ (57)
B = $\displaystyle -2\frac{\omega_{\rm p}^{2}R^{2}(z_{\rm c})}{V_{\rm A}}-\frac{\omega_{\rm b}^{2}}{V_{\rm A}
}\left(x_{1}^{2}(t)+\,\,y_{1}^{2}(t)\right)$  
    $\displaystyle +\frac{\omega_{\rm b}^{2}}{V_{\rm A}}\left[ x_{1}^{2}(t)-\,\,y_{1}^{2}(t)\right]
\cos(2\eta(t))$  
    $\displaystyle +\frac{\omega_{\rm b}}{V_{\rm A}}\left(x_{1}(t)\frac{{\rm d}x_{1}(t)}{{\rm d}t}-\,\,y_{1}
(t)\frac{{\rm d}y_{1}(t)}{{\rm d}t}\right) \sin(2\eta(t))$  
    $\displaystyle +\omega_{\rm b}\,\left[ y_{1}(t)-\,\,x_{1}(t)\right] \frac{{\rm d}R}{{\rm d}z_{\rm c}}\sin
(\eta(t)+\psi(t))$  
    $\displaystyle +\frac{\omega_{\rm p}R(z_{\rm c})}{V_{\rm A}}\left( \frac{{\rm d}...
...)}{{\rm d}t}-\,\,\frac
{{\rm d}y_{1}(t)}{{\rm d}t}\right) \sin(\eta(t)+\psi(t))$  
    $\displaystyle -\omega_{\rm b}\left[ x_{1}(t)+\,\,y_{1}(t)\right] \frac{{\rm d}R}{{\rm d}z_{\rm c}}\sin
(\eta(t)-\psi(t))$  
    $\displaystyle -\frac{\omega_{\rm p}R(z_{\rm c})}{V_{\rm A}}\left( \frac{{\rm d}...
...)}{{\rm d}t}+\,\,\frac
{{\rm d}y_{1}(t)}{{\rm d}t}\right) \sin(\eta(t)-\psi(t))$  
    $\displaystyle -2\frac{\omega_{\rm b}\omega_{\rm p}R(z_{\rm c})}{V_{\rm A}}\left[ x_{1}(t)+\,\,y_{1}
(t)\right] \cos(\eta(t)-\psi(t))$  
    $\displaystyle +\left( \frac{{\rm d}x_{1}(t)}{{\rm d}t}+\,\,\frac{{\rm d}y_{1}(t)}{{\rm d}t}\right) \frac
{{\rm d}R}{{\rm d}z_{\rm c}}\cos(\eta(t)-\psi(t))$  
    $\displaystyle +2\frac{\omega_{\rm b}\omega_{\rm p}R(z_{\rm c})}{V_{\rm A}}\left[ x_{1}(t)-\,\,y_{1}
(t)\right] \cos(\eta(t)+\psi(t))$  
    $\displaystyle +\left( \frac{{\rm d}x_{1}(t)}{{\rm d}t}-\,\,\frac{{\rm d}y_{1}(t)}{{\rm d}t}\right) \frac
{{\rm d}R}{{\rm d}z_{\rm c}}\cos(\eta(t)+\psi(t)),$ (58)
$\displaystyle {\rm and}$      
C = $\displaystyle \omega_{\rm p}^{2}R^{2}(z_{\rm c})+\frac{\omega_{\rm b}^{2}}{2}\left( x_{1}
^{2}(t)+\,\,y_{1}^{2}(t)\right)$  
    $\displaystyle +\frac{1}{2}\left( \frac{{\rm d}x_{1}(t)}{{\rm d}t}\right) ^{2}+\frac{1}{2}\left(
\frac{{\rm d}y_{1}(t)}{{\rm d}t}\right) ^{2}-v_{\rm c}^{2}$  
    $\displaystyle +\frac{1}{2}\cos(2\eta(t))\times$  
    $\displaystyle \left( \left( \frac{{\rm d}x_{1}(t)}{{\rm d}t}\right) ^{2}-\left(...
...ght) ^{2}-\omega_{\rm b}^{2}\left( x_{1}^{2}(t)-\,\,y_{1}^{2}(t)\right)
\right)$  
    $\displaystyle -\omega_{\rm b}\left( x_{1}(t)\frac{{\rm d}x_{1}(t)}{{\rm d}t}-\,\,y_{1}(t)\frac{{\rm d}y_{1}
(t)}{{\rm d}t}\right) \sin(2\eta(t))$  
    $\displaystyle -\omega_{\rm p}R(z_{\rm c})\left( \frac{{\rm d}x_{1}(t)}{{\rm d}t}-\,\,\frac{{\rm d}y_{1}(t)}
{{\rm d}t}\right) \sin(\eta(t)+\psi(t))$  
    $\displaystyle +\omega_{\rm p}R(z_{\rm c})\left( \frac{{\rm d}x_{1}(t)}{{\rm d}t}+\,\,\frac{{\rm d}y_{1}(t)}
{{\rm d}t}\right) \sin(\eta(t)-\psi(t))$  
    $\displaystyle +\omega_{\rm b}\omega_{\rm p}R(z_{\rm c})\left[ x_{1}(t)+\,\,y_{1}(t)\right] \cos
(\eta(t)-\psi(t))$  
    $\displaystyle -\omega_{\rm b}\omega_{\rm p}R(z_{\rm c})\left[ x_{1}(t)-\,\,y_{1}(t)\right] \cos
(\eta(t)+\psi(t)).$ (59)

References

 
Copyright ESO 2001