Multiplanet disc interactions in binary systems

ABSTRACT

We investigate the evolution of a multiplanet–disc system orbiting one component of a binary star system. The planet–disc system is initially coplanar but misaligned to the binary orbital plane. The planets are assumed to be giants that open gaps in the disc. We first study the role of the disc in shaping the mutual evolution of the two planets using a secular model for low initial tilt. In general, we find that the planets and the disc do not remain coplanar, in agreement with previous work on the single planet case. Instead, the planets and the disc undergo tilt oscillations. A high-mass disc between the two planets causes the planets and the disc to nodally precess at the same average rate but they are generally misaligned. The amplitude of the tilt oscillations between the planets is larger while the disc is present. We then consider higher initial tilts using hydrodynamical simulations and explore the possibility of the formation of eccentric Kozai–Lidov (KL) planets. We find that the inner planet’s orbit undergoes eccentricity growth for a large range of disc masses and initial misalignments. For a low disc mass and large initial misalignment, both planets and the disc can undergo KL oscillations. Furthermore, we find that sufficiently massive discs can cause the inner planet to increase its inclination beyond 90° and therefore to orbit the binary in a retrograde fashion. The results have important implications for the explanation of very eccentric planets and retrograde planets observed in multiplanet systems.

accretion, accretion discs, hydrodynamics, planets and satellites: formation, binaries: general

1 INTRODUCTION

The discovery of thousands of exoplanets presents a variety of exotic configurations of their orbits around the host star. Horch et al. (2014), Deacon et al. (2016), Matson et al. (2018), and Ziegler et al. (2018) estimated that roughly 50 per cent of these planets might be hosted by binary systems. Planet formation in a binary system is particularly worth investigating in order to explain high inclination, eccentric close-in planets that have been observed.

A large fraction of observed exoplanet orbits is highly misaligned with respect to the spin of the host star (Triaud et al. 2010; Winn et al. 2010; Winn & Fabrycky 2015; Triaud 2018). The misalignment is measured with the Rossiter–McLaughlin effect (McLaughlin 1924; Rossiter 1924). There are three main mechanisms that have been proposed to explain the observed misalignments. First, a planet might lie on a misaligned orbit as a result of the interaction with another object, either a star or a planet (e.g Takeda, Kita & Rasio 2008). Secondly, the misalignment might arise because the protoplanetary disc feels a torque exerted by a binary companion (e.g. Batygin 2012). Finally, the magnetic field of the star may induce changes in the orientation of the spin of the star relative to the disc angular momentum (Lai, Foucart & Lin 2011). In this work, we consider the second effect.

Among the planets that show projected spin–orbit misalignments, a significant fraction shows apparent retrograde motion (Triaud et al. 2010). For instance, observations of the Rossiter–McLaughlin effect confirmed the prediction made by Anderson et al. (2010) on the nature of the retrograde orbit of WASP-17b.

According to the current commonly accepted scenario, planets form in a protoplanetary disc. Giant planets are believed to be embedded in the disc around the host star in the earliest stages of their formation (Pollack et al. 1996; Hubickyj, Bodenheimer & Lissauer 2005; Papaloizou & Nelson 2005). When the planet grows to the mass of Neptune, it undergoes runaway-gas accretion and its tidal torque opens a gap in the gas disc (Lin & Papaloizou 1986; D’Angelo, Henning & Kley 2002; Bate et al. 2003). In this scenario, these planets form outside the snow line because the presence of ice enhances the surface density of the planetesimals (Lecar et al. 2006). They are also required to form within the lifetime of the gas disc, roughly a few million years (Haisch, Lada & Lada 2001; Pascucci et al. 2006; Fedele et al. 2010). Terrestrial planets are believed to form from agglomeration of Moon to Mars sized planetary embryos and planetesimals (Chambers 2001; O’Brien, Morbidelli & Levison 2006). Therefore, while these kind of small planets are most likely to remain coupled with the gaseous disc, giant planets interact more strongly with the disc. The study of these kind of interactions is essential in order to explain the presence of giant planets on inclined and/or eccentric orbits close to the host star.

Understanding this kind of system is key for explaining observations such as the Kepler-56 planetary system (Huber et al. 2013). This is a red giant star hosting two close-in planets of Neptune and Saturn size, respectively, on highly inclined orbits with respect to the spin of the star. Misalignments between the planet orbits and the host star spin are believed to be caused by the torque exerted from a wide orbiting companion star. Radial velocity measurements revealed indeed a third object in Kepler-56.

Once a planet is massive enough to form a gap in the gas disc, the binary torque dominates that of the disc. Picogna & Marzari (2015), Lubow & Martin (2016), and Martin et al. (2016) found that it is more likely for the planet and the disc to undergo inclination oscillations on long time-scales and not to remain coplanar. The phase angles of the components that orbit the primary star can either be in a librating or circulating state. In the librating state, the phase difference between different components is limited. Therefore, the mean precession rates for the librating objects are locked and they evolve together. If instead the phase angles have different mean precession rates, the system is circulating and the components evolve independently from one another.

In this paper, we extend the work done by Martin et al. (2016) and Lubow & Martin (2016) on a single planet–disc system to investigate multiplanet systems. We study the orbital evolution of two giant planets with a gas disc. In Section 2, we first use analytical calculations to investigate the role that the accretion disc has in shaping the planets orbital evolution for low initial misaligment. Then, in Section 3 we use smoothed particle hydrodynamical (SPH) simulations to explore higher initial misalignment and the possibility of planets undergoing Kozai–Lidov oscillations (Kozai 1962; Lidov 1962). In principle, if the two planets are circulating, the inner planet is likely to undergo KL oscillations if its inclination exceeds the critical value i_crit = 39°. On the contrary, if the distance between the planets is small enough, there is a chance that the outer planet suppresses the KL mechanism acting on the inner planet. From the results presented in Lubow & Martin (2016), it is clear that the disc mass is a crucial parameter to determine whether the planets evolve together or separately. We therefore investigate in this paper the role of the disc mass in shaping the planets evolution and determining the possibility of KL oscillations of the inner planet. We draw our conclusions in Section 4.

2 ANALYTICAL CALCULATIONS: SECULAR EVOLUTION

The system under investigation is composed of a primary star surrounded by a gas disc with two planets embedded in it and a companion star orbiting the primary. In the secular analytical calculations, we assume that both planets already carved a gap in the disc so that the objects orbiting the primary are in the following order: inner planet, inner disc, outer planet, outer disc, and companion star.

The six bodies interact only via gravitational forces, the discs are taken to be non-viscous, rigid, and flat. Since there is no dissipation within the disc, the amplitude of the inclination oscillations between the discs and the planets remains constant with time in the case where there is only one planet interacting with one disc. In the multiplanet–disc system, the oscillations become more chaotic. The orbital planes of the planets and the two discs are initially coplanar but misaligned with respect to the binary orbital plane.

The binary is assumed to be circular, equal mass with total mass M = M₁ + M₂ and separation a. We put two planets with mass |$M_{\rm p,i}=0.001\, M$| around the primary star located at orbital radii |$d_1=0.05\, a$| and |$d_2=0.15\, a$|⁠, respectively.

We adopt the same procedure outlined in Lubow & Martin (2016) describing the time-dependent tilt of each object relative to the binary plane with the complex variable W(t) = l_x(t) + il_y(t), where the unit angular momentum vector as a function of time is |$\boldsymbol {l}(t)=(l_{\rm x}(t),l_{\rm y}(t),l_{\rm z}(t))$|⁠. Assuming that the tilts are small, we solve the secular evolution equations for the gravitational interaction between the different objects in the system (Lubow & Ogilvie 2001). According to this model, the torques are evaluated in the small angle approximation. As Lubow & Martin (2016) pointed out, the analytical results are therefore accurate for initial tilts up to about 20°. Furthermore, this approach assumes the orbits to remain circular therefore it cannot be used to study systems undergoing KL oscillations. We investigate higher inclination systems using SPH simulations in Section 3.

The interaction between the components of the system is described by a linear model with coupling coefficients. The evolution equations for the planets tilt W_p1, W_p2, and the two discs W_d1, W_d2 are

$$\begin{eqnarray} J_{\rm p1} \frac{\mathrm{ d}W_{\rm p1}}{\mathrm{ d}t} &=& {} iC_{\rm p1d1}(W_{\rm d1} - W_{\rm p1}) + iC_{\rm p1p2}(W_{\rm p2} - W_{\rm p1}) \nonumber \\ && +\, iC_{\rm p1d2}(W_{\rm d2} - W_{\rm p1}) - iC_{\rm p1s}W_{\rm p1} , \end{eqnarray}$$

(1)

$$\begin{eqnarray} J_{\rm p2} \frac{\mathrm{ d}W_{\rm p2}}{\mathrm{ d}t} &=& {} iC_{\rm p2d1}(W_{\rm d1} - W_{\rm p2}) + iC_{\rm p1p2}(W_{\rm p1} - W_{\rm p2}) \nonumber \\ && +\, iC_{\rm p2d2}(W_{\rm d2} - W_{\rm p2}) - iC_{\rm p2s}W_{\rm p2} , \end{eqnarray}$$

(2)

$$\begin{eqnarray} J_{\rm d1} \frac{\mathrm{ d}W_{\rm d1}}{\mathrm{ d}t} &=& {} iC_{\rm p1d1}(W_{\rm p1} - W_{\rm d1}) + iC_{\rm d1p2}(W_{\rm p2} - W_{\rm d1}) \nonumber \\ && +\, iC_{\rm d1d2}(W_{\rm d2} - W_{\rm d1}) - iC_{\rm d1s}W_{\rm d1} , \end{eqnarray}$$

(3)

$$\begin{eqnarray} J_{\rm d2} \frac{\mathrm{ d}W_{\rm d2}}{\mathrm{ d}t} &=& {} iC_{\rm d2d1}(W_{\rm d1} - W_{\rm d2}) + iC_{\rm d2p2}(W_{\rm p2} - W_{\rm d2}) \nonumber \\ && +\, iC_{\rm d1d2}(W_{\rm d1} - W_{\rm d2}) - iC_{\rm d2s}W_{\rm d2} , \end{eqnarray}$$

(4)

where p1, d1, p2, d2, and s denote, respectively, the inner planet, the inner disc, the outer planet, the outer disc, and the companion star with J being the respective angular momentum. The coupling coefficients describing the interaction of the various components of the system have the same form as the ones in Lubow & Martin (2016).

In order to understand the role of the disc in this type of architecture, we first consider only the evolution of the two planets misaligned with respect to the binary orbital plane, i.e. we ignore the presence of the disc. The coupling coefficient between the planets is given by

$$\begin{eqnarray} C_{\rm p1p2} = G\, M_{\rm p1}M_{\rm p2}\, K(d_1,d_2) , \end{eqnarray}$$

(5)

and the interaction between the planets and the companion star has the same form as equation (5) (see equation 7 in Lubow & Martin 2016). The kernel K(d₁, d₂) has the same form as in Lubow & Martin (2016). The left-hand panel of Fig. 1 shows the inclination and phase angles evolution of the two planets. The two objects oscillate away from each other and from the phase angle evolution we can see that each planet precesses at its own rate. Looking at the phase portrait in the right-hand panel of Fig. 1, we see that the planets are circulating, meaning that they are evolving independently from one another. The inclination of each component relative to the binary orbital plane is defined as i_j(t) = |W_j(t)| and the phase angle is ϕ_j(t) = tan ⁻¹(ℑ[W_j(t)]/ℜ[W_j(t)]), where |$j=\, \rm p1,\, p2,\, d1,\, d2,\, s$|⁠. The relative tilt between the i, j components is defined as Δi/i₀ = |W_i − W_j|, while the nodal phase difference is Δϕ = ϕ_i − ϕ_j.

Figure 1.

Left-hand panel: secular tilt (upper) relative to the binary orbital plane and phase (lower) evolution of a planet–planet system that orbits around one component of a circular binary. The time is in units of the orbital period of the binary P_b. Initially, the two planets are coplanar and lie on a misaligned orbit with inclination i₀ with respect to the binary orbital plane. The blue line represents the inner planet, while the green line is the outer planet. Right-hand panel: Phase portrait of the relative planet–planet tilt Δi/i₀ versus nodal phase difference Δϕ.

We then add the two portions of the accretion disc in the calculations and investigate their effect on the planets evolution. We assume that the mass contained in each portion of the disc is the same and the discs extend from |$R_{\rm in,1}=0.06\, a$| to |$R_{\rm out,1}=0.125\, a$| and from |$R_{\rm in,2}=0.175\, a$| to |$R_{\rm out,2}=0.25\, a$|⁠, respectively. The discs aspect ratio at the inner edge is H/R(R_{in, i}) = 0.035, the sound speed has the form c_s ∝ (R/R_{in, i})^−q with q = 0.75, the surface density profile in the two discs is Σ₁ ∝ (R/R_{in, 1})^−p and Σ₂ ∝ (R/R_{in, 2})^−p with p = 1.5. The index i = 1, 2 indicates the two discs. The constant of proportionality for the surface densities is determined by the mass contained in each disc, which is half of the total disc mass. The interaction between the planets and the discs is described as

$$\begin{eqnarray} C_{{\rm p}i{\rm d}j} = 2\pi \int _{R_{{\rm in},j}}^{R_{{\rm out},j}} G\, M_{{\rm p}i}\, R \Sigma _j(R)\, K(R,d_{\rm i})\, \mathrm{ d}R, \end{eqnarray}$$

(6)

where i = 1, 2 indicates the inner and outer planet and j = 1, 2 indicates the inner and outer disc, respectively. We also introduce a new coupling coefficient that describes the interaction between the inner and outer disc in the form

$$\begin{eqnarray} C_{\rm d1d2} = (2\pi)^2 \int _{R_{\rm in,2}}^{R_{\rm out,2}} \int _{R_{\rm in,1}}^{R_{\rm out,1}} r_1\, r_2\, \Sigma _1(r_1)\, \Sigma _2(r_2)\, K(r_1,r_2)\, \mathrm{ d}r_1\, \mathrm{ d}r_2, \nonumber \\ \end{eqnarray}$$

(7)

where x represent the radius inside the inner disc and y is the radius within the outer disc.

Fig. 2 shows the results for the same system but including the inner and outer discs. The total disc mass is |$M_{\rm d}=0.001\, M$|⁠. The upper left-hand panel of Fig. 2 shows the evolution of the inclination of the inner (blue) and outer (green) planet and of the inner (yellow) and outer (red) disc. The lower left-hand panel of Fig. 2 shows that the two planets (the blue and green line) precess independently. The right-hand panel of Fig. 2 shows the phase portrait of the relative planet–planet (blue), inner disc–outer planet (yellow), and outer planet–outer disc (green) tilt. We see that the inner planet is circulating with the outer planet but it is librating with the inner disc (see the blue and yellow lines in the lower left-hand panel). The outer planet is librating with the outer disc (the green line in the phase portrait). Therefore, a low-mass disc with a total mass of |$M_{\rm d}=0.001\, M$| leaves the planets free to circulate and precess each one with its own frequency around the primary star. However, the amplitude of the tilt oscillations relative to the binary orbital plane is significantly increased by the presence of the disc.

Figure 2.

Left-hand panel: secular tilt relative to the binary orbital plane (upper) and phase (lower) evolution of a planet–disc–planet–disc system that orbits around one component of a circular binary. The time is in units of the orbital period of the binary P_b. Initially, the two planets and the two discs are coplanar and lie on a misaligned orbit with respect to the binary orbital plane. The blue, yellow, green, and red lines represent the inner planet, the inner disc, the outer planet, and the outer disc, respectively. The total mass of the disc is |$M_{\rm d}=0.001\, M$|⁠. Right-hand panel: Phase portrait of the relative planet–planet (blue), inner disc–outer planet (yellow), and outer planet–outer disc (green) tilt Δi/i₀ versus nodal phase difference Δϕ.

Fig. 3 shows the results for the same system but with a higher total disc mass |$M_{\rm d}=0.002\, M$|⁠. The lower left-hand panel of Fig. 3 shows that the two planets (the blue and green line) precess at the same rate together with the portion of the disc between them (the yellow line). The outer disc (the red line) behaves independently with respect to the three inner components. However, the outer planet librates with the outer disc for very brief intervals during the evolution of the system, as seen in the closed green lines in the right-hand panel of Fig. 3. Also, in that panel, we see that the inner planet is librating with the outer planet (the blue line) and the outer planet librates with the inner disc (the yellow line). Therefore, increasing the disc mass results in librating solutions of the two planets and the portion of the disc between them. For a higher mass disc, the relative peak inclination drops (the blue line in right-hand panel of Fig. 3), as expected from the action of the secular resonance (see fig. 2 in Lubow & Martin 2016). The amplitude of the tilt oscillations of the planets relative to the binary orbital plane is again much larger as a result of the presence of the disc.

Figure 3.

Same as Fig. 2 but with total disc mass |$M_{\rm d}=0.002\, M$|⁠.

We see from Figs 1–3 that the presence of the two discs leads also to an increase in the mutual inclination between the planets. When the disc is present, the inclination of the inner planet relative to the binary orbital plane increases by roughly a factor 2 with respect to the initial inclination. Therefore, even if the same system starts with a i₀ = 20° inclination, it is likely that the planet reaches an inclination above the critical angle for the KL mechanism to operate. Since the secular model is valid for small tilts (i.e. up to 20°), we need to explore the higher inclination regime with hydrodynamical simulations.

3 MULTIPLANET DISC HYDRODYNAMICAL SIMULATIONS

In order to further study the behaviour of the planets embedded in the accretion disc, including the viscous evolution of the disc and higher initial inclinations, we perform smoothed particle hydrodynamics (SPH) numerical simulations using the code phantom (Lodato & Price 2010; Price & Federrath 2010; Nixon 2012; Price 2012; Nixon, King & Price 2013; Price et al. 2017). Planets embedded in accretion discs around one component of a binary star system have been previously studied with phantom (e.g. Lubow & Martin 2016; Martin et al. 2016). We perform SPH simulations using N = 5 × 10⁵ particles initially. The resolution of the simulation depends on N, the viscosity parameter α, and the disc scale height H. The Shakura & Sunyaev (1973) viscosity parameter is modelled by adapting the artificial viscosity according to the approach of Lodato & Price (2010). The main implication of this treatment is that in order to provide a constant α in the disc, the disc scale height H must be uniformly resolved. This is achieved by choosing both a surface density profile and a temperature profile that ensures this condition (Lodato & Pringle 2007), as discussed next.

Disc evolution strongly depends on the value of the protoplanetary discs viscosity that is not well known. In our simulations, we model disc viscosity using the artificial viscosity parameter α_AV (Lodato & Price 2010). The main consequence of modelling disc viscosity through the artificial viscosity parameter is that there is a lower limit below which a physical viscosity is not resolved in SPH (α_AV ≈ 0.1). Viscosities smaller than this value can produce disc spreading that is independent of the value of α_AV (Bate, Bonnell & Price 1995; Meru & Bate 2012).

3.1 Simulation set–up

If we simply place the formed planets into the disc (with no gap at their orbital location), the two planets and the disc can lose the initial coplanarity before the two gaps are carved. Furthermore, the planets can increase their mass significantly before carving the gaps and we would lose control over the mass of the planet when an equilibrium gap is established. Thus, we follow the methods described in Lubow & Martin (2016) and find the surface density of a disc with an equilibrium gap already established at the start of the simulation.

We consider an initial configuration with two planets embedded in an accretion disc that orbits the primary star. Both planets and the disc are coplanar and lie in the same plane as the binary. The equal mass binary has total mass M = M₁ + M₂ = 1 and a circular orbit in the x–y plane with separation a. We put two planets with mass |$M_{\rm p}=0.001\, M$| around the primary star located at |$d_1=0.05\, a$| and |$d_2=0.15\, a$|⁠, respectively. The two stars and the two planets are all modelled as sink particles. We choose the accretion radius for the particle removal to be |$R_{\rm acc}=0.0125\, a$| for the stars and 10 per cent of the Hill radius around each planet (⁠|$0.0044\, a$| and |$0.0013\, a$|⁠, respectively, for the inner and outer planet). The simulation results do not vary significantly for smaller values of these accretion radii. The mass and angular momentum of any particle that enters the sink radius are added to that of the sink particle.

We start with a very low initial disc mass of M_d,in = 4 × 10⁻⁶M and an initial disc radial extent from |$R_{\rm in}=0.0125\, a$| to |$R_{\rm out}=0.225\, a$|⁠. The disc then slightly expands outwards reaching |$R_{\rm out}=0.25\, a$| that corresponds to the tidal truncation radius of a disc in a coplanar binary system with separation a (Paczynski 1977).

The initial surface density profile of the disc is set as |$\Sigma \propto (R/R_{\rm in})^{-3/2}(1-\sqrt{R_{\rm in}/R})$| between |$R_{\rm in}=0.0125\, a$| and |$R_{\rm out}=0.225\, a$|⁠. With this density distribution and a sound speed given by c_s ∝ (R/R_{in, i})^−q with q = 0.75, the disc is uniformly resolved. The circumprimary disc aspect ratio at the inner edge is chosen to be H/R(R_in) = 0.035, which is a typical value for protoplanetary discs. The viscosity parameter is chosen to be α = 0.01.

The system evolves in the plane of the binary (i.e. i = 0°) until the planets carve two gaps in the disc. We find that a reasonably stable equilibrium gap structure is established in roughly 50 outer planet orbits, or 4 binary orbital periods. Fig. 4 shows the disc surface density profile after four binary orbits. The portion of the disc that lies inside the inner planet is accreted on to the star so the resulting system orbiting the primary star is composed of an inner planet, a disc, an outer planet, and an outer disc. The amount of mass contained in each disc is roughly half of the total disc mass.

Figure 4.

Surface density profile in log scale at a time of four binary orbits for a coplanar disc with initial mass of M_d,in = 4 × 10⁻⁶M. The surface density Σ is normalized by M/a², the binary mass divided by the square of its separation. The simulations start with the two planets embedded in the disc without a gap, as discussed in Section 3.1.

The disc density profile shape does not depend on the disc mass if the disc is locally isothermal and non-selfgravitating and if the planets are on fixed orbits. Note that planet migration occurs on a much longer time-scale than the time-scales considered in this work. Therefore, in the rest of our simulations, we rescale this density distribution obtained after four binary orbital periods to achieve the desired disc mass for each simulation. We then tilt the planets and the two portions of the disc all by the same angle so that they are mutually coplanar but misaligned to the binary orbital plane. The initial surface density distribution for the simulations is the equilibrium for a coplanar disc. The upper panels in Fig. 5 show the initial set-up for a simulation that is tilted by 40° to the binary orbital plane. The upper middle panel shows that the two planets and the disc are all initially aligned with each other, but tilted with respect to the binary orbital plane.

Figure 5.

Simulation with initial disc mass |$M_{\rm d}=0.002\, M$| and tilt 40°. The colour of the gas denotes the column density with yellow being about two orders of magnitude larger than blue. The red points show the primary star (located at the centre of each panel) and the two planets. The binary companion is not shown on this scale but is located along the positive x-axis. The left-hand panels show the view looking down on the binary orbital plane, the x−y plane, the middle panels show the x−z plane, and the right-hand panels show the y−z plane. The upper panels show the initial conditions and the lower panels show the system after a time of |$5\, P_{\rm orb}$|⁠.

3.2 Low inclination system

We first investigate the evolution of a system with initial total disc mass |$M_{\rm d}=0.001\, M$| and initial inclination i = 10°. The results are shown in Fig. 6. We evaluate the inclination i and phase angle ϕ for the inner and outer disc using a weighted density average from |$R_{\rm in,1}=0.06\, a$| to |$R_{\rm out,1}=0.125\, a$| and from |$R_{\rm in,2}=0.175\, a$| to |$R_{\rm out,2}=0.25\, a$|⁠, respectively.

Figure 6.

Left-hand panel: Evolution of the inclination (upper panel), phase angle (middle panel), and eccentricity (bottom panel) of the components that orbit the primary star: inner planet (the blue line), inner disc (yellow), outer planet (green), and outer disc (red) as a function of time in units of the binary orbital period P_b for a hydrodynamic simulation. Right-hand panel: Phase portrait of the relative planet–planet tilt Δi/i₀ versus nodal phase difference Δϕ. The system is initially coplanar, inclined by i₀ = 10° with respect to the binary plane and with total disc mass |$M_{\rm d}=0.001\, M$|⁠.

The two planets that orbit the primary are circulating with respect to one another, in agreement with the secular model with the same parameters. The inner disc (yellow) looks like it might be librating with the outer planet (green). However the inner disc is accreted after |$30\, P_{\rm b}$| therefore it is hard to make predictions on their mutual evolution. While in the analytical calculations we find the inner planet to be librating with the portion of the disc between the two planets, in our simulations this seems not to be the case and the two components are circulating. The simulations suggest that the analytic model may have underestimated the gap size around the planets. However, the behaviour of the planets with respect to each other (circulating) is the same in both the analytic model and the simulations. In the SPH simulations, the gap size evolves depending on the relative inclination of the planet and the disc. This effect is not taken into account in the analytic model. We show the results up to |$30\, P_{\rm b}$| because at that time the disc between the planets is essentially accreted. We discuss the rapid accretion and how it affects our simulations further in Section 4.

The right-hand panel of Fig. 6 shows the phase portrait of the relative planet–planet tilt. We can see that this is consistent with the right-hand panel of Fig. 2 that shows that the two planets are circulating with each other. The increase in relative tilt between the two planets is in good agreement with the secular theory.

The critical disc mass required for the planets to librate is higher in the SPH simulation compared to the linear theory. The reason is that in the simulations the disc is partially accreted on to the planets as the system evolves, while the analytical calculations do not include accretion of mass. Also, the surface density profile is assumed to be a power law for both discs in the calculations in Section 2, while in the numerical simulations the profile is determined by the planets carving gaps in the disc. Furthermore, the disc outer radius might be underestimated in the secular model because it does not take into account the viscous spreading of the disc.

3.3 KL oscillations of the inner planet

In order to explore the possibility of the inner planet undergoing eccentricity growth through KL oscillations, we increase the initial misalignment of the system with respect to the plane of the binary to i = 40°. If the initial inclination is higher, the portion of the disc between the two planets is likely to survive for longer since the misalignment between the planets and the disc becomes large enough to slow accretion of the disc material on to the planets.

The evolution of the inclination, phase angle, and eccentricity of the four components orbiting the primary star as a function of time are shown in Fig. 7 for initial total disc mass |$M_{\rm d}=0.001\, M$| and in Fig. 8 for disc mass |$M_{\rm d}=0.002\, M$|⁠. For both simulations, the four components around the primary all undergo precession on different time-scales initially. This can also been seen in the lower panels of Fig. 5 that show the higher mass disc at a time of |$t=5\, P_{\rm orb}$|⁠. The inner and outer discs are clearly misaligned to each other and to the binary orbit. Later on in the simulation, the inclination of the outer planet to the disc becomes sufficiently high that the torque from the outer planet is no longer strong enough to keep the gap open. Material flows through the gap and there is just one disc rather than two distinct components.

Figure 7.

Same as Fig. 6 but with initial inclination of i₀ = 40° and with total disc mass |$M_{\rm d}=0.001\, M$|⁠. Note that the yellow line is not distinguishable in the phase plot after roughly 50 orbits because the gap is filled and so the two parts of the disc behave in the same way.

Figure 8.

Same as Fig. 6 but with initial inclination of i₀ = 40° and with total disc mass |$M_{\rm d}=0.002\, M$|⁠. Note that the yellow line is not distinguishable in the phase plot after roughly 50 orbits because the gap is filled and so the two parts of the disc behave in the same way.

We can clearly see in Figs 7 and 8 that the inner planet, represented by the blue line, undergoes KL oscillations of eccentricity and inclination. For a lower mass disc, the inner planet reaches an inclination above the critical value on a slightly longer time-scale and this results in a delay in the first KL oscillation. Furthermore, we can see that a higher disc mass results in larger amplitude KL oscillations.

The KL oscillation period changes as the system evolves and we can see from Fig. 8 that the inner planet is starting a second KL cycle with a different periodicity at roughly |$80\, P_{\rm b}$|⁠. The comparison between the results obtained using two different disc masses shows that the period of the first KL oscillation is shorter if the mass of the disc is larger, in agreement with the analytical calculations that show shorter tilt oscillations periods for larger disc masses (see also right-hand panel of fig. 2 in Lubow & Martin 2016).

Starting the multiplanet–disc system at a higher inclination results in libration between the two discs after roughly |$60\, P_{\rm b}$| for the lower mass disc (Fig. 7) and about |$20\, P_{\rm b}$| for the higher mass disc (Fig. 8). The misalignment between the outer planet and the two discs becomes large enough so that the planet torque is weaker and the disc material can replenish the gap carved by the outer planet. In the bottom panel in Fig. 7), the disc between the planets starts to librate with the inner planet and becomes slightly eccentric at roughly |$60\, P_{\rm b}$|⁠. Then, the disc material replenishes the gap and the inner disc starts to librate, as it circularizes, with the outer disc. The same mechanism occurs on a shorter time-scale for the more massive case (Fig. 8).

We can see from Figs 7 and 8 that if the two planets are circulating the inner planet undergoes KL oscillations, while the outer planet remains on a circular orbit and does not affect the evolution of the inner one.

We further increase the mass of the disc to |$M_{\rm d}=0.02\, M$| in order to understand the role of the second planet on the KL oscillations of the inner planet. The evolution might be significantly different if the disc is massive enough for the two planets to start librating together. In this scenario, the outer planet could cause damping of the KL oscillations of the inner planet.

Fig. 9 shows the evolution of the inclination, phase angle, and eccentricity (left-hand panel) of each component that orbits the primary star and the phase portrait (right-hand panel) of the two planets. We can see from the right-hand panel that the two planets start to librate at the very beginning of the simulation when the mass of the disc is high enough. Furthermore, since the mutual inclination between the outer planet and the two discs increases very rapidly with time, the gap is replenished and therefore the disc starts to behave as a whole and undergoes two KL oscillations after roughly |$70\, P_{\rm b}$|⁠.

Figure 9.

Left-hand panel: same as Fig. 8 but with initial inclination of i₀ = 40° and with total disc mass |$M_{\rm d}=0.02\, M$|⁠. Right-hand panel: Phase portrait of the relative planet–planet tilt Δi/i₀ versus nodal phase difference.

The higher mass disc does not seem to affect the possibility of having KL eccentricity growth of the inner planet, infact, it strengthens its effect. The inner planet reaches a much higher eccentricity and inclination value compared to the lower disc mass cases. We can see from the left-hand panel of Fig. 9 that the inner planet inclination increases beyond 90°, meaning that the planet is orbiting in a retrograde fashion before coming back to the prograde orbit.

3.4 KL oscillations of the outer planet

We have also investigated the possibility for the outer planet to undergo KL oscillations. From Figs 7–9, we can see that the outer planet inclination decreases dramatically during the first tens of binary orbits. The higher the disc mass, the faster this decrease towards i = 0° occurs. Therefore, the system would have to start at higher inclinations than those explored so far in order for the outer planet inclination to exceed the critical angle for the KL mechanism to operate.

We ran a simulation with a disc mass of |$M_{\rm d}=0.001\, M$| starting at a 60° inclination. The results are shown in Fig. 10. The outer planet is indeed able to undergo eccentricity and inclination oscillations and to maintain a moderate eccentricity (e ≈ 0.5) for the whole duration of the simulation. At the same time, the two parts of the disc undergo KL oscillations as well since the initial inclination is well above the critical angle (Martin et al. 2014; Fu, Lubow & Martin 2015a; Lubow & Ogilvie 2017; Franchini, Martin & Lubow 2019). Since the discs become very eccentric, they accrete more efficiently on to the central star and are completely drained after |$40\, P_{\rm b}$|⁠.

Figure 10.

Same as Fig. 7 but starting at i₀ = 60°.

Fig. 11 shows the disc–binary–planet at the start of the simulation and at a time |$t=17\, P_{\rm orb}$|⁠. The outer planet is sufficiently misaligned to the disc that material flows through the planet gap. The disc itself is seen to be highly eccentric as a result of its KL oscillations.

Figure 11.

Same as Fig. 5 except for a simulation with initial disc mass |$M_{\rm d}=0.001\, M$| and tilt 60°. The upper panels show time t = 0, while the lower panels show time |$t=17\, P_{\rm orb}$|⁠.

4 CONCLUSIONS

In this paper, we investigated the interaction between two giant planets embedded in a protoplanetary disc that orbits one component of a binary system. The planets and the disc are initially coplanar but misaligned with respect to the binary orbital plane. We first investigated the effect of the disc in determining whether the planets librate together or precess independently (i.e. circulate) using the secular model (Lubow & Martin 2016). The results show that the presence of the disc leads to larger mutual inclinations between the two planets. Furthermore, for less massive discs (⁠|$M_{\rm d}\le 0.001\, M$|⁠) the planets circulate while for higher disc masses the planets precess together around the primary star.

Since the secular model assumes small tilts and its results are therefore accurate up to about 20°, we then used SPH simulations to investigate the possibility of KL oscillations of the planets starting with the planets and the disc misaligned by i = 40°. In particular, we find that the inner and outer planets, as well as the disc, can undergo eccentricity growth depending on the initial inclination of the planets–disc system and the disc mass. Furthermore, more massive discs lead to larger oscillation amplitudes of the inner planet eccentricity and inclination.

We confirmed the prediction of the secular model with our numerical simulations in terms of how the disc influences the planets behaviour. However, the critical disc mass for the planets to librate together is very different in the two approaches because of the presence of viscosity and therefore accretion of mass on to the planets in the SPH simulations. If the disc is as massive as |$M_{\rm d}=0.02\, M$|⁠, we find that the two planets do start to librate at the very beginning of the simulations but they switch to circulation when the disc mass decreases from the initial value. The inner planet reaches an eccentricity close to unity and transitions to a retrograde orbit before its oscillations are damped and it returns to orbit the primary star prograde. Therefore, the interaction between two planets and the massive disc where they are embedded might explain the formation of observed retrograde planets such as WASP-17b (Anderson et al. 2010).

One limitation of our simulations is that the disc disperses on a short time-scale. This is a consequence of several factors. First, the accretion radii of the planets and the central star cause material to accrete faster than it otherwise should. Secondly, there may be some overlap in the planet gaps that would cause the disc between the two planets to be depleted on a short time-scale (e.g. Morbidelli & Crida 2007). This may not be an important effect in the case of the highly misaligned systems that we have considered here. A more massive disc would result in larger tilt oscillations and a wider range of conditions for mutual libration between the planets than we found here. For the types of systems we considered, a higher disc mass would also mean that the inner planet would be more likely to undergo KL oscillations, while the outer planet would be less likely.

Note that if we want to further increase the mass of the accretion disc we would need to consider the effect of the disc self-gravity. Fu, Lubow & Martin (2015b) showed for an equal mass binary that the disc mass required to quench KL disc oscillations is a few per cent of the mass of the central star. Therefore, we expect that if we increase the mass of the disc up to, say, |$M_{\rm d}=0.05\, M$| the KL oscillations of the disc are likely to be suppressed.

ACKNOWLEDGEMENTS

We thank Giovanni Dipierro for useful comments and discussion. We thank Daniel Price for providing the phantom code for SPH simulations and acknowledge the use of splash (Price 2007) for the rendering of the figures. We acknowledge support from NASA through grant NNX17AB96G. Computer support was provided by UNLV’s National Supercomputing Center.