Abstract

We investigate the evolution of low mass (Md/Mb = 0.005) misaligned gaseous discs around eccentric supermassive black hole (SMBH) binaries. These are expected to form from randomly oriented accretion events on to a SMBH binary formed in a galaxy merger. When expanding the interaction terms between the binary and a circular ring to quadrupole order and averaging over the binary orbit, we expect four non-precessing disc orientations: aligned or counter-aligned with the binary, or polar orbits around the binary eccentricity vector with either sense of rotation. All other orientations precess around either of these, with the polar precession dominating for high eccentricity. These expectations are borne out by smoothed particle hydrodynamics simulations of initially misaligned viscous circumbinary discs, resulting in the formation of polar rings around highly eccentric binaries in contrast to the coplanar discs around circular binaries. Moreover, we observe disc tearing and violent interactions between differentially precessing rings in the disc significantly disrupting the disc structure and causing gas to fall on to the binary with little angular momentum. While accretion from a polar disc may not promote SMBH binary coalescence (solving the ‘final-parsec problem’), ejection of this infalling low-angular momentum material via gravitational slingshot is a possible mechanism to reduce the binary separation. Moreover, this process acts on dynamical rather than viscous time-scales, and so is much faster.

1 INTRODUCTION

It is widely accepted that most massive galaxies host a supermassive black hole (SMBH) at their centre. Galaxy mergers, expected from the hierarchical structure growth scenario based on the Λ cold dark matter (ΛCDM) cosmological model, then result in the formation of SMBH binaries (Begelman, Blandford & Rees 1980). Both SMBHs of such a binary sink towards the galactic centre due to dynamical friction and form a hard binary (Merritt 2001). However, most such SMBHs appear to be single rather than binary SMBHs, implying that SMBH binaries quickly coalesce and merge. One process driving further binary shrinking is slingshot interactions with stars in the ‘loss cone’: those on orbits intersecting with the binary (Saslaw, Valtonen & Aarseth 1974). Since the slingshot mechanism ejects the stars, the ‘loss cone’ needs to be replenished in a relatively short time-scale in order to shrink the binary all the way down to separations ≲10−2 pc where gravitational waves are expected to drive coalescence. For spherical collisionally relaxed stellar systems, it is thought that the slingshot mechanism stalls well before reaching this separation, resulting in what is known as ‘the final-parsec problem’ (Milosavljević & Merritt 2003; Berczik, Merritt & Spurzem 2005).

Potential stellar dynamical solutions have been sought for gas-poor systems. Berczik et al. (2006) studied SMBH binaries evolution in realistic triaxial rotating galaxies and found that the galaxies supply stars on centrophilic orbits refilling the loss cone at a high enough rate to prevent the SMBH binary from stalling and that complete coalescence is achieved in less than 10 Gyr. Khan et al. (2013) found that for axisymmetric galaxies with axis ratio c/a ≤ 0.8 the hardening rate is 25 times faster than for spherical galaxies. Self-consistent N-body simulations of merging galaxies containing SMBH found binary hardening rates much higher than idealized spherical models and sufficient to shrink the binary to the gravitational wave coalescence regime (Khan, Just & Merritt 2011; Gualandris & Merritt 2012).

Interactions with circumbinary gas discs may change the evolution of the SMBH binary. The mass of such a disc is uncertain and depends on its formation. A galaxy-rich merger can channel large amounts of gas towards the centre. If this gas can cool efficiently and avoid fragmentation and substantial star formation, a disc with mass comparable to that of the binary may form.

For a prograde disc, spiral density waves in the disc driven by the outer Lindblad resonance with the binary transport angular momentum away from the binary (Goldreich & Tremaine 1979). This mechanism is efficient if the disc reaches very close to the binary, so that it occupies the resonance, and is sufficiently massive for the angular momentum absorption not to result in an expansion of the inner disc edge (Escala et al. 2005; MacFadyen & Milosavljević 2008). For a disc with Md = 0.2 Mb, Cuadra et al. (2009) found that binary orbital decay can stall because the disc expands due to absorption of angular momentum from the binary, severely slowing further angular momentum exchange (see also Lodato et al. 2009).

Apart from the classical density-wave mechanism, the infall of gas from the inner edge of the disc into the cavity can be important (Roedig et al. 2012; Roedig & Sesana 2014). The binary may either eject such infalling gas via a gravitational slingshot whereby losing angular momentum and energy, or capture it on to an accretion disc around either component, which adds to the binary angular momentum. The binary evolution is determined by the competition between these two effects and it remains unclear, which one wins in the long term.1

For a retrograde coplanar disc, the lack of orbital resonances allows the disc to extend to small radii. This enables the binary to accrete or capture material with negative angular momentum (Nixon et al. 2011a). If Md ∼ Mb, this may suffice to achieve coalescence (Roedig & Sesana 2014).

In reality, discs with mass in excess of their aspect ratio times the binary mass are gravitationally unstable and hence, due to the short cooling time in these discs, fragment and form stars much faster than binary coalescence (Gammie 2001; Goodman 2003; Levin 2007). The numerical treatment of fragmentation, star formation, and stellar feedback is extremely challenging. In all of the aforementioned simulations with such massive discs, these processes have simply been suppressed (by assuming slow cooling which prevents star formation), overestimating the efficiency of disc-driven binary coalescence. Although star formation will rob the disc of a significant amount of gas, the newly formed stars may still contribute to binary orbital decay (e.g. Sesana, Haardt & Madau 2007, 2008), though less so than the gas owing to the lack of an efficient dissipation mechanism to reduce their pericentres.

A more likely scenario than binary coalescence driven by the interaction with a single massive disc is the repeated interaction with low-mass discs resulting from the infall and tidal disruption of molecular clouds on to the binary. Such discs are expected to have masses 105–6 M typical of molecular clouds, small compared to the typical mass 106–9 M of a SMBH binary. Nixon et al. (2011a) studied retrograde discs of this type and found that they are very efficient in reducing the binary angular momentum through accretion of gas with negative angular momentum on to the secondary black hole. This enhanced accretion on to the secondary black hole increases the binary's eccentricity, decreasing the pericentre distance in the process, and coalescence is achieved when a mass comparable to the secondary black hole has been accreted.

If accretion events in galactic nuclei are chaotic and randomly oriented (King & Pringle 2006, 2007; King, Pringle & Hofmann 2008), we expect the formation of misaligned circumbinary discs around SMBH binaries. In the case of a circular SMBH binary, the interaction between the misaligned disc and the binary is similar to Lense–Thirring precession on an accretion disc around a spinning black hole (Bardeen & Petterson 1975; Pringle 1992; Scheuer & Feiler 1996). King et al. (2005) showed that the induced differential precession will cause a misaligned disc to counter-align with the black hole spin provided:
(1)
where |$\boldsymbol{J}_{\mathrm{d}}$| and |$\boldsymbol{J}_{\mathrm{h}}$| are the disc and black hole angular momenta, respectively, and θ is the angle between them. The disc will co-align with the black hole spin if this relation is not satisfied. Nixon, King & Pringle (2011b) showed that the same analysis applies to the case of a misaligned disc around a binary (though the precession rate is slightly different). Thus, if counter-alignment is stable (Nixon 2012), this mechanism can provide a solution to the final-parsec problem by supplying retrograde discs to achieve coalescence.

Recently, Nixon, King & Price (2013) performed 3D hydrodynamical simulations of circumbinary discs around a circular binary for various tilt angles θ. In addition to co- and counter-alignment, they found that in many cases the discs is torn into distinct rings which precess almost independently (Nixon et al. 2012). The precessing rings, which have partially opposed angular momentum, may interact causing partial cancellation of their angular momenta and thus gas infall close to the binary. This disc tearing significantly increases the accretion rate and may play an important role in promoting the binary final coalescence.

Those studies considered the case of a circular binary interacting with a circumbinary disc, when disc precession is only around the pole of the binary plane. In this study, we consider the more general situation of an eccentric binary. For a SMBH binary formed via a galaxy merger, we expect high eccentricities in many cases (Aarseth 2003; Khan et al. 2011; Wang et al. 2014). Moreover, retrograde accretion on to a circular binary naturally results in eccentricity growth as discussed earlier. One important effect of binary eccentricity is to make the time averaged binary potential triaxial rather than axisymmetric as for a circular binary. Previous studies have shown that misaligned discs in triaxial galaxies can precess around both the major and the minor axes (Steiman-Cameron & Durisen 1984; Thomas, Vine & Pearce 1994).

The paper is organized as follows. In Section 2 we present analytic results for a simple orbit-averaged model for the binary–disc interaction up to quadrupole order. Section 3 describes the set-up of our 3D hydrodynamical simulations, results of which are presented in Section 4 and discussed in Section 5. Finally, we summarize and conclude in Section 6.

2 BINARY–DISC QUADRUPOLE INTERACTION

The dynamics of a circumbinary gaseous ring orbiting an eccentric binary is not analytically treatable, even without considering any dissipation. However, useful insight can be obtained by (i) truncating the binary gravitational potential at quadrupole order, (ii) assuming that the ring is circular, (iii) time-averaging over the binary orbit, and (iv) neglecting dissipation. Assumptions (i) and (ii) are valid as long as the ring is sufficiently distant from the binary, while assumption (iii) requires that orbital resonances between ring and binary are not important.

The monopole of the gravitational interaction results in Keplerian motion of the ring around the binary centre of mass, while the quadrupole describes the lowest order deviation of the binary from a central point mass.

Recently, Naoz et al. (2013) have used Hamiltonian perturbation theory to obtain the equations for the secular evolution of a hierarchical triple up to octopole order. For a circular outer binary, their results are equivalent to the situation of a circular circumbinary ring. We now summarize the relevant relations (obtained in Appendix A with Newtonian dynamics, but otherwise equivalent to those of Naoz et al.) in terms of vectors rather than orbital elements to describe the system.

The binary is parametrized by its mass ratio qm2/m1 ≤ 1, total mass M = m1 + m2, semimajor axis a, specific angular momentum |$\boldsymbol{h}$|⁠, and eccentricity vector |$\boldsymbol{e}$|⁠. Let |$\boldsymbol{R}\equiv \boldsymbol{x}_1-\boldsymbol{x}_2$| the instantaneous binary separation vector, then
(2)
and
(3)
The vector |$\boldsymbol{e}$| is conserved for the binary orbit and points from the centre of mass to periapse (hence is always orthogonal to |$\boldsymbol{h}$|⁠). Its magnitude is the orbital eccentricity and is related to that of |$\boldsymbol{h}$| by
(4)

2.1 Ring evolution

The circular circumbinary ring is parametrized by its mass m, radius r, and pole |$\hat{\boldsymbol{l}}$|⁠. The latter is the unit vector in direction of the ring's specific angular momentum |$\boldsymbol{l}$|⁠, which has amplitude |$l=\sqrt{G{(}M{+m)}r}$|⁠. The ring radius must satisfy r > a(1 + e)/(1 + q) for the quadrupole approximation to be valid (and in order to avoid collision with a binary component). Note that the tilt angle θ of the ring with respect to the binary satisfies |$\cos \theta =\hat{\boldsymbol{l}}{\cdot }\hat{\boldsymbol{h}}$|⁠.

The quadrupole interaction energy between binary and ring, averaged over both the binary orbit and the ring, is
(5)
with |$\omega =\sqrt{G{(}M{+m)}/r^3}$| the orbital frequency of the ring, in agreement with equation (22) of Naoz et al. (2013). The time-averaged binary quadrupole torques the ring according to
(6)
with the vector
(7)
From equation (6) we have |$\boldsymbol{l}\cdot \dot{\boldsymbol{l}}=0$|⁠, i.e. |$\dot{l}=0=\dot{r}$| and the ring is merely precessing (this is no longer true at octopole and higher order, when |$\dot{l}\ne 0$|⁠; see Naoz et al. 2013).
Since |$\omega m r^2\boldsymbol{\Theta }=\mathrm{\partial} {\langle {E_{\mathrm{br}}}\rangle }/\mathrm{\partial} \hat{\boldsymbol{l}}$|⁠, equation (6) implies
(8)
Thus (in the assumed approximation), no energy is exchanged between binary and a single ring, but only angular momentum (if the binary interacts with several rings, the individual interaction energies with each ring are no longer conserved).

2.2 Binary evolution

The torque of the binary from the ring can be worked out analogously to that of the ring from the binary. After averaging over the binary orbit, we obtain
(9)
In particular, the total angular momentum, |$M\boldsymbol{h}+m\boldsymbol{l}$|⁠, is conserved at quadrupole order. For the case m ≪ M considered here, the orientation |$\hat{\boldsymbol{h}}$| only varies slightly even if the disc orientation |$\hat{\boldsymbol{l}}$| undergoes large changes.

For a circular binary |$\boldsymbol{\Theta }$| is parallel to |$\hat{\boldsymbol{h}}$| such that |$\dot{\boldsymbol{h}}{\cdot }\boldsymbol{h}=0$|⁠, i.e. |$h=|\boldsymbol{h}|$| is conserved and the binary is merely precessing (with an amplitude that is smaller than that of the disc by a factor m/M). This fact together with conservation of total angular momentum was the basis of the analysis by Nixon et al. (2011b).

For an eccentric binary, the evolution of |$\boldsymbol{h}$| is not simply a precession and h not conserved. Instead, we find
(10)
with |$\Omega =\sqrt{GM/a^3}$| the binary orbital frequency. Thus, h remains unchanged only if e = 0 (circular binary), or if |$\hat{\boldsymbol{l}}$| is perpendicular to either |$\hat{\boldsymbol{e}}$| or |$\hat{\boldsymbol{k}}$|⁠, i.e. if either |$\hat{\boldsymbol{e}}$| or |$\hat{\boldsymbol{k}}$| are in the ring plane. Otherwise, h oscillates, since |$\hat{\boldsymbol{l}}\cdot \hat{\boldsymbol{k}}$| oscillates around zero under the ring precession.
The change of the eccentricity vector is
(11)
and the corresponding change in eccentricity
(12)
in agreement with equation (A34) of Naoz et al. (2013), but also with equation (10) in conjunction with equation (4). In addition to the precession of the orbital plane and the oscillation of the eccentricity (both already described by equation 10), the binary also undergoes apsidal precession with rate
(13)
which is prograde for near-planar disc orientations (when |$|\hat{\boldsymbol{l}}\cdot \hat{\boldsymbol{h}}|\sim 1$|⁠), but retrograde for near-polar discs (when |$|\hat{\boldsymbol{l}}\cdot \hat{\boldsymbol{e}}|\sim 1$|⁠).

2.3 Ring precession

The rates of change of the directions of the binary and ring angular momenta satisfy
(14)

and a similar relation holds for |$|\mathrm{d}\hat{\boldsymbol{e}}/\mathrm{d}t|$|⁠. Thus, as long as m ≪ M, the binary orientation changes only very little and/or much more slowly than that of the ring (except for extreme binary eccentricities when |$Mh\propto \sqrt{1-e^2}$| can be small). We therefore consider in this subsection the limit m/M → 0 when the binary orientation and eccentricity are constant.

Then, equation (8) implies that the ring precesses along curves of constant 〈Ebr〉. Isolated minima and maxima of 〈Ebr〉 denote stable, non-precessing ring orientations. In the presence of dissipation (due to viscosity in the disc), these orientations are attractors, i.e. the dissipative damping of the precession eventually aligns the pole |$\hat{\boldsymbol{l}}$| of the ring with the extrema of 〈Ebr〉 (Steiman-Cameron & Durisen 1984). For any e < 1, the orientations |$\hat{\boldsymbol{l}}=\pm \hat{\boldsymbol{h}}$| are isolated minima of 〈Ebr〉 and correspond to coplanar ring orientations either co- or counter-rotating with the binary.

For a circular binary (e = 0), these are the only stable orientations, but all polar orbits (θ = 90°) maximize 〈Ebr〉. |$\boldsymbol{\Theta }$| is parallel to |$\boldsymbol{h}$| and ring precession is circular: |$\hat{\boldsymbol{l}}$| describes a circle around either of the stable orientations, see also the left-hand plot in Fig. 1. The precession rate is lower than the orbital frequency by the factor 3qa2cos θ/4r2(1 + q)2. This is the situation previously studied by Nixon et al. (2011b). We now turn to the more general case of an eccentric binary.

Precession paths for the direction $\hat{\boldsymbol{l}}$ of the angular momentum of a dissipationless circular ring of negligible mass orbiting a binary with eccentricity e as indicated. The binary orbits counterclockwise in the plane perpendicular to its specific angular momentum $\boldsymbol{h}$ with periapse in the direction $\hat{\boldsymbol{e}}$. For a circular binary (e = 0, left), $\hat{\boldsymbol{l}}$ always precesses around $\boldsymbol{h}$ in a retrograde sense. For eccentric binaries, prograde polar precession (blue) around $\boldsymbol{e}$, the long axis of the time-averaged binary potential, is also possible. The regions of polar and azimuthal precession are separated by two great circles (black). The four ring orientations $\hat{\boldsymbol{l}}=\pm \hat{\boldsymbol{h}}$ and $\hat{\boldsymbol{l}}=\pm \hat{\boldsymbol{e}}$ are stable (non-precessing), while the orientations $\hat{\boldsymbol{l}}=\pm \hat{\boldsymbol{k}}$ are unstable. Dissipation would damp the precession and eventually align the ring with one of the four stable orientations. In the case of a massive ring, the binary orbit evolves too: the vectors $\boldsymbol{h}$ and $\boldsymbol{e}$ oscillate and precess, and $\boldsymbol{e}$ and $\hat{\boldsymbol{k}}$ rotate around $\boldsymbol{h}$.
Figure 1.

Precession paths for the direction |$\hat{\boldsymbol{l}}$| of the angular momentum of a dissipationless circular ring of negligible mass orbiting a binary with eccentricity e as indicated. The binary orbits counterclockwise in the plane perpendicular to its specific angular momentum |$\boldsymbol{h}$| with periapse in the direction |$\hat{\boldsymbol{e}}$|⁠. For a circular binary (e = 0, left), |$\hat{\boldsymbol{l}}$| always precesses around |$\boldsymbol{h}$| in a retrograde sense. For eccentric binaries, prograde polar precession (blue) around |$\boldsymbol{e}$|⁠, the long axis of the time-averaged binary potential, is also possible. The regions of polar and azimuthal precession are separated by two great circles (black). The four ring orientations |$\hat{\boldsymbol{l}}=\pm \hat{\boldsymbol{h}}$| and |$\hat{\boldsymbol{l}}=\pm \hat{\boldsymbol{e}}$| are stable (non-precessing), while the orientations |$\hat{\boldsymbol{l}}=\pm \hat{\boldsymbol{k}}$| are unstable. Dissipation would damp the precession and eventually align the ring with one of the four stable orientations. In the case of a massive ring, the binary orbit evolves too: the vectors |$\boldsymbol{h}$| and |$\boldsymbol{e}$| oscillate and precess, and |$\boldsymbol{e}$| and |$\hat{\boldsymbol{k}}$| rotate around |$\boldsymbol{h}$|⁠.

For e > 0, the orientations |$\hat{\boldsymbol{l}}=\pm \hat{\boldsymbol{e}}$| are maxima of 〈Ebr〉, corresponding to polar rings (with opposite senses of rotation) around |$\hat{\boldsymbol{e}}$|⁠. For e < 1, |$\hat{\boldsymbol{l}}=\pm \hat{\boldsymbol{k}}$| are saddle points and correspond to polar rings around
(15)
the intermediate axis of the time-averaged binary potential. These latter ring orientations are unstable, i.e. small deviations will result in precession around any of the four stable orientations. For 0 < e < 1, ring precession is never circular: |$\hat{\boldsymbol{l}}$| describes a curve elongated towards the unstable orientations, rather than a circle. Azimuthal and polar precessions are retrograde and prograde, respectively. See Fig. 1 for a visualization of the precession paths.
The regions of polar and azimuthal precession are separated by the contours of 〈Ebr〉 passing through the saddle points. These separatrices are circular and shown as black in Fig. 1. The fraction of ring orientations undergoing polar precession are
(16)
At small e, this grows linearly (⁠|${\propto }\sqrt{20}e/\pi$|⁠) with eccentricity (see Fig. 2). Azimuthal and polar precession are equally likely for e = 6−1/2 ≈ 0.408.
Percentage of ring orientations undergoing polar precession as a function of binary eccentricity.
Figure 2.

Percentage of ring orientations undergoing polar precession as a function of binary eccentricity.

3 SIMULATION SET-UP

We perform a set of 3D smoothed particle hydrodynamics (SPH; Gingold & Monaghan 1977; Lucy 1977) simulations of geometrically thin accretion discs with different initial misalignment around an eccentric binary. We use a range of different binary eccentricities e = 0, 0.3, 0.6, and 0.9. The disc set-up is very similar to that used by Nixon et al. (2013): the disc is initially flat and extends from an inner radius of 2a to an outer radius of 8a with an inner thickness H/R = 0.01. We use a disc viscosity coefficient (Shakura & Sunyaev 1973) α = 0.1 which we set up using an appropriate SPH artificial viscosity coefficient αAV corresponding to our resolution (Lodato & Price 2010). All simulations start with four million SPH particles, while the binary is modelled using two equal mass sink particles with accretion radius of 0.05a. The disc initial surface density follows the profile Σ ∝ R−3/2, and we use a locally isothermal equation of state with sound speed cs ∝ R−3/4. These choices ensure a uniform vertical resolution (and hence uniform physical viscosity; see Lodato & Pringle 2007). We assume a disc mass of Md/Mb = 0.005 < H/R, which ensures that disc self-gravity is not important (we do not include gas self-gravity in our simulations, but we do self-consistently include the backreaction from the gas on the binary). The simulations were performed using our own code (Dehnen & Aly 2012), which implements an SPH scheme very similar to that used by Nixon et al. (2013), and we verified that our results for e = 0 agree with theirs. The disc has initial angular momentum direction
(17)
where ϕ is the twist angle of the disc. We ran a total of 118 simulations for e = 0, 0.3, 0.6, and 0.9; θ = 0°, 10°, 30°, 45°, 60°, 80°, 90°, 100°, 120°, 135°, 150°, 170°, and 180°; ϕ = 0°, −45°, and −90°. In the next section ϕ will be taken to be zero whenever it is not specified, we discuss the effects of varying ϕ separately.

We point out that the choice of a locally isothermal equation of state implies that the disc instantly radiates away all the heat gained from the viscous dissipation and shocks. This is justified if the cooling time is much shorter than the precession time. When this assumption does not hold, the disc thickness will increase and will be more able to resist breaking. We leave more advanced thermodynamic treatment to future investigation.

4 SIMULATION RESULTS

Fig. 3 shows snapshots after ≈95 binary orbits for the 20 simulations with initial tilt angles θ = 30°, 45°, 60°, 80°, and 90° initial binary eccentricities e = 0, 0.3, 0.6, and 0.9. As expected from the analysis in Section 2, the disc precesses around |$\hat{\boldsymbol{h}}$| for circular and low-eccentricity binaries, while for very high eccentricities the precession is predominantly around |$\hat{\boldsymbol{e}}$|⁠. In almost all cases, the disc breaks into distinct rings, which interact with each other and, depending on the details of each case, result in either co-, counter-, or polar alignment of the disc. In some cases the interaction between independently precessing such rings is very violent and disruptive, leading to ejection of gas. We now visit each possible outcome in detail.

Density rendering of the simulations after t = 600 (≈95 binary orbits) for different binary eccentricities and initial misalignment angles as indicated. The projections are along the intermediate binary axis $\boldsymbol{k}$ with the angular momentum vector $\boldsymbol{h}$ pointing upwards and the eccentricity vector $\boldsymbol{e}$ to the right.
Figure 3.

Density rendering of the simulations after t = 600 (≈95 binary orbits) for different binary eccentricities and initial misalignment angles as indicated. The projections are along the intermediate binary axis |$\boldsymbol{k}$| with the angular momentum vector |$\boldsymbol{h}$| pointing upwards and the eccentricity vector |$\boldsymbol{e}$| to the right.

4.1 Polar alignment

Nixon et al. (2011b) showed that for a circular binary, where the binary induced precession is only around |$\hat{\boldsymbol{h}}$|⁠, the disc eventually co- or counter-aligns with the binary orbital plane depending on the disc angular momentum and its initial misalignment angle. Our analysis in Section 2 suggests that for an eccentric binary the precession will be around |$\hat{\boldsymbol{h}}$| or |$\hat{\boldsymbol{e}}$|⁠. In the latter case, dissipation results in polar alignment.

Fig. 4 shows the time evolution of the tilt profiles θ(R) for two simulations with initial θ = 60° but either e = 0 (top) or 0.9 (bottom). For the circular binary case, the inner part of the disc eventually co-aligns with the binary (θ → 0°), as expected. In contrast, for the highly eccentric binary we see the disc aligning in a polar configuration with respect to the binary angular momentum vector (θ → 90°). The variations in θ as function of both time and radius are caused by the precession around |$\hat{\boldsymbol{e}}$|⁠.

Evolution of tilt profiles for e = 0 (top) and 0.9 (bottom) discs with initial tilt θ = 60° at t = 0, 100, 200, 300, 400, 500, and 600 (see legend) in code units.
Figure 4.

Evolution of tilt profiles for e = 0 (top) and 0.9 (bottom) discs with initial tilt θ = 60° at t = 0, 100, 200, 300, 400, 500, and 600 (see legend) in code units.

This is more evident from Fig. 5, where we plot the orientation |$\hat{\boldsymbol{l}}$| of the angular momentum in annuli of the disc at t = 100 for four simulations with different binary eccentricity but identical initial disc orientation at θ = 60°. We see that the discs in our simulations closely follow the predicted precession paths especially in the outer parts of the disc. The inner parts of the disc, which have higher precession rates, dissipate faster and start to align with the |$\hat{\boldsymbol{h}}$| or |$\hat{\boldsymbol{e}}$| as expected.

Projections of angular momenta of the radially binned disc in our simulations compared to the analytical precession paths shown in Fig. 1. Solid circles represent seven radial bins of the disc ranging from R = 1 (light grey) to 8 (dark grey) at t = 100 for simulations of different eccentricity (as indicated) but the same initial tilt θ = 60°. Obviously, the inner disc precesses faster, nicely following the theoretical precession paths. The innermost parts of the discs around eccentric binaries start to align with the stable polar orientation.
Figure 5.

Projections of angular momenta of the radially binned disc in our simulations compared to the analytical precession paths shown in Fig. 1. Solid circles represent seven radial bins of the disc ranging from R = 1 (light grey) to 8 (dark grey) at t = 100 for simulations of different eccentricity (as indicated) but the same initial tilt θ = 60°. Obviously, the inner disc precesses faster, nicely following the theoretical precession paths. The innermost parts of the discs around eccentric binaries start to align with the stable polar orientation.

4.2 Violent ring interactions

Our simulations starting from discs misaligned to both the |$\hat{\boldsymbol{h}}$| and |$\hat{\boldsymbol{e}}$| show rather violent gas dynamics. The radially differential binary torque tears the disc and causes the formation of separate rings. These rings are mutually misaligned and start to interact with each other, presumably because they gained some eccentricity from interactions with the binary. The ring interactions cause partial cancellation of angular momentum and hence a significant increase in the accretion rate. This is identical to the picture reported by Nixon et al. (2013) for circumbinary discs around circular binaries.

However, for very high eccentricities we find disc tearing to be much more violent and lead to a different evolution from that for circular and low-eccentricity binaries. There are two reasons for this difference: first the precession rate increases with eccentricity; second, the low-angular-momentum gas resulting from the interactions and falling on to the binary will align to polar orientation rather than a prograde or retrograde orientation as in the case of a near-circular binary. This allows this highly eccentric low-angular-momentum gas to come very close to the binary without suffering a lot of accretion. This non-circular gas in the central zone interacts with the outer disc further increasing its orbital eccentricity, throwing more gas to the centre, and promoting more interaction. This run away effect is shown in Fig. 6 and can also be seen in the bottom left-hand panel of Fig. 3. Eventually, this process sends an increasing amount of gas plunging on to the binary on almost radial orbits that can reach very close to the binary, avoiding significant accretion, and receiving energy kicks from one of the binary components in a manner very similar to the slingshot mechanism. Some of this gas will get ejected producing outward, almost radial, streams that can act as a possible observational signature of a highly eccentric SMBH binary.

Density rendering (left-hand panels) and particle plots coloured by eccentricity magnitude (right-hand panels) of five snapshots of the e = 0.9, θ = 150° run at times (from top to bottom) t = 0, 200, 400, 600, and 800.
Figure 6.

Density rendering (left-hand panels) and particle plots coloured by eccentricity magnitude (right-hand panels) of five snapshots of the e = 0.9, θ = 150° run at times (from top to bottom) t = 0, 200, 400, 600, and 800.

Fig. 6 shows density rendering (left-hand panels) and particle plots coloured by eccentricity magnitude (right-hand panels) of five snapshots for the simulation with e = 0.9 and initially θ = 150° at times t = 0, 200, 400, 600, and 800. We can see that the amount of chaotic gas resulting from the ring interactions keeps increasing during the simulation. The outward streams of gas resulting from the slingshot mechanism are very clear.

4.3 Precession rate

In order to provide a quantitative comparison between the predictions of our analytical model in Section 2 and the results obtained from the simulations, we plot in Fig. 7 the analytical precession rate Θ derived from equation (7) for a disc with an initial misalignment of θ = 60° around binaries with four different eccentricities along with the equivalent precession computed from the simulation and averaged over 10 binary orbits starting from t = 50. We find that the simulations agree quite well with the predicted precession rate at radii ≳2.5R/a. For discs around eccentric binaries, we observe oscillations on binary orbital time-scales and a good agreement is only found when the precession rate is averaged over a few binary orbits. We note that only a modest agreement is to be expected since our model ignores dissipative effects, contributions from higher than quadrupole order, and orbital resonances.

Comparison between the disc precession rate measured from the simulation (solid) with initial θ = 60° and e = 0, 0.3, 0.6, and 0.9 and our analytical model (dotted) of equation (7). For each case we plot the dominant component of Θ, i.e. Θh for e = 0, 0.3 and Θe for e = 0.6, 0.9.
Figure 7.

Comparison between the disc precession rate measured from the simulation (solid) with initial θ = 60° and e = 0, 0.3, 0.6, and 0.9 and our analytical model (dotted) of equation (7). For each case we plot the dominant component of Θ, i.e. Θh for e = 0, 0.3 and Θe for e = 0.6, 0.9.

4.4 Non-zero initial disc twist angle

So far, all the results shown here are for twist angles ϕ = 0°, i.e. initially the disc line of nodes with the binary plane is the |$\hat{\boldsymbol{k}}$| direction: |$\hat{\boldsymbol{l}}$| is tilted towards |$\hat{\boldsymbol{e}}$|⁠. In general, however, we should expect any disc orientation, i.e. non-zero ϕ.

In Fig. 8 we present snapshots for simulations with binary eccentricity e = 0.9, initial twist angles ϕ = 0°, 45°, and 90°, and initial tilt angles θ = 30°, 45°, 60°, 80°, and 90°. For ϕ = 90°, we only observe azimuthal precession, akin to the circular binary-induced precession. This confirms our prediction since for that case |$\hat{\boldsymbol{l}}$| and |$\hat{\boldsymbol{e}}$| are always orthogonal, causing the first term in equation (7) to vanish, and we are left with only azimuthal precession. For ϕ = 45°, we find the same trend discussed earlier, i.e. experiencing either polar or azimuthal precession, or violent ring interaction. In Fig. 9 we compare the precession paths of all three ϕ values for the case of e = 0.9 and θ = 60° to the analytical contours. Similar to Fig. 5, we see the simulations closely follow the analytical contours apart from the innermost parts where disc breaking and alignment are dominant. This strongly suggests that our analysis above for ϕ = 0° carries over to the general case.

Density rendering of the e = 0.9 simulations at t = 600 (≈95 orbits of the binary) projected on the x–z plane with different values for ϕ and θ as indicated in the figure.
Figure 8.

Density rendering of the e = 0.9 simulations at t = 600 (≈95 orbits of the binary) projected on the xz plane with different values for ϕ and θ as indicated in the figure.

Precession paths for e = 0.9 and θ = 60° at t = 100 with three different ϕ values: ϕ = 0° (circles), 45° (rectangles), and 90° (triangles).
Figure 9.

Precession paths for e = 0.9 and θ = 60° at t = 100 with three different ϕ values: ϕ = 0° (circles), 45° (rectangles), and 90° (triangles).

5 IMPLICATIONS FOR THE FINAL-PARSEC PROBLEM

The solution suggested by Nixon et al. (2011a) to the final-parsec problem requires the binary to accrete negative angular momentum from a retrograde disc, which gradually increases the binary eccentricity until coalescence is achieved via energy losses to gravitational radiation at pericentre. Nixon et al. (2011b) showed that for a circular binary counter-alignment of randomly oriented accretion events can provide a continuous supply of the required retrograde discs. In particular, they showed that for cases where Jb > 2Jd, all accretion events with initial misalignment of θ > 90° will result in a retrograde disc, i.e. roughly half of randomly oriented accretion events will counter-align with the binary as long as the binary dominates the angular momentum of the system.

Our results somewhat change this picture. As the binary eccentricity increases (due to retrograde accretion as in Nixon et al. or earlier stellar dynamical processes) disc counter- (and co-) alignment becomes ever less likely at the expense of polar alignment. The subsequent accretion of such polar discs merely rotates the angular momentum vector of the binary presumably hardly affecting the binary eccentricity. Thus simply retrograde gas accretion appears less viable a solution to the final parse problem.

There are, however, still several ways gas can solve this problem. First, a single massive retrograde accretion event may, in principle, supply enough negative angular momentum to complete the binary merger. However, for a massive disc self-gravity becomes important, likely causing clumping and star formation, which reduces the amount of gas that can be accreted. Moreover, a single massive retrograde accretion event may be not be sufficiently likely to explain the coalescence of all SMBH binaries (which form with each major merger of massive galaxies).

A more intriguing possibility involves more violent gas dynamics. We showed that, in many cases, the disc does not smoothly align, instead the strong differential precession (in particular for misaligned discs around eccentric binaries) leads to tearing of the disc into separate mutually misaligned rings. In the inner disc close to the binary, the gravity of the binary causes these rings to become eccentric such that they inevitably interact with each other and with the outer disc. These interactions cause further eccentricity growth on a dynamical time-scale and eventually result in plunging gas infall. Some of this infalling gas will be accreted by either binary component. This will change the binary angular momentum, but may not reduce its absolute value, depending on the orientation and in contrast to the situation with predominantly retrograde accretion (Nixon et al. 2013).

If the infalling gas evades this fate, it will most likely get ejected from the binary via a slingshot interaction. This in turn reduces the binary separation in much the same way as the ejection of penetrating stars, thus exactly as required to solve the final-parsec problem. Indeed, we find in our simulations which undergo violent gas dynamics not only significant gas accretion but also a shrinking of the binary orbit.

Clearly, this violent interaction and accretion processes are rather complex and chaotic and certainly not well resolved or adequately modelled in our simulations. Nonetheless, what our simulations quite clearly show is that such violent gas dynamical processes are inevitable if the gas is initially misaligned with the binary, in particular if the binary is eccentric. We leave a more detailed investigation into the binary orbital evolution in this chaotic environment for a future study.

6 CONCLUSIONS

We have studied the interaction of an eccentric binary with a gaseous disc initially misaligned with the binary angular momentum. Such a configuration should occur naturally from the infall and subsequent circularization of gas into the inner few parsec of a merger remnant still hosting a SMBH binary (e.g. Dunhill et al. 2014). The binary exerts a torque on the disc, resulting in disc precession and, due to viscous dissipation, in eventual alignment of the disc with the binary. In the case of a circular binary, this alignment is always coplanar, resulting either in a progradely or retrogradely rotating circumbinary disc (Nixon et al. 2011b).

We find that in the general case of an eccentric binary polar alignment also occurs, when disc angular momentum is aligned with the binary periapse or apoapse direction. The binary torque on the disc can be quite accurately understood analytically assuming an orbit-averaged binary potential to quadrupole order (see Section 2). The fraction of initial disc orientations which give rise to polar alignment grows with binary eccentricity, reaching 0.5 at e ≈ 0.4. The precession paths (neglecting dissipation) are not circular, but elongated towards the intermediate axis of the orbit-averaged binary.

The prospect of accretion on to the binary from a polar instead of coplanar disc impedes the solution (proposed by Nixon et al. 2011a) to the final-parsec problem for coalescing a SMBH binary. In that picture consecutive randomly oriented accretion events lead to the formation of either prograde or retrograde coplanar circumbinary discs. While accretion from the former is largely suppressed (by orbital resonances as discussed in the Introduction), accretion from the latter reduces the binary angular momentum and drives it to larger eccentricities. However, at large eccentricities polar disc orientations dominate, when accretion (not resolved in our simulations) has presumably little effect on the binary orbit (since the accreted angular momentum is perpendicular to that of the binary). Thus, eccentricity growth via accretion is likely to be significantly reduced well before gravitational wave emission can take over as driver for coalescence.

However, in many of our simulations, in particular for larger binary eccentricity and stronger initial misalignment, the disc does not smoothly align, but is torn into separate mutually misaligned rings. This process was already reported by Nixon et al. (2013) for a circular binary and can be understood by the radially differential binary torque, which overcomes the adhesive effect of gas viscosity. The prominence of tearing with binary eccentricity and initial disc misalignment can be understood as consequence of the larger binary torque in these cases.

The subsequent evolution of these gas rings can be rather chaotic and is not quite adequately modelled in our simulations. However, some basic results appear to be robust. The innermost rings are sufficiently perturbed by the binary to acquire some orbital eccentricity. This in turn inevitably leads to interactions between the rings, resulting in partial cancellation of their angular momenta. This process is more prominent in more eccentric binaries, because the stronger binary torque results in larger mutual misalignment between adjacent rings. The cancellation of angular momentum of the rings will increase their eccentricity, providing a positive feedback loop and hence a runaway process, eventually resulting in gas plunging on to the central binary. This material may be accreted on to either hole, but when coming from a near-polar orientation, this will hardly help with the final-parsec problem, as explained above.

Alternatively, the infalling gas, which for a highly eccentric binary can come much closer to the binary whilst avoiding accretion, may get ejected from the binary via a gravitational slingshot interaction with one of its components. This also helps to solve the final-parsec problem, though this time by reducing its semimajor axis. This is similar to the stellar dynamical process of shrinking the binary orbit via ejection of stars penetrating into the binary. The difference is that the total amount of stars in the ‘loss cone’ (whose orbit carries them into inner parsec) is limited and cannot be easily refilled, while gas being dissipative and collisional by nature may provide a better agent. This is particularly the case at the parsec scale where the SMBH dominates the dynamics and by its gravitational torques shepherds some gas into the loss cone.

We thank the referee for useful suggestions. Research in Theoretical Astrophysics at Leicester is supported by an STFC rolling grant. The calculations presented in this paper were performed using the ALICE High Performance Computing Facility at the University of Leicester. Some resources on ALICE form part of the DiRAC Facility jointly funded by STFC and the Large Facilities Capital Fund of BIS. CN acknowledges support provided by NASA through the Einstein Fellowship Program, grant PF2-130098. We used splash (Price 2007) visualization code for some of the figures in this paper.

Einstein Fellow.

1

This effect was present in the simulations of Cuadra et al. (2009) but did not effectuate significant binary evolution. On the other hand, based on an extrapolation to 50 times longer than actually modelled, Roedig & Sesana (2014) claim efficient binary shrinking. However, since the infall of gas depends on the disc structure at its inner edge, this result is very sensitive to the thermodynamical treatment. Roedig et al. (2012), for example, found that for isothermal instead of adiabatic gas with an imposed standard β cooling prescription, the binary orbital decay can be significantly reduced.

REFERENCES

Aarseth
S. J.
Ap&SS
,
2003
, vol.
285
pg.
367
Bardeen
J. M.
Petterson
J. A.
ApJ
,
1975
, vol.
195
pg.
L65
Begelman
M. C.
Blandford
R. D.
Rees
M. J.
Nature
,
1980
, vol.
287
pg.
307
Berczik
P.
Merritt
D.
Spurzem
R.
ApJ
,
2005
, vol.
633
pg.
680
Berczik
P.
Merritt
D.
Spurzem
R.
Bischof
H.-P.
ApJ
,
2006
, vol.
642
pg.
L21
Cuadra
J.
Armitage
P. J.
Alexander
R. D.
Begelman
M. C.
MNRAS
,
2009
, vol.
393
pg.
1423
Dehnen
W.
Aly
H.
MNRAS
,
2012
, vol.
425
pg.
1068
Dunhill
A. C.
Alexander
R. D.
Nixon
C. J.
King
A. R.
MNRAS
,
2014
, vol.
445
pg.
2285
Escala
A.
Larson
R. B.
Coppi
P. S.
Mardones
D.
ApJ
,
2005
, vol.
630
pg.
152
Gammie
C. F.
ApJ
,
2001
, vol.
553
pg.
174
Gingold
R. A.
Monaghan
J. J.
MNRAS
,
1977
, vol.
181
pg.
375
Goldreich
P.
Tremaine
S.
ApJ
,
1979
, vol.
233
pg.
857
Goodman
J.
MNRAS
,
2003
, vol.
339
pg.
937
Gualandris
A.
Merritt
D.
ApJ
,
2012
, vol.
744
pg.
74
Khan
F. M.
Just
A.
Merritt
D.
ApJ
,
2011
, vol.
732
pg.
89
Khan
F. M.
Holley-Bockelmann
K.
Berczik
P.
Just
A.
ApJ
,
2013
, vol.
773
pg.
100
King
A. R.
Pringle
J. E.
MNRAS
,
2006
, vol.
373
pg.
L90
King
A. R.
Pringle
J. E.
MNRAS
,
2007
, vol.
377
pg.
L25
King
A. R.
Lubow
S. H.
Ogilvie
G. I.
Pringle
J. E.
MNRAS
,
2005
, vol.
363
pg.
49
King
A. R.
Pringle
J. E.
Hofmann
J. A.
MNRAS
,
2008
, vol.
385
pg.
1621
Levin
Y.
MNRAS
,
2007
, vol.
374
pg.
515
Lodato
G.
Price
D. J.
MNRAS
,
2010
, vol.
405
pg.
1212
Lodato
G.
Pringle
J. E.
MNRAS
,
2007
, vol.
381
pg.
1287
Lodato
G.
Nayakshin
S.
King
A. R.
Pringle
J. E.
MNRAS
,
2009
, vol.
398
pg.
1392
Lucy
L. B.
AJ
,
1977
, vol.
82
pg.
1013
MacFadyen
A. I.
Milosavljević
M.
ApJ
,
2008
, vol.
672
pg.
83
Merritt
D.
ApJ
,
2001
, vol.
556
pg.
245
Milosavljević
M.
Merritt
D.
Centrella
J. M.
AIP Conf. Ser. Vol. 686, The Astrophysics of Gravitational Wave Sources
,
2003
New York
Am. Inst. Phys.
pg.
201
Naoz
S.
Farr
W. M.
Lithwick
Y.
Rasio
F. A.
Teyssandier
J.
MNRAS
,
2013
, vol.
431
pg.
2155
Nixon
C. J.
MNRAS
,
2012
, vol.
423
pg.
2597
Nixon
C. J.
Cossins
P. J.
King
A. R.
Pringle
J. E.
MNRAS
,
2011a
, vol.
412
pg.
1591
Nixon
C. J.
King
A. R.
Pringle
J. E.
MNRAS
,
2011b
, vol.
417
pg.
L66
Nixon
C.
King
A.
Price
D.
Frank
J.
ApJ
,
2012
, vol.
757
pg.
L24
Nixon
C.
King
A.
Price
D.
MNRAS
,
2013
, vol.
434
pg.
1946
Price
D. J.
Publ. Astron. Soc. Aust.
,
2007
, vol.
24
pg.
159
Pringle
J. E.
MNRAS
,
1992
, vol.
258
pg.
811
Roedig
C.
Sesana
A.
MNRAS
,
2014
, vol.
439
pg.
3476
Roedig
C.
Sesana
A.
Dotti
M.
Cuadra
J.
Amaro-Seoane
P.
Haardt
F.
A&A
,
2012
, vol.
545
pg.
A127
Saslaw
W. C.
Valtonen
M. J.
Aarseth
S. J.
ApJ
,
1974
, vol.
190
pg.
253
Scheuer
P. A. G.
Feiler
R.
MNRAS
,
1996
, vol.
282
pg.
291
Sesana
A.
Haardt
F.
Madau
P.
ApJ
,
2007
, vol.
660
pg.
546
Sesana
A.
Haardt
F.
Madau
P.
ApJ
,
2008
, vol.
686
pg.
432
Shakura
N. I.
Sunyaev
R. A.
A&A
,
1973
, vol.
24
pg.
337
Steiman-Cameron
T. Y.
Durisen
R. H.
ApJ
,
1984
, vol.
276
pg.
101
Thomas
P. A.
Vine
S.
Pearce
F. R.
MNRAS
,
1994
, vol.
268
pg.
253
Wang
L.
Berczik
P.
Spurzem
R.
Kouwenhoven
M. B. N.
ApJ
,
2014
, vol.
780
pg.
164

APPENDIX A: BINARY–DISC QUADRUPOLE INTERACTION

Here, we give the details of the analysis leading to the results reported in Section 2. Our results, obtained via Newtonian dynamics, agree with the more general results of Naoz et al. (2013), obtained via Hamiltonian perturbation theory, for a circular ring and ignoring octopole terms.

The three unit vectors |$\hat{\boldsymbol{h}}$|⁠, |$\hat{\boldsymbol{e}}$|⁠, and |$\hat{\boldsymbol{k}}=\hat{\boldsymbol{h}}\times \hat{\boldsymbol{e}}$| are conserved under the binary motion and are mutually orthogonal, such that
(A1)
The binary components are at positions
(A2)
with
(A3)
Here, η is the eccentric anomaly, which is related to the mean anomaly ℓ via
(A4)
such that dℓ = (1 − ecos η)dη and an orbit average becomes |$\langle \cdot \rangle =(2\pi )^{-1}\int _0^{2\pi }\cdot \,(1-e\cos \eta )\mathrm{d}\eta$|⁠. When orbit-averaging RiRj, the cross term between |$\hat{\boldsymbol{e}}$| and |$\hat{\boldsymbol{k}}$| averages to zero and
(A5)
(A6)
where the second form follows from eliminating |$\hat{k}_i\hat{k}_{\!j}$| in favour of |$\hat{h}_i\hat{h}_{\!j}$| with the help of the identity (A1). From this result, we can work out the orbit-averaged trace-free specific quadrupole moment of the binary as
(A7)
(A8)
(A9)
We will also need the orbit average
(A10)
with |$\Omega =\sqrt{GM/a^3}$| the binary orbital frequency.

A1 Ring evolution

A ring particle at position |$\boldsymbol{r}$| experiences the orbit-averaged quadrupole potential of the binary:
(A11)
Averaging over the ring, we obtain
(A12)
such that the trace-free specific quadrupole moment of the ring is
(A13)
The quadrupole interaction energy between binary and ring, averaged over the binary orbit and the ring, is then
(A14)
with interaction tensor |$\boldsymbol{\mathsf {\Theta }}\equiv \boldsymbol{\mathsf {q}}\cdot \boldsymbol{\mathsf {Q}}$|⁠, which has components
(A15)
Taking its trace in equation (A14) we obtain equation (5). The orbit-averaged torque on the ring also involves the interaction tensor. Using index notation, we have
(A16)
Note that only the antisymmetric part of |$\boldsymbol{\mathsf {\Theta }}$| contributes to the torque. Inserting equation (A15), we find that we can write this as |$\dot{\boldsymbol{l}} = \boldsymbol{\Theta } \times \boldsymbol{l}$| with the vector
(A17)
where |$\omega =\sqrt{G(M+m)/r^3}$| is the ring angular frequency.

A2 Binary evolution

The torque of the binary from the ring can be worked out analogously to that of the ring from the binary. The quadrupole potential due to the ring at the binary is
(A18)
Adding the torque from each binary component and averaging over the orbit, we find
(A19)
In particular, the total angular momentum, |$M\boldsymbol{h}+m\boldsymbol{l}$|⁠, is conserved at quadrupole order. Together with the precession of the ring, this implies that the evolution of |$\boldsymbol{h}$| is not simply a precession and that |$h=|\boldsymbol{h}|$| is in general not conserved. Instead, we find
(A20)
Thus, h remains unchanged only for a circular binary. |$\dot{h}=0$| for e = 0 or if |$\hat{\boldsymbol{l}}$| is perpendicular to either |$\boldsymbol{e}$| or |$\hat{\boldsymbol{k}}$|⁠. Otherwise, h oscillates, since |$\hat{\boldsymbol{l}}\cdot \hat{\boldsymbol{k}}$| oscillates around zero under the ring precession.
The change of the eccentricity vector is
(A21)
with
(A22)
Inserting (A22) into (A21) and orbit averaging (using equation A10), we find
(A23)
(A24)
(A25)
From this, we obtain
(A26)
in agreement with equation (A34) of Naoz et al. (2013), but also with equation (A20) and |$(1-e^2)\dot{h}=-he\dot{e}$| (from equation 4).

A3 Ring precession

If the mass of the ring is negligible compared to that of the binary, we can approximate the binary orientation as fixed and the vectors |$\hat{h}$|⁠, |$\hat{e}$|⁠, and |$\hat{k}$| as constants. In this case, the evolution of the ring orientation allows some further analytical treatment.

Since |$\boldsymbol{\Theta }$| is parallel to |$\mathrm{\partial} {\langle {E_{\mathrm{br}}}\rangle }/\mathrm{\partial} {\hat{\boldsymbol{l}}}$|⁠, precession is along the lines of constant 〈Ebr〉. This gives the equation
(A27)
with constant C for the precession paths. C = 0 corresponds to the contour of 〈Ebr〉 through the unstable orientations |$\hat{\boldsymbol{l}}=\pm \hat{\boldsymbol{k}}$|⁠. Hence, this contour separates the regions of polar and azimuthal precession. The right-hand side (rhs) of equation (A27) can be written as |$(\hat{\boldsymbol{l}}\cdot \boldsymbol{u}_{1})\;(\hat{\boldsymbol{l}}\cdot \boldsymbol{u}_{2})$| with
(A28)
Thus, the separatrices are great circles with poles |$\boldsymbol{u}_{1,2}$|⁠. The fraction of ring orientations undergoing polar precession is
(A29)
At small e, this grows linearly (⁠|${\propto }\sqrt{20}e/\pi$|⁠) with eccentricity. Azimuthal and polar precession are equally likely for |$\hat{\boldsymbol{u}}_1\cdot \hat{\boldsymbol{u}}_2=0$|⁠, which occurs at e = 6−1/2 ≈ 0.408.
If the constant C in equation (A27) is positive (negative), we have azimuthal (polar) precession. This equation for |$\hat{\boldsymbol{l}}$| has the parametric solutions for the precession paths (see also Fig. 1):
(A30)
for 0 < C ≤ 1 − e2 (azimuthal precession), and
(A31)
for −5e2 ≤ C < 0 (polar precession). In either case, the third component of |$\hat{\boldsymbol{l}}$| follows from the normalization condition |$|\hat{\boldsymbol{l}}|=1$|⁠.
The instantaneous precession rate |$|\mathrm{d}\hat{\boldsymbol{l}}/\mathrm{d}t|=|\boldsymbol{\Theta }\times \hat{\boldsymbol{l}}|$| varies along the precession paths. It is minimal at the largest value of |$|\hat{\boldsymbol{l}}\cdot \hat{\boldsymbol{k}}|$| along the precession path and drops to zero for |$\hat{\boldsymbol{l}}=\hat{\boldsymbol{k}}$| (the unstable orientations). The instantaneous precession rate becomes maximal at |$\hat{\boldsymbol{l}}\cdot \hat{\boldsymbol{k}}=0$| (at zero twist ϕ), when
(A32)

Thus, the maximum precession rate is much larger for highly eccentric than for circular binaries. The largest variation of the precession rate occurs for precession paths close to the separatrices.