On the evolution of a twisted thin accretion disc in eccentric inclined supermassive binary black holes

We propose a model of a twisted accretion disc around a Kerr black hole interacting with a secondary black hole of a smaller mass on an inclined eccentric orbit. We use parameters of the system, which may be appropriate for the so-called 'precessing massive' model of OJ 287. We calculate expressions for torque exerted on the disc by the secondary and a contribution of the secondary to the apsidal precession of disc elements by a double averaging procedure over the periods of the secondary and the disc elements. These expressions are used at all scales of interest, including the ones inside the binary orbit. We calculate numerically the evolution of the disc tilt and twist assuming a flat initial configuration. We consider the disc aspect ratio $h/r=10^{-3}$, a rather large viscosity parameter $\alpha=0.1$ and several values of the primary rotational parameter, $\chi$. We find that, after a few periods of Lense-Thirring precession of the orbit, the disc relaxes to a quasi-stationary configuration in the precessing frame with a non-trivial distribution of the disc inclination angle, $\beta$, over the radial scale. We propose an analytic model for this configuration. We show that the presence of the twisted disc leads to multiple crossings of the disc by the secondary per one orbital period, with time periods between the crossings being different from the flat disc model. Our results should be taken into account in the modelling of OJ 287. They can also be applied to similar sources.


INTRODUCTION
Supermassive binary black holes (SBBHs) forming as a result of the inevitable process of galactic mergers could lead to many striking observational effects, see e.g.Ferrarese & Ford (2005), Komossa (2006) and Dotti, Sesana & Decarli (2012) for general reviews.Typically, candidates to SBBHs are identified as quasars and active galactic nuclei exhibiting either quasi-periodic activity or some peculiar spectral features in the ultraviolet, optical, and infrared bands, see e.g.Graham et al (2015) and Eracleous et al (2012).Some of the SBBH candidates can also be resolved by the present and planned mm-VLBI missions, see Malinovsky & Mikheeva (2023).Among all SBBH candidates, the most prominent place, perhaps, belongs to the famous BL Lac object OJ 287 that demonstrates quasi-periodic activity with a period ∼ 11 − 12yr over the time span of a century, which may be attributed to the presence of a SBBH, Sillanpaa et al (1988).
As was noted, for example, by Lehto & Valtonen (1996), since the outbursts from OJ 287 often come in pairs, with individual outbursts in a pair separated by time periods order of or slightly larger than a year, a natural explanation of this activity would be to associate the outbursts with collisions between a smaller black hole (a perturber) orbiting a larger primary black hole on an eccentric orbit and an accretion disc around the primary.In this model, the orbital plane and the disc plane are misaligned.Since the time between the outbursts in a pair is much shorter than that between two successive pairs, it is natural to assume that the perturber-accretion disc collisions happen near the periastron of the orbit.Implementation of this model has led to rather extreme parameters of the system, namely, the primary mass being as large as 20 billion solar masses and the perturber with mass of the order of 100 million solar masses having quite an eccentric orbit with eccentricity e ≈ 0.7 and semi-major axis a order of a hundred gravitational radii, see e.g Valtonen et al (2008) and Valtonen et al (2023) and references therein for a more recent discussion.
Additionally, as we have mentioned above, the orbit should be significantly inclined with respect to the disc plane, so that its inclination angle well exceeds the disc opening angle.The model is referred hereafter to as PM (precessing massive) model of OJ 287.
The unusual parameters of the system naturally have led to a number of questions regarding validity of the model.For example, it was pointed out a long time ago (e.g.Ivanov, Papaloizou & Polnarev (1999)) that, when a binary of unequal masses interacts with an accretion disc and the mass ratio, q, of the smaller component to the larger one is small enough, the low mass component is rather quickly dragged to the disc plane.This happens when the perturber mass, m, is smaller than a typical disc mass enclosed within the orbit.On the other hand, in the opposite case of the light disc, it A twisted disc in binary black holes 3 is deformed in such a way that it aligns with the orbital plane of the binary.The alignment time is again expected to be rather short.Both processes should operate at scales much larger than those required to explain the activity of OJ287 and, therefore, a system with the value of semi-major axis appropriate for OJ 287 is expected to have orbit and the disc, which are aligned.Moreover, in certain cases, a cavity in the disc inside the orbit is expected to be formed.Therefore, the existence of a system with the required properties seems unlikely from a theoretical point of view, at least, at first glance.
Moreover, the suggested model has also been criticised from the observational point of view.
For example, it was claimed recently in Komossa et al (2023) that an outburst of OJ 287 expected to occur in October 2022 according to the model did not actually happen.Based on this fact and using some empirical relationships, these authors proposed to revise the mass of the primary to a much more modest value ∼ 10 8 M ⊙ .In its turn, the criticism of Komossa et al (2023) has prompted Valtonen et al (2023) to modify the original PM model by slightly changing the arrival times of outbursts in such a way that it is still able to explain the observations.Here we note that in order to model the accretion disc, the authors of the original PM model used systems of noninteracting particles (see Valtonen et al (2023) and references therein), which may be inadequate for the description of hydrodynamical effects in the disc.It was also assumed that the disc is globally flat and resides in a plane, which coincides with the equatorial plane of the primary.
Clearly, it is important to examine these assumptions in order to either validate or refute this model, for example, in favour of the models associated with jet precession (see, e.g.Britzen et al (2023)).Of a special importance would be, for example, to check whether or not the disc has a flat shape or it acquires a more complicated geometrical form due to secular interactions with the perturber.
As we discuss below, that the orbital evolution time is determined by the emission of gravitational waves in PM model alleviates the significance of the theoretical criticism.Namely, it may then be shown that the evolution time of the bound orbit can be significantly shorter than a typical viscous time of the accretion disc provided that the disc is geometrically thin.That means that an extended cavity in the disc cannot be formed before the final merger of the binary and the perturber always has enough material when it intersects the disc to produce hot outflows and the corresponding flaring activity, see e.g.Ivanov, Igumenshchev & Novikov (1998).However, the timescale for alignment of the disc with the orbital plane is expected to be smaller than the viscous time provided that the Shakura-Sunyaev parameter α is small, see Papaloizou & Pringle (1983).
Therefore, the disc may be able , in principle, to have enough time to align with the orbital plane P. B. Ivanov and V. V. Zhuravlev thus changing the basic assumptions of the model relying on perturber-disc collisions as the source of outbursts.This is one of the questions we are going to address in this Paper.
In this work we apply the theory of twisted discs to find qualitative properties of the disc shape on time scales, that are longer than the orbital period at the spatial scales of interest, but shorter than, or of the order of a characteristic time of orbital evolution due to emission of gravitational waves, T GW .We assume that at the scale of interest, the disc is razor thin, with a typical ratio of the disc's half-thickness h to its radius r (the disc opening angle) order of 10 −3 (see e.g.Shakura & Sunyaev (1973) and Ivanov, Igumenshchev & Novikov (1998)) .We employ our model of a fully relativistic non-stationary twisted disc proposed in Zhuravlev & Ivanov (2011) and developed and tested against GRMHD simulations in Morales & Teixeira et al (2014) and Zhuravlev et al (2014) to see what kind of geometrical distortions could be caused by a perturber with the parameters of the system appropriate to PM model and how it can influence predictions of the model concerning the times of outbursts associated with times of crossing of the disc plane by the perturber.In order to bring the problem to a treatable form and be able to use our formalism, we make a number of assumptions.1) We assume that both inclinations of the orbital plane β b and the disc plane β(r, t) are formally small, which allows us to split equations describing dynamical properties of the disc into a 'background part' and 'twisted perturbations'.For the background part, we use the Novikov & Thorne (1973) model and assume that the primary rotational parameter is small enough to treat the background space-time as that of a Schwarzschild black hole.We also neglect potentially important effects associated with the disc's self-gravity.
2) In this Paper we do not make attempts to find a quantitative agreement between our approach and the observations.Basic parameters of the model such as the primary and the perturber masses, orbital semi-major axis and eccentricity are specified according to the values proposed in Valtonen et al (2023).Accordingly, we neglect the evolution of the orbit due to emission of gravitational waves.The disc is initially flat, and we numerically evolve the dynamical equations describing our model up to the end time order of T GW .However, we fully account for effects caused by the Einstein apsidal precession and Lense-Thirring nodal precession of the binary orbit and of the disc rings.
3) Since we are interested in time scales longer than the orbital time scales of the binary and of the disc, we double average terms in the perturbing gravitational potential, which are responsible for apsidal and nodal precession of slightly perturbed disc rings, over the periods of the binary der consideration.Although a similar procedure was frequently used in the case of circumbinary and circumprimary discs in protoplanetary systems (see e.g.Zanazzi & Lai (2018) and Abod et al (2022) for the case of eccentric circumbinary discs) we are not aware of an extensive use of this procedure1 in the entire range of r.Since we use this procedure even in the region between orbital periastron and apoastron, its justification is not relied upon the usual assumption of a large mismatch between the orbital frequencies.Rather, we assume that gas elements of the disc remain on slightly perturbed circular orbits of almost the same inclination over the averaging times since the perturbing potential is proportional to the small mass ratio q.In this way we neglect the potential influence of resonances of various kinds on the evolution of the disc rings2 .The influence of resonances can, in principle, be done separately, see e.g.Artymowicz (1994).

4
) In this Paper we neglect the possibility of the disc's breaking (see e.g.Young et al (2023) and references therein) and the influence of non-linear effects on warp propagation (see e.g.Ogilvie (1999), Ogilvie (2006) and Ogilvie & Latter (2013)).5) Finally, for simplicity, we restrict our attention to the case of rather large viscosity parameter α = 0.1, which is less numerically demanding from the computational point of view.As we are going to show, the use of a large α could lead to rather non-trivial consequences.
It is shown in the Paper, that the disc quickly acquires a non-trivial twisted shape, which is quasi-stationary in the frame precessing with the orbit due to the Lense-Thirring effect.Its typical inclination angle varies rather significantly over the scale of the order of the binary semi-major axis, with a typical variation amplitude comparable to β b .We provide a preliminary numerical study of the conditions of intersections of the perturber with the twisted disc, and we show that both their number per one orbital period and the intersection times differ significantly from the flat disc case.Thus, the effects determined by the non-trivial disc shape should be taken into account in the modelling of the flaring activity of OJ 287 in the framework of PM model (see e.g.Dey et al (2018) and Zwick & Mayer (2023) and the references therein).They can, in principle, either refute the model or lead to a significant modification of the model parameters 3 .The Paper is organised as follows.In Section 2 we introduce basic notations, definitions, and quantities used in the following Sections as well as discuss various physical processes playing a role in our study.In Section 3 we apply our approach to the calculation of the nodal and apsidal precessions.In Section 4 we present our numerical model, discuss the results of computations, and propose a simple analytical model to clarify our numerical results.Finally, in Section 5 we conclude and discuss some possible further developments of this study.

Basic notations
We consider a binary black hole with two unequal masses M and m ≪ M orbiting around each other on a slowly evolving eccentric Keplerian orbit with its semi-major distance being significantly larger than the gravitational radii of both components.The heavier and lighter components are referred hereafter to as the primary component and the perturber, respectively, and their mass ratio q = m M is assumed to be of the order of 10 −2 as suggested by the precessing binary model of OJ 287 (see e.g Valtonen et al (2023) and references therein).It is assumed that there is a thin accretion disc around the primary component and the orbital plane is inclined with respect to the plane of unperturbed disc.

Characteristic temporal and spatial scales
In our problem there are parameters and characteristic scales of two different origins, the ones associated with the binary and the ones associated with the disc.

Characteristic scales associated with the binary
We follow Valtonen et al (2023) for estimates of tentative parameters and associated characteristic scales of the binary.Namely, we assume that the primary's mass, M, could be as large as 2 • 10 10 M ⊙ , while mass of the secondary, m, is about 1.5 • 10 8 M ⊙ .This gives the mass ratio q = 7.5 • 10 −3 .We use M * = M/2 • 10 10 M ⊙ and q * = q/7.5 • 10 −3 .Valtonen et al (2023) assume the orbital period in the observer frame P obs ≈ 12yr, which gives the orbital period in the source frame P orb ≈ 9yr for the redshift z = 0.3 corresponding to OJ 287.Given these values, it is straightforward to evaluate the orbital semi-major axis, a.It is convenient to express it in units of the characteristic gravitational scale r g = GM c 2 and we use hereafter tilde for various spatial scales expressed in units of r g .We have

Characteristic scales and times associated with the disc
We use the standard Shakura & Sunyaev (1973) α-model of the disc.In the framework of this model it can be shown that the disc opening angle δ = h/r, where h is the disc halfthickness is nearly constant and very small, δ ≈ 10 −3 .We use the expression for this angle provided by Ivanov, Igumenshchev & Novikov (1998), hereafter IIN, which is rescaled to take into account that typically we consider larger masses of this primary and larger accretion rate where α * = α/0.1,ṁ = Ṁ /0.1 ṀE , Ṁ is accretion rate and ṀE is its Eddington value defined as in INN. Background quantities of the disc, such as its surface density, accretion rate and δ evolve on a "viscous" time scale where Ω K = GM r 3 is the Keplerian angular frequency of a disc element on a circular orbit of radius r.At scales r ∼ a, T ν is larger than any time scale of our problem.That means that, in general, in regions where the perturber does not strongly influence the disc, it is possible to consider stationary values of the background quantities provided by Shakura & Sunyaev (1973) solution.
Note, however, that this assumption breaks down in the region between the orbital periastron and apoastron, where the perturber intersects the disc plane.
When α > δ the disc tends to align with a symmetry plane of a problem on smaller scales characterised by the corresponding characteristic alignment radius, r al .In our case, the situation is more complicated, since there are two such planes -the black hole equatorial plane and the orbital plane, and, accordingly, two characteristic radii r al,BH and r al,orb .We use the result of Ivanov & Illarionov (1997) for r al,BH , where δ * = δ/10 −3 .For an estimate of r al,orb we use the result of Ivanov, Papaloizou & Polnarev (1999), hereafter IPP, which was obtained under the assumption of a circular orbit and considering radial scales much larger than the semi-major axis.Both these assumptions break down in our case, but, we are going, nonetheless, to use this result for a crude estimate.We have One can see that both characteristic scales can be comparable, and, therefore, one can expect relaxation of the twisted disc to a more non-trivial quasi-stationary configuration than normally expected in the situation when either the Lense-Thirring effect or perturbations of the disc by the binary dominate.
When α > δ the disc tends to relax to a stationary configuration during characteristic time Papaloizou & Pringle (1983).Using (9) we have Comparing ( 7) and ( 12) one can see that when α ∼ 0.01 a and e are changing faster than the relaxation time of the twisted disc.In this situation, on a time scale of the order of T GW effects determined by hydrodynamical interactions in the disc are expected to be less important than the torques on the disc from the side of the binary and black hole.Since this could lead to some new dynamical effects this case is chosen below for our numerical and analytic studies.

A possibility of gap formation
Even when the orbital plane is inclined with respect to the disc plane, there are several physical processes which could lead to the formation of a gap in the disc in the vicinity of the orbit.
Firstly, the perturber intersects the disc twice per orbital period at some radius r int .Each time, the disc material inside the so-called 'accretion radius' measured from the intersection point, is removed (see INN, IPP).This results in the outflow rate of the disc material, Ṁout order of ∼ 2πΣr 2 acc /P orb , where Σ is the disc surface density near the intersection point.Provided that the rate of inflow, Ṁin , of the disc material in the vicinity of the intersection point is smaller than Ṁout a depression of the surface density near r int is formed.As was argued in e.g.INN when Ṁin is associated with the accretion rate for an unperturbed disc even a perturber with a rather small mass can form such a depression.This effect was later observed in numerical simulations, see e.g Xiang-Gruess & Papaloizou (2014) and Hayasaki, Saito & Mineshige (2013) for the inclined case.However, in our case the orbital parameters change with time due to emission of gravitational waves with a characteristic time, which is much smaller than the viscous time T ν .That means that, in our case in the frame comoving with the semi-major axis, the inflow rate can be estimated as ∼ 2πΣr 2 int /T GW , where we take into account that r int ∼ a.In this case, from the condition Ṁout < Ṁin we obtain where we use ( 7).We see that this inequality is marginally satisfied for our chosen value q = 7.5 • 10 −3 and, therefore, the formation of the depression due to the mechanical withdrawal of the disc material is unlikely.
Secondly, when a perturber corotates with the disc, an extended gap can be formed due to the presence of the Lindblad resonances.This effect is expected to be very efficient for a very thin disc and quite eccentric binary, as in the situation we consider, see e.g.Artymowicz & Lubow (1994).
Although the torque transferred to the disc by waves launching at the resonances gets smaller when orbital inclination increases, this effect is expected to be insignificant in a situation when orbital inclination is of the order of or smaller than orbital eccentricity, see e.g.Artymowicz (1994).
However, in our case a characteristic size of the orbit shrinks due to emission of gravitational waves on a time scale, which is much smaller than the viscous time T ν .This allows to suggest that an extended gap is likely to be absent in our case, see e.g.Artymowicz & Lubow (1994) 4 .
We conclude that the gas dynamical effects do not significantly change the surface density distribution.The effects related to excitation of waves in the disc are unlikely to produce an extended gap.On the other hand, we have checked numerically that a possible thin gap practically doesn't change our results.Therefore, we neglect the possibility of gap formation for our analysis below and assume that the basic properties of the disc are given by the standard Novikov & Thorne 1973 solution.

Coordinate systems and associated expressions
We introduce two Cartesian coordinate systems (X h , Y h , Z h ) and (X b , Y b , Z b ) centered at the position of the primary component and associated with the equatorial plane of the central black hole and the orbital plane, respectively.It is implied that the coordinates (X b , Y b , Z b ) are rotated with respect to the coordinates (X h , Y h , Z h ) by Euler angles β b , γ b , Ψ b , where the choice of elementary rotations is the same as in Ivanov & Illarionov (1997) and Demianski & Ivanov (1997).Namely, we adopt the following convention for a choice of rotational angles.First, X and Y axes are rotated by the angle γ b about the Z axis.Then, we perform rotation by the angle β b about the new X axis, and, then rotation about the new Z axis by the angle Ψ b , see Fig. 1 of Ivanov & Illarionov (1997) for a graphical representation.
It is also implied that the periastron of the binary is situated at the X b axis.Thus, the angles β b , γ b and Ψ b define inclination of the orbital plane with respect to the black hole's equatorial plane, a position of the nodal line with respect to X h axis and a position of the apsidal line with respect to the nodal line.
The position vector of the perturber has only two components: D x (t) and D y (t) in the coordinate system associated with the orbit, which are given by the standard Keplerian expressions where a is the semi-major axis, e is the eccentricity and E is the eccentric anomaly, which is related to time t as The effects of General Relativity lead to a secular evolution of the orbital elements.In particular, the Einstein precession and the Lense-Thirring effect lead to the evolution of Ψ b and γ b , respectively, with the corresponding evolution time scales t E and t L−T , following from equations ( 4) and ( 6) above, while emission of gravitational waves causes the evolution of a and e, with the corresponding evolution time scale t GW given by equation ( 7).From these equations it follows that t E ≪ t L−T ≪ t GW .Note that we neglect any contribution to the evolution of orbital elements due to interaction of the binary with the disc.This is justified in the case when a typical disc mass is much smaller than m, as in the case of our system.
Similarly, a disc ring can be described by the 'twisted coordinates' (r, Ψ, ξ) (e.g.Petterson (1978) and also Zhuravlev & Ivanov (2011)) and references therein) with the help of two Euler angles β(r, t) and γ(r, t) that are associated with the black hole equatorial plane.

Concrete expressions for the coordinate transformations
In what follows, we assume that the inclination angles β and β b are small.In this situation one usually employs transformation laws between different coordinate systems in the linear approximation over the inclination angles, see e.g.II.However, in our case, this leads to formally divergent expressions for apsidal and nodal precession frequencies of the disc's rings due to the presence of the perturber, which we calculate below, see the next Section.This happens when r p r r a , where r p and r a are apoastron and periastron of the perturber's orbit.Since we would like to have expressions for these quantities, which are finite everywhere apart from some radius, where the perturber crosses the disc, where such divergences are physically motivated, we have to take into account terms quadratic over the inclination angles when considering transformations between different coordinate systems.
More concretely, the divergences arise due to the fact, that in the linear approximation over β and β b the distance between the perturber and a particular gas element situated at the disc midplane , where X d , Y b , Z b are coordinates of a particular element of the disc corresponding to ξ = 0, is formally the same as in the case when the perturber's orbital plane and the disc midplane coincide.
Clearly, in this situation, ∆r(ξ = 0) is equal to zero at the time when pertuber's distance to the origin of our coordinate systems, D, is the same as r.This happens twice per orbital period, P orb .In order to circumvent this obstacle, we consider the terms quadratic in β and β b in the transformation laws.This results in finite distances at the times of the closest approach apart from a particular value of r = r int belonging to the interval (r p , r a ), which corresponds to the intersection of the orbital plane and the disc midplane.It is assumed below that the disc material must be removed by the action of the perturber's gravitational field close to r in , and, therefore, our expression for the precessional frequencies remain finite at all potentially interesting radial scales.
In order to calculate ∆r 0 and other quantities interesting for us, we express Cartesian coordinates (X b , Y b , Z b ) in terms of the twisted coordinates associated with the equatorial plane, (r, Ψ, ξ) using the rotational matrices defined as in II.It is then convenient to represent , where X 0 , Y 0 do not depend on the inclination angles, X 1 and Y 1 are linear over the inclination angles and X 2 and Y 2 are quadratic over them.
Note that X 0 , Y 0 and X 2 , Y 2 are even with respect to ξ, while X 1 and Y 1 are odd.
We easily obtain the following relations from the expressions provided in Ivanov & Illarionov (1997) where © 2010 RAS, MNRAS 000, 1-37 and Additionally, we have 2.4.2Calculation of a concrete expression for the distance between the perturber and a gas element at the disc's midplane Setting ξ = 0 in ( 23) and substituting the result and ( 18) in ( 17) we obtain where where and the angle Ψ ef f is defined by the conditions It follows from ( 27) and ( 28) that when β ≪ β b we have Finally, we would like to take into account that disc elements with ∆r 0 < r acc , where r acc is given by ( 13), are removed from the disc and, accordingly, do not contribute to torque exerted on the disc from the side of perturbed.Taking into account again that this is important only when ∆r 0 is small, this can be done by redefinition of In this way, we finally obtain 2.5 The perturbing potential In general, the Newtonian gravitational potential of a binary system, with masses M (primary) and m (secondary) in a coordinate frame with the origin at the primary position, has the form (see e.g.Larwood & Papaloizou (1997)) where r is the radius vector of a particular gas element in the accretion disc, D is the vector pointing from the central black hole to the companion, with magnitude equal to the distance between the components, and we take into account that in the coordinate frame associated with the binary its orbit lies in the plane (X b , Y b ).Clearly, the last term in the expression for φ p is determined by the non-inertial character of the chosen coordinate frame.In what follows we are going to use below only quantities doubly averaged over the azimuthal angle Ψ and the binary orbital period P orb .One can easily see that after the averaging over the azimuthal angle the non-inertial term disappears.It is, therefore, neglected hereafter.

Calculation of nodal evolution
In order to calculate the nodal evolution of a disc's rings due to the presence of a perturber we need only the m = 1 term in Fourier decomposition of φ p over the angle Ψ, which is odd with respect to coordinate ξ, φ 1 p .We also average φ 1 p over binary's orbital period We are going to denote by bar any other other orbital averaged quantities.Let us also remind that we are interested only in quantities, which are linear in the inclination angles.
Substituting (18-23) in (30) we obtain where ∆r 0 is given by eq. ( 24), we use the correction ∝ β 2 b only in denominators in (31), determines the value of the perturbing potential in the orbital plane.It is not important for the calculation of nodal precession, but it will be used for the calculation of apsidal precession later in Section 3.2.Using (18-23) it is easy to see that X 0 X 1 + Y 0 Y 1 + Z 2 /2 = 0 in the linear approximation over the inclination angles.Therefore, we need to evaluate and orbit average the m = 1 term of the expression It is easy to see that this expression depends on Ψ only through cos(Ψ − Ψ * ), and so does the m = 1 term, φ 1 p .Therefore, for φ 1 p we obtain from ( 33) where where 2r 2 , where we take into account again that the correction is only important when D ≈ r.Note that (35) can be expressed in terms of complete elliptic integrals.Taking into account that D x is odd and D y is even with respect to time reversal, we obtain from (34) where and we take into account that terms, even with respect to time reversal, disappear after averaging over the orbit to obtain the last equality.
We substitute the explicit expressions for X 1 and Y 1 (20) in ( 35) and represent the result in the form where where δ = γ − γ b .
To perform the orbital averaging we need to integrate (37) over the orbit: φ1 p = 1 From the expression (37) it follows that in order to find φ1 p we need to evaluate two quantities Using the quantities K 1 and K 2 from the expression (37), it follows that orbit average perturbation of the potential has the form Note that it follows from the definition of I (see eq. ( 35)), that it has the dimension of the inverse cube of a distance, I ∝ r −3 .Accordingly, it is seen from ( 41), that K 1 and K 2 have the dimension of the inverse square, and, therefore, φ1 had the correct dimension of the inverse distance, φ1 ∝ r −1 .
Remembering that the perturbing gravitational potential is −Gmφ p , we can easily find the averaged force acting on a disc ring in the direction perpendicular to the ring plane, F ξ , as When the influences of hydrodynamical interactions and the relativistic effects on motion of the ring are neglected, and the influence of the perturber is treated as perturbation, the ring remains to be circular, but the Euler angles β and γ are changing.It is easy to find their evolution from the relation ∂U ∂Ψ = 1 2rΩ K F ξ , where we remind that Ω K = GM r 3 , and U = β sin Ψ − β γ cos Ψ, see e.g.Demianski & Ivanov (1997) and references therein.From this relation we have It is useful to calculate the corresponding time derivative of a complex quantity W = βe iγ .
Using eqns.( 39), ( 40) and ( 43) we easily obtain where . Solid, and dashed curves represent Ω 1 and Ω 2 in units of Ω K as functions of r, while the dotted curve is the Lense-Thirring frequency (in units of Ω K ) of a disc element on a circular orbit.Since Ω 2 changes its sign, only the absolute value is shown.Note that it is positive in both asymptotic limits.Solid and dashed curves with larger values of the argument correspond to smaller β b .We show β b = 0.01 labeled as 1± just above the curves with ± corresponding to Ω 1,2 , respectively.β b = 0.1 and 0.5 are labeled as 2± and 3± just below the curves with the same sign convention as in the previous case.Ψ ef f = 0 for all curves.0.01, 0.1 and 0.5 in Fig. 1.It is seen that the curves corresponding to β b = 0.01 and 0.1 are close to each other.This is because when β b is small, a contribution of the term proportional to β b to ∆ b in eq. ( 29) plays a smaller role than the term determined by the accretion radius, which depends only on the mass ratio.It is seen from eq. ( 44) that the term proportional to Ω 1 causes precession of a disc ring about the axis perpendicular to the orbital plane.In principal, when the second term proportional to Ω 2 is also taken into account, the ring dynamics could be non-trivial.As was discussed in e.g.Abod et al (2022), Aly et al (2015) and Zanazzi & Lai (2018) when precession of Ψ b is solely determined by the binary, under certain circumstances, it is possible to obtain configurations of free particles and particles in a disc precessing about the direction pointing towards the periastron of the orbit.However, in our case, the Einstein precession of Ψ b appears to be much faster than the one determined by the binary.The term proportional to Ω 2 in this situation is averaged out and causes little effect on the evolution of our system, see below.
3.1.1Expressions for K 1 and K 2 in asymptotic limits r/a ≪ 1 and r/a ≫ 1 In both limits r/a ≪ 1 and r/a ≫ 1, it is easy to evaluate K 1 and K 2 analytically.For that, we note that in these limits we can set ∆ 2 b = 0 and that the integral J defined in ( 35) is approximately equal to 3πx 2 in both limits.Using these facts, it can be shown that integrals entering (41) are elementary when r ≫ a and they can be shown to be reduced to Legendre polynomials of ǫ = √ 1 − e 2 when P. B. Ivanov and V. V. Zhuravlev r ≪ a.In the first case, we obtain while in the second case, we have Substituting these expressions in (42) we obtain when r ≪ a.

The rate of apsidal precession
In order to calculate the apsidal advance rate, ν, at first we calculate the angular frequency of circular motion, perturbed by the presence of the companion, Ω, according to the relation and epicyclic frequency, κ, which follows from the relation Then, the apsidal advance rate is determined by their difference, ν = Ω − κ.Assuming that the perturbing potential is small in absolute value, we can easily obtain where the bar stands again for the double average over the azimuthal angle and the orbit.
For our purposes, it is sufficient to use in (52) the doubly averaged value of the perturbing potential in the orbital plane (32).Thus, we neglect any contribution to the apsidal advance rate due to the difference between orbital plane and the plane of the disc.This value can be represented in the form where J(x) is given in (35).Note that we assume that the disc material always moves in the positive direction, while the motion of the perturber can be both in positive and negative directions.In the latter case, (52) should be multiplied by minus.
We show the absolute value of ν in Fig. 2 for e = 0.7, Ψ ef f = 0 and β b = 0.01, 0.1 and 0.5 together with the Einstein apsidal precession frequency.It is seen that Einstein's precession dominates at radial scales of interest.Therefore, for systems with the considered parameters, the contribution of gravitational potential of the perturber is rather insignificant.

NUMERICAL AND ANALYTIC ANALYSIS OF THE BEHAVIOUR OF THE TWISTED DISC IN THE CASE OF RELATIVELY LARGE VISCOSITY
As we have mentioned above, in this Paper we consider only the discs with a relatively large viscosity, α ∼ 0.1.This is because they are expected to evolve in manner, which is qualitatively different from the low viscosity case, since their alignment time T al is smaller than the system life time T GW .We also fix the orbital parameters of the binary to the values introduced above and consider two values of the binary's inclination β b = 0.1 and 0.5.P. B. Ivanov and V. V. Zhuravlev

A simple qualitative analytic model
Since we are going to consider the disc at times less than T al , in a zero approximation, it appears to be reasonable to neglect hydrodynamical interactions in the disc at all and assume that the disc is evolving only due to the torques exerted by the binary and black hole where we remind that W = β(t, r)e γ (t,r) , is the Lense-Thirring frequency and Ẇb is given by eq. ( 44).Since Ẇb depends on γ b and Ψ b as well as on β and γ through the dependencies of the frequencies Ω 1 and Ω 2 on the disc's Euler angles, eq. ( 54) is a rather complicated non-linear equation with time dependent coefficients.To make it analytically treatable, at first we average Ω 1 and Ω 2 over the angle Ψ ef f , thus considering Ω1,2 = 1 2π 2π 0 dΨ ef f Ω 1,2 and set β = β b in the resulting expressions.Secondly, we neglect the term proportional to Ω 2 in (44).
After making these simplifications, eq. ( 54) becomes elementary with the solution subject to the initial condition W(r, t = 0) = 0, where Ω b LT is given by (5).Equation ( 55) has a resonance at r = r r , such as ∆Ω(r r ) = 0.At radii r ∼ r r , the assumptions leading to (55) are broken and a more accurate treatment is required, see below.
When |∆Ω|t ≫ 1 the last term in the brackets in (56) has a large phase.Since ∆Ω is a function of r, this leads to very sharp variations of this term over r at a given time, which are supposed to be smeared out by the action of hydrodynamical interactions.We average over these variations, thus neglecting this term.In the end, we have an expression for a twisted disc, which is precessing with the Lense-Thirring frequency Ω b LT and is approximately stationary in the precessing frame, with the angle β varying with r as We compare below the expression (57) with the results of numerical simulations.
4.1.2A more accurate approach to the calculation of the quasi-stationary shape of the disc in the precessing frame Now let us consider the disc shape near r r .It follows from eq. ( 57) that β formally diverges at this radius.In order to remove this divergency, we take into account hydrodynamical interactions between the disc's rings.For that we use an appropriate modification of the simplest form of dynamical equations describing twisted discs proposed in Demianski & Ivanov (1997), see their equations ( 38) and ( 39).Note that these equations are approximately valid in the limits r ≫ 1 and α ≪ 1.They contain the variable W and an auxiliary complex variable A. Assuming that the disc is stationary in the frame precessing with Lense-Thirring frequency of the binary, Ω b LT , these variables can be represented as W = Ŵ(r)e iΩ b LT and A = Â(r)e iΩ b LT .Instead of the last term in eq. ( 39) of Demianski & Ivanov (1997) we use eq.( 55).In eq. ( 38) of Demianski & Ivanov (1997) we neglect the last term in the bracket on r.h.s, express Â in terms of Ŵ and substitute the result in eq. ( 39).The resulting equation can be further simplified in a region of |r − r r | ≪ r r .We represent ∆Ω in this region as where Ω * = −r r d∆Ω dr (r = r r ), and introduce two 'renormalised' dimensionless frequencies ω * and ω 1 , according to the rule It is evident from their definitions that both ω * and ω 1 do not depend on r.
In this way, we obtain It is easy to see from ( 60) that it can be reduced to a standard form of inhomogeneous Airy equation by rescaling Ŵ and x.Since ω * and ω 1 are proportional to δ −2 ∼ 10 6 they are expected to be quite large.That means that when |x| ≫ x * , where x * ≪ 1 is specified below, the first term on r.h.s. of ( 60) can be neglected and we have Using ( 58) and ( 59) it is evident that (57) tends to (61) in the limit x → 0. Thus, we are looking for a solution to (60), which must tend to (61) in the case |x| ≫ x * .
One can see that such a solution can be expressed in terms of Scorer's function Hi(z) (see e.g.Olver et al (2010) for its definition and properties) of a complex argument Note that Hi(z) admits the integral representation (e.g.Olver et al (2010)) which allows for a very efficient numerical evaluation of (62).The scale x * can be defined by the condition that the absolute value of the argument in ( 62) is equal to one.From that we have , where ∆Ω(r) and Ω1 (r) are understood to be functions of an arbitrary r.When r → r r x → x by its construction.The expression is reduced to (62) when |x| ≪ 1 and it is reduced to (57) when |x| ≫ x * .It provides an approximate analytic solution to our problem valid for any value of radius.We compare our analytic expression ( 64) with the results of our numerical work below, see Section 4.2.3.

Specific details of our numerical approach
The hydrodynamical model of a twisted disc which includes the gravitomagnetic action of the primary black hole is considered here employing the dynamical equations derived in Zhuravlev & Ivanov (2011).The equations (60-61) of Zhuravlev & Ivanov (2011) describe the geometrically thin relativistic accretion disc around a slowly rotating black hole.These equations are linear with respect to the disc inclination angle and based on the background model of the vertically isothermal flat disc as introduced by the Novikov-Thorne solution.It is assumed that the disc aspect ratio is determined by the Thomson scattering, which yields its profile according to equation (37) of Zhuravlev & Ivanov (2011).In the results shown below, it is assumed that the characteristic aspect ratio entering equation (37) of Zhuravlev & Ivanov (2011), δ * = 0.001.The standard form of the viscosity prescription is given by ( 35) of Zhuravlev & Ivanov (2011), which introduces the dimensionless viscosity parameter α.
In this Paper, we incorporate the new terms entering the right-hand side of equation ( 44) into the right-hand side of the twist equation ( 61) of Zhuravlev & Ivanov (2011).Also, the correction to epicyclic frequency due to the gravitational influence of the secondary black hole, see equation ( 52), is taken into account in equation ( 60) of Zhuravlev & Ivanov (2011), which describes the velocity perturbations induced by the disc twist.
The dynamical equations generalised in this way are numerically advanced using an explicit second-order scheme.The dynamical variables W and B6 are set on the spatial grids uniform in √ r with the grid nodes for B shifted with respect to the grid nodes for W by a half step.The quantities K 1 , K 2 and ν representing the influence of the pertuber are tabulated first on the separate three-dimensional mesh of its arguments, which are β ef f , ψ ef f and r.Next, they are evaluated in the grid nodes of the scheme by the successive linear interpolations along r, ψ ef f and β ef f .Since the disc shape and the binary orbit evolve, this procedure is repeated at each slice.The time step is limited by a value inverse to the largest absolute value of the frequency being the solution to local dispersion relation derived from the dynamical equations.We performed additional tests of the scheme in order to check that the convergence of the numerical solution is achieved for the case α = 0.1 and, additionally, for α = 0.01.For that, the spatial grid step must be significantly smaller than the disc thickness in the vicinity of the inner boundary in the case α = 0.01.However, it may even exceed the inner disc thickness in the case α = 0.1.Thus, we find that the case of smaller viscosity is far more demanding for the scheme.The regularity boundary conditions, B = 0 and dW/dr = 0, are imposed close to the inner boundary of the disc, where the surface density of the standard accretion disc vanishes.Note that in the limit of slow primary black hole rotation, the inner boundary of a twisted disc is formally set to its Schwarzschild value, r = 6GM/c 2 , see the analysis in Zhuravlev & Ivanov (2011).
4.2.2 Results of numerical calculations of the case with α = 0.1 We consider numerical runs with χ = ±0.5, ±0.25, and β b = 0.1, 0.5.We evolve all runs up to the end times order of or larger than the evolution time T GW .
In Figs. 3 and 4 we show ratios β/β p on r for different times t, for the case |χ = 0.5|, for β b = 0.1 and β b = 0.5, respectively.Solid (labeled as 1), dashed (labeled as 2), dot dashed (labeled as 3) and dot dot dashed (labeled as 4) curves correspond to t , respectively.Note that all labels are always in the bottom  LT t the evolution is practically harmonic, with the phase of W 1 differing from that of W 2 by π/2 and the time period approximately equal to 2πΩ b LT −1 .We have checked that the same type of evolution is observed for other radii and other calculation parameters.This confirms that the disc is indeed precessing with its precessional frequency equal to the Lense-Thirring frequency of the binary as suggested by our analytic model.The case |χ| = 0.25 is qualitatively similar to the case |χ| = 0.5.Therefore, we show only the dependencies of β on r for particular moments of time together with the analytical curves represented again by dotted lines.These are shown in Figs. 6 and 7, which are analogous to Figs. 3 and 4. As seen from these figures, the agreement between analytical and numerical results looks slightly better for the smaller absolute value of χ.

A more detailed comparison of our analytic and numerical results
Since the expression for Ω 1 given by equation ( 45) depends both on r and the angle Ψ ef f , in order to compare our results with the analytic expression (64), we remind that we average numerically calculated profiles of Ω 1 over this angle in order to obtain Ω1 entering (64).We use the case β b = 0.1 and χ = 0.5 for the comparison, the result is shown in Fig. 8 for the case β b = 0.1 together with the corresponding result of numerical computation taken at some late computational time when the numerical distribution of β is almost stationary.As seen from this figure, the numerical and analytic curves almost coincide at radii larger than the radii corresponding to the maximum of β.These maximal values agree with the accuracy order of 10 per cent.Also, the radius corresponding to the maximum value of the analytic curve is shifted towards smaller values of r/r g by the same amount.Overall, we see that our rather simple analytic approach gives quite a reasonable approximation to the numerical result.

4.3
The influence of the disc tilt and twist on timing of intersections between the perturber and the disc Now let us briefly consider how possible conditions of the intersections between the perturber and the disc are modified by the presence of the disc's tilt and twist.For that, we look for the intersections of the disc by the pertuber by comparing the vertical coordinate of the perturber in the coordinate system associated with the primary, Z b h , calculated along its orbit, with the corresponding height of the disc, Z d h .Z d h (t) is strictly defined as the vertical coordinate of a point on the disc belonging to the line perpendicular to the equatorial plane, which crosses the orbit at the point with the coordinate Z b h (t).Since both quantities are taken along the orbit, it is clear that the perturber crosses the disc when and only when Z b h (t) = Z d h (t).It is also evident that when the disc is flat the condition of intersection is just Z b h (t) = 0.The instant location of the perturber with respect to the ascending node of its orbit is evaluated as Ψ = Ψ b + n b , where n b is the current true anomaly of the perturber.The latter is determined numerically employing an iterative solution of the Keplerian equation ( 16) and the standard rela-P.B. Ivanov and V. V. Zhuravlev tion between the true and the eccentric anomalies.We assume below that n b (t = 0) = 0, and that n b = 0 corresponds to the passage of periastron.In the linear approximation over the angle β b , we have where r b (n b ) is the current distance between the holes.
Similarly, in the linear approximation over β, we have where Ψ tw is the position angle of the disc element located above (or below) the perturber with respect to the ascending node of the disc ring specified by the angles β(r b ) and γ(r b ).In the same approximation we can use which determines Z d h together with equation ( 66).LT , respectively, assuming that Ψ b (t = 0) = π/4 and γ b (t = 0) = 0.For our orbital parameters, it gives the minimal time interval of the order of one year between the successive intersections of the disc in the flat disc case, which is broadly consistent with the PM model of OJ 287.Since we are going to discuss only the principal importance of the disc tilt and twist, this appears to be sufficient for our study.
In Figs. 9 we show the dependencies of Z b h and Z d h on time during one orbital period, for the disc model with β b = 0.1 and |χ| = 0.5 evolved for a sufficiently long time to the quasi-stationary state.We see that both the number of crossings and the time passed between particular crossings of the disc differ very significantly from the flat disc case.Formally, we have six crossings of the disc instead of two in the flat case.Note, however, that certain crossings happen very close to each other, as e.g. the ones close to Ω b K ∆t ≈ 3 in the prograde case, these particular pairs could merge if we take into account e.g. the finite thickness of the disc or the fact that the perturber is expected to expel the disc material within its accretion radius from the orbit.We further illustrate this case in Fig. 10 showing the evolution of Z b h and Z d h on the plane (r h , Z h ), where r h is the distance from the primary.One can see from this figure that the section of this disc (r h , Z d h ) has rather non-trivial from, with a 'butterfly-like' feature close to r h /a = 0.5.This feature is explained using the results of Section 4.1.From these results it follows that close to the resonance radius r r the disc's nodal angle γ rotates by π.The same effect is also observed in our numerical calculations.This rotation leads to a very sharp variation of the disc height along the orbit when the perturber's distance is close to r r .We confirm that our main conclusion that the number of disc's crossings per period is significantly modified is robust in Figs.11-13.In Fig. 11 we illustrate the case with the same value of black hole rotational parameter |χ| = 0.5, but for the orbit with β b = 0.5, while in Fig. 12 and 13 we show the results for the case |χ| = 0.25, for β b = 0.1 and β b = 0.5, respectively.One can see that the qualitative situation remains the same, although the numbers of crossings and the crossing times differ from case to case.Additionally, on the left panel of Fig. 11 we show the case of a disc model with the same parameters, but with Ω 2 set artificially to zero.This corresponds to averaging of our disc model over a time period order of Ω −1 E .It is seen again that that the situation remains very similar to the general case.Finally, in Fig. 14 we show a three-dimensional picture of the disc, the perturber's orbit and their intersections.One can clearly see a 'spiral rim' at radii close to the resonance radius in this Fig.

CONCLUSIONS AND DISCUSSION
In this Paper we considered effects caused by interaction of an eccentric inclined binary black hole with a very thin (h/r ∼ 10 −3 ) accretion disc.The orbital parameters of the binary were assumed to be the same as in the so-called 'precessing massive' (PM) model of OJ 287 (see e.g.Valtonen et al (2023) and references therein).Namely, it was assumed that the binary consists of two black holes with masses M ≈ 2 • 10 10 M ⊙ and m ≈ 10 8 M ⊙ on a very eccentric orbit with eccentricity e = 0.7 and the semi-major axis a ≈ 60 GM c 2 corresponding to the orbital period in the observer's frame ≈ 12yr.A distinctive feature of such a system is that its evolution time due to the emission of gravitational waves T GW is smaller than the disc relaxation time T ν , which suggests that the disc is unable to react to the presence of the binary by e.g.opening a gap inside its orbit.Moreover, when the disc viscosity parameter α > 0.01, the evolution time is smaller than the time T al associated with alignment of the disc to a quasi-stationary twisted configuration.It is this case that has been studied in this Paper in some detail.
We calculated the torque exerted on the disc by an inclined binary, assuming that the mass ratio q is small and one can use a perturbation theory.Technically, we double averaged the expression for the gravitational force responsible for this torque over the binary period and over a period of a particular disc element on a slightly perturbed circular orbit.The same procedure is used to obtain a correction to the apsidal precession rate of the disc elements.We assumed that the double averaging procedure can be used to describe the evolution of the disc over characteristic time exceeding the binary orbital period and typical orbital periods of the disc elements.It was also assumed that the perturbing component does not produce torques on the disc inside its 'accretion radius', since © 2010 RAS, MNRAS 000, 1-37 gas inside this radius is likely to be removed from the system.In general, the corresponding torque and the apsidal precession rate are obtained numerically, but simple asymptotic expressions can be found for the disc's elements with radii much smaller or much larger than a.We have checked that these expressions agree with similar expressions obtained in other studies.We used the torque term and the apsidal precession rate together with the standard expressions for these quantities induced by a relatively slowly rotating central black hole (primary) to calculate numerically the evolution of the disc tilt and twist from the flat configuration initially aligned with the equatorial plane of the primary.We set α = 0.1 for our runs and used four values of the primary rotational parameter χ = ±0.25 and χ = ±0.5,where the plus (minus) sign corresponds to prograde (retrograde ) rotation of the primary with respect to the direction of motion of the disc's elements.We found that in all cases the disc relaxes to a non-trivial quasi-stationary configuration in the frame rotating with the Lense-Thirring frequency of the orbit Ω b LT .The corresponding relaxation time of the order of the corresponding period P b LT = 2π/Ω b LT is much smaller than T GW .After the period of relaxation, the disc's inclination angle β is almost stationary at scales order of a.It is shown that β significantly varies inside the orbit.When α = 0.1, its value in the vicinity of its periastron is much smaller than the orbital inclination angle β b .At radii between the periastron and apoastron of the orbit, β becomes a few times larger than β b while approaching back to β b in the vicinity of apoastron.
Our preliminary analysis of the conditions of the crossing times made in Section 4.3 shows that there are, in general, more intersections of the perturber with the twisted disc per one orbital period that in the flat disc case.This happens in part due to the large disc's tilt at scales of the order of the perturber's orbit, and, in part, due to a large disc's twist near 'resonance' radius r r , where the precession frequency of a disc element both due to the actions of the binary and the primary coincides with Ω b LT .Close to this resonance the disc's nodal angle γ rotates by π.Clearly, our results must be taken into account in the modelling of OJ 287 and similar sources.
In Section 4.1.1we proposed an elementary analytic model for the distribution of the disc inclination over r, which agrees with the numerical results surprisingly well apart from scales close to r r .A more accurate, albeit more involved treatment was considered in Section 4.1.2.It describes the quasi-stationary shape of the disc in the frame precessing with the Lense-Thirring frequency of the binary at all radii, including the ones close to the resonance.We compared the analytic and numerical results in Section 4.2.3, and showed that those are in good agreement.It is important to stress that one of the reasons of the agreement is that the term proportional to Ω 2 in (44) is not very important.It is proportional to e 2iΨ b , and Ψ b evolves on a relatively fast timescale of Einstein's apsidal precession.However, the disc tilt and twist mainly evolve on a longer timescale order of (Ω b LT ) −1 , so that a contribution of the term proportional to Ω 2 is averaged out.Therefore, we can assume that the timescale of averaging the potential perturbations in Section 3 can be shorter than, but of the order of (Ω b LT ) −1 .
Our analytic results can be helpful in calculating some approximate shapes of the disc without the use of numerical integration of the equations for the disc's tilt and twist for other parameters of the system.For that, one can numerically evaluate the distribution of Ω 1 using the results obtained in Section 3.1 for any desired parameters of the binary and the disc and use results of Section 4.1 to find the disc shape.This approach can, for example, be used in an accurate modelling of light curves of OJ 287 in the framework of the PM model, allowing for the possibility of multiple disc crossings per orbital period found in this work.Of course, our approach is valid only for systems with sufficiently large values of α ≫ δ.
In this Paper we did not make attempts to compare our model with observations.This requires a detailed survey of the parameter space which deserves a separate study with careful consideration of a lot of other factors influencing the determination of the light curve (see Dey et al (2018) and Zwick & Mayer (2023) for a state-of-the-art modelling of the timing of OJ287 flares).Note, however, that it seems to be clear that taking into account of our results may lead to a significant revision of the parameters of the system in the framework of PM model.Alternatively, it can provide arguments in favour of other models, such as the models based on the jet activity (e.g.Britzen et al ( 2023)) or, for example, the models of eccentric binary black holes of comparable masses immersed in a circumbinary disc, e.g.Hayasaki, Mineshige & Sudou (2007) and Hayasaki, Mineshige & Ho (2008).
The numerical calculations should also be extended to smaller values of α.As we briefly discussed in the text, smaller values of α turn out to be demanding from the numerical point of view, and, therefore, numerical simulations of the discs with lower viscosity are left for a future This problem is beyond the scope of this work and is left for the future.
Of course, our model should be checked with either SPH or grid based three-dimensional simulations.For that, one does not need to consider very thin discs or very long simulation times.
and the angle φ orb is defined by the conditions Dx D = cos φ orb and Dy D = sin φ orb .The correction ∆ 2 b is expected to be important only when Ψ ≈ Ψ * and r ≈ D. Therefore, we can set D = r and Ψ = Ψ * in the expression (25), and rewrite it in the form

W
b = β b e iγ b and * stands for the complex conjugate.We show Ω 1 and the absolute value of Ω 2 together with Lense-Thirring frequency of a circular disc element 5 for e = 0.7 and β b =

Figure 2 .
Figure2.The absolute value of ν in units of Ω K is shown as a function r/a together with the Einstein precession frequency in the same units represented by the dashed curve.Curves with larger values of the argument correspond to smaller values of β b .Note that at large values of r/a ν is negative.We label the curves corresponding to β b = 0.01 and 0.1 by 1 and 2 just above the curves, and the curve corresponding to β b = 0.5 is labeled as 3 just below the curve.

Figure 3 .Figure 4 .Figure 5 .
Figure 3. Left panel.The numerical dependencies of β/β b on r for different times as well as the analytic dependency (57), see the text for a description of particular curves.The case β b = 0.1 and χ = 0.5 is shown.Right panel.The same as the left panel, but for the retrograde black hole rotation χ = −0.5.

Figure 6 .
Figure6.The same as Fig.3, but for smaller |χ| = 0.25.Also note that the dot dot dashed curves correspond to a slightly smaller t = 8620Ω b K

Figure 8 .
Figure 8.The ratio of the inclination angle β to the binary inclination angle β b is shown as a function of r.The solid and dashed curves represent the result of a numerical calculation and our analytic expression (64), respectively.Note we use χ = 0.5, rr ≈ 0.5a, Ω 1 (rr) ≈ 2.9 • 10 −3 Ω b K Z d h and Z b h are shown for a particular orbital period in Figs.(9-13).The instant values of Ψ b and γ b are obtained according to the uniform apsidal and nodal precessions of the binary orbit with the frequencies Ω E and Ω b

Figure 9 .Figure 10 .
Figure 9. Left panel.The solid curve shows the height of the disc above the equatorial plane of the primary black hole at the location of the perturber, Z d h , as a function of time; the dashed curve shows the simultaneous height of the perturber, Z b h .Both quantities are expressed in units of a.The dotted line Z h = 0 shows the equatorial plane of the primary black hole.The disc parameters are β b = 0.1 and χ = 0.5 and the disc evolution time is approximately 9 • 10 3 Ω b K −1 .Right panel.The same as the left panel, but for the retrograde rotation of the primary χ = −0.5.Note that, for convenience, we shift the time origin in our figures in such a way, that the orbit periastron corresponds to Ω b K ∆t ≈ 3.

Figure 11 .
Figure 11.The same as Fig.9, but for the case β b = 0.5.Additionally, here we plot Z d d corresponding to the model of disc evolution with Ω 2 artificially set to zero by the dot dashed curve on the left panel.

Figure 14 .
Figure14.A three-dimensional image of disc and its intersections with the binary orbit.The trajectory of the secondary is shown by the black line with an arrow indicating the direction of the binary orbit.The disc parameters are β b = 0.5 and χ = 0.5 and the disc evolution time is close to that used in Figs.9-13.Note that for this particular orbit its apsidal and nodal lines are close to each other.Top panel.The observer's line of sight is in the equatorial plane of the primary while it has an angle −45 degrees with respect to X h -axis.Bottom panel.The observer's line of sight is inclined by 45 degrees with respect to the equatorial plane of the primary while its projection onto this plane has an angle −20 degrees with respect to X h -axis.
work.It is interesting to mention that the low viscosity discs can, in principle, have a qualitatively different behaviour.Namely, the presence of the resonance discussed in Section 4.1 suggests a possibility of launching warp waves near the resonance by a mechanism similar to the well-knownGoldreich & Tremaine (1979) mechanism of launching density waves near Lindblad resonances.