© A. Roisin et al. 2022

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 Sun-like stars are part of multiple star systems (Duquennoy & Mayor 1991; Raghavan et al. 2010). Until now, more than 100 S-type 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 planet-disk 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., Xiang-Gruess & 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 Zeipel-Lidov-Kozai (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 self-gravity 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 S-type 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 Type-II 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 long-term 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 set-up of the N-body 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 long-term 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 S-type 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 semi-major 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 well-known N-body 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 type-II 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 Type-II regime given by (both in the disk-dominated case and planet-dominated cases, see Sotiriadis et al. 2017, for more details): $${\tau }_{\text{mig}}=\frac{2}{3}{\alpha }^{-1}{h}^{-2}{\text{Ω}}_{\text{pl}}^{-1}\times \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 Shakura-Sunyanev viscosity parameter (Shakura & Sunyaev 1973), h = 0.05 the disk aspect ratio, ${\text{Ω}}_{\text{pl}}^{-1}$ the orbital frequency of the planet, a_pl the semi-major axis of the planet, m_pl the mass of the planet, r_pl the distance of the planet to the host star, and Σ ∝ r⁻⁰⁵ 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 gas-free inner cavity using an hyperbolic tangent function $\tanh \left(\frac{r-R_{\text{in}}}{\text{Δ}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.025m_A}{\left(a_e{\left[i_{\text{pl}}+\frac{m_{\text{pl}}}{3}\right]}^{-2b_e}+C_e{i}_{\text{pl}}^{-2d_e}\right)}^{-\frac{1}{2}}\hfill \\ \hfill & +12.65\frac{m_{\text{pl}}{m}_{\text{ld}}}{m_A^2}{e}_{\text{pl}}\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.025m_A}(a_i{i}_{\text{pl}}^{-2b_i}\exp \left(\frac{-{\left(i_{\text{pl}}/gi\right)}^2}{2}\right)\hfill \\ \hfill & {+c_i{\left[\frac{i_{\text{pl}}}{40}\right]}^{-2d_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=80e_{\text{pl}}^2\exp \left\{-e_{\text{pl}}^2\frac{m_{\text{pl}}}{0.26}\right\}{15}^{m_{\text{pl}}}\left(20+11m_{\text{pl}}-m_{\text{pl}}^2\right),$$(4) $$b_e=0.3m_{\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(2-3e_{\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(2-3e_{\text{pl}}\right)\left(5+e_{\text{pl}}^2{\left[m_{\text{pl}}+2\right]}^3\right)},$$(10) $$d_i=-3+2e_{\text{pl}},$$(11)

and $$g_i=\sqrt{\frac{3m_{{}_{\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}}\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{\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}\sin \left({\theta }_{\text{pl}}\right)}{r_{\text{pl}}^2+{\tilde{r}}^2+2r_{\text{pl}}\tilde{r}\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}=-\sin \left({\theta }_{\text{pl}}\right)\cos \left({\phi }_{\text{pl}}\right)\frac{\partial \Phi }{\partial r_{\text{pl}}}-\cos \left({\theta }_{\text{pl}}\right)\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}=-\sin \left({\theta }_{\text{pl}}\right)sin\left({\phi }_{\text{pl}}\right)\frac{\partial \Phi }{\partial r_{\text{pl}}}-\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}=-\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}}\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)\cot \left({\theta }_{\text{pl}}\right)\sqrt{\frac{\tilde{r}}{r_{\text{pl}}\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-2r_{\text{pl}}\tilde{r}\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{\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)\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 N-body 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)=2i_B\left|sin\left(2\pi f_{d}t\right)\right|\\ {\text{Ω}}_d\left(t\right)=\frac{\pi }{2}+{\text{Ω}}_B-\left(2\pi f_{d}t\mathrm{mod}\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}{1m_{\odot }}\right),$$(25)

with a_B and m_B being the semi-major 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}\exp \left(-t/T_0\right)<{10}^{-9}{m}_0/\text{yr,}$$(26)

where m₀ denotes the mass of the central star, m_disk,0 the disk mass at t = 0 (fixed here to 8 m_Jup), and T₀ = 2.8 × 10⁵ 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 planet-planet 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 j-th planet, respectively, and mj the mass of the j-th 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 Laplace-Lagrange 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₁/a₂ 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+4c_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+4c_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 }{1-3q\sqrt{\alpha }/b_{3/2}^{\left(1\right)}\left(\alpha \right)}\ll 1$$(32)

with q = m₂/m₁. 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.

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 mean-motion resonances (MMRs) via disk migration. We then evolve the simulations for another 100 Myr, in order to study the long-term 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⁻³, 0.1, 0.3, and 0.5) are considered here. The semi-major 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 long-term 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⁻³ and i_B = 50° and planetary masses of m₁ = 1.44 M_Jup and m₂ = 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 θ₁ = 2λ₂ − λ₁ − ϖ₁ and θ₂ = 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 non-resonant 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 aB = 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 eB = 10−3 and iB = 50° for the binary companion and a1 = 6.7 au, m1 = 1.44 MJup, a2 = 18.9 au, and m2 = 3.55 MJup for the planets.

Fig. 2 Same as Fig. 1, but for the initial parameters eB = 0.5, iB = 50°, a1 = 5.7 au, m1 = 2.77 MJup, a2 = 18.15 au, and m2 = 2.99 MJup.

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 planet-planet 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 long-term 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 eB = 0.1, iB = 50°, a1 = 6 au, m1 = 4.11 MJup, a2 = 12.74 au, and m2 = 3.47 MJup.

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 two-planet 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 two-planet 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 one-planet 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 two-planet 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.

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 two-planet 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.

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 inclination-type resonance. This mechanism is known to operate at large eccentricities in two-planet 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⁵ yr when the system enters an inclination-type 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{Ω}}_1$ with the evolution of the mutual inclination (in logarithmic scale). We see that the inclination increases in correlation with the libration of the inclination-type 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 spin-orbit 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 inclination-type resonance during the disk phase. The initial parameters are eB = 0.5, iB = 0.001°, a1 = 6.86 au, m1 = 1.64 MJup, a2 = 17.1 au, and m2 = 3.79 MJup.

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 long-term 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 semi-major 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 semi-major 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 semi-major axes of the binary companion considered here. As expected from Eqs. (27), (30), and (31), we see that the smaller the semi-major 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. 8 Same as Fig. 4, but for the different semi-major axes of the binary companion.

Fig. 9 Same as Fig. 5, but for aB = 250 au.

Fig. 10 Eccentricity and pericenter argument distributions at the end of the simulations for single planet systems (left panels) and two-planet systems (right panels) when aB = 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 long-term 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 one-planet 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 two-planet 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 one-planet 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 one-planet 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 disk-induced precession cannot counteract it. Note that some planets may also have been expelled from the disk’s midplane due to inclination-type resonances, as illustrated in Fig. 6. For binaries with inclinations larger than 40° (blue color), the distribution of the mutual disk-planet 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 semi-major 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 semi-major 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 aB = 250 au (in the Laplace plane reference frame), during the disk phase. The initial parameters (with respect to the disk plane) are eB = 0.5 and iB = 50° for the binary companion and a1 = 10.22 au, m1 = 1.39 MJup, a2 = 25.33 au, and m2 = 1.07 MJup 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 aB = 250 au. Blue (resp. orange) bins are for binary orbital inclinations above (resp. below) 40°.

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 aB = 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 semi-majoraxis, for different semi-majoraxes 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 iB 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 self-gravity. 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 self-gravity 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 Laplace-Lagrange theory to estimate the nodal precession of the disk, we studied in Fig. 14 the uniformity of the disk for the different semi-major 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 self-gravity 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 self-gravity 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 Laplace-Lagrange 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 semi-majoraxes considered for the binary companion. The differences observed for a_B = 250 au mainly arise during the long-term 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 two-planet 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 quasi-periodic orbits and orbits close to stable periodic orbits, and diverges with time for irregular orbits. In practice, the simulations were made with the symplectic Wisdom-Holman 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 semi-major axes of the two planets. The two-planet 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 two-planet 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 long-term 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 Semi-major axis of the outer planet with respect to the one of the inner planet for the two-planet systems formed in our simulations with the binary at aB = 250 au (left panel) and at aB = 1000 au (right panel). The color code refers to the MEGNO value.

7 Comparison with Observations

In this section we compare our results with the 187 S-type 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 S-type 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 Jupiter-like planets with a mass higher than 0.5 MJup. As noticed in Kaib et al. (2013), planets in S-type 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 two-planet systems are weakly perturbed during their long-term 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 quasi-circular 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 S-type 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 multi-planet 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 long-term 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 disk-planets 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⁵−10⁶ 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 disk-planet 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 S-type 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 planet-planet numerical experiments, and finally to include tides and spin dynamics to properly assess the distribution of spin-orbit 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 MIS-project). Computational resources have been provided by the Consortium des Equipements de Calcul Intensif, supported by the FNRS-FRFC, the Walloon Region, and the University of Namur (Conventions No. 2.5020.11, GEQ U.G006.15, 1610468 et RW/GEQ2016).

References

All Tables

All Figures

Fig. 1 Typical evolution of a system of two giant planets with a wide binary at aB = 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 eB = 10−3 and iB = 50° for the binary companion and a1 = 6.7 au, m1 = 1.44 MJup, a2 = 18.9 au, and m2 = 3.55 MJup for the planets.

In the text

	Fig. 2 Same as Fig. 1, but for the initial parameters eB = 0.5, iB = 50°, a1 = 5.7 au, m1 = 2.77 MJup, a2 = 18.15 au, and m2 = 2.99 MJup.
In the text

	Fig. 3 Same as Fig. 1, but for the initial parameters eB = 0.1, iB = 50°, a1 = 6 au, m1 = 4.11 MJup, a2 = 12.74 au, and m2 = 3.47 MJup.
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 inclination-type resonance during the disk phase. The initial parameters are eB = 0.5, iB = 0.001°, a1 = 6.86 au, m1 = 1.64 MJup, a2 = 17.1 au, and m2 = 3.79 MJup.
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 semi-major axes of the binary companion.
In the text

	Fig. 9 Same as Fig. 5, but for aB = 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 two-planet systems (right panels) when aB = 250 au.
In the text

	Fig. 11 Typical evolution of a system of two giant planets with a binary at aB = 250 au (in the Laplace plane reference frame), during the disk phase. The initial parameters (with respect to the disk plane) are eB = 0.5 and iB = 50° for the binary companion and a1 = 10.22 au, m1 = 1.39 MJup, a2 = 25.33 au, and m2 = 1.07 MJup 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 aB = 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 aB = 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 semi-majoraxis, for different semi-majoraxes 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 iB is fixed to 60°.
In the text

	Fig. 15 Semi-major axis of the outer planet with respect to the one of the inner planet for the two-planet systems formed in our simulations with the binary at aB = 250 au (left panel) and at aB = 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