Multiple black hole system in 4C31.61 (2201+315)

Modeling trajectories of radio components ejected by the nucleus of 4C31.61 (2201+315) and observed by very long baseline interferometry (VLBI) in the frame of the MOJAVE survey suggests that they are ejected from three different origins that possibly host three different supermassive black holes. These origins correspond to three stationary components, one of which one is the VLBI core. Most of the mass of the nucleus is associated with a supermassive binary black hole system whose separation is ≈0.3 milliarc second, that is, a distance of ≈1.3 parsec and the mass ratio is ≈2. In contrast, the mass ratio with respect to the third black hole is ≈1/100. The three origins lie within 0.6 milliarc second, or a distance of ≈2.6 parsec. Based in this structure of the nucleus, we explain the variations observed in the astrometric coordinate time series obtained from VLBI geodetic surveys. This study shows that it is possible to exploit large MOJAVE-like VLBI databases to propose more insights into the structure of the extragalactic radio sources that are targeted by VLBI in geodetic and astrometry programs.


Introduction
Very long baseline interferometry (VLBI) observations of nuclei of extra galactic radio sources show that the ejection of VLBI components does not follow a straight line, but wiggles. These observations suggest a precession of the accretion disk that can be explained if the nuclei contain either a spinning single black hole (BH) or a binary black hole (BBH) system. We developed a minimization method by first modeling the ejection of a VLBI component by a spinning single BH, which produces the precession of the accretion disk. Then we modeled the ejection of a VLBI component by a BBH system, taking the precession of the accretion disk and the motion of the BHs into account.
When a nucleus contains a BBH system, the two BHs can eject VLBI components, and both can be detected using VLBI observations. In this case, we observe the core and one stationary component, and we observe two families of trajectories. This is, for instance, the case for 3C 279 (Roland et al. 2013) and 1928+738 (Roland et al. 2015). However, it is not always the case, and in some cases, only one BH can be detected using VLBI observations, the other being "inactive".
In this article, we use for the Hubble constant H o ≈ 73 km s −1 Mpc −1 and D L ≈ 1460 Mpc for the luminosity distance of 2201+315 and model the ejection of the VLBI components of 4C31.61 ≡ 2201+315 using data from the Monitoring Of Jets in Active galactic nuclei with VLBA Experiments (MOJAVE), that is, the 15 GHz observations and the definitions of the different components of Table 3 of Lister et al. (2019). The analysis of MOJAVE data with our minimization method suggests that the nucleus of 2201+315 contains three BHs associated with the stationary components C0 (the core), C3, and C9 (see Fig. 1). Most of the mass of the nucleus is associated with the BBH system C0-C9, whose separation is R C0−C9 = 0.296 +0.066 −0.020 mas, corresponding to a projected linear distance = 1.3 +0.26 −0.09 pc, and the mass ratio is M C9 /M CO = 2 ± 0.5. The mass of the third BH is M C3 /(M C0 + M C9 ) ≈ 1/100. The distance between C0 and C3 is R C0−C3 = 0.557 +0.124 −0.124 mas, which translates into ≈2.4 pc. Calling G, the center of gravity of the system C0-C9, the distance between component C3 and the center of gravity G is R G−C3 = 0.369 +0.124 −0.124 mas, that is, ≈1.6 pc (see Table 1). We show that component C6 is ejected from the core C0, component C7 is ejected from the stationary component C3, and component C13 is ejected from the stationary component C9. We deduce the characteristics of the BBH system constituted by the core C0 and the stationary component C9. We also deduce the characteristics of the third BH that is associated with the stationary component C3. Finally, knowing the structure of the nucleus of 2201+315, we discuss the influence of the variations in the ratio of the flux densities of the three BHs and of the ejection of a new VLBI component on the coordinate time series obtained by the analysis of geodetic VLBI observations (Gattano et al. 2018).  The analysis of MOJAVE data with our minimization method suggests that the nucleus of 2201+315 contains three BHs. G is the center of gravity of the BBH system C0-C9 and L1 is the first Lagrange point of the BBH system C0-C9; this is the point where the gravitational forces of the BHs C0 and C9 are equal. In order to obtain a stable solution, the radii of the accretion disks around C0 and C9 have to be smaller than the distances L1-C0 and L1-C9, respectively.

Properties of the BBH system solution
To model the ejection of a VLBI component, we assumed that the VLBI component corresponds to the relativistic ejection of an electron-positron (e − e + ) plasma referred to as the beam. It is surrounded by a subrelativistic electron-proton (e − p) plasma called the jet. The relativistic beam is responsible for the formation of VLBI components, their synchrotron radio emission, and inverse Compton emission in the UV, X-rays, and γ-rays. The subrelativistic electron-proton jet that carries most of the mass and kinetic energy ejected by the BH is responsible for the formation of kiloparsec jets, hot spots, and extended lobes. The magnetic field in the beam and the mixing layer between the beam and the jet is parallel to the flow, and the magnetic field in the jet rapidly becomes toroidal. This model is called the twofluid model ( Fig. 2; see also Roland et al. 2013, and references therein). The beam can propagate along the magnetic lines and is stable if the bulk Lorentz factor, γ c of the ejected VLBI component is γ c < 30 and if the magnetic field in the beam and in the mixing layer is grater than a critical value B c . The beam and the jet are perturbed by the precession of the accretion disk and the motion of the BHs.
When in addition to the radio, optical observations are available that peak in the light curve, this optical emission can be modeled as the synchrotron emission of a point source ejected in the perturbed beam, see Britzen et al. (2001) and Lobanov & Roland (2005). This short burst of very energetic, relativistic e ± plasma is followed immediately by a very long burst of less energetic relativistic e ± plasma. This long burst is modeled as Fig. 2. Two-fluid model. The outflow consists of an e − − e + plasma, the beam, which moves at a highly relativistic speed and is surrounded by an e − − p plasma, and of the jet, which moves at a mildly relativistic speed.
an extended structure along the beam and is responsible for the VLBI radio emission.
In the case of ejection by a spinning single BH, we found in general that the surface χ 2 of the solution that corresponds to the inclination angle is convex, indicating that there is no stable solution. However, we found that it becomes concave when the ejection is caused by a BBH system. For details on the model geometry, the perturbation of the VLBI ejected component, and the coordinates of the ejected VLBI component, see Appendix A.
When the two BHs eject VLBI components, we observe two families of trajectories. We find by fitting the ejection of VLBI components of the two families that if one family is characterized by the mass ratio M 1 /M 2 = a, where M 1 and M 2 are the masses of the two BHs, the second family is characterized by the inverse mass ratio 1/a, showing the consistency of the BBH model. The solution of the fit is not unique and shows a degeneracy; the degeneracy parameter being V a , the propagation speed of the perturbation of the jet and the beam. This means that χ 2 (V a ) = constant, that is, χ 2 (V a ) does not depend on the value of V a when V a varies. It also means that there is a possible range of values for the BBH system period T b and the BBH system mass M 1 + M 2 .

Properties of the radio source 2201+315
The source has been observed using the Very Long Baseline Array (VLBA) in the frame of the MOJAVE survey 1 at 15 GHz (Lister et al. 2019) and at 43 and 86 GHz (Cheng et al. 2018). The positions of the components J1, J2, and J3 detected at 43 and 86 GHz are shown in Fig. 8. The position of 2201+315 provided by the MOJAVE survey is shown in Figs. 3 and 4 with the stacked 15 GHz map and a plot of separation versus time, with the component names as provided by the MOJAVE survey. These observations reveal the stationary radio components C3 and C9, whose flux densities are variable with time. We show that these components can be linked with the BHs that eject VLBI components. Source 2201+315 has been observed during more than 20 years by MOJAVE, and an important property of these components is that they are detected during different periods of time.
Radio source 2201+315 2 is associated with a quasar whose active galactic nucleus (AGN) class is low-spectral peak  <10 14 Hz (LSP), and it is a quasar whose fractional linear optical polarization is consistently below 3% (LQP). Its redshift is z ≈ 0.2947 according to Ho & Kim (2009). Hutchings & Morris (1995) reported that Canada France Hawaii Telescope (CFHT) 3.6 m optical telescope images of 2201+315 indicated an overall elliptical shape, extended normal to the radio axis, but with a radial color gradient (redder with increasing distance), and some irregular knots within the host galaxy. The total magnitudes (R15.6 and V16.6) are close to those quoted in Hewitt & Burbidge (1993) and Hutchings & Neff (1992), and the small differences are probably due to differing extensions of the host galaxy that were measured. Healey et al. (2008) reported an Rmag of M R ≈ 14.33. The Gaia DR2 catalog (Brown et al. 2018) gives a Gmag of 15.3969, a GBPmag of 15.4428, and a GRPmag of 14.9321 (Brown et al. 2018). Gmag (wavelength range 330-1050 nm), GBPmag (blue photometer, wavelength range 330-680 nm) and GRPmag (red photometer, wavelength range 640-1050 nm) are Gaia magnitudes defined in Carrasco et al. (2016). The source has also been detected by the Fermi/LAT -γ-ray space telescope (Acero et al. 2015). It is a strong X-ray AGN observed by Swift (Maselli et al. 2010).

Stationary component C3
The stationary component C3 is detected at 15 GHz between 1995 and 2011 (Lister et al. 2019). The mean position of C3 is X C3 = −0.288 ± 0.055 mas and Y C3 = −0.477 ± 0.111 mas (see Fig. 1), and the distance between cores C0 and C3 is R C3−C0 = 0.557 ± 0.124 mas (see Fig. 4). The flux density of C3 at 15 GHz is variable, but weaker than or equal to the flux density of the core (see Fig. 10). We show that component C3 is responsible for the ejection of component C7.

Stationary component C9
The stationary component C9 has been identified and observed between 2009 and 2011 (Lister et al. 2019). However, the identification of the first points of other components indicates that the first two points of component C10, the first three points of component C15, and component C16, which are located at the same position as component C9, are probably associated with component C9 (see Figs. 4 and 5). This means that component C9 is detected between 2009 and 2017.
The mean position of C9 is X C9 = −0.082 ± 0.025 mas and Y C9 = −0.284 ± 0.061 mas (see Figs. 1 and 5) and the distance between cores C0 and C9 should be R C0−C9 = 0.296±0.066 mas. However, the error on the distance between cores C0 and C9 is not symmetrical because of the condition R disk,C0 < D L1−C0 (see Eq. (3) and Appendix B.2.1), and we adopt R C0−C9 = 0.296 +0.066 −0.020 mas (see Fig. 4) for the distance between the core and C9. A better determination of the distance between C0 and C9 and its error bars could be achieved during the fit of component C6 (see Appendix B.2.1), but this is beyond the scope of this article. Its flux density at 15 GHz is variable but often higher than the flux density of the core (see Fig. 12). It is detected at 15 GHz between 2009 and 2017.
We show that component C13 is ejected from the stationary component C9. Components C10, C15, and C17 are probably also ejected by the BH associated with component C9.

Components used to fit the model
We fit the precession model and the BBH model parameters using components C6, C7, and C13, that is, components that on the one hand, have enough observed points to have a welldefined trajectory, and on the other hand, have observed points in the first two or three mas to precisely define the origin of the component. We therefore did not fit components C1, C2, and C4, which have no observations in the first 3 mas (C1) and in the first 2 mas (C2 and C4). Component C5, which has a poorly defined trajectory, has not been fit either.
The observations at 15 GHz that we used for this study correspond to 37 epochs and are reported in Table 3 of Lister et al. (2019). We used the model fit data of Lister et al. (2019) to fit the coordinates X(t) and Y(t) of components C6, C7, and C13 using the precession and the BBH model. The results of the fits of components C6, C13, and C7 are given in Sects. 3.4, 3.5, and 3.6, respectively. The positions of the stationary components C0, C9, and C3, the positions of the moving components C6, C7, and C13, and the trajectories obtained from the fit of the BBH system are shown in Fig. 7.

Stability condition
We describe an important constraint that can be used to fit an ejected VLBI component below. To obtain a stable solution in a BBH system, for instance, the BBH system C0-C9 (see Fig. 1), the radii of the accretion disks around C0 and C9, that is, R disk,C0 and R disk,C9 , must be smaller than the distances D L1−C0 and D L1−C9 , respectively. Here L1 is the first Lagrange point of the BBH system (the point where the gravitation forces of C0 and C9 are equal).
Calling T b the orbital period and T p the precession period of the accretion disk, we can calculate the mass of the ejecting BH M C0 , T b , and T p for each value of V a the propagation speed of the perturbation along the jet and the beam based on the knowledge of the mass ratio M C0 /M C9 and the ratio T p /T b .
The rotation period of the accretion disk around C0, T disk,C0 , is given by (Britzen et al. 2001) When we assume that the mass of the accretion disk is M disk,C0 M C0 , the radius of the accretion disk R disk,C0 is and we must have the stability condition where D L1−C0 is the distance between C0 and the first Lagrange point L1, and R C0−C9 is the separation of the BBH system C0-C9. The radius of the accretion disk does not depend on V a , which is the propagation speed of the perturbation along the beam and the jet.

Separation of the BBH systems
At the VLBI core, the VLBI jet becomes transparent to a given synchrotron frequency. The position of the VLBI core does not correspond to the position of the supermassive BH (SMBH), and the distance between the VLBI core and the SMBH depends on the observational frequency. At higher frequency, the VLBI core is closer to the SMBH; see, for instance, Marscher et al. (2008). Although the exact dependence of the distance of the VLBI core to the SMBH with the frequency of observation is not known, the distance between the two VLBI cores that are detected at the same frequency provides a good estimate of the distance between the two SMBH when the ejection directions of Fig. 6. Top: the distance between the two VLBI cores that are detected at the same frequency provides a good estimate of the distance between the two SMBH, BH1 and BH2, if ejection directions of the VLBI jets are equal to within few degrees. Bottom: when the jets are ejected in opposite directions or in very different directions, then the distance between the two VLBI cores is very different than the distance between the two SMBH.
the VLBI jets are equal to within few degrees (see Fig. 6). This effect will produce an additional error on the distance determination of the two BH. If the distance between the core and the BH is about 25 µas and the angle between the two jets is 3 • , the additional error is ≈2 × 25 × sin(3 • ) = 2.6 µas. When the jets are ejected in opposite directions or in very different directions, the distance between the two VLBI cores is very different than the distance between the two SMBH (see Fig. 6).

Parameter ranges we explored for the fit
In this section, we provide the parameter ranges we explored to fit the VLBI components C6, C7, and C13. In order to cover a wide range of possible inclination angles, we explored the range 2 • ≤ i o ≤ 40 • . For the mass ratio of the BBH system, we explored the range 10 −7 ≤ M 1 /M 2 ≤ 10, where M 1 is the mass of the BH ejecting the VLBI component. For the ratio T p /T b , we explored the range 1 ≤ T p /T b ≤ 10 4 . For the propagation speed of the perturbation, we explored the range 0.001 c ≤ V a ≤ 0.45 c (we limit ourselves to nonrelativistic hydrodynamics in this model).

Introduction
We studied the three different models given in Table 2. We calculated the curve χ 2 (i o ), which corresponds to the χ 2 of the model of the VLBI component ejection for a given inclination angle i o . The details of the calculations corresponding to the three models are given in Appendix B. In this section, we present the results of fitting the ejection of component C6 from the origin C0 by the BBH system C0-C9.
We found that component C6 cannot be ejected by a single spinning BH but must be ejected by a BBH system.

BBH system C0-C9 : Origin C0
The BBH system is constituted by core C0 and the stationary component C9. Its separation is R C0−C9 = 0.3 mas. Component C6 is ejected by the VLBI core, that is, component C0 (see Fig. 7).

Model
Origin Concave: Solution rejected Table 3. Parameter ranges for the BBH system that ejects C6.
≈3.1 × 10 10 M stable solutions. Moreover, the solution with M C0 /M C9 ≈ 0.5 has the smallest χ 2 min . Thus, component C6 is ejected from core C0 of the BBH system C0-C9, which is characterized by The corresponding solution is characterized by χ 2 min ≈ 36.9, the inclination angle is i o ≈ 6.7 • , the ratio T p /T b is ≈5.5, and the angle between the accretion disk and the rotation plane of the BBH system is Ω(C0) ≈ 3.1 • . The bulk Lorentz factor of the VLBI component is γ c ≈ 9.6, and the mean apparent speed Using the parameters of the solution, we gradually varied V a between 0.001 c and 0.45 c. The function χ 2 (V a ) remained constant, indicating a degeneracy of the solution (see Fig. B.6). We deduced the variation range of the C0-C9 parameters in the BBH system that eject C6. They are given in Table 3.
Based on the knowledge of the mass ratio M C0 /M C9 and the ratio T p /T b , we calculated the mass of the ejecting BH M C0 , the orbital period T b , and the precession period T p for each value of V a . We found that the radius of the accretion disk around C0 does not depend on V a and is R disk,C0 = 0.106 mas = 0.46 pc.

Introduction
We studied the four different cases given in Table 4. We calculated the curve χ 2 (i o ), which corresponds to the χ 2 of the model of the VLBI component ejection for a given inclination angle i o . Details of the calculations corresponding to the four models are given in Appendix C. In this section, we present the results of fitting the ejection of component C13 from the origin C9 by the BBH system C0-C9.

BBH system C0-C9 : Origin C9
The BBH system is constituted by core C0 and the stationary component C9. Its separation is R C0−C9 = 0.3 mas. Component C13 is ejected by the stationary component C9 (see Fig. 7).
We found that the mass ratio M C9 /M C0 is M C9 /M C0 > 1 (see Fig. C.5). Using the stability condition, R disk,C9 < D(L1−C9), we found that solutions with the mass ratio M C9 /M C0 ≥ 3 have R disk,C9 ≥ D(L1−C9) and do not correspond to stable solutions. Moreover, the solution with M C9 /M C0 = 2 has the smallest χ 2 min . Thus, component C13 is ejected from the stationary component C9 of the BBH system C0-C9, which is  The corresponding solution is characterized by χ 2 min ≈ 9.3, and the inclination angle is i o ≈ 5.8 • . The ratio T p /T b is ≈10.2, and the angle between the accretion disk and the rotation plane of the BBH system is Ω(C9) ≈ 3.3 • . The bulk Lorentz factor of the VLBI component is γ c ≈ 7.3, and the mean apparent speed Using the parameters of the solution, we gradually varied V a between 0.001 c and 0.45 c. The function χ 2 (V a ) remained constant, indicating a degeneracy of the solution (see Fig. C.8). We deduced the range of variation of the BBH system parameters ejecting C13. They are given in Table 5.
We found that component C13 cannot be ejected by a single spinning BH but must be ejected by a BBH system.
We found that the radius of the accretion disk around C9 does not depend on V a and is R disk,C9 = 0.140 mas = 0.61 pc. The ratio

Introduction
We studied the three different models given in Table 6. We calculated the curve χ 2 (i o ), which corresponds to the χ 2 of the model of the VLBI component ejection for a given inclination angle i o . The details of the calculations corresponding to the three models are given in Appendix D. In this section, we present the results of fitting the ejection of component C7 from the origin C3 by the BBH system (C0+C9)-C3.
We found that component C7 cannot be ejected by a single spinning BH but must be ejected by a BBH system.
We found that the mass ratio Fig. D.1). We found that solutions A101, page 5 of 16 A&A 634, A101 (2020) Table 6. Different models investigated for the fit of C7.

Model
Origin The corresponding solution is characterized by χ 2 min ≈ 38.4, and the inclination angle is i o ≈ 6.9 • . The ratio T p /T b is ≈10.7, and the angle between the accretion disk and the rotation plane of the BBH system is Ω(C0) ≈ 3.8 • . The bulk Lorentz factor of the VLBI component is γ c ≈ 7.1, and the mean apparent speed is V ap (C7) ≈ 6 c. The ejection time of the VLBI component is Using the parameters of the solution, we gradually varied V a between 0.001 c and 0.45 c. The function χ 2 (V a ) remained constant, indicating a degeneracy of the solution (see Fig. D.3). We deduced the parameter range of variation of the BBH system that ejects C7. They are given in Table 7.
We found that the radius of the accretion disk around C3 does not depend on V a and is R disk,C3 = 0.020 mas = 0.09 pc. The ratio

BH system characteristics in 2201+315
The inclination angle of 2201+315 is 5.8 • ≤ i o ≤ 6.9 • . The kiloparsec-scale radio map (Cooper et al. 2007) shows two extended lobes, and the total extension of the source is about 75 arcsec. When we assume i o ≈ 6.5 • , the total extension of the extended lobes is about 2.9 Mpc. The structure of the nucleus is shown in Fig. 1. Table 7. Parameter ranges for the BBH system that ejects C7.
To obtain the final range of the parameters of the BH system C9-C0-C3, we intersected the ranges of the BBH system parameters we found after the fits of C13, C6, and C7. The characteristics of the BH system C9-C0-C3 are given in Table 8 where Ω is the angle between the accretion disk and the rotation plane of the BHs.
The BH system and the component trajectories of C6, C7, and C13 using the BBH model we obtained in Sects. 3.4, 3.6, and 3.5 are shown in Fig. 7.

Discussion and conclusion
The analysis of MOJAVE model fit data with our minimization method suggests that none of components C6, C7, and C13 can be ejected by a single spinning BH, but that they are ejected by BBH systems and that the nucleus of 2201+315 contains three BHs. The three BHs are associated with the VLBI core C0 and the stationary components C9 and C3 as defined by Lister et al. (2019). This source constitutes another example of a nucleus that contains several BHs that are detected in radio and eject VLBI components.
Stationary VLBI components are frequently observed (Jorstad et al. 2017) and are generally assumed to be associated with stationary or recollimation shocks (see Mizuno et al. 2015, Martí et al. 2016, Hada et al. 2018 and references therein), but in the two-fluid model when the relativistic e − −e + beam dissipates into the subrelativistic e − −p + jet, a VLBI component is observed that moves subrelativistically. This appears as a quasi stationary VLBI component (see, e.g., component 1 of 1532+016; Lister et al. 2019), and finally the stationary components can also be associated with BHs that eject VLBI components. The stationary components are not necessarily all associated with BHs or with recollimation shocks. To explain the MOJAVE observations of 2201+315, we do not need stationary or recollimation shocks.
As indicated in Roland et al. (2015), a BBH system can produce three perturbations of the VLBI ejection by the precession A101, page 6 of 16 Fig. 8. Positions of components J1, J2, and J3 detected at 43 and 86 GHz by Cheng et al. (2018). The VLBI components C6, C13, and C7 are shown.
of the accretion disk, the motion of the two BHs around the center of gravity of the BBH system, and the possible motion of the BBH system around either a third BH or another BBH system. This third perturbation produces a change in the VLBI jet direction of between 5 to 15 mas. It is observed for instance in the cases of 1928+738 (Roland et al. 2015 When we observe a change in VLBI jet direction between 5 to 15 mas, for example, the slow rotation of the BBH system, has to be modeled and the VLBI data have then to be corrected for this slow rotation, so that finally the ejection of the corrected data can be modeled using a BBH system. This has been done by Roland et al. (2015) in the case of 1928+738, which contains two BBH systems. In the case of 2201+315, the mass of the third BH is far lower than that of the two others, and the ejection of components C6 and C13 can be modeled using the BBH system C0-C9 and the ejection of component C7 using the BBH system (C0+C9)-C3.
As indicated in the introduction, when the nucleus contains a BBH system, the two BHs can eject VLBI components and both can be detected using VLBI observations, but it is not necessarily always the case. Roland et al. (2015) showed that 1928+738 contains two BBH systems: the two BHs of the first BBH system are detected (associated with a stationary component) and eject VLBI components; only one BH of the second BBH system ejects VLBI components during the observations; no BH of the second BBH system is associated with a stationary component.
The source 2201+315 has been observed between 1995 and 2018 in the MOJAVE survey. While core C0 is detected during the whole period, component C3 is detected between 1995 and 2010 and component C9 is detected between 2009 and 2017. Because the stationary components C3 and C9 correspond the place in the VLBI jet that becomes transparent to a given synchrotron frequency (see Sect. 3.2), the disappearance of component C3 in 2013 means that at least the ejection of the relativistic e − e + plasma (the beam) stopped in 2010; this does not mean that the ejection of the subrelativistic e − p plasma (the jet) also stopped in 2010. Component C9 started to be detected in 2009, which means that the ejection of the relativistic e − e + plasma (the beam) started in 2010.
Source 2201+315 has been observed at 43 and 86 GHz by Cheng et al. (2018). They detected three components J1, J2, and J3. Their positions are shown in Fig. 8. Component J1 probably belongs to the family of trajectories defined by C13 and has been ejected from the stationary component C9. It is just as probable that component J3 belongs to the family of trajectories defined by C6 and has been ejected from core C0 or is ejected from the stationary component C3. For component J2, the situation is less clear because of the shift that is observed with the family of trajectories defined by C13, but it can probably be associated with the family of trajectories defined by C13 and has been ejected from the stationary component C9.
An interesting corollary of our study resides in the importance of detecting and separating possible stationary components and determining their flux density ratios, and thereby explaining the position variations measured during geodetic VLBI sessions. Radio source 2201+315 has been intensively observed in the framework of the permanent geodetic VLBI program devoted to absolute astrometry (Fey et al. 2015). In geodetic VLBI analysis, the coordinates of the "radio center" can be estimated on average over observing sessions of 24 h each, and session-wise coordinate time series can be obtained in this way. The radio center is not located at the brightest component of the radio source but can rather be seen as a barycenter of the brightness distribution mitigated by the effect of the network geometry. The accuracy of coordinates as determined by geodetic VLBI is typically a few dozen microarcseconds. This allows recording signatures of subtle changes in the source structure as small variations in the coordinates. Radio source 2201+315 was observed in more than 1000 sessions since 1979, mainly between 1990 and 2010. In Fig. 9 each dot represents the position of the radio center determined in one session of 24 h (see Gattano et al. 2018, for more details). During some periods, dramatic changes in the coordinate time series can be observed. These variations can be explained by two types of variability that occur in the nucleus, that is, the ejection of a new VLBI component whose flux density can be higher than the flux densities of the two stationary components, and in the case of a BBH, the variations in the ratio of the flux densities of the VLBI components associated with two BHs.
It is generally observed that changes in astrometric position correspond to the direction of the VLBI jet (Moór et al. 2011). However, if the direction of the BBH system is different from the direction of the VLBI jet, these two types of variations will produce changes in astrometric position in two different directions. During 1997During -2005 of the core, S ν (C0), is higher than the flux density of component C3, S ν (C3). However, the flux densities of the core and C3 can change rapidly and significantly, and consequently, there are periods where the ratio S ν (C0)/S ν (C3) can be S ν (C0)/S ν (C3) ≈ 1. In 1998-1999, the flux density of the core was similar to the flux density of C3, and the flux density of the nucleus was dominated by the flux density of component C5, which was close to C3 (see Figs. 10 and 11). Thus the mean astrometric position deduced at 8 GHz is shifted to the south. Between 2001 and 2002, the flux density of the nucleus was dominated by the flux density of the core, C0, which increased with time, and the mean astrometric position deduced at 8 GHz is shifted to the north. However, in 2003, the flux density of the core, C0, was similar to the flux density of C3 and the flux density of the nucleus was dominated by the flux density of component C6, which was close to C3 (see Figs. 10 and 11). Thus the mean position deduced at 8 GHz is shifted to the south and we observe a quick change in the astrometric position time series. During the periods 2009 to 2012 and 2014 to 2015, the flux density of component C9 was higher than the flux density of the core (see Fig. 12). During these periods, the astrometric position is shifted along the direction C0-C9, which is different from the VLBI jet direction. Thus the astrometric position of 2201+315 moves with time in two different directions. The mean astrometric position of 2201+315 is located between the stationary components C0, C9, and C3 and does not correspond to any of their positions.
This result shows the importance of determining the structure of nuclei of extragalactic radio sources and especially of quasars using VLBI observations at frequencies ≥15 GHz in order to detect and separate possible stationary components and to determine their flux density ratio. This will allow explaining the time shifts of the position observed during geodetic VLBI observations and Gaia observations. This also shows the importance of large MOJAVE-like VLBI databases for improving the understanding of the apparent displacement of the radio center measured by geodetic VLBI and thereby in the maintenance of the celestial reference frames. Before fitting components C6, C7, and C13, we recall the geometry, parameters, and basic equations of the model, which have been developed in Britzen et al. (2001), Lobanov & Roland (2005), and Roland et al. (2008Roland et al. ( , 2013. Details can be found in Roland et al. (2008Roland et al. ( , 2013.

A.1. Geometry of the model
We call Ω the angle between the accretion disk and the orbital plane (XOY) of the BBH system. The component is ejected in a cone (the precession cone) with its axis in the Z OZ plane and its opening angle is Ω. We assumed that the line of sight is in the plane (YOZ) and that it forms an angle i o with the axis Z OZ (see Fig. A.1). The axis η corresponds to the mean ejection direction of the VLBI component projected on a plane perpendicular to the line of sight, so that the plane perpendicular to the line of sight is the plane (ηOX). We call ∆Ξ the rotation angle in the plane perpendicular to the line of sight to transform the coordinates η and X into coordinates N (north) and W (west), which are directly comparable with the VLBI observations. We have The sign of the coordinate W was changed from Roland et al. (2008) to use the same definition as in the VLBI observations.

A.2. General perturbation of the VLBI ejection
For VLBI observations, the origin of the coordinates is BH 1, that is, the BH that ejects the VLBI components. For the sake of simplicity, we assumed that the two BHs have circular orbits, that is, e = 0. Therefore, the coordinates of the moving components in the frame of reference where BH 1 is considered the ejection origin are (Roland et al. 2008) where R o (z) is the amplitude of the precession perturbation, given by R o (z) = R o z c (t)/(a + z c (t)), with a = R o /(2 tanΩ); ω p is ω p = 2π/T p , where T p is the precession period, and k p is defined by k p = 2π/T p V a , where V a is the speed of the propagation of the perturbations; ω b is ω b = 2π/T b , where T b is the BBH system period and k b is defined by k b = 2π/T b V a ; T d is the characteristic time of the damping of the perturbation, and x 1 and y 1 are given by We define with R bin the distance between the two BHs as the separation of the BBH system. It is (A.7) The differential equation governing the evolution of z c (t) can be obtained by defining the speed of the component, where v c is related to the bulk Lorentz factor by v c /c = The coefficients A, B, and C are calculated in Appendix A of Roland et al. (2008). Equation (A.9) admits two solutions corresponding to the jet and the counter-jet.

A.3. Coordinates of the VLBI component
Solving Eq. (A.9), we determine the coordinate z c (t) of a pointsource component that is ejected relativistically in the perturbed beam. Then, using Eqs. (A.3) and (A.4), we can find the coordinates x c (t) and y c (t) of the component. In addition, for each point of the trajectory, we can calculate the derivatives dx c /dt, dy c /dt, and dz c /dt and then deduce cos θ, δ c , S ν , and t obs (see Roland et al. 2013).
After calculating the coordinates x c (t), y c (t), and z c (t), they can be transformed into the w c (t) (west) and n c (t) (north) coordinates using Eqs. (A.1) and (A.2).
As explained in Britzen et al. (2001), Lobanov & Roland (2005), and Roland et al. (2008), the radio VLBI component has to be described as an extended component along the beam. We call n rad the number of points (or integration steps along the beam) for which we integrate to model the component. The coordinates W c (t) and N c (t) of the VLBI component are then .11) and can be compared with the observed coordinates of the VLBI component, which correspond to the radio peak intensity coordinates provided by model-fitting during the VLBI data reduction process.

A.4. Model parameters
In this section, we list the possible free model parameters. The parameter n rad is known when the size of the VLBI component is known. In this article we fit components C6, C7, and C13 assuming that their projected size on the sky is = 0.1 mas. This means that practically, the problem we have to solve is a 15 free-parameter problem.
We have to investigate the different possible scenarios with regard to the sense of the rotation of the accretion disk and the sense of the orbital rotation of the BBH system. These possibilities correspond to ±ω p (t − z/V a ) and ±ω b (t − z/V a ). Because the sense of the precession is always opposite to the sense of the orbital motion (Katz 1997), we study the two cases denoted by +− and −+, where we have ω p (t − z/V a ), −ω b (t − z/V a ) and −ω p (t − z/V a ), ω b (t − z/V a ), respectively.

A.5. Method for solving the problem
This method is a practical one that provides solutions, but the method is not unique and does not guarantee that all possible solutions are found. We calculate the projected trajectory on the plane of the sky of an ejected component and determine the parameters of the model to simultaneously produce the best fit with the observed west and north coordinates. The parameters we found minimize where χ 2 (W c (t)) and χ 2 (N c (t)) are the χ 2 calculated by comparing the VLBI observations with the calculated coordinates W c (t) and N c (t) of the component. For instance, to determine the inclination angle that provides the best fit, we minimize χ 2 t (i o ). The concave parts of the surface χ 2 (i o ) contain a minimum. We can find solutions without a minimum; they correspond to the convex parts of the surface χ 2 (i o ) and are called mirage solutions.
To illustrate the properties of the surface χ 2 (i o ), we plot in Fig. A.2 a possible example of a profile of the solution χ 2 (i o ). Figure A.2 shows two possible solutions for which χ 2 (Sol1) ≈ χ 2 (Sol2). Solution 2 is more robust than solution 1, that is, it is the deepest one, and we adopted this solution. There are two possible solutions for which χ 2 (Sol1) ≈ χ 2 (Sol2). They correspond to the concave parts of the surface χ 2 (i o ). However, solution 2 is more robust than solution 1, i.e., it is the deepest one, and we adopted this solution.
We define the robustness of the solution as the square root of the difference between the smallest maximum close to the minimum and the minimum of the function χ 2 . A solution of robustness 3 is a 3 σ solution, that is, 3 σ ⇔ ∆χ 2 = 9.
The main difficulties we have to solve are the following: to find all possible solutions, to eliminate the mirage solutions, to find the most robust solutions.
One parameter allows us to determine the possible solutions. This fundamental parameter is the ratio T p /T b , where T p and T b are the precession period of the accretion disk and the binary period of the BBH system, respectively.
To find the possible solutions, for a given inclination angle and for several values of T p /T b (say T p /T b = 10, 30, 100. . . ), we calculate the function χ 2 (M 1 /M 2 ), for 10 −7 ≤ M 1 /M 2 ≤ 100 in a first step, and then for a given value of M 1 /M 2 determined previously, we calculate χ 2 (T p /T b ) for 1 ≤ T p /T b ≤ 10 000. The determination of M 1 /M 2 does not depend on the choice of the inclination angle. Because we investigate a wide range for the parameters M 1 /M 2 and T p /T b , we expect to be able to find all possible solutions.
In a second step, for different values of M 1 /M 2 found previously, assuming that T p /T b is a free parameter, we calculate χ 2 (i o ). When the solution is found, it is not unique, but is a family of solutions. The solution shows a degeneracy, and we show below that the parameter for studying the degeneracy, that is, to find the range of parameters that provide the family of solutions, is V a , the propagation speed of the perturbation along the beam.
Generally, for any value of the parameters, the surface χ 2 (λ), where λ can be any of the free parameters of the system, i o , M 1 /M 2 , or T p /T b , for instance. is convex and does not present a minimum. Moreover, when we are in the convex part of the surface χ 2 (λ), one of the important parameters of the problem can diverge. The important parameters of the problem that can diverge are: the bulk Lorentz factor of the e ± beam, which has to be γ c ≤ 30. This limit is imposed by the stability criterion for the propagation of the relativistic beam in the subrelativistic e − −p jet; the radius of the accretion disk, which has to be smaller than the distance between the distance between the BH and the first Lagrange point (see Sect. 3.1); the ratio T p /T b , which has to be >1, when T p /T b = 1 the radius of the accretion disk becomes equal to the separation of the BBH system; and the total mass of the BBH system.
We studied the three different models given in Table 2. We calculated the curves χ 2 (i o ), which corresponds to the χ 2 of the model of the ejection of the VLBI component for a given inclination angle i o .

B.1. Precession model: Origin C0
In this section we model the ejection of the VLBI component C6 assuming a single spinning BH, that the VLBI ejection is perturbed by the precession of the accretion disk, and that the ejection origin is core C0. The direction for the precession rotation is chosen to be +.
We found that the possible range for the inclination angle is limited to 2 • ≤ i o ≤ 11.6 • . The bulk Lorentz factor γ c exceeds 30 and diverges when i o > 11.6 • . The curve γ c (i o ) is shown in Fig. B.1.

B.2. BBH C0-C9 model: Origin C0
In this section we model the ejection of the VLBI component C6 assuming a BBH system that consists of core C0 and the stationary component C9. Its separation is R C0−C9 = 0.3 mas, the VLBI ejection is perturbed by the precession of the accretion disk and the motion of the BHs, and the ejection origin is core C0. The directions for the precession rotation and the BBH rotation are chosen to be + and −, respectively.
Following the method given in Appendix A.5, we calculated the function χ 2 (M C0 /M C9 ) for i o = 8 • and T p /T b = 10. We found that the mass ratio is 0.1 ≤ M C0 /M C9 < 1 (see Fig. B.3), that is, component C6 ejected from core C0 is ejected by the smallest BH of the BBH system C0-C9. It is possible to constrain and determine the mass ratio M C0 /M C9 when the radius of the accretion disk around C0, R disk,C0 is smaller than the distance D L1−C0 to obtain a stable solution (see Sect. 3.1 and Eqs. (2) and (3)).
Then we calculated χ 2 (T p /T b ) and χ 2 (i o ) for various values of M C0 /M C9 and found that solutions with M C0 /M C9 ≥ 0.6 have R disk,C0 ≥ D L1−C0 and do not correspond to stable solutions. Moreover, we found that for solutions with M C0 /M C9 ≤ 0.5 the solution with M C0 /M C9 ≈ 0.5 has the smallest χ 2 min , that is, the solution that provides the best fit of the VLBI component C6 corresponds to the BBH system C0-C9 characterized by R C0−C9 = 0.3 mas and M C0 /M C9 = 0.5. The corresponding solution is characterized by i o ≈ 6.7 • , T p /T b ≈ 5.5 and Ω ≈ 3. We found that the curve χ 2 (i o ) is convex. It does not have a minimum, that is, there is no stable solution. The curve χ 2 (i o ) is shown in Fig. B.2.

B.2.1. Separation of the BBH system C0-C9
We showed in Sect. 2.4 that from the mean position of C9, that is, X C9 = −0.082 ± 0.025 mas and Y C9 = −0.284 ± 0.061 mas, the distance between the core C0 and C9 should be R C0−C9 = 0.296 ± 0.066 mas. However, the solution corresponding to R C0−C9 = 0.300 mas is characterized by R disk,C0 /D L1−C0 = 0.85, and when we reduce R C0−C9 , we will have R disk,C0 /D L1−C0 > 1. We found that we must have R C0−C9 ≥ 0.280 mas, thus the distance between C0 and C9 is R C0−C9 = 0.296 +0.066 −0.020 mas, that is, the errors are asymmetrical.  When the solution is known, we can determine the best distance between C0 and C9 and its error. We can choose 30 points, for example, whose positions are close to the C9 position within the errors ∆X = 0.025 mas and ∆Y = 0.061 mas. For each point, we calculate its new distance to C0 and the new offset of the VLBI data, and then we calculate the function χ 2 (i o ). Then we will be able to deduce the best value of the distance between the two BHs and its error. However, this is beyond the scope of the article.

B.2.2. Solution family
Using the parameters of the solution, we gradually varied V a between 0.001 c and 0.45 c. The function χ 2 (V a ) remained constant, indicating a degeneracy of the solution (see Fig. B.6). The corresponding ranges of the parameters of the BBH system are given in Sect. 3.4.

B.3. BBH model C0-C3: Origin C0
In this section we model the ejection of the VLBI component C6 assuming the BBH system consists of core C0 and the stationary component C3. Its separation is R C0−C3 = 0.56 mas, the VLBI ejection is perturbed by the precession of the accretion disk and the motion of the BHs, and the ejection origin is C0. The directions for the precession rotation and the BBH rotation are chosen to be + and −, respectively. We searched for the solution with R C0−C3 = 0.56 mas and M C0 /M C3 = 0.5 to compare it with the solution found previously in Appendix B.2. The curve χ 2 (i o ) is shown in Fig. B.7. It . The curve χ 2 (i o ) has a a minimum, but the corresponding solution has a χ 2 min larger than the χ 2 min of the solution of the fit of C6 assuming it is ejected by the BBH system C0-C9 characterized by R C0−C9 ≈ 0.3 mas and M C0 /M C9 = 0.5 (see Fig. B.4). Moreover, it is a very weak solution, i.e., its robustness is ≤ 1 (see Appendix A.5) and is rejected. has a minimum, but the corresponding solution has a χ 2 min ≈ 39.1 larger than the χ 2 min ≈ 36.9 of the solution of the fit of C6 assuming it is ejected by the BBH system C0-C9 characterized by R C0−C9 = 0.3 mas and M C0 /M C9 = 0.5 (see Fig. B.4). Moreover, it is a very weak solution, that is, its robustness is ≤ 1 (see Appendix A.5) and is rejected.  We found that the curve χ 2 (i o ) is convex, it does not have a minimum, that is, there is no stable solution. The curve χ 2 (i o ) is shown in Fig. C.2.

C.2. Precession model: origin C9
In this section we mode the ejection of the VLBI component C13 assuming a single spinning BH, the VLBI ejection is perturbed by the precession of the accretion disk, and the ejection origin is the stationary component C9. The direction for the precession rotation is chosen to be +.
We found that the possible range for the inclination angle is limited to 2 • ≤ i o ≤ 14.7 • . The bulk Lorentz factor γ c exceeds 30 and diverges when i o > 14.7 • . The curve γ c (i o ) is shown in Fig. C.3.
We found that the curve χ 2 (i o ) is convex, it does not have a minimum, that is, there is no stable solution. The curve χ 2 (i o ) is shown in Fig. C.4.

C.3. BBH model C0-C9: Origin C9
In this section we model the ejection of the VLBI component C13 assuming that the ejection origin is the stationary component C9 and that it is ejected by the same BBH system that has ejected C6, that is, the BBH system C0-C9, whose separation is R C0−C9 = 0.3 mas. The VLBI ejection is perturbed by the precession of the accretion disk and the motion of the BHs, and the ejection origin is core C0. The directions for the preces- sion rotation and the BBH rotation are chosen to be + and −, respectively. First, we calculated the function χ 2 (M C0 /M C9 ) for i o = 8 • and T p /T b = 10. We found that the mass ratio is M C9 /M C0 ≥ 1.0 (see Fig. C.5), that is, component C13 ejected from the stationary component C9 is ejected by the largest BH of the BBH system. It is possible to constrain and determine the mass ratio M C9 /M C0 when the radius of the accretion disk around C9, R disk,C9 is smaller than the distance D L1−C9 to obtain a stable solution (see Sect. 3.1).
Then we calculated χ 2 (T p /T b ) and χ 2 (i o ) for various values of the mass ratio M C9 /M C0 , that is, M C9 /M C0 = 1, 2, and 3. The curves corresponding to χ 2 (i o ) are shown in Fig. C.6. We found that solutions with M C9 /M C0 ≥ 3.0 have R disk,C9 ≥ D L1−C0 and do not correspond to stable solutions. The curves (R disk,C9 /D L1−C9 )(i o ) corresponding to the mass ratios M C9 /M C0 = 1, 2, and 3 are shown in Fig. C.7. Moreover, we found that for solutions with M C9 /M C0 ≤ 2.0 the solution with M C9 /M C0 = 2.0 has the smallest χ 2 min , that is, the solution that provides the best fit of the VLBI component C13 corresponds to the BBH system C0-C9 with R C0−C9 = 0.3 mas and M C9 /M C0 = 2.0 ± 0.5 assuming that component C13 is ejected from the stationary component C9. The solution is characterized by by the inclination angle i o ≈ 5.8 • , the ratio T p /T b ≈ 10.2, and the angle Ω ≈ 3.3 • .
Using the parameters of the solution, we gradually varied V a between 0.001 c and 0.45 c. The function χ 2 (V a ) remained constant, indicating a degeneracy of the solution (see Fig. C.8). The corresponding ranges of the parameters of the BBH system are given in Sect. 3.5.

C.4. BBH model C0-C9: Origin C0
In this section we model the ejection of the VLBI component C13 assuming that the ejection origin is core C0 and that it is ejected by the same BBH system that has ejected C6, that is, the BBH system C0-C9, whose separation is R C0−C9 = 0.3 mas, and the mass ratio is M C0 /M C9 = 0.5. The VLBI ejection is perturbed by the precession of the accretion disk and the motion of the BHs, and the ejection origin is core C0. The directions for the precession rotation and the BBH rotation are chosen to be + and −, respectively. First, we calculated the function χ 2 (M C0 /M C9 ) for i o = 8 • and T p /T b = 100. We found that the mass ratio is 1.0 ≤ M C0 /M C9 ≤ 10 (see Fig. C.9), that is, component C13 ejected from the core is ejected by the largest BH of the BBH system. This result contradicts the result found for the mass ratio M C0 /M C9 ≈ 0.5. during the fit of component C6 ejected from the core C0 of the same BBH system C0-C9 (see Appendix B.2). This result shows already that component C13 cannot be ejected from core C0 of the BBH system C0-C9.  is no solution, and thus component C13 cannot be ejected from the core C0 of the BBH system C0-C9 characterized by R C0−C9 = 0.3 mas and the mass ratio M C0 /M C9 = 0.5.

Appendix D: Fit of component C7
We studied the three different models given in Table 6. We calculated the curves χ 2 (i o ).

D.1. Precession model: Origin C0
In this section we model the ejection of the VLBI component C7 assuming a single spinning BH, that the VLBI ejection is perturbed by the precession of the accretion disk, and that the ejection origin is core C0. The direction for the precession rotation is chosen to be +. We found that the possible range for the inclination angle is limited to 2 • ≤ i o ≤ 14.4 • . The bulk Lorentz factor γ c exceeds 30 and diverges when i o > 14.4 • .
We found that the curve χ 2 (i o ) is convex, it does not have a minimum, that is, there is no stable solution. The curve χ 2 (i o ) is shown in Fig. D.2.

D.2. BBH model (C0+C9)-C3: Origin C3
In this section we model the ejection of the VLBI component C7 assuming that it is ejected by the BBH system formed by C0+C9 and C3, whose separation is R G−C3 = 0.370 mas (see Fig. 1). The VLBI ejection is perturbed by the precession of the accretion disk and the motion of the BHs, and the ejection origin is component C3. The directions for the precession rotation and the BBH rotation are chosen to be + and −, respectively.  have a minimum, however, and when i o becomes smaller than 4.8 • , the radius of the accretion disk diverges and becomes larger than the separation of the BBH system (see Fig. D.2). We adopted the solution with M C3 /(M C0 + M C9 ) = 0.01 and R G−C3 = 0.37 mas, which corresponds to the best solution. The corresponding curve χ 2 (i o ) is shown in Fig. D.2. The solution is characterized by an inclination angle i o ≈ 6.9 • , T p /T b ≈ 10.7, and Ω ≈ 3.8 • . The value at the minimum χ 2 min ≈ 38.4 is lower than the value of the minimum of χ 2 (i o ) when component C7 is ejected by the BBH system C0-C9 from the core C0 (see Fig. D.4).
Using the parameters of the solution, we gradually varied V a between 0.001 c and 0.45 c. The function χ 2 (V a ) remained constant, indicating a degeneracy of the solution (see Fig. D.3). The corresponding ranges of the parameters of the BBH system are given in Sect. 3.6.

D.3. BBH model C0-C9: origin C0
In this section we model the ejection of the VLBI component C7 assuming that it is ejected by the same BBH system that has ejected C6, that is, the BBH system C9-C0, whose separation is R C0−C9 = 0.3 mas, and the mass ratio is M C0 /M C9 = 0.5. The VLBI ejection is perturbed by the precession of the accretion disk and the motion of the BHs, and the ejection origin is core C0. The directions for the precession rotation and the BBH rotation are chosen to be + and −, respectively.
First, we calculated the function χ 2 (M C0 /M C9 ) for i o = 8 • and T p /T b = 10. The solutions corresponding to the ejection of C7 are characterized by M C0 /M C9 < 1.0.
The curve χ 2 (i o ) corresponding to the mass ratio M C0 /M C9 = 0.5 is shown in Fig. D.4. The curve χ 2 (i o ) has a minimum, χ 2 min ≈ 41.8 at i o ≈ 5.0 • . The solution is characterized by T p /T b ≈ 9.3 and Ω ≈ 1.2 • . These values should be the same as those found for the solution corresponding to the ejection of C6 by C0, that is, (T p /T b ) C6 ≈ 5.5 and Ω C6 ≈ 3.1 • . When C0 ejects different components, they are characterized by the same geometrical parameters T p /T b and Ω. This solution is therefore not the correct one and is rejected. Its χ 2 min ≈ 41.8 is larger than the χ 2 min ≈ 38.4 of the solution obtained for the BBH system (C0+C9)-C3 (see Appendix D.2).