Issue 
A&A
Volume 664, August 2022



Article Number  A87  
Number of page(s)  15  
Section  Celestial mechanics and astrometry  
DOI  https://doi.org/10.1051/00046361/202243344  
Published online  11 August 2022 
Resonance capture and longterm evolution of planets in binary star systems
NaXys, Department of Mathematics, University of Namur,
61 Rue de Bruxelles,
5000
Namur, Belgium
email: arnaud.roisin@unamur.be
Received:
16
February
2022
Accepted:
13
June
2022
Aims. The growing population of planets discovered in orbit around one stellar component of a binary star raises the question of the influence of the binary companion on the formation process of planetary systems. The aim of this work is to study the impact of a binary companion on the evolution of twoplanet systems during both the TypeII migration phase and their longterm evolution after the dissipation of the protoplanetary disk.
Methods. We used the symplectic integrator SyMBA, modified to include a wide binary companion. We also included the TypeII migration of giant planets during the protoplanetary disk phase with suitable eccentricity and inclination damping as well as the gravitational potential acting on the planets due to the disk and the nodal precession of the disk induced by the binary companion. We considered various inclinations, eccentricities, and separations of the binary companion.
Results. Disk migration allows for the formation of planet pairs in meanmotion resonances despite the presence of the binary companion. When the binary separation is wide (1000 au), the timescale of the perturbations that it raises on the planets is longer than the disk’s lifetime and resonant pairs are routinely formed in the 2:1, 5:2, and 3:1 commensurabilities. Provided the planetplanet interaction timescale is smaller than the timescale of binary perturbations, these systems can remain in resonance long after the disk has dissipated. When the binary separation is smaller (250 au), only planets in the 2:1 resonance tend to remain in a resonant state and more chaotic evolutions are observed, as well as more ejections. After those ejections, the remaining planet can become eccentric due to the perturbations from the binary companion in addition, for strongly inclined binary companions, captures in the von ZiepelLidovKozai resonance can occur. Whereas in systems with two planets, this mechanism is quenched by planetplanet interactions. Our simulations reveal that the interplay between planetdisk, planetplanet, and planetbinary interactions can lead to the formation of resonant pairs of planets which remain stable over timescales much longer than the disk’s lifetime.
Key words: planetdisk interactions / planetstar interactions / binaries: general / planets and satellites: dynamical evolution and stability / planets and satellites: formation
© A. Roisin et al. 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
1 Introduction
Planetary system formation under the influence of a binary companion is an important question since it is estimated that about half of the Sunlike stars are part of multiple star systems (Duquennoy & Mayor 1991; Raghavan et al. 2010). Until now, more than 100 Stype planets (also called circumprimary planets) have been discovered and their eccentricities are slightly higher than around single star systems (Kaib et al. 2013).
A binary companion plays an important role in the formation process. Firstly, it has a considerable impact on the protoplanetary disk. Regarding the disk formation, it was observed that disks in close binary star systems are less present and if they are, they tend to be less massive than in single star systems (Kraus & Ireland 2012; Franchini et al. 2019; Chachan et al. 2019). The binary companion can induce a truncation in the disk, as shown by Artymowicz & Lubow (1994) and Savonije et al. (1994). The binary companion also influences the shape of the disk (Terquem & Bertout 1993; Papaloizou & Terquem 1995) and could even lead to the warping of the disk if the self gravity or pressure forces are not strong enough to maintain its uniformity (Batygin et al. 2011; Zanazzi & Lai 2017). Furthermore, the disk precesses nodally due to the presence of the binary companion (Batygin et al. 2011; Roisin et al. 2021).
Secondly, the binary companion also influences the planet embedded in the protoplanetary disk. The planetdisk interactions were extensively studied for close binaries by means of 3D hydrodynamical simulations where a decoupling between the motion of the planet and the disk was generally observed (e.g., XiangGruess & Papaloizou 2014; Picogna & Marzari 2015; Lubow & Martin 2016; Martin et al. 2016). Regarding wide binaries, it was shown that the disk gravitational potential acting on the planet and the damping forces exerted by the disk on the planet tend to keep the latter in the midplane of the former (e.g., Roisin et al. 2021).
One key effect affecting the planets in binary star systems is the von ZeipelLidovKozai (ZLK) resonance which consists of coupled eccentricity and inclination variations for the planets (von Zeipel 1910; Lidov 1962; Kozai 1962). This effect influences the formation of giant planets in binary stars in multiple ways. It has been shown that the binary companion can generate ZLK cycles in the disk (Martin et al. 2014; Fu et al. 2015a), which could be suppressed by the selfgravity of the disk (Batygin et al. 2011; Rafikov & Silsbee 2015). Moreover, for sufficiently inclined disks, ZLK instabilities can arise (Lubow & Ogilvie 2017; Zanazzi & Lai 2017). Inclined binary companions also influence planets migrating in the protoplanetary disk via the ZLK effect (Roisin & Libert 2021), but this effect is counteracted by the gravitational potential of the disk acting on the planets in the case of wide binaries (Roisin et al. 2021). Finally, it has been widely shown that the ZLK mechanism, when associated with tidal friction, could lead to the formation of hot Jupiters in binary star systems (e.g., Fabrycky & Tremaine 2007).
Most of the planet formation studies focus on close binary systems, due to the strong interaction between the planets and the binary companion. However a wide binary could also in some case have a significant impact on the Stype planetary system architecture. This work aims at extending the studies on single planet systems made by Roisin & Libert (2021) and Roisin et al. (2021) by considering the influence of a wide binary companion on the resonance captures of two giant planets migrating in the TypeII regime during the protoplanetary disk phase where the gravitational influence of the disk on the planets and the mutual gravitational influence between the planets also act. The longterm evolution of the systems after the dispersal of the disk will be thoroughly investigated. We also focus our attention on the ZLK resonance with the wide binary companion and how it could play a role during the two stages mentioned above.
The paper is organized as follows. Section 2 describes the setup of the Nbody simulations carried out in this work. Section 3 presents an overview of the relevant timescales of the problem. Section 4 presents the results of the simulations for a wide binary companion at 1000 au by outlining the system architectures as well as the different dynamical effects observed during and after the disk phase, while a similar analysis is performed in Sect. 5 for closer binary stars. The longterm stability of the systems found in our simulations is assessed in Sect. 6. In Sect. 7, we discuss our results in light of the observational data on Stype planets. Finally our conclusions are given in Sect. 8.
2 Methods
In this work, we follow the evolution of two planets perturbed by a binary companion during and after the disk phase. Different semimajor axis a_{B}, eccentricity e_{B}, and inclination i_{B} values are considered for the binary. The masses of both stars are 1 M_{⊙} in all the simulations.
To run those simulations, the code used in this work is the wellknown Nbody SyMBA integrator, which handles close encounters (Duncan et al. 1998). It was modified to include a wide binary companion following Chambers et al. (2002; see Roisin & Libert 2021, for more details). The code also considers the gravitational influence of the protoplanetary disk in which the planets are embedded, namely the typeII migration of the planets (Ivanov et al. 1999; Nelson et al. 2000; Crida & Morbidelli 2007), through the acceleration (Papaloizou & Larwood 2000) $${a}_{\mathrm{mig}}=\frac{{v}_{\text{pl}}}{{\tau}_{\text{mig}}},$$(1)
with v_{pl} the velocity of the migrating planet and τ_{mig} the timescale for the TypeII regime given by (both in the diskdominated case and planetdominated cases, see Sotiriadis et al. 2017, for more details): $${\tau}_{\text{mig}}=\frac{2}{3}{\alpha}^{1}{h}^{2}{\text{\Omega}}_{\text{pl}}^{1}\times \mathrm{max}\left\{1,\frac{{m}_{\text{pl}}}{\left(4\pi /3\right)\Sigma \left({r}_{\text{pl}}\right){r}_{\text{pl}}^{2}}\right\},$$(2)
with α = 0.005 the classical value for the ShakuraSunyanev viscosity parameter (Shakura & Sunyaev 1973), h = 0.05 the disk aspect ratio, ${\text{\Omega}}_{\text{pl}}^{1}$ the orbital frequency of the planet, a_{pl} the semimajor axis of the planet, m_{pl} the mass of the planet, r_{pl} the distance of the planet to the host star, and Σ ∝ r^{−05} the surface density profile of the disk. Unless otherwise stated, in the simulations, we fix the disk inner and outer edges to a_{in} = 0.05 AU and a_{out} = 30 AU. The code also adopts a smooth transition in the gasfree inner cavity using an hyperbolic tangent function $\mathrm{tanh}\left(\frac{r{R}_{\text{in}}}{\text{\Delta}r}\right)$ where ∆r = 0.001 AU, following Matsumoto et al. (2012).
Moreover, the inclination and eccentricity damping induced on the planets by the disk is included through the formulas deduced from the 3D hydrodynamical simulations in Bitsch et al. (2013): $$\begin{array}{ll}\frac{\text{d}e}{\text{d}t}\left({m}_{\text{pl}},{e}_{\text{pl}},{i}_{\text{pl}}\right)=\hfill & \frac{{m}_{\text{ld}}}{0.025{m}_{A}}{\left({a}_{e}{\left[{i}_{\text{pl}}+\frac{{m}_{\text{pl}}}{3}\right]}^{2{b}_{e}}+{C}_{e}{i}_{\text{pl}}^{2{d}_{e}}\right)}^{\frac{1}{2}}\hfill \\ \hfill & +12.65\frac{{m}_{\text{pl}}{m}_{\text{ld}}}{{m}_{A}^{2}}{e}_{\text{pl}}\mathrm{exp}\left({\left[\frac{{i}_{\text{pl}}}{{m}_{\text{pl}}}\right]}^{2}\right),\hfill \\ \frac{\text{d}i}{\text{d}t}\left({m}_{\text{pl}},{e}_{\text{pl}},{i}_{\text{pl}}\right)=\hfill & \frac{{m}_{\text{ld}}}{0.025{m}_{A}}({a}_{i}{i}_{\text{pl}}^{2{b}_{i}}\mathrm{exp}\left(\frac{{\left({i}_{\text{pl}}/gi\right)}^{2}}{2}\right)\hfill \\ \hfill & {+{c}_{i}{\left[\frac{{i}_{\text{pl}}}{40}\right]}^{2{d}_{i}})}^{\frac{1}{2}},\hfill \end{array}$$(3)
with e_{pl} and i_{pl} the planetary eccentricity and inclination, respectively, m_{A} the mass of the central star, and with the coefficients defined as $${a}_{e}=80{e}_{\text{pl}}^{2}\mathrm{exp}\left\{{e}_{\text{pl}}^{2}\frac{{m}_{\text{pl}}}{0.26}\right\}{15}^{{m}_{\text{pl}}}\left(20+11{m}_{\text{pl}}{m}_{\text{pl}}^{2}\right),$$(4) $${b}_{e}=0.3{m}_{\text{pl}},$$(5) $${c}_{e}=450+{2}^{{m}_{\text{pl}}},$$(6) $${d}_{e}=1.4+\frac{\sqrt{{m}_{\text{pl}}}}{6}.$$(7) $${a}_{i}=1.5\times {10}^{4}\left(23{e}_{\text{pl}}\right){m}_{\text{pl}}^{3},$$(8) $${b}_{i}=1+\frac{{m}_{\text{pl}}{e}_{\text{pl}}^{2}}{10},$$(9) $${c}_{i}=\frac{1.2\times {10}^{6}}{\left(23{e}_{\text{pl}}\right)\left(5+{e}_{\text{pl}}^{2}{\left[{m}_{\text{pl}}+2\right]}^{3}\right)},$$(10) $${d}_{i}=3+2{e}_{\text{pl}},$$(11)
and $${g}_{i}=\sqrt{\frac{3{m}_{{}_{\text{pl}}}}{{e}_{{}_{\text{pl}}}+0.001}.}$$(12)
To include the damping in the code, we converted it into an acceleration according to the formulas (Papaloizou & Larwood 2000) $${a}_{\text{ecc}}=2\frac{\left({v}_{\text{pl}}\cdot {r}_{\text{pl}}\right){r}_{\text{pl}}}{{\Vert {r}_{\text{pl}}\Vert}^{2}{\tau}_{\text{ecc}}},$$(13) $${a}_{\text{inc}}=2\frac{\left({v}_{\text{pl}}\cdot k\right)k}{{\tau}_{\text{inc}}},$$(14)
where v_{pl} and r_{pl} are the velocity and the position of the planet, respectively, k is the unitary vertical vector, $${\tau}_{\text{ecc}}=\left\frac{{e}_{\text{pl}}}{\text{d}e/\text{d}t}\right,$$(15)
and $${\tau}_{\text{inc}}=\left\frac{{i}_{\text{pl}}}{\text{d}i/\text{d}t}\right.$$(16)
In Eq. (3), the time is expressed in orbital period and the inclination damping in degrees per orbit. We refer to Sotiriadis et al. (2017) for more details.
The code also includes the disk gravitational potential acting on the planets expressed in spherical coordinates as (Terquem & Ajmia 2010; Huré & Hersant 2011) $$\Phi \left({r}_{\text{pl}},{\phi}_{\text{pl}},{\theta}_{\text{pl}}\right)=2G{\displaystyle {\int}_{{R}_{\text{in}}}^{{R}_{\text{out}}}\Sigma \left(\tilde{r}\right)\sqrt{\frac{\tilde{r}}{{r}_{\text{pl}}\mathrm{sin}\left({\theta}_{\text{pl}}\right)}}}kK\left(k\right)\text{d}\tilde{r},$$(17)
where, (r_{pl}, φ_{pl}, θ_{pl}) is the position of the planet in spherical coordinates and K(k) is the complete elliptic integral of first kind: $$K\left(k\right)={\displaystyle {\int}_{0}^{\pi /2}\frac{\text{d}\tilde{\phi}}{\sqrt{1{k}^{2}{\mathrm{sin}}^{2}\left(\tilde{\phi}\right)}}}$$(18)
with k between 0 and 1 and equal, in our case, to $$k=2\sqrt{\frac{{r}_{\text{pl}}\tilde{r}\mathrm{sin}\left({\theta}_{\text{pl}}\right)}{{r}_{\text{pl}}^{2}+{\tilde{r}}^{2}+2{r}_{\text{pl}}\tilde{r}\mathrm{sin}\left({\theta}_{\text{pl}}\right)}}.$$(19)
We converted this potential in terms of acceleration in cartesian coordinates through (Roisin et al. 2021) $$\begin{array}{l}\frac{{d}^{2}{x}_{\text{pl}}}{d{t}^{2}}=\mathrm{sin}\left({\theta}_{\text{pl}}\right)\mathrm{cos}\left({\phi}_{\text{pl}}\right)\frac{\partial \Phi}{\partial {r}_{\text{pl}}}\mathrm{cos}\left({\theta}_{\text{pl}}\right)\mathrm{cos}\left({\phi}_{\text{pl}}\right)\frac{1}{{r}_{{}_{\text{pl}}}}\frac{\partial \Phi}{\partial {\theta}_{{}_{\text{pl}}}}\hfill \\ \frac{{d}^{2}{y}_{\text{pl}}}{d{t}^{2}}=\mathrm{sin}\left({\theta}_{\text{pl}}\right)sin\left({\phi}_{\text{pl}}\right)\frac{\partial \Phi}{\partial {r}_{\text{pl}}}\mathrm{cos}\left({\theta}_{\text{pl}}\right)sin\left({\phi}_{\text{pl}}\right)\frac{1}{{r}_{{}_{\text{pl}}}}\frac{\partial \Phi}{\partial {\theta}_{{}_{\text{pl}}}}\hfill \\ \frac{{d}^{2}{z}_{\text{pl}}}{d{t}^{2}}=\mathrm{cos}\left({\theta}_{\text{pl}}\right)\frac{\partial \Phi}{\partial {r}_{\text{pl}}}+sin\left({\theta}_{\text{pl}}\right)\frac{1}{{r}_{{}_{\text{pl}}}}\frac{\partial \Phi}{\partial {\theta}_{{}_{\text{pl}}}},\hfill \end{array}$$(20)
where (x_{pl}, y_{pl}, z_{pl}) is the cartesian coordinate of the planet and where $$\begin{array}{ll}\frac{\partial \Phi}{\partial {r}_{\text{pl}}}=\hfill & \frac{G}{{r}_{\text{pl}}}{\displaystyle {\int}_{{R}_{\text{in}}}^{{R}_{\text{out}}}\Sigma \left(\tilde{r}\right)\sqrt{\frac{\tilde{r}}{{r}_{\text{pl}}\mathrm{sin}\left({\theta}_{\text{pl}}\right)}}k\left[K\left(k\right)\frac{{\tilde{r}}^{2}{r}_{\text{pl}}^{2}}{{R}^{2}}E\left(k\right)\right]d\tilde{r}}\hfill \\ \frac{\partial \Phi}{\partial {r}_{\text{pl}}}=\hfill & G{\displaystyle {\int}_{{R}_{\text{in}}}^{{R}_{\text{out}}}\Sigma \left(\tilde{r}\right)\mathrm{cot}\left({\theta}_{\text{pl}}\right)\sqrt{\frac{\tilde{r}}{{r}_{\text{pl}}\mathrm{sin}\left({\theta}_{\text{pl}}\right)}}k[K\left(k\right)}\hfill \\ \hfill & \frac{{\tilde{r}}^{2}+{r}_{\text{pl}}^{2}}{{R}^{2}}E\left(k\right)]d\tilde{r}\hfill \end{array}$$(21)
with ${R}^{2}={r}_{\text{pl}}^{2}+{\tilde{r}}^{2}2{r}_{\text{pl}}\tilde{r}\text{\hspace{0.17em}}\mathrm{sin}\left({\theta}_{\text{pl}}\right)$ and E(k) the complete elliptic integral of the second kind: $$E\left(k\right)={\displaystyle {\int}_{0}^{\pi /2}\sqrt{1{k}^{2}{\mathrm{sin}}^{2}\left(\tilde{\phi}\right)}d\tilde{\phi}.}$$(22)
We note that all the accelerations presented in this section have been added in the code preserving the symmetry of SyMBA as prescribed in Lee & Peale (2002): $$f\left({t}_{k+1}\right)={E}_{\text{disk}}\left(\frac{\tau}{2}\right)\mathrm{exp}\left(T\tau \right){E}_{\text{disk}}\left(\frac{\tau}{2}\right)f\left({t}_{k}\right),$$(23)
where T is the linear operator associated with the Hamiltonian H of the Nbody system without the disk and E_{disk} (·) is the evolution due to the disk.
Moreover, we considered the nodal precession of the disk induced by the binary companion as in Roisin et al. (2021), which leads to a uniform precession of the disk’s ascending node and a variation of the disk inclination with amplitude 2i_{B} with respect to the initial plane of the disk: $$\{\begin{array}{c}{i}_{d}\left(t\right)=2{i}_{B}\leftsin\left(2\pi {f}_{d}t\right)\right\\ {\text{\Omega}}_{d}\left(t\right)=\frac{\pi}{2}+{\text{\Omega}}_{B}\left(2\pi {f}_{d}t\text{\hspace{0.17em}}\mathrm{mod}\text{\hspace{0.17em}}\pi \right)\end{array},$$(24)
where i_{d} and Ω_{d} are the inclination and longitude of the ascending node of the disk, i_{B} and Ω_{β} are the inclination and longitude of the ascending node of the binary companion in the initial plane of the disk, and $${f}_{d}=0.2145{\left(\frac{{a}_{\text{out}}}{\text{1au}}\right)}^{3/2}{\left(\frac{{a}_{B}}{\text{1au}}\right)}^{3}\left(\frac{{m}_{B}}{1{m}_{\odot}}\right),$$(25)
with a_{B} and m_{B} being the semimajor axis and mass of the binary companion, respectively.
In our simulations, the disk mass decreases exponentially with time. We consider the disk as being fully dissipated and thus we neglect its influence when $$\frac{\text{d}{m}_{\text{disk}}\left(t\right)}{\text{d}t}=\frac{{m}_{\text{disk},0}}{{T}_{0}}\mathrm{exp}\left(t/{T}_{0}\right)<{10}^{9}{m}_{0}/\text{yr,}$$(26)
where m_{0} denotes the mass of the central star, m_{disk,0} the disk mass at t = 0 (fixed here to 8 m_{Jup}), and T_{0} = 2.8 × 10^{5} yr. More information about the code can also be found in Roisin et al. (2021).
As commonly adopted, we applied an inward migration to the outer planet only (this approach favors convergent migration), while the eccentricity and inclination damping due to the disk is applied on both planets. For the integration, we used a time step of 0.005 yr for the disk phase and 0.01 yr for the dynamical evolution after the dispersal of the disk. In both cases, a recursive reduction of the timestep is implemented in case of close encounters, as described in Duncan et al. (1998). The simulations were carried out for 100 Myr. In the following, subscripts 1 and 2 will refer to the inner planet and outer planet, respectively.
3 Timescale Analysis
There are several timescales relevant to our study: the disk’s lifetime, the timescale of the planetplanet interactions, and the timescale of the precession induced by the binary (i.e., the ZLK resonance timescale). Since we use them to interpret the outcomes of our simulations throughout the paper, we review them first.
All effects inducing a precession of the pericenter of a planet faster than the one induced by the ZLK resonance act to suppress it. It is the case with, for example, general relativity, tidal and rotational bulges, and, as is the case here, gravitational interaction with another planet (Wu & Murray 2003; Fabrycky & Tremaine 2007; Takeda et al. 2008). The timescale associated with the ZLK resonance induced by a binary companion, τ_{zlk}, can be approximated by (Kiseleva et al. 1998) $${\tau}_{\text{ZLK}}=\frac{2}{3\pi}\frac{{P}_{B}^{2}}{{P}_{j}}{\left(1{e}_{B}^{2}\right)}^{3/2}\frac{{m}_{0}+{m}_{j}+{m}_{B}}{{m}_{B}},$$(27)
where P_{B} and Pj are the orbital periods of the binary companion and the jth planet, respectively, and mj the mass of the jth planet (j = 1, 2). If the ZLK timescale if longer than the disk’s lifetime, then the binary companion does not significantly perturb the orbits of the planets while they migrate in the disk.
The timescale of the pericenter precession due to the gravitational interaction between two planets, denoted τ_{pl}, can be estimated by using the LaplaceLagrange secular theory. It can be approximated by means of the eigenvalues of the following matrix (Brouwer & Clemence 1961; Murray & Dermott 1999; Takeda et al. 2008): $$A=\left(\begin{array}{cc}{c}_{1}& {c}_{0}{c}_{1}\\ {c}_{0}{c}_{2}& {c}_{2}\end{array}\right)$$(28)
with $$\begin{array}{l}{c}_{0}=\frac{{b}_{3/2}^{\left(2\right)}\left(\alpha \right)}{{b}_{3/2}^{\left(1\right)}\left(\alpha \right)}\hfill \\ {c}_{1}=\frac{1}{4}{n}_{1}\frac{{m}_{2}}{{m}_{0}+{m}_{1}}{\alpha}^{2}{b}_{3/2}^{\left(1\right)}\left(\alpha \right)\hfill \\ {c}_{2}=\frac{1}{4}{n}_{2}\frac{{m}_{1}}{{m}_{0}+{m}_{2}}{\alpha}^{2}{b}_{3/2}^{\left(1\right)}\left(\alpha \right),\hfill \end{array}$$(29)
where ${b}_{3/2}^{\left(1\right)}\left(\alpha \right)$ and ${b}_{3/2}^{\left(2\right)}\left(\alpha \right)$ are the Laplace coefficients of the first and second kind, respectively, α = a_{1}/a_{2} the planetary semimajor axis ratio, and n_{j} the planetary mean motions (j = 1, 2). The expressions of the eigenvalues are (Zhou & Sun 2003): $${g}_{+}=\frac{1}{2}\left({c}_{1}+{c}_{2}+\sqrt{{\left({c}_{1}{c}_{2}\right)}^{2}+4{c}_{0}^{2}{c}_{1}{c}_{2}}\right)$$(30) $${g}_{}=\frac{1}{2}\left({c}_{1}+{c}_{2}\sqrt{{\left({c}_{1}{c}_{2}\right)}^{2}+4{c}_{0}^{2}{c}_{1}{c}_{2}}\right).$$(31)
The timescale of the pericenter precession τ_{pl} for each planet is thus given by 2π/g_{+} and 2π/g_{−} with the highest value associated with the more massive planet. This approximation is accurate under the condition presented in Takeda et al. (2008), namely $$\frac{\alpha}{13q\sqrt{\alpha}/{b}_{3/2}^{\left(1\right)}\left(\alpha \right)}\ll 1$$(32)
with q = m_{2}/m_{1}. In our work, the condition is not always fulfilled but even when it is not the case, the formula will still give us an indication about the order of magnitude of the precession timescale.
Initial parameters for the simulations.
4 Wide Binary with a_{B} = 1000 au
In this section we start by considering a binary companion at a_{B} = 1000 au. At such a separation, the timescale of the perturbations arising from the binary is longer than the disk’s lifetime, and the binary will have a minimal impact on the migrating planets. We use this setup to smoothly form planets in meanmotion resonances (MMRs) via disk migration. We then evolve the simulations for another 100 Myr, in order to study the longterm impact of the binary on the planets, with emphasis on the resonant dynamics. In the next section, closer binaries will be considered, which will have a visible impact on the planets as early as in the migration phase.
Regarding the initial parameters of the simulations, the masses and orbital elements of the two giant planets and the binary companion are given in Table 1. For the binary, eight different inclinations (with respect to the initial plane of the disk) ranging from nearly 0° to 70°, as well as four different eccentricities (10^{−3}, 0.1, 0.3, and 0.5) are considered here. The semimajor axes of the planets are chosen randomly between 5 and 7 au for the inner planet and between 12 and 21 au for the outer planet, in order to allow captures in different MMRs during the migration. The closest initial location of the two planets is beyond the 2:1 MMR ((12/7)^{3/2} ≃ 2.24), while the furthest one is beyond the 8:1 MMR ((21/5)^{3/2} ≃ 8.61). The planetary masses, eccentricities, and inclinations are chosen randomly as m ~ U[1; 5]m_{Jup}, e ~ U[0.001; 0.01], and i ~ U[0.01°; 0.1°]. For each combination of the fixed parameters of the binary companion, we drew 100 different random initial parameters for the planets.
4.1 Typical Examples
We start by presenting several typical evolutions observed in the simulations. For each case, we discuss both the disk phase and the longterm dynamical evolution after the dispersal of the disk. The evolutions are shown in the Laplace plane reference frame which is more suitable for the study of the ZLK resonance since the latter can easily be identified in this reference frame from the libration of the planetary pericenter argument (e.g., Libert & Tsiganis 2009). The Laplace plane is the invariant plane orthogonal to the total angular momentum of the system that almost corresponds, in our case, to the plane of the binary star. In this particular reference plane, the disk does not change in inclination but only precesses as shown in Roisin et al. (2021; see their Sect. 2.3.2).
A first example is shown in Fig. 1, for a binary companion with e_{B} = 10^{−3} and i_{B} = 50° and planetary masses of m_{1} = 1.44 M_{Jup} and m_{2} = 3.55 M_{Jup}. The evolution during the disk phase is shown in the left panels. We observe that the outer planet migrates inward and is captured at about 0.25 Myr in the 2:1 MMR with the inner planet, as confirmed by the libration of both resonant angles θ_{1} = 2λ_{2} − λ_{1} − ϖ_{1} and θ_{2} = 2λ_{2} − λ_{1}  ϖ_{2} (last panel), with λ = ω + Ω + M and ϖ = ω + Ω being the mean longitude and the longitude of the pericenter, respectively. After the capture in MMR, the two planets start migrating together and the eccentricity of the inner planet increases while the outer planet is kept on a quasi circular orbit (Lee & Peale 2002, see, e.g.). The eccentricity increase is limited due to the eccentricity damping by the disk and is followed by a slight increase of the mutual inclination between the planets which is rapidly damped by the disk. The disk keeps the planets within its midplane during the whole disk phase (see the mutual inclination between the disk and the planets in the sixth left panel), despite the influence of the inclined binary companion. Regarding the dynamical evolution of the system after the dispersal of the disk (right panels of Fig. 1), we observe that the MMR is preserved during the whole simulation. We also notice the linear evolutions of the planetary longitudes of the ascending nodes which, coupled with the constant inclinations, indicate the nodal precession of the planets due to the binary companion. Moreover, there is no evidence of a ZLK resonance with the highly inclined binary companion (i.e., no libration of the pericenter arguments around 90° or 270°). This is due to the strong gravitational interaction between the planets.
A second system is displayed in Fig. 2. Its evolution is similar to the one of the previous system during the disk phase, except now the planets get captured in the 3:1 MMR instead (left panels). However, after the disk phase (right panels), a scattering event rapidly occurs leading to the ejection of the outer planet. The remaining planet is captured in a ZLK resonance with the binary companion, as shown by the libration of the pericenter argument around 270° in the fourth right panel of Fig. 2 as well as the coupled variation in the planetary eccentricity and inclination.
In some cases, a destabilization of the system takes place during the disk phase, as illustrated in Fig. 3. At about 0.25 Myr, while the planets evolve in the 3:1 MMR and pursue their migration, a scattering event occurs. The outer planet is scattered to a larger separation and does not evolve in the disk anymore. The evolution of both planets rapidly becomes chaotic. They alternate between ZLK resonant states and nonresonant evolutions (see the evolution of the planetary arguments of pericenters), until the disruption of the system at about 6 Myr. The instability comes from the nearly equal strength of the ZLK effect induced by the binary and the gravitational interaction between the planets, as will be shown in Sect. 5.
Fig. 1 Typical evolution of a system of two giant planets with a wide binary at a_{B} = 1000 au (in the Laplace plane reference frame), during the disk phase (left panels) and after the disk dispersal (right panels). The initial parameters (with respect to the disk plane) are e_{B} = 10^{−3} and i_{B} = 50° for the binary companion and a_{1} = 6.7 au, m_{1} = 1.44 M_{Jup}, a_{2} = 18.9 au, and m_{2} = 3.55 M_{Jup} for the planets. 
Fig. 2 Same as Fig. 1, but for the initial parameters e_{B} = 0.5, i_{B} = 50°, a_{1} = 5.7 au, m_{1} = 2.77 M_{Jup}, a_{2} = 18.15 au, and m_{2} = 2.99 M_{Jup}. 
4.2 Results from a Large Ensemble of Simulations
We now study the system configurations found in the 3200 simulations performed with a_{B} = 1000 au. In order to interpret these simulations, it is useful to consider the relevant timescales. The ZKL and planetplanet interaction timescales are shown in Fig. 4, calculated with the system parameters found at the dispersal of the disk. As expected, for nearly all the systems, the ZLK timescale is much longer than the timescale associated with the gravitational interaction between the two planets, which explains that the influence of the binary is very limited in the simulations. In addition, the ZKL timescale is longer than the disk’s lifetime, indicating that the binary had little influence over the planets during the migration process. As a result, here we focus on the longterm dynamics after the disk phase, which a particular emphasis on the interplay between planetary MMRs and binary perturbations.
The distributions of the final orbital parameters are displayed in Fig. 5. To make it clearer, we present separately the systems for which the inclination of the binary companion is strictly below 40° (orange color) and the other systems for which the binary companion is highly inclined (i_{B} ≥ 40°) and which are thus likely to experience a ZLK resonance with the binary companion (blue color). Let us note that we also run the same 3200 simulations without the binary companion for comparison (yellow color). The distributions are shown at the dispersal of the disk in the left panels and at the end of the simulations in the right panels.
Fig. 3 Same as Fig. 1, but for the initial parameters e_{B} = 0.1, i_{B} = 50°, a_{1} = 6 au, m_{1} = 4.11 M_{Jup}, a_{2} = 12.74 au, and m_{2} = 3.47 M_{Jup}. 
4.2.1 MMR Captures
We first look at the period ratios between the planets (top panels of Fig. 5) to study the MMR captures during the disk phase. Two resonant commensurabilities are more prevalent, namely the 2:1 MMR followed by the 3:1 MMR, although several systems are also trapped close to the 5:2 MMR. Overall, the effect of the wide binary companion is small, as expected from the timescale analysis. We see that the systems with a binary companion (blue and orange) are captured slightly more often in the 2:1 MMR and less captured in the 3:1 MMR than the systems without a binary companion (yellow curve). No significant influence of the inclination of the binary companion on the resonance capture is visible. The left and right panels look very similar, which leads us to the conclusion that the systems captured in a MMR during the disk phase stay trapped for a long time (100 Myr in this case), even under the influence of the binary companion. This shows that for the wide binary case investigated here, the chaotic evolutions after the disk phase are rather limited in twoplanet resonant systems.
To confirm the previous results, we studied the evolution of the resonant angles for each system. In practice, the planets are considered locked in a MMR if at least one of the resonant angles librates with an amplitude of maximum 270° during the 0.2 Myr before the time considered, namely the dissipation of the disk or the end of the simulation. The results regarding the MMR captures are gathered in Tables 2 and 3. As previously observed, there are slightly more twoplanet systems captured in the 2:1 MMR for systems with a binary companion (53.6% against 49.7% at the dispersal of the disk, see Table 2) while it is the contrary for the 3:1 MMR (28.2% against 34.4%) and the 5:2 MMR (0.66% against 1.2%). Let us note that roughly the same percentages of planetary ejections at the dispersal of the disk are observed for single star systems and binary star systems (see the two last rows of Table 2). In Table 3, we see that the percentages of final oneplanet systems have increased, in particular for binary star systems due to chaotic events occurring after the dispersal of the disk. Examples of chaotic evolutions after the disk phase were shown in Figs. 2 and 3. A study of the stability of the twoplanet systems formed here is done in Sect. 6.
Fig. 4 Timescale of the pericenter precession due to the gravitational interaction between the planets versus the ZLK timescale due to the binary companion. The green and red colors refer to the inner planet and outer planet, respectively. 
Fig. 5 Normalized distributions of final orbital elements at the dispersal of the disk (left panel) and at the end of the simulation (right panel). From top to bottom: period ratio, eccentricity, mutual inclination, and argument of the pericenter (in the Laplace plane reference frame). The color code refers to the presence of a binary companion and its inclination value for the systems with a companion. 
Percentages of the different system configurations at the dissipation of the disk.
4.2.2 Eccentricities
The eccentricity distribution of the different systems is shown in the panels of the second row of Fig. 5. We observe that the single star systems with two planets have slightly higher eccentricities (orange and blue curves above the yellow curve for eccentricities below 0.25 while the yellow curve is above the two others for higher values). This observation is coherent with the percentages of resonance captures. The proportion of twoplanet systems which are found in MMR is indeed more important for single star systems (see Tables 2 and 3) and those resonance captures are generally followed by an increase of the planetary eccentricities. Again, the similar profile between the curves in the left panel (at the dispersal of the disk) and the right panel (at the end of the simulation) leads us to say that the wide binary companion has a limited influence on the system even after the dispersal of the disk.
Percentages of the different system configurations at the end of the simulation.
4.2.3 Inclinations
In the third row panels of Fig. 5, we show the distribution of the mutual inclination between the planets. We observe mutual inclination values slightly higher for systems including a binary companion. However, most of the mutual inclinations are below 5°. Only 3% of the systems have a mutual inclination above 5° in binary star systems. These higher mutual inclinations are generally produced by an inclinationtype resonance. This mechanism is known to operate at large eccentricities in twoplanet systems (Thommes & Lissauer 2003; Libert & Tsiganis 2009; Teyssandier & Terquem 2014; Sotiriadis & Libert 2020), but can also be present at small to moderate eccentricities in systems with more than two planets (Libert et al. 2018). An example of a highly mutually inclined system is shown in Fig. 6. We first observe a temporary capture in the 3:1 MMR, before a second capture in the 5:2 MMR. The eccentricities increase to moderate values and the inclinations suddenly grow at 9 × 10^{5} yr when the system enters an inclinationtype resonance. To better identify the resonant process leading to the inclination increase, we display at the bottom panel of Fig. 6 the evolution of the angle ${\theta}_{{i}_{1}^{2}}=5{\lambda}_{2}2{\lambda}_{1}{\varpi}_{1}2{\text{\Omega}}_{1}$ with the evolution of the mutual inclination (in logarithmic scale). We see that the inclination increases in correlation with the libration of the inclinationtype resonant angle ${\theta}_{{i}_{1}^{2}}$.
The planets mostly remain in the midplane of the disk during the disk phase thanks to the gravitational and damping forces of the disk. We note that the planets tend to be slightly more misaligned with respect to the disk in systems with a binary companion but this misalignment is not significant (less than 4°).
In Fig. 7, we display the planetary inclinations with respect to the initial plane of the disk. We observe in the left panel showing the inclinations at the dispersal of the disk, that for binary star systems the planetary inclinations are proportional to the inclination of the binary companion. This is coherent with the precession of the disk due to the binary companion, which induces a periodic variation with an amplitude equal to 2i_{B}. Since the planets tend to stay coupled with the disk, they follow its nodal precession, which explains the inclination distribution in the left panel of Fig. 7. On the right panel, we observe that the inclination is amplified during the dynamical phase. This is a consequence of the nodal precession induced by the binary companion acting this time on the coupled planets. Moreover, on the right panel, we observe that the normalized frequencies decrease with the inclination. This effect comes from the amplitude of the nodal precession. Indeed, we have 400 systems for each inclination of the binary, whose planetary inclinations are uniformly spread between 0 and 2i_{B}, leading to an accumulation around the smallest values.
Assuming that the disk angular momentum vector is originally aligned with the stellar spin and that the stellar spin does not change direction over time (although see Storch et al. 2014, for conditions where this assumption might break down), and observing that the planets mostly remain in the disk’s midplane throughout the migration, we can use the disk’s inclination as a proxy for measuring spinorbit angles. The nodal precession induced by the binary companion thus creates a misalignment between the spin axis of the primary star and the angular momentum of the planets.
Fig. 6 Typical evolution of two giant planets (in a wide binary star system) experiencing an inclinationtype resonance during the disk phase. The initial parameters are e_{B} = 0.5, i_{B} = 0.001°, a_{1} = 6.86 au, m_{1} = 1.64 M_{Jup}, a_{2} = 17.1 au, and m_{2} = 3.79 M_{Jup}. 
Fig. 7 Inclination distribution at the dissipation of the disk (left panel) and at the end of the simulation (right panel). The inclinations are given in the initial disk plane reference frame. 
4.2.4 ZLK Resonance
To further analyze the influence of the binary companion, in particular the ZLK mechanism, we now focus on the pericenter argument of the planets. The bottom panels of Fig. 5 display the pericenter argument distributions at the dissipation of the disk (left panel) and at the end of the simulation (right panel). We do not notice any specific pattern associated with the ZLK resonance (i.e., accumulations around 90 or 270°). Regarding the evolution after the disk phase (right panel), the gravitational interaction between the two planets tends to overcome the effect of the binary companion and the establishment of a ZLK resonant evolution in highly inclined binary systems. Note that we also numerically investigated the libration of the pericenter arguments and observed that only the systems suffering from an ejection (i.e., systems with one remaining planet) possibly hold a planet locked in the ZLK resonance.
Overall, our simulations indicate that for planetary systems formed in very wide binaries, resonant captures still occur frequently and are robust against perturbations from the binary. In particular, on the longterm evolution, resonant interactions seem to quench ZLK oscillations.
5 Closer Binary with a_{B} = 250 au
As seen in the previous section, in case of a wide binary companion at 1000 au, the influence of the binary companion on planet pairs embedded in the gas disk is very limited and the resonant evolution of the pairs is preserved well after the dispersal of the disk. In this section, we investigate the case of closer binary stars.
5.1 Second Set of Simulations
Based on the observations related to Fig. 4, we aimed to conduct a parametric study to identify parameter values leading to an increase of the ratio between the precession timescale pl and the ZLK resonance timescale τ_{LK}. Following Eq. (27), to reduce the ZLK timescale, we decreased the semimajor axis of the binary companion by adopting the values 250, 500, and 750 au. Note that when the binary companion is closer to the primary star, our disk model suffers from some limitations that will be discussed in detail in Sect. 5.3. We also increased the initial semimajor axes of the planets by multiplying their range of values by 1.5 in order to increase their orbital periods. A last change consisted in decreasing the lower bound on the planetary masses to 0.65 M_{Jup}. All the other parameters were kept unchanged, as given in Table 1. For each of the four aB values and for the 32 choices of the binary eccentricity and inclination, we randomly drew 50 different initial conditions for the planetary parameters. In total, we ran 8000 simulations (single star simulations included).
We now study the system configurations resulting from this second set of simulations at the dispersal of the disk. In Fig. 8, we reproduce the same plot as in Fig. 4 for the different semimajor axes of the binary companion considered here. As expected from Eqs. (27), (30), and (31), we see that the smaller the semimajor axis of the binary, the closer the systems to the dashed line of equal timescales. As a result, the impact of the binary companion will be stronger for simulations with a_{B} = 250 au and we therefore chose to focus on this binary separation in the remainder of this section.
Fig. 10 Eccentricity and pericenter argument distributions at the end of the simulations for single planet systems (left panels) and twoplanet systems (right panels) when a_{B} = 250 au. 
5.2 Results
The results are displayed in Fig. 9. We observe that systems with a binary companion (cyan color) are more eccentric than the ones without companion (yellow color), especially at the end of the simulations where highly eccentric planetary orbits are observed. They mainly result from scattering events and ejections during the longterm evolution. As shown in Fig. 10, for the simulations with a binary companion ending with two planets, almost all the eccentricities are below 0.2 (right panel), while most of the oneplanet systems present very eccentric configurations (left panel). Note that in Fig. 10 we also see that the inclination of the binary companion influences the planetary eccentricities, since we observe that the majority of the highly eccentric systems have a binary companion with an inclination higher than 40° (i.e., close to the critical value of 39.23° for the ZLK mechanism.).
Again, the mutual inclinations between the planets are small, mainly below 4° (panels of the third row of Fig. 9). Moreover, the distributions of the pericenter arguments shown in the bottom panels of Fig. 9 are nearly uniform, but we note this time a possible evidence of the ZLK resonance. Indeed, if we present separately the single planet systems and the twoplanet systems formed at the end of the simulations, we see in Fig. 10 that the distribution of the systems ending with one planet only presents two peaks centered at 90° and 270° associated with the ZLK resonance (bottom left panel). A study of the libration of the pericenter arguments for each system individually leads us to the observation that no system with two planets end up with one of the planets in the ZLK resonance with the binary companion. For the oneplanet systems, the percentage of systems locked in the ZLK resonance with an inclined binary companion is ~30%, similarly to the observation made in Roisin et al. (2021). However, even if the system is not evolving in the ZLK resonance, its dynamics is still influenced by the ZLK islands which can explain the high planetary eccentricities observed here. For a dynamical study of the oneplanet case, we refer to Roisin & Libert (2021); Roisin et al. (2021).
Contrary to the case where a_{B} = 1000 au, the inclination of the planets with respect to the disk was modified during the disk phase for a_{B} = 250 au, as a consequence of the much shorter timescale of the binary perturbations. This can be observed in the typical evolution displayed in Fig. 11. In Fig. 12 we show the inclination of the planets with respect to the disk at the time of the dispersal of the disk. For binaries with inclinations less than 40° (red color), about half of the planets remain in the disk’s midplane, while the other half show inclinations ranging up to 90°, with a peak around 40°. This distribution could be in part caused by the nodal precession which induces inclinations ranging from 0 to 2i_{b}. Indeed, some of these planets may have been separated from the disk’s midplane due to the differential precession induced by the binary, which occurs on a timescale short enough that the damping forces and diskinduced precession cannot counteract it. Note that some planets may also have been expelled from the disk’s midplane due to inclinationtype resonances, as illustrated in Fig. 6. For binaries with inclinations larger than 40° (blue color), the distribution of the mutual diskplanet inclinations broadens and an excess is observed at inclinations near 30° and 140°. After being ejected from the disk’s midplane, these planets may be undergoing ZLK cycles driven by the disk, as shown by Terquem & Ajmia (2010), for which the critical angle can be lower than 40° (Teyssandier et al. 2013).
Regarding the MMRs, a variation in the initial semimajor axes of the planets leads to a modification of their migration rates and thus possibly to different MMR captures. In the new set of simulations with smaller a_{B} values, almost all the resonant systems are trapped in the 2:1 MMR. During the disk phase, the evolutions of the systems are very similar to the ones with a wide binary companion described in Sect. 4. At the end of the simulation, a lot of ejections are reported for the simulations with a_{B} = 250 au, as can be observed in Table 4 which shows the system configurations at the end of the simulation for the four considered binary separations for comparison. In Fig. 13, we observe that, when a_{B} = 250 au, the systems ending with a single planet at the end of the simulation (purple curve) have a mean ratio between τ_{pl} and τ_{LK} higher than the one for systems with two planets (green curve). This confirms that systems for which both timescales are roughly similar experience more planetary ejections. For all the other values of the binary semimajor axis, Table 4 shows that the influence of the binary companion on the system architecture is still quite limited.
In conclusion, we observed that a binary companion close enough to strongly influence the planets (a_{B} = 250 au) will tend to destabilize the systems and lead to the ejection of a planet in about one third of the cases. The ejection will possibly induce high eccentricity and ZLK resonance capture for the remaining planet.
Fig. 11 Typical evolution of a system of two giant planets with a binary at a_{B} = 250 au (in the Laplace plane reference frame), during the disk phase. The initial parameters (with respect to the disk plane) are e_{B} = 0.5 and i_{B} = 50° for the binary companion and a_{1} = 10.22 au, m_{1} = 1.39 M_{Jup}, a_{2} = 25.33 au, and m_{2} = 1.07 M_{Jup} for the planets. 
Fig. 12 Histogram of the inclination of the planets with respect to the disk, at the dissipation of the disk for simulations with a binary separation a_{B} = 250 au. Blue (resp. orange) bins are for binary orbital inclinations above (resp. below) 40°. 
Percentages of the different system configurations at the end of the simulation for the second set of simulations.
Fig. 13 Normalized frequency of the logarithm of the ratio between the timescale of the gravitational interaction between the planets and the ZLK timescale, for the simulations with a_{B} = 250 au. The color code refers to the number of planets in the system at the end of the simulation. 
Fig. 14 Time evolution of the inclination of the disk as a function of the semimajoraxis, for different semimajoraxes of the binary (in the initial disk plane reference frame). The line type refers to different times: the solid line stands for t = 0.5 Myr, the dashed line for t = 0.75 Myr and the dotted line for t = 1.123 Myr. The inclination of the binary i_{B} is fixed to 60°. 
5.3 Caveats of the Simulations with a_{B} = 250 au
In this section we discuss the possible limitations of our disk model for the parameters of the second set of simulations. Indeed, we have assumed here that the disk is uniform due to its selfgravity. This hypothesis is valid for wide binary stars (e.g., Roisin et al. 2021). Closer binary companions as the ones considered in the second part of our work could have the effect to overcome the disk selfgravity and lead to deformations of the disk. By splitting the disk in several massive rings adjacent to one another and making them evolve using the LaplaceLagrange theory to estimate the nodal precession of the disk, we studied in Fig. 14 the uniformity of the disk for the different semimajor axes of the binary companion considered here. More details about this technique can be found in Murray & Dermott (1999); Levison & Morbidelli (2007); Batygin et al. (2011); Roisin et al. (2021). Note that we fix the inclination of the binary companion to i_{B} = 60° arbitrarily, but the results are similar regardless of the chosen inclination. We see that the disk remains uniform for a binary companion at 750 and 1000 au. On the contrary, the disk uniformity is not maintained for 250 and 500 au. This does not mean that the disk is not uniform in this case, but the selfgravity is not sufficient alone to ensure it. Hydrodynamical effects such as radial pressure force and viscous diffusion (e.g., Papaloizou & Terquem1995; Larwood et al. 1996) could play a role and keep the disk uniform. For instance, it has been shown that if the wave crossing time in the disk is shorter than the precession time due to the binary companion, the rigidity of the disk could be maintained (Zanazzi & Lai 2018a). Investigations of such effects is left for future work.
Moreover, the disk could experience ZLK cycles due to the presence of the binary companion (see, e.g., Martin et al. 2014; Fu et al. 2015b; Zanazzi & Lai 2017). Batygin et al. (2011) showed that the cycles could be suppressed by the selfgravity of the disk if the disk is massive enough. This assumption will eventually break down in our simulations since we implemented an exponential decrease of the disk mass. Inaddition, ZLK cycles in the disk could lead to eccentric orbits and the hypothesis of the LaplaceLagrange theory requiring annuli without intersections could not be verified in this case.
Nevertheless, most of the results of the second set of simulations were very similarforall the semimajoraxes considered for the binary companion. The differences observed for a_{B} = 250 au mainly arise during the longterm evolution after the dispersal of the disk. Since this phase does not suffer from the limitations in the disk model presented in this section, we therefore conclude that despite those limitations, our results are quite robust.
6 Stability
To assess the stability of the twoplanet systems formed in our simulations, we used the Mean Exponential Growth factor of Nearby Orbits (MEGNO) chaos indicator which is based on the evolution of one deviation vector (Cincotta et al. 2003). This indicator differentiates efficiently the orbits since it converges to 0 for stable periodic orbits, to 2 for quasiperiodic orbits and orbits close to stable periodic orbits, and diverges with time for irregular orbits. In practice, the simulations were made with the symplectic WisdomHolman integrator WHFast implemented in REBOUND (Rein & Tamayo 2015) and we started to follow the evolution of the MEGNO indicator just after the disk phase.
A visual representation of the results is given in Fig. 15 which shows for each system the semimajor axes of the two planets. The twoplanet systems with a binary companion at 250 au are shown in the left panel and the ones when a_{B} = 1000 au are shown in the right panel. The green dashed lines indicate the different MMRs. The high density of points close to the lines confirms that the systems gather around the 2:1 MMR for closer binaries, while captures in higher order resonances like the 3:1 and 5:2 MMRs are possible forwiderbinaries. Moreover, as previously noted, fewer twoplanet systems are formed when the binary companion is at 250 au due to the largest proportion of planetary ejections.
The color code of Fig. 15 refers to the MEGNO value. Note that the color bar upper value is fixed to 10 since we stopped the simulation when the MEGNO reached this value. We observe that in the wide binary case, many systems close to the 2:1 and 3:1 MMRs have a regular motion while the evolution of most of the systems near the 5:2 MMR is irregular. The exact percentages of systems with regular (defined here as a MEGNO value smaller than 2.5) and irregular (MEGNO value higher than 2.5) orbits is given in Table 5. On the contrary, for close binaries, only 22% of the systems near the 2:1 MMR present a regular evolution, showing the significant impact of a close binary companion on the longterm evolution of planetary systems. As expected, we observe in Table 5 that the percentage of regular systems is slightly lower for highly inclined binary stars (above 40°).
Fig. 15 Semimajor axis of the outer planet with respect to the one of the inner planet for the twoplanet systems formed in our simulations with the binary at a_{B} = 250 au (left panel) and at a_{B} = 1000 au (right panel). The color code refers to the MEGNO value. 
Percentages of regular and irregular evolutions as given by the MEGNO indicator.
7 Comparison with Observations
In this section we compare our results with the 187 Stype planets in binary star systems presented in the Open Exoplanet Catalogue of Rein (2012). In particular, the comparison of the cumulative eccentricities is made in Fig. 16. The eccentricity distribution of Rein (2012) for Stype planets is indicated by the green curve, while the one of all the detected planets is indicated by the black curve. Note that we only considered Jupiterlike planets with a mass higher than 0.5 MJup. As noticed in Kaib et al. (2013), planets in Stype systems have larger eccentricities. We observe that for the systems with a binary companion at 1000 au (blue curve) and without binary companion (yellow curve) found in our simulations, nearly all the eccentricities are below ~0.2, since the twoplanet systems are weakly perturbed during their longterm evolution. On the contrary, high eccentricities are found in simulations with a close binary companion at 250 au (orange curve) and can notably be explained by the influence of the ZLK resonance on the planetary evolution after the ejection of one of the planets (see Roisin & Libert 2021, for more details). We also added in Fig. 16 the systems found in Roisin et al. (2021) for simulations of single planets migrating in the disk with a binary companion at 1000 au (purple curve). While binary companions with i_{B} < 40° do not excite the initial quasicircular orbits of the planets, high eccentricities can be reached in the presence of a highly inclined binary companion (cyan curve). As a result, we showed that the high eccentricities observed for Stype planets in Fig. 16 can be caused by highly inclined wide binary companions in case of single planet systems and close binary companions in case of multiplanet systems.
Moreover, several observations of misalignment between the stellar spin axis and the orbital angular momentum of hot Jupiters have been reported (Wright et al. 2011; Albrecht et al. 2012). A common scenario to explain this misalignment is the action of tidal friction during the ZLK cycles due to an highly inclined binary companion (e.g., Naoz & Fabrycky 2014). In this work we observed misalignment induced by the nodal precession of the planets due to their interaction with the binary companion, regardless of the inclination of the binary companion. This effect was extensively studied during the disk phase (see, e.g., Batygin & Adams 2013; Lai 2014; Spalding & Batygin 2014; Zanazzi & Lai 2018b) but it continues even after the dissipation of the disk (Boué & Fabrycky 2014). The nodal precession effect could lead to possibly strong misalignment even for planets further away of the star than the hot Jupiters and could help explaining the observed misalignment, alongside other mechanisms. Note that in this work we considered the spin axis of the star as fixed and parallel to the disk angular momentum at the beginning of the simulation. Of course, it is a strong assumption since the disk would most likely precess before the formation of the giant planets and would be misaligned. Moreover, the spin axis of the host star could also be chaotic and not constant, as shown by Storch et al. (2014) for a primary star under the influence of a hot Jupiter for example.
Fig. 16 Cumulative eccentricity of the planets found in our simulations and in the observations (Rein 2012). See text for more details. 
8 Conclusions
In this paper, we studied the impact of a binary companion on the migration of two giant planets in a protoplanetary disk, followed by the longterm evolution of the systems after the disk vanished. The eccentricity and inclination damping induced by the disk, its gravitational potential, and the nodal precession induced by the binary companion on the disk were all considered here. A timescale analysis hinted at the existence of two main patterns in the disk migration phase. First, for very wide binaries (a_{B} = 1000 au), the nodal precession induced by the binary occurs on a timescale longer than the disk’s lifetime. In this case, damping forces from the disk cause the planets to remain close to its midplane while they migrate and the entire diskplanets system precesses rigidly very slowly around the binary. Resonant pairs are routinely formed in these simulations, mostly at the 2:1, 5:2, and 3:1 commensurabilities. In addition, the strong gravitational interaction between the resonant planets acts to reduce the impact of the binary companion after the dissipation of the disk, preserving their resonant structure, and shielding them from undergoing ZLK oscillations. Second, for closer binaries (a_{B} = 250 au), the nodal precession timescale induced by the binary is in the range 10^{5}−10^{6} yr, which is similar to the disk’s lifetime. In this case, the damping forces and gravitational potential from the disk are not always able to maintain the planets within its midplane, and some primordial diskplanet inclinations can be generated, especially for highly inclined binary companions. At the end of these simulations, the fraction of systems in MMRs has dropped significantly compared to simulations with a_{B} = 1000 au, while the fraction of unstable systems with ejections has increased. Those ejections usually leave the remaining planets with high eccentricity. This could give an additional explanation for the existence of many Stype planets detected with high eccentricity values.
Our study shows that the existence of a binary companion can lead to significant misalignment between the planets and their natal protoplanetary disk. This suggests that misalignment between planetary orbital angular momentum and stellar spin could be generated early on in the system’s lifetime, while the disk is still present. A natural continuation of this work would be to consider binary companions whose parameter distributions closely match observations, to increase the number of planets as is done in planetplanet numerical experiments, and finally to include tides and spin dynamics to properly assess the distribution of spinorbit misalignments.
Acknowledgements
The authors would like to warmly thank the referee for providing insightful comments which improved this work. We also would like to warmly thank A. Crida, A. Lemaître, S. Raymond and K. Tsiganis for useful discussions. The work of A. Roisin was supported by a F.R.S.FNRS Research Fellowship and the work of J. Teyssandier was supported by a F.R.S.FNRS Postdoctoral Research Fellowship. This work was also supported by the Fonds de la Recherche Scientifique – FNRS under Grant No. F.4523.20 (DYNAMITE MISproject). Computational resources have been provided by the Consortium des Equipements de Calcul Intensif, supported by the FNRSFRFC, the Walloon Region, and the University of Namur (Conventions No. 2.5020.11, GEQ U.G006.15, 1610468 et RW/GEQ2016).
References
 Albrecht, S., Winn, J. N., Johnson, J. A., et al. 2012, ApJ, 757, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Artymowicz, P., & Lubow, S. H. 1994, ApJ, 421, 651 [Google Scholar]
 Batygin, K., & Adams, F. C. 2013, ApJ, 778, 169 [NASA ADS] [CrossRef] [Google Scholar]
 Batygin, K., Morbidelli, A., & Tsiganis, K. 2011, A&A, 533, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bitsch, B., Crida, A., Libert, A.S., & Lega, E. 2013, A&A, 555, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Boué, G., & Fabrycky, D. C. 2014, ApJ, 789, 111 [CrossRef] [Google Scholar]
 Brouwer, D., & Clemence, G. M. 1961, Methods of Celestial Mechanics (New York: Academic Press) [Google Scholar]
 Chachan, Y., Booth, R. A., Triaud, A. H. M. J., & Clarke, C. 2019, MNRAS, 489, 3896 [NASA ADS] [Google Scholar]
 Chambers, J. E., Quintana, E. V., Duncan, M. J., & Lissauer, J. J. 2002, AJ, 123, 2884 [NASA ADS] [CrossRef] [Google Scholar]
 Cincotta, P. M., Giordano, C. M., & Simó, C. 2003, Phys. D Nonlinear Phenomena, 182, 151 [NASA ADS] [CrossRef] [Google Scholar]
 Crida, A., & Morbidelli, A. 2007, MNRAS, 377, 1324 [Google Scholar]
 Duncan, M. J., Levison, H. F., & Lee, M. H. 1998, AJ, 116, 2067 [Google Scholar]
 Duquennoy, A., & Mayor, M. 1991, A&A, 248, 485 [NASA ADS] [Google Scholar]
 Fabrycky, D., & Tremaine, S. 2007, ApJ, 669, 1298 [NASA ADS] [CrossRef] [Google Scholar]
 Franchini, A., Martin, R. G., & Lubow, S. H. 2019, MNRAS, 485, 315 [NASA ADS] [CrossRef] [Google Scholar]
 Fu, W., Lubow, S. H., & Martin, R. G. 2015a, ApJ, 807, 75 [NASA ADS] [CrossRef] [Google Scholar]
 Fu, W., Lubow, S. H., & Martin, R. G. 2015b, ApJ, 813, 105 [NASA ADS] [CrossRef] [Google Scholar]
 Huré, J.M., & Hersant, F. 2011, A&A, 531, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ivanov, P. B., Papaloizou, J. C. B., & Polnarev, A. G. 1999, MNRAS, 307, 79 [Google Scholar]
 Kaib, N. A., Raymond, S. N., & Duncan, M. 2013, Nature, 493, 381 [Google Scholar]
 Kiseleva, L. G., Eggleton, P. P., & Mikkola, S. 1998, MNRAS, 300, 292 [NASA ADS] [CrossRef] [Google Scholar]
 Kozai, Y. 1962, AJ, 67, 591 [Google Scholar]
 Kraus, A. L., & Ireland, M. J. 2012, ApJ, 745, 5 [Google Scholar]
 Lai, D. 2014, MNRAS, 440, 3532 [NASA ADS] [CrossRef] [Google Scholar]
 Larwood, J. D., Nelson, R. P., Papaloizou, J. C. B., & Terquem, C. 1996, MNRAS, 282, 597 [NASA ADS] [Google Scholar]
 Lee, M. H., & Peale, S. J. 2002, ApJ, 567, 596 [Google Scholar]
 Levison, H., & Morbidelli, A. 2007, Icarus, 189, 196 [NASA ADS] [CrossRef] [Google Scholar]
 Libert, A.S., & Tsiganis, K. 2009, MNRAS, 400, 1373 [NASA ADS] [CrossRef] [Google Scholar]
 Libert, A.S., Sotiriadis, S., & Antoniadou, K. I. 2018, Celest. Mech. Dyn. Astron., 130, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Lidov, M. L. 1962, Planet. Space Sci., 9, 719 [Google Scholar]
 Lubow, S. H., & Martin, R. G. 2016, ApJ, 817, 30 [NASA ADS] [CrossRef] [Google Scholar]
 Lubow, S. H., & Ogilvie, G. I. 2017, MNRAS, 469, 4292 [NASA ADS] [CrossRef] [Google Scholar]
 Martin, R. G., Nixon, C., Lubow, S. H., et al. 2014, ApJ, 792, L33 [Google Scholar]
 Martin, R. G., Lubow, S. H., Nixon, C., & Armitage, P. J. 2016, MNRAS, 458, 4345Matsumoto, Y., Nagasawa, M., & Ida, S. 2012, Icarus, 221, 624 [Google Scholar]
 Murray, C. D., & Dermott, S. F. 1999, Solar System Dynamics (Cambridge: Cambridge University Press) [Google Scholar]
 Naoz, S., & Fabrycky, D. C. 2014, ApJ, 793, 137 [Google Scholar]
 Nelson, R. P., Papaloizou, J. C. B., Masset, F., & Kley, W. 2000, MNRAS, 318, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Papaloizou, J. C. B., & Larwood, J. D. 2000, MNRAS, 315, 823 [Google Scholar]
 Papaloizou, J. C. B., & Terquem, C. 1995, MNRAS, 274, 987 [NASA ADS] [Google Scholar]
 Picogna, G., & Marzari, F. 2015, A&A, 583, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rafikov, R. R., & Silsbee, K. 2015, ApJ, 798, 70 [Google Scholar]
 Raghavan, D., McAlister, H. A., Henry, T. J., et al. 2010, ApJS, 190, 1 [Google Scholar]
 Rein, H. 2012, ArXiv eprints [arXiv: 1211.7121] [Google Scholar]
 Rein, H., & Tamayo, D. 2015, MNRAS, 452, 376 [Google Scholar]
 Roisin, A., & Libert, A. S. 2021, A&A, 645, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Roisin, A., Teyssandier, J., & Libert, A.S. 2021, MNRAS, 506, 5005 [NASA ADS] [CrossRef] [Google Scholar]
 Savonije, G. J., Papaloizou, J. C. B., & Lin, D. N. C. 1994, MNRAS, 268, 13 [NASA ADS] [CrossRef] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
 Sotiriadis, S., & Libert, A.S. 2020, MNRAS, 495, 192 [NASA ADS] [Google Scholar]
 Sotiriadis, S., Libert, A.S., Bitsch, B., & Crida, A. 2017, A&A, 598, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Spalding, C., & Batygin, K. 2014, ApJ, 790, 42 [NASA ADS] [CrossRef] [Google Scholar]
 Storch, N. I., Anderson, K. R., & Lai, D. 2014, Science, 345, 1317 [NASA ADS] [CrossRef] [Google Scholar]
 Takeda, G., Kita, R., & Rasio, F. A. 2008, ApJ, 683, 1063 [NASA ADS] [CrossRef] [Google Scholar]
 Terquem, C., & Ajmia, A. 2010, MNRAS, 404, 409 [NASA ADS] [Google Scholar]
 Terquem, C., & Bertout, C. 1993, A&A, 274, 291 [NASA ADS] [Google Scholar]
 Teyssandier, J., & Terquem, C. 2014, MNRAS, 443, 568 [Google Scholar]
 Teyssandier, J., Terquem, C., & Papaloizou, J. C. B. 2013, MNRAS, 428, 658 [NASA ADS] [CrossRef] [Google Scholar]
 Thommes, E. W., & Lissauer, J. J. 2003, ApJ, 597, 566 [Google Scholar]
 von Zeipel, H. 1910, Astron. Nachr., 183, 345 [Google Scholar]
 Wright, J. T., Fakhouri, O., Marcy, G. W., et al. 2011, PASP, 123, 412 [NASA ADS] [CrossRef] [Google Scholar]
 Wu, Y., & Murray, N. 2003, ApJ, 589, 605 [Google Scholar]
 XiangGruess, M., & Papaloizou, J. C. B. 2014, MNRAS, 440, 1179 [NASA ADS] [CrossRef] [Google Scholar]
 Zanazzi, J. J., & Lai, D. 2017, MNRAS, 467, 1957 [Google Scholar]
 Zanazzi, J. J., & Lai, D. 2018a, MNRAS, 477, 5207 [NASA ADS] [CrossRef] [Google Scholar]
 Zanazzi, J. J., & Lai, D. 2018b, MNRAS, 478, 835 [NASA ADS] [CrossRef] [Google Scholar]
 Zhou, J.L., & Sun, Y.S. 2003, ApJ, 598, 1290 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Percentages of the different system configurations at the dissipation of the disk.
Percentages of the different system configurations at the end of the simulation for the second set of simulations.
All Figures
Fig. 1 Typical evolution of a system of two giant planets with a wide binary at a_{B} = 1000 au (in the Laplace plane reference frame), during the disk phase (left panels) and after the disk dispersal (right panels). The initial parameters (with respect to the disk plane) are e_{B} = 10^{−3} and i_{B} = 50° for the binary companion and a_{1} = 6.7 au, m_{1} = 1.44 M_{Jup}, a_{2} = 18.9 au, and m_{2} = 3.55 M_{Jup} for the planets. 

In the text 
Fig. 2 Same as Fig. 1, but for the initial parameters e_{B} = 0.5, i_{B} = 50°, a_{1} = 5.7 au, m_{1} = 2.77 M_{Jup}, a_{2} = 18.15 au, and m_{2} = 2.99 M_{Jup}. 

In the text 
Fig. 3 Same as Fig. 1, but for the initial parameters e_{B} = 0.1, i_{B} = 50°, a_{1} = 6 au, m_{1} = 4.11 M_{Jup}, a_{2} = 12.74 au, and m_{2} = 3.47 M_{Jup}. 

In the text 
Fig. 4 Timescale of the pericenter precession due to the gravitational interaction between the planets versus the ZLK timescale due to the binary companion. The green and red colors refer to the inner planet and outer planet, respectively. 

In the text 
Fig. 5 Normalized distributions of final orbital elements at the dispersal of the disk (left panel) and at the end of the simulation (right panel). From top to bottom: period ratio, eccentricity, mutual inclination, and argument of the pericenter (in the Laplace plane reference frame). The color code refers to the presence of a binary companion and its inclination value for the systems with a companion. 

In the text 
Fig. 6 Typical evolution of two giant planets (in a wide binary star system) experiencing an inclinationtype resonance during the disk phase. The initial parameters are e_{B} = 0.5, i_{B} = 0.001°, a_{1} = 6.86 au, m_{1} = 1.64 M_{Jup}, a_{2} = 17.1 au, and m_{2} = 3.79 M_{Jup}. 

In the text 
Fig. 7 Inclination distribution at the dissipation of the disk (left panel) and at the end of the simulation (right panel). The inclinations are given in the initial disk plane reference frame. 

In the text 
Fig. 8 Same as Fig. 4, but for the different semimajor axes of the binary companion. 

In the text 
Fig. 9 Same as Fig. 5, but for a_{B} = 250 au. 

In the text 
Fig. 10 Eccentricity and pericenter argument distributions at the end of the simulations for single planet systems (left panels) and twoplanet systems (right panels) when a_{B} = 250 au. 

In the text 
Fig. 11 Typical evolution of a system of two giant planets with a binary at a_{B} = 250 au (in the Laplace plane reference frame), during the disk phase. The initial parameters (with respect to the disk plane) are e_{B} = 0.5 and i_{B} = 50° for the binary companion and a_{1} = 10.22 au, m_{1} = 1.39 M_{Jup}, a_{2} = 25.33 au, and m_{2} = 1.07 M_{Jup} for the planets. 

In the text 
Fig. 12 Histogram of the inclination of the planets with respect to the disk, at the dissipation of the disk for simulations with a binary separation a_{B} = 250 au. Blue (resp. orange) bins are for binary orbital inclinations above (resp. below) 40°. 

In the text 
Fig. 13 Normalized frequency of the logarithm of the ratio between the timescale of the gravitational interaction between the planets and the ZLK timescale, for the simulations with a_{B} = 250 au. The color code refers to the number of planets in the system at the end of the simulation. 

In the text 
Fig. 14 Time evolution of the inclination of the disk as a function of the semimajoraxis, for different semimajoraxes of the binary (in the initial disk plane reference frame). The line type refers to different times: the solid line stands for t = 0.5 Myr, the dashed line for t = 0.75 Myr and the dotted line for t = 1.123 Myr. The inclination of the binary i_{B} is fixed to 60°. 

In the text 
Fig. 15 Semimajor axis of the outer planet with respect to the one of the inner planet for the twoplanet systems formed in our simulations with the binary at a_{B} = 250 au (left panel) and at a_{B} = 1000 au (right panel). The color code refers to the MEGNO value. 

In the text 
Fig. 16 Cumulative eccentricity of the planets found in our simulations and in the observations (Rein 2012). See text for more details. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.