TY - JOUR
AU - Nelson, Richard, P.
AB - Abstract We present the results of hydrodynamic simulations examining migration and growth of planets embedded in self-gravitating circumbinary discs. The binary star parameters are chosen to mimic those of the Kepler-16, -34 and -35 systems; the aim of this study is to examine the role of disc mass in determining the stopping locations of migrating planets at the edge of the cavity created by the central binary. Disc self-gravity can cause significant shrinkage of the cavity for disc masses in excess of 5–10 × the minimum mass solar nebula model. Planets forming early in the disc lifetime can migrate through the disc and stall at locations closer to the central star than is normally the case for lower mass discs, resulting in closer agreement between simulated and observed orbital architecture. The presence of a planet orbiting in the cavity of a massive disc can prevent the cavity size from expanding to the size of a lower mass disc. As the disc mass reduces over long time-scales, this indicates that circumbinary planet systems retain memory of their initial conditions. Our simulations produce planetary orbits in good agreement with Keper-16b without the need for self-gravity; Kepler-34 analogue systems produce wide and highly eccentric cavities, and self-gravity improves the agreement between simulations and data. Kepler-35b is more difficult to explain in detail due to its relatively low mass, which results in the simulated stopping location being at a larger radius than that observed. accretion, accretion discs, hydrodynamics, methods: numerical, planets and satellites: formation, planet–disc interactions 1 INTRODUCTION With the current tally of Kepler circumbinary planets standing at 11 (Kepler-16b,1 Kepler-34b and Kepler-35b,2 Kepler-38b,3 Kepler-47b,c,4 and d,5 Kepler-64b,6 Kepler-413b,7 Kepler-453b8 and Kepler-1647b9), this class of object is one of the most interesting outcomes of the now two decade old search for planets around stars other than the Sun. They form part of an exoplanet catalogue that contains planets with very diverse orbital and physical characteristics orbiting within a broad variety of stellar systems. Interest in circumbinary planets pre-dates their discovery, with their theorized presence prompting a number of studies into their formation and dynamical evolution. One of the key findings from this period, which has been validated by observations in recent years, is the work of Holman & Wiegert (1999). They found a limit for dynamical stability around short-period binaries. This critical limit depends on the mass and orbital properties of the binary. The majority of the Kepler circumbinary planets lie close to this limit within their respective systems. It is assumed that this, as well as the mutual co-planarity of the planet and binary orbital planes, is a fingerprint of the planets in these systems having formed in a common circumbinary disc. Two general scenarios of planet formation theory have been used to try and explain the observed positions of the planets in these systems. The first is that the planets formed in situ from material in close vicinity of the binary. In this case planet-forming material must be brought together under the strong influence of the gravitational field of the binary. Previous studies have shown numerous disruptive effects to the goal of bringing planetesimals together in a manner which results in mass growth: N-body simulations show excitation of planetesimal eccentricity leading to relative velocities disruptive to accretion and the formation of planetary bodies (Paardekooper et al. 2012; Meschiari 2012a,b; Lines et al. 2014; Bromley & Kenyon 2015); differential pericentre-alignment of eccentric planetesimals of different sizes leads to corrosive collisions (Scholl, Marzari & Thébault 2007); and gravitational interactions with asymmetric features in the gas disc, and the global eccentric mode, leading to large impact velocities of planetesimals (Marzari, Thébault & Scholl 2008; Kley & Nelson 2010). The second, and the focus of this series of papers, is that the protoplanetary cores formed in the quiescent exterior of the disc, where the disruptive influence of the binary is negligible, and then moved inwards to their observed position through disc-driven migration – either Type-I (Ward 1997; Tanaka, Taku & Ward 2002) or Type-II (Lin & Papaloizou 1986; Nelson et al. 2000). Whether this period of migration and mass growth occurs early or late in the lifetime of these circumbinary discs is unclear. Whilst this scenario solves the problems with in situ formation we must now answer how the migrating planets stop in the inner disc. Again, it was prior to the discovery of the Kepler circumbinary planets that this question was answered. The influence of the binary exerts a tidal torque on the inner circumbinary disc which sweeps material away from the binary, creating a central cavity that can act as a barrier to migration. The radial extent of this cavity depends on the mass and orbital properties of the binary, as well as disc parameters (Artymowicz & Lubow 1994). More recent work has shown that the interaction of the binary with this feature leads to an asymmetric, eccentric, precessing disc (Pierens & Nelson 2013; Pelupessy & Portegies Zwart 2013; Kley & Haghighipour 2014; Mutter, Pierens & Nelson 2017). The directly imaged circumbinary disc in the GG Tau system shows an inner cavity (Dutrey, Guilloteau & Simon 1994). Paper I in this series (Mutter et al. 2017), examined the impact of disc self-gravity and disc mass on the evolution and structure of circumbinary discs in analogue Kepler-16, -34 and -35 systems. Self-gravity has already been examined in low-mass circumbinary discs (Marzari et al. 2009), where it was discounted as an unimportant factor in the disc evolution. As pointed out in Lines et al. (2015), even at low-mass, disc self-gravity can modify the precession frequencies associated with low-frequency global eccentricity modes (Papaloizou 2002). Disc masses ranging from 1 to 20 × MMSN (where MMSN refers to the minimum mass solar nebula; Hayashi 1981) were examined in Paper I, where we confirmed these findings for low-mass discs. However we found for more massive discs, corresponding to 10 × MMSN and beyond, that self-gravity can dramatically alter the evolution and structure of the disc. In addition to the central eccentric cavity, a series of eccentric modes were found to develop at larger orbital radii in all three of the Kepler circumbinary systems simulated, arising from the tidal field of the binary. Furthermore, we see that the radial size of the cavity decreases in more massive discs. Fig. 1 shows the surface density profiles of discs in the Kepler-16 system, for three disc masses: 1, 10 and 20MMSN once they have reached a pseudo-steady state – i.e. the cavity size evolution has stopped and the eccentricity of the disc oscillates around a constant value. In this paper, we present the results of simulations examining a number of migration and accretion scenarios for planet cores, in the full-set of evolved self-gravitating discs from Paper I, including the massive models which show additional eccentric features. While the long-term aim of our work on circumbinary discs is to produce simulation outcomes that fit the Kepler data, because of the relative simplicity of our disc models in this study we set ourselves the less ambitious goal of examining how disc mass and self-gravity influence the final orbital elements of planets that form and migrate in circumbinary discs. Figure 1. Open in new tabDownload slide Changed plot 2D surface density image contours from self-gravitating circumbinary disc simulations around the Kepler-16 system from Paper I. These are all taken at a point when the disc has reached a pseudo-steady state. The panel on the left shows a typical low-mass (1 MMSN) self-gravitating disc in this system. The middle shows a 10MMSN case, and the right-hand panel shows the most massive 20MMSN model. Note the additional eccentric features in the massive discs. The central eccentric cavity common to circumbinary discs can be seen around the binary in all three cases. The positions of the binary stars are marked with the ‘x’ (primary) and ‘+’ (secondary) symbols. Figure 1. Open in new tabDownload slide Changed plot 2D surface density image contours from self-gravitating circumbinary disc simulations around the Kepler-16 system from Paper I. These are all taken at a point when the disc has reached a pseudo-steady state. The panel on the left shows a typical low-mass (1 MMSN) self-gravitating disc in this system. The middle shows a 10MMSN case, and the right-hand panel shows the most massive 20MMSN model. Note the additional eccentric features in the massive discs. The central eccentric cavity common to circumbinary discs can be seen around the binary in all three cases. The positions of the binary stars are marked with the ‘x’ (primary) and ‘+’ (secondary) symbols. The evolution of giant planets in evolved self-gravitating circumbinary discs has not been studied before, however their counterparts in non-self-gravitating discs have been. The interaction of giant migrating planets with circumbinary discs was first studied by Nelson (2003). This study showed that Jovian-mass planets generally migrated into the central cavity, where they were captured into a 4:1 mean-motion resonance with the binary. These giant planets often underwent close-encounters with the binary, with scattering events ejecting them from the system. Lighter, Saturnian-mass planets underwent stable migration to the disc cavity where they then remained in stable orbits (Pierens & Nelson 2008a). Less massive planets undergo Type-I migration until they are halted at the inner cavity edge by a strong positive co-rotation torque, counteracting the Lindblad torque – see Pierens & Nelson (2007, 2008a,b) for more details. The techniques developed in these works were then applied to a number of the newly discovered Kepler circumbinary systems, in attempts to explain and recreate the orbits of their planets (Pierens & Nelson 2013; Kley & Haghighipour 2014, 2015). Pierens & Nelson (2013), henceforth referred to as PN13, had difficulty recreating both the semimajor axes and eccentricities for the observed planets, with a range of disc parameters, under the assumption of an isothermal equation of state. Kley & Haghighipour (2014, 2015), referred to from now on as KH14 and KH15, included a more realistic equation of state and radiation effects, as well as the role of multiplanet migration and interaction. These works also had difficulty in recreating all the observed properties of the Kepler circumbinary planets, although with a little more success than PN13. Using 3D SPH simulations Dunhill & Alexander (2013) argued the near-circular orbit of Kepler-16b hints that it formed in a massive disc in which the orbit of the planet is heavily damped by the disc. Understanding the physics and parameters which affect the environment in which the planets form and evolve – the circumbinary disc – is key to understanding the final, observed states of these intriguing systems. Our motivation in Paper I and this work is to probe the early dynamical history of circumbinary discs – as we increase the disc mass we effectively examine earlier and earlier times in the system's history. We aim to address the questions: Does a high-mass disc leave a fingerprint on the planet population if circumbinary planets form early? Is this erased by the transition to a low-mass disc as the system evolves? Does the epoch when planets form, accrete gas, and migrate affect their final orbital configuration or mass? Using the results from Paper I as a starting point, we examine the impact of self-gravity and disc mass on migration and accretion scenarios for protoplanetary cores, in systems intended to mimic the Kepler-16, -34 and -35 circumbinary systems. The scenarios are carried out in evolved self-gravitating discs with masses equivalent to 1, 2, 5, 10 and 20 × MMSN, in each of the binary systems. These different disc masses are proxies for different eras in the lifetime of the disc, so we can answer the questions raised above. To simulate the evolution of the system from a high- to low-mass disc state we also carry out simulations where the orbital evolution of the planet is tracked as the disc mass is exponentially dissipated. The outline for this paper is as follows. Section 2 describes the physical model and initial conditions used in our simulations. Section 3 looks at the results of the orbital migration of protoplanetary cores in the whole range of evolved self-gravitating discs from Paper I. Section 4 examines the orbital evolution of gas accreting cores. Section 5 contains results from our investigation into the impact of disc dissipation on planetary core migration, and final halting position. Our results from Sections 3– 5 are summarized and discussed in Section 6. 2 NUMERICAL SETUP In this section, we outline the extensions to our numerical model from Paper I, which deal with the interaction of the planet with the disc. For a full description of the equations pertaining to the evolution of the binary-disc system please refer to Section 2.1 of Paper I. The simulations are conducted in a reference frame based on the centre of mass of the binary system. The stellar orbital elements maintain constant values appropriate to the various Kepler systems that we are studying, and the N-body system comprising the stars and planet is evolved such that the centre of mass of this system is non-accelerating. Given that we are working in a frame centred on the binary centre of mass, an indirect term is required in the equations of motion that accounts for the acceleration of the binary centre of mass, as outlined below. 2.1 Equations of motion 2.1.1 Disc evolution The equations of motion of the binary-disc system are detailed in Section 2.1.1 of Paper I in two-dimensional polar coordinates (R, ϕ) with the origin kept at the centre of mass of the binary. Table 1. Binary and planet parameters. . Kepler-16 . Kepler-34 . Kepler-35 . MA (M⊙) 0.690 1.048 0.888 MB (M⊙) 0.203 1.021 0.809 mp (MJ) 0.333 0.220 0.127 qb = MB/MA 0.294 0.974 0.912 qp = mp/M⋆ 3.54 × 10−4 1.01 × 10−4 7.13 × 10−5 ab (au) 0.224 0.228 0.176 ap (au) 0.705 1.090 0.603 eb 0.159 0.521 0.142 ep 0.007 0.182 0.042 Reference (Doyle et al. 2011) (Welsh et al. 2012) . Kepler-16 . Kepler-34 . Kepler-35 . MA (M⊙) 0.690 1.048 0.888 MB (M⊙) 0.203 1.021 0.809 mp (MJ) 0.333 0.220 0.127 qb = MB/MA 0.294 0.974 0.912 qp = mp/M⋆ 3.54 × 10−4 1.01 × 10−4 7.13 × 10−5 ab (au) 0.224 0.228 0.176 ap (au) 0.705 1.090 0.603 eb 0.159 0.521 0.142 ep 0.007 0.182 0.042 Reference (Doyle et al. 2011) (Welsh et al. 2012) Open in new tab Table 1. Binary and planet parameters. . Kepler-16 . Kepler-34 . Kepler-35 . MA (M⊙) 0.690 1.048 0.888 MB (M⊙) 0.203 1.021 0.809 mp (MJ) 0.333 0.220 0.127 qb = MB/MA 0.294 0.974 0.912 qp = mp/M⋆ 3.54 × 10−4 1.01 × 10−4 7.13 × 10−5 ab (au) 0.224 0.228 0.176 ap (au) 0.705 1.090 0.603 eb 0.159 0.521 0.142 ep 0.007 0.182 0.042 Reference (Doyle et al. 2011) (Welsh et al. 2012) . Kepler-16 . Kepler-34 . Kepler-35 . MA (M⊙) 0.690 1.048 0.888 MB (M⊙) 0.203 1.021 0.809 mp (MJ) 0.333 0.220 0.127 qb = MB/MA 0.294 0.974 0.912 qp = mp/M⋆ 3.54 × 10−4 1.01 × 10−4 7.13 × 10−5 ab (au) 0.224 0.228 0.176 ap (au) 0.705 1.090 0.603 eb 0.159 0.521 0.142 ep 0.007 0.182 0.042 Reference (Doyle et al. 2011) (Welsh et al. 2012) Open in new tab The first extension we make to the system described in Paper I consisting of a close binary system surrounded by a self-gravitating disc is the addition of a massive, interacting planetary core. The planet is free to interact with the disc, and vice versa. This results in the potential felt by the disc (equation 4 in Paper I) having two additional terms, represented by Φp and Φi: \begin{equation} \Phi = \Phi _\mathrm{SG} + \sum _{k=1}^{2}\Phi _{\mathrm{s},k} + \Phi _\mathrm{p} + \Phi _\mathrm{i}. \end{equation} (1) The first two terms in this equation are those created by the disc itself through self-gravity, and the two binary stars (with indices s). Their form is described in equations (5) and (7) in Paper I. The form of the potential created by the planet, of mass mp is as follows: \begin{equation} \Phi _\mathrm{p} = -\frac{Gm_\mathrm{p}}{\sqrt{R^2 + R^2_\mathrm{p} - 2RR_\mathrm{p}\cos \left(\phi - \phi _\mathrm{p}\right) + \epsilon ^2}}. \end{equation} (2) ε is a softening length used to avoid singularities in the calculation of the planet's potential; it takes a value equal to 0.4H in this work, where H is the disc thickness. Readers of Paper I will note that this is the same prescription as used for the smoothing length used for the calculation of the disc's self-gravitating potential. The term represented by Φi is the indirect term resulting from the acceleration of the binary centre of mass due to the gravity of the planet.10 2.1.2 Orbital evolution Table 1 contains the best-fitting observed binary and planetary orbital and mass parameters of the Kepler-16, -34 and -35 circumbinary planetary systems, as quoted in Doyle et al. (2011) and Welsh et al. (2012). One of the eventual goals of this work is to recreate the observed state of the Kepler circumbinary systems. To minimize the initial parameter space of this work, and as noted in Paper I, the binaries’ orbital parameters remain fixed throughout our simulations. The orbital evolution of the binary system is therefore independent of the disc and planet system; we hope to revisit the back-reaction of massive self-gravitating discs on the binary in a later work. In Paper I we discuss the drawbacks to this approach, where in our most massive systems – where the disc mass is comparable to the mass of one of the binary stars contained within a notional radius of 30 au – significant modification of the binary will occur if back-reaction is allowed. See Paper I for a detailed description of these problems, in the initial set-up of our simulations and in the discussion. The equation of motion for the binary stars remains unchanged from equation (8) in Paper I. The equation of motion for a planet of mass mp, interacting with the binary and disc system is as follows: \begin{equation} \frac{\text{d}^2\boldsymbol {R}_\mathrm{p}}{\text{d}t^2} = -\sum _{k=1}^{2}\frac{GM_{\mathrm{s},k}\left(\boldsymbol {R}_\mathrm{p}-\boldsymbol {R}_{\mathrm{s}, k}\right)}{|\boldsymbol {R}_\mathrm{p}-\boldsymbol {R}_{\mathrm{s}, k}|^3} + \boldsymbol {f}_\mathrm{dp} - \boldsymbol {f}_\mathrm{i}, \end{equation} (3) where |$\boldsymbol {f}_\mathrm{dp}$| is the force acting on the planet from the disc, and is given by \begin{equation} \,\boldsymbol {f}_\mathrm{dp} = \int _S\frac{\Sigma (\boldsymbol {R})\mathrm{d}\boldsymbol {R}}{\sqrt{R^2 + R_\mathrm{p}^2 - 2RR_\mathrm{p}\cos \left(\phi - \phi _\mathrm{p}\right) + \epsilon ^2}}. \end{equation} (4) The term |$\boldsymbol {f}_\mathrm{i}$| represents the acceleration of the binary centre of mass by the gravity of the planet. 2.2 Hydrodynamic model The hydrodynamic set-up used in this work follows that in Paper I, and builds on the disc-binary results. The simulation work-load for the results presented here, and in Paper I, was split across two separate numerical codes, after comparing test simulations to verify that the results agreed. These codes were FARGO-ADSG and GENESIS. FARGO-ADSG is an updated version of the widely used FARGO code (Masset 1999), which includes the calculation of disc self-gravity as well as an adiabatic equation of state (Baruteau & Masset 2008a,b). GENESIS uses an advection scheme based on the van Leer (1977) monotonic transport algorithm to solve the disc equations, and contains the FARGO time stepping upgrade, as well as a module to calculate self-gravity. In both codes, the binary and planetary orbits are evolved using a fifth-order Runge–Kutta integrator scheme (Press et al. 1992). These codes were used to run 2D hydrodynamic simulations in the plane of the binary's orbit. The calculations presented here use a grid resolution of NR × Nϕ = 550 × 550 cells. The radial grid-spacing between Rin and Rout is logarithmic, as required by the self-gravity calculations (refer to Table 3 for the values used). This has the added benefit of having a finer grid in the inner region of the disc closest to the binary. The azimuthal grid is equally spaced between [0, 2π]. All the disc models use a kinematic alpha-prescription to model turbulence in the disc (Shakura & Sunyaev 1973), where α = 1 × 10−3, and a constant disc aspect ratio, H/R = 0.05, which gives rise to a locally isothermal set-up. The following computational units are used: the total mass of the binary M* = MA + MB = 1, the gravitational constant G = 1 and the radius R = 1 is equivalent to 1 au. To present the results of simulations we use the binary orbital period, |$P_\mathrm{b}=2\pi \sqrt{GM_\star /a^3_\mathrm{b}}$|⁠, as the unit of time. For the migration scenarios where our giant planet cores are allowed to accrete gas from the disc, we follow the prescription of Kley (1999). Accretion is modelled by removing a fraction of the gas within the planet's Hill sphere, RHill = ap(mp/3M*)1/3, from the disc and adding the equivalent mass to that of the planet. The rate at which gas is removed from the Hill sphere is determined by the accretion time-scale, tacc = f tdyn, where tdyn is the orbital period of the planet, and f is an adjustable factor. For the simulations where there is no accretion this corresponds to f = 0. 2.3 Initial conditions The initial conditions used to set up the simulations of the disc models used in this work are detailed in Section 2.3 of Paper I. This section will instead focus on the procedure used to initialize the planet cores in the migration, gas accretion and disc dissipation scenarios presented here. Fig. 2 summarizes the ‘initial’ conditions the planets are inserted into – the pseudo-steady state, azimuthally averaged surface density profiles for the 1, 2, 5, 10 and 20MMSN models in the Kepler-16, -34 and -35 systems. These snapshots are taken at t = 6000 Pb, once the disc has reached a pseudo-equilibrium. Figure 2. Open in new tabDownload slide Azimuthally averaged surface density profile results of all self-gravitating disc models in the Kepler-16, -34 and -35 systems from Paper I. One can see the central cavity in all disc models, and the density spikes associated with the additional eccentric features in the most massive 10 and 20MMSN models. These profiles are calculated once the discs have reached pseudo-steady state – 5000 Pb in the Kepler-16 system, and 6000 Pb in the Kepler-34 and -35 systems. Figure 2. Open in new tabDownload slide Azimuthally averaged surface density profile results of all self-gravitating disc models in the Kepler-16, -34 and -35 systems from Paper I. One can see the central cavity in all disc models, and the density spikes associated with the additional eccentric features in the most massive 10 and 20MMSN models. These profiles are calculated once the discs have reached pseudo-steady state – 5000 Pb in the Kepler-16 system, and 6000 Pb in the Kepler-34 and -35 systems. In our first set of simulations we launch protoplanetary cores, on initially circular orbits, in the outer regions of the evolved discs from Paper I and allow them to interact with the discs. The initial mass of the core in each system is chosen so that qp,0 = mp,0/M* = 6 × 10−5. If M* = 1 M⊙ this is equivalent to a 20 M⊕ core. We release the cores into the outer region of the disc, where the surface density profile is unperturbed by the binary, and the disc eccentricity is negligible. Referring to the profiles in Fig. 2 this lies at 2 au in the low-mass Kepler-16 and -35 models, and around 2.5 au in the low-mass Kepler-34 system. The situation in the high-mass systems is a little more complicated due to the additional eccentric features in the outer disc. Starting the planets at an initial starting position beyond 4 au – exterior to any strong eccentric features – means the time needed to migrate into the inner disc is too long. However, we speculated in Paper I that the migrating planets could interact with these additional features to produce interesting behaviour, therefore we did not want to place the planets too close to the binary. In the high-mass discs, we used an approach which placed the planet beyond the first additional feature, but not too far out in the outer disc. Table 2 summarizes the starting semimajor axes of the cores in all our models. FARGO-ADSG modifies the initial planet semimajor axes with an initial self-gravity ‘boost’, which increases the starting position by ≈10 per cent in the most massive 20MMSN models. Table 2. Starting semimajor axes of protoplanetary cores in each of the Kepler-16, -34 and -35 system models. . ap,0 (au) . . Kepler-16 . Kepler-34 . Kepler-35 . 1MMSN 2.0 2.5 2.5 2MMSN 2.0 2.5 2.5 5MMSN 2.0 2.5 2.5 10MMSN 3.0 2.5 2.0 20MMSN 3.0 2.5 2.5 . ap,0 (au) . . Kepler-16 . Kepler-34 . Kepler-35 . 1MMSN 2.0 2.5 2.5 2MMSN 2.0 2.5 2.5 5MMSN 2.0 2.5 2.5 10MMSN 3.0 2.5 2.0 20MMSN 3.0 2.5 2.5 Open in new tab Table 2. Starting semimajor axes of protoplanetary cores in each of the Kepler-16, -34 and -35 system models. . ap,0 (au) . . Kepler-16 . Kepler-34 . Kepler-35 . 1MMSN 2.0 2.5 2.5 2MMSN 2.0 2.5 2.5 5MMSN 2.0 2.5 2.5 10MMSN 3.0 2.5 2.0 20MMSN 3.0 2.5 2.5 . ap,0 (au) . . Kepler-16 . Kepler-34 . Kepler-35 . 1MMSN 2.0 2.5 2.5 2MMSN 2.0 2.5 2.5 5MMSN 2.0 2.5 2.5 10MMSN 3.0 2.5 2.0 20MMSN 3.0 2.5 2.5 Open in new tab 2.4 Boundary conditions In Paper I, we carried out a fairly exhaustive investigation into the impact of inner boundary conditions on the structure of circumbinary discs around close binary systems. This was motivated by similar discussions in Marzari et al. (2009), PN13, KH14 and Lines et al. (2015), which found discrepancies between different outflow choices, or between the same boundary condition in different systems. For a full description of the investigation and its results we direct the reader to Sections 2.4 and 3 in Paper I. Again, we found no boundary choice which was consistent between the three systems presented here. The choice of inner boundary condition is a way to simulate how much mass flows from the inner disc, out of the boundary and on to the central binary. The Open, Closed or Viscous boundary conditions tested previously all guess at how mass flows through the eccentric cavity, and on to the stars. Our findings from these initial investigations prompted us to develop a way to treat the inner boundary which resulted in a more physically realistic treatment of the material accreting on to the binary. This required us to shrink the size of the inner edge of the disc domain, so the binary is partially embedded in the computational domain rather than sitting entirely interior to the inner boundary. Decreasing the radius of the inner edge, Rin, increases the computational runtime of the simulations, so a compromise is struck between accuracy and speed. We found this balance came when at least 70 per cent of the Roche lobe area of the least massive star was contained in the computational grid at all times. The various values for Rin in the three binary systems simulated here are given in Table 3. The outer boundary of the disc is treated the same in each simulation. A value of Rout = 5 au is used, with an Open outflow boundary condition. A brief description of the Open and Viscous boundary conditions is given below: Open – material is allowed to freely leave the disc i.e. outflow. No inflow is allowed. A zero-gradient condition is set in both υR and Σ. Viscous – this is a limiting condition to stop the inner disc from emptying of gas too quickly. Material in the innermost cells is given a radial velocity, υR = βυR(Rin), where υR(Rin) = −3ν/2Rin, is the viscous drift velocity and β is a free factor (Pierens & Nelson 2008a). We follow previous works which use this condition and set β = 5. Table 3. Inner and outer boundary conditions. . Kepler-16 . Kepler-34 . Kepler-35 . Rin (au) 0.090 0.040 0.056 Rin BC Viscous Rout (au) 5.0 Rout BC Open . Kepler-16 . Kepler-34 . Kepler-35 . Rin (au) 0.090 0.040 0.056 Rin BC Viscous Rout (au) 5.0 Rout BC Open Open in new tab Table 3. Inner and outer boundary conditions. . Kepler-16 . Kepler-34 . Kepler-35 . Rin (au) 0.090 0.040 0.056 Rin BC Viscous Rout (au) 5.0 Rout BC Open . Kepler-16 . Kepler-34 . Kepler-35 . Rin (au) 0.090 0.040 0.056 Rin BC Viscous Rout (au) 5.0 Rout BC Open Open in new tab As we describe in Paper I, the Viscous outflow condition at the inner disc radius tries to model the accretion flow on to the central star(s). It acts as compromise between an unphysical reflecting boundary, and an Open boundary, which empties the inner disc of material at too fast a rate. Our models do not (fully) resolve the circumstellar discs which would form in the Roche lobes of the primary and secondary stars however if we assume that these discs evolve and accrete on to the central objects on the viscous time-scale, the Viscous boundary condition ‘feeds’ these discs at a self-consistent rate. How the azimuthal velocity is treated at the inner radial edge of the disc is also modified. Usually in hydrodynamical codes the viscous stress is maintained by setting υϕ to the sub-Keplerian orbital velocity at the locations Rin and Rout. At the inner boundary the potential created by the binary is extremely non-Keplerian. We therefore set a zero-gradient condition for the azimuthal velocity at this location. 3 MIGRATION OF PROTOPLANETARY CORES In this section, we present the results of simulations examining the migration of protoplanetary cores in evolved self-gravitating discs around the Kepler-16, -34 and -35 binary systems. We insert a non-accreting core, with mass ratio qp,0 = 6 × 10−5, into each of the 1, 2, 5, 10 and 20MMSN discs from Paper I, once the disc has reached a pseudo-steady state. In the Kepler-16 system this is at 5000 Pb, whilst the discs reach this state after 6000 Pb in the Kepler-34 and -35 systems. At this point in the simulation the inner eccentric cavity has a stable precession frequency, and in the high-mass discs, models which were shown to exhibit additional eccentric features at any point, have done so. The core mass used in these simulations lies in the regime where Type-I migration is rapid, but is not massive enough to open a gap in the disc i.e. the gap-opening criteria of Crida, Morbidelli & Masset (2006) is not met. 3.1 Kepler-16 Our results for the orbital evolution of protoplanetary cores in self-gravitating discs around Kepler-16 are summarized in Fig. 4. The upper panel shows the evolution of the cores’ semimajor axes, until a pseudo-steady orbit is reached. In addition to the 1–20MMSN models being shown on this plot, several other quantities are plotted. The red dotted line is the semi-empirical critical semimajor axis for stable orbits around Kepler-16 (Holman & Wiegert 1999); the green dashed line is the best-fitting observed semimajor axis for Kepler-16b from Doyle et al. (2011); the blue dotted line (with label acore) shows the final values for the non-self-gravitating results with comparable disc and core properties from PN13; and the grey dashed lines show the locations of the 5:1–9:1 mean-motion resonances with the binary – locations which have been shown to lead to eccentricity growth, leading to ejections or scattering with the binary (Kostov et al. 2014; Kley & Haghighipour 2015; Kostov et al. 2016). The middle panel shows the evolution of the core eccentricity results for the low-mass (1–5MMSN) disc models, with the high-mass results plotted in the bottom panel for clarity. We can see that like the disc evolution models in Paper I, the evolution of the protoplanetary cores in the low- and high-mass discs can be separated into two distinct regimes of behaviour. In the low-mass discs the cores migrate inwards, albeit with increased rates in the more massive models (the migration rate scales moderately superlinearly with the surface density at the planets’ location; Baruteau & Masset 2008b), from their initial starting position until they reach 1.2 au. This location corresponds well with the surface density peak in the material bounding the tidally truncated inner cavity (Fig. 2), a result fully expected from previous work. As can be seen in Fig. 4 the low-mass results agree extremely well with those from PN13, but not with the observed state of Kepler-16b. This result is slightly unexpected as the disc cavity size seen in our models is somewhat smaller than those in PN13 due to our more realistic treatment of the inner disc boundary. To explain this we must examine what dictates the halting position of these protoplanetary cores in circumbinary discs. From prior work (Pierens & Nelson 2007, 2013; Kley & Haghighipour 2014), we know the stopping behaviour of planets across a range of planetary masses. In the Type-I regime, Earth-like planets are stopped by the growth of a strong positive co-rotation torque which counteracts the influence of the negative Lindblad torque (Masset 2006; Pierens & Nelson 2007). These two torques balance each other when the surface density gradient is sufficiently positive. For more massive Saturn-like planets, a different stopping mechanism operates. If the planetary eccentricity is large enough, a torque reversal can be induced – at apoapse the planet orbits amongst material in the outer disc that is locally travelling faster than itself. When this material overtakes the planet it is focused by the planet's gravity, leading to a positive torque. The reverse of this occurs at periapse, leading to a negative torque from the inner disc (Pierens & Nelson 2008a,b, 2013). When at a cavity edge, the inner torque is naturally smaller in magnitude than the outer torque, leading to a net positive torque arising from this effect. As can be seen in the first panels of Figs 2 and 3, and the middle panel of Fig. 4, the migration of the protoplanets starts to slow when the planetary eccentricity reaches a significant level, ep ≈ 0.13. This coincides when the local disc cell eccentricity and planetary eccentricity are comparable. These findings lead us to the same conclusion as PN13, that for the protoplanetary core mass used here, it is the torque reversal induced by significant planetary eccentricity which halts migration. Comparing our results with those in PN13, the fact that the planet's stopping location is essentially the same in that study and this one, in spite of the different size of the cavity, arises because of differences in the planetary eccentricity and the structure of the cavity (eccentricity and surface density profile). Fig. 5 shows the orbit of the core at the exterior edge of the cavity, where the protoplanet's eccentricity is high enough to induce a torque reversal. Figure 3. Open in new tabDownload slide Azimuthally averaged cell eccentricity profile results of all self-gravitating disc models in the Kepler-16, -34 and -35 systems from Paper I. The eccentric inner cavity can be seen in all models, with the additional eccentricity bumps associated with the eccentric features in the high-mass discs. Again, these profiles are taken once the models have reached pseudo-steady state. Figure 3. Open in new tabDownload slide Azimuthally averaged cell eccentricity profile results of all self-gravitating disc models in the Kepler-16, -34 and -35 systems from Paper I. The eccentric inner cavity can be seen in all models, with the additional eccentricity bumps associated with the eccentric features in the high-mass discs. Again, these profiles are taken once the models have reached pseudo-steady state. Figure 4. Open in new tabDownload slide The top panel shows the evolution of qp,0 = 6 × 10−5 protoplanetary cores’ semimajor axes in evolved self-gravitating discs around the Kepler-16 system. The grey dashed lines in this plot show the positions of n:1 mean-motion resonances with the binary, which have been shown to be unstable to planetary orbits (Nelson 2003; Kostov et al. 2013; Kley & Haghighipour 2014, 2015). The red dotted line shows the location of acrit from Holman & Wiegert (1999). The middle and bottom panels show these cores’ eccentricity evolution in the low- and high-mass disc models, respectively. The green dotted lines in these plots show the values of ap and ep of the observed planet from Doyle et al. (2011), and the blue dotted lines show the final values of simulation with the same disc parameters and comparable core mass from PN13. Figure 4. Open in new tabDownload slide The top panel shows the evolution of qp,0 = 6 × 10−5 protoplanetary cores’ semimajor axes in evolved self-gravitating discs around the Kepler-16 system. The grey dashed lines in this plot show the positions of n:1 mean-motion resonances with the binary, which have been shown to be unstable to planetary orbits (Nelson 2003; Kostov et al. 2013; Kley & Haghighipour 2014, 2015). The red dotted line shows the location of acrit from Holman & Wiegert (1999). The middle and bottom panels show these cores’ eccentricity evolution in the low- and high-mass disc models, respectively. The green dotted lines in these plots show the values of ap and ep of the observed planet from Doyle et al. (2011), and the blue dotted lines show the final values of simulation with the same disc parameters and comparable core mass from PN13. Figure 5. Open in new tabDownload slide A snapshot of the 2D surface density profile in the 1MMSN model once the planet has reached a pseudo-steady-state orbit at the edge of the eccentric cavity. The instantaneous orbit of the planet is shown by the black dashed ellipse, the red dashed circle shows the location of the critical stability limit. The disc and planet system look nearly identical in the 2 and 5MMSN models around the same system. The planet's inwards migration has been halted by a strong positive co-rotation torque balancing the negative Lindblad torque, once the eccentricity of the planet attains a value ≈ed(ap). Figure 5. Open in new tabDownload slide A snapshot of the 2D surface density profile in the 1MMSN model once the planet has reached a pseudo-steady-state orbit at the edge of the eccentric cavity. The instantaneous orbit of the planet is shown by the black dashed ellipse, the red dashed circle shows the location of the critical stability limit. The disc and planet system look nearly identical in the 2 and 5MMSN models around the same system. The planet's inwards migration has been halted by a strong positive co-rotation torque balancing the negative Lindblad torque, once the eccentricity of the planet attains a value ≈ed(ap). In the high-mass disc regime we observe planetary evolution behaviour not seen in the low-mass discs, or previous work on this topic. Whilst the cores still migrate inwards, in the 10MMSN model the core briefly halts at 2.2 au. During the period when it is trapped at this location its eccentricity steadily grows from 0.05 to a maximum of 0.4. When the core is then released its eccentricity is quickly damped and it migrates into the inner disc, halting at a location in good agreement with the low-mass results, ap = 1.1 au. Unlike the low-mass models the core's eccentricity is slowly damped by the large amount of gas in its vicinity, down to a value in good agreement with the observed value of Kepler-16b. Whilst the core in the 20MMSN model does not show signs of this trapping immediately, its inwards migration is halted at 1.5 au, a location significantly exterior to the prior results. Examining Figs 6 and 7, we can start to explain this behaviour. The discovery of the additional eccentric features in high-mass self-gravitating circumbinary discs prompted us to theorize that they could act as planet traps. In the Type-I regime, the positive surface density gradient creates a strong co-rotation torque which could counteract the Lindblad torque or, for more massive planets, the excited eccentricity in these regions could excite the eccentricity of the body sufficiently to induce a torque reversal. Whilst this process requires the planet to have a non-negligible eccentricity, it also requires there to be a surface density gradient across the extremes of the orbit. At apocentre it should find itself in an area of high surface density, and a low surface density at pericentre. In the discs which we obtain in these models, this can be achieved by the planet and disc eccentricities not being exactly equal, or a misalignment between the respective line of nodes. In this case, the planet is on a less eccentric orbit than the surrounding disc material. Figure 6. Open in new tabDownload slide Evolution of the 2D and azimuthally averaged 1D surface density profiles of the Kepler-16 10MMSN model, over the course of the protoplanet's migration from its initial starting point, to its final location. The instantaneous orbit of the planet is shown by the black dashed ellipse in the 2D plots; the red dashed circle shows the location of the critical stability limit. The black line in the 1D profiles is the core's actual distance from the binary CoM, the orange dashed line is ap and the inner and outer dashed red lines show the location of rper = ap(1 − ep) and rapo = ap(1 + ep), respectively. Figure 6. Open in new tabDownload slide Evolution of the 2D and azimuthally averaged 1D surface density profiles of the Kepler-16 10MMSN model, over the course of the protoplanet's migration from its initial starting point, to its final location. The instantaneous orbit of the planet is shown by the black dashed ellipse in the 2D plots; the red dashed circle shows the location of the critical stability limit. The black line in the 1D profiles is the core's actual distance from the binary CoM, the orange dashed line is ap and the inner and outer dashed red lines show the location of rper = ap(1 − ep) and rapo = ap(1 + ep), respectively. Figure 7. Open in new tabDownload slide A 2D surface density profile of the evolved 20MMSN disc model around the Kepler-16 binary. In this snapshot a protoplanetary core has been interacting with the disc for 1000 Pb. The instantaneous orbit of the core is shown by the black dashed ellipse. The core can be seen to be orbiting between the inner cavity, and the first eccentric feature. The planet's orbit crosses this feature, and because it has a period of precession ranging in the hundreds of binary orbits (see Paper I), the planet will cross and interact with this feature repeatedly. The eccentric feature is strong enough to trap the planet at this exterior position. Figure 7. Open in new tabDownload slide A 2D surface density profile of the evolved 20MMSN disc model around the Kepler-16 binary. In this snapshot a protoplanetary core has been interacting with the disc for 1000 Pb. The instantaneous orbit of the core is shown by the black dashed ellipse. The core can be seen to be orbiting between the inner cavity, and the first eccentric feature. The planet's orbit crosses this feature, and because it has a period of precession ranging in the hundreds of binary orbits (see Paper I), the planet will cross and interact with this feature repeatedly. The eccentric feature is strong enough to trap the planet at this exterior position. Despite this prediction, we see two different end results in the 10 and 20MMSN models. In the least massive of these cases the core migrates inwards until it reaches the first eccentric feature. At this location it halts. Repeated interaction with this highly eccentric feature lead to the planet's own eccentricity being excited. This can be seen in the second panel of Fig. 6. However, the planet's eccentricity becomes so high (ep = 0.4), the planet's pericentre position decreases until it interacts with the material bounding the inner eccentric cavity. A similar process to the initial trapping then occurs, however the planet's orbit is circularized by the far less extended inner feature. The core's semimajor axis shrinks until the orbit is moderately eccentric, which matches that of the observed Kepler-16b relatively well, at the location of the inner cavity. The core in the 20MMSN model can get trapped at the first outer eccentric cavity because the eccentric feature is more tightly localized due to the disc's stronger self-gravity. Therefore, the planet's orbit does not take it into close proximity of the strong inner eccentric feature, and it remains trapped between the inner and first outer eccentric features. 3.2 Kepler-34 As a result of the variety in the evolved disc structures in the Kepler-34 models, we see a large range of results for the migration of protoplanetary cores in these evolved discs. In Paper I, we found that as we increase the disc mass from 1 to 20MMSN, the size of the initially very eccentric, extended cavity in the least massive disc gradually decreases – self-gravity acts to compact the scale of the system. As can be seen from Fig. 8 a similar pattern can be seen in the final stopping positions of the migrating cores, where the halting of migration occurs because ep increases and induces a torque reversal. Figure 8. Open in new tabDownload slide The top panel shows the evolution of qp,0 = 6 × 10−5 protoplanetary cores’ semimajor axes in evolved self-gravitating discs in the Kepler-34 system. The middle and bottom panels show these cores’ eccentricity evolution in the low- and high-mass disc models, respectively. Figure 8. Open in new tabDownload slide The top panel shows the evolution of qp,0 = 6 × 10−5 protoplanetary cores’ semimajor axes in evolved self-gravitating discs in the Kepler-34 system. The middle and bottom panels show these cores’ eccentricity evolution in the low- and high-mass disc models, respectively. In the low-mass discs there is a clear trend for ap to decrease from ap ≈ 1.95 to 1.5 au, and for ep to increase from ep = 0.275 to 0.3 as the disc mass rises from 1 to 5MMSN. Looking at the disc eccentricity distributions in the middle panel of Fig. 3, for a given radius one obtains a smaller value for the average disc cell eccentricity for larger disc masses – the core therefore has to migrate further through the disc so that ep ≈ ed. The large ep seen for the cores in these discs means the strength of the co-rotation torque would be greatly diminished (Fendyke & Nelson 2014), therefore making torque reversal the dominant mechanism for halting migration. The low-mass discs, especially the 1MMSN model, do not match the results obtained in PN13, and show poor agreement with the observed Kepler-34b. The discs in these models tend to have large, highly eccentric cavities compared to the equivalent models in PN13. Our more realistic treatment of the inner disc boundary, allowing for a more accurate capturing of angular momentum flux through the disc due to the binary, is the likely explanation for this. Whilst the low-mass disc results do not agree well with past results or the observed state of the Kepler-34 planetary system, the 10MMSN and 20MMSN models agree relatively well with the planetary orbital elements quoted by Welsh et al. (2012). The final stopping positions in the 10MMSN and 20MMSN systems, 1.2 and 1.0 au, respectively, bracket the observed value of ap ≈ 1.1 au due to the compacting of the system as disc-mass and self-gravity increase. The trend for ep to increase as the disc-mass increases is reversed in the high-mass regime, possibly due to the disc-mass in the vicinity of the planet providing significant damping. We see no evidence of the core being trapped in the outer disc in either model. Examining the last two panels in Fig. 9, we can see that although additional eccentric features are present in the outer disc, in both the 10 and 20MMSN models, they are far less well-defined than those in Kepler-16. These washed-out features are not strong enough to halt the inwards migration of the cores in this system. Figure 9. Open in new tabDownload slide 2D surface density profiles showing the structure of the binary-disc-core system once migration has halted for the 1–20MMSN disc models around the Kepler-34 binary. The planetary orbit is shown by the black dashed line, with the red dashed line showing the critical stability limit. For each of these models it can be seen that the planet migrates into the inner disc, where it then interacts with the strong eccentric feature. The stopping position can be seen to decrease as the disc mass increases. Clear evidence of pericentre-alignment between the core and the eccentric cavity can be seen. Figure 9. Open in new tabDownload slide 2D surface density profiles showing the structure of the binary-disc-core system once migration has halted for the 1–20MMSN disc models around the Kepler-34 binary. The planetary orbit is shown by the black dashed line, with the red dashed line showing the critical stability limit. For each of these models it can be seen that the planet migrates into the inner disc, where it then interacts with the strong eccentric feature. The stopping position can be seen to decrease as the disc mass increases. Clear evidence of pericentre-alignment between the core and the eccentric cavity can be seen. A common feature that can be seen for all the disc models in this system can be observed in Fig. 9. One can clearly see that the orbits of the cores in each system are aligned with the precessing eccentric inner cavity. An examination of the evolution of the planet's longitude of pericentre, alongside that of the mean disc longitude of pericentre shows this as well (Fig. 10). As the planet migrates into the inner disc, the phase and period of precession both evolve into lockstep with that of the inner disc cavity (which the local calculation of ωd traces). The planet and eccentric feature precess with each other, in a pericentre-aligned fashion, a behaviour previously seen for full mass planets in non-self-gravitating discs in the Kepler-34 system (Kley & Haghighipour 2015). This is not true in the most massive disc model presented here, 20MMSN. In this case, the precession of the planet and disc are half a precession period out of phase such that their eccentric orbits are anti-aligned. In Paper I, we give an explanation for the global and local calculations of the disc eccentricity and longitude of pericentre. Whilst the global calculation takes into account all the material in the disc between Rin and Rout, the local calculation only takes into account material up to and just beyond the position of the surface density peak associated with the inner cavity. This procedure ignores the effect of exterior eccentric, precessing material in the outer disc. Figure 10. Open in new tabDownload slide Evolution of longitude of pericentre of the disc-planet system in the 1, 10 and 20MMSN disc mass models around Kepler-34. Both the global and local calculations of the disc longitude of pericentre are included (see the text for a description of the differences between these two calculations. In the 1–10MMSN discs, the planet evolves into a state where the phase and period of circulation of its precession match that of the inner eccentric cavity of the disc. The two match when the planet has halted its migration at the inner cavity. Figure 10. Open in new tabDownload slide Evolution of longitude of pericentre of the disc-planet system in the 1, 10 and 20MMSN disc mass models around Kepler-34. Both the global and local calculations of the disc longitude of pericentre are included (see the text for a description of the differences between these two calculations. In the 1–10MMSN discs, the planet evolves into a state where the phase and period of circulation of its precession match that of the inner eccentric cavity of the disc. The two match when the planet has halted its migration at the inner cavity. We note that a second separate simulation of a protoplanetary core released at 3 au in the Kepler-34 20MMSN disc model was undertaken to examine migration from a larger radius. This location corresponds to a radius between the second and third additional eccentric features in the disc. Whilst these features are relatively weak, they still alter the surface density profile of the disc. These regions of positive surface density gradient are sufficient to hamper any inwards migration of the planet, but insufficient to excite sufficient eccentricity for it to escape. The forces acting on the planet at this outer position which normally result in inwards migrations are overcome by the small perturbations in surface density – whilst at the starting position of the first 20MMSN core the rate of inwards migration is greater. These weak features could play an important role in the early stages of planet formation, trapping large numbers of planetesimals, boulders or pebbles, but with low eccentricity, providing a reservoir for protoplanetary core creation. For clarity we have not included this simulation in the plots for this section, but we will discuss its further evolution in subsequent sections. 3.3 Kepler-35 The results from protoplanetary core migration in the evolved Kepler-35 disc models look very similar to those from the Kepler-16 models, with minor changes caused by the differences in evolved disc structure. This result was expected due to the similarity in evolution and final structure results from Paper I for these low-eccentricity binaries. To recap, in the low-mass discs, the cores migrate inwards through the disc until they reach a location where their eccentricity is excited enough (ep ≈ ed(ap)) to induce a torque reversal. This location corresponds to the edge of the eccentric cavity. The location of this edge, which is easily identifiable as the peak in the surface density, lies at a smaller radius in the Kepler-35 system than in Kepler-16, due to its lower binary eccentricity. The final semimajor axis for the cores in the 1–5MMSN models is 0.9 au – slightly smaller than the previous PN13 result – and ep oscillates around 0.11 for all three models (see the upper and middle panels of Fig. 11). The mean value of ep results from a balance between the highly eccentric disc pumping up the eccentricity and the surrounding material damping the eccentricity. None of the final planetary orbital elements obtained in the low-mass regime are in good agreement with those quoted in Welsh et al. (2012) for Kepler-35b. Figure 11. Open in new tabDownload slide The top panel shows the evolution of qp,0 = 6 × 10−5 protoplanetary cores’ semimajor axes in evolved self-gravitating discs in the Kepler-35 system. The middle and bottom panels show these cores’ eccentricity evolution in the low- and high-mass disc models, respectively. Figure 11. Open in new tabDownload slide The top panel shows the evolution of qp,0 = 6 × 10−5 protoplanetary cores’ semimajor axes in evolved self-gravitating discs in the Kepler-35 system. The middle and bottom panels show these cores’ eccentricity evolution in the low- and high-mass disc models, respectively. The results from the high-mass models also show the same evolutionary history, with slightly different final values for ap and ep, as the Kepler-16 high-mass models. The 10MMSN model shows evidence of trapping by the m = 1 eccentric mode at 1.6–1.9 au, where its eccentricity gets rapidly excited to 0.5. This highly eccentric orbit then brings the pericentre close enough to the inner cavity to allow the planet to be captured by the large amount of material skirting the boundary. This material damps the orbit of the planet, decreasing the semimajor axis (ap = 0.9 au) and eccentricity (ep = 0.03), to a near-circular orbit (see upper and lower panels of Fig. 11, and Fig. 12). Figure 12. Open in new tabDownload slide Evolution of pericentre distance of the protoplanetary cores in the 10MMSN and 20MMSN disc models. Of particular note is the decrease between 3000 and 5000 Pb. During this time the core has been trapped at the location of the first additional eccentric feature. Whilst its eccentricity is being excited, the pericentre distance is decreasing, to the point where it reaches the location of the inner eccentric cavity. Whilst the pericentre distance remains relatively constant, the eccentricity and semimajor axis are damped by this massive feature. Figure 12. Open in new tabDownload slide Evolution of pericentre distance of the protoplanetary cores in the 10MMSN and 20MMSN disc models. Of particular note is the decrease between 3000 and 5000 Pb. During this time the core has been trapped at the location of the first additional eccentric feature. Whilst its eccentricity is being excited, the pericentre distance is decreasing, to the point where it reaches the location of the inner eccentric cavity. Whilst the pericentre distance remains relatively constant, the eccentricity and semimajor axis are damped by this massive feature. The 20MMSN model, matching the evolution of the core in the Kepler-16 20MMSN disc, migrates inwards through the disc – keeping a low eccentricity, ep ≈ 0.08 – until it is trapped at the first extra eccentric feature, with a final semimajor axis, ap = 1.4 au. The eccentricity damping provided by the disc is sufficient that it remains at this location. Both the 10 and 20MMSN discs produce cores whose final eccentricity is in good agreement with that of Kepler-35b, but the simulated semimajor axes are too large. 4 MIGRATION OF ACCRETING CORES The simulations that have been presented so far in this paper all adopted a fixed mass for the planetary cores, corresponding to a mass ratio between the planet and central binary of qp = 6 × 10−5. The actual mass ratios for the observed systems are all larger than this by various factors (see Table 1), and so we now consider what happens to the orbital elements if the planets accrete gas and achieve their observed masses while migrating. The results shown above indicate that the planets considered so far normally halt their migration at a location that is too far from the binary to provide good agreement with the observations, so we examine whether or not the stopping orbital radii decrease as we increase the planet masses to their observed values. Only the 10MMSN and 20MMSN cases for the Kepler-34 system produce final stopping radii that agree well with the observations, and this occurs because of the dramatic shrinking of the cavity size in this case for discs where self-gravity is important. For the simulations presented in this section we undertook accretion scenarios for the evolved binary-disc-planet Kepler-16 and Kepler-34 systems. The initial core mass ratio used throughout this work means that the protoplanetary core in the Kepler-35 models is within ≈20 per cent of the observed planet mass, so we did not simulate gas accretion in this case. The accretion routine of Kley (1999) was used to grow the mass of the protoplanetary cores (qp, 0 = 6 × 10−5) to that of the observed planet mass in the specific system (qp = 3.54 × 10−4 and 1.01 × 10−4 in the Kepler-16 and -34 systems, respectively). This prescription removes a portion of the gas from the Hill sphere and adds its mass to that of the planet. The accretion time-scale, i.e. the time in which the Hill sphere is emptied of gas, is determined as a fraction of the dynamic time-scale of the planet, tacc = f tdyn. The variable constant f is tuned for the Kepler-16 simulations so that the planet reaches its final mass over 5000 Pb. We use this approach to inhibit the growth of the planet. A constant value is used (f = 0.01−1) for all the disc-mass models in the Kepler-34 system, as the final planet mass is relatively low. The issue worth noting with this set-up, in relation to a realistic comparison with the masses of the observed circumbinary planets, is that when accretion is turned on in these simulations the planet finds itself at a location with a wealth of material. Even conservative estimates for the accretion time-scale lead to rapid mass growth. If planets in circumbinary discs only reach a gas-accretion phase when they are already at the cavity edge, it would be logical to assume that the planet could quickly grow to Jovian mass unless gas accretion is very slow indeed, or occurs at the end of the disc lifetime. Our simulations apply to the former possibility, but it is worth noting that circumbinary systems are self-selecting because too much gas accretion leads to the formation of a Jovian-mass planet, and these tend to be much more unstable due to dynamical interaction with the central binary (Nelson 2003). Even if circumbinary planets grow to be of Jovian-mass close to the cavity edge, we are unlikely to see them as they have a significant probability of being ejected from the system. The accretion scenarios that we consider here are run from the point in the simulations from the last section when the planet has reached a pseudo-steady orbit. 4.1 Kepler-16 Kepler-16b, with mp ≃ 0.3 MJup is the most massive of the three circumbinary planets that we consider in this work (Doyle et al. 2011). Using the gap-opening criteria of Crida et al. (2006) which states that for a given set of disc parameters, a planet of mass ratio, q will open a gap if \begin{equation} 1.1\left(\frac{q}{h^3}\right)^{-1/3} + \frac{50\alpha h^2}{q} \le 1, \end{equation} (5) one can see that for the viscous stress parameter and disc aspect ratio used in these simulations, the core will significantly alter the surface density profile of the disc when it approaches its final mass – as can be seen in Fig. 13. In the low-mass discs the core slowly migrates from its initial stopping position at ap = 1.1 au further inwards to ap = 0.75 au, between the 6:1 and 7:1 MMRs with the binary (top panel of Fig. 14). As the planets migrate into the cavity evacuated by the binary, they carve out this cavity further – opening one side of a gap. This process destroys the eccentric cavity, as the planet's mass dominates – resulting in a decrease in eccentricity ep = 0.12 → 0.03. One anomaly in these results is the rapid outward migration of the core in the 5MMSN model at 4000 Pb (Fig. 14). This is accompanied by a sharp decrease in ep. Examining the evolution of the planet at this epoch, it can be seen that as the planet's mass grows it appears to interact with the 8:1 MMR with the binary. These n:1 MMR locations have been shown to be unstable to planetary orbits because they excite the eccentricity (Nelson et al. 2000), and in this case the planet is scattered out. It is not ejected and is able to migrate back into the inner disc, avoiding further scattering events. Figure 13. Open in new tabDownload slide 1D snapshots of the surface density profile in the 1MMSN (top panel) and 10MMSN (bottom panel) models in Kepler-16, once accretion on to the protoplanetary core is allowed. In the low-mass disc the core has grown from its initial mass to the final, observed mass of Kepler-16b between the first and second profiles. It can be seen that whilst the planet has migrated further into the inner disc, it has also carved out a pseudo-gap in the inner edge of the cavity. In the massive disc, whilst the planet does not seem to open a gap, it does significantly alter the surface density profile over the whole radial extent. Figure 13. Open in new tabDownload slide 1D snapshots of the surface density profile in the 1MMSN (top panel) and 10MMSN (bottom panel) models in Kepler-16, once accretion on to the protoplanetary core is allowed. In the low-mass disc the core has grown from its initial mass to the final, observed mass of Kepler-16b between the first and second profiles. It can be seen that whilst the planet has migrated further into the inner disc, it has also carved out a pseudo-gap in the inner edge of the cavity. In the massive disc, whilst the planet does not seem to open a gap, it does significantly alter the surface density profile over the whole radial extent. Figure 14. Open in new tabDownload slide Evolution of accreting protoplanetary cores’ semimajor axes in evolved self-gravitating discs in the Kepler-16 system. The middle and bottom panels shows these cores’ eccentricity evolution in the low- and high-mass disc models, respectively. Figure 14. Open in new tabDownload slide Evolution of accreting protoplanetary cores’ semimajor axes in evolved self-gravitating discs in the Kepler-16 system. The middle and bottom panels shows these cores’ eccentricity evolution in the low- and high-mass disc models, respectively. The planet in the 10MMSN disc alters the surface density similarly to the planets in the low-mass discs. It can be seen in the bottom panel of Fig. 13 that it does not open such a deep gap at the cavity edge. The presence of the planet leads to the destruction of the additional eccentric features in the outer disc. Any planet forming and evolving in the outer disc in a multiplanet formation and migration scenario (see Kley & Haghighipour 2015) would have a very different migration pathway to the first planet. During the accretion phase of the simulation, slow inwards migration occurs to ap = 1.0 au, whilst the eccentricity of the core's orbit falls. The next 10 000 Pb is spent at this distance, after which it migrates further into the inner disc, where it reaches ap = 0.75 au. For the remainder of the simulation lifetime it has an eccentricity ep = 0.04, in relatively good agreement with Kepler-16b. However, between 2.5 × 104 and 3 × 104 Pb the core seems to undergo a similar scattering event as that seen in the 5MMSN disc – oscillating around the 7:1 MMR and consequently scattering out. This scattering and subsequent inwards migration seems to happen repeatedly over the course of the simulation. The mass of the disc is sufficient to maintain a significant eccentricity of its own, and excite that of the planet. The evolution of the accreting core in the most massive 20MMSN disc is even more disruptive. Initially, when its mass starts to grow, in the first few 1000 Pbs of the simulation, it escapes the outer planet trap and migrates into the inner disc (ap = 0.9 au). During this phase it maintains a significant eccentricity (ep ≈ 0.1) – because of this its orbit enters the n:1 MMR region. It spends the remainder of the disc undergoing repeated scattering and migration events, where its eccentricity dramatically rises to ≈0.3 and then circularizes. If this continues, the core could in principal enter the critical stability limit during one of these events and be ejected from the system – although we have not yet seen this happen. Fig. 15 we can see that the structure of the circumbinary disc has been significantly altered by the growth of the planet to its observed mass. The opening of the gap, as well as strong spiral wakes launched at the Lindblad resonances with the planet act to destroy the eccentric cavity, making it more circular. This is clear in the top panel of Fig. 16; during the first 10 000 Pb of the simulation, when the planet is accreting mass from the disc and migrating slowly inwards, the eccentricity of the disc decreases to 0.01. In the most massive discs, the growing and migrating planet disrupts the eccentric features in the exterior disc as well as the inner cavity – leading again to a decrease in ed (bottom panel of Fig. 16). The erratic changes in the orbit of the planets in the 10 and 20MMSN models also leads to corresponding fluctuations in the disc eccentricity. When the planet is on a wider, more eccentric orbit, the eccentricity of the disc can also grow. Figure 15. Open in new tabDownload slide Surface density map of the Kepler-16 2MMSN disc model once the core has grown to its observed mass and reached its final semimajor axis. The near-circularity of the planetary orbit is clear, whilst the complex interplay between the density waves launched from the Lindblad resonances with the planet, and those launched by the binary in the self-gravitating disc, has destroyed the eccentric inner cavity. A gap has been opened with a circular cavity interior to it. Figure 15. Open in new tabDownload slide Surface density map of the Kepler-16 2MMSN disc model once the core has grown to its observed mass and reached its final semimajor axis. The near-circularity of the planetary orbit is clear, whilst the complex interplay between the density waves launched from the Lindblad resonances with the planet, and those launched by the binary in the self-gravitating disc, has destroyed the eccentric inner cavity. A gap has been opened with a circular cavity interior to it. Figure 16. Open in new tabDownload slide Comparison of planetary and disc (global) eccentricity evolution once gas accretion has started (at t = 0 Pb) in the 2 and 10MMSN models around the Kepler-16 system. The growth of the core to its observed mass in the first 5000 Pb of the scenario results in significant alteration of the disc eccentricity profile, as well as the average disc eccentricity. As the planet grows and it approaches the gap-opening regime, a decrease in disc eccentricity can be seen – in parallel to the decrease in planetary eccentricity. In the 10MMSN model the erratic changes in ap beyond 2.5 × 104 Pbs are also accompanied by a significant growth in ed. Figure 16. Open in new tabDownload slide Comparison of planetary and disc (global) eccentricity evolution once gas accretion has started (at t = 0 Pb) in the 2 and 10MMSN models around the Kepler-16 system. The growth of the core to its observed mass in the first 5000 Pb of the scenario results in significant alteration of the disc eccentricity profile, as well as the average disc eccentricity. As the planet grows and it approaches the gap-opening regime, a decrease in disc eccentricity can be seen – in parallel to the decrease in planetary eccentricity. In the 10MMSN model the erratic changes in ap beyond 2.5 × 104 Pbs are also accompanied by a significant growth in ed. In summary, we find that increasing the planet's mass in the low-mass discs leads to further inwards migration, and final orbital elements that are in rather good agreement with the observed values. Gas accretion in the high-mass discs, however, leads to repeated interactions with the binary that cause the orbits of the planets to change erratically. 4.2 Kepler-34 Increasing the protoplanet's mass from qp,0 = 6 × 10−5 to the observed mass of Kepler-34, qp = 1 × 10−4 – an increase a little over 60 per cent – results in little change of orbital parameters. This is unsurprising as the core is still in the Type-I planet migration regime in our disc models, and according to equation (5), is not capable of sufficiently disturbing the surface density distribution to open a gap. This lack of significant activity was apparent after a relatively short simulation time (≈2000 Pb), where after a period of relaxation the system reaches a pseudo-steady state. In Fig. 17, one can see a slight outward migration of the planets in the 2MMSN and 5MMSN models, associated with a circularization of the orbit. A lack of change in the low-mass models mean that there is still poor agreement with the observed configuration of Kepler-34b. The semimajor axes of the cores in this system are too large (ap = 1.7–2 au), with eccentricities which are too excited (ep = 0.25–0.275) – although good agreement with the observed eccentricity is obtained for the 5MMSN core model where ep is oscillating around 0.18. The 10 and 20MMSN model cores also show little change when accretion is switched on, apart from a slight decrease in eccentricity in the 10MMSN case due to more efficient damping by the disc. Figure 17. Open in new tabDownload slide Evolution of accreting protoplanetary cores’ semimajor axes in evolved self-gravitating discs in the Kepler-34 system. The middle and bottom panels shows these cores’ eccentricity evolution in the low- and high-mass disc models, respectively. The difference between the two cores in 20MMSN models is the initial starting position. It can be seen that a small increase in planet mass is sufficient to alter the surface density in the outer disc, so that the previously trapped core can migrate inwards. Figure 17. Open in new tabDownload slide Evolution of accreting protoplanetary cores’ semimajor axes in evolved self-gravitating discs in the Kepler-34 system. The middle and bottom panels shows these cores’ eccentricity evolution in the low- and high-mass disc models, respectively. The difference between the two cores in 20MMSN models is the initial starting position. It can be seen that a small increase in planet mass is sufficient to alter the surface density in the outer disc, so that the previously trapped core can migrate inwards. The second 20MMSN model run in the Kepler-34 system which is initially released further out in the disc but is trapped close to its starting position – mentioned at the end of Section 3.2 – shows the most dramatic response to accreting mass. The increased core mass is sufficient for it to escape the region of weak eccentric features, created by the self-gravitating disc response to the binary potential, in the outer disc. It quickly migrates through the disc, finally reaching ap = 1.0 au, the same as the first 20MMSN model presented and in good agreement with the observed value of ap, although both models have small values of ep. In this system it is especially hard to produce a planet so close-in with a non-negligible eccentricity that matches the observations. Increasing ep, hence lowering the pericentre distance, further increases the risk of destabilizing encounters with the n:1 MMR region and the chance of a catastrophic ejection event. This, along with post-disc dissipation evolution with the binary, may be the reason why we may yet to observe a very close-in circumbinary planet (like Kepler-16b or -35b) with a significant eccentricity like that of Kepler-34b. 5 DISC DISSIPATION IMPACT ON MIGRATING CORES From the beginning of this investigation we have been using the disc-mass as a proxy for the age of the circumbinary disc. It is logical to assume that when the disc first forms into a stable entity around the central binary it is at its most massive, and over the course of its lifetime loses mass due to a number of different processes – accretion on to the central binary, loss from photo-evaporative and/or magnetized winds from the surface of the disc, etc. Whilst we have simulated the disc structure and evolution at different eras throughout its lifetime, we have not investigated the effects of transitioning from a high-mass environment to that of a low-mass one. The dichotomy of results from Paper I suggest that the additional eccentric features seen in the outer disc will disperse as the disc-mass and the strength of self-gravity decrease. Without a sustaining action, the viscous forces in the disc will dissipate these eccentric features. As the strength of self-gravity diminishes in the disc we would also expect the compactness of the system to relax back to that seen in the least-massive 1MMSN disc – the eccentric cavity will increase in size, especially those seen in the Kepler-34 system. The surface density profile will alter as the disc relaxes and we would expect the planet to migrate outwards with the cavity. In those discs where the planets are halted by the counteracting of the Lindblad torque by the positive co-rotation torque, the planet may be able to stay at this stable stopping location whilst the disc relaxes. We are not investigating the mechanisms and physics which dictate disc mass-loss and dispersal – these are topics of ongoing research – and perhaps deserve their own work in the context of massive self-gravitating discs. Instead, and as a computation time saving exercise, we use a simple exponential decay to dissipate the mass of the disc: \begin{equation} \Sigma ^{n+1}_{ij} = \Sigma ^{n}_{ij} \exp {\left(-\frac{\Delta t}{\tau }\right)}, \end{equation} (6) where Σij is the cell surface density value, Δt = tn + 1 − tn, is the time between successive time levels n and n + 1, and τ is the decay time constant. In each disc this value is chosen so that after 5000 Pb the total disc mass md will have decreased from its initial value down to the equivalent 1MMSN model in that system. For reference, to reach a 1MMSN disc from a 10MMSN or 20MMSN mass disc in the Kepler-16 system, a time constant of τ = 1520 or 1175 is used, respectively. This decay length is sufficiently large that the dynamical time-scales associated with the disc and planet are much smaller, and can therefore respond to any changes in the disc. Once the disc has reached a total disc-mass equivalent to the initial 1MMSN disc mass, the dissipation mechanism is stopped, to allow the disc and planet to reach a pseudo-steady state on time-scales of a few ×10 000 binary orbits. This procedure is started in the disc once – similarly to the previous subsection looking into accretion scenarios – the initial binary-disc-protoplanet systems from Section 3 have reached quasi-steady state. This allows us to track the response of planets, trapped at the cavity edge or by eccentric rings, to the diminishing disc mass and relaxation or dissipation of eccentric features. During this procedure, and in the post-dissipation evolution of the planets, we consider non-accreting cores. 5.1 Kepler-16 The similarity of the results from the 1–5MMSN models seen in Section 3.1 for the Kepler-16 system prompted us to only carry out disc dispersal simulations for the most massive 10MMSN and 20MMSN models. Fig. 18 shows the response of the planets’ semimajor axes (top panel) and eccentricities (bottom panel) to the disc dispersal, which occurs during the first 5000 Pb of these plots. A significant amount of post-dissipation evolution of the cores can be seen, especially in the 20 →1MMSN model. Figure 18. Open in new tabDownload slide Evolution of protoplanetary cores’ semimajor axes in dissipating self-gravitating discs, in the Kepler-16 system (top panel). The bottom panel shows the core eccentricity evolution in the 10MMSN and 20MMSN disc models. Figure 18. Open in new tabDownload slide Evolution of protoplanetary cores’ semimajor axes in dissipating self-gravitating discs, in the Kepler-16 system (top panel). The bottom panel shows the core eccentricity evolution in the 10MMSN and 20MMSN disc models. As previously detailed, the core in the 10MMSN disc is orbiting at the edge of the inner eccentric cavity when migration halts. When dissipation starts to occur, the semimajor axis of the core increases from 1 to 1.2 au – the position of Rmax, or the cavity edge in the least massive 1MMSN model. The eccentricity of the core also increases in this period, increasing from around 0.05 to 0.1. This increase in eccentricity, due to reduced damping by the disc, increases the positive torque contribution from the outer disc, even as it relaxes due to dissipation. The balance between reduced eccentricity damping and the diminishing influence of the positive torque from the outer disc, dictates whether the core migrates inwards or outwards as the disc dissipates. After this initial period of outwards migration, this balance inverts. The semimajor axis decreases – past the initial stopping distance – further into the inner disc as the eccentricity drops – reaching a final orbit with ap = 0.9 au and ep ≈ 0.025. This model gives better final agreement with the observed Kepler-16 system than the low-mass models. When dissipation starts, the core's small ap and non-negligible ep mean its pericentre distance lies around 0.8 au. The core retains this small value during dissipation, and as the eccentricity is damped by the disc. The final semimajor axis corresponds to a location between the 8:1 and 9:1 MMR with the binary – the core's low eccentricity however keeps it clear of interaction with these destabilizing regions. Examining the surface density profile and planetary orbit in Fig. 19, one can see the similarity to the low-mass Kepler-16 discs (Fig. 5). The eccentric features in the outer disc have dissipated and the inner eccentric cavity has relaxed, to a size in good agreement with the 1–5MMSN models. The difference in the shape of the planetary orbit is also clear, a more circular orbit inside, rather than tracing the outside edge of the cavity, is attained. Figure 19. Open in new tabDownload slide Surface density map of the Kepler-16 10MMSN disc model once the disc has undergone mass dissipation, and the core has reached its final orbital position. Whilst the eccentric cavity has relaxed from the very compact initial 10MMSN state to a cavity with a larger extent, the core has managed to retain a close-in, circular orbit. Figure 19. Open in new tabDownload slide Surface density map of the Kepler-16 10MMSN disc model once the disc has undergone mass dissipation, and the core has reached its final orbital position. Whilst the eccentric cavity has relaxed from the very compact initial 10MMSN state to a cavity with a larger extent, the core has managed to retain a close-in, circular orbit. A very different evolution is seen in the 20MMSN model. The core in this model is trapped at the location of the first eccentric feature in the outer disc when the process of dissipation begins. During dissipation the planet migrates outwards from ap = 1.4 to 2.2 au. After the first 2000 Pb of disc dissipation, migration reverses and the core migrates into the inner disc, reaching a final semimajor axis of 1.1 au. Examining the azimuthally averaged surface density profiles in Fig. 20 this evolutionary history can be explained. Comparing the three plots several important things can be extracted; the first being that between the first two panels the surface density profile has relaxed to one resembling a 10MMSN profile. The tightly wound eccentric features in the 20MMSN model dissipate outwards in the disc, but one relatively strong eccentric feature at 2.1 au can still be seen. If dissipation stopped here, we might see a migration scenario much like the 10MMSN model from Section 3.1, where the planet gets trapped, but then subsequently escapes. Dissipation does continue however, and the core is free to immediately migrate into the inner disc, as the eccentric features are destroyed. The last panel shows the planet at its final location of ap = 1.1 au. This value and the surface density profile are very similar to the final state of the 1MMSN model from Section 3.1, and its final ep ≈ 0.075 is in good agreement with observations. To summarize, the core – still trapped by the eccentric feature – migrates outwards as it follows the dissipating perturbation, until the eccentricity of the core diminishes enough so that the positive torque contribution from the outer disc stops. The net negative torque migrates the core inwards towards the central cavity, where the eccentricity increases again, inducing another torque reversal, halting migration in the inner disc. Figure 20. Open in new tabDownload slide 1D surface density profiles are shown here during the evolution of the 20MMSN disc, and its embedded core. All line types and colours carry the same meaning as previous plots of this type. The first panel shows the system at the beginning of disc dissipation, the second is 2000 Pb orbits through this process, and the final subplot shows the system once the disc has relaxed and the planet has reached its final stopping position. The middle panel's surface density distribution looks very similar to those seen from 10MMSN, whilst the last panel looks very similar to a 1MMSN model. The planet initially migrates outwards as it follows the dissipating eccentric feature. Once this feature has completely dissipated, and the surface density gradient in the outer disc is negative once more, the planet then migrates into the inner disc. It finally stops at the edge of the inner cavity, ap =1.1 au, in good agreement with the results from core migration in the 1–5MMSN models presented earlier. Figure 20. Open in new tabDownload slide 1D surface density profiles are shown here during the evolution of the 20MMSN disc, and its embedded core. All line types and colours carry the same meaning as previous plots of this type. The first panel shows the system at the beginning of disc dissipation, the second is 2000 Pb orbits through this process, and the final subplot shows the system once the disc has relaxed and the planet has reached its final stopping position. The middle panel's surface density distribution looks very similar to those seen from 10MMSN, whilst the last panel looks very similar to a 1MMSN model. The planet initially migrates outwards as it follows the dissipating eccentric feature. Once this feature has completely dissipated, and the surface density gradient in the outer disc is negative once more, the planet then migrates into the inner disc. It finally stops at the edge of the inner cavity, ap =1.1 au, in good agreement with the results from core migration in the 1–5MMSN models presented earlier. 5.2 Kepler-34 In contrast to Kepler-16, in the Kepler-34 system we undertook disc dissipation in the 2–20MMSN models, reducing their disc mass to 1MMSN. With these models, we would expect the core semimajor axes and eccentricities to converge on the values reached by the core in the 1MMSN model, as the disc mass dissipates. Examining the top panel of Fig. 21, the evolution of ap, we see this is not the case. Whilst there is evidence of slight outwards migration as a result of the disc relaxing, they do not migrate significantly to ap ≈ 2 au – the location of the 1MMSN model core. The 2MMSN core shows little change, 5MMSN migrates outwards to 1.6 au where it halts, and the high-mass discs all converge to 1.4 au. As the eccentric features in the outer disc dissipate, the core in the second 20MMSN model is able to escape the outer disc and migrate into the inner disc, where it halts at 1.4 au, close to the stopping radius of the other 20MMSN run with reducing disc mass, and the corresponding 10MMSN case. Comparing the surface density maps in Fig. 22, the lack of agreement between the basic 1MMSN migration scenario and the disc-dissipated 20→1MMSN models is clear. Whilst there is some evidence of the disc relaxing during its dissipation, the cavity in the latter model is still more tightly bound around the central binary, as a result the planet is in a much closer orbit than expected. It appears that the presence of the planet in the inner cavity interferes with the relaxation of the disc and prevents it from relaxing to the configuration expected from the 1MMSN run. It is for this reason that we achieve a smaller stopping radius for the planets when the disc mass transitions from high to low mass, and indicates that the history of the system influences the final stopping location of the planet. Figure 21. Open in new tabDownload slide Evolution of protoplanetary cores’ semimajor axes in dissipating self-gravitating discs in the Kepler-34 system. The middle and bottom panels shows these cores’ eccentricity evolution in the low- and high-mass disc models, respectively. Figure 21. Open in new tabDownload slide Evolution of protoplanetary cores’ semimajor axes in dissipating self-gravitating discs in the Kepler-34 system. The middle and bottom panels shows these cores’ eccentricity evolution in the low- and high-mass disc models, respectively. Figure 22. Open in new tabDownload slide Surface density maps of the Kepler-34 1MMSN system once the planet has reached its final orbit (left), and the final state of the 20→1MMSN dissipation model (right), once the disc has relaxed and the planet has migrated outwards to the orbit seen here. Despite the mass being equal in these discs, the cavity in the dissipated case is still much smaller and less eccentric. Figure 22. Open in new tabDownload slide Surface density maps of the Kepler-34 1MMSN system once the planet has reached its final orbit (left), and the final state of the 20→1MMSN dissipation model (right), once the disc has relaxed and the planet has migrated outwards to the orbit seen here. Despite the mass being equal in these discs, the cavity in the dissipated case is still much smaller and less eccentric. 5.3 Kepler-35 Similar results to Kepler-16 in Section 3.1 prompted a similar approach for running disc dissipation scenarios in the Kepler-35 systems; ignoring the low-mass models which show consistent results and focusing on the high-mass models which show the most variation, both with each other and the low-mass cases. The cores in the 10MMSN and 20MMSN models start in much the same positions as those in the Kepler-16 models; the 10MMSN core on a close-in orbit (ap = 0.9 au) with a low eccentricity, and the 20MMSN core trapped in the outer disc by the first additional eccentric feature. Examining Fig. 23 the similarity continues – the cores follow the same migration pathway as their counterpart cores in the Section 5.1 simulations. The replication of the same evolutionary scenarios in different mass-ratio binary systems, suggest that the mechanisms observed in the above sections are relatively robust. For low-mass, isothermal discs, with the same structure, the zero torque location should be the same – hence planets in these discs halt migration at the same location. The core in the dissipating 20MMSN model reaches the same semimajor axis as the core in the 1MMSN model, in our first set of simulations, because when it reaches the inner disc the structure is the same because the disc has already relaxed. On the other hand, the core in 10MMSN has already reached the inner disc, with a different disc structure. As the disc dissipates, the planet also has an impact on the final disc structure, which shifts the zero torque location – inwards in this case. The tendency for the cores in the 10MMSN models to converge on to shorter period orbits between the 8:1 and 9:1 MMRs, resulting in better agreement with the observations than the other models when the planet mass is kept constant, is likely a fingerprint of the high-mass disc structure at 10MMSN, which when it disperses allows the planets to achieve shorter period orbits with low eccentricities. Figure 23. Open in new tabDownload slide Evolution of protoplanetary cores’ semimajor axes in dissipating self-gravitating discs in the Kepler-35 system. The bottom panels show these cores’ eccentricity evolution in the high-mass disc models. Figure 23. Open in new tabDownload slide Evolution of protoplanetary cores’ semimajor axes in dissipating self-gravitating discs in the Kepler-35 system. The bottom panels show these cores’ eccentricity evolution in the high-mass disc models. 6 SUMMARY AND DISCUSSION This is the second paper in a series that examines the influence of disc self-gravity on the evolution of gaseous circumbinary discs, and on the evolution of planets that are embedded in those discs. The focus of Paper I was on the evolution of the discs alone. Several disc masses, ranging between 1 and 20MMSN equivalent discs, were used to probe the evolution of disc structure throughout the lifetime of a circumbinary disc under the influence of self-gravity. The main results to emerge from this study were that self-gravity leads to two important effects: (i) the size of the tidally truncated, eccentric inner disc cavity that forms tends to be smaller for larger disc masses as self-gravity compacts the system scale; (ii) additional precessing eccentric modes emerge at large orbital radii in discs where self-gravity is important. In this paper, we use the end-points of the simulations from Paper I as initial conditions for simulations that examine the orbital evolution of embedded planets, with the binary parameters having been chosen to correspond to the Kepler-16, -34 and -35 planet-hosting systems. Most of the simulations that we present assume that the planet-binary mass ratio is fixed at q = 6 × 10−5. The aim of this work is to examine whether or not self-gravity can improve the level of agreement between the migration stopping locations of planets in the simulations and their currently observed orbital radii. In addition to examining the influence of disc mass, we also examined how the results changed when allowing planets to accrete gas so that they reached their observed masses while migrating, and the influence of allowing the disc mass to decrease with time such that high-mass discs transition to become low-mass discs after the planets have migrated to their stopping radii. We summarize and discuss the results for the different binary systems below. 6.1 Kepler-16 We found that the cavity size in this case only changes significantly when the disc mass exceeds 10MMSN. The migration of planetary cores of fixed mass in the lower mass discs resulted in them stopping close to the edge of the cavity, but with semimajor axes and eccentricities that were too large compared to the observations (ap ∼ 1.15 au and ep ∼ 0.11 versus observed values of 0.705 au and 0.007, respectively). The stopping location, however, was found to agree well with our previous work presented in PN13 that adopted different boundary conditions. Migration in the 20MMSN disc resulted in the planet being halted by one of the additional eccentric features further out in the disc that acted as a planet trap, so in spite of the disc cavity being significantly smaller in this case, the planet was unable to reach the cavity such that it could park closer to the central binary as required by the observations. Allowing the planets to accrete gas so that they reach the mass inferred from observations (this requires the mass ratio to grow from 6 × 10−5 to 3.54 × 10−4) resulted in much better agreement with observations for the low-mass discs. Here, the planet grows in excess of the gap forming mass, and this allows it to push deeper into the tidally truncated cavity. Furthermore, the growth of the planet causes the eccentricity of the central cavity (and the other eccentric features) to diminish significantly, and this leads to the eccentricity of the planet orbits reducing significantly. For the lowest mass disc we obtain ap = 0.78 au and ep = 0.03, which agrees rather well with the observed values for Kepler-16b, giving us confidence that the formation and evolution scenario that we are exploring in this work is probably the correct one. The evolution of the accreting planets in the higher mass 10MMSN and 20MMSN discs did not result in such good agreement with observations. Here, the planets have their eccentricities excited by the eccentric disc modes to values that cause them to interact more strongly with the central binary, leading to a sequence of scattering events that send them out into the disc and then back again over the full run times of the simulations. Finally, allowing the disc mass to decrease for the heavy 10MMSN and 20MMSN discs, while keeping the planet-binary mass ratio =6 × 10−5 caused the planets to end up orbiting closer to the star than when the disc masses were at their initial values. In particular, the reduction in disc mass causes the additional eccentric features in the disc to dissipate, and this allows the planet in the 20MMSN case to migrate inwards. The level of agreement with observations in these cases, however, is not as good as that obtained by allowing the planet masses to increase to their observed values. We conclude that for the Kepler-16 system, self-gravity of the disc does not provide a positive contribution to obtaining agreement between observations and theoretical predictions. 6.2 Kepler-34 The binary system in Kepler-34 has an eccentric orbit, and this leads to the formation of a wide and highly eccentric cavity when the disc mass is low. For large disc masses, however, self-gravity causes the disc cavity to shrink substantially, and this has a strong influence on the orbital evolution of embedded planets. In the low-mass discs the planets migrate inwards and stop at the edge of the cavity, which is too far from the binary for the stopping location to agree with the observations (ap  ∼  1.8 au and ep  ∼ 0.28 versus the observed values 1.09 au and 0.182, respectively). In the high-mass cases, however, we find that the additional eccentric features that form in the disc are somewhat weaker than in the Kepler-16 run described above, and consequently the planets can normally migrate all the way to the central cavity in these cases. We find that the semimajor axes and eccentricities of the planets in the 10MMSN and 20MMSN discs straddle the observed values for Kepler-34b, indicating that self-gravity in this case provides the possibility of obtaining much better agreement with the observations. Switching on gas accretion makes very little difference to the results of these simulations because the final mass of Kepler-34b is only 60 per cent larger than the initial mass that we start with. Allowing the disc mass to decrease inevitably leads to the cavity sizes of the most massive discs increasing as the influence of self-gravity is diminished. Interestingly, however, we find that the presence of the planet prevents the cavity from relaxing to the size expected for a lower mass disc, and instead the amount of expansion observed is relatively modest. (The time to establish the approximately steady-state cavity configuration, in the absence of planets, is typically ∼3500 binary orbits. We have run our simulations for longer than this to ensure that we have achieved a quasi-steady state.) Although the planets in these more massive discs no longer show such good agreement with observations once the disc mass has diminished, they provide much better agreement than those planets that form and migrate in low-mass discs. This leads us to conclude that self-gravity can have a positive impact on obtaining agreement between simulations and the observations of Kepler-34b because of the rather dramatic influence that it has on the cavity size, and also because the system retains memory of its larger initial disc mass when the mass of the disc is slowly decreased. Formation of a planet in a heavy disc, followed by its migration and then rapid disc removal would seem to provide one way in which agreement with observations could be obtained for this system. 6.3 Kepler-35 As mentioned previously in this paper, the similar eccentricity of the Kepler-16 and -35 binaries leads to very similar outcomes both in terms of disc structure and orbital evolution. Mass growth of Kepler-35b from the initial planet-binary mass ratio of 6 × 10−5 was not considered in this paper because the final mass is only 20 per cent larger than the initial mass. One consequence of this lower planet mass is that growth to a gap forming object that can push further into the inner cavity is difficult to invoke so that good agreement between the simulation outcomes and the observations of Kepler-35b can be obtained, in contrast to the situation with Kepler-16b. The final orbital radii of Kepler-35b analogues were always too large by a factor of 1.5 compared to the observed values. Allowing the planet to be in the partial gap forming regime, such that it might push deeper into the cavity, can probably only be achieved by a significant reduction in the disc pressure scaleheight. Given that the two stars in Kepler-35 are more massive and hotter than in the Kepler-16, it is not immediately obvious why the disc should be cooler in this case. Fitting this system using simulations therefore remains an unsolved problem, and will require a more sophisticated treatment of the disc thermodynamics to examine whether or not Kepler-35b could have been in the gap forming regime when the protoplanetary disc was present. Assuming that the scenario we have explored in this paper, namely that planets form at large orbital radii in circumbinary discs, and then migrate inwards to be stopped near the cavity edge, is the correct one, then we can use the Kepler-16, -34 and -35 systems to constrain models of planet formation and protoplanetary disc dynamics. The simplicity of our models means that there is still a large amount of work that needs to be done to achieve this goal. In this work, gas accretion on to planets to their final masses and disc dissipation scenarios have been carried out separately. Combining the two, akin to the method in Pierens & Nelson (2013), in self-gravitating discs might help to fit the orbital properties of the planets but also shed light on the era in which the planets may have accreted their masses. Whilst disc-mass, and the influence of self-gravity can significantly alter disc structure, other physics will also play important roles. The inclusion of an adiabatic equation of state with radiative physics has been investigated in circumbinary systems, along with the effect on planet migration (Kley & Haghighipour 2014, 2015). The non-uniform, time dependent, radiation field produced by the two stars, however, has not yet been explored in combination with a more realistic thermal treatment of the disc. 3D effects are also likely to be important. These include, but are not limited to: disc warping when the disc and orbit plane of the central binary are misaligned (Larwood & Papaloizou 1997), which can in turn lead to the development of a parametric instability in the disc that may be a source of hydrodynamic turbulence (Ogilvie & Latter 2013); the development of eccentric modes in the disc leading to parametric instability and hydrodynamic turbulence (Papaloizou 2005; Barker & Ogilvie 2014) and the formation of Spiral Wave Instabilities existing in a disc that is tidally forced by a binary system, can lead to a parametric instability and hydrodynamic turbulence (Bae et al. 2016). Equally importantly will be the inclusion of MHD, since the underlying angular momentum transport mechanism operating in circumbinary discs is likely to be of magnetic origin (e.g. Balbus & Hawley 1991; Bai & Stone 2013). Simulations carried out in 3D will allow us to examine these and other effects. Finally, both theory and observations indicate that planets do not normally form in isolation, so the evolution of multiplanet systems (Kley & Haghighipour 2015) may provide better agreement with at least a subset of circumbinary planet observations. Acknowledgments The authors thank the referee for the comments and suggestions in their helpful report. Computer time for the simulations performed with GENESIS was provided by HPC resources of Cines under the allocation A0010406957 made by GENCI (Grand Equipement National de Calcul Intensif). Those performed with FARGO utilized: Queen Mary's MidPlus computational facilities, supported by QMUL Research-IT and funded by EPSRC grant EP/K000128/1; and the DiRAC Complexity system, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment is funded by BIS National E-Infrastructure capital grant ST/K000373/1 and STFC DiRAC Operations grant ST/K0003259/1. DiRAC is part of the National E-Infrastructure. This research was also supported in part by the National Science Foundation under Grant No. NSF PHY-1125915. 1 " Doyle et al. (2011). 2 " Welsh et al. (2012). 3 " Orosz et al. (2012a). 4 " Orosz et al. (2012b). 5 " Welsh et al. (2015a). 6 " Kostov et al. (2013). 7 " Kostov et al. (2014). 8 " Welsh et al. (2015b). 9 " Kostov et al. (2016). 10 " In principle, we should also include an additional term from the disc acting on the centre of mass of the binary. Extensive tests were undertaken in Paper I which demonstrated that including this term did not change the results. REFERENCES Artymowicz P. , Lubow S. H., 1994 , ApJ , 421 , 651 Crossref Search ADS Bae J. , Nelson R. P., Hartmann L., Richard S., 2016 , ApJ , 829 , 13 Crossref Search ADS Bai X.-N. , Stone J. M., 2013 , ApJ , 769 , 76 Crossref Search ADS Balbus S. A. , Hawley J. F., 1991 , ApJ , 376 , 214 Crossref Search ADS Barker A. J. , Ogilvie G. I., 2014 , MNRAS , 445 , 2637 Crossref Search ADS Baruteau C. , Masset F., 2008a , ApJ , 672 , 1054 Crossref Search ADS Baruteau C. , Masset F., 2008b , ApJ , 678 , 483 Crossref Search ADS Bromley B. C. , Kenyon S. J., 2015 , ApJ , 806 , A98 Crossref Search ADS Crida A. , Morbidelli A., Masset F., 2006 , Icarus , 181 , 587 Crossref Search ADS Doyle L. R. et al. , 2011 , Science , 333 , 1602 Crossref Search ADS PubMed Dunhill A. C. , Alexander R. D., 2013 , MNRAS , 435 , 2328 Crossref Search ADS Dutrey A. , Guilloteau S., Simon M., 1994 , A&A , 286 , 149 Fendyke S. M. , Nelson R. P., 2014 , MNRAS , 437 , 96 Crossref Search ADS Hayashi C. , 1981 , Prog. Theor. Phys. Suppl. , 70 , 35 Crossref Search ADS Holman M. J. , Wiegert P. A., 1999 , AJ , 117 , 621 Crossref Search ADS Kley W. , 1999 , MNRAS , 303 , 696 Crossref Search ADS Kley W. , Haghighipour N., 2014 , A&A , 564 , A72 (KH14) Crossref Search ADS Kley W. , Haghighipour N., 2015 , A&A , 581 , A20 (KH15) Crossref Search ADS Kley W. , Nelson R. P., 2010 , Planets in Binary Star Systems . Springer Netherlands , Dordrecht , . 135 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Kostov V. B. , McCullough P. R., Hinse T. C., Tsvetanov Z. I., Hébrard G., Díaz R. F., Deleuil M., Valenti J. A., 2013 , ApJ , 770 , 52 Crossref Search ADS Kostov V. B. et al. , 2014 , ApJ , 787 , A14 Crossref Search ADS Kostov V. B. et al. , 2016 , ApJ , 827 , 86 Crossref Search ADS Larwood J. D. , Papaloizou J. C. B., 1997 , MNRAS , 285 , 288 Crossref Search ADS Lin D. N. C. , Papaloizou J. C. B., 1986 , ApJ , 309 , 846 Crossref Search ADS Lines S. , Leinhardt Z. M., Paardekooper S.-J., Baruteau C., Thebault P., 2014 , ApJ , 782 , L11 Crossref Search ADS Lines S. , Leinhardt Z. M., Baruteau C., Paardekooper S.-J., Carter P. J., 2015 , A&A , 582 , A5 Crossref Search ADS Marzari F. , Thébault P., Scholl H., 2008 , ApJ , 681 , 1599 Crossref Search ADS Marzari F. , Scholl H., Thébault P., Baruteau C., 2009 , A&A , 508 , 1493 Crossref Search ADS Masset F. , 1999 , A&AS , 141 , 165 Crossref Search ADS Masset F. S. , 2006 , ApJ , 652 , 730 Crossref Search ADS Meschiari S. , 2012a , ApJ , 752 , A72 Crossref Search ADS Meschiari S. , 2012b , ApJ , 761 , L7 Crossref Search ADS Mutter M. M. , Pierens A., Nelson R. P., 2017 , MNRAS , 465 , 4735 Crossref Search ADS Nelson R. P. , 2003 , MNRAS , 345 , 233 Crossref Search ADS Nelson R. P. , Papaloizou J. C. B., Masset F., Kley W., 2000 , MNRAS , 318 , 18 Crossref Search ADS Ogilvie G. I. , Latter H. N., 2013 , MNRAS , 433 , 2420 Crossref Search ADS Orosz J. A. et al. , 2012a , Science , 337 , 1511 Crossref Search ADS Orosz J. A. et al. , 2012b , ApJ , 758 , A87 Crossref Search ADS Paardekooper S.-J. , Leinhardt Z. M., Thébault P., Baruteau C., 2012 , ApJ , 754 , L16 Crossref Search ADS Papaloizou J. C. B. , 2002 , A&A , 388 , 615 Crossref Search ADS Papaloizou J. C. B. , 2005 , A&A , 432 , 743 Crossref Search ADS Pelupessy F. I. , Portegies Zwart S., 2013 , MNRAS , 429 , 895 Crossref Search ADS Pierens A. , Nelson R. P., 2007 , A&A , 472 , 993 Crossref Search ADS Pierens A. , Nelson R. P., 2008a , A&A , 478 , 939 Crossref Search ADS Pierens A. , Nelson R. P., 2008b , A&A , 483 , 633 Crossref Search ADS Pierens A. , Nelson R. P., 2013 , A&A , 556 , A134 (PN13) Crossref Search ADS Press W. H. , Teukolsky S. A., Vetterling W. T., Flannery B. P., 1992 , Numerical Recipes in C. The Art of Scientific Computing , 2nd edn . University Press , Cambridge Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Scholl H. , Marzari F., Thébault P., 2007 , MNRAS , 380 , 1119 Crossref Search ADS Shakura N. I. , Sunyaev R. A., 1973 , A&A , 24 , 337 Tanaka H. , Taku T., Ward W. R., 2002 , ApJ , 565 , 1257 Crossref Search ADS van Leer B. , 1977 , J. Comput. Phys. , 23 , 276 Crossref Search ADS Ward W. , 1997 , Icarus , 126 , 261 Crossref Search ADS Welsh W. F. et al. , 2012 , Nature , 481 , 475 Crossref Search ADS PubMed Welsh W. F. , Orosz J., Quarles B., Haghighipour N., 2015a , ESS Meeting , 3 , 402.01 Welsh W. F. et al. , 2015b , ApJ , 809 , 26 Crossref Search ADS © 2017 The Authors Published by Oxford University Press on behalf of the Royal Astronomical Society
TI - The role of disc self-gravity in circumbinary planet systems – II. Planet evolution
JO - Monthly Notices of the Royal Astronomical Society
DO - 10.1093/mnras/stx1113
DA - 2017-08-21
UR - https://www.deepdyve.com/lp/oxford-university-press/the-role-of-disc-self-gravity-in-circumbinary-planet-systems-ii-planet-Da2laqNRmj
SP - 4504
VL - 469
IS - 4
DP - DeepDyve
ER -