THE EVOLUTION OF PLANET–DISK SYSTEMS THAT ARE MILDLY INCLINED TO THE ORBIT OF A BINARY COMPANION

# The Evolution of Planet–disk Systems That Are Mildly Inclined to the Orbit of a Binary Companion

## Abstract

We determine the evolution of a giant planet–disk system that orbits a member of a binary star system and is mildly inclined with respect to the binary orbital plane. The planet orbit and disk are initially mutually coplanar. We analyze the evolution of the planet and the disk by analytic means and hydrodynamic simulations. We generally find that the planet and the disk do not remain coplanar unless the disk mass is very large or the gap that separates the planet from the disk is very small. The relative planet–disk tilt undergoes secular oscillations whose initial amplitudes are typically of order the initial disk tilt relative to the binary orbital plane for disk masses 1% of the binary mass or less. The effects of a secular resonance and the disk tilt decay enhance the planet–disk misalignment. The secular resonance plays an important role for disk masses greater than the planet mass. At later times, the accretion of disk gas by the planet causes its orbit to evolve towards alignment, if the disk mass is sufficiently large. The results have several implications for the evolution of massive planets in binary systems.

accretion, accretion disks – binaries: general – hydrodynamics – planetary systems: formation

## 1 Introduction

Disks around young stars in binary systems may be misaligned with respect to their binary orbital planes. A misaligned disk in a binary system is expected to evolve towards coplanarity due to tidal dissipation associated with turbulent viscosity (Papaloizou & Terquem, 1995; Bate et al., 2000; Lubow & Ogilvie, 2000; King et al., 2013). The alignment process may occur on relatively short timescales in binaries whose separations are small and thus the tidal torques and associated dissipation are strong. The alignment may also be a consequence of initial conditions in the binary formation process. Observations show that the stellar rotation and binary orbital axes are better aligned in closer systems with binary separations (Hale, 1994). These observations provide indirect evidence for binary–disk alignment during the star formation stage in closer binaries and misalignment for wider binaries.

There is some direct observational evidence of disk misalignment with respect to the binary orbital planes in wider binary systems (e.g., Jensen et al., 2004; Skemer et al., 2008; Roccatagliata et al., 2011). For the case of the young binary system HK Tau with a projected separation of about 350 AU, circumstellar disks are observed around each component with one disk edge–on and the other more face–on (Stapelfeldt et al., 1998). Although the orbit of the binary is not known, at least one of the disks must be substantially misaligned to the binary orbital plane. Recent ALMA observations of this system by Jensen & Akeson (2014) suggest that the misalignment between the two disks is . Williams et al. (2014) observed a wide binary (with a projected separation of 440 AU) in Orion and found the misalignment between the projected disk rotation axes to be about . These results suggest that wide binary star systems do not form directly from a single large corotating primordial structure. Instead, they may be subject to smaller scale effects, such as turbulence (Offner et al., 2010; Tokuda et al., 2014; Bate, 2012) that could result in a lack of correlation between the rotational axes of the accreting gas associated with the two stars.

Less direct evidence of noncoplanarity comes from the existence of extrasolar planets whose orbits are tilted with respect to the spin axis of the central star (e.g. Albrecht et al., 2012; Huber et al., 2013; Lund et al., 2014; Winn & Fabrycky, 2014). If such planets reside in binary star systems, these observations suggest that the planets may have formed in disks that are misaligned with the binary orbital plane (e.g., Bate et al., 2010; Batygin et al., 2011; Batygin, 2012).

A misaligned disk will undergo nodal precession due to the torques caused by the companion star. For typical protostellar disk parameters, the disk remains nearly flat and undergoes little warping as it precesses about the binary orbital axis (Larwood et al., 1996). A misaligned disk whose inclination angle with respect to the binary orbital plane is between about and can additionally undergo Kozai–Lidov oscillations (Martin et al., 2014b; Fu et al., 2015). These oscillations cause the disk inclination and eccentricity to vary in time. In this paper we restrict our attention to cases where the disk inclination angle is sufficiently small that these oscillations do not occur and the disk remains circular.

We analyze the orbital evolution of a giant planet that interacts with the gaseous disk in which it forms. Such a planet will open a gap in a disk due to tidal forces and will undergo migration as the gas accretes towards the central star (Lin & Papaloizou, 1986). The accretion is driven by viscous forces in the disk. The planet orbit and disk are taken to be initially coplanar and slightly misaligned with the respect to the binary orbital plane. In this paper, we explore the evolution of the planet and disk. A study along these lines was recently carried out by Xiang-Gruess & Papaloizou (2014) who concluded that coplanarity between the orbital planes of the planet and the disk can be maintained. Very recently, Picogna & Marzari (2015) reported results on SPH simulations with the same initial conditions that ran over longer timescales and with higher resolution. They found that coplanarity is not maintained.

A giant planet–disk system that orbits a single star (not in a binary) is subject to effects of misalignment due to mean motion resonances (Borderies et al., 1984; Lubow, 1992). The misalignment can be suppressed by damping effects associated with disk viscosity (Lubow & Ogilvie, 2001). In this paper, we concentrate on the effects of the binary companion that can bring about misalignment.

For sufficiently low mass disks, we expect that planet–disk coplanarity cannot be maintained, since the precessional torque on the planet caused by the disk is weak compared to the precessional torque on the disk provided by the binary companion. The disk and planet could then precess nearly independently. On the other hand, for sufficiently high mass disks, we would expect that coplanarity can be maintained, since the disk torque on the planet dominates. We explore a range of disk parameters in order to understand the conditions under which coplanarity breaks down. We follow the tilt evolution in time by means of an analytic model of the secular evolution and by SPH simulations.

In Section 2 we consider analytically the secular evolution of a planet and a rigid nonviscous outer disk that are initially on mutually coplanar orbits, but misaligned with respect to the orbit of a binary companion. As a test of the secular theory, in Section 3, we consider a binary system with two planets orbiting one of the stars and compare the analytic secular evolution to that of 4–body simulations. In Section 4 we describe the results of 3D smoothed particle hydrodynamic (SPH) simulations that model a viscous disk that interacts with a planet in a binary system. The planet and disk are initially mutually coplanar, but misaligned with respect to the binary orbital plane. We compare the results with the secular theory of nonviscous disks. Section 5 contains a discussion and Section 6 summarizes our results.

## 2 Secular Evolution of Planet–Disk–Binary Systems

### 2.1 Secular Equations

In this section we consider the secular evolution of a planet and a disk that orbit one member of a binary system. In particular, we model a system that orbits a central star consisting of a planet, a disk that lies exterior to the planet, and a companion star. The orbital planes of the planet and the disk are initially aligned, but misaligned with respect to the binary orbital plane. The planet, the disk, and the binary interact through gravitational forces. The disk interior to the planet provides little torque on the other components in our SPH simulations described in Section 4 and we ignore its effects here. The changes in the orbit of the binary due to its interactions with the planet and disk are expected to be small, since binary angular momentum is much larger than the angular momentum of the planet or disk. The binary orbit is then taken to be fixed. In addition, we assume the binary orbit is circular and Keplerian. The disk is taken to be nonviscous (nondissipative), rigid, and flat (does not warp).

Consider a Cartesian coordinate system in the inertial frame centered on star 1 with the plane defined as the binary orbital plane. The axis is parallel to the line joining the stars at the initial time. We describe the time dependent tilt of the planet and disk relative to the plane with the complex variable notation , where the unit angular momentum vector is denoted by and is the time. We assume that the tilts are small, .

The angular frequency is assumed to be Keplerian for the planet and disk, so that , where is the mass of the star 1 that lies at the disk center and is the distance from the center of star 1. The disk extends from a radius out to and has angular momentum

 Jd=2π∫RoutRinΣ(R)R3Ω(R)dR, (1)

where is the surface density distribution of the disk. The angular momentum of the planet with mass and orbital radius from the primary is given by

 Jp=Mpa2pΩ(ap). (2)

The binary companion star 2 has mass and orbital radius .

We apply the secular evolution equations for the gravitational interactions between slightly misaligned components in Lubow & Ogilvie (2001). In this model, the torques between the components of the system (planet, disk, and binary) are evaluated in the small angle approximation for their relative tilts. In Section 3, we consider a binary system with two planets orbiting one of the stars and compare the analytic secular evolution to that of 4–body simulations. We find that the analytic results are accurate for initial tilt angles up to about , due in part to the small angle approximation that limits the accuracy at larger initial tilt angles. In addition, the components are assumed to remain on circular orbits. Consequently, this formalism cannot be used to analyze planet-disk systems undergoing Kozai-Lidov oscillations, where the relative tilts are large and the planet orbit and disk may acquire substantial eccentricity.

The interaction between the components and is described by a linear model through coupling coefficients denoted by . The evolution equations for the planet tilt and the disk tilt are given by

 JpdWpdt=iCpd(Wd−Wp)−iCpsWp (3)

and

 JddWddt=iCpd(Wp−Wd)−iCdsWd, (4)

where subscripts , and refer to the planet, disk and companion star 2, respectively. The first term on the right-hand side of Equation (3) is the horizontal torque (along the plane of the binary) on the planet due to the disk and the second term is the horizontal torque on the planet due to the binary companion star. The first term on the right-hand side of Equation (4) is the horizontal torque on the disk due to the planet and the second term is the horizontal torque on the disk due to the binary companion.

The coupling coefficients are given by

 Cpd=2π∫RoutRinGMpRΣ(R)K(R,ap)dR, (5)
 Cds=2π∫RoutRinGM2RΣ(R)K(R,a)dR, (6)

and

 Cps=GMpM2K(ap,a), (7)

where the symmetric kernel, with units of inverse length, is given by

 K(Rj,Rk)=RjRk4π∫2π0cosϕdϕ(R2j+R2k−2RjRkcosϕ)3/2. (8)

The kernel contains a singularity as . This singularity is resolved by the finite thickness of the disk . We assume in this paper that the separation between the components is greater than and consequently we do not smooth the kernel.

Consider a time-periodic vector of tilts with frequency

 W(t)=~Wexp(iωt) (9)

where

 ~W=(~Wp~Wd). (10)

The tilts satisfy the matrix equation

 M~W=0, (11)

where

 M=(ωJp+Cpd+Cps−Cpd−CpdωJd+Cpd+Cds). (12)

The frequency can be shown to be real.

Although we assume that the actual tilts are small, for simplicity in this linear problem, we set to solve for the eigenvectors and the eigenfrequencies. The contributions of these eigenvectors to the tilt evolution are determined by initial conditions. We denote the two eigenvectors by for and 2. We have determined these eigensolutions analytically. We apply these analytic solutions to determine the numerical results that we describe below.

With these two eigensolutions for the matrix equation, we can solve any initial value problem. We determine the contributions of the eigenvectors through equation

 c1~W1+c2~W2=~W0 (13)

by solving for constants and , where and and are the initial tilts for the orbits of the planet and disk, respectively, and where denotes the transpose of the vector from row to column form. The evolution of the tilts of the planet-disk components is then given by

 W(t)=c1~W1eiω1t+c2~W2eiω2t, (14)

where describes the tilt evolution of the planet and disk, respectively.

We consider systems with the initial tilts given by This vector represents a planet and a disk whose orbital planes are initially aligned with each other, but are slightly misaligned with respect to the orbital plane of the binary with small . The inclination of each component relative to the binary orbital plane is given by

 ij(t)=|Wj(t)|, (15)

for . The linear model assumes that tilts and are small. The phase angle is given by

 ϕj(t)=tan−1(Im(Wj(t))Re(Wj(t))). (16)

We apply these equations in the following subsections to determine the evolution of some planet-disk-binary systems.

### 2.2 System Parameters

We choose a set of standard system parameters and consider the effects of variations from these values. The parameter values for the disk are selected as plausible and are not tuned to quantitatively match the results of the SPH simulations that are described in Section 4. The results of the secular models demonstrate trends that explain properties of the simulations that are due to gravitational torques.

The surface density of the disk is assumed to follow , where is a constant defined through the constraint that the mass of the disk, , is given. The disk is assumed to extend only exterior to the orbit of the planet. The standard disk mass is taken to be , where is the binary mass.

As discussed in Section 2.1, the binary orbit is taken to be circular. The value of the disk outer radius is taken to be that is about equal to the tidal truncation radius of a disk in a coplanar equal mass binary system (see Paczynski, 1977). The tidal truncation radius for a noncoplanar configuration may be somewhat larger due to the weakening of the tidal torques (Lubow et al., 2015). The planet is taken to have mass and orbital radius . The planet and disk orbit star 1 in an equal mass binary with .

The clearance between the orbit of planet and the inner edge of the disk is determined by the size of the gap. The gap size depends on the level of disk viscosity, planet mass, etc. Also, the ”gap” region is not completely clear of gas. Using the gap density profile corresponding to Bate et al. (2003) for a planet–to–star mass ratio of (a planet orbiting a star), we compute a value of in Equation (5) with . We then determine the effective disk inner radius that gives the same value of , but for a disk with an empty gap region and a sharp disk inner edge. This procedure yields in an effective disk inner radius For a planet–to–star ratio of , as we consider here (), the gap size is expected to be somewhat larger, perhaps by a factor of , based on scaling by the Hill radius, and we adopt a fiducial value of .

### 2.3 Effect of Varying the Disk Mass

Fig. 1 shows the evolution of the inclinations and phase angles in a system with the standard parameters described in Section 2.2, except that we consider variations to the disk mass . Since the secular evolution model is linear, the tilts scale with , and so we plot .

As seen in Fig. 1, the inclinations of the planet and the disk for the smaller disk mass of oscillate away from each other, while each precesses on different timescales. The planet precesses more slowly than the disk because it is farther away from the companion star. For the larger disk mass of , the precession of the planet is clearly affected by its interaction with the disk that causes its precession rate to be close to that of the disk that is somewhat longer than in the left panel. For both disk masses, the inclination of the planet is generally larger than its initial value, while the inclination of the disk is generally smaller than its initial value. The planet and disk spend very little time in a coplanar configuration. Instead, we find generally that .

The level of planet–disk misalignment undergoes periodic oscillations. We measure this misalignment by the amplitude of these oscillations. The left hand panel of Fig. 2 plots the normalized oscillation amplitudes of planet–disk tilt differences, as a function of disk mass. The right hand panel plots the oscillation period as a function of disk mass. The amplitude is determined in two different ways. One set of lines (solid and short–dashed) plots the oscillation amplitude (the maximum value over time) of the difference in the planet and disk tilts relative to the binary orbital plane , normalized by the initial tilt . The amplitude in this case is defined by

 Ab=maxt(|Wp(t)|−|Wd(t)|). (17)

The other set of lines (long-dashed and dashed–dotted) plots the oscillation amplitude of the planet tilt relative to the disk , normalized by the initial tilt . The amplitude in this case is defined by

 Apd=maxt|Wp(t)−Wd(t)|. (18)

This quantity depends on the relative phasing of the planet and disk tilts. For disk masses beyond the peaks in the curves, both ways of measuring the tilt differences coincide because the planet and disk mean precession rates are locked.

In the left panel of Fig. 1, the maximum inclination difference is and this occurs at time of about , corresponding to an oscillation period of about . These values correspond to the oscillation amplitude and oscillation period that lie on the solid lines in the left and right panels of Fig. 2, respectively, for a disk mass of .

The combined planet-disk angular momentum is generally not conserved. In the left panel of Fig. 1, the planet and disk precess independently due to the torque caused by the binary. While precessing at different rates, their individual angular momentum vectors vary differently in time. In the right panel, they precess together while undergoing angular momentum changes due to the binary torque.

We might expect that with increasing disk mass the planet becomes more aligned with the disk. However, there are peaks of in the left panel of Fig. 2. It is somewhat surprising that the level planet-disk misalignment can increase with increasing disk mass. In Section 2.7 we show that this peak is associated with a secular resonance. For a disk inner radius of in Fig. 1, the peak occurs near . For disk masses smaller than this value, the planet and disk precess independently. The tilt oscillation frequency is equal to the difference in the nodal precession frequencies between the disk and planet. However, for increasing disk masses greater than this value, the precessions of the planet and disk become locked and the inclination differences and the oscillation periods decrease. For a smaller disk inner radius (smaller gap), the planet and disk interact more strongly and reach peak values in Fig. 2 at a lower disk mass. This effect is described further in Section 2.4. Note however that except for small values of the disk mass, the inclination difference is generally comparable to or greater than the initial inclination, , for disk masses up to . This result suggests that departures from coplanarity are significant in many cases.

Fig. 3 plots the relative disk–planet tilt as a function of nodal phase difference for the model plotted in Fig. 2 with disk inner radius for different disk masses. For a disk mass of , the planet is only slightly affected by the presence of the disk and they precess independently at different rates, resulting in a fully circulating phase difference. varies in this case because the disk undergoes small amplitude tilt oscillations due to its interaction with the planet. When the planet and disk phase angles are aligned, the relative tilt is zero. When the difference in the phase angles is , the relative tilt is maximum.

For a larger disk mass of , Fig. 3 shows that the interactions lead the system to undergo stronger tilt oscillations with larger planet–disk relative tilts, while still fully circulating in phase. When the disk mass becomes large enough for them begin to precess together (), the phase difference is limited and the system is librating. The mean precession rates for the planet and for the disk over a libration cycle are then equal, unlike the fully circulating case at lower disk masses. The planet and disk mean precession rates are then locked. At this disk mass, the amplitude of relative tilt oscillations is maximum, as discussed above. At the times of both maximum and minimum relative tilt , the disk and planet have the same phase, . For increasingly larger disk masses, there is a decrease in the amplitude of the oscillations and the range of phase . However, the relative tilts can be significant compared to the initial tilt with respect to the binary for disk masses up to one or more percent of the binary mass.

### 2.4 Effect of Varying the Disk Inner Radius

For the higher disk mass especially, the evolution of the disk–planet system is somewhat sensitive to the radius of the inner edge of the disk for a fixed planet orbital radius. The reason is that the kernel defined Equation (8) varies as for , resulting in a variation of in the disk–planet coupling coefficient .

For a smaller disk inner radius (smaller gap), the planet and disk interact more strongly and reach peak relative tilt values in Fig. 2 at a lower disk mass. The left panel of Fig. 4 shows the amplitude of the planet–disk relative inclination oscillations, , as a function of disk inner radius with disk mass . The right panel shows the oscillation period. If the inner edge of the disk is very close to the planet, then the oscillation amplitude may be very small. There is a peak in the oscillation amplitude. If the disk inner radius is smaller than the inner radius at the peak, then the mean precession rates of the planet and the disk are locked. But for a larger inner disk radius, they are more disconnected and precess independently. When the planet and disk do not interact with each other, their maximum relative tilt is twice the initial tilt of that occurs when they are separated in phase by . (As noted below Equation (8), we assume that the gap size is larger than the disk thickness in evaluating . This plot is then valid for cases where .) But unless the gap is small, less than about , the planet and a larger mass disk can be substantially misaligned with misalignment angle greater than about half of their initial tilts relative to the binary orbital plane.

### 2.5 Effect of the Varying the Orbital Radius of the Planet

In Fig. 4 we include plots of the relative disk–planet inclination for a planet at a smaller orbital radius of for different disk inner radii. The planet–disk relative inclinations at the two planet orbital radii are very similar. This plot then shows that the orbital radius of the planet does not strongly affect the amplitude of the oscillations. The more important factor is the disk gap size.

### 2.6 Effect of Varying the Disk Outer Radius

The disk outer radius is controlled by a competition between the disk viscous torques that act to expand the outer parts of the disk with tidal torques that act to truncate the disk. As discussed in Section 2.2, for the standard model we choose the outer disk radius to be near the tidal truncation radius of a disk in a coplanar binary system (Paczynski, 1977). However, in a misaligned binary, the disk can be somewhat larger because the tidal torque that truncates the disk is reduced (Lubow et al., 2015). In Fig. 5 we plot the amplitude of the oscillations of the planet orbital tilt relative to the disk as a function of disk outer radius for the standard model described in Section 2.2. A larger disk experiences a stronger binary torque that leads to greater planet–disk misalignment. The results show that the amplitude of the relative tilt oscillations is a sensitive function of the disk outer radius. However, the relative tilts are substantial over the plotted range of disk outer radii.

### 2.7 Role of Secular Resonances

Secular resonances are known to play an important role in the dynamics of small mass objects, such as asteroids, in multi-planet systems (Murray & Dermott, 1999). Secular resonances occur when there is a matching of the precession frequency of a test particle (e.g., an asteroid) with the frequency of one of the precessional modes of the planetary system. The frequencies involved are much lower than in the mean-motion resonance case, i.e., Lindblad resonances. At such a resonance, the motion of a test particle can be strongly driven by the planets, resulting in a high orbital inclination (for a nodal resonance) or eccentricity (for an apsidal resonance). In the context of this paper, we are concerned with nodal secular resonances. An example of such resonance is the so-called resonance caused by Saturn and Jupiter driving strong vertical motions in the asteroid belt, resulting in a narrow gap at the radius where a secular resonance condition is satisfied among the nodal precession frequencies.

We analytically determined the properties of the peak value of as a function of disk mass, such as is seen in the left panel of Fig. 2. Using analytic expressions for and , we find that this the maximum occurs when

 ωds=ωpd(1+JpJd)+ωps, (19)

where is the precession rate of the disk caused by the binary companion star, is the precession rate of the planet caused by the disk, and is the precession rate of the planet caused by the binary companion star. In this equation, and are linear functions of disk mass , while the other parameters are independent of disk mass. The maximum value of the planet–disk tilt oscillation amplitude for varying disk mass is given by

 Amaxi0=√1+JdJp, (20)

where is the value of the disk angular momentum that occurs when Equation (19) is satisfied.

For , the resonance full-width at half-maximum in terms of disk mass is given by

 δMdMd=4√3√JpJd, (21)

where and are the values of the disk angular momentum and disk mass that occurs when Equation (19) is satisfied. This resonance width estimate assumes that the gap size is maintained in taking this low limit. While this assumption will not likely hold, this estimate is made for comparison with the resonance, as described below.

Equation (19) is essentially a secular resonance condition. If we consider the planet to be a very small mass object so that , then Equation (19) reduces to the resonance condition where we regard the planet as an asteroid that interacts with Jupiter and Saturn, instead of the disk and companion star. is regarded as the sum of precession rates due to the mutual torques involving the dominant components, Jupiter and Saturn. Since the angular momentum of an asteroid is very small compared the angular momentum of Jupiter or Saturn (implying in effect ), we see from Equations (20) and (21) that the inclination change is very large and is confined to a narrow range of perturbing mass at fixed position that translates into a narrow range of radii at fixed mass.

The case of a disk-planet system is somewhat different, since the giant planet mass is not extremely small compared to the disk mass. For the upper curve plotted in left panel of Fig. 2, near the peak value of , where , the angular momentum ratio is . As a result, the value of is not very large and the dimensionless resonance width in terms of disk mass is not small. Consequently, the effects of the resonance extend over a broad range of disk mass.

We define a dimensionless measure of the closeness to the resonance condition (19) (the detuning parameter) by

 D(Md)=ωpd(1+Jp/Jd)+ωpsωds−1, (22)

where and are linear functions of . For positive (negative), the relative nodal phasing of the planet and disk is librating (circulating). For negative, the planet orbit inclination relative to the disk is of order , as a consequence of the circulation. For large and positive, the planet inclination relative to the disk is much less than and alignment occurs. For of order unity, the planet inclination relative to the disk is of order , as a consequence of the secular resonance.

In Fig. 6, we plot as a solid line the case of the dashed and overlapping solid line in the left panel of Fig. 2. The dashed line plots the resonance condition (22) and the dotted line plots the predicted maximum value of in Equation (20) with evaluated at the resonance disk mass (although the plot extends over different disk masses). As expected, the peak misalignment occurs for a disk mass where and therefore where Equation (19) is satisfied. The peak value agrees with the predicted value. In this case, the resonance condition is satisfied for a disk mass that is about twice the planet mass.

The nonmontonic behavior of as a function of disk mass is then a consequence of the secular resonance. That is, the increase in relative planet-disk tilt with increasing disk mass is due to the effects of the secular resonance driven by the disk and binary companion that causes misalignment as measured by to have order unity values ( or larger) in the cases we have considered. The results of this section suggest that maintaining coplanarity by means of gravitational torques between the planet and the disk appears to only be possible if the mass of the disk is very large compared to that of the planet, or if the disk extends very close to the planet. Smaller mass planets open a smaller gap in the disk and are more likely to remain coplanar.

## 3 Two Planets in a Binary System

In this section we consider a binary star system with two misaligned planets around one star. We are essentially replacing the disk in Section 2 with an outer planet. By doing so, we can easily carry out simulations that can be compared with the secular theory of Section 2 to test its accuracy. Of course, a disk is an extended body that can lie closer to both the binary companion and the inner planet than can a point mass object such as a planet. This test then involves a somewhat different situation than having a disk. But it can give a general estimate of the accuracy of the secular theory for the range of parameters involved in the planet–disk case.

### 3.1 Secular Interactions

We adopt an equal mass binary system with and planets and with masses orbiting at radii and , respectively. We apply the methods described in Section 2, but replace the disk with a point mass with mass at radius and relabel the planet as planet . Thus, we replace the angular momentum of the disk component with

 J2=Mp2a2p2Ω(ap2), (23)

replace the coupling coefficient by

 Cp1p2=GMp1Mp2K(ap1,ap2), (24)

and replace by

 Cp2s=GMp2M2K(ap2,a). (25)

The planet orbits are initially coplanar but misaligned with respect to the binary orbital plane. We consider three different initial binary misalignments of , , and . Because the secular model equations are linear, each is just a scaled version of the others. Fig. 7 shows the tilt and phase angle evolution for each planet. The orange line corresponds to the inner planet and the red line the outer planet .

### 3.2 4–body simulations

We describe the numerical simulations where we integrate the gravitational force equations in time. We work in the inertial frame and fix the orbits of the stars. As in the planet–disk secular model, we do not allow the orbit of the binary to evolve under the influence of the other objects. We apply a Cartesian coordinate system with the origin at the binary center of mass, the axis along the line joining the initial positions of the two stars, and the axis along the binary rotation axis. Star 1 that is central to the two planets has position

 r1(t)=a1(cost,sint,0), (26)

where is its binary orbital radius. The companion star 2 has position

 r2(t)=−a2(cost,sint,0), (27)

where is its orbital radius. For an equal mass binary, we have that . The position of the inner planet is denoted by and the outer planet is denoted by . We solve the differential equations for the position of the inner planet

 d2rp1dt2+GM1(rp1−r1)|rp1−r1|3+GM2(rp1−r2)|rp1−r2|3+GMp2(rp1−rp2)|rp1−rp2|3=0 (28)

and for the outer planet

 d2rp2dt2+GM1(rp2−r1)|rp2−r1|3+GM2(rp2−r2)|rp2−r2|3+GMp1(rp2−rp1)|rp2−rp1|3=0. (29)

Planets and begin on the –axis at radii and . from star 1. The planets are given an initial velocity in the plane such that they are on a circular Keplerian but tilted orbits about the central star 1. The resulting evolution of the inclination of the two planets is shown in the black and blue lines in Fig. 7.

The numerical simulations and the secular evolution predicted in the previous section show very similar behavior for low inclination angles in Fig. 7. Thus, we conclude from this section that the secular evolution model in Section 2 provides an approximately accurate description of tilt evolution for small tilts .

## 4 Hydrodynamical Simulations

In this Section we describe 3D hydrodynamical simulations to model planet–disk–binary systems. We use the smoothed particle hydrodynamics (SPH; e.g., Price 2012) code phantom (Price & Federrath, 2010; Lodato & Price, 2010; Nixon, 2012; Nixon et al., 2012, 2013). The binary star system, disk, and planet parameters are summarized in Table 1. We consider an equal mass binary star system, with total mass that has a circular orbit in the - plane with separation . The mass of the planet is and its initial distance from the central star is . We choose the accretion radius about each star to be and about the planet to be . We have found that the simulation results do not vary significantly for smaller values of these accretion radii. With these values we are able to run the simulations faster than with smaller values. Particles that fall within this radius are removed from the simulation and their mass and momentum are added to the accreting object.

The disk is locally isothermal with sound speed and at . With these parameters, and are constant over the disk (Lodato & Pringle, 2007). The Shakura & Sunyaev (1973) parameter is taken to be 0.05. We implement the disk viscosity by adapting the SPH artificial viscosity according to the procedure described in Lodato & Price (2010), using and . The disk is resolved with shell–averaged smoothing length per scale height .

In order to simulate a planet–disk system, we must choose an appropriate initial disk density distribution. In the case of coplanar systems, we can choose a disk profile with or without an initial gap. As the system evolves, a reasonably stable equilibrium gap structure is produced typically in less than 100 planet orbits (e.g., Bate et al., 2003). In case of a planet–disk system that is misaligned with respect to a binary, the initial conditions need to be more accurately established because the disk could precess substantially before an equilibrium gap is established. That is, an initially aligned planet–disk system may not retain this desired initial relative alignment until the equilibrium gap is established. A further problem is that if we start without a gap, the planet can gain substantial mass before gap opening and thereby prevent us from easily controlling the planet mass when an equilibrium gap is established.

To mitigate these problems, we have developed a procedure in which we first simulate a coplanar planet–disk–binary system with a very low mass initial disk, . The disk starts without a gap. The disk initially extends from to . The outer disk radius is close to the tidal truncation radius of a disk in a coplanar binary system (see Paczynski, 1977). The initial surface density follows the power law . We evolve this coplanar simulation for 10 binary orbital periods. By this time, a stable disk gap structure is created while the planet mass and orbit remain nearly unchanged. Fig. 8 shows the resulting surface density profile.

It can be shown that for a locally isothermal, non-selfgravitating disk (as we assume), and fixed orbits of the binary and planet, the disk density profile shape is independent of disk mass. We then rescale the density distribution obtained after 10 binary orbital periods to achieve the desired disk mass. We also tilt both the planet and disk relative to the binary orbital plane by the desired angle so that they are mutually coplanar. This planet–disk configuration then serves as the initial conditions for the subsequent simulation of a planet–disk system that is misaligned with respect to the binary orbital plane.

After the initial coplanar simulation completes in 10 binary orbital periods, much of the disk that lies interior to the orbit of the planet has been accreted onto the central star. The resulting disk lies primarily exterior to the planet orbit. Due to the approximations made in the establishing the initial conditions for the tilted system, the system makes some initial readjustment to the tilt and increased disk mass. The tidal forces are weaker for a tilted disk than a coplanar one (Lubow et al., 2015). Consequently, the disk expands outwards somewhat. However, since the tilt is small, this should not be a large readjustment. Another effect is that the increase in the disk mass results in outward gravitational forces on the planet that cause its orbit to initially expand slightly.

We examine disks with three different masses, , and , for the initially tilted disk. The simulations start with SPH particles in the initial coplanar (pre-tilted) configuration.

We first consider the evolution of a small mass disk with a total mass of . The left panel of Fig. 9 shows the evolution of the system in time. The upper two graphs show the inclination and phase angle evolution for the disk and the planet. The disk remains largely unwarped (flat) during the evolution. We plot the disk evolution at a representative radius of . The black lines show the evolution of the planet and the blue lines the disk. After the first oscillation, does not return to zero but remains approximately periodic. Unlike the oscillations in the secular model, the oscillations in the SPH simulations are damped through the disk viscosity.

In comparing these results to the left panel of Fig. 1, we see that the periods of the tilt oscillations in the SPH simulation and secular model have similar values of about . There is some difference in the tilt amplitudes . They are likely due to differences in our parameter choices for the secular model, such as the disk outer radius. The red line in the middle graph plots the relative inclination between the planet and the disk, . The relative inclination of the planet to the disk reaches values of about 1.8 times or that is similar to the value of about as implied by in Fig. 3 for this disk mass. The fourth graph shows the mass of each component and the bottom graph shows the orbital semi–major axis of the planet. There is little accretion or migration of the planet caused by this small mass disk.

The right panel of Fig. 9 plots the relative disk–planet tilt as a function of nodal phase difference with the individual phases plotted in the second graph of the left panel. The system starts at and The phase wraps from to , as the system evolves. The system is circulating, rather than librating. In both the secular model and the SPH simulation, the relative phasing between the planet and disk is fully circulating (see Fig. 3 and right panel of Fig 9).

Fig. 10 shows the evolution of planet–disk system with a higher initial disk mass of . The tilt oscillations are larger with the higher disk mass, as was seen for the secular model in Fig. 1. There is more migration of the planet with the higher mass disk, but the amount of migration over the course of the simulation is small. The planet gains about 50% more mass by the end of the simulation which is about 1/3 of the amount of gas that has been lost from the disk. Gas giant planets in this general mass range in coplanar planet–disk systems typically accrete most of the gas that flows past their orbits (Lubow & D’Angelo, 2006). As seen in Fig. 1 and Fig. 10, the secular model and the initial oscillation in the simulations show a similar peak value of . On the other hand, the inclination of the disk reaches lower values than predicted by the secular model.

In the right panel of Fig. 1 , the secular model predicts libration at this disk mass, rather than the circulation found in the SPH simulations. It is clear from Fig. 10 that the phasing of the planet and disk is not locked, even at early times. This comparison indicates that peak planet–disk misalignment, as seen in Fig. 6, occurs at a higher value of the disk mass than the secular model suggests, higher than . For , the planet–disk misalignment at higher disk masses is even stronger than estimated by the secular model with the adopted standard parameters. This shift could be a consequence, for example, of underestimating the disk outer radius as for the secular model compared with possibly larger values, as suggested by Fig. 8. The SPH results indicate that the transition from circulation to libration occurs at a disk mass of about . Such a shift would in effect change the plot of in Fig. 6 from passing through zero at to passing through zero at , about 5 times the planet mass, where would be maximum.

At later times, additional differences from the secular model are likely due to accretion of gas by planet as well as the disk tilt decay by dissipation in the simulations that are not included in the secular model. The disk viscosity can also result in the decay of the planet–disk relative tilt (Lubow & Ogilvie, 2001). The planet gains both mass and momentum from the disk as it accretes gas. Consequently, its orbit tilt becomes more aligned with the disk that is becoming more coplanar with the binary at later times. However, over the time span of the simulations, the tilt of the planet orbit relative to the disk is generally of order the initial tilt .

Fig. 11 shows results for a somewhat higher disk mass of . As a consequence of this higher disk mass, the planet orbital radius initially expands more than in the lower disk mass cases, the planet gains about an additional of its initial mass by the end of the simulation. The disk tilt decays to a nearly coplanar configuration with the binary orbital plane. Again, the planet gains both mass and momentum from the disk as it accretes gas. Consequently, its orbit tilt becomes more aligned with the disk that is nearly coplanar with the binary at later times.

The right panel of Fig. 11 plots the relative disk–planet tilt as a function of nodal phase difference . The system starts at and After the initial time, does not return to zero and is not very periodic. The system is initially librating and then transitions to circulating at late stages of the evolution when the disk has lost most of its initial mass. Essentially, the system is passing through the cycles plotted in Fig. 3 with decreasing disk mass over time, starting at a loop for high disk mass, reaching the near flat-top at the edge of libration, and breaking to circulation at low disk mass.

The disk surface density evolution is plotted in Fig. 12. The solid line plots the initial surface density distribution that is obtained by rescaling the density distribution obtained from the evolved coplanar configuration (plotted in Fig. 8) so that the initial disk mass is . The density distributions at later times show some further outward expansion along with a decrease in disk mass. The plots show that the deep disk gap about the planet that is located near is maintained over time.

## 5 Discussion

We have found in Section 2.7 that a secular resonance plays an important role in the misalignment process. Secular resonances have been previously considered in the context of disk–planet systems. Ward (1981) showed that secular resonances caused by planets could have swept across portions of the early solar system. The sweeping is due the gravitational effects of the gaseous disk, even if the planets do not migrate. The reason for the sweeping is that even a minimum-mass solar nebula can have an important influence on the relevant precession rates. The role of the disk is to modify the precession rates. As the nebula disperses, the precession rates vary, along with the resonance locations. Solid bodies, such as the terrestrial planets, can be driven into significantly eccentric and inclined orbits as these resonances pass through their orbital locations.

In another study, Lubow & Ogilvie (2001) examined the response of a gaseous disk to secular driving by two planets that are misaligned with respect to a disk. Because the precessional forcing period is long compared to the sound crossing time, the disk response is global in the form of a large scale warp. In these previous models the precessional modes that drive the resonances are due to compact objects (planets). In the secular resonance described in Section 2.7, the disk and binary companion star drive the secular resonance.

Also, unlike the case examined by Ward (1981) the angular momentum of the smallest mass object in the system, the planet in this paper, is not extremely small compared to the angular momentum of all other objects in the system. This property moderates the degree of tilt misalignment produced by the resonance considered in this paper, but also broadens its range of influence.

Xiang-Gruess & Papaloizou (2014) considered the evolution of a tilted disk with a planet in a binary and found approximate planet–disk coplanarity was maintained over the binary orbital periods that they simulated. Their system has an initial disk mass of (where is the binary mass) that comes closest to our model with disk mass at seen in Figs. 1 and 10 that shows significant misalignment over that time interval. Even higher disk mass models than considered by Xiang-Gruess & Papaloizou (2014) show substantial misalignment over this timescale (see Figs. 2 and 11). In particular, Fig. 2 shows that significant misalignment occurs to disk masses of . Their simulated disk does not have a fully cleared gap, but instead has a partial density depression at the orbit of the planet, as is seen in their Figs. 2 and 7. Our simulations begin with a planet fully embedded in a disk without a gap. During an initial adjustment phase, the disk evolves to have a clear gap near the orbit of the planet, as seen here in Fig. 12. For the parameters of this model, the standard gap opening criterion is satisfied (Lin & Papaloizou, 1986). As shown in Section 2.4, substantial misalignment does not occur for a system without a significant gap. Picogna & Marzari (2015) also found in SPH simulations that substantial misalignment occurs within planet-disk systems. Their simulations involve initial tilts relative to binary orbital plane of and , as were also simulated by Xiang-Gruess & Papaloizou (2014).

We have assumed that the binary orbit is circular, while binaries typically have eccentricities . The secular model described in Section 2.1 could be adapted to account for the binary eccentricity by modifying the coupling coefficients, while assuming the planet orbit and disk remain circular. Another effect of binary eccentricity is to truncate the nearly coplanar disk to a smaller radius in terms of binary semi-major axis (Artymowicz & Lubow, 1994). For moderate eccentricity, we may expect qualitatively similar behavior to the circular case, but do not pursue the analysis in this paper.

The results of this work have several implications for giant planet gas accretion. If a planet forms in a massive misaligned disk, it will remain coplanar with the precessing disk until it becomes massive enough to open a gap in the disk. Once a gap opens, the torque from the disk on the planet becomes weaker. Thus, the disk and the planet orbit may not remain coplanar. The planet–disk misalignment is enhanced by the effects of the secular resonance and also the effects of disk dissipation that cause disk tilt decay towards the binary orbital plane. If the planet gains substantial mass from the disk after its tilt has decayed, the planet will become more aligned with the disk and binary. Also, if terrestrial planets form at a late stage from the remains of the disk, then they may not be aligned to giant planet orbit. Instead they are likely to be more closely aligned to the binary orbital plane.

## 6 Summary

We have explored the evolution of planet–disk systems that orbit a member of a binary star system. The planet, disk, and binary interact through gravitational forces. The planet is taken to have an initial mass of 0.1% of the binary mass that is large enough mass to open a gap in the disk. The planet orbit and disk are initially coplanar and mildly inclined () relative to the orbital planet of the binary.

The planet–disk system undergoes secular oscillations. Over the course of the oscillations, the planet and disk generally have a level of misalignment that is comparable to the initial planet–disk tilt relative to the binary orbital plane for the parameters we considered. The misalignment is aided by the effects of a secular resonance and the decay of the disk tilt to the binary orbit plane. At later times, the planet orbit can evolve towards alignment with the disk, if the planet has gained a substantial amount of mass from a disk that has become nearly aligned with the binary orbital plane. This tendency toward alignment is due to the advection of disk momentum by the planet.

We determined the tilt evolution of the planet and disk by means of secular theory and SPH simulations. The secular model describes the general properties of the gravitational interactions between the planet, disk, and binary that are found SPH simulations. Since the secular model parameters were not tuned to match the simulations, the quantitative agreement is approximate.

In Section 2, we apply the secular theory to a planet–disk system in a binary for a nondissipative disk that lies external to the orbit a planet with a clearance that depends on the gap size. The tilts of the planet and disk undergo oscillations as the objects precess. One would expect that at very small disk mass, the disk and planet precess independently with substantial misalignment between the orbital plane of the planet relative to the disk plane. In addition one expects that for very high disk mass, the planet precession becomes locked to that of the disk with a smaller relative tilt. These expectations are realized in the secular model (see Fig. 2).

However, the amplitude of the relative planet–disk tilt oscillations does not vary monotonically with disk mass. As seen in Figs. 2 and 3, for small disk masses, the relative planet-disk tilt oscillation amplitude increases with disk mass and reaches a peak value at a disk mass for which the planet–disk interactions are just strong enough for them to begin to precess together in a mean sense (librate). That is, precession rate locking (in a mean sense, i.e., libration) then does not guarantee coplanarity. Just the opposite occurs. When mean precession rate locking sets in as a function of increasing disk mass, the planet–disk misalignment is largest. We attribute this effect to a secular nodal resonance driven by the disk and binary companion as described in Section 2.7. This peak planet–disk relative tilt occurs at disk masses that are several times the planet mass. The resonance is broad and enhances the misalignment for higher disk masses. Substantial relative inclinations between the planet orbit and the disk, of order the initial planet–disk tilts relative to the binary orbital plane, are possible for outer disk masses of the binary mass.

By means of SPH simulations, we analyzed this process with a dissipative (viscous) disk. The planet has an initial mass of of the binary mass. It advects mass and momentum from the disk and can migrate. The planet is initially embedded in a disk without gap. We describe the evolution after an initial disk adjustment period of discussed in Section 4. Several aspects of the SPH simulations agree with the secular model. However, there are some differences. The disk dissipation causes a decay of the disk tilt to the binary orbital plane. Consequently, the disk tilt does not rise back to its initial value as occurs in the oscillations of the secular model. Advection of disk momentum by the planet from a sufficiently large disk may cause the planet’s orbit to evolve towards alignment at later times, if it has gained substantial mass from a disk that has become nearly aligned with the binary orbital plane. At a higher disk mass, 0.6% of the binary mass, the planet–disk system evolves from nodal libration to circulation as the disk mass decreases (Fig. 11), as discussed in Section 4.

In the simulations, there is generally a substantial misalignment between the orbit of the planet and the disk that is comparable to the initial tilt of the system (relative to the binary orbital plane). The peak planet–disk misaligment occurs for a disk mass that is about 5 times the planet mass and significant misalignment extends to higher disk masses, largely due to the effects of the secular resonance. Picogna & Marzari (2015) also found in SPH simulations that substantial misalignment occurs within planet-disk systems that orbit a member of a binary.

In this work, we have considered only mild initial misalignments between the planet–disk system and the binary orbital plane. Recently, we found that substantially misaligned disks (tilts between about and with respect to the binary orbital plane) can undergo coherent Kozai–Lidov tilt and eccentricity oscillations (Martin et al., 2014a, b; Fu et al., 2015). In a future paper, we will investigate the evolution of planet–disk systems with such initial misalignments.

We thank the referee for useful comments and informing us about the paper by Picogna & Marzari (2015). We thank Jim Pringle for helpful conversations. SHL acknowledges support from NASA grant NNX11AK61G. We acknowledge the use of SPLASH (Price, 2007) for the rendering of the figures. Computer support was provided by UNLV’s National Supercomputing Center. This work also utilised the Janus supercomputer, which is supported by the National Science Foundation (award number CNS–0821794), the University of Colorado Boulder, the University of Colorado Denver, and the National Center for Atmospheric Research. The Janus supercomputer is operated by the University of Colorado Boulder.

### Footnotes

1. affiliationmark:
2. affiliationmark:
3. affiliationmark:
4. affiliationmark:

### References

1. Albrecht, S., Winn, J. N., Johnson, J. A., Howard, A. W., Marcy, G. W., Butler, R. P., Arriagada, P., Crane, J. D., & et al. 2012, ApJ, 757, 18
2. Artymowicz, P., & Lubow, S. H. 1994, ApJ, 421, 651
3. Bate, M. R. 2012, MNRAS, 419, 3115
4. Bate, M. R., Bonnell, I. A., Clarke, C. J., Lubow, S. H., Ogilvie, G. I., Pringle, J. E., & Tout, C. A. 2000, MNRAS, 317, 773
5. Bate, M. R., Lodato, G., & Pringle, J. E. 2010, MNRAS, 401, 1505
6. Bate, M. R., Lubow, S. H., Ogilvie, G. I., & Miller, K. A. 2003, MNRAS, 341, 213
7. Batygin, K. 2012, Nature, 491, 418
8. Batygin, K., Morbidelli, A., & Tsiganis, K. 2011, A&A, 533, A7
9. Borderies, N., Goldreich, P., & Tremaine, S. 1984, ApJ, 284, 429
10. Fu, W., Lubow, S. H., & Martin, R. G. 2015, ApJ, in press
11. Hale, A. 1994, AJ, 107, 306
12. Huber, D., Carter, J. A., Barbieri, M., Miglio, A., Deck, K. M., Fabrycky, D. C., Montet, B. T., Buchhave, L. A., & et al. 2013, Science, 342, 331
13. Jensen, E. L. N. & Akeson, R. 2014, Nature, 511, 567
14. Jensen, E. L. N., Mathieu, R. D., Donar, A. X., & Dullighan, A. 2004, ApJ, 600, 789
15. King, A. R., Livio, M., Lubow, S. H., & Pringle, J. E. 2013, MNRAS, 431, 2655
16. Kozai, Y. 1962, AJ, 67, 591
17. Larwood, J. D., Nelson, R. P., Papaloizou, J. C. B., & Terquem, C. 1996, MNRAS, 282, 597
18. Lidov, M. L. 1962, Planet. Space Sci., 9, 719
19. Lin, D. N. C. & Papaloizou, J. 1986, ApJ, 309, 846
20. Lodato, G. & Price, D. J. 2010, MNRAS, 405, 1212
21. Lodato, G. & Pringle, J. E. 2007, MNRAS, 381, 1287
22. Lubow, S. H. 1992, ApJ, 398, 525
23. Lubow, S. H., & D’Angelo, G. 2006, ApJ, 641, 526
24. Lubow, S. H., Martin, R. G., & Nixon, C. 2015, ApJ, 800, 96
25. Lubow, S. H. & Ogilvie, G. I. 2000, ApJ, 538, 326
26. —. 2001, ApJ, 560, 997
27. Lund, M. N., Lundkvist, M., Silva Aguirre, V., Houdek, G., Casagrande, L., Van Eylen, V., Campante, T. L., Karoff, C., & et al. 2014, A&A, 570, A54
28. Martin, R. G., Nixon, C., Armitage, P. J., Lubow, S. H., & Price, D. J. 2014a, ApJL, 790, L34
29. Martin, R. G., Nixon, C., Lubow, S. H., Armitage, P. J., Price, D. J., Doğan, S., & King, A. 2014b, ApJL, 792, L33
30. Murray, C. D., & Dermott, S. F. 1999, Solar system dynamics by Murray, C. D., 1999
31. Nixon, C., King, A., & Price, D. 2013, MNRAS, 434, 1946
32. Nixon, C., King, A., Price, D., & Frank, J. 2012, ApJl, 757, L24
33. Nixon, C. J. 2012, MNRAS, 423, 2597
34. Offner, S. S. R., Kratter, K. M., Matzner, C. D., Krumholz, M. R., & Klein, R. I. 2010, ApJ, 725, 1485
35. Paczynski, B. 1977, ApJ, 216, 822
36. Papaloizou, J. C. B. & Terquem, C. 1995, MNRAS, 274, 987
37. Picogna, G. & Marzari, F. 2015, A&A, 583, id.A133, 8 pp.
38. Price, D. J. 2007, Pasa, 24, 159
39. —. 2012, Journal of Computational Physics, 231, 759
40. Price, D. J. & Federrath, C. 2010, MNRAS, 406, 1659
41. Roccatagliata, V., Ratzka, T., Henning, T., Wolf, S., Leinert, C., & Bouwman, J. 2011, A&A, 534, A33
42. Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 24, 337
43. Skemer, A. J., Close, L. M., Hinz, P. M., Hoffmann, W. F., Kenworthy, M. A., & Miller, D. L. 2008, ApJ, 676, 1082
44. Stapelfeldt, K. R., Krist, J. E., Ménard, F., Bouvier, J., Padgett, D. L., & Burrows, C. J. 1998, ApJL, 502, L65
45. Tokuda, K., Onishi, T., Saigo, K., Kawamura, A., Fukui, Y., Matsumoto, T., Inutsuka, S.-i., Machida, M. N., & et al. 2014, ApJL, 789, L4
46. Williams, J. P., Mann, R. K., Di Francesco, J., Andrews, S. M., Hughes, A. M., Ricci, L., Bally, J., Johnstone, D., & Matthews, B. 2014, ApJ, 796, 120
47. Winn, J. N. & Fabrycky, D. C. 2014, ArXiv e-prints
48. Xiang-Gruess, M. & Papaloizou, J. C. B. 2014, MNRAS, 440, 1179
103645