© ESO, 2015

1. Introduction

According to observational surveys, about half of solar-type stars reside in multiple stellar systems (Raghavan et al. 2010) with the frequency declining to roughly 30% for less massive M stars (Lada 2006). It is presumed that this frequency is higher among young stars (Reipurth 2000) and the subsequent decay may be due to dynamical instability or gravitational encounters with other stars. The presence of circumstellar disks around both components of binaries does not seem to be strongly affected by the gravitational perturbations of the companion. Spatially resolved observations of disks in binaries in the Orion nebula cluster (Daemgen et al. 2012) suggest that the percentage of circumstellar disks around individual components of binary systems is about 40%, i.e. only slightly lower than that for single stars (roughly 50%). On the one hand, the secondary star perturbations may lead to disk truncation (Artymowicz & Lubow 1994, 1996; Bate 2000), reducing its lifespan. On the other hand, the presence of a circumbinary disk may feed the truncated circumstellar disks through gas streams so that their dissipation timescale does not vary significantly when compared to single stars (Monin et al. 2007), and they have a chance to form planets.

Until now, 97 exoplanetary systems have been discovered around multiple stars (Rein 2012a)¹, triggering statistical comparisons with planets around single stars (Desidera & Barbieri 2007; Roell et al. 2012). The analysis indicates that planetary masses are higher for stellar separations smaller than 100 au, while for wider binaries the physical and orbital parameters appear very similar. In effect, planet formation in close binary configurations is a complex problem (see Thebault & Haghighipour 2014, for a review), while it is expected that for large separations the effects of the binary companion on the planet growth may be less relevant. However, if the orbit of the companion star is significantly inclined in relation to the initial protostellar disk plane, which is possibly coplanar with the star equator, this may lead to important dynamical consequences for the final orbit of the planet, in particular for its inclination in relation to the star equator. In fact, it has been suggested (Batygin 2012) that the evolution of the protoplanetary disk under the perturbations of the binary companion may be responsible for the observed spin-orbit misalignment of some exoplanets. According to Triaud et al. (2010) and Albrecht et al. (2012), about 40% of hot Jupiters have orbits that are significantly tilted in relation to the equatorial plane of the star. A fundamental requirement of the model of Batygin (2012) is that the planet resides within the disk during its evolution. In this way it will continue its migration by tidal interaction with the disk and, at the same time, follow the inclination evolution of the disk. Once the disk is dissipated its migration stops, while the evolution of the inclination continues, although in a different fashion, until the binary companion is possibly stripped away leaving the planet on an inclined orbit.

Previous studies seem to suggest that a giant planet is indeed forced to evolve within its birth disk even in the presence of external perturbing forces. Any relative inclination between the planet and the disk plane, induced either by a planet-planet scattering event or by resonances (Thommes & Lissauer 2003), is quickly damped by the disk, which leads to a realignment of the planet’s orbit with the disk plane (Cresswell et al. 2007; Marzari & Nelson 2009; Bitsch & Kley 2011; Rein 2012b; Teyssandier et al. 2013). However, the precession rate of the planetary orbit that results from its interaction with the binary star is faster than for the disk. As a consequence, if the damping of the disk on the planet orbit is not strong enough to keep it within the disk plane, the planet will evolve independently under the secular perturbations of the secondary star and develop a significant relative inclination in relation to the disk. The assumed Type II migration will not occur in this case since the planet is not embedded in the disk and cannot open a gap.

In a recent paper, Xiang-Gruess & Papaloizou (2014) investigate the evolution of a disk and massive planet under the influence of an inclined binary companion. According to their smoothed particle hydrodynamic (SPH) simulations, the planet and disk maintain approximate coplanarity during the evolution of the system, and the planet quickly migrates towards the star and stops in a close orbit. These findings seem to confirm the model described in Batygin (2012) for small separations between binary components. To test these findings further, we performed numerical simulations with a set-up similar to that adopted in Xiang-Gruess & Papaloizou (2014), and used two different SPH codes: VINE (Wetzstein et al. 2009; Nelson et al. 2009) and PHANTOM (Price & Federrath 2010; Lodato & Price 2010). Our runs extend over a much longer timescale, and we adopted a higher resolution. The simulations show that the initial coplanarity is maintained only for a limited amount of time and that the planet detaches from the disk plane once and for all as soon as the disk completes a quarter of its precession period. The secular perturbations of the binary companion in an inclined orbit overcome the damping force of the disk and the planet evolves independently of the disk. A significant mutual inclination develops between the disk and the planet orbital plane, but this does not halt the orbital migration. When the planet crosses the disk, it is still affected by the gas via dynamical friction and its inward drift continues, even if at a slower rate compared to what is computed in the initial stages of the system evolution by Xiang-Gruess & Papaloizou (2014).

In Sect. 2 we study the theoretical model of the precession rates in a star+disk+planet+star system and compare evolution timescales. In Sect. 3 we describe the model setup for the numerical simulations in detail. Then, in Sects. 4 and 5 we show the results of the long-term high resolution models for different binary inclinations, and in Sect. 6 we discuss our results and their implications.

2. Relevant timescales

In this section, we compare the timescales of precession of the angular momentum of a planet, computed within a pure N-body problem, to those estimated for the disk. This may provide important clues as to the forces that try to separate the planet from the disk and provide better insight into the results of our SPH simulations. To compute the precession timescales, we consider a planet of mass M_p, initially on a circular orbit around a star of mass M_⋆ at a semi-major axis a, with a distant binary companion of mass M_b, semimajor axis a_b, inclination i_B, and eccentricity e_b, which we assume to be 0 at the beginning. The perturbations of the companion star are particularly strong for inclinations higher than i_B ~ 39° when the Kozai-Lidov mechanism (Kozai 1962; Lidov 1962) forces wide oscillations of both the inclination and eccentricity of the planet whose values are tied by the conservation of ${\mathit{L}}_{\mathrm{z}}\mathrm{=}\sqrt{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\mathrm{)}}\cos \mathrm{\left(}{\mathit{i}}_{\mathrm{p}}\mathrm{\right)}$. The maximum eccentricity attainable by the planet is ${\mathit{e}}_{\max }\mathrm{\simeq }\sqrt{\mathrm{1}\mathrm{-}\mathrm{5}\mathit{/}\mathrm{3}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{lb}}^{\mathrm{0}}}$. However, short-range forces like general relativity and tidal distortions (which we ignore here) and the dumping effect of the protoplanetary disk tend to reduce this value.

The Kozai cycles occur at a characteristic period given by (Kiseleva et al. 1998) $\begin{array}{ccc}{\mathit{P}}_{\mathrm{KZ}}& \mathrm{\simeq }& \frac{{\mathit{P}}_{\mathrm{b}}^{\mathrm{2}}}{{\mathit{P}}_{\mathrm{p}}}{\mathrm{(}\mathrm{1}\mathrm{-}{{\mathit{e}}_{\mathrm{b}}^{\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{3}\mathit{/}\mathrm{2}}\frac{{\mathit{M}}_{\mathit{\star }}\mathrm{+}{\mathit{M}}_{\mathrm{b}}}{{\mathit{M}}_{\mathrm{b}}}\\ & \mathrm{=}& \frac{\mathrm{2}\mathit{\pi }}{{\mathrm{\Omega }}_{\mathrm{p}}}{\left(\frac{{\mathit{a}}_{\mathrm{b}}}{\mathit{a}}\right)}^{\mathrm{3}}\left(\frac{{\mathit{M}}_{\mathit{\star }}}{{\mathit{M}}_{\mathrm{b}}}\right){\mathrm{(}\mathrm{1}\mathrm{-}{{\mathit{e}}_{\mathrm{b}}^{\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{3}\mathit{/}\mathrm{2}}\end{array}$(1)(see also Ford et al. 2000), where P_p and P_b are the orbital periods of the planet and binary. The dynamical evolution of the planet can also be viewed from the point of view of the orbital angular momentum. According to Storch et al. (2014), the planetary orbital angular momentum vector precesses around the binary axis () at a rate that, in the absence of tidal dissipation, is approximately given by ${\mathrm{\Omega }}_{\mathrm{pb}}\mathrm{\simeq }\frac{\mathrm{3}\mathit{\pi }}{\mathrm{2}}{\mathit{P}}_{\mathrm{KZ}}^{-1}\cos {\mathit{\theta }}_{\mathrm{pb}}^{\mathrm{0}}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{0}}^{\mathrm{2}}}\left[\mathrm{1}\mathrm{-}\mathrm{2}\left(\frac{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{0}}^{\mathrm{2}}}{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}}\right)\frac{{\sin }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{pb}}^{\mathrm{0}}}{{\sin }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{pb}}}\right]\mathit{,}$(2)where e₀ is the initial eccentricity of the planet, θ_pb is the inclination of the planetary orbit compared to the binary orbit, and ${\mathit{\theta }}_{\mathrm{pb}}^{\mathrm{0}}$ is its initial value. The precession of the orbital angular momentum translates approximately into an oscillation of the planet’s inclination in relation to the initial plane, which is the disk plane. As such, the timescales estimated by Eqs. (1) and (2) are then comparable.

The disk is also affected by the binary companion gravity and the disk axis precesses around the binary axis , assuming that it satisfies the condition for a coherent behaviour as a solid body, with a period given by $\begin{array}{ccc}{\mathit{P}}_{\mathrm{d}}^{\mathrm{prec}}& \mathrm{\simeq }& {\mathit{P}}_{\mathrm{out}}{\left(\frac{{\mathit{a}}_{\mathrm{b}}}{{\mathit{r}}_{\mathrm{out}}}\right)}^{\mathrm{3}}\left(\frac{{\mathit{M}}_{\mathit{\star }}}{{\mathit{M}}_{\mathrm{b}}}\right)\frac{\mathrm{1}}{\mathit{K}\cos {\mathit{\theta }}_{\mathrm{db}}}\\ & \mathrm{\simeq }& \mathrm{2}\mathit{\pi }{\left(\frac{{\mathit{a}}_{\mathrm{b}}^{\mathrm{6}}}{\mathit{G}{\mathit{M}}_{\mathit{\star }}{\mathit{r}}_{\mathrm{out}}^{\mathrm{3}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\left(\frac{{\mathit{M}}_{\mathit{\star }}}{{\mathit{M}}_{\mathrm{b}}}\right)\frac{\mathrm{1}}{\mathrm{3}\mathit{/}\mathrm{8}\cos {\mathit{\theta }}_{\mathrm{db}}}\end{array}$(3)(see Bate et al. 2000; Lai 2014), where r_out and P_out are the radius and the Keplerian period of the outer disk edge, respectively, and θ_db is the inclination of the disk compared to the orbital plane of the binary. The relative precession rate is given by ${\mathrm{\Omega }}_{\mathrm{db}}\mathrm{\simeq }\mathrm{-}\frac{\mathrm{3}}{\mathrm{8}}{\left(\frac{\mathit{G}{\mathit{M}}_{\mathit{\star }}{\mathit{r}}_{\mathrm{out}}^{\mathrm{3}}}{{\mathit{a}}_{\mathrm{b}}^{\mathrm{6}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\left(\frac{{\mathit{M}}_{\mathrm{b}}}{{\mathit{M}}_{\mathit{\star }}}\right)\cos {\mathit{\theta }}_{\mathrm{db}}\mathit{.}$(4)Even in this case, the precession can be translated into a periodic oscillation of the disk inclination in relation to the initial plane, which can be assumed to have been the equatorial plane of the primary star.

If we now compare the two precession rates, those of the planet and the disk, we find that the initial ratio between them, under the assumption that the planet formed in the midplane of the disk (${\mathit{\theta }}_{\mathrm{pb}}\mathrm{=}{\mathit{\theta }}_{\mathrm{pb}}^{\mathrm{0}}\mathrm{=}{\mathit{\theta }}_{\mathrm{db}}$) in a circular orbit (e = e₀ = 0), is given by $\begin{array}{ccc}{\left(\frac{{\mathit{P}}_{\mathrm{d}}^{\mathrm{prec}}}{{\mathit{P}}_{\mathrm{p}}^{\mathrm{prec}}}\right)}_{\mathrm{ini}}\mathrm{=}{\left(\frac{{\mathrm{\Omega }}_{\mathrm{pb}}}{{\mathrm{\Omega }}_{\mathrm{db}}}\right)}_{\mathrm{ini}}& \mathrm{\simeq }& \mathrm{4}\mathit{\pi }{\mathit{P}}_{\mathrm{KZ}}^{-1}{\left(\frac{{\mathit{a}}_{\mathrm{b}}^{\mathrm{6}}}{\mathit{G}{\mathit{M}}_{\mathit{\star }}{\mathit{r}}_{\mathrm{out}}^{\mathrm{3}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\frac{{\mathit{M}}_{\mathit{\star }}}{{\mathit{M}}_{\mathrm{b}}}\\ & \mathrm{=}& \mathrm{2}{\left(\frac{\mathit{a}}{{\mathit{r}}_{\mathrm{out}}}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{=}\mathrm{2}{\left(\frac{\mathit{a}}{{\mathit{a}}_{\mathrm{b}}\mathit{/}\mathrm{3}}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathit{,}\end{array}$(5)where, in the last step, we assume that the binary truncates the disk maximum radius to one-third of the binary separation (because of tidal interaction), which is a good approximation for binary separations less than ~300 au (Artymowicz & Lubow 1994).

For t ≠ t₀ and e₀ = 0 the ratio evolves as $\frac{{\mathit{P}}_{\mathrm{d}}^{\mathrm{prec}}}{{\mathit{P}}_{\mathrm{p}}^{\mathrm{prec}}}\mathrm{\simeq }\mathrm{2}\frac{\cos {\mathit{\theta }}_{\mathrm{lb}}^{\mathrm{0}}}{\cos {\mathit{\theta }}_{\mathrm{db}}}{\left(\frac{\mathit{a}}{{\mathit{a}}_{\mathrm{b}}\mathit{/}\mathrm{3}}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}\left[\mathrm{2}\left(\frac{\mathrm{1}}{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}}\right)\frac{{\sin }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{lb}}^{\mathrm{0}}}{{\sin }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{lb}}}\mathrm{-}\mathrm{1}\right]\mathit{.}$(6)If we consider a typical case of an equal-mass binary, made of solar-type stars on a circular orbit and separated by 100 au, and a planet orbiting the primary with a = 5 au, and embedded in a disk with an outer radius of 30 au (see Eq. (5)), then the two periods differ by about a factor of ten. This is a large difference and the disk damping must overcome the tendency of the planet to evolve on a slower timescale before moving out of the disk plane.

3. Model description

In this section we briefly describe the two numerical algorithms that we use to model the system planet+disk+binary companion and outline the differences that might affect our long-term results.

3.1. SPH codes and model set-up

We used two different SPH codes to model the evolution of a Jupiter-sized planet that is embedded in the disk that surrounds the primary star of an inclined binary system:

• VINE (Wetzsteinet al. 2009; Nelsonet al. 2009), which is a hybrid N-body/SPH code and has been updated to improve momentum and energy conservation, as described in Picogna & Marzari (2013);

• PHANTOM (Price & Federrath 2010; Lodato & Price 2010), a modern SPH code that models the massive bodies as sink particles (Bate et al. 1995).

Both these codes solve the hydrodynamic equations in the Lagrangian approach by replacing the fluid with a set of particles (see Price 2012, for a review).

Equation of state

A locally isothermal equation of state, similar to the one described in Pepliński et al. (2008), is adopted in all simulations: ${\mathit{c}}_{\mathrm{s}}\mathrm{=}\frac{{\mathit{h}}_{\mathrm{s}}{\mathit{r}}_{\mathrm{s}}{\mathit{h}}_{\mathrm{p}}{\mathit{r}}_{\mathrm{p}}}{\mathrm{\left[}\mathrm{\right(}{\mathit{h}}_{\mathrm{s}}{\mathit{r}}_{\mathrm{s}}{\mathrm{)}}^{\mathit{n}}\mathrm{+}\mathrm{(}{\mathit{h}}_{\mathrm{p}}{\mathit{r}}_{\mathrm{p}}{\mathrm{)}}^{\mathit{n}}{\mathrm{]}}^{\mathrm{1}\mathit{/}\mathit{n}}}\sqrt{{\mathrm{\Omega }}_{\mathrm{s}}^{\mathrm{2}}\mathrm{+}{\mathrm{\Omega }}_{\mathrm{p}}^{\mathrm{2}}}\mathit{,}$(7)where r_s, r_p are the distances of the fluid element from the primary star and the planet, respectively, and h_s, h_p are the circumstellar and circumplanetary disk aspect ratios. Moreover, Ω_s and Ω_p are the circular Keplerian orbits of a fluid element around the star and planet, and n = 3.5 is a non-dimensional parameter chosen to smoothly join the equation of state near the planet (c_s = h_pr_pΩ_p for r_p ≪ r_s) with that of a flat circumstellar disk (c_s = h_sr_sΩ_s).

For the simulations in this work we chose a disk aspect ratio of h_s = 0.037 and a circumplanetary scale height of h_s = 0.6.

Accretion

VINE treats the stars and planet as N bodies. The mass accretion has been modelled for the stars inside a radius of 0.5 au, but not for the planet.

In the simulations performed with the PHANTOM code, the planet and stars are all modelled as Lagrangian sink particles (Bate et al. 1995). A sink particle is evolved as an SPH particle but it only experiences the gravitational force. A gas particle that comes within its accretion radius, R_acc = 0.075 au, can be accreted if

• it is inside the Hill sphere of the sink particle;

• its specific angular momentum is less than what is required to form a circular orbit at the accretion radius;

• it is bound more to the candidate sink particle than to any other sink particle.

Viscosity

An artificial viscosity term is introduced in SPH codes to correctly model shock waves that inject entropy into the flow over distances that are much shorter than a smoothing length and to simulate the evolution of viscous disks. The term broadens the shock over a small number of smoothing lengths and correctly resolves it, ensuring at the same time that the Rankine-Hugoniot equations are satisfied. In this way it prevents discontinuities in entropy, pressure, density, and velocity fields.

The implementation of the viscosity term is slightly different in the two numerical codes: VINE adopts the standard formulation introduced by Monaghan & Gingold (1983) where, for approaching particles (), an artificial viscosity term is introduced in the momentum and energy equations as ${\mathrm{\Pi }}_{\mathit{ab}}\mathrm{=}\{\begin{array}{c}\\ & {{v}}_{\mathit{ab}}\mathrm{·}{{r̂}}_{\mathit{ab}}\mathrm{\le }\mathrm{0}\\ & {{v}}_{\mathit{ab}}\mathrm{·}{{r̂}}_{\mathit{ab}}\mathit{>}\mathrm{0}\mathit{,}\end{array}$(8)where μ_ab is a velocity divergence term ${\mathit{\mu }}_{\mathit{ab}}\mathrm{=}\frac{\mathit{h}{{v}}_{\mathit{ab}}\mathrm{·}{{r}}_{\mathit{ab}}}{{{r}}_{\mathit{ab}}^{\mathrm{2}}\mathrm{+}{\mathit{\eta }}^{\mathrm{2}}{\mathit{h}}^{\mathrm{2}}}\mathit{,}$(9)with η² = 0.01. This implementation conserves total linear and angular momentum, and the viscosity vanishes for a rigid body rotation. The symbol α is the linear term (bulk viscosity) that dissipates kinetic energy as particles approach each other to reduce subsonic velocity oscillations following a shock. The symbol β is the quadratic term (von Neumann-Richtmyer-like viscosity) that converts kinetic energy to thermal energy, therby preventing mutual penetration of particles in shocks.

PHANTOM uses a more general formulation of dissipative terms (Monaghan 1997). It is built on an analogy with Riemann solvers, where the dissipative terms of conservative variables (density, specific momentum, and energy), which experience a jump across a shock front, are multiplied by eigenvalues similar to signal velocities (v_sig). The viscosity term for the momentum equation becomes ${\mathrm{\Pi }}_{\mathit{ab}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathit{\alpha }{\mathit{v}}_{\mathrm{sig}}{{v}}_{\mathit{ab}}\mathrm{·}{{r̂}}_{\mathit{ab}}}{{\mathit{\rho }}_{\mathit{ab}}}\mathit{,}$(10)where the signal velocity is given by ${\mathit{v}}_{\mathrm{sig}}\mathrm{=}\{\begin{array}{c}\\ & {{v}}_{\mathit{ab}}\mathrm{·}{{r̂}}_{\mathit{ab}}\mathrm{\le }\mathrm{0}\\ & {{v}}_{\mathit{ab}}\mathrm{·}{{r̂}}_{\mathit{ab}}\mathit{>}\mathrm{0.}\end{array}$(11)For the energy evolution, an additional signal velocity v_sig,u is adopted, and the viscosity term is defined as ${\left(\frac{\mathrm{d}\mathit{u}}{\mathrm{d}\mathit{t}}\right)}_{\mathrm{diss}\mathit{,a}}\mathrm{=}\mathrm{-}\sum _{\mathit{b}}\frac{{\mathit{m}}_{\mathit{b}}}{{\mathit{\rho }}_{\mathit{ab}}}\left[\frac{\mathrm{1}}{\mathrm{2}}\mathit{\alpha }{\mathit{v}}_{\mathrm{sig}}\mathrm{(}{{v}}_{\mathit{ab}}\mathrm{·}{{r̂}}_{\mathit{ab}}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}{\mathit{\alpha }}_{\mathrm{u}}{\mathit{v}}_{\mathrm{sig}\mathit{,}\mathrm{u}}{\mathit{u}}_{\mathit{ab}}\right]{{r̂}}_{\mathit{ab}}\mathrm{·}{\mathrm{\nabla }}_{\mathit{a}}{\mathit{W}}_{\mathit{ab}}\mathit{,}$(12)where the signal velocity for the energy jumps is chosen to be as in Wadsley et al. (2008).

To properly compare these two different approaches and derive a Shakura & Sunyaev-like viscosity value for the different simulations, we adopt the relations given by Meru & Bate (2012)${\mathit{\alpha }}_{\mathrm{SS}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{20}}\mathit{\alpha }\frac{\mathit{h}}{\mathit{H}}\mathrm{+}\frac{\mathrm{3}}{\mathrm{35}\mathit{\pi }}\mathit{\beta }\left(\frac{\mathit{h}}{\mathit{H}}\right)2$(13)for the Monaghan & Gingold (1983) formalism, and ${\mathit{\alpha }}_{\mathrm{SS}}\mathrm{=}\frac{\mathrm{31}}{\mathrm{525}}\mathit{\alpha }\frac{\mathit{h}}{\mathit{H}}\mathrm{+}\frac{\mathrm{9}}{\mathrm{70}\mathit{\pi }}\mathit{\beta }\left(\frac{\mathit{h}}{\mathit{H}}\right)2$(14)for the Monaghan (1997) one, where h is the averaged smoothing length, and H is the disk scale height.

The initial average value of α_SS for our VINE run is 0.02, similar to that in Xiang-Gruess & Papaloizou (2014). The α_SS in the PHANTOM run is about twice as large, with an average initial value of 0.04.

Initial conditions

The scenario we modelled includes a binary system composed of two equal-mass stars m_p = m_s = 1 M_⊙. One of them (which we call primary) harbours a protoplanetary disk with a Jupiter mass planet M_p = 1 M_J embedded in it, which is on an orbit with an initial semi-major axis a_p = 5 au. Different inclinations between the two components of the binary have been studied in a reference frame centred on the primary star. The disk is modelled with 600 000 SPH particles, approximately three times the number used by Xiang-Gruess & Papaloizou (2014), and it extends from 0.5 to 30 au, with a surface density profile $\mathrm{\Sigma }\mathrm{=}{\mathrm{\Sigma }}_{\mathrm{0}}{\mathit{r}}^{-0.5}\mathit{,}$(15)where Σ₀ is defined so that the total disk mass is M_D = 0.01 M_⊙.

The choice of the initial parameters are such that the disk will precess as a solid body. In fact, according to Papaloizou & Terquem (1995) and Larwood et al. (1996), the condition for this behaviour is ${\mathit{P}}_{\mathrm{out}}\mathit{/}{\mathit{P}}_{\mathrm{d}}\mathrm{\le }\mathit{H}\mathit{/}\mathit{R}\mathit{.}$(16)In our initial set-up, this condition is fully satisfied since H/R is 0.037, while the ratio between ω_p/ Ω_d is approximately 0.01. This is also confirmed by the calculations of Xiang-Gruess & Papaloizou (2013). As a reference frame, we adopt the initial plane of the disk and planet orbit while the binary orbit has different initial inclinations defined in relation to this plane. The initial semi-major axis of the companion star is a_b = 100 au, and it is set on a circular orbit.

4. Initial mutual inclination of 45°

The model where the inclination between the binary orbit and the initial disk+planet plane is 45° is assumed as the standard model. The simulation with VINE was halted after 11 600 orbits of the planet (125 000 yr), while that with PHANTOM was stopped after about half that period (60 000 yr). The run with VINE required approximately eight months of CPU on a 24-processor node, while the run with PHANTOM required four months on a 32-processor machine. In both simulations we find that, in the first 6000 yr, the behaviour is similar to that as described in Xiang-Gruess & Papaloizou (2014), but soon after, the dynamical evolution substantially changes. The inclination of the planet i_p, which initially closely follows that of the disk i_d, effectively departs from i_d and grows at a lower rate. In Fig. 2, both i_p and i_d are shown as a function of time for the run with VINE and PHANTOM. While the periodic oscillations of i_d continue with an almost unaltered frequency of about 2 × 10⁴ yr, the planet inclination i_p grows on a much longer timescale, reaching a maximum after about 8 × 10⁴ yr. There are clear indications of interactions between the planet and the disk, particularly when the two inclinations are comparable, e.g. after 2.5 × 10⁴, 5 × 10⁴. At these times the planet crosses the disk and interacts with its gas particles.

There are some minor differences between the VINE and the PHANTOM models that are possibly related to the different handling of viscosity and accretion on the planet. A circumplanetary disk develops in the initial phase of the run with VINE when i_p ~ i_d, but later on when the planet departs from the disk and periodically crosses it with high mutual inclination, it is dissipated. In both simulations the dynamical behaviour of the planet shows clear competition between the gravitational force of the binary and the interaction with the disk, but the binary secular perturbations finally dominate. The inclination oscillations of the disk show a slow damping in both simulations (Fig. 2), a behaviour already suggested by Martin et al. (2014). Meanwhile, the inclination of the planet roughly follows a pure N-body trend, shown as a dotted light blue line in Fig. 2, even if i_p is strongly perturbed by the disk. The detachment between the planet and the disk can also be seen in Fig. 1, where a 3D rendering of the system is displayed at an optical depth of ~2.

Fig. 1 3D rendering of the disk + planet system in the VINE run. The plot is created through a ray-tracing process, from an observer point of view placed at ~150 au from the disk. The disk is shown at an optical depth of ~2, defined as the cross section per SPH particle mass.

Fig. 2 Evolution of the disk and planet inclination with respect to the initial disk+planet plane in the two runs of VINE and PHANTOM, respectively. The dashed magenta (VINE) and green (PHANTOM) lines show the planet inclination ip, while the continuous black (VINE) and red (PHANTOM) lines illustrate the disk inclination id. After about 6000 yr, the planet evolution departs from that of the disk in both simulations and ip grows at a slower rate compared to id. The light blue line illustrates the evolution of the planet inclination in a pure N-body problem star-planet-companion star.

As discussed in Sect. 2, the distinct evolution of the planet and disk is related to their different precession timescales around the binary axis, which translates into the inclination evolution shown in Fig. 2. As an additional indication of the different evolution of the planet in relation to the disk, in Fig. 3 we illustrate the evolution of the precession angle of the disk L_D and planet L_p angular momentum around that of the binary star L_B, which is calculated in relation to the initial orbital plane (Xiang-Gruess & Papaloizou 2014; Larwood et al. 1996) $\cos {\mathit{\beta }}_{\mathrm{p}}\mathrm{=}\frac{{{L}}_{\mathrm{D}}\mathrm{\times }{{L}}_{\mathrm{B}}}{\mathrm{|}{{L}}_{\mathrm{D}}\mathrm{\times }{{L}}_{\mathrm{B}}\mathrm{|}}\mathrm{·}{u}\mathit{,}$(17)where u is the fixed unit reference vector for the initial orbital plane. The evolution timescale of the planet precession rate is ~10 times slower than that of the disk, confirming the theoretical prediction.

Fig. 3 Evolution of the disk and planet precession angle in relation to the initial disk+planet plane, as defined in Eq. (17). The dashed green line shows the planet precession angle βp, while the continuous red line illustrates the disk’s one βd.

To test whether a potential warping of the disk as a result of binary perturbations is fully responsible for the decoupling of the planet from the disk in our models (Terquem 2013), we checked the disk shape every 300 planetary orbits. In Fig. 4 we show the disk evolution at different times, and there is no evidence of marked warping. In addition, we performed a simulation where a smaller disk was considered.

In Fig. 5 the PHANTOM code was used to model the evolution of a disk with the same mass as our standard disk but extending only to 15 au in radial distance. Because it is less affected by the binary perturbations that extend well within the tidal truncation radius (Artymowicz & Lubow 1994), the warping for such a disk is expected to be negligible. Even here, the planet decouples from the disk in a way very similar to that shown in Fig. 2, albeit on slightly different timescales because of the disk’s different configuration and density.

Initially, the smaller disk has a superficial density that is about six times higher than that of our standard disk. As a consequence, when the planet decouples from the disk, the repeated crossing of the disk plane led to a different dynamical evolution since the frictional force depends on the local disk density. The increase of only 25% in the precession period of the disk, which initially extends over 15 au instead of the expected factor 2 (on the basis of Eq. (3)), is due to the radial spreading of the disk caused by the strong binary perturbations. The standard disk with an initial outer radius of about 30 au cannot spread beyond the tidal truncation limit, and as a consequence, its outer radius remains constant. On the other hand, the smaller disk has room to expand in response to the binary gravitational perturbations, extending after one cycle (about 2 × 10⁴ yr) beyond 25 au.

The outcomes of the numerical simulations rule out warping as being responsible for the decoupling of a planet from its disk in our dynamical configuration since, according to (Terquem 2013), a severe warping is needed to affect the planet inclination. This mechanism may, however, play an important role in closer binaries, increasing the decoupling rate between planet and disk.

Fig. 4 Disk shape shown with a time interval of 300 planetary orbits, when the inclination of the planet is decoupling from that of the planet. Disk warping is negligible.

Fig. 5 Comparison of the planet and disk inclination evolution in the standard case, where the disk extends to 30 au and a model with a smaller disk extends to 15 au. No warping is expected in the latter case, but the planet decouples from the disk as in our standard case with a larger disk.

In Fig. 6, the planet migration is significantly slower compared to that found by Xiang-Gruess & Papaloizou (2014), Fig. 3. Apart from the initial fast migration rate and the planet still embedded in the disk, when the planet does detach from the disk plane, it is the friction that develops during the periodic crossing of the disk by the planet that dominates the semi-major axis evolution and migration. According to Teyssandier et al. (2013), the planet experiences a significant friction while crossing the disk due to aerodynamic drag and to a dynamical drag from the gravitational scattering of the disk particles. The latter dominates since the ratio between the two friction forces can be expressed as (Teyssandier et al. 2013) $\frac{{\mathit{F}}_{\mathrm{aer}}}{{\mathit{F}}_{\mathrm{dyn}}}\mathrm{~}\frac{\mathrm{1}}{\mathrm{8}\ln \left(\frac{\mathit{H}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{{\mathit{R}}_{\mathrm{p}}}\right)}{\left(\frac{{\mathit{R}}_{\mathrm{p}}}{\mathit{a}}\right)}^{\mathrm{2}}{\left(\frac{{\mathit{M}}_{\mathrm{s}}}{{\mathit{M}}_{\mathrm{p}}}\right)}^{\mathrm{2}}\mathit{,}$(18)where M_p and R_p are the mass and radius of the planet, respectively, and a is its semi-major axis. In our initial configuration where a = 5 au, H(r) = h_sa ~ 0.19 and, assuming a radius of the planet R_p ~ 2R_J ~ 0.001 au, we get a ratio lower than 0.02. In this situation the drag force is mostly due to the particle scattering during the planet crossing. There may also be an angular momentum exchange between the planet and the disk, when the planet approaches the disk, but it probably has less effect on the long-term drift.

The migration rate appears to be more irregular in the run with VINE compared with PHANTOM. This different behaviour is possibly due to the mass accretion onto the planet in the PHANTOM run (the planet is modelled as a sink particle), which is not implemented in the VINE run. In the latter, the planet cannot grow and gas is temporarily trapped and removed around it during the disk crossing, possibly leading to a noisy migration. In the run with PHANTOM, the planet has reached a mass of 1.6M_J after about 4 × 10⁴ yr and then the mass growth is reduced.

Fig. 6 Evolution of the semi-major axis of the planet in the VINE and PHANTOM runs, respectively. The planet migrates because of the friction with the gas as it periodically crosses the disk plane.

5. Initial mutual inclination of 60°

Even in the scenario with i_B = 60° we have run two different simulations, one with VINE and one with PHANTOM. In Fig. 7 we show the evolution of the disk and planet inclinations over a short period compared to the previous simulations, where i_B = 45°. However, it is already clear that the planet abandons the plane of the disk and evolves independently under the action of the companion star in both simulations. The disk influence also appears to be a perturbing force in this scenario. The planet inclination, as in the i_B = 45°, follows that of the disk only for a few 10³ yr, after which the planet detaches from the disk following an N-body behaviour while the disk perturbs the dynamics, especially when the planet crosses the disk. In the two models the planet evolution is different, and this could be ascribed to the diverse way of handling the gas trapped around the planet and to the different behaviour of the disks in each case. After an initial quick increase of the disk inclination, the growth halts at about i_d = 90°, and damped oscillations are afterwards observed in both cases. However, the increase in inclination appears more irregular in the run with VINE, and stronger damped oscillations are observed immediately after the first maximum. In the run with PHANTOM, the evolution appears more regular and the damping more progressive. This different behavior may be ascribed to the different viscosity that seems to play a more relevant role on the evolution of the disk inclination when i_B is higher.

The damped oscillations recall, to some extent, those observed in Martin et al. (2014). However, the parameters of our model significantly differ from theirs so the exact same behaviour is not expected. In particular, Martin et al. (2014) use: 1) different temperature and density profiles; 2) their viscosity parameters are set to have a constant value of α_SS all over the disk; and 3) their disk is ten times less massive than ours. As a consequence, the evolution of the disk inclination is only qualitatively similar, but in both scenarios, damped oscillations are observed. An initial hint of damped oscillations of the disk inclination is also seen in the model with i_B = 45° (Fig. 2) even if, as already observed, it appears feeble.

A potential contribution to the strong damping in the i_B = 60° case is possibly related to the exchange of disk mass between the two stars. After about 2 × 10⁴ yr from the beginning of the simulation, the secondary star gains a gaseous disk whose mass is about 5% that of the primary disk. The mass exchange is significantly lower in the i_B = 45° case, where the acquired disk has a mass of only 0.05% of the initial disk mass. As a consequence, the viscosity of the disk also plays a role in an indirect way by forcing the mass exchange for higher values of i_B , which in turn leads to a damping of the disk inclination oscillations.

Fig. 7 Evolution of the disk and planet inclination with respect to the initial disk+planet plane in the case with iB = 60° in the VINE and PHANTOM runs, respectively. The symbols are the same as in Fig. 2.

Fig. 8 Evolution with time of the disk in the PHANTOM model. Some mass leaves the primary disk and forms a ring around the secondary star. This behaviour was already observed in Picogna & Marzari (2013). No significant warping is observed.

Even in this case we tested the potential warping of the disk by plotting the disk evolution every 300 orbits of the planet. In this more inclined configuration we observe some mass transfer from the disk around the primary to the secondary star where a ring of gas slowly builds up. It is interesting to note that the disk eccentricity is higher than in the case where i_B = 45° and comparable to what has been observed in Martin et al. (2014). This is why the mass flux from the primary towards the secondary star is more marked. In any case, the amount of mass extending out of the primary disk is not enough to cause a significant gravitational perturbation on the planet.

6. Discussion and conclusions

We have shown that a Jupiter-like planet that is embedded in a circumstellar disk will evolve almost independently of the disk in response to the perturbations of a misaligned binary stellar companion for values of the semi-major axis a_B around 100 au. After an initial coupled evolution, and for inclinations of the binary plane of the order of 45° and 60°, the planet dynamically detaches from the disk, and its inclination oscillates on a significantly longer timescale compared to that of the disk. This occurs because the N-body secular perturbations of the companion star dominate the damping force of the disk, which tends to drag the planet back into its median plane. Our simulations have three times the resolution of those performed by Xiang-Gruess & Papaloizou (2014) and were performed with two different, up-to-date, SPH codes, VINE and PHANTOM.

Our results show that one of the fundamental requirements in the model of Batygin (2012), which is used to explain the observed spin-orbit misalignment of some exoplanets, is not met for small binary separations and a reasonable choice of physical parameters of the disk and planet. However, the migration of the planet appears to occur at a significant rate even if the orbit of the planet and the disk plane are not aligned. This is the result of the friction with the gas as planet repeatedly crosses the disk. As a consequence, the inward drift of the planet results not from Type I/II migration, as argued for bigger separations by Crida & Batygin (2014), but from the dynamical friction with the disk. This could possibly lead to similar conclusions to those in Batygin (2012), which concern the long-term evolution of the planet spin axis and semi-major axis migration, but the physical mechanism would be different. The formation of hot Jupiters with their spin misaligned in relation to the stellar equator would, in this case, not occur because the planet follows the disk in its precession, except through a combination of N-body Kozai cycles and repeated crossings of the disk plane by the inclined planet that evolves out of the disk. Even if the disk oscillations are damped over a long time to a given value of inclination, the continuous crossing of the disk by the planet, following Kozai-Lidov cycles, would cause a significant migration.

A deeper exploration of the parameter space is of course needed to have a more complete picture of the systems’ behaviour. We have only explored one particular step of the system evolution, while previous steps also need to be investigated; for example, when does the planet begin to evolve independently from its birth disk? At which stage of its growth? And for which value of mass?

We have shown in our simulations that a Jupiter-mass planet exits almost immediately from the disk, but does this also occur for smaller mass planets? Possibly not, since the final infall of gas on the planet core is a quick process that occurs on a shorter timescale than the decoupling between the planet and disk inclination. However, the problem of growing a planet core in a disk that is perturbed by an inclined binary is a critical process. The condensation of dust into larger pebbles and, eventually, kilometer-sized planetesimals may be strongly perturbed, particularly when the planetesimals decouple from the gas. For small inclinations between the circumprimary gas disk and the binary orbital plane, Xie & Zhou (2009) and Xie et al. (2010) show that planetesimal accretion may be fast because size segregation among planetesimals favours low velocity impacts.

For larger inclinations, the nodal longitude randomisation forced by the companion star may lead to the dispersion of the planetesimal disk, thereby reducing the chances of planet formation. Marzari et al. (2009b), however, show that for binary semi-major axes larger than 70 au the nodal longitude randomisation timescale is longer than the planetesimal-accretion process and thus core formation might occur. As a consequence, planet formation might occur within the disk, even in the presence of an inclined binary, up until the planet reaches a Jupiter mass or even larger. After that, it would decouple from the disk and follow an independent evolution.

Additional simulations are needed to explore the problem more fully, and different initial values need to be sampled for the planet mass, the binary orbital elements, and the disk parameters (such as mass, viscosity, temperature distribution, density profile) which can influence the propagation of perturbations. Running more numerical models may help us to understand not only the response of the disk to the binary perturbations in the long term and the possible occurrence of damping, as well as the dynamical evolution of the planet. The problem with these types of simulations is the large amount of CPU required to run an SPH code in this configuration.

Even if the migration of the planet in this scenario is the result of the friction with the gas of the disk, it remains an unresolved problem as to how to drive the planet very close to the star. Isothermal simulations of disks in binaries (Kley et al. 2008; Marzari et al. 2009a) show the formation of an elliptical hole in the density distribution close to the primary star in the long term. This could evidently halt the planet migration once the planet moves within the hole where either Type I/II migration or, as we suggest in this paper, friction with the gas during the crossing of the disk would stop. However, more recent simulations with radiative disks do not show the formation of such an internal low-density region (Müller & Kley 2012; Marzari et al. 2012; Picogna & Marzari 2013). This is an additional angle that should be investigated to better understand the dynamical evolution of inclined hot Jupiters in misaligned binary star systems.

Acknowledgments

We thank an anonymous referee for his useful comments and suggestions, and Daniel Price for having given us access to the PHANTOM code. Many of our plots were made with the SPLASH software package (Price 2007). G. Picogna acknowledges support through the German Research Foundation (DFG) grant KL 650/21 within the collaborative research programme, “The first 10 Million Years of the Solar System”. Some simulations were performed on the bwGRiD cluster in Tübingen, which is funded by the Ministry for Education and Research of Germany and the Ministry for Science, Research and Arts of the state of Baden-Württemberg, and the cluster of the Forschergruppe FOR 759, “The Formation of Planets: The Critical First Growth Phase”, funded by the DFG.

References

All Figures

	Fig. 1 3D rendering of the disk + planet system in the VINE run. The plot is created through a ray-tracing process, from an observer point of view placed at ~150 au from the disk. The disk is shown at an optical depth of ~2, defined as the cross section per SPH particle mass.
In the text

Fig. 2 Evolution of the disk and planet inclination with respect to the initial disk+planet plane in the two runs of VINE and PHANTOM, respectively. The dashed magenta (VINE) and green (PHANTOM) lines show the planet inclination ip, while the continuous black (VINE) and red (PHANTOM) lines illustrate the disk inclination id. After about 6000 yr, the planet evolution departs from that of the disk in both simulations and ip grows at a slower rate compared to id. The light blue line illustrates the evolution of the planet inclination in a pure N-body problem star-planet-companion star.

In the text

	Fig. 3 Evolution of the disk and planet precession angle in relation to the initial disk+planet plane, as defined in Eq. (17). The dashed green line shows the planet precession angle βp, while the continuous red line illustrates the disk’s one βd.
In the text

	Fig. 4 Disk shape shown with a time interval of 300 planetary orbits, when the inclination of the planet is decoupling from that of the planet. Disk warping is negligible.
In the text

	Fig. 5 Comparison of the planet and disk inclination evolution in the standard case, where the disk extends to 30 au and a model with a smaller disk extends to 15 au. No warping is expected in the latter case, but the planet decouples from the disk as in our standard case with a larger disk.
In the text

	Fig. 6 Evolution of the semi-major axis of the planet in the VINE and PHANTOM runs, respectively. The planet migrates because of the friction with the gas as it periodically crosses the disk plane.
In the text

	Fig. 7 Evolution of the disk and planet inclination with respect to the initial disk+planet plane in the case with iB = 60° in the VINE and PHANTOM runs, respectively. The symbols are the same as in Fig. 2.
In the text

	Fig. 8 Evolution with time of the disk in the PHANTOM model. Some mass leaves the primary disk and forms a ring around the secondary star. This behaviour was already observed in Picogna & Marzari (2013). No significant warping is observed.
In the text