https://orcid.org/0000-0002-5187-7561
ABSTRACT
A test particle orbit around an eccentric binary has two stationary states in which there is no nodal precession: coplanar and polar. Nodal precession of a misaligned test particle orbit centres on one of these stationary states. A low-mass circumbinary disc undergoes the same precession and moves towards one of these states through dissipation within the disc. For a massive particle orbit, the stationary polar alignment occurs at an inclination less than 90°, which is the prograde-polar stationary inclination. A sufficiently high angular momentum particle has an additional higher inclination stationary state, the retrograde-polar stationary inclination. Misaligned particle orbits close to the retrograde-polar stationary inclination are not nested like the orbits close to the other stationary points. We investigate the evolution of a gas disc that begins close to the retrograde-polar stationary inclination. With hydrodynamical disc simulations, we find that the disc moves through the unnested crescent shape precession orbits and eventually moves towards the prograde-polar stationary inclination, thus increasing the parameter space over which circumbinary discs move towards polar alignment. If protoplanetary discs form with an isotropic orientation relative to the binary orbit, then polar discs may be more common than coplanar discs around eccentric binaries, even for massive discs. This has implications for the alignment of circumbinary planets.
1 INTRODUCTION
Binary star systems are commonly observed (e.g. Ghez, Neugebauer & Matthews 1993; Duchêne & Kraus 2013). They are known to form in turbulent molecular clouds (McKee & Ostriker 2007), where the accretion process may be chaotic (Bate, Bonnell & Bromm 2003; Bate 2018). This results in a high likelihood that the protoplanetary disc that forms is misaligned with respect to the orbital plane of the binary (Monin et al. 2007; Bate, Lodato & Pringle 2010; Bate 2018). Circumbinary disc observations suggest that misalignments are common (e.g. Chiang & Murray-Clay 2004; Winn et al. 2004; Capelo et al. 2012; Kennedy et al. 2012; Brinch et al. 2016; Aly, Lodato & Cazzoletti 2018). Circumbinary discs have been observed with polar inclinations of 90° to the binary orbit. HD98800 and V773 Tau B have polar circumbinary gas discs (Kennedy et al. 2019; Kenworthy et al. 2022), while 99 Her is a polar circumbinary debris disc (Kennedy et al. 2012; Smallwood et al. 2020). The formation and evolution of planets around binaries may be altered by the torque from the binary (e.g. Nelson 2000; Mayer et al. 2005; Boss 2006; Martin et al. 2014; Fu, Lubow & Martin 2015a, b, 2017; Franchini, Martin & Lubow 2019). Giant planets that form in a misaligned disc may not remain coplanar to the disc (Picogna & Marzari 2015; Lubow & Martin 2016; Martin et al. 2016). We seek to understand the evolution of misaligned circumbinary discs in order to better understand the observed properties of exoplanets.
The Kepler Mission has so far detected about 20 circumbinary planets (e.g. Doyle et al. 2011; Orosz et al. 2012; Welsh et al. 2012; Kostov et al. 2016) and the TESS Mission recently discovered one more (Kostov et al. 2021). All of the circumbinary planets detected so far are nearly coplanar to the binary orbital plane. However, this may be a selection effect due to the small orbital period of the Kepler binaries (Czekala et al. 2019). Longer orbital period binaries are expected to host planets with a wide range of inclinations. Terrestrial circumbinary planets that form through core accretion in the absence of gas may only form in polar or coplanar configurations (Childs & Martin 2021a, b, 2022), but gas giants that form within the gas disc may form with the same initial inclination as the gas disc (e.g. Pierens & Nelson 2018). Planets with large misalignments are much more difficult to detect than those that orbit in the binary orbital plane because their transits are much rarer (Schneider 1994; Martin & Triaud 2014, 2015; Martin 2017; Chen & Kipping 2022). However, other detection methods are possible, such as eclipse timing variations (Murat Esmer et al. 2022). Polar planets (those in an orbit close to perpendicular to the binary orbital plane) may be distinguished from coplanar planets through eclipse timing variations of the binary (Zhang & Fabrycky 2019). Two circumbinary planets have been found around eccentric orbit binaries. Kepler-34b has a mass of 0.22MJ and orbits an eclipsing binary star system (Kepler-34) that has an orbital eccentricity of 0.52 (Welsh et al. 2012; Kley & Haghighipour 2015).
The evolution of low-mass circumbinary discs is fairly well understood. A misaligned circumbinary disc around a circular orbit binary experiences uniform nodal precession with constant tilt in the absence of dissipation (e.g. Nixon 2012; Facchini, Lodato & Price 2013; Foucart & Lai 2013; Lodato & Facchini 2013). The angular momentum vector of the disc precesses about the angular momentum vector of the binary. We call this a circulating solution. In this work, we focus on protoplanetary discs that are in the wave-like regime in which the disc aspect ratio is much larger than the Shakura & Sunyaev (1973a) viscosity parameter, H/R > α. For a sufficiently warm and radially narrow disc, the disc precesses as a solid body (e.g. Larwood & Papaloizou 1997). In these systems, dissipation in the disc leads to alignment with the binary orbital plane (Papaloizou & Terquem 1995; Lubow & Ogilvie 2000; Nixon, King & Pringle 2011; Nixon 2012; Facchini, Lodato & Price 2013; Foucart & Lai 2013, 2014; Lodato & Facchini 2013).
The evolution around an eccentric orbit binary is more complex. For a sufficiently high initial tilt, a low-mass disc precesses about the eccentricity vector of the binary (rather than the orbital rotation axis of the binary; Verrier & Evans 2009; Farago & Laskar 2010; Doolin & Blundell 2011; Aly et al. 2015). We call this a librating solution, where the orbit is precessing about the prograde-polar stationary state. In the low-mass disc case, the prograde-polar stationary state is aligned with the binary eccentricity vector. Dissipation leads to polar alignment where the circumbinary disc is aligned to the eccentricity vector of the binary and perpendicular to the binary orbital plane (Martin & Lubow 2017; Lubow & Martin 2018; Martin & Lubow 2018; Zanazzi & Lai 2018; Cuello & Giuppone 2019).
The evolution of a massive circumbinary gas disc around an eccentric binary has not yet been fully explored. In the absence of dissipation within the disc, the dynamical evolution of a radially narrow circumbinary gas disc is qualitatively similar to that of a massive particle orbit around the binary. Thus, in order to understand the evolution of a massive disc, in Section 2, we first explore the dynamics of massive particle orbits (see also Chen et al. 2019). In this case, the binary feels the gravitational force of the planet, which causes the binary orbit to evolve. The eccentricity vector of the binary precesses, the binary orbit tilts, and the magnitude of its eccentricity oscillates. The prograde-polar stationary aligned state occurs at a lower level of misalignment (Farago & Laskar 2010; Zanazzi & Lai 2018; Chen et al. 2019; Martin & Lubow 2019). For librating orbits, the angular momentum vector precesses around the prograde-polar stationary state. For a gas disc, dissipation causes these oscillations to damp and the disc settles at the prograde-polar stationary tilt angle.
In the massive particle case, there is a third type of nodal precession orbit, which we call the crescent orbits (Martin & Lubow 2019; Chen et al. 2019; see Section 2.3 for more details). These occur for sufficiently massive particles that begin closer to retrograde than to prograde. These crescent orbits are not nested like the circulating and librating orbits, meaning that they do not have a common centre for the precession. For a sufficiently high angular momentum particle, there is an additional retrograde-polar stationary inclination in the region of these crescent orbits. The particle precession in a crescent orbit is not centred around a stationary state. As we show in this paper, a gas disc on such orbits often evolves towards the prograde-polar stationary tilt angle. In Section 3, we consider the long-term evolution of a gas disc that begins in this regime, close to the retrograde-polar stationary inclination. Finally, in Section 4, we discuss the implications of our results and state our conclusions.
2 THREE-BODY SIMULATIONS
In this section, before running our circumbinary disc simulations, we first consider the evolution of a three-body system in order to understand the dynamics of massive circumbinary particle orbits. The parameters we vary are the initial inclination of the third body/planet orbit with respect to the inner binary orbital plane, i0, and the angular momentum ratio of the third body compared to the inner binary, j. Note that in this paper, j refers to the angular momentum ratio between both planet-binary and disc-binary for the three-body systems and circumbinary disc systems, respectively. For the three-body simulations, we follow the methods outlined in Chen et al. (2019). In Section 3, we compare these results with circumbinary disc simulations that have roughly the same initial parameters, j and i0.
2.1 Three-body simulation set-up
We use the N-body simulation package, rebound, with the WHfast integrator that is a second-order symplectic Wisdom Holman integrator with 11th order symplectic correctors (Rein & Tamayo 2015). We solve the gravitational equations for a planet orbiting around a binary star system. The primary mass in the binary is set to be |$M_{1} = 0.9M$|, while the secondary mass is |$M_{2} = 0.1M$|, where M = M1 + M2 is the total mass of the binary. The initial eccentricity of the binary orbit is set to eb = 0.8, and the initial semimajor axis of the binary is ab.
This table lists the simulation parameters and outcomes. Each line consists of both an N-body planet simulation and an SPH disc simulation with the same combination of j and i0. The first column describes the simulation name. The second column is the initial inclination of the planet/disc. The third column is the angular momentum ratio of the planet/disc to the binary. The fourth column is the analytically calculated retrograde-polar stationary inclination (see equation 5). The fifth column is the mass of the planet. The sixth column describes whether the orbit of the planet is retrograde circulating (CR), librating (L) or in a crescent orbit (CRESCENT). The seventh column is the mass of the disc. The eighth column describes the orbit type of the disc prior to its 30 per cent mass-loss cut-off. The ninth column describes the disc break radius. The 10th column is the time at which the circumbinary disc simulation has lost 30 per cent of its initial mass.
| Name | i0 (°) | j | irs(°) | Mp/M | Planet orbit | Md/M | Disc orbit | Break radius/ab | Time (Tb) |
|---|---|---|---|---|---|---|---|---|---|
| Low-j-120 | 120 | 0.5 | – | 0.006 | L | 0.01 | L | – | 2110 |
| Low-j-130 | 130 | 0.5 | – | 0.006 | L | 0.01 | L | – | 1540 |
| Low-j-140 | 140 | 0.5 | – | 0.006 | L | 0.01 | L | – | 1070 |
| Low-j-150 | 150 | 0.5 | – | 0.006 | L | 0.01 | L | – | 850 |
| Low-j-160 | 160 | 0.5 | – | 0.006 | L | 0.01 | L | 6 | 1560 |
| Low-j-170 | 170 | 0.5 | – | 0.006 | L | 0.01 | L | – | 2520 |
| Mid-j-100 | 100 | 1.5 | 139 | 0.018 | L | 0.03 | L | – | 2880 |
| Mid-j-110 | 110 | 1.5 | 139 | 0.018 | CRESCENT | 0.03 | L | – | 2320 |
| Mid-j-120 | 120 | 1.5 | 139 | 0.018 | CRESCENT | 0.03 | CRESCENT → L | – | 1830 |
| Mid-j-130 | 130 | 1.5 | 139 | 0.018 | CRESCENT | 0.03 | CRESCENT | – | 1410 |
| Mid-j-140 | 140 | 1.5 | 139 | 0.018 | irs | 0.03 | |$i_{\rm rs}\, \rightarrow$| L | – | 1140 |
| Mid-j-150 | 150 | 1.5 | 139 | 0.018 | CRESCENT | 0.03 | CRESCENT (?) | – | 870 |
| Mid-j-160 | 160 | 1.5 | 139 | 0.018 | CR | 0.03 | CR | 2.5 | 1000 |
| Mid-j-170 | 170 | 1.5 | 139 | 0.018 | CR | 0.03 | CR | – | 2620 |
| High-j-100 | 100 | 2.5 | 129 | 0.031 | L | 0.05 | L | – | 2380 |
| High-j-110 | 110 | 2.5 | 129 | 0.031 | CRESCENT | 0.05 | CRESCENT → L | – | 2310 |
| High-j-120 | 120 | 2.5 | 129 | 0.031 | CRESCENT | 0.05 | CRESCENT | – | 2230 |
| High-j-130 | 130 | 2.5 | 129 | 0.031 | irs | 0.05 | |$i_{\rm rs}\, \rightarrow$| L | – | 4200 |
| High-j-140 | 140 | 2.5 | 129 | 0.031 | CRESCENT | 0.05 | CRESCENT → L | – | 1980 |
| High-j-150 | 150 | 2.5 | 129 | 0.031 | CR | 0.05 | CR | – | 2020 |
| High-j-160 | 160 | 2.5 | 129 | 0.031 | CR | 0.05 | CR | 2.5 | 1110 |
| High-j-170 | 170 | 2.5 | 129 | 0.031 | CR | 0.05 | CR | – | 2620 |
| Name | i0 (°) | j | irs(°) | Mp/M | Planet orbit | Md/M | Disc orbit | Break radius/ab | Time (Tb) |
|---|---|---|---|---|---|---|---|---|---|
| Low-j-120 | 120 | 0.5 | – | 0.006 | L | 0.01 | L | – | 2110 |
| Low-j-130 | 130 | 0.5 | – | 0.006 | L | 0.01 | L | – | 1540 |
| Low-j-140 | 140 | 0.5 | – | 0.006 | L | 0.01 | L | – | 1070 |
| Low-j-150 | 150 | 0.5 | – | 0.006 | L | 0.01 | L | – | 850 |
| Low-j-160 | 160 | 0.5 | – | 0.006 | L | 0.01 | L | 6 | 1560 |
| Low-j-170 | 170 | 0.5 | – | 0.006 | L | 0.01 | L | – | 2520 |
| Mid-j-100 | 100 | 1.5 | 139 | 0.018 | L | 0.03 | L | – | 2880 |
| Mid-j-110 | 110 | 1.5 | 139 | 0.018 | CRESCENT | 0.03 | L | – | 2320 |
| Mid-j-120 | 120 | 1.5 | 139 | 0.018 | CRESCENT | 0.03 | CRESCENT → L | – | 1830 |
| Mid-j-130 | 130 | 1.5 | 139 | 0.018 | CRESCENT | 0.03 | CRESCENT | – | 1410 |
| Mid-j-140 | 140 | 1.5 | 139 | 0.018 | irs | 0.03 | |$i_{\rm rs}\, \rightarrow$| L | – | 1140 |
| Mid-j-150 | 150 | 1.5 | 139 | 0.018 | CRESCENT | 0.03 | CRESCENT (?) | – | 870 |
| Mid-j-160 | 160 | 1.5 | 139 | 0.018 | CR | 0.03 | CR | 2.5 | 1000 |
| Mid-j-170 | 170 | 1.5 | 139 | 0.018 | CR | 0.03 | CR | – | 2620 |
| High-j-100 | 100 | 2.5 | 129 | 0.031 | L | 0.05 | L | – | 2380 |
| High-j-110 | 110 | 2.5 | 129 | 0.031 | CRESCENT | 0.05 | CRESCENT → L | – | 2310 |
| High-j-120 | 120 | 2.5 | 129 | 0.031 | CRESCENT | 0.05 | CRESCENT | – | 2230 |
| High-j-130 | 130 | 2.5 | 129 | 0.031 | irs | 0.05 | |$i_{\rm rs}\, \rightarrow$| L | – | 4200 |
| High-j-140 | 140 | 2.5 | 129 | 0.031 | CRESCENT | 0.05 | CRESCENT → L | – | 1980 |
| High-j-150 | 150 | 2.5 | 129 | 0.031 | CR | 0.05 | CR | – | 2020 |
| High-j-160 | 160 | 2.5 | 129 | 0.031 | CR | 0.05 | CR | 2.5 | 1110 |
| High-j-170 | 170 | 2.5 | 129 | 0.031 | CR | 0.05 | CR | – | 2620 |
This table lists the simulation parameters and outcomes. Each line consists of both an N-body planet simulation and an SPH disc simulation with the same combination of j and i0. The first column describes the simulation name. The second column is the initial inclination of the planet/disc. The third column is the angular momentum ratio of the planet/disc to the binary. The fourth column is the analytically calculated retrograde-polar stationary inclination (see equation 5). The fifth column is the mass of the planet. The sixth column describes whether the orbit of the planet is retrograde circulating (CR), librating (L) or in a crescent orbit (CRESCENT). The seventh column is the mass of the disc. The eighth column describes the orbit type of the disc prior to its 30 per cent mass-loss cut-off. The ninth column describes the disc break radius. The 10th column is the time at which the circumbinary disc simulation has lost 30 per cent of its initial mass.
| Name | i0 (°) | j | irs(°) | Mp/M | Planet orbit | Md/M | Disc orbit | Break radius/ab | Time (Tb) |
|---|---|---|---|---|---|---|---|---|---|
| Low-j-120 | 120 | 0.5 | – | 0.006 | L | 0.01 | L | – | 2110 |
| Low-j-130 | 130 | 0.5 | – | 0.006 | L | 0.01 | L | – | 1540 |
| Low-j-140 | 140 | 0.5 | – | 0.006 | L | 0.01 | L | – | 1070 |
| Low-j-150 | 150 | 0.5 | – | 0.006 | L | 0.01 | L | – | 850 |
| Low-j-160 | 160 | 0.5 | – | 0.006 | L | 0.01 | L | 6 | 1560 |
| Low-j-170 | 170 | 0.5 | – | 0.006 | L | 0.01 | L | – | 2520 |
| Mid-j-100 | 100 | 1.5 | 139 | 0.018 | L | 0.03 | L | – | 2880 |
| Mid-j-110 | 110 | 1.5 | 139 | 0.018 | CRESCENT | 0.03 | L | – | 2320 |
| Mid-j-120 | 120 | 1.5 | 139 | 0.018 | CRESCENT | 0.03 | CRESCENT → L | – | 1830 |
| Mid-j-130 | 130 | 1.5 | 139 | 0.018 | CRESCENT | 0.03 | CRESCENT | – | 1410 |
| Mid-j-140 | 140 | 1.5 | 139 | 0.018 | irs | 0.03 | |$i_{\rm rs}\, \rightarrow$| L | – | 1140 |
| Mid-j-150 | 150 | 1.5 | 139 | 0.018 | CRESCENT | 0.03 | CRESCENT (?) | – | 870 |
| Mid-j-160 | 160 | 1.5 | 139 | 0.018 | CR | 0.03 | CR | 2.5 | 1000 |
| Mid-j-170 | 170 | 1.5 | 139 | 0.018 | CR | 0.03 | CR | – | 2620 |
| High-j-100 | 100 | 2.5 | 129 | 0.031 | L | 0.05 | L | – | 2380 |
| High-j-110 | 110 | 2.5 | 129 | 0.031 | CRESCENT | 0.05 | CRESCENT → L | – | 2310 |
| High-j-120 | 120 | 2.5 | 129 | 0.031 | CRESCENT | 0.05 | CRESCENT | – | 2230 |
| High-j-130 | 130 | 2.5 | 129 | 0.031 | irs | 0.05 | |$i_{\rm rs}\, \rightarrow$| L | – | 4200 |
| High-j-140 | 140 | 2.5 | 129 | 0.031 | CRESCENT | 0.05 | CRESCENT → L | – | 1980 |
| High-j-150 | 150 | 2.5 | 129 | 0.031 | CR | 0.05 | CR | – | 2020 |
| High-j-160 | 160 | 2.5 | 129 | 0.031 | CR | 0.05 | CR | 2.5 | 1110 |
| High-j-170 | 170 | 2.5 | 129 | 0.031 | CR | 0.05 | CR | – | 2620 |
| Name | i0 (°) | j | irs(°) | Mp/M | Planet orbit | Md/M | Disc orbit | Break radius/ab | Time (Tb) |
|---|---|---|---|---|---|---|---|---|---|
| Low-j-120 | 120 | 0.5 | – | 0.006 | L | 0.01 | L | – | 2110 |
| Low-j-130 | 130 | 0.5 | – | 0.006 | L | 0.01 | L | – | 1540 |
| Low-j-140 | 140 | 0.5 | – | 0.006 | L | 0.01 | L | – | 1070 |
| Low-j-150 | 150 | 0.5 | – | 0.006 | L | 0.01 | L | – | 850 |
| Low-j-160 | 160 | 0.5 | – | 0.006 | L | 0.01 | L | 6 | 1560 |
| Low-j-170 | 170 | 0.5 | – | 0.006 | L | 0.01 | L | – | 2520 |
| Mid-j-100 | 100 | 1.5 | 139 | 0.018 | L | 0.03 | L | – | 2880 |
| Mid-j-110 | 110 | 1.5 | 139 | 0.018 | CRESCENT | 0.03 | L | – | 2320 |
| Mid-j-120 | 120 | 1.5 | 139 | 0.018 | CRESCENT | 0.03 | CRESCENT → L | – | 1830 |
| Mid-j-130 | 130 | 1.5 | 139 | 0.018 | CRESCENT | 0.03 | CRESCENT | – | 1410 |
| Mid-j-140 | 140 | 1.5 | 139 | 0.018 | irs | 0.03 | |$i_{\rm rs}\, \rightarrow$| L | – | 1140 |
| Mid-j-150 | 150 | 1.5 | 139 | 0.018 | CRESCENT | 0.03 | CRESCENT (?) | – | 870 |
| Mid-j-160 | 160 | 1.5 | 139 | 0.018 | CR | 0.03 | CR | 2.5 | 1000 |
| Mid-j-170 | 170 | 1.5 | 139 | 0.018 | CR | 0.03 | CR | – | 2620 |
| High-j-100 | 100 | 2.5 | 129 | 0.031 | L | 0.05 | L | – | 2380 |
| High-j-110 | 110 | 2.5 | 129 | 0.031 | CRESCENT | 0.05 | CRESCENT → L | – | 2310 |
| High-j-120 | 120 | 2.5 | 129 | 0.031 | CRESCENT | 0.05 | CRESCENT | – | 2230 |
| High-j-130 | 130 | 2.5 | 129 | 0.031 | irs | 0.05 | |$i_{\rm rs}\, \rightarrow$| L | – | 4200 |
| High-j-140 | 140 | 2.5 | 129 | 0.031 | CRESCENT | 0.05 | CRESCENT → L | – | 1980 |
| High-j-150 | 150 | 2.5 | 129 | 0.031 | CR | 0.05 | CR | – | 2020 |
| High-j-160 | 160 | 2.5 | 129 | 0.031 | CR | 0.05 | CR | 2.5 | 1110 |
| High-j-170 | 170 | 2.5 | 129 | 0.031 | CR | 0.05 | CR | – | 2620 |
2.2 Low-j three-body system
We first consider three-body simulations with a relatively low third body angular momentum of j = 0.5 < jcr = 0.91. The prograde-polar stationary inclination calculated with equation (4) is is = 82°. There is no retrograde-polar stationary inclination in this case since the angular momentum ratio is lower than the critical defined by equation (5). The two left-hand panels in Fig. 1 show the results of the Low-j three-body simulations. The top left plot shows the inclination phase plot, i cos ϕ−i sin ϕ. In this low-mass case, all of the orbits shown in the phase plot are librating, meaning the particle angular momentum vector is precessing about the generalized prograde-polar inclination. All of the orbits begin at their respective initial inclinations and carry out their paths in a counterclockwise direction. The coloured triangle is located at the start of the path and points to the initial direction. The lower left panel shows the eccentricity phase plot eb cos ϕ − eb sin ϕ. All of the results in the eccentricity phase plot begin at eb = 0.8, and sin ϕ = 1 and carry out their paths in a clockwise direction. In all cases, the binary eccentricity oscillates and initially increases from its initial value.
The upper panels show the i cos ϕ−i sin ϕ phase plots, while the lower panels show the eb cos ϕ−eb sin ϕ phase plane. The left-hand panels are for the three-body systems (see Section 2), while the right-hand panels are for the circumbinary disc systems (see Section 3). The dashed lines represent any part of the system that is simulated past the point where the disc has lost |$30 {{\ \rm per\ cent}}$| of its initial mass (see Fig. 6). All simulations, both disc and planet, have j = 0.5 initially. The three-body system has a planet with mass |$M_{\rm p} = 0.006M$| and initial planet distance from system centre of mass of |$r=20 \, a_{\rm b}$|. The initial mass of the disc is |$M_{\rm d} = 0.01M$|. Additionally, the pink dot on each of the right plots indicates the point where that circumbinary disc broke apart at time |$T=2600 \, T_{\rm b}$|.
2.3 Mid-j three-body system
We now consider a more massive planet such that the angular momentum ratio j = 1.5 > jcr = 0.91. The prograde-polar stationary inclination calculated with equation (4) is is = 73° and the retrograde-polar stationary inclination is irs = 139°. The left two panels of Fig. 2 show the results of the Mid-j three-body simulations. The top phase plot shows that for initial inclination i0 ≲ 100°, the orbit is librating about the prograde-polar stationary inclination. For i0 ≳ 160°, the orbits are in retrograde-polar circulation, meaning that the planet angular momentum vector is precessing about the negative of the binary angular momentum vector. In the approximate initial inclination range, 110–150° we see crescent orbits. These crescent orbits are not nested on a common centre like the circulating or librating orbits. The Mid-j-140 three-body simulation does not show significant variation from its initial inclination in the phase plot because of its proximity to the retrograde-polar stationary inclination at irs = 139°. Additionally, as seen by the coloured triangles, the three-body systems with i0 < irs precess counterclockwise (prograde), while those with i0 > irs precess clockwise (retrograde).
Same as Fig. 1 except j = 1.5. The three-body systems have a planet with mass |$M_{\rm p} = 0.018M$| and initial planet distance from system centre of mass |$r=20 \, a_{\rm b}$|. The initial mass of the disc is |$M_{\rm d} = 0.03M$|. Additionally, the pink dot on each of the right plots indicates the point where that circumbinary disc broke apart at |$T=240 \, T_{\rm b}$|.
Looking at the lower left panel in Fig. 2, we see that the binary eccentricity oscillates in all cases, but the initial direction varies depending upon the initial inclination. If i0 < irs, then the eccentricity initially increases, while for i0 > irs, the eccentricity initially decreases. The initial horizontal direction of the paths in the inclination phase plot correlates with the initial vertical direction of the eccentricity plots. The orbital paths that have rightward-facing triangles have an initially decreasing eb, while orbits with leftward-facing triangles have an initially increasing eb.
To better understand what the phase plots represent in three dimensions, Fig. 3 connects the angular momentum evolution to the phase plots in Fig. 2. From this depiction, we are able to see the true shapes of these orbital paths without the alterations that are made when mapping the spherical path on to a Cartesian coordinate system.
Precession paths in 3D for the Mid-j simulations. Each row corresponds to a different initial inclination, i0, that is highlighted in the inclination phase plot on the right that is the same as the upper left panel in Fig. 2. The left and middle columns show 3D visual representations of the three-body simulations. The left column is viewed from azimuth |$\, =0^{\circ }$| and elevation |$\, =-5^{\circ }$|, while the middle column is viewed from azimuth |$\, =50^{\circ }$| and elevation |$\, =-20^{\circ }$|. The two yellow stars in the centre are the binary stars that are shown at apastron, the dotted coloured line and circle outside of the grey sphere are the third body with its trajectory, the black arrow is the angular momentum vector of the third body, and the solid coloured line on the sphere is the path of the planet’s angular momentum vector for a complete precession. The black and grey dotted lines on the sphere show the plane in which the binary orbits.
2.4 High-j three-body system
We now consider our most massive planet with the highest angular momentum ratio, j = 2.5 > jcr = 0.91. The prograde-polar stationary inclination calculated with equation (4) is is = 70°, and the retrograde-polar stationary inclination is irs = 129°. The left two panels of Fig. 4 show the results of these High-j three-body simulations. The phase plot behaviour is similar to the Mid-j simulations, but with some of the orbital phenomenon shifted downwards by about 10°. For i0 ≲ 100°, the orbit is librating about is. For i ≳ 150°, the orbits are in retrograde circulation. We find the crescent orbit regime to be approximately within 110–140°. Similar to the behaviour of Mid-j-140, we find that High-j-130 does not deviate significantly from its initial inclination due to its proximity to the retrograde-polar stationary inclination at irs = 139°. When looking at the eccentricity plots, we find the similar behaviours as described in Section 2.3.
Same as Fig. 1 except j = 2.5. The three-body systems have a planet mass |$M_{\rm p} = 0.031M$| and initial planet distance from system centre of mass |$r=20\, a_{\rm b}$|. The initial mass of the disc is |$M_{\rm d} = 0.05M$|.
3 CIRCUMBINARY DISC SIMULATIONS
In this section, we explore hydrodynamical simulations of a circumbinary gas disc with the same initial properties (i0 and j) as the three-body simulations in the previous section in order to understand the long-term evolution of a massive circumbinary disc.
3.1 Circumbinary disc simulation set-up
One detriment to note in our SPH simulations is that the flow in the central gap region is not well resolved by the code in the intrinsically 3D flows. That flow takes the form of rapid low-density gas streams (e.g. Artymowicz & Lubow 1996; Mösta, Taam & Duffell 2019; Muñoz, Miranda & Lai 2019). This causes some uncertainty in the binary evolution.
3.2 Low-j circumbinary disc
Here, we describe the Low-j hydrodynamical simulations of circumbinary discs, where j = 0.5 < jcr = 0.91. The right two plots of Fig. 1 show the orbital phase and eccentricity results of these systems. The top right plot shows the i cos ϕ−i sin ϕ plane and the lower right plot shows the eccentricity phase plot eb cos ϕ−eb sin ϕ. When compared to the equivalent three-body system (left), we see general agreement in the phase plots for i0 ≲ 160°. While the three-body orbits are exactly periodic, the disc simulations include dissipation that causes the disc to gradually align towards the prograde-polar stationary inclination in a spiral-like motion. The upper panel of Fig. 6 shows the mass of each disc in time. The lines are dashed after the disc mass has lost 30 per cent of the initial value. This corresponds to the dashed lines in the phase plots in Fig. 1.
Snapshots from the Low-j-160 broken-disc system. The top row shows the initial conditions of the circumbinary disc from three different viewing angles. The middle row shows the time where we classify the disc as broken and begin analysing the data from the inner and outer disc separately. Note that |$t=2600\, T_{\rm b}$| is where the red and blue lines meet in Fig. 7. The bottom row shows the simulation at a later time where the two discs are fully separated, and the inner disc has gone nearly polar.
The disc mass evolution of all simulations. The dashed lines represent where each simulation has lost more than |$30 {{\ \rm per\ cent}}$| of its initial mass.
The higher initial inclinations, where i ≳ 160°, show somewhat different behaviour. For Low-j-160, the disc initially follows a similar evolution to the particle, but then begins to move quickly towards prograde-polar alignment. This behaviour can be seen in the pink line in Fig. 1, which experienced a disc break at a time of |$t = 2600 \, T_{\rm b}$|. The disc breaks into two disjoint rings at a break radius of around |$R = 6 \, a_{\rm b}$|. This is visually displayed in Fig. 5. Disc breaking occurs when the radial communication in the disc occurs on a longer time-scale than the precession time-scale (Papaloizou & Terquem 1995; Larwood et al. 1996; Nixon, King & Price 2013). Protoplanetary discs are typically in the wave-like regime since H/R ≫ α (Papaloizou & Pringle 1983) and bending waves, that propagate at about half the sound speed, communicate the warp (Papaloizou & Terquem 1995). Also, Fig. 7 shows the inclination phase plot for just this simulation. The pink line here is the same as in the upper right panel of Fig. 1. The blue and red lines show the inner and outer disc, respectively, plotted. The inner disc is analysed at |$R=5 \, a_{\rm b}$|, while the outer disc is set at |$R=25 \, a_{\rm b}$|. The outer disc continues much on the original path, while the inner disc quickly spirals to prograde-polar. The combination of these two paths is why we observe a sharp, but looping path towards prograde-polar alignment. As seen in Fig. 5, the radius of the disc break grows over time. Thus, the inclination of the total disc is dominated by the inclination of the inner disc at larger times.
The i cos ϕ−i sin ϕ phase plot for the Low-j-160 broken disc. The blue line represents the inner disc at |$R=5\, a_{\rm b}$|, while the red line represents the outer disc, set at |$R=25\, a_{\rm b}$| after the disc has broken at time of |$t=2600\, T_{\rm b}$|. The pink line is the same as shown in Fig. 1.
Similar to the three-body paths, we find that the initial direction of the orbital paths (denoted with the coloured triangles) corresponds with the initial direction of the eccentricity plots. The orbital paths with rightward-facing triangles have an initially decreasing eb.
3.3 Mid-j circumbinary disc
The right two panels of Fig. 2 show the phase plots of the Mid-j circumbinary disc simulations with j = 1.5 > jcr = 0.91. When compared to the equivalent Mid-j three-body systems (left), we see general agreement for the highest inclination case of the Mid-j-170 system. For i0 ≲ 110°, we see librating precession and the disc spirals towards prograde-polar alignment. The Mid-j-110 system in the three-body case shows a crescent orbit, but this is not seen for the disc. However, in the approximate range of 120–150°, the disc begins at first on a crescent-type orbit. In this range, the disc moves towards libration and prograde-polar alignment. For Mid-j-120 and Mid-j-130, we see at least one complete phase of a crescent orbit before the disc begins librating. It is important to note that some of these transitions into librating orbits occur after the disc has lost more than |$30 {{\ \rm per\ cent}}$| of its initial mass, and thus direct comparison to the three-body system shown may not be applicable. The value for the angular momentum ratio j is decreasing during the simulation because mass is lost from the disc and accreted on to the binary (see the middle panel of Fig. 6). However, for the Mid-j-120 simulation, the mass remains high while the orbit transitions to libration.
For Mid-j-140 and Mid-j-150, we observe the beginning of a descent to prograde-polar alignment, which differs greatly from the respective three-body simulations. Note that Mid-j-140 begins very close to the retrograde-polar stationary inclination irs and we see that the initial direction of this orbital path is opposite to that of the respective three-body orbit. We find that most initial inclinations seem to ignore the retrograde-polar stationary inclination irs in order to evolve towards a librating orbit around prograde-polar alignment is.
Finally, for i = 160°, we again observe the disc breaking at a time of |$t=300 \, T_{\rm b}$|. We can see the moment of the disc break in the pink inclination and eccentricity paths where there is a small bump and pink dot to indicate the time of the disc break. The approximate radius of this disc break was at |$R=2.5\, a_{\rm b}$|.
For i ≲ 130°, we find that the initial direction of the orbital paths, denoted with the coloured triangles, corresponds with the initial direction of the eccentricity plots. However, for Mid-j-140 and Mid-j-150 we find the eccentricities initially increasing unlike their respective three-body paths. For Mid-j-160, we see an initial decrease followed by an increase in a counterclockwise rotation unlike the clockwise rotation in the three-body path. For Mid-j-170, we see general agreement with both the eccentricity and inclination phase plots.
3.4 High-j circumbinary system
Here, we describe simulations of circumbinary discs with j = 2.5 > jcr = 0.91. The right two plots of Fig. 4 show the results of the High-j circumbinary disc simulations. When compared to the equivalent three-body system, we see general agreement for only the i0 = 100° and i0 = 170° systems. For High-j-100, we see a librating orbit spiraling towards prograde-polar alignment. For High-j-110 and High-j-120, we see at least one phase of crescent orbits that descend into librating orbits around prograde-polar alignment. While High-j-120 has lost a significant amount of mass before it begins librating, High-j-110 has not (see the lower panel of Fig. 6). Specifically for High-j-120, we see at least two crescent orbits before descending to prograde-polar libration. For High-j-130, we observe the beginning of a descent to prograde-polar alignment, which differs greatly from the respective three-body simulation. Note that High-j-130 is the closest to the retrograde-polar stationary inclination irs for j = 2.5, and we see that the initial direction of this orbital path is opposite to the respective three-body simulation.
For i ≲ 120°, we find that the initial direction of the orbital paths (denoted with the coloured triangles) corresponds with the initial direction of the eccentricity plots (as discussed in Section 2.4). The High-j-140 and High-j-150 eccentricities initially decrease like their respective three-body paths. For High-j-160, we see an initial decrease followed by an increase in a counterclockwise rotation unlike the clockwise rotation in the three-body path. Similar behaviour can be seen in the pink line in Fig. 1. This system also experienced a disc break at a time of |$t=270\, T_{\rm b}$|. The disc breaks into two disjoint rings at a break radius of around |$R=2.5\, a_{\rm b}$|. For High-j-170, we see general agreement with both the eccentricity and inclination phase plots.
4 CONCLUSIONS
We have investigated the evolution of a massive circumbinary disc in order to examine the effect of the disc mass on the probability of polar alignment. We consider a disc that begins with an inclination that is closer to retrograde alignment than prograde alignment to the binary orbit. We have found that at least initially a circumbinary disc displays similar behaviour to a three-body simulation with the same angular momentum ratio of the outer body to the inner binary. Dissipation and radial communication across the disc lead to differences later on. For a narrow range of initial disc inclinations around i = 160°, the disc communication time-scale is longer than the nodal precession time-scale and this leads to breaking, where the disc has two disjoint rings that can precess independently. Qualitatively, the behaviour is not affected by the disc break.
With three-body simulations, we have explored the dynamics of a massive third body orbiting an eccentric binary. There are three different types of nodal precession for massive bodies. First, the third body may be circulating, meaning that it precesses about the binary angular momentum vector (prograde circulation) or the negative of the binary angular momentum vector (retrograde circulation). Secondly, if the initial inclination is close to the generalized prograde-polar stationary inclination, it may be librating, meaning that it precesses about the generalized prograde-polar state. For inclinations between the retrograde circulation and the libration, the particle may be in a crescent-type orbit in the inclination phase plot. These orbits do not have a common centre. There are two stationary inclinations for which there is no nodal precession of the orbit. The prograde-polar stationary inclination lies at the centre of the librating orbits. The retrograde-polar stationary inclination exists for sufficiently high angular momentum of the third body and lies within the crescent orbits.
We have found that the dissipation within the circumbinary disc leads to (prograde-)polar alignment for all orbits that begin in the crescent orbit regime. Even if a disc begins very close to the retrograde-polar stationary inclination, it does not stay there. Therefore, discs are not expected to be aligned to the retrograde-polar stationary inclination. Similarly, we do not expect planets to be aligned to the retrograde-polar stationary inclination after the disc has dissipated. Previous models for the probability of polar alignment are based on a massless disc (Aly et al. 2015; Zanazzi & Lai 2018). However, we have shown that the phase space for which a retrograde disc moves towards polar alignment is larger than would be predicted by just the librating orbits. A massive and close to retrograde circumbinary disc eventually moves towards the prograde-polar stationary inclination. This result has implications for the formation of planets around an eccentric orbit binary since polar discs may form from a wide range of initial disc misalignments.
The existence of crescent-type precession orbits and the retrograde-polar stationary inclination relies on a large disc angular momentum (see the critical angular momentum ratio given in equation 5). This corresponds to a disc mass of a few per cent or more of the binary mass for a disc that extends to a radius of about 10 times the binary separation that we have considered here. We now have a number of observational tracers of disc mass that show that disc masses are typically up to a few per cent of a solar mass (Anderson et al. 2022a, b); some may be up to |$0.2\, \rm M_\odot$| (McClure et al. 2016), even larger than those considered in this work. For these disc masses, the range of initial inclinations that evolve to a polar orientation may be large. For example, in Fig. 2, discs with initial inclination in the range 60–150° evolve towards polar. Turbulence in giant molecular clouds is the expected way that these discs can form initially misaligned (Bate 2018). If the initial orientation of discs relative to the binary is isotropic, then polar discs may be more common than coplanar discs around eccentric binaries even when then disc is massive. Furthermore, as the disc dissipates and its mass decreases, polar alignment becomes even more likely (see Fig. 1 where inclinations ≳20° evolve towards polar).
ACKNOWLEDGEMENTS
We thank an anonymous referee for useful comments that improved the manuscript. Computer support was provided by UNLV’s National Supercomputing Center. We acknowledge support from NASA through grants 80NSSC19K0443 and 80NSSC21K0395. SHL acknowledges visitor support from the Simons Foundation. This research was supported in part by the National Science Foundation under grant no. NSF PHY-1748958. Simulations in this paper made use of the rebound code, which can be downloaded freely from http://github.com/hannorein/rebound. Visualizations of phantom simulations were created through Splash from Price (2007).
DATA AVAILABILITY
The SPH simulation results in this paper can be reproduced using the phantom code (Astrophysics Source Code Library identifier ascl.net/1709.002). The N-body simulation results can be reproduced with the rebound code (Astrophysics Source Code Library identifier ascl.net/1110.016). The data underlying this article will be shared on reasonable request to the corresponding author.
