© ESO 2021

1 Introduction

About half of the Sun-like stars are part of multiple-star systems (Duquennoy & Mayor 1991; Raghavan et al. 2010). Since the first discoveries of giant planets in binary stars during the past decade (e.g., GL86; Queloz et al. 2000, γ Cephei; Campbell et al. 1988; Hatzes et al. 2003), their number continues to increase, and today more than 100 S-type planets (also known as circumprimary planets) are known that move around one stellar component (Schwarz et al. 2016). Most of these are found in wide binaries with separations larger than 500 AU. Giant planets in wide binaries have slightly higher eccentricities than around single stars (Kaib et al. 2013) and some reported eccentricities are even higher than 0.85.

The question of the formation and long-term evolution of these planets remains open, since the stellar companion can strongly affect the planet formation process. Extensive studies have analyzed the long-term stability and habitability of giant planets detected in binary star systems (e.g., Rabl & Dvorak 1988; Holman & Wiegert 1999; Haghighipour 2010; Bazsó et al. 2017). Planet formation has principally been investigated for close binaries in which giant planets have been detected despite a very strong perturbation of the stellar companion. For S-type planets in wide binaries, owing to a higher stellar separation, it is reasonable to expect that the binary companion would have a more limited (but still significant) impact on the planet formation process.

The different stages of planet formation can be affected by a binary companion (see, e.g., Thébault & Haghighipour 2015, for a review). Firstly, during the formation of the protoplanetary disk, the disk can be truncated by the stellar companion (Artymowicz & Lubow 1994; Savonije et al. 1994). Observations have shown that disks in close binaries tend to be less frequent and less massive than around single stars (Kraus & Ireland 2012). Secondly, the intermediate stage of kilometer-sized planetesimal accretion is extremely sensitive to stellar companion perturbations (Thébault et al. 2006). The formation of planets in close binaries or in the presence of a massive and inclined distant perturber (inducing the Lidov-Kozai excitations of the protoplanetary disk; e.g., Martin et al. 2014; Fu et al. 2015; Zanazzi & Lai 2017) is still a theoretical issue, since even moderate dynamical perturbations can generate high-velocity impacts and thus inhibit the formation of planetary embryos (Lissauer 1993). However, Batygin et al. (2011) showed that disk self-gravity can be sufficient to suppress the Lidov-Kozai oscillations and maintain orbital coherence and planetary growth in the presence of a massive and inclined distant perturber. The decisive role that disk gravity plays in planetesimal dynamics was also confirmed by Rafikov & Silsbee (2015). Lastly, it was shown that in the late stage of planetary accretion, from Lunar-sized embryos to fully formed planets, the embryo accretion region roughly corresponds to the stability domain (e.g., Quintana et al. 2007). Haghighipour & Raymond (2007) analyzed the formation of terrestrial planets in close binaries, assuming the existence of a giant planet already formed further out in the disk and a hundred planetary embryos (Moon- to Mars-sized bodies) formed in the region between the primary and the giant planet. These authors showed that planet formation can be efficient in moderately close binary planetary systems, and the final assembly of the planets and their water content strongly depend on the parameters of the stellar companion.

Regarding the formation of giant planet systems in binaries, much attention was directed toward the formation of hot Jupiters in binary star systems via Lidov-Kozai cycles and tidal friction (e.g., Fabrycky & Tremaine 2007). The planet-planet scattering mechanism was considered for Jupiter-mass planets orbiting the central star of a close binary system (Marzari et al. 2005). Several hydrodynamical simulations of giant planets evolving in gas disks were performed for binary systems. Without providing an exhaustive list, we cite the following works. The evolution of planets embedded in circumbinary disks was considered for close binary stars, for instance, in Kley & Haghighipour (2014) for Kepler-38. Xiang-Gruess & Papaloizou (2014), Picogna & Marzari (2015), and Lubow & Martin (2016) studied the evolution of a giant planet-disk system in the presence of an inclined binary companion and found that a substantial misalignment between the orbit of the planet and the disk generally occurs. Martin et al. (2016) studied the formation of giant planets in Lidov-Kozairesonance with a highly misaligned binary companion, starting from a coplanar planet-disk system configuration.

In the present work, we aim to realize a dynamical study of the influence of a wide binary companion on the evolution of a single giant planet migrating in the protoplanetary disk in the Type II regime. We particularly focus on how the Lidov-Kozai resonance articulates with the planetary migration when the stellar companion is highly inclined. Using a symplectic N-body code including eccentricity and inclination damping formulae to properly model the influence of the disk, we carry out 3200 numerical simulations with different initial eccentricity and inclination values of the binary companion. We follow the evolution of the giant planet when embedded in the protoplanetary disk and pursue the analysis well after the dispersion of the disk.

The paper is organized as follows. In Sect. 2, the code and initial settings of the simulations are described. In Sect. 3, we present the results of the simulations by focusing on typical planetary evolutions as well as the final parameter distributions for the planet. In order to dynamically characterize the planetary evolutions, the results are discussed using a quadrupolar Hamiltonian approach in Sect. 4. Finally, our conclusions are given in Sect. 5.

2 Setup of the simulations

In this work, we focus on the evolution of a giant planet embedded in a protoplanetary disk undergoing Type II migration around one of the two stars of a wide binary star system (S-type planet). We denote with subscripts 0 for the central star, B for the binary companion, and pl for the planet. The code used for the simulations and the initial parameters of our study are described in the following.

For our study we adapted the SyMBA N-body code (Duncan et al. 1998) to wide binaries by following the strategy of Chambers et al. (2002), based on a new set of coordinates (wide-binary coordinates). This last set induces a splitting of the Hamiltonian of the three-body problem composed of two stars and a planet, into three terms suitable for symplectic integration. Working in a Hamiltonian formalism, the family of symplectic integrators uses the symplectic structure of the problem to keep the Hamiltonian well conserved for a long-time scale (i.e., the energy error remains bounded), even when using a large time step. The evolution of the three terms can be performed separately via the Baker-Campbell-Hausdorff formula (Bourbaki 1972), upon which symplectic schemes such as that of Laskar & Robutel (2001) have been built. However these schemes require that one part of the Hamiltonian is far larger than the other parts. The splitting that is classically used for the solar system (i.e., a term for the Keplerian motion and two terms depending on the position and momentum, respectively, for the mutual perturbation between the planets) is not suitable for binary systems because the perturbation is no longer smaller than the Keplerian part. The wide-binary coordinates of Chambers et al. (2002) are defined such that all the large terms of the Hamiltonian can be incorporated into a single part of the Hamiltonian, making a symplectic integration possible with a leapfrog or SABA₁ (Laskar & Robutel 2001) integrator.

The N-body code also includes the Type II migration of the giant planet caused by the angular momentum exchange with the protoplanetary gas disk (Goldreich & Tremaine 1979), mimicked by a suitable Stokes-type drag force added in the equations of motion. We used the same approach as in Sotiriadis et al. (2017), to which we refer for more details. The migration of a planet in Type II regime is on the same timescale as the viscous accretion time but, when the mass of the planet is comparable to the mass of the material in its vicinity, the migration rate scales with the ratio of the planetary mass over the local disk mass (Ivanov et al. 1999; Nelson et al. 2000; Crida & Morbidelli 2007). Thus the timescale for Type II regime adopted in this work is written as $${\tau }_{\text{mig}}=\frac{2}{3}{\alpha }^{-1}{h}^{-2}{\Omega }_{\text{pl}}^{-1}\times \text{max}\left\{1,\frac{m_{\text{pl}}}{(4\pi /3)\Sigma (r_{\text{pl}})r_{\text{pl}}^2}\right\},$$(1)

where α = 0.005 is the classical value for the Shakura-Sunyanev viscosity parameter (Shakura & Sunyaev 1973), h = 0.05 the disk aspect ratio, ${\Omega }_{\text{pl}}^{-1}=2\pi a_{\text{pl}}^{(-3/2)}$ the orbital frequency of the planet, a_pl the semimajor axis of the planet, m_pl the mass of the planet, and Σ ∝ r^−0.5 the surface density profile of the disk. In the simulations, we fixed the disk inner and outer edges to R_in = 0.05 AU and R_out = 30 AU. The local disk mass is considered as the mass of the disk between 0.2a_pl and 2.5a_pl. The code also includes a smooth transition in the gas-free inner cavity using an hyperbolic tangent function $\text{tan}h\left(\frac{r-R_{\text{in}}}{\Delta r}\right),$ where Δr = 0.001 AU, following Matsumoto et al. (2012).

Regarding the damping of the planetary eccentricity and inclination caused by the disk, we adopted the damping formulae of Bitsch et al. (2013) achieved by three-dimensional hydrodynamical simulations of protoplanetary disks with embedded high-mass planets. The analytic formulae depend on the local mass disk, planetary mass, eccentricity, and inclination evaluated throughout the integration. These formulae are conceived for planets between 1 M_J and 10 M_J. We note that it is rather unclear if a highly inclined planet continues to migrate on a viscous accretion timescale. Nevertheless, we adopted this recipe since, when initially embedded in the disk, the planet only occasionally reached high inclination with respect to the disk plane during the protoplanetary disk phase (only a few simulations are concerned and only for a very brief period of time). No inclination damping is applied when i_pl < 0.5°.

Following Papaloizou & Larwood (2000), the orbital migration and the eccentricity and inclination damping effects are added to the acceleration of the planet as follows: $$a_{\text{disk}}=-\frac{v_{\text{pl}}}{{\tau }_{\text{mig}}}-2\frac{(v_{\text{pl}}\text{ cdotp}{r}_{\text{pl}})r_{\text{pl}}}{r_{\text{pl}}^2{\tau }_{\text{ecc}}}-2\frac{(v_{\text{pl}}\text{ cdotp}k)k}{{\tau }_{\text{inc}}},$$(2)

where r_pl is the position of the planet, v_pl the velocity of the planet, k the unit vector in the vertical direction, τ_mig the timescale for Type-II regime (see Eq. (1)), and τ_ecc and τ_inc the timescales for eccentricity and inclination damping (issued from Bitsch et al. (2013), see their formulae (11)–(16)), respectively. To keep the symmetry of the symplectic algorithm, the damping forces are implemented in the code as in Lee & Peale (2002) as follows:

$$f(t_{k+1})=E_{\text{disk}}\left(\frac{\tau }{2}\right)\exp \left(T\tau \right)E_{\text{disk}}\left(\frac{\tau }{2}\right)f(t_k),$$

where T is the linear operator associated to the Hamiltonian of the N-body problem and $E_{\text{disk}}(\cdot )$ stands for the evolution due to the disk. The time step of our simulations is fixed to 0.001 yr during the migration phase and 0.01 yr after the dispersal of the disk, since the evolution of the system is then characterized by secular effects. The simulations are followed for 100 Myr.

The initial disk mass in our simulations is fixed to 8 M_Jup and decreases exponentially through the evolution of the system, with a dispersal time of ~ 1 Myr. More specifically, the disk is considered as dissipated and the interaction is neglected when $$\text{d}{m}_{\text{disk}}(t)/\text{d}t=m_{\text{disk}}(0)\exp (-t/T_0)/T_0<{10}^{-9}{m}_0/\text{yr},$$

where m₀ denotes the mass of the central star and T₀ = 2.8 × 10⁵ yr. The self-gravity of the disk is not considered in this work.

The initial orbital elements and masses for the planet and the binary stars adopted in the simulations are summarized in Table 1. These are expressed in the disk plane reference frame. The wide binary companion has a mass of 1 M_⊙ and is located at 500 AU. The planetary mass m_pl, eccentricity e_pl, inclination i_pl, and argument of the pericenter ω_pl follow a uniform distribution. For the binary companion, we considered four values of the eccentricity e_B (namely 10⁻⁵, 0.1, 0.3, and 0.5) as well as eight values of the inclination i_B (from ~0° up to 70°). In total, by picking randomly 100 values for each uniform distribution, we carried out 3200 simulations. The dynamical evolutions observed in the simulations are presented in the next section.

3 Results

In this section, we describe in detail the long-term evolution of the planets as well as their orbital characteristics at the end of the simulations. Particular attention is given to the capture in the Lidov-Kozai resonance in the case of highly inclined binary companions.

Fig. 1 Typical evolution of a Lidov-Kozai resonant system (libration of the pericenter argument in the invariant plane reference frame). The initial parameters are mpl = 2.02 MJup, ωpl = 16°, eB = 0.3, and iB = 70°.

3.1 Lidov-Kozai resonance

In our simulations, a strong influence of the wide binary companion is observed on the evolution of the planet embedded in the protoplanetary disk, especially for highly inclined binary companions for which the Lidov-Kozai mechanism comes into play (Lidov 1962; Kozai 1962). When referred to the invariant plane orthogonal to the total angular momentum of the system (also called the Laplace plane), the planet in a Lidov-Kozai resonant state exhibits coupled variations in eccentricity and inclination as well as a libration of the pericenter argument around ± 90°. We note that the large orbital variations associated with the Lidov-Kozai resonance take place in a coherent way such that the long-term stability of the system is assured. A typical example of a Lidov-Kozai evolution is shown in Fig. 1, for an inclination i_B = 70° of the binary companion with respect to the disk plane. This inclination value is higher than the critical value $\arccos \sqrt{3/5}\approx {39.23}^{°}$ introduced by Kozai (1962) for the restricted three-body problem. After the migrating phase, the planet is rapidly captured in a Lidov-Kozai resonance, as indicated by the libration of the pericenter argument observed in the last panel of Fig. 1, and high eccentricity and inclination variations are observed.

It is important to note that the libration of the pericenter argument in the invariant plane reference frame is the proper criteria for the identification of the Lidov-Kozai resonant regime among the simulations (Libert & Tsiganis 2009). High eccentricity and inclination variations can also be observed when the planet is not in a Lidov-Kozairesonant state, as shown in Fig. 2. While the planet does experience large variations in eccentricity and inclination, the argument of the pericenter circulates (see the last panel of Fig. 2). For this reason, all the dynamical evolutions displayed in the following are presented in the invariant plane reference frame.

Note that, in practice, a planet will be considered here as locked in the Lidov-Kozai resonance if a libration of the argument of the pericenter is observed at the end of the simulation during three Lidov-Kozai cycles (for the calculation of the period, see e.g. Ivanov et al. 1997) or, in case it exceeds the timescale of the simulation, during the last 8 × 10⁷ yr.

Fig. 2 Typical evolution of a planet that undergoes high eccentricity and inclination variations, while not in a Lidov-Kozai resonance (circulation of the pericenter argument in the invariant plane reference frame). The initial parameters are mpl = 4.34 MJup, ωpl = 13°, eB = 10−5, and iB = 70°.

3.2 Finaldistributions

In this section, we describe the orbital parameters of the planets at the end of the simulations (i.e., 100 Myr). We recall that all the planetary orbital parameters are given in the invariant plane reference frame. Figure 3 shows the planetary inclination as a function of the semimajor axis and the eccentricity. The color code indicates if the planet is in a Lidov-Kozai resonance state at the end of the evolution.

As seen in the left panel of Fig. 3, all the planets end up with a semimajor axis that is smaller than 10 AU as a result of the Type II migration of the planets during the protoplanetary disk phase. The planets found locked in a Lidov-Kozairesonant state (red stars) all have an inclination (with respect to the invariant plane) that is higher than ~ 40°, which comes as no surprise considering the critical value of 39.23° for the establishment of the Lidov-Kozai resonance. Nevertheless, the highly inclined planets are not all in a Lidov-Kozai resonant state at the end of the simulations, as shown by the color code.

The critical inclination value is even more obvious in the right panel of Fig. 3, where the planetary inclinations are presented jointly with the eccentricities. While a pileup around low eccentricities is observed at low planetary inclinations, high eccentricities are only associated with inclinations above ~ 40°. However, although the Lidov-Kozai resonance could contribute to the excitation of the planetary eccentricities, we note that only a small fraction of the planets evolves in a Lidov-Kozai resonance at the end of the simulations. The percentages of final Lidov-Kozai resonant evolutions found among the 3200 simulations are given in Table 2, for each initial inclination value of the binary companion. Let us remark that the percentages of resonance captures are slightly lower when the simulations do not include the disc phase. It seems to indicate that the migration in the disc promotes the establishment of Lidov-Kozai resonant sytems. Additionally, we note that the values in eccentricity and inclination displayed in Fig. 3 follow particular curves (green lines) whose trend will be discussed in Sect. 4.

Since the argument of the pericenter plays a key role in the Lidov-Kozai mechanism, we therefore pay particular attention tothis angle in the following. Figure 4 shows the normalized distribution of the argument of the pericenter of the planets at the end of the simulations. We clearly observe two pileups around 90° and 270°. Interestingly, the highest peaks are found at the borders of the two pileups, namely for pericenter arguments close to 45°, 135°, 225°, and 315°. In the left panel of Fig. 4, the systems in a Lidov-Kozai resonance belong to the pileups at 90° and 270°, as expected. Inthe right panel of Fig. 4, the color code of the histogram refers to the initial inclination values of the binary companion i_B. We can deduce from the plot that the pileups are formed by the highly inclined binary companions (i_B ≥ 40°), since the normalized frequencies of the pericenter arguments for planets with a binary companion initially below 40° are almost uniformly distributed.

Finally, in Fig. 5, we present the inclination of the planet as a function of the argument of the pericenter, at the end of the simulation, with the same color codes as in Fig. 4. We note the existence of two arc-shaped curves. The Lidov-Kozai resonant systems mainly lie around these two curves, as observed in the left panel. In the right panel, we see that the two arc-shaped curves mostly correspond to planets with highly inclined binary companions. We note that the accumulations observed at the borders of the two pileups in the histograms of Fig. 4 can be connected with the particular shape of these curves. The edges of the arc-shaped curves are nearly parallel to the y-axis and thus more systems are gathered at the edges rather than in the central part of the curves (i.e., around 90° and 270°) when considering parts of the x-axis of the same length.

3.3 Initial conditions leading to Lidov-Kozai resonance

We might wonder whether it is possible to identify the initial planetary orbital elements leading to a Lidov-Kozai resonant state for the planet. As can be deduced from Table 2, approximately 36% of the planets with a binary companion inclined at more than 40° end up trapped in a Lidov-Kozai resonance at the end of the simulation. A statistical analysis of the initial parameters of these systems was performed, but no trend was observed. In particular, the initial values of the pericenter argument of the planets captured in the Lidov-Kozai resonance follow an uniform distribution. This is further illustrated in the evolutions of Figs. 1 and 2, which have close initial values of the arguments of pericenter but different resonant behaviors.

Moreover, we carried out a study of the specific time of the capture in the Lidov-Kozai resonance. We present in Fig. 6 the same plot as in the left panel of Fig. 5, but for different times. The top panel shows the planetary inclination versus the pericenter argument values at half of the disk lifetime (0.55 Myr), the middle panel at the dispersion of the disk (1.1 Myr) and the bottom panel at the end of the simulation (100 Myr). We clearly see the accumulation of the planets around the arc-shaped curves through time. Interestingly, the planets firstly gather in the left part of the curves (i.e., around a pericenter argument of 45° and 135°) during the disk phase and spread along the curves in a second phase.

Fig. 3 Inclination of the planets as a function of the semimajor axis (left panel) and the eccentricity (right panel) at the end of the simulations (100 Myr). The planets evolving in the Lidov-Kozai resonance are drawn in red.

Fig. 4 Normalized distribution of the pericenter argument of the planet, at the end of the simulation (100 Myr). The color codes refer to planets evolving in the Lidov-Kozai resonance in the left panel (red color) and to the initial inclination of the binary companion iB in the right panel.

Fig. 5 Inclination of the planet as a function of the argument of the pericenter, at the end of the simulation (100 Myr). The color codes refer to planets evolving in the Lidov-Kozai resonance in the left panel (red color) and to the initial inclination of the binary companion iB in the right panel.

Fig. 6 Same as Fig. 5, for different times: 0.55 Myr (half of the disk lifetime) in the top panel, 1.1 Myr (dispersion of the disk) in the middle panel, and 100 Myr (end of the simulation) in the bottom panel.

4 Dynamical analysis

In this section, we give a dynamical interpretation of the arc-shaped curves observed in Fig. 5. To describe the dynamics of the migrating planets, we draw phase portraits of a simplified Hamiltonian formulation and follow the trajectory of the planets on this portrait until their final location.

Following the work of Bataille et al. (2018), we used a quadrupolar Hamiltonian approach (i.e., the dominating terms of the secular Hamiltonian) that has one degree of freedom and states as follows: $$\begin{array}{ll}\mathcal{H}=\hfill & -\left(5-3\frac{h^2}{{\cos }^2{i}_{\text{pl}}}\right)(1-3{\cos }^2(i_{\text{pl}}))\hfill \\ \hfill & +15\left(1-\frac{h^2}{{\cos }^2{i}_{\text{pl}}}\right)(1-{\cos }^2{i}_{\text{pl}})\cos (2{\omega }_{pl}).\hfill \end{array}$$

Since the angular momentum of the wide binary orbit is much greater than that of the planetary orbit, the dimensionless angular momentum of the planet $h=\sqrt{1-e_{\text{pl}}^2}\cos i_{\text{pl}}$ is nearly constant and is usually referred to as the Kozai constant. For a given value of the Kozai constant, phase portraits of the dynamics can be drawn, consisting of level curves of the Hamiltonian in the plane (ω_pl, e_pl) (or equivalently (ω_pl, i_pl)).

For a low value of the inclination of the binary companion, the Hamiltonian curves in the plane (ω_pl, e_pl) are nearly straight lines. An example of such an evolution is shown in Fig. 7, for an inclination of 30° for the companion star. In the left panel, we see that no significant eccentricity and inclination variation is observed during the long-term evolution of the planet. This evolution is overplotted on the phase portrait given in the right panel of Fig. 7. The time evolution is indicated by the color scale, while the initial planetary parameters are denoted with a magenta star symbol. The value of the Kozai constant h used to compute the phase portrait is fixed to the mean value of h after the migration phase. Thus, at the beginning (dark blue points), the planetary evolution can depart from the Hamiltonian curves represented on the phase portrait (since the system is dissipative owing to the migration), but after the migration phase the trajectory of the planet evolves on a fixed level curve of the analytic model. In Fig. 7, the planet follows a quasi-straight line of constant Hamiltonian, as expected.

For systems with a highly inclined companion (i_B ≥ 40°), the phase space shows the Lidov-Kozai resonance islands that strongly influence the dynamical behavior of the planet. Two typical dynamical outcomes are observed in the simulations and these are both illustrated in Fig. 8. In the left panel, we show the phase portrait for the system of Fig. 1, for which a capture around the Lidov-Kozai equilibrium around 90° takes place. The evolution of the system of Fig. 2 is represented on the phase space of the right panel. No capture in the Lidov-Kozai resonance is observed, and the planet evolves around the two islands with its pericenter argument circulating. Finally, we note that in both cases the planets spend time in libration or circulation around the Kozai equilibria. This gives an explanation for the arc-shaped curves observed in the final distribution of the simulations in Fig. 5.

Having a careful look to the two phase spaces of Fig. 8, one could expect that the points along the horizontal lines inside the two arc-shaped curves in Fig. 5 are associated with Lidov-Kozai resonant motion (red color). The blue points on the horizontal lines correspond to planets that alternate between phases of Lidov-Kozai resonant captures and phases of circulation around the islands. An example of such an evolution is represented on Fig. 9. The planet is first captured around the island centred on ω = 90° at the dispersion of the disk, then after some cycles leaves the resonance to circulate around the islands, before being captured around the second island centred on ω = 270° and startingto alternate between circulation and resonance capture.

Finally, we show that the trend observed in the right panel of Fig. 3 can be explained by the conservation of the Kozai constant (green curves). Indeed, for a given value of the Kozai constant, the ranges of possible eccentricity and inclination values for the planet are fixed by the relation $h=\sqrt{1-e_{\text{pl}}^2}\cos i_{\text{pl}}$. Here the Kozai constant is evaluated at the dispersal of the disk and has been approximated by the quantity cos i_B, since the disk tends to circularize the planetary orbit and to keep it in the disk through eccentricity and inclination damping.

Fig. 7 Left panel: typical evolution of a planet in the presence of a binary companion with an inclination lower than 40°. Right panel: phase portrait of the evolution. The initial parameters are mpl = 4.6 MJup, ωpl = 152°, eB = 0.1, and iB = 30°.

Fig. 8 Phase portraits of different planetary evolutions in the presence of a binary companion at iB = 70°. Left: Lidov-Kozai resonant evolution from Fig. 1. Right: non Lidov-Kozai resonant evolution from Fig. 2.

5 Discussion and conclusions

In this work, we studied the dynamical influence of a wide binary companion on the migration of a single giant planet in the protoplanetary disk. The simulations were made using the SyMBA symplectic N-body integrator modified for binary star systems and eccentricity and inclination damping formulae for the disk influence. By varying the eccentricity and inclination of the binary companion as well as the physical and orbital parameters of the planet, we carried out 3200 numerical simulations and carefully analyzed the orbital configuration of the planets after 100 Myr.

As a consequence of the Lidov-Kozai mechanism, the influence of the wide binary companion is strong for highly inclined companions even during the disk phase. We showed that the orbital parameters of the planets at the end of the simulations gather around two arc-shaped curves in the planes (ω_pl, e_pl). We noticed that only about 36% of the planets with a highly inclined binary companion (i_B ≥ 40°) end up in a Lidov-Kozai resonant state (with the libration of the pericenter argument). Nevertheless, the nonresonant evolutions are strongly influenced by the Lidov-Kozai resonance and present high eccentricity and inclination variations associated with a circulation around the Lidov-Kozai islands. Using an analytical quadrupolar Hamiltonian approach, we realized phase portraits in which the different dynamical evolutions can easily be observed.

It is important to note that in this work we neglected the influence of the companion star on the disk, in particular, the induced nodal precession of the disk. To study the robustness of our results, and in particular the impact that nodal precession could have, we ran several tests by considering different migration rates and eccentricity and inclination damping rates as well as a location of the planet initially outside of the disk. The results of these tests are summarized in Fig. 10 which shows, for the different tests, the same analysis as in Fig. 5, namely, the planetary inclination as a function of the argument of the pericenter at the end of the simulation, in which the color code indicates the initial inclination of the binary companion i_B. In the top panels, the planets are initially inside of the disk plane (i_pl = 0°), while the planets have an initial inclination of i_pl = 20° with respect to the disk plane in the bottom panels. Three modifications for the planetary migration are considered. In the left panels, the migration rates and eccentricity and inclinations damping rates are unchanged. In the middle panels, the Type II migration rate is scaled by the factor cos (i_pl). In the right panels, no eccentricity and inclination damping is included. Thus, in total, six tests were realized, each one consisting in 3200 simulations, as previously detailed. As observed in Fig. 10, the modifications performed have no significant impact on the results and the arc-shaped curves are present in each test.

These extra simulations make us confident in the robustness of the results presented in this work. We obtained similar outcomes when changing the planetary inclination and, therefore, the mutual inclination between the planet and the disk. Therefore, we believe that the addition of the nodal precession of the disk caused by the binary star, inducing a change in the mutual inclination, does not significantly modify the long-term evolution of the planets.

For this study, we did not include additional effects due to the disk and the binary that are likely to have a significant impact on the evolution of the planet. Aside from the nodal precession of the disk caused by the binary star, we did not consider the aspidal and nodal precessions caused by the gravitational potential of the disk, which could dominate over the precession induced by the binary star (e.g., Zanazzi & Lai 2018). Also, the general relativity and the stellar oblateness could compete or suppress the Lidov-Kozai effect. The influence of these additional effects on our results is left for future work.

We believe that since it is currently estimated that about half of the Sun-like stars are part of multiple-star systems, the complete understanding of the influence of a wide binary companion on the planetary eccentricities and inclinations is very important to fully explain the observational parameter distributions of the exoplanets.

Fig. 9 Left panel: typical evolution of a planet alternating between resonance capture and circulation. Right panel: phase portrait of the evolution. The initial parameters are mpl = 4.74 mJup, ωpl = 35°, eB = 0.1, and iB = 60°.

Fig. 10 Same as Fig. 5 for different migration and damping rates. Rates are shown unchanged in the left panels, Type II migration rate scaled by the factor cos (ipl) indicated in the middle panels, and no eccentricity and inclination damping considered in the right panels. Two initial inclination values of the planet with respect to the disk plane are also shown (ipl is fixed to 0° in the top panels and 20° in the bottom panels).

Acknowledgements

The authors would like to warmly thank J. Teyssandier, K. Tsiganis and A. Crida for useful discussions as well as the anonymous referee for his/her comments that improve the present document. The work ofA. Roisin is supported by a F.R.S.-FNRS research fellowship. Computational resources have been provided by the Consortium des Equipements de Calcul Intensif, supported by the FNRS-FRFC, the Walloon Region, andthe 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 Lidov-Kozai resonant system (libration of the pericenter argument in the invariant plane reference frame). The initial parameters are mpl = 2.02 MJup, ωpl = 16°, eB = 0.3, and iB = 70°.
In the text

	Fig. 2 Typical evolution of a planet that undergoes high eccentricity and inclination variations, while not in a Lidov-Kozai resonance (circulation of the pericenter argument in the invariant plane reference frame). The initial parameters are mpl = 4.34 MJup, ωpl = 13°, eB = 10−5, and iB = 70°.
In the text

	Fig. 3 Inclination of the planets as a function of the semimajor axis (left panel) and the eccentricity (right panel) at the end of the simulations (100 Myr). The planets evolving in the Lidov-Kozai resonance are drawn in red.
In the text

	Fig. 4 Normalized distribution of the pericenter argument of the planet, at the end of the simulation (100 Myr). The color codes refer to planets evolving in the Lidov-Kozai resonance in the left panel (red color) and to the initial inclination of the binary companion iB in the right panel.
In the text

	Fig. 5 Inclination of the planet as a function of the argument of the pericenter, at the end of the simulation (100 Myr). The color codes refer to planets evolving in the Lidov-Kozai resonance in the left panel (red color) and to the initial inclination of the binary companion iB in the right panel.
In the text

	Fig. 6 Same as Fig. 5, for different times: 0.55 Myr (half of the disk lifetime) in the top panel, 1.1 Myr (dispersion of the disk) in the middle panel, and 100 Myr (end of the simulation) in the bottom panel.
In the text

	Fig. 7 Left panel: typical evolution of a planet in the presence of a binary companion with an inclination lower than 40°. Right panel: phase portrait of the evolution. The initial parameters are mpl = 4.6 MJup, ωpl = 152°, eB = 0.1, and iB = 30°.
In the text

	Fig. 8 Phase portraits of different planetary evolutions in the presence of a binary companion at iB = 70°. Left: Lidov-Kozai resonant evolution from Fig. 1. Right: non Lidov-Kozai resonant evolution from Fig. 2.
In the text

	Fig. 9 Left panel: typical evolution of a planet alternating between resonance capture and circulation. Right panel: phase portrait of the evolution. The initial parameters are mpl = 4.74 mJup, ωpl = 35°, eB = 0.1, and iB = 60°.
In the text

Fig. 10 Same as Fig. 5 for different migration and damping rates. Rates are shown unchanged in the left panels, Type II migration rate scaled by the factor cos (ipl) indicated in the middle panels, and no eccentricity and inclination damping considered in the right panels. Two initial inclination values of the planet with respect to the disk plane are also shown (ipl is fixed to 0° in the top panels and 20° in the bottom panels).

In the text