Orbital alignment of circumbinary planets that form in misaligned circumbinary discs: the case of Kepler-413b
Arnaud Pierens Richard P. Nelson


A Framework for Evaluating 
 6-DOF Object Trackers
    Mathieu Garon 0000-0003-1811-4156, Denis Laurendeau 0000-0003-2858-5955 and Jean-François Lalonde 0000-0002-6583-2364
    December 30, 2023
=========================================================================================================================


 Although most of the circumbinary planets detected by the Kepler spacecraft are on orbits that are closely aligned with the binary orbital plane, the systems Kepler-413 and Kepler-453 exhibit small misalignments of ∼ 2.5^∘. One possibility is that these planets formed in a circumbinary disc whose midplane was inclined relative to the binary orbital plane. Such a configuration is expected to lead to a warped and twisted disc, and our aim is to examine the inclination evolution of planets embedded in these discs. We employed 3D hydrodynamical simulations that examine the disc response to the presence of a modestly inclined binary with parameters that match the Kepler-413 system, as a function of disc parameters and binary inclinations. The discs all develop slowly varying warps, and generally display very small amounts of twist. Very slow solid body precession occurs because a large outer disc radius is adopted. Simulations of planets embedded in these discs resulted in the planet aligning with the binary orbit plane for disc masses close to the minimum mass solar nebular, such that nodal precession of the planet was controlled by the binary. For higher disc masses, the planet maintains near coplanarity with the local disc midplane. Our results suggest that circumbinary planets born in tilted circumbinary discs should align with the binary orbit plane as the disc ages and loses mass, even if the circumbinary disc remains misaligned from the binary orbit. This result has important implications for understanding the origins of the known circumbinary planets.





accretion, accretion discs –
                planet-disc interactions–
                planets and satellites: formation –
                hydrodynamics –
                methods: numerical




§ INTRODUCTION

To date, 11 circumbinary planets orbiting both components of a close binary system have been discovered by the Kepler spacecraft.  With the exception of Kepler-413b (Kostov et al. 2015) and Kepler-453b (Welsh et al. 2015) whose orbital planes are misaligned by ∼ 2.5^∘ with respect to the binary orbital planes, most of these planets have orbits that are closely aligned to those of the central binaries, supporting the idea that these planets formed in a circumbinary disc that was close to being coplanar with the central binary. 

A number of the Kepler circumbinary planets are on orbits that are close to the region of dynamical instability near the central binary (Holman & Wiegert 1999). Two scenarios have been put forward to explain the presence of planets in these orbits. In-situ formation, involving the accretion of planetesimals at the observed locations of the planets, has been suggested. However, various studies have highlighted the difficulty of building circumbinary planets in this way, due to the planetesimals being excited onto highly eccentric orbits (Paardekooper et al. 2012; Meschiari 2012 a,b; Lines et al. 2014; Bromley & Kenyon 2015), size-dependent pericentre alignment of planetesimals (Scholl et al. 2007)  or gravitational interactions with non-axisymmetric eccentric features in the disc (Marzari et al. 2008; Kley & Nelson 2010) all contributing to mutual collisions being disruptive rather than accretionary. Studies involving pebble accretion have yet to be undertaken, but one might expect that the large velocity perturbations induced in the gas and solid component of the disc in the inner regions would render planet formation via pebble accretion an inefficient process there.

 An alternative, more plausible,  possibility is that these planets formed further out in the disc, in a more accretion-friendly environment, and migrated to their present location through Type I (Ward et al. 1997; Tanaka et al. 2002) and/or Type II migration (e.g. Lin & Papaloizou 1986; Nelson et al. 2000).  Previous studies of the interaction of  planets with circumbinary discs (Nelson 2003; Pierens & Nelson 2007, 2008a,b, 2013; Kley & Haghighipour 2014, 2015) have indeed shown that migrating circumbinary planets naturally park at the edge of the central cavity formed by the binary. 
 
The fact that the circumbinary planets detected so far evolve on orbits that are coplanar with the binary is at first sight consistent with the general expectation that the angular momentum vectors of the circumbinary disc and binary should be primordially aligned because both originate from a single cloud. We note, however, that in the context of the Kepler mission, it is of course easier to detect planets whose orbits are aligned with the binary (Winn et al. 2015), since even a small misalignment can prevent the planet from transiting for long periods of time (precession does, however, lead to times when even a misaligned planet transits one or both stars). Hence, the lack of misaligned circumbinary planets may simply be due to selection bias. A primordial misalignment between the disc midplane and the binary orbit, that could in principle lead to the formation of a misaligned circumbinary planet, might arise due to turbulence the cloud (Padoan & Nordlund 2002),  leading to infall of gas whose angular momentum vector varies with time (Bate et al. 2003). This would  be consistent with the observation of inclined circumbinary discs, such as around the pre-main-sequence star binary KH 15D (Chiang et al. 2004), or the binary protostar IRS 43 (Brinch et al. 2016). 

The misalignment between the circumbinary disc and orbit plane of the binary leads to disc warping (Larwood & Papaloizou 1997), to a degree that depends strongly on the disc aspect ratio, h, and viscosity parameter, α. In the case where  h>α and the warp propagation time is shorter than the differential precession timescale, the disc undergoes rigid precession with a small warp, with a precession frequency that tends to  zero if the disc is arbitrarily large (Larwood & Papaloizou 1997). Stronger warps develop in the diffusive regime for warp propagation (i.e. when α>h) and/or if the warp propagation timescale is longer than the differential precession timescale (Larwood et al 1996; Fragner & Nelson 2010). 

Such a disc warp may  have a significant impact on the inclination evolution of embedded planets. Terquem (2013) has shown that secular interaction with a warped disc can make the planetary orbit precess, and  can even lead to very high inclinations when the warp is strong enough. For a disc that is inclined relative to a binary companion, the outcome of the secular interaction between the disc and an embedded planet depends mainly on the disc mass and the ability for the planet to open a gap in the disc. Coplanarity between the disc midplane and the planetary orbit can be maintained if the disc is massive enough (Xiang-Gruess & Papaloizou 2014; Lubow & Martin 2016),  whereas the planet tends to decouple from the disc in low-mass discs or when it opens a gap in the disc (Picogna & Marzari 2015). In the context of circumbinary discs, Foucart & Lai (2013) have shown that the back-reaction torque exerted by a warped disc on the binary tends to align the disc midplane with the binary. Because the alignment timescale increases with the binary semi-major axis,  they suggest that wide binaries with separations ≳ 1 AU can have misaligned circumbinary discs whereas close binaries should be aligned with their circumbinary disc. The implication is that, provided they formed in the late phase of the circumbinary disc, the circumbinary planets should evolve on orbits which are aligned with the binary, consistently with what is currently observed.

In this paper, we employ 3D hydrodynamical simulations to investigate whether or not, for the Kepler-413b  circumbinary planet, the small misalignment between the planetary and binary orbits could have resulted from migration in a circumbinary disc whose midplane is initially misaligned with respect to the binary orbit. The first issue that we address is the hydrodynamic response of the disc to the precessional torque exerted by the binary, as a function of disc parameters and binary inclination. We show that the disc tends towards a quasi-steady state displaying essentially no precession, and generally develops a small warp (i.e. one in which the rate of change of the disc inclination with radius is small), such that the outer disc remains in its original orbit plane, but the inner disc settles towards a plane that is modestly inclined with respect to the binary plane. We use a subset of these relaxed disc models as initial conditions for simulations that examine the orbital evolution of planets that are subject to gravitational interactions with the disc and binary. We find that coplanarity between the planet orbit and the disc midplane is maintained for high mass discs. As the disc mass is decreased, however, the planet decouples gravitationally from the disc and becomes dominated by the binary, which drives precession out of the disc midplane. Tidal interaction between the disk and planet leads to damping of the planet inclination, and the planet evolves towards coplanarity with the binary, as this provides a plane of symmetry for the system. Hence, we suggest that as an inclined circumbinary disc containing a planet on a short period evolves, such that its mass decreases due to accretion onto the central stars or perhaps through mass loss in a wind, a time will be reached when the planet will start to align with the binary. Given sufficient time, it will settle into the binary plane. Hence, even misaligned circumbinary discs may give rise to short period planets on orbits that are coplanar with the binary. 

This paper is organized as follows. In Section 2 we describe the hydrodynamical model and the initial conditions that are used in the simulations. In Section 3 we present the results of simulations of circumbinary discs that are modestly misaligned with respect to the binary orbit. The orbital evolution of embedded planets that are initially coplanar with the disc is then studied in Section 4. Finally, we discuss our results and draw our conclusions in Section 5.



§ THE HYDRODYNAMIC MODEL



 §.§ Numerical setup

The simulations presented in this paper were performed using FARGO3D (Benitez-Lamblay & Masset 2016). This code uses finite differences, and an advection scheme based on the monotonic transport algorithm (Van Leer 1977).  The FARGO algorithm (Masset 2000) is employed to avoid time step limitations due to the Keplerian velocity at the inner edge of the disc.
 
We adopt locally isothermal disc models, and solve the following equations for the conservation of mass and momentum in spherical coordinates (r,θ,ϕ) with the origin of the frame located at the centre of mass of the binary,  
where r is the radial distance from the origin, θ is the polar angle measured from the z-axis and ϕ is the 
azimuthal coordinate starting from the x-axis:
 
    ∂ρ/∂ t+∇·(ρ v)=0,

 
    ρ(∂ v/∂ t+ v·∇ v)=-∇ P-ρ∇Φ +∇·Π,

where ρ is the density, v the velocity, P the pressure, Π the viscous stress tensor and Φ is the gravitational potential, which can be written as:
 
    Φ=-GM_1/|𝐫-𝐫_1|-GM_2/|𝐫-𝐫_2|,

where M_1, M_2 are the masses, and r_1, r_2  the radius vectors of the primary and secondary stars respectively. We do not include disc self-gravity in this work. Previous work (Mutter et al. 2017) for flat 2D discs showed that for disc masses ≲ 5 MMSN (where MMSN refers to the minimum mass solar nebula model of Hayashi (1981)), the disc structure is only weakly affected by self-gravity. Here, our reference disc mass corresponds to 1 MMSN, so development of the disc structures reported by Mutter et al.(2017) should not arise. Self-gravity, however, can lead to changes in the precession rates of discs, and this could in principle modify the disc structure. We leave exploration of this effect to future work.

When a planetary companion is included, it contributes to the gravitational potential through the expression 
    Φ_p=-Gm_p/(|𝐫-𝐫_𝐩|^2+b^2)^1/2

where m_p is the planet mass, r_p is the planet radius vector and b is the smoothing length which corresponds approximately to the cell diagonal size. 

Computational units are chosen such that  the total mass of the binary is M_ bin=1, the gravitational constant G=1, and the radius r=1 in the computational domain corresponds to the binary semi-major axis for the  Kepler-413 system (a_bin∼ 0.1 AU, see Table 1). When presenting simulation results, we use the binary orbital period T_bin=2π√(a_bin^3/GM_bin) as the unit of time. The computational domain extends from R_in=1.5a_bin to R_out=80a_bin. Here, employing a large disc model is important  because the disc precession frequency strongly depends on the value for R_out. For an effectively infinite disc that is inclined with respect to the binary orbital planet, we would expect  the disc to maintain a steady warped structure without any global precession, and which globally aligns to the binary orbit on a viscous time scale. Because of the large disc radial extent, we employ a logarithmic grid using 574 radial grid cells. In the azimuthal direction the simulation domain extends from 0 to 2π with 700 uniformly spaced grid cells.  In the meridional direction, the simulation domain covers  E[5+γ_bin/h_max] disc pressure scale heights above and below the disc midplane, where E[x] denotes the nearest integer to the real x,  γ_bin is the binary inclination and h_max is the aspect ratio at the outer edge of the disc.



 §.§ Initial conditions

In this work, we focus on locally isothermal discs for which the aspect ratio does not vary with time and is given by h=h_0(r/r_0)^f,  where f is the disc flaring index and r_0=1. Most of the models have constant 
aspect ratio with h_0=0.05, but we also considered a flared disc model with f=0.25 and h_0=0.01. This particular model was designed in a such a way that a planet with mass equal to that of Kepler-413b (see Table 1) would be able to carve a gap in the disc when located at a_p∼ 3.5 a_bin.

The effective kinematic viscosity ν is modelled using the standard alpha prescription  ν=α c_s H (Shakura & Sunyaev 1973) where  c_s 
the sound speed and α the viscosity parameter. In this work we explore how changing α modifies the disc structure, and adopt the values α=10^-5, 4× 10^-3 and 10^-1.

We consider circumbinary disc whose midplane start off in the equatorial plane of the computational grid, with a binary whose orbital plane is inclined.  The initial density and azimuthal velocity profiles satisfy hydrostatic equilibrium and are given by:

    ρ(r,θ)=f_gapρ_0(r/r_0)^-ξ[sin (θ) ]^1/h^2+2f-(ξ+1)

and 

    v_ϕ(r,θ)=[GM/r(1+(2f-(ξ+1))h^2)]^1/2,


 with ρ_0  the density at r_0=1 and with ξ=1+p+f, where p=3/2 and β=1-2f are the power-law indexes for the surface density and temperature, respectively. ρ_0 is set to a value corresponding to a 1 MMSN disc, namely containing 2% of the mass of the binary within 30 AU. f_gap is a gap-function used to initiated the disc with an inner cavity (created by the binary), and is given by:

    f_gap=(1+exp[-r-r_gap/0.1r_gap])^-1,

 where r_gap=2.5a_bin is the analytically estimated gap size (Artymowicz & Lubow 1994).
The binary semi-major axis is set to a_bin=1.0 while the binary eccentricity e_bin and mass ratio q_bin are set to match those of Kepler-413b (see Table 1).  Regarding the binary inclination γ_bin, we explored four cases with values  0^∘≤γ_bin≤ 8^∘. In each case, the ascending node of the binary is located on the x-axis. The gravitational back-reaction from the disc onto the binary is not included, so the binary orbit remains fixed. Including the effect of the torque exerted by the disc onto the binary would be expected to lead to a decrease of the binary semi-major axis and evolution of its eccentricity (Pierens & Nelson 2007; Fleming & Quinn 2017). On time scales that depend on the disc mass, alignment of the orbital angular momentum of the binary with that of the disc would also arise (e.g. Foucart et Lai 2013). In total, we considered a set of 7 different disc models that are characterized by their values for α, γ_bin, h_0, f. Details of parameters for each run can be found in Table 2.







 §.§ Boundary conditions

We employ a closed boundary condition at the outer edge of the disc, and a viscous outflow condition (Pierens & Nelson 2008) at the inner edge, for which  the radial velocity in the ghost zones is set to  v_r=β v_r(R_in), where v_r(R_in)=-3ν/2R_in is the 
gas drift velocity due to viscous evolution and β is a free parameter which is set to β=5.  In a recent study, Mutter et al. (2016) showed that the radial structure of a circumbinary disc (e.g. surface density profile, tidally truncated cavity size) is robust regarding the choice of the inner boundary condition only when the inner edge of the computational domain is small enough for the secondary star to be completely embedded within the computational domain. The computational expense incurred by 3D hydrodynamic simulations prevents us from  adopting such a small value for the inner radius because of the time step constraint that it would impose, so the structure of the tidally truncated cavity may be impacted by our  choice for the inner boundary condition (Kley & Haghighipour 2014; Mutter et al. 2016). This issue is clearly of importance if trying to predict the final stopping location of a migrating planet, but in this study we are focussing on how the planet inclination evolution varies with disk mass. The precise form of the radial density distribution should not be so important when addressing that issue. Outflow boundary conditions are used at the meridional boundaries, for which all quantities in the ghost zones have the same values as in the first active zones, except the meridional velocity whose value is set to 0 if it is directed towards the disc midplane to prevent inflow of material.



§ WARP EVOLUTION IN CIRCUMBINARY DISCS



We now describe the warp evolution for each of the misaligned circumbinary disc models we considered.  In the following, the inclination angle, i_d, between the angular momentum vectors of disc annuli located at radius r and the binary is defined as:

    cos(i_d(r))= L_D· L_B/| L_B|| L_D|

where L_B is the binary orbital angular momentum vector and
 L_D is the angular momentum of disc annuli located at radius r:
 
    L_D=∫_0^2π∫_π/2-Δθ^π/2+Δθ l r^2sin(θ)drdθ dϕ,

 where dr is the radial grid spacing and l=ρ ( r× v) the angular momentum density. Similarly, we define the precession angle 
 of the disc angular momentum vector around L_B as:
 
    cos(β_d(r))= L_D× L_B/| L_B× L_D|· u,

 where u is a reference unit vector in the orbital plane of the binary, defined to be u= e_x. 
 
A disc annulus located at r undergoes precession at a frequency given approximately by the 
free particle precession rate due to the torque exerted by the binary (e.g. Nixon et al. 2013):

    Ω_prec(r)=-3η/4a_b^2/r^2Ω(r)cos(i_D(r))

where η=M_1M_2/M_bin^2 and Ω is the Keplerian frequency. Differential precession of neighbouring annuli induces a twist that causes their midplanes to become misaligned.  This generates pressure-induced horizontal shear motions which, in the case with h>α,  are responsible for bending waves propagating  across the disc at speed c_s/2. In the case h<α, however, the horizontal shear resulting from misalignment of disc midplanes is damped locally and a warp that is built up in the disc will evolve diffusively in that case. 




 §.§ Evolution for a fiducial run: Model 1


Figure <ref> displays the evolution of the disc inclination and precession angles, i_D and β_D,  as a function of radius for Model 1 at  equally spaced times t=1, 100, 200, ...  up to the final run time t=10^3 T_bin.  Warp evolution takes place in the wave-like regime since  α<h. It is clear that the system evolves towards a quasi-steady state with very small warp (i.e. changes in inclination occur smoothly over a large range of radii), due to bending waves communicating information about the local warp structure efficiently across the disc. At t=10^3   T_bin, we see that an equilibrium has been reached at r ∼ 3.5 a_bin, the current location of Kepler-413b,  while the outer parts of the disc are still evolving because of the very long differential precession times there. It is interesting to notice that at this location, the inclination angle at equilibrium i_D∼ 2.5^∘ is very close to the orbital misalignment of Kepler-413b relative to the binary orbital plane (see Table 1). 

Here, we note that the fact that the disc is distorted up to 
∼ 50  a_bin at t=10^3   T_bin is not  a consequence of the bending waves propagating from the inner disc but  rather simply results from disc precession.  Given that the warp propagates at speed c_s/2, we would indeed expect the outer regions of the disc located at ∼ 50  a_bin to respond to the changing inclination of the inner disc, due to transient bending waves, after  a time longer than ∼ 10^3   T_bin.   In fact, our results  suggest that  only a small amount of precession is required for the disc to respond by setting up internal stresses to ensure uniform precession. The consequence is that  on a timescale corresponding to the differential precession time in the outer disc, we would expect a steady-state warped structure to develop across the disc, where the binary torques that try to twist up and warp the disc are counter balanced by internal stresses involving pressure and viscosity. On even longer timescales, we would expect the disc and binary orbit planes to align. For a fixed binary orbit, as considered here, where realignment occurs because the inner disc brings the outer disc into alignment with the binary orbit plane, this is likely to occur on the global viscous evolution time of the disc once the steady-state warped structure has been established. 
 
 In the wave-like regime for warp propagation, we note that a significant warp may eventually be established in the innermost regions where the differential precession induced by the binary is stronger, provided that the precession  timescale there t_prec∼Ω_prec^-1 is shorter than the warp propagation timescale t_w=2r/c_s. Using Eq. <ref>, we find t_w<t_prec for r > 3.6 a_bin. As discussed later in the paper (see Sect. <ref>), 
this corresponds roughly to the tidal truncation radius  such that we  would expect the whole disc to remain smooth and only be modestly warped, consistently with the simulations.  Also, we see in the bottom panel of Fig. <ref> that the system evolves towards a very slowly precessing configuration with almost no twist (i.e. β_D is ∼ constant), which corresponds to what would be expected for an infinite disc with t_w<t_prec everywhere. In fact, when the latter condition is fullfilled, a disc with outer radius R_out should precess as a rigid body with precession rate given by (Larwood & Papaloizou 1997):
 
    ω_prec=-3/4Gη a_bin^2cos(γ_bin)∫_R_in^R_outΣ r^-2dr/∫_R_in^R_outΩΣ r^3dr,

 which we can write, for a surface density profile Σ∝ r^-3/2:
  
    ω_prec≃ -3/10ηcosγ_bin(a_bin/R_out) (a_bin/R_in)^5/2Ω_bin
  
  For a finite disc with R_out=80a_bin, the previous relation gives  ω_prec∼ 10^-4Ω_bin. The results of the simulations, showing a very small precession rate,  are 
 clearly in very good agreement with such an expectation. These numbers, however, also demonstrate the extreme difficulty of simulating circumbinary systems with realistic outer disc radii over run times that allow the long-term, global, quasi-steady evolution to be examined. Nonetheless, as our simulations demonstrate, the more rapid evolution of the inner disc regions towards a quasi-steady state allows analysis of the dynamics there to be captured at reasonable computational cost.
 


 §.§ Dependence on binary inclination: Models 2, 3, 4


We now examine how changing the binary inclination, γ_bin, influences the evolution. Although not shown here, we find that for Models 2 and 3, that have γ_bin=2.5^∘ and 8^∘, respectively, the inclination and precession angles again reach constant values in the innermost parts of the disc at t∼ 1000 T_bin.  For these models, the inclination and precession angles at equilibrium are plotted as a function of radius  in Fig. <ref>. The weak dependence of  the differential precession rate on the binary inclination (see Eq. <ref>) implies that for the range of values that we consider (2.5≤γ_bin≤ 8^∘), the relaxed structure of the disc should not vary too much from one model to another,  and this is indeed what we observe.  Fig. <ref> reveals that (i)  the final value for the precession angle does not depend on γ_bin, and (ii) a very small warp is established with i_D∼γ_bin/2 at  r=3.5 a_bin which is about the tidal truncation radius.

The azimuthally-averaged surface density profiles at t=1000 T_bin are displayed in Fig. <ref>, and indicate that the size of the tidally cleared cavity and location of the maximum surface density are similar for all models.




For the disc parameters that we adopted,  it is somewhat encouraging that the location of the density peak in each model almost coincides with the current location of Kepler-413b, since it is  expected that the inward migration of a protoplanetary core will be stopped at the tidally truncated cavity  due to the action of corotation torques (Pierens & Nelson 2007). 







 §.§ Dependence on α viscosity: Models 5, 6

Models 5 and 6 both have h=0.05 and γ_bin=2.5^∘, but  the viscous stress parameter α is different in each model. Model 5  has α=10^-4  such that warp propagation occurs via bending waves, whereas Model 6 has α=0.1 so warps propagate diffusively because α>h. The disc inclination and precession angles are shown in Fig. <ref>.  Model 5 develops a very small warp  whose  structure at the end of the simulation is very similar to that of Model 2 (see Fig. <ref>). This  is not surprising given that both models are in the wave regime with the same value for  the warp wave speed  ∼ c_s/2.  As for Model 5, the disc in Model 6 which has a larger viscosity also attains a state with almost no global precession at the end of the simulation,  but is found to be much more twisted, with a value for the twist between the inner and outer parts of the disc of ∼ 90 degrees. This is 
in good agreement with the results of Fragner & Nelson (2010) who studied discs in external binary systems, and found that discs with large viscosity tend to develop significant twists before achieving a state of solid 
body precession. As demonstrated by these authors, this arises because when warp communication is less efficient, the disc needs to be significantly distorted before internal stresses  are able to  counterbalance the effect of the differential precession induced by the binary.  Compared to Models 2 and 5, the amplitude of the warp established in Model 6 is also slightly higher, with the inner and outer regions of the disc having relative inclination of ∼ 1.8 degrees.  This can be also observed in Fig. <ref> which shows column density plots for a particular viewing geometry for Models 5 and 6 at t=700 T_bin. For both models,  we would again expect alignment between the disc and binary orbital plane to occur on a time corresponding to  the global viscous evolution timescale of the disc once a steady warped structure has been established, corresponding to ∼ 4.5× 10^8 and ∼ 4.5× 10^5 T_bin for models 5 and 6, respectively.

 





 §.§ Evolution for a flared disc model: Model 7


We now discuss the circumbinary disc evolution for Model 7, which corresponds to a flared disc with aspect ratio h=0.01(r/r_0)^0.25. The aspect ratio is clearly smaller than for Model 1, such that we expect the tidally truncated inner cavity to be deeper for Model 7 (e.g. Pierens & Nelson 2013). The surface density profile at t=1000 T_bin is shown in the upper panel of Fig. <ref>.  The density peak is located slightly further out at R∼ 3.8 for this run compared to Model 1, consistent with expectations, and the cavity edge is clearly steeper.

 Using Eq . <ref> for the free particle precession rate, we  estimate t_w<t_prec=Ω_prec^-1 in regions where r≳ 6 a_bin, which is slightly larger than the tidal truncation radius. The implication is that we expect a small amount of warping to be established in the innermost disc regions.  For this model, the disc inclination and precession angles are plotted as functions of time and radius in the middle and lower panels of Fig. <ref>, respectively. Similarly to the previous models that we have considered, a quasi-steady configuration is reached at t∼ 1000 T_bin interior to r∼ 40  a_bin.  This quasi-steady state is also characterised by a small degree of warping near R∼ 5-6 a_bin and a small twist, as expected. At t=1000 T_bin, the values at R=3.5  a_bin for the disc inclination and precession angles are  i_D∼ 1.7^∘ and β∼ 300^∘ respectively, which  are close to those that we have obtained in Model 1. The outer edge of where the disc is being distorted is also very similar to what has been found in Model 1. This confirms that this location is set by differential precession rather than by bending waves, since these propagate much  slower in Model 7.  

The main 
 difference here will be that a planet with the same mass as Kepler-413b will open a deep gap in the Model 7 disc, whereas it will at most create a small density depression and  experience Type I migration in Model 1. In the next section, we discuss the dynamical evolution of planets that are initially locally aligned with the circumbinary disc.










§ DYNAMICS OF PLANETS EMBEDDED IN WARPED CIRCUMBINARY DISCS



 §.§ Evolution for two fiducial cases: Models 1 and 7

To examine the evolution of planets that are initially locally aligned with the circumbinary disc, we focussed on the Models 1 and 7, which have h=0.05 and h=0.01(r/r_0)^0.25, respectively. The viscous stress parameter α=4× 10^-3 in both models. We inserted a non-accreting planet with mass ratio 
 q ≡ m_p/M_bin=1.6× 10^-4 and semi-major axis a_p=3.5 a_bin on a circular orbit after each of the discs had evolved for 1000 T_bin. Hence, the planet is located close to the tidal truncation radius at the beginning of this phase of the simulations.
 
We begin discussion of the simulation results by considering discs with masses equivalent to 1 MMSN. We discuss how the results depend  on disc mass in the next section. The left panel of Fig. <ref> shows the evolution of the planet semi-major axis and eccentricity for the two models. The criterion for deep gap formation (Crida et al. 2006) is

    P=1.1(q/h^3)^-1/3+50α h^2/q≤ 1

For Model 1, P∼ 4 such that gap formation should not occur and migration will be in the Type I regime. In that case, previous work has shown that the planet can become trapped at the edge of the inner cavity for two different reasons: i) halting of migration through the action of corotation torques in a region where the surface density increases with radius (Pierens & Nelson 2007); ii) stopping of inward migration once the planet eccentricity becomes high enough to induce a Lindblad torque reversal (Pierens & Nelson 2013). Inspection of Fig. <ref> reveals that migration is halted at the cavity edge, as expected. The horseshoe libration timescale t_lib=8π a_p /3 Ω_p x_s is significantly shorter than the viscous timescale across the horseshoe region t_ν=x_s^2/ν , where x_s=1.2√(q/h^3) is the  half-width of the horseshoe region (Paardekooper et al. 2010), such that corotation torques are partially saturated and therefore significantly attenuated. The implication is that stalling of migration likely occurs because of Lindblad torque reversal. To explain this effect, Pierens & Nelson (2013) showed that a planet on a significantly eccentric orbit has smaller angular velocity at apocentre than neighbouring disc material, resulting in the outer disc exerting a positive torque on the planet. The opposite happens at pericentre where the inner disc exerts a negative torque. Hence, the torque on the planet oscillates from negative to positive values as the planet moves from apocentre to pericentre and back again. The  time averaged total torque, however, can equal zero resulting in migration being halted.  Figure <ref> displays the total torque together with the planet radial position over a few orbital periods, and shows that the torque is indeed positive (resp. negative) when the planet is at apocentre (resp. pericentre), but with the orbit-averaged torque almost cancelling as expected for Lindblad torque reversal. As noted by Pierens & Nelson (2008), the slight phase shift between the two curves is due to the time needed for the planet to create an inner (resp. outer) wake in the disc at apocentre (resp. pericentre). 
 
 Model 7 resulted in the planet migrating inwards over the duration of the run.  This is because the planet opens a gap in the disc. The gap opening criterion given by Eq. (<ref>) gives P∼ 0.6, such that deep gap formation is expected, and is confirmed by   inspection of Fig. <ref> which shows a projection of the column density through the disc for a viewing angle (θ,ϕ)=(70^∘,80^∘) at time t=1000 T_bin after the planet was inserted. The inner disc is depleted so inward migration is driven by the outer disc, pushing the planet into the cavity created by the binary.  Due to the very long run times required by the simulations, the final fate of the system remains uncertain. It is possible that the planet will approach 6:1 mean motion resonance with the binary at later times. Previous work has shown different outcomes for long term simulations involving gap forming planets. These include ejection of the planet from the system through interaction with the binary, scattering into the outer disc by the binary such that the planet repeatedly migrates to the inner cavity and is scattered out again, and slow inwards migration that stalls when the density of gas in the cavity becomes very small (Nelson 2003; Pierens & Nelson 2008; Pierens & Nelson 2013).
 
We notice that for the two models that we have considered here, the planet eccentricity reaches a constant value but we do not obtain particularly good agreement with observational data. Model 1 overestimates the value for e_p (e_p∼ 0.15 at the end of the run) whereas the value 
for e_p inferred from Model 7 (e_p∼ 0.05) is much smaller that the observed value e_p∼ 0.1. We note that high values of the eccentricities are generally driven by non-axisymmetric features in the disc and not by the central binary in these types of simulations (Kley & Nelson 2010). Presumably a model with parameters somewhere between those of Model 1 and Model 7 could produce a better fit.
 
The right panel of Fig. <ref> shows the time evolution of the planet inclination i_p and precession angle  β_p of the planet angular momentum vector around the angular momentum of the binary. These are defined by: 

    cos(i_p)= L_p· L_B/| L_B|| L_p|

and
 
    cos(β_p)= L_p× L_B/| L_B× L_p|· u

where L_p is the orbital angular momentum vector of the planet. 
Since the disc angular 
momentum is much higher than that of the planet, changes in the global disc inclination resulting from 
secular interaction with the planet is expected to be small. Although not shown here, we checked that the 
disc inclination i_D and precession angles β_D remain almost constant over the duration of the simulation (∼ 1000 T_bin). Looking at the right panel of Fig. <ref>, it is immediately obvious that coplanarity between the disc and planet is not maintained for a 1 MMSN disc,
and the planet orbit tends to align with the binary orbit. Since (i_p-i_d)^2 ∝ (β_p-β_d)^2 (Fragner, Nelson & Kley 2011), the nodal precession of the planet induced by the binary, which occurs significantly faster than the disc precession, induces the relative inclination between the planet and local disc midplane to grow at early times,  leading to the initial 
decrease in i_p that can be observed in the upper right panel of Fig. <ref>. i_p then oscillates with period T_osc∼ 300 T_bin which, for r∼3.4, is close to the free particle precession period 2π/Ω_prec(r), where Ω_prec(r) is given by Eq. <ref>. As expected, the relative tilt i_p-i_d (not shown here) is at a minimum when the disc and planet precession angles are in phase whereas it is at a maximum when |β_p-β_d |=π. For both Models 1 and 7, the long term evolution involves the planet inclination relative to the binary progressively decreasing. This occurs because the binary causes the planet to continuously precess relative to the disc. When the planet and disc planes are misaligned the passage of the planet through the disc occurs at high velocity, and the disc provides a source of dynamical friction on the planet's orbit. In the limit that the rapid precession of the planet is driven by the binary, with the disc only making a negligible contribution by virtue of it having a small mass, the binary orbit plane provides a plane of symmetry to which the planetary orbit damps due to the disc-induced dynamical friction.











 
 




 §.§ Dependence on disc mass

To examine how the disc mass M_D affects the inclination evolution of embedded planets, we performed two additional Model 1 simulations with disc masses equivalent to M_D=2 and 20 MMSN.  We note in passing that the model with disc mass equivalent to 20 MMSN is used for illustrative purpose only, since it would be gravitationally unstable outside 
R∼ 60 a_bin. Fig. <ref> shows the planet inclination and precession angles as a function of time for M_D=1, 2, and 20 MMSN.  As expected, the planet tries to maintain coplanarity with the disc as the disc mass is increased, since the gravitational acceleration of the disc tries to pull the planet back towards the midplane, in competition with the binary which tries to drive precession of the planet relative to the disc, and hence away from the disc midplane. For M_D≥ 2 MMSN, we see that the inclination and precession angle exhibit damped oscillations, leading i_p to reach a constant value that increases with disc mass. For M_D=2 MMSN, i_p∼ 1.2^∘ once a quasi steady state is reached whereas for M_D=20 MMSN, the planet remains almost fully aligned with the disc.  For both cases, the precession angle is also observed to librate around a value that is very close to that of the disc (see Fig. <ref>),  whereas it circulates in the case with M_D=1 MMSN. In summary, our results show that in the limit of high disc mass, the precession of the planet relative to the disc is prevented from occurring because disc gravity causes the planet to remain orbiting very close to the disc midplane. In a low mass disc, precession is completely dominated by the central binary, and the damping of the disc-planet relative motion as the planet passes through the disc during periods of high mutual inclination leads to the planet orbit damping towards the binary orbit plane. For intermediate disc masses, the planet orbit achieves a stationary state with a precession angle that is just offset from that of the non-precessing disc, and a non-zero orbital inclination that is smaller than that of the disc.



 §.§ A secular model with inclination damping, migration and disc dispersal

In the context of secular theory, libration of the precession angle is expected for sufficiently massive discs that induce onto the planet a precessional torque that is stronger than the one induced by the binary. The inclination damping effect, however, is due to resonant planet-disc interactions which tend to damp the planet orbit into the local disc plane. For very high mutual inclinations, a model of inclination damping based on dynamical friction can also be formulated. 

The evolution of the disc-planet-binary system under the influence of secular gravitational interactions, migration and inclination damping can examined in a simplified formalism by considering the evolution equation for the planet tilt W_p=l_p,x+i l_p,y including effects resulting from  i) secular interaction with both the binary and the disc, and ii) resonant interaction with the disc. Following Lubow & Martin (2016), the secular evolution equation for the planet  is given by:

    J_pdW_p/dt=i C_pd(W_d-W_p)-iC_ps W_p,

where J_p=m_pa_p^2Ω_p is the planet angular momentum, W_D=l_d,x+i l_d,y is the disc tilt,  and C_pd (resp. C_ps) is a coupling coefficient between the planet and the disc (resp. binary).  These are given by:

    C_pd=2π∫_R_in^R_outGm_pRΣ(R)K(R,a_p)dR

and

    C_ps=Gm_pM_2K(a_p,a_bin),

where the kernel K is given by:

    K(R_i,R_j)=R_iR_j/4π∫_0^2πcos(ϕ)dϕ/(R_i^2+R_j^2-2R_iR_jcos(ϕ))^3/2.

Regarding the evolution equation for the disc, we only include the term corresponding to the secular interaction with the planet so that it reads: 

    J_pdW_d/dt=i C_pd(W_p-W_d)

where

    J_d=2π∫_R_in^R_outΣ(R)R^3Ω(R) dR

is the disc angular momentum. The effect of the planet on the disc is expected to be small because of the small angular momentum of the planet compared to the disc, but is included for completeness. The effect of the binary on a disc with large radius is also negligible in terms of driving its precession, as demonstrated by our simulations.

The effect of the resonant interaction with the disc 
can be incorporated by adding an additional term  (Lubow & Ogilvie 2001) to the secular evolution equation for the planet, which then becomes:

    J_pdW_p/dt=i C_pd(W_d-W_p)-iC_ps W_p-J_pW_p-W_d/τ_inc

where τ_inc is the inclination damping timescale which is set to (Tanaka & Ward 2004):

    τ_inc=t_wave/0.544

with:

    t_wave=q^-1Σ_p a_p^2/Mh^4Ω_p^-1

where Σ_p is the disc surface density at the location of the planet. We solved the secular evolution equations for the disc and planet,  assuming that the disc and planet are initially coplanar 
with i_p=i_d=2.5^∘ and β_p=β_d=0^∘. We  adopted disc parameters corresponding to the initial disc parameters for Model 1, and used the surface density profile at equilibrium to compute the initial planet–disc coupling coefficient C_pd in Eq. <ref>. To model gas disc dissipation due to accretion onto the star and photoevaporation, this coefficient and the disc angular momentum J_D are forced to decay exponentially with an e-folding time t_dis, for which we adopt a reference value of  10^5 T_bin. The planet is also allowed to migrate on a timescale t_m=J_p/2J̇_p, where J̇_p is the disc torque acting on the planet which is given by the sum of the Lindblad torque and fully unsaturated horseshoe drag. In that case, J̇_p is given by (Paardekooper et al. 2010):

    J̇_p/Γ_0=-(2.5-0.5β-0.1σ)(0.4/b/h)^0.71-1.4β(0.4/b/h)^1.26+
    
    1.1
    (3/2-σ)(0.4/b/h)

where β=1 and σ are the negatives of the power-law indices for the temperature and surface density, respectively, and where:

    Γ_0=(q/h)^2Σ_pa_p^4Ω_p^2


In order to investigate how this secular model  compares with the results of the previously presented hydrodynamical simulations, we first consider the case with  t_dis=10^8 T_bin so that the disc mass remains almost constant over the considered evolution times of
10^5 T_bin, and we also assume a_p=3.45 initially such that the migration is almost stalled due to the action of the corotation torque.
Fig. <ref> shows the planet inclination i_p as a function of the precession angle β_p (the phase portrait) for initial disc masses ranging from 0.1 to 20 MMSN. The upper panel corresponds to the phase portrait that is 
obtained when the resonant term is not taken into account. For M_D≤ 2 MMSN, the precession angle  β_p circulates, while it librates for M_D=20 MMSN, resulting in a very small change in the planet inclination. This is consistent with the expectation that the planet precession becomes locked to the disc precession for massive discs (Fragner et al 2011; Lubow & Martin 2016). Not surprisingly, the change in i_p is also  modest for very  low-mass discs, and the precessional torque 
exerted by the binary makes β_p circulate in that case with precession frequency given by Eq. <ref>.

The lower panel of Fig. <ref> corresponds to the case where resonant effects are considered. The oscillations 
for the tilt and precession angles are clearly damped in that case, consistently with what is observed in 
Fig. <ref>.  For M_D≥2 MMSN, the final values for i_p and β_p are in very good agreement with the  results of hydrodynamical simulations. For M_D=1 MMSN, however, the hydrodynamical simulation shows circulation 
of the precession angle and continuous decrease in the planet inclination whereas in the  lower panel 
of Fig. <ref>, the precession angle is already librating due to damping effects after one synodic precession period. Nevertheless, the case corresponding to M_D=0.1 MMSN confirms that for sufficiently low mass discs, the resonant interaction with the disc causes the planet orbit to become almost perfectly aligned with respect to the binary orbit. 

The time evolution of the planet semi-major axis, a_p, and inclination, i_p, for a disc dissipation timescale of  t_dis=10^5 T_bin and a_p=6 × a_bin initially is presented in Fig. <ref> to illustrate how the migration of a planet in an inclined, dissipating disc around a central binary might evolve.  The idea here is not to fit any particular system, but to illustrate general effects.  As expected, coplanarity with the disc is  approximately maintained at early times, although there is a slight decrease in the relative planet-disc tilt due to the increase in the precessional torque exerted by the binary as the planet migrates inward. The planet is then observed to reach the edge of the truncated cavity at t∼ 10^4 T_bin where it remains trapped, while the planet–binary relative inclination decreases as the disc disperses. At t=10^5 T_bin, the disc mass is M_D∼ 1.8 MMSN and i_p∼ 1.2^∘, which is in good agreement with the results  of the hydrodynamical simulations (see Fig. <ref>). The subsequent evolution involves  continuous damping of the planet towards the binary plane as the disc mass decreases, which is again consistent with the results of both the hydrodynamical simulations and the secular model  with 
t_dis=10^8 T_bin (see Fig. <ref>).


 
  
 
   
 
 

§ DISCUSSION AND CONCLUSION

We have presented the results of 3D hydrodynamical simulations that examine the migration and orbital evolution of planets orbiting around a system with stellar and planetary parameters intended to match the Kepler-413 circumbinary planet system.  We began our study by considering the evolution of circumbinary discs with large outer radii in the absence of embedded planets, assuming that the disc equatorial planes were initially moderately misaligned with respect to the orbital plane of the binary.  Our aim was to examine the effect of the tidal and
precessional torque exerted by the binary on the disc structure as a function of disc parameters. Most of our models have constant aspect ratio h=0.05, and we studied how the disc response varies as a function of the binary inclination γ_bin and viscosity parameter α. Binary inclinations of γ_bin=2.5, 5, 8 ^∘  were studied and we considered values for the viscosity parameter of α=10^-4, 4× 10^-3, 10^-2. 

Our results confirm that in the wave-like regime for warp propagation, corresponding to h>α, a small slowly-varying warp is set up in the disc due to bending waves communicating efficiently across the full radius of the disc.  The inner regions of the disc models are torqued towards alignment with the binary (with final tilt angles here being about half of the tilt at the outer edge of the disc). In that case, the disc achieves a quasi-steady state in the inner regions with r ≲ 10 a_bin, where a_bin is the binary semi-major axis, after ∼ 1000 binary orbital periods.  This quasi-equilibrium state is characterized by the disc maintaining a very small warp and twist, in addition to exhibiting almost no precession due to its large radial extent. 
The inner regions of the discs are tidally truncated, such that a low density inner cavity is formed, and are torqued towards alignment with the binary, with final tilt angles here being about half of the tilt at the outer edge of the disc. The disc structure is found to depend only weakly on the binary inclination. An equilibrium configuration is also achieved in the diffusive warp propagation regime with α> h. In agreement with previous work (e.g. Larwood et al. 1996; Fragner & Nelson 2010), we find that  in this regime the disc tends to be more twisted, simply because less efficient warp propagation requires a more distorted disc structure to be established before internal stresses counterbalance the effect of the binary precessional torque. 

Having obtained quasi-equilibrium (or at least slowly varying) disc configurations, we then examined the dynamical evolution of planets that are initially embedded in the disc on orbits that are locally coplanar with the disc midplane and close the the outer edge of the tidally truncated inner cavity. We adopted a planet-binary mass ratio equivalent to that of Kepler-413b. We first considered a circumbinary disc model with constant aspect ratio h=0.05, viscosity parameter α=4× 10^-3 and whose inclination with respect to the orbital plane of the binary is ∼ 2.5^∘ at the initial position of the planet (the outer edge of the disc was inclined by 5^∘). Consistently with previous work, we find stalling of migration once the eccentricity of the planet becomes high enough to induce a Lindblad torque reversal. For a disc mass M_d equivalent to 1 MMSN, we also find that coplanarity with the circumbinary disc is not maintained. The rapid nodal precession induced by the binary causes the planet orbit to precess with respect to the disc, and hence to periodically become highly inclined with respect to it as the relative precession angle reaches π. Damping of the relative motion due to dynamical friction (or equivalently due to inclination damping due to resonant disc-planet interactions) causes the planet orbit to settle towards the binary orbit plane.

For  higher disc masses with M_d≥ 2 MMSN, however, the disc gravity starts to become important and changes the evolution dramatically. For M_d =20 MMSN, the disc gravity is strong enough to force the planet to remain orbiting very close to the disc midplane. For M_d =2 MMSN, after a transient period during which the inclination and precession angles of the planetary orbit are observed to undergo damped oscillations, the planet settles into an orbit with a nodal precession angle very similar to that of the local disc annuli, and with an inclination angle that is smaller that the disc, but not fully aligned with the binary orbit plane. 

We also considered a flared disc model with a small aspect ratio h=0.01 at the initial position of the planet, for which the planet is expected to open a gap. This resulted in a similar mode of evolution for i_p for a disc mass of 1 MMSN. Because of the gap formation by the planet,  the orbital evolution was different in that case but was in good agreement with previous work, showing continuous inward migration deep inside the inner cavity until the planet enters a MMR with the binary, or the amount of gas in the vicinity of the planet becomes too small to make the planet migrate further in. 
 
Previous work has demonstrated the plausibility of a formation scenario for the Kepler circumbinary planets involving core formation far away from the binary, followed by migration to the edge of the tidally truncated cavity. Gas accretion could have occurred during the migration or after the planet has reached the cavity edge. In the context of this scenario, our results indicate that it is likely that a core formed at large distance from the binary would maintain coplanarity with the outer disc due to the weak precessional torque exerted by the binary. This coplanarity, however, would ultimately be lost once the disc mass starts to significantly decrease during one of the following scenarios: i) during disc dissipation after the planet has reached the edge of the cavity; ii) in the course of migration if the core formed in a relatively low-mass disc. 

The main result of this paper is that when we consider the long term evolution of misaligned circumbinary discs with embedded planets, during which the disc mass decreases due to accretion onto the star and/or because a photoevaporative or magneto-centrifugal wind removes mass, the planet will tend to align with the binary orbit plane even if the disc-binary system remains misaligned up until the end of the disc life time. Hence, it is not necessary for the disc to align with the binary in order for the planet to lie within the binary orbit plane. The picture presented in this paper may therefore represent an alternative scenario to that presented in Foucart & Lai (2013), who suggest that realignment of binaries and circumbinary discs should occur rapidly. If our results are correct, and other processes do not intercede to misalign circumbinary planets, then these planets should have orbits that are preferentially close to being aligned with the binaries, and far from being isotropically distributed. This has important implications for assessing the occurrence rates of circumbinary planets, since they have been estimated to be more common than planets around single stars in their inclination distribution is close to being isotropic (Armstrong et al. 2014). We have shown that for a single planet orbiting a pair of stars in a circumbinary disc, even if the disc inclinations are themselves isotropic (an unlikely scenario perhaps), the planetary orbits will not be. 

 Regarding the Kepler-413 and Kepler-453 systems, which contain planets with modest inclinations estimated to be ∼ 2.5^∘, it is possible that rapid gas disc removal may have prevented these systems from fully aligning with the central binary.  Other mechanisms, such as planet-planet scattering (Rasio & Ford 1996; Chatterjee et al. 2008) or nodal precession due to interaction with a distant inclined companion might also be invoked to interpret the slight misalignments observed.  We also note that we have only considered
circumbinary discs that are slightly misaligned with respect to the binary orbit, and the case of higher disc-binary tilts needs to be investigated. For larger misalignment angles i≳ 20^∘, tilt oscillations between the planet and the disc can induce the inclination between the planet and the disc to increase above 
  40^∘ because the planet precesses with respect to the disc, resulting in the formation a Kozai-Lidov planet (Martin et al. 2016). Considering larger binary inclinations will allow us to examine this and other effects in the contex of circumbinary planets, and will be the subjects of future work. 


§ ACKNOWLEDGMENTS

Computer time for this study was provided by the computing facilities MCIA (Mésocentre de Calcul Intensif Aquitain) of the Universite de Bordeaux and by HPC resources of Cines under the allocation A0010406957 made by GENCI (Grand Equipement National de Calcul Intensif).  

[Armstrong et al.(2014)]2014MNRAS.444.1873A Armstrong, D. J., Osborn, H. P., Brown, D. J. A., et al. 2014, MNRAS, 444, 1873
[Artymowicz & Lubow(1994)]1994ApJ...421..651A Artymowicz, P., & Lubow, S. H. 1994, ApJ, 421, 651 
[Bate et al.(2003)]2003MNRAS.339..577B Bate, M. R., Bonnell, I. A., & Bromm, V. 2003, MNRAS, 339, 577 
[Batygin(2012)]2012Natur.491..418B Batygin, K. 2012, Nature, 491, 418
[Benítez-Llambay & Masset(2016)]2016ApJS..223...11B Benítez-Llambay, P., & Masset, F. S. 2016, ApJS, 223, 11
[Brinch et al.(2016)]2016ApJ...830L..16B Brinch, C., Jørgensen, J. K., Hogerheijde, M. R., Nelson, R. P., & Gressel, O. 2016, ApJL, 830, L16
[Bromley & Kenyon(2015)]2015ApJ...806...98B Bromley, B. C., & Kenyon, S. J. 2015, ApJ, 806, 98
[Chatterjee et al.(2008)]2008ApJ...686..580C Chatterjee, S., Ford, E. B., Matsumura, S., & Rasio, F. A. 2008, ApJ, 686, 580-602
[Chiang & Murray-Clay(2004)]2004ApJ...607..913C Chiang, E. I., & Murray-Clay, R. A. 2004, ApJ, 607, 913
[Crida et al.(2006)]2006Icar..181..587C Crida, A., Morbidelli, A., & Masset, F. 2006, Icarus, 181, 587
[Fleming & Quinn(2017)]2017MNRAS.464.3343F Fleming, D. P., & Quinn, T. R. 2017, MNRAS, 464, 3343
[Foucart & Lai(2013)]2013ApJ...764..106F Foucart, F., & Lai, D. 2013, ApJ, 764, 106
[Fragner & Nelson(2010)]2010A   A...511A..77F Fragner, M. M., & Nelson, R. P. 2010, A& A, 511, A77 
[Fragner et al.(2011)]2011A   A...528A..40F Fragner, M. M., Nelson, R. P., & Kley, W. 2011, A&A, 528, A40
[Kley & Nelson(2010)]2010ASSL..366..135K Kley, W., & Nelson, R. P. 2010, Planets in Binary Star Systems, 366, 135 
[Kley & Haghighipour(2014)]2014A   A...564A..72K Kley, W., & Haghighipour, N. 2014, A& A, 564, A72
 [Kley & Haghighipour(2015)]2015A   A...581A..20K Kley, W., & Haghighipour, N. 2015, A& A, 581, A20
[Kostov et al.(2014)]2014ApJ...784...14K Kostov, V. B., McCullough, P. R., Carter, J. A., et al. 2014, ApJ, 784, 14
[Larwood et al.(1996)]1996MNRAS.282..597L Larwood, J. D., Nelson, R. P., Papaloizou, J. C. B., & Terquem, C. 1996, MNRAS, 282, 597 
[Larwood & Papaloizou(1997)]1997MNRAS.285..288L Larwood, J. D., & Papaloizou, J. C. B. 1997, MNRAS, 285, 288
[Lin & Papaloizou(1986)]1986ApJ...309..846L Lin, D. N. C., & Papaloizou, J. 1986, ApJ, 309, 846
[Lines et al.(2014)]2014ApJ...782L..11L Lines, S., Leinhardt, Z. M., Paardekooper, S., Baruteau, C., & Thebault, P. 2014, ApJL, 782, L11
[Lubow & Ogilvie(2001)]2001ApJ...560..997L Lubow, S. H., & Ogilvie, G. I. 2001, ApJ, 560, 997
[Lubow et al.(2015)]2015ApJ...800...96L Lubow, S. H., Martin, R. G., & Nixon, C. 2015, ApJ, 800, 96 
[Lubow & Martin(2016)]2016ApJ...817...30L Lubow, S. H., & Martin, R. G. 2016, ApJ, 817, 30
[Masset(2000)]2000A   AS..141..165M Masset, F. 2000, A&AS, 141, 165
[Masset et al.(2006)]2006ApJ...642..478M Masset, F. S., Morbidelli, A., Crida, A., & Ferreira, J. 2006, ApJ, 642, 478
[Marzari et al.(2008)]2008ApJ...681.1599M Marzari, F., Thébault, P., & Scholl, H. 2008, ApJ, 681, 1599-1608
[Meschiari(2012)]2012ApJ...752...71M Meschiari, S. 2012, ApJ, 752, 71
[Meschiari(2012)]2012ApJ...761L...7M Meschiari, S. 2012, ApJL, 761, L7
[Mutter et al.(2017)]2017MNRAS.465.4735M Mutter, M. M., Pierens, A., & Nelson, R. P. 2017, MNRAS, 465, 4735
[Nelson et al.(2000)]2000MNRAS.318...18N Nelson, R. P., Papaloizou, J. C. B., Masset, F., & Kley, W. 2000, MNRAS, 318, 18
[Nixon et al.(2013)]2013MNRAS.434.1946N Nixon, C., King, A., & Price, D. 2013, MNRAS, 434, 1946
[Ogilvie & Latter(2013)]2013MNRAS.433.2420O Ogilvie, G. I., & Latter, H. N. 2013, MNRAS, 433, 2420 
[Padoan & Nordlund(2002)]2002ApJ...576..870P Padoan, P., & Nordlund, Å. 2002, ApJ, 576, 870
[Paardekooper et al.(2010)]2010MNRAS.401.1950P Paardekooper, S.-J., Baruteau, C., Crida, A., & Kley, W. 2010, MNRAS, 401, 1950
[Paardekooper et al.(2012)]2012ApJ...754L..16P Paardekooper, S.-J., Leinhardt, Z. M., Thébault, P., & Baruteau, C. 2012, ApJL, 754, L16
[Picogna & Marzari(2015)]2015A   A...583A.133P Picogna, G., & Marzari, F. 2015, A& A, 583, A133
[Pierens & Nelson(2007)]2007A   A...472..993P Pierens, A., & Nelson, R. P. 2007, A& A, 472, 993
[Pierens & Nelson(2008)]2008A   A...478..939P Pierens, A., & Nelson, R. P. 2008, A& A, 478, 939
[Pierens & Nelson(2008)]2008A   A...482..333P Pierens, A., & Nelson, R. P. 2008, A& A, 482, 333
[Pierens & Nelson(2008)]2008A   A...483..633P Pierens, A., & Nelson, R. P. 2008, A&A, 483, 633
[Pierens & Nelson(2013)]2013A   A...556A.134P Pierens, A., & Nelson, R. P. 2013, A& A, 556, A134
[Rasio & Ford(1996)]1996Sci...274..954R Rasio, F. A., & Ford, E. B. 1996, Science, 274, 954
[Scholl et al.(2007)]2007MNRAS.380.1119S Scholl, H., Marzari, F., & Thébault, P. 2007, MNRAS, 380, 1119
[Shakura & Sunyaev(1973)]1973A   A....24..337S Shakura, N. I., & Sunyaev, R. A. 1973, A& A, 24, 337
[Tanaka et al.(2002)]2002ApJ...565.1257T Tanaka, H., Takeuchi, T., & Ward, W. R. 2002, ApJ, 565, 1257
[Tanaka & Ward(2004)]2004ApJ...602..388T Tanaka, H., & Ward, W. R. 2004, ApJ, 602, 388
[Terquem(2013)]2013MNRAS.435..798T Terquem, C. 2013, MNRAS, 435, 798
[van Leer(1977)]1977JCoPh..23..276V van Leer, B. 1977, Journal 
of Computational Physics, 23, 276 
[Ward(1997)]1997Icar..126..261W Ward, W. R. 1997, Icarus, 126, 261
[Welsh et al.(2015)]2015ApJ...809...26W Welsh, W. F., Orosz, J. A., Short, D. R., et al. 2015, ApJ, 809, 26
[Winn & Fabrycky(2015)]2015ARA   A..53..409W Winn, J. N., & Fabrycky, D. C. 2015, ARAA, 53, 409
[Xiang-Gruess & Papaloizou(2014)]2014MNRAS.440.1179X Xiang-Gruess, M., & Papaloizou, J. C. B. 2014, MNRAS, 440, 1179

