PLANET FORMATION IN STELLAR BINARIES. I. PLANETESIMAL DYNAMICS IN MASSIVE PROTOPLANETARY DISKS

Results of radial velocity surveys demonstrate that ∼20% of exoplanets reside in stellar binaries (Desidera & Barbieri 2007). While most of these binaries have wide separation between stellar components (hundreds of AU), some of them are systems with relatively short binary periods of ∼100 yr. One of the best examples of such a binary is γ Cephei (Hatzes et al. 2003), which consists of two stars with masses of M_p = 1.6 M_☉ and M_s = 0.41 M_☉ with a semi-major axis of a_b = 19 AU (orbital period P_b = 58 yr) and an eccentricity of e_b = 0.41. The planet with the projected mass M_plsin i = 1.6 M_J is in orbit around the primary with a semi-major axis of a_pl ≈ 2 AU. Several more planetary systems within tight (a_b ≈ 20 AU) binaries are known at present (Chauvin et al. 2011), including the terrestrial planet around our stellar neighbor α Cen (Dumusque et al. 2012; see Hatzes 2013).

For a long time, theorists struggled to explain the origin of planets in such S-type (according to Dvorak's 1982 classification) systems. The issue lies in the strong dynamical excitation to which any object in a binary is subjected. Gravitational perturbations due to the eccentric companion are expected to adversely affect planet formation already at the stage of planetesimal growth. As first shown by Heppenheimer (1978), companion perturbations drive planetesimal eccentricities to high values, easily approaching 0.1 at 2 AU from the primary. Planetesimals would then be colliding at relative speeds of a couple km s⁻¹, which is much higher than the escape speed from the surface of even a 100 km object (about 100 m s⁻¹). As a result, collisions should lead to planetesimal destruction rather than growth.

A number of possibilities have been explored to at least alleviate this problem. In particular, Marzari & Scholl (2000) studied the dissipative effects of gas drag as the means of damping the relative velocities of planetesimals. These authors have shown that for a circular disk in secular approximation, gas drag induces an alignment of planetesimal orbits such that the periapses of small objects strongly affected by gas drag tend to cluster around 3π/2 with respect to the binary apsidal line. This was originally thought (Marzari & Scholl 2000; Thébault et al. 2004) to assist planetesimal agglomeration since the relative velocities of colliding bodies are reduced by such orbital phasing. However, it was subsequently recognized (Thébault et al. 2008) that the reduction of the relative velocity caused by apsidal alignment is effective only for planetesimals of similar sizes. Objects of different sizes still collide at high speeds, which complicates their growth.

These studies have generally arrived at the same conclusion— the difficulty of planetesimal accretion—despite the different ways in which the gas disk and its interaction with planetesimals was treated. While the early calculations (Thébault et al. 2004, 2006, 2008, 2009; Thébault 2011) typically assumed disk properties to be described by some (semi-)analytic axisymmetric models, several recent studies followed the properties and evolution of gas disks in binaries using direct hydrodynamical simulations (Paardekooper et al. 2008; Kley et al. 2008; Marzari et al. 2009a, 2012; Regály et al. 2011; Müller & Kley 2012; Picogna & Marzari 2013). One of the most important aspects of disk physics that the latter captures is the development of non-axisymmetry in the surface density distribution of the gaseous disk. This non-axisymmetry emerges under the gravitational perturbation of the binary companion, predominantly in the form of non-zero eccentricity of the fluid trajectories (Marzari et al. 2012). Another phenomenon is the disk precession, which sometimes develops in simulations and has a subsequent effect on planetesimal dynamics.

An entirely different way of lowering planetesimal eccentricities in binaries has been pursued by Rafikov (2013, hereafter R13), who demonstrated that planetesimal eccentricities can be considerably lower than previously thought by properly accounting for the gravity of a massive axisymmetric gaseous disk in which planetesimals form. The non-Keplerian potential of the disk drives rapid precession of planetesimal orbits, suppressing the companion's driving of their eccentricity.

Note that massive protoplanetary disks must have been quite natural in γ Cep-like systems since all of the known systems (with the exception of α Cen) harbor Jupiter-like planets with M_plsin i = (1.6–4) M_J (Chauvin et al. 2011). It is natural to expect the parent disk mass to exceed the planet mass by at least a factor of several (this number is very uncertain but is believed to be ∼10 for the Minimum Mass Solar Nebula), which does not make an assumption of a (0.01–0.1) M_☉ disk unreasonable despite the fact that sub-millimeter surveys find very low fluxes of thermal dust emission in young binaries with semi-major axes of several tens of AU (Harris et al. 2012).

As discussed above, the assumption of a purely axisymmetric disk may be too simplistic since simulations indicate that protoplanetary disks in binaries often develop significant eccentricities. To that effect, Silsbee & Rafikov (2015, hereafter SR15) presented the first investigation of the secular excitation of planetesimal eccentricities by the simultaneous action of the gravitational perturbations due to both the eccentric gaseous disk and the companion star. They showed that the non-axisymmetric gravitational field of such a disk excites planetesimal eccentricity (in addition to the excitation produced by the companion) and usually does not allow it to drop below the disk eccentricity, which may be rather high as suggested by some simulations (Okazaki et al. 2002; Paardekooper et al. 2008; Kley & Nelson 2008). This would again suppress planetesimal growth. On the other hand, SR15 outlined several ways in which this issue can be alleviated, for example, if the gaseous disk is precessing rapidly or if its own self-gravity is capable of reducing disk eccentricity to low levels. At the same time, SR15 did not include gas drag in their calculations, eliminating the possibility of planetesimal apsidal alignment and their eccentricity damping by drag.

Our current work builds upon the results of previous studies by exploring planetesimal dynamics in disks coplanar with the binary under the combined effects of (1) gravitational perturbations due to the eccentric gaseous disk in which planetesimals are embedded, (2) the gravity of the eccentric companion, and (3) gas drag. While we do not model disks in binaries using hydrodynamical simulations, we still capture their main features, namely their non-axisymmetry and the possibility of precession. Our results are then used in a companion paper (Rafikov & Silsbee 2015, hereafter Paper II) to explore the details of planet formation in binaries.

We thereby extend the existing semi-analytical studies in which only gas drag (and not disk gravity) was accounted for (Thébault et al. 2004, 2006, 2008, 2009; Leiva et al. 2013). We also go beyond the works of R13 and SR15 in which gas drag was neglected and only the gravitational effects of the gaseous disk and binary companion were considered. In addition, we extend the study of Beaugé et al. (2010) devoted to exploring planetesimal dynamics in eccentric, precessing disks by accounting for the gravitational potential of such a disk.

This paper is structured as follows. We discuss our general setup in Section 2 and then derive equations for the evolution of orbital elements of planetesimals in eccentric disks in Section 3. A prescription for the gas-drag-induced eccentricity evolution is described in Section 4. Solutions to the equations of planetesimal dynamics in non-precessing and precessing disks are presented in Sections 5 and 6 (as well as Appendices B and C) correspondingly. The diversity of planetesimal dynamical behavior is discussed in Section 7. We derive the relative velocity distribution of objects of different sizes in Section 8. We provide an extensive discussion of our dynamical results and their applications in Section 10. We compare different approximations for treating planetesimal dynamics in Section 10.1, briefly discuss limitations of this work in Section 10.2, and summarize our main conclusions in Section 11. Finally, some of our analytical derivations use the local approximation for treating elliptical motion, which is reviewed in Appendix A.

Our general setup is similar to that explored in SR15. We consider an elliptical disk around a primary star in a binary with semi-major axis a_b, eccentricity e_b, and component masses M_p (primary) and M_s (secondary). We define μ ≡ M_s/(M_p + M_s) and ν ≡ M_s/M_p. Binary, disk, and planetesimal orbits within it are assumed to be coplanar. This distinguishes our work from many other studies focused on the effects of Lidov–Kozai oscillations (Lidov 1962; Kozai 1962) on planetesimal dynamics in systems with inclined companions (Marzari et al. 2009a; Batygin et al. 2011; Zhao et al. 2012).

Non-axisymmetric disk structure is described via the non-zero disk eccentricity, which is a viable approximation given that simulations tend to show the prevalence of the m = 1 azimuthal harmonic of the disk shape distortion (Marzari et al. 2012). Fluid elements in a disk follow elliptical trajectories with an eccentricity e_g(a_d) that is a function of the semi-major axis a_d of a particular ellipse. All of them have the primary star of the binary as a focus. For simplicity, we assume all fluid elliptical trajectories have aligned apsidal lines, uniquely determining disk orientation via a single parameter ϖ_d—the angle between the disk and binary apsidal lines. The latter is assumed to be fixed in space as the precession of the binary under the gravity of the disk is slower than all other processes. The assumption of apsidal alignment does not affect the qualitative features of planetesimal dynamics and can be easily relaxed using the results of Statler (2001).

Because gas moves on ellipses, its surface density generally varies along the trajectory (Statler 2001; Ogilvie 2001). To obtain the gas surface density Σ(r_d, ϕ_d) at a point in a disk with polar coordinates (r_d, ϕ_d), we specify the gas surface density at periastron for each elliptical trajectory Σ_p(a_d) as a function of the semi-major axis of the corresponding ellipse a_d. SR15 show how this and the knowledge of e_g(a_d) can be used to derive Σ(r_d, ϕ_d) everywhere in the disk. In this work, following SR15, we assume a simple power-law dependence for both e_g and Σ_p:

$$\begin{eqnarray} &&\Sigma _p(a_d)=\Sigma _0\left(\frac{a_{\rm out}}{a_d}\right)^p,\quad e_g(a_d)=e_0\left(\frac{a_{\rm out}}{a_d}\right)^q, \end{eqnarray} \tag{ 1 }$$

where a_out is the semi-major axis of the outermost elliptical trajectory of the disk, and Σ₀ and e₀ are the values of Σ_p and e_g at a_out. The gravity of the companion truncates the disk at this outer radius a_out, which for eccentric binaries with e_b = 0.4 is about (0.2–0.3)a_b (Artymowicz & Lubow 1994; Regály et al. 2011). Unless stated otherwise, we will be using a_out = 5 AU in this work.

In all calculations of this paper we will be using a disk model with p = 1 and q = −1. Some motivation for singling out these particular values of p and q for circumstellar disks in binaries has been provided in R13 and SR15.

The total disk mass $M_d\approx 2\pi \int ^{a_{\rm out}}_{a_{\rm in}}\Sigma _p(a_d)a_d da_d$ enclosed within a_out can be used to express Σ_p as

$$\begin{eqnarray} \Sigma _p(a_d) & = & \frac{2-p}{2\pi }\frac{M_d}{a_{\rm out}^2} \left(\frac{a_{\rm out}}{a_d}\right)^p\nonumber \\ & \approx & 3\times 10^3 {\rm g\, cm}^{-2}M_{d,-2}a_{\rm out,5}^{-1}a_{d,1}^{-1} \end{eqnarray} \tag{ 2 }$$

(the numerical estimate is for p = 1), where M_{d, −2} ≡ M_d/(10⁻² M_☉), a_{out, 5} ≡ a_out/(5 AU) and a_{d, 1} ≡ a_d/AU. In Equation (2) we neglected disk ellipticity and assumed p < 2, so that most of the disk mass is concentrated in its outer part.

We are interested in the dynamics of planetesimals orbiting the primary within the disk and coplanar to it. We characterize their orbits with the semi-major axis a_p, the eccentricity e_p, and the apsidal angle (with regard to the binary apsidal line) ϖ_p. The orbital evolution of planetesimals is treated in a secular approximation, i.e., neglecting short-term gravitational perturbations (Murray & Dermott 1999). We also assume e_p ≪ 1 as well as e_g ≪ 1 and introduce, for convenience, the planetesimal eccentricity vector e_p = (k_p, h_p) = e_p(cos ϖ_p, sin ϖ_p).

In this work, we fully account for gravitational perturbations due to both the binary companion and the eccentric disk using the approach advanced in SR15. For the disk properties described by Equation (1), SR15 calculated an analytic expression for the planetesimal disturbing function accounting for the gravity of both the disk and the secondary. They then derived a set of Lagrange equations (see their Equations (16) and (17)) describing the evolution of e_p under the influence of gravitational forces alone.

In addition, we take into account the effects of gas drag on the secular evolution of planetesimal eccentricity. Drag-induced dissipation also results in non-conservation of energy and evolution of a_p. However, to zeroth order, we can neglect this since the radial inspiral of planetesimals usually occurs on much longer timescales than their eccentricity evolution (Adachi et al. 1976). As a result, we can concentrate on the behavior of e_p at fixed a_p and determine the relative velocities of planetesimals and their collisional outcomes.

Gas drag introduces additional terms in the eccentricity evolution equations of SR15, which we re-write in the following form:

$$\begin{eqnarray} &&\frac{dh_p}{dt}=Ak_p+B_b+B_d\cos \varpi _d(t)+\dot{h}_p^{\rm drag}, \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} &&\frac{dk_p}{dt}=-Ah_p-B_d\sin \varpi _d(t)+\dot{k}_p^{\rm drag}. \end{eqnarray} \tag{ 4 }$$

Here A = A_b + A_d is the planetesimal precession rate. It is contributed both by the gravity of the secondary (A_b) and the disk (A_d), with

$$\begin{eqnarray} A_b & = & \frac{\nu }{4}n_p \alpha _b^2 b_{3/2}^{(1)}(\alpha _b)\approx \frac{3}{4}n_p\nu \left(\frac{a_p}{a_b}\right)^3 \\ & \approx & 5.9\times 10^{-4}{\rm yr}^{-1}\nu \frac{M_{p,1}^{1/2}}{a_{b,20}}a_{p,1}^{3/2}, \nonumber \end{eqnarray} \tag{ 5 }$$

where $n_p=(GM_p/a_p^3)^{1/2}$ is the planetesimal mean rate, M_{p, 1} ≡ M_p/M_☉, a_{p, 1} ≡ a_p/AU, a_{b, 20} ≡ a_b/(20 AU), $b_s^{(j)}(\alpha)$ is the standard Laplace coefficient (Murray & Dermott 1999), α_b ≡ a_p/a_b and the approximation in Equation (5) works for α_b ≪ 1. The disk contribution is

$$\begin{eqnarray} A_d & = & 2\pi \frac{G \Sigma _p(a_p)}{a_p n_p}\psi _1= (2-p)\psi _1 n_p\frac{M_d}{M_p}\left(\frac{a_p}{a_{\rm out}} \right)^{2-p} \nonumber\\ & \approx & -6.3\times 10^{-3}{\rm yr}^{-1}a_{p,1}^{-1/2} \frac{M_{d,-2}}{M_{p,1}^{1/2}a_{\rm out,5}}, \end{eqnarray} \tag{ 6 }$$

where the numerical estimate is for p = 1 so that ψ₁ = −0.5 (SR15). Dimensionless coefficients of order unity ψ₁ and ψ₂ (see Equation (8)) have been calculated in SR15 and are functions of the disk model and the distance of the planetesimal orbit from the disk edges. One can see that for reasonable assumptions about the disk mass (M_d ∼ 10⁻² M_☉), the planetesimal precession rate at 1 AU is dominated by the disk gravity.

Eccentricity excitation by the binary (B_b) and the disk (B_d) is described by

$$\begin{eqnarray} B_b & = &{ -} \frac{\nu }{4}n_p\alpha _b^2 b_{3/2}^{(2)}(\alpha _b)e_b\approx -\frac{15}{16}n_p\nu \left(\frac{a_p}{a_b}\right)^4e_b, \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} B_d & = & \pi \frac{G \Sigma _p(a_p)}{a_p n_p}e_d(a_p)\psi _2, \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} & = & \frac{2-p}{2}\psi _2 e_g(a_p)n_p\frac{M_d}{M_p}\left(\frac{a_p}{a_{\rm out}}\right)^{2-p}, \end{eqnarray} \tag{ 9 }$$

with the latter explicitly depending on the local value of the disk eccentricity e_g(a_p).

Note that ϖ_d in Equations (3) and (4) is not necessarily constant—it can be an explicit function of time, allowing one to treat the case of a precessing disk.

The terms $\dot{h}_p^{\rm drag}$ and $\dot{k}_p^{\rm drag}$, absent in the original version of Equations (3) and (4) in SR15, represent the effect of gas drag on the eccentricity evolution; they are derived in Section 4. The main goal of this work is to see how their introduction affects planetesimal dynamics.

Next we derive the expressions for the drag-induced eccentricity evolution terms $\dot{h}_p^{\rm drag}$ and $\dot{k}_p^{\rm drag}$ applicable to the case of an eccentric disk.

Because of our assumption of small eccentricities for both gas and planetesimals, it is reasonable to employ the local (or guiding center) approximation. This approach is often used in studies of planetesimal and galactic dynamics (Binney & Tremaine 2008) and forms a basis of the so-called Hill approximation (Hénon & Petit 1986; Hasegawa & Nakazawa 1990). Local approximation considers planetesimal motion in a local Cartesian x−y reference system aligned with the radial and azimuthal directions, respectively. The main features of this approximation are reviewed in Appendix A. In particular, Equations (A8) describe how k_p and h_p evolve under the effect of external force F.

In our case, F is the drag force arising because of the motion of planetesimals with respect to gas. Adachi et al. (1976) gives the following expression for quadratic drag force appropriate for rapidly moving objects with sizes larger than the mean free path of gas molecules:

$$\begin{eqnarray} &&{\bf F}=-\frac{C_D}{2}\pi d_p^2\rho _g v_r{\bf v}_r, \end{eqnarray} \tag{ 10 }$$

where C_D is a constant drag coefficient taken to be 0.5 throughout this paper, d_p is the particle size, and ρ_g is the gas density. The relative particle–gas velocity v_r is given by Equations (A5)–(A7) with relative particle–gas eccentricity components

$$\begin{eqnarray} &&h_r=h_p-h_g,\quad k_r=k_p-k_g, \end{eqnarray} \tag{ 11 }$$

and e_g = (k_g, h_g) = e_g(cos ϖ_d, sin ϖ_g) being the local value of the gas eccentricity vector. Using these expressions, we obtain the force components F_x and F_y:

$$\begin{eqnarray} F_x & =& {-}\frac{3C_D}{8}m_p D v_r^a (k_r\sin n_p t - h_r\cos n_p t), \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray} F_y & =&{ -}\frac{3C_D}{16}m_p D v_r^a (k_r\cos n_p t + h_r\sin n_p t), \end{eqnarray} \tag{ 13 }$$

where m_p is the planetesimal mass, and the relative velocity $v_r^a$ is given by Equation (A7). The prefactor D is given by

$$\begin{eqnarray} &&D=n_p\frac{\Sigma _g}{\rho _p d_p}\frac{r}{h}, \end{eqnarray} \tag{ 14 }$$

with ρ_p being the particle bulk density and h = c_s/n_p being the disk scale height (c_s = (kT_g/μ)^1/2).

Now we plug the expressions for F_x, F_y into the first two Equations (A8) and then average them in time t over the planetesimal orbital period (this is the secular, i.e., time-averaged approximation). One can easily see that to get the result to lowest-order in e_r we do not need to keep the terms O(e_r, e_d, e_p) in the expression for ρ_g. As a result, we find

$$\begin{eqnarray} && \dot{k}_p^{\rm drag}=-\frac{3C_D}{4\pi }{\rm E}\left(\frac{\sqrt{3}}{2} \right)D k_r e_r, \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray} && \dot{h}_p^{\rm drag}=-\frac{3C_D}{4\pi }{\rm E}\left(\frac{\sqrt{3}}{2} \right)D h_r e_r, \end{eqnarray} \tag{ 16 }$$

where $E(\sqrt{3}/2)\approx 1.211$ is a complete elliptic integral, and $e_r^2=k_r^2+h_r^2$.

We can rewrite Equations (15) and (16) in the following form:

$$\begin{eqnarray} &&\dot{k}_p^{\rm drag}=-\frac{k_p-k_g}{\tau _d},\quad \dot{h}_p^{\rm drag}=-\frac{h_p-h_g}{\tau _d}, \end{eqnarray} \tag{ 17 }$$

where the eccentricity damping time

$$\begin{eqnarray} \tau _d & = & \frac{4\pi }{3C_D{\rm E}\left(\sqrt{3}/2\right)}D^{-1}e_r^{-1} \nonumber\\ & \approx & 600 {\rm yr } C_D^{-1} \frac{a_{\rm out,5}a_{p,1}}{M_{p,1}^{1/2}M_{d,-2}} \frac{h/r}{0.1}\frac{10^{-2}}{e_r}d_{p,1}. \end{eqnarray} \tag{ 18 }$$

Here d_{p, 1} ≡ d_p/(1 km) and numerical estimate is for p = 1 and ρ_p = 3 g cm⁻³; in the case of quadratic drag law (10) τ_d depends on k_p and h_p through e_r; see Equations (11).

Results of Section 4 allow us to understand the behavior of e_p. For simplicity, we start by considering the case of a non-precessing disk, i.e., ϖ_d =const. Even in this case Equations (3) and (4) with the quadratic drag terms (15) and (16) cannot generally be solved analytically because of the τ_d dependence on e_r.

However, it can be easily shown that solutions of these equations inevitably converge to a steady-state form—the free eccentricity, which depends on initial conditions (R13; SR15), damps out and e_p converges to the forced eccentricity vector (Beaugé et al 2010). This is illustrated in Figure 1 where we solve evolution equations numerically. It is clear that starting with arbitrary initial conditions and after initial (sometimes oscillatory) evolution, k_p and h_p do converge to the same steady state values (depending only on the disk parameters and planetesimal size), which are given by Equations (22) and (23) derived below. This point is additionally illustrated in Figure 2(a) where we plot the trajectory of e_p as it evolves in h_r–k_r coordinates. There one can clearly see e_p converging to a fixed point solution, in an oscillatory fashion for large planetesimals, and exponentially for small objects, which rapidly couple to the gas disk.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Damping of the memory of initial conditions can also be demonstrated by solving Equations (3) and (4) analytically in the simplified but qualitatively similar case of a linear drag law, when τ_d is independent of h_p and k_p. Such a solution is presented in Appendix B for the general case of a precessing gaseous disk. The non-precessing disk solution is obtained by setting $\dot{\varpi }_d=0$. This clearly demonstrates the convergence of e_p to a time-independent, forced value.

SR15 have demonstrated that in the absence gas drag, under the action of only the gravity of the disk and the companion star, the steady state (forced) eccentricity is given by

$$\begin{eqnarray} {\bf e}_p^{\rm n/drag} & = & \left\lbrace \begin{array}{@{}l@{}}k_p^{\rm n/drag}\\ h_p^{\rm n/drag} \end{array} \right\rbrace = {\bf e}_{b}+{\bf e}_{d}, \end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray} {\bf e}_{b} & = & \left\lbrace \begin{array}{@{}l@{}}k_b\\ h_b \end{array} \right\rbrace = -\frac{B_b}{A} \left\lbrace \begin{array}{@{}l@{}}1 \\ 0 \end{array} \right\rbrace, \end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray} {\bf e}_{d} & = & \left\lbrace \begin{array}{@{}l@{}}k_d\\ h_d \end{array} \right\rbrace = -\frac{B_d}{A} \left\lbrace \begin{array}{@{}l@{}}\cos \varpi _d \\ \sin \varpi _d \end{array} \right\rbrace, \end{eqnarray} \tag{ 21 }$$

where e_b and e_d are forced eccentricity vectors due to the secondary and disk gravity, respectively. Note that the accuracy of the analytical expression (20) for the binary contribution is known to worsen (beyond the ∼10% level) when a_p/a_b ≳ 0.1 (Thébault et al. 2006; Barnes & Greenberg 2006). More refined calculations of e_b are possible (Veras & Armitage 2007; Giuppone et al. 2011) but for the purposes of this work it is sufficient to use Equation (20).

With the gas drag included, the behavior of e_p changes. To determine the steady-state values of k_p and h_p and analyze their properties, we use the prescription (17), set to zero time derivatives in the left-hand sides of Equations (3) and (4) and solve the resulting algebraic system with respect to h_p and k_p. We find as a result

$$\begin{eqnarray} k_p & = & k_b+k_d + \frac{ (k_g-k_b-k_d) - (h_g-h_d)(A\tau _d) }{1+(A\tau _d)^2}, \end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray} h_p & = & h_d+ \frac{ (h_g - h_d) + (k_g-k_b-k_d)(A\tau _d)}{1+(A\tau _d)^2}, \end{eqnarray} \tag{ 23 }$$

where k_b, k_d, h_b, h_d are defined in Equations (20) and (21). These asymptotic results are valid even if τ_d is a function of e_r —in that case they simply represent two implicit relations for h_p and k_p.

Solutions (22) and (23) can be re-written in vectorial form as

$$\begin{eqnarray} {\bf e}_p & = & \left\lbrace \begin{array}{@{}l@{}}k_p\\ h_p \end{array} \right\rbrace ={\bf e}_{f,b}+{\bf e}_{f,d}, \end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray} {\bf e}_{f,b} & = & k_b\frac{(A\tau _d)}{1+(A\tau _d)^2} \left\lbrace \begin{array}{@{}l@{}}(A\tau _d)\\ -1 \end{array} \right\rbrace, \end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray} {\bf e}_{f,d} & = & \left[\frac{e_g^2+\tau _d^2 B_d^2}{1+(A\tau _d)^2}\right]^{1/2} \left\lbrace \begin{array}{@{}l@{}}\cos \left(\varpi _d+\phi \right)\\ \sin \left(\varpi _d+\phi \right) \end{array} \right\rbrace, \end{eqnarray} \tag{ 26 }$$

where the phase shift ϕ is given by

$$\begin{eqnarray} &&\cos \phi =\frac{e_g- A B_d\tau _d^2}{\left(e_g^2+\tau _d^2 B_d^2\right)^{1/2} \left[1+(A\tau _d)^2 \right]^{1/2}}. \end{eqnarray} \tag{ 27 }$$

In the limit of vanishing drag, Aτ_d → ∞, one finds ϕ → π and the solution (24)–(26) reduces to the non-drag result with no free eccentricity (19)–(21); see SR15.

In the limit of strong drag (Aτ_d → 0) in a circular disk (i.e., e_g = 0) and no disk gravity (i.e., A_d = B_d = 0) one finds h_p/k_p → −∞. This means that in this case planetesimal apsidal lines cluster around ϖ_p = 3π/2, in agreement with Marzari & Scholl (2000). Also, |e_p| → B_bτ_d directly depends on planetesimal size, which implies that in this limit, planetesimals of different sizes collide with non-zero speeds even despite their apsidal alignment (Thébault et al. 2008).

Expressions (24)–(26) clearly show that e_p can be split into two distinct components: a contribution e_{f, b} due to the gravity of the binary and a contribution e_{f, b} related to both the gravitational and gas drag effects of the disk. It is also clear that after reaching steady state, planetesimal orbits are generally aligned with neither the disk (ϖ_p ≠ ϖ_d) nor the binary (ϖ_p ≠ 0).

5.1. Relative Particle–Gas Eccentricity

In the case of quadratic drag (10), we can further analyze eccentricity behavior. Using Equations (22) and (23) we express relative particle–gas eccentricity as

$$\begin{eqnarray} &&e_r=|{\bf e}_p - {\bf e}_g|= e_c\frac{\left(A\tau _d\right)}{\sqrt{1+\left(A\tau _d\right)^2}}, \end{eqnarray} \tag{ 28 }$$

where we introduced a characteristic eccentricity $e_c=|{\bf e}_p^{\rm n/drag} - {\bf e}_g|$ given by

$$\begin{eqnarray} e_c & \equiv & \left[(h_g-h_d)^2+ (k_g-k_b-k_d)^2\right]^{1/2} \nonumber \\ & = & \frac{\left[(Ae_g+B_d)^2+B_b^2+2\cos \varpi _dB_b(Ae_g+B_d)\right]^{1/2}}{|A|}.\nonumber \\ \end{eqnarray} \tag{ 29 }$$

Plugging this expression for e_r into Equation (18), one obtains the following bi-quadratic equation for (Aτ_d):

$$\begin{eqnarray} &&\left(A\tau _d\right)^4=\left(\frac{d_p}{d_c}\right)^2 \left[\left(A\tau _d\right)^2+1\right], \end{eqnarray} \tag{ 30 }$$

where we have introduced a characteristic planetesimal size d_c defined as

$$\begin{eqnarray} &&d_c \equiv \frac{3C_D{\rm E} (\sqrt{3}/2)}{4\pi }\frac{n_p}{|A|}\frac{\Sigma _g}{\rho _p} \frac{r}{h}e_c. \end{eqnarray} \tag{ 31 }$$

All our subsequent results can be formulated completely in terms of e_c and d_p/d_c, underscoring the significance of these variables. A detailed discussion of the characteristic values and general behavior of e_c and d_c is provided in Sections 7.1 and 7.3.

Solving Equation (30), one finds

$$\begin{eqnarray} &&\left|A\tau _d\right|=\frac{d_p}{d_c}\left[\frac{1}{2}+ \sqrt{\frac{1}{4}+\left(\frac{d_c}{d_p}\right)^2}\right]^{1/2}, \end{eqnarray} \tag{ 32 }$$

i.e., that |Aτ_d| is a function of d_r/d_c only.

Plugging Equation (32) into Equation (28), one also finds the general expression for the relative particle–gas eccentricity

$$\begin{eqnarray} &&e_r=e_c\frac{d_p}{d_c} \left[\sqrt{\frac{1}{4}+\left(\frac{d_c}{d_p}\right)^2} -\frac{1}{2}\right]^{1/2} \end{eqnarray} \tag{ 33 }$$

valid for arbitrary d_p/d_c.

We illustrate the behaviors of |Aτ_d| and e_r given by Equations (32) and (33) in Figure 3. This reveals the meaning of the characteristic size d_c: objects with d_p ∼ d_c have |Aτ_d| ∼ 1, i.e., their stopping time due to gas drag is comparable to their orbital precession period and their relative eccentricity with respect to gas is e_r ∼ e_c.

Zoom In Zoom Out Reset image size

It is instructive to further explore general solutions (32), (33) valid for arbitrary d_p/d_c in the two limits covered next.

5.2. Small Objects, d_p ≲ d_c — Strong Drag (|Aτ_d|  ≲  1)

In the limit of strong gas drag we expect the damping time τ_d to be very short and |Aτ_d| ≪ 1, so that the gas-particle velocity differential is rapidly reduced to zero. According to Equation (32), this regime is valid for small objects with d_p ≲ d_c, when

$$\begin{eqnarray} &&\left|A\tau _d\right|\approx \left(d_p/d_c\right)^{1/2}\lesssim 1. \end{eqnarray} \tag{ 34 }$$

From Equation (28), the relative particle–gas eccentricity is

$$\begin{eqnarray} &&e_r \approx \left|A\tau _d\right| e_c\approx e_c\left(d_p/d_c\right)^{1/2} \end{eqnarray} \tag{ 35 }$$

to leading order in (Aτ_d).

Equations (22) and (23) become

$$\begin{eqnarray} k_p &\rightarrow& k_g + [(h_d-h_g)(A\tau _d) \nonumber \\ & & - (k_g-k_b-k_d)(A\tau _d)^2 ], \end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray} h_p & \rightarrow& h_g + [(k_g-k_b-k_d)(A\tau _d) \nonumber \\ && + (h_d-h_g)(A\tau _d)^2 ]. \end{eqnarray} \tag{ 37 }$$

Here brackets encompass the leading order subdominant terms, compared to the zeroth order terms outside brackets.

It is clear from these asymptotic expressions that in the case of strong drag, the eccentricity vector of planetesimals tends toward the eccentricity vector of the gas, e_p → e_g. It is only weakly sensitive to gravitational perturbations due to either the companion or the disk. Thus, to leading order, the value of the eccentricity vector is independent of particle size (which enters only through τ_d).

5.3. Big Objects, d_p ≳ d_c — Weak Drag (|Aτ_d|  ≳  1)

In the opposite limit of weak drag or long damping time |Aτ_d| ≫ 1 valid for large objects with d_p ≳ d_c, Equation (32) yields

$$\begin{eqnarray} &&\left(A\tau _d\right)\approx d_p/d_c\gtrsim 1, \end{eqnarray} \tag{ 38 }$$

while the relative particle–gas eccentricity is

$$\begin{eqnarray} &&e_r \approx e_c; \end{eqnarray} \tag{ 39 }$$

see Equation (28). Thus, in the weak drag regime, e_r saturates at the value independent of the size of the object.

Equations (22) and (23) reduce in this limit to

$$\begin{eqnarray} k_p &\rightarrow &k_b+k_d + [(h_d-h_g)(A\tau _d)^{-1} \nonumber \\ & & + (k_g-k_b-k_d)(A\tau _d)^{-2} ], \end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray} h_p& \rightarrow & h_g + [(k_g-k_b-k_d)(A\tau _d)^{-1} \nonumber \\ && - (h_d-h_g)(A\tau _d)^{-2} ]. \end{eqnarray} \tag{ 41 }$$

Again, terms in brackets are subdominant compared to the leading terms (outside brackets).

This solution shows that in the limit of weak drag ${\bf e}_p\rightarrow {\bf e}_p^{\rm n/drag}$, i.e., the behavior of the particle eccentricity vector is determined predominantly by the gravitational effects of the secondary and the disk. Thus, e_p is again almost independent of the particle size.

So far we have assumed the orientation of the disk to be fixed in the binary frame. However, some simulations find disks in binaries to precess (e.g., Marzari et al 2009b; Müller & Kley 2012). We now study how planetesimal dynamics change in the case of a disk uniformly precessing at a constant rate of $\dot{\varpi }_d$. Figure 2(b) displays the evolution of e_p for the same parameters as in panel (a) of that figure, but in a disk precessing at the rate $\dot{\varpi }_d=A$. One can see that the main difference compared to the non-precessing case is that in the long run e_p converged to limit cycle behavior (Beaugé et al. 2010) rather than to a fixed point, as in panel (a). The sizes and shapes of the asymptotic limit cycles depend on both the planetesimal size d_p and the disk precession rate $\dot{\varpi }_d$, as discussed in detail in Appendix C and shown in Figure 9. This certainly complicates planetesimal dynamics.

To gain additional insights, in Appendix B, we derive a full time-dependent solution for e_p in a precessing disk for the case of linear gas drag, when τ_d is independent of the relative particle–gas eccentricity e_r. This solution fully accounts for the gravitational and gas drag effects of the precessing disk as well as for the gravity of the binary companion.

We use this solution as a basis for understanding planetesimal dynamics in a precessing disk in the more complicated but realistic case of quadratic gas drag. This regime, which does not admit a general analytical solution even for the long-term behavior, is explored in Appendix C. There we show that planetesimal dynamics with drag law (10) depend on the relative role played by the binary companion, as described in the next section.

6.1. Strong Binary Perturbation Case

The results in Appendices B and C show that whenever binary gravity dominates e_p excitation and the condition

$$\begin{eqnarray} &&|(A-\dot{\varpi }_d)e_g+B_d |\lesssim |B_b| \end{eqnarray} \tag{ 42 }$$

is fulfilled, planetesimal dynamics proceed as if the disk were not precessing: neither the gas eccentricity e_g nor the eccentricity driven by disk gravity e_d, Equation (21), is significant compared to the forced eccentricity due to the binary e_b = B_b/A (note that both the binary and disk gravity contribute to A).

In this case, e_p is close to the relative planetesimal-gas eccentricity e_r and is approximately constant. As a result, the planetesimal orbit maintains a roughly fixed orientation with respect to the binary orbit and

$$\begin{eqnarray} &&k_p\approx k_b\frac{\left(A\tau _d\right)^2}{1+\left(A\tau _d\right)^2},\quad h_p\approx -k_b\frac{\left(A\tau _d\right)}{1+\left(A\tau _d\right)^2}, \end{eqnarray} \tag{ 43 }$$

with k_b defined by Equation (20). Planetesimal orbits are aligned with the binary (ϖ_p → 0) for |Aτ_d| → ∞ (weak drag), but in the case of strong drag, |Aτ_d| → 0, the planetesimal apsidal line points at ϖ_d = 270°, which agrees with Marzari & Scholl (2000) despite the disk precession.

Interestingly, even though gas eccentricity e_g does not appear in these expressions (and neither does the precession rate $\dot{\varpi }_d$, at the lowest order), the effect of the gas drag is explicitly present via the non-trivial τ_d dependence. Thus, our precessing disk results obtained in the limit (42) apply equally well to planetesimal dynamics in a purely axisymmetric (e_g = 0) gaseous disk, extending the results of R13 to the case of non-zero gas drag—note that A in Equations (43) and in the definition of k_b is the full precession rate due to both the binary and the disk.

The value of e_r in the regime (42) is given by Equations (28) and (33) with d_c and |Aτ_d| computed using e_c ≈ e_b = |B_b/A| (i.e., Equation (29) in the limit B_d → 0, e_g → 0); see Equations (31) and (32).

6.2. Weak Binary Perturbation Case

In the opposite case of weak driving of e_p by the binary companion, we combine solutions (C5) and find the relative particle–gas eccentricity to be

$$\begin{eqnarray} &&e_r=\left|e_c^{\rm pr}\right|\frac{\left|A-\dot{\varpi }_d\right|\tau _d}{\sqrt{1+\left(A-\dot{\varpi }_d\right)^2\tau _d^2}}, \end{eqnarray} \tag{ 44 }$$

replacing Equation (28) in the case of the precessing disk. Here we defined the characteristic eccentricity as

$$\begin{eqnarray} &&e_c^{\rm pr}=-\frac{B_d}{A-\dot{\varpi }_d}-e_g, \end{eqnarray} \tag{ 45 }$$

which, according to SR15, is the relative particle–gas forced eccentricity in the no drag (τ_d → ∞) and no binary (B_b → 0) case. As $\dot{\varpi }_d\rightarrow 0$, one finds $|e_c^{\rm pr}|\rightarrow e_c$ given by Equation (29) with k_b = h_b = 0; also, Equation (44) reduces to the non-precessing disk result (28).

Plugging this expression for e_r into Equation (18) one finds

$$\begin{eqnarray} &&\left|A-\dot{\varpi }_d\right|\tau _d= \frac{d_p}{d_c^{\rm pr}}\left[\frac{1}{2}+ \sqrt{\frac{1}{4}+\left(\frac{d_c^{\rm pr}}{d_p}\right)^2}\right]^{1/2}, \end{eqnarray} \tag{ 46 }$$

with a new characteristic planetesimal size

$$\begin{eqnarray} &&d_c^{\rm pr} \equiv \frac{3C_D{\rm E} (\sqrt{3}/2)}{4\pi }\frac{n_p}{\left|A-\dot{\varpi }_d\right|}\frac{\Sigma _g}{\rho _p} \frac{r}{h}\left|e_c^{\rm pr}\right|. \end{eqnarray} \tag{ 47 }$$

These expressions are different from Equations (31) and (32) in that they use $A-\dot{\varpi }_d$ instead of A and $|e_c^{\rm pr}|$ instead of e_c. It is then clear that whenever a precessing disk dominates the planetesimal dynamics Equation (33) also holds provided that we replace $d_c\rightarrow d_c^{\rm pr}$ and $e_c\rightarrow |e_c^{\rm pr}|$. The same is true for our asymptotic results on e_p behavior presented in Sections 5.2 and 5.3 if we also take k_b → 0.

In the limit $\dot{\varpi }_d\rightarrow 0$, the value of e_f reduces to e_{f, d} given by Equation (26), but when $|\dot{\varpi }_d|\gg |A|$, rapid disk precession suppresses excitation of planetesimal eccentricity by the disk gravity, i.e., the first term in Equation (45).

It is worth noting that the results of Appendix B for the case of linear drag suggest that neglecting binary gravity in the case of a precessing disk might require a condition different from the direct opposite to the constraint (42). Indeed, the asymptotic solution (B5) for the relative eccentricity of planetesimals in the case of weak drag ($\tau _{d,1},\tau _{d,2}\gg |A-\dot{\varpi }_d|^{-1}$) shows that the term proportional to k_b can be neglected only when

$$\begin{eqnarray} &&|(A-\dot{\varpi }_d)e_g+B_d|\gtrsim |B_b| \left(\frac{A-\dot{\varpi }_d}{A}\right)^2, \end{eqnarray} \tag{ 48 }$$

which is a more stringent criterion whenever $|\dot{\varpi }_d|\gg |A|$. The same constraint may be needed in the case of quadratic drag. However, in practice, one often finds $|\dot{\varpi }_d|\lesssim |A|$; see Paper II in which case Equation (48) is just the opposite of the condition (42).

Results of Section 5 demonstrate that the steady state value of the eccentricity vector e_p is fully determined by just two key parameters—the characteristic eccentricity e_c and the critical planetesimal size d_c; see Equation (33). The eccentricity e_c sets the overall scale of e_p, while d_c is the planetesimal size at which planetesimal coupling to gas changes from weak to strong. We now explore the behavior of these variables as a function of system parameters to elucidate some important features of planetesimal dynamics.

7.1. Behavior of e_c

In Figures 4(a) and (b) we show e_c computed for the γ Cep system at 2 AU—the semi-major axis of its planet—as a function of disk mass M_d and eccentricity e₀, for two disk orientations—aligned (ϖ_d = 0) and anti-aligned (ϖ_d = π) with the apsidal line of the binary.

Zoom In Zoom Out Reset image size

One can immediately see a feature common to both panels—a narrow valley of high e_c (white because of saturation at high e_c) at almost constant M_d. It appears because at this value of disk mass, A_d = −A_b and A = 0, giving rise to a secular resonance. According to Equations (19)–(21) and (29), e_c gets driven to high values as A → 0. This resonance has been previously discussed in R13 and SR15.

Equations (5) and (6) predict that at a given distance from the primary a_p, this resonance occurs for the disk mass

$$\begin{eqnarray} M_{d,A=0} &=& M_s\frac{3}{4(2-p)|\psi _1|} \left(\frac{a_p}{a_b}\right)^{1+p} \left(\frac{a_{\rm out}}{a_b}\right)^{2-p}\nonumber \\ \ms & \approx & 1.5\times 10^{-3}\,M_\odot \frac{M_s}{0.4 \,M_\odot } \frac{a_{\rm out,5}}{a_{b,20}^3}a_{p,2}^2, \end{eqnarray} \tag{ 49 }$$

where a_{out, 5} ≡ a_out/(5 AU), a_{p, 2} ≡ a_p/(2 AU), and a_{b, 20} ≡ a_b/(20 AU). This estimate agrees with Figures 4(a) and (b) for the γ Cep parameters and a disk with p = 1 and ψ₁(p = 1) = −0.5 (SR15).

The existence of this resonance is independent of the relative disk–binary orientation because planetesimal precession rates A_b and A_d are determined by the axisymmetric components of the binary and disk gravitational potentials. For this reason, M_{d, A = 0} is the same for all disk orientations. To the right of the secular resonance, disk gravity dominates the planetesimal precession rate and suppresses e_c if the disk eccentricity is small (R13).

At high disk eccentricity, typically e₀ ≳ 0.05, this suppression vanishes because for large M_d ≳ 10⁻³ M_☉, disk gravity starts to dominate e_p excitation. This statement is true above the blue line |B_b| = |B_d| in Figures 4(a) and (b) (the origin of the low-e_c band at small M_d and high e₀ in Figure 4(a) is discussed in Section 7.2). Further increase of the disk mass in this region does not affect e_c because the planetesimal dynamics switch to the so-called DD regime (SR15) in which e_p(a_p) ≈ |ψ₂/ψ₁|e_g(a_p), independent of M_d. As a result, high e_g leads to high e_p.

In Figures 4(c) and (d) we explore the dependence of e_c on the distance from the binary a_p and M_d for two different values of the disk eccentricity e₀ = 0.1 and 0.01. Here we look only at an aligned disk case. Again, an obvious feature of these maps is the secular resonance around the blue dashed curve for |A_d| = |A_b|, where e_c is very large and collisional growth is impossible. In Figure 4(d) there is also a "valley" of low e_c to the right of the blue line |B_b| = |B_d|, whose origin is discussed in Section 7.2.

These maps make it clear that e_c becomes independent of M_d (at a given separation a_p) when the disk mass becomes large enough. This is a direct consequence of the planetesimal dynamics switching into the DD regime (SR15) when both the eccentricity excitation and apsidal precession of planetesimals are dominated by the disk gravity with a negligible contribution from the binary companion. In the high-M_d regime, e_c decreases as a_p decreases. This is a consequence of our adopted disk model, in which e_d∝a_p and the fact that e_c∝e_d in the DD regime.

7.2. Valley of Stability in Aligned Disks

Figures 4(a) and (b) show that irrespective of the disk orientation, e_c is low for high M_d ≳ 10⁻² M_☉ and small disk eccentricity, e₀ ≲ 10⁻². Outside this corner of phase space, e_c is much higher, which makes planetesimal growth problematic there. At the same time, in the case of an aligned disk (ϖ_d = 0), low values of e_c are also possible in a narrow "valley" stretching toward high e₀ and low M_d. Since this feature may have interesting implications for planet formation in binaries (see Paper II for details), we discuss its origin in more detail.

Equation (29) implies that in an aligned disk h_g = h_d = 0 so that

$$\begin{eqnarray} &&e_c\approx |k_g-k_d-k_b|=\left|\frac{B_b+B_d+A e_g}{A}\right|. \end{eqnarray} \tag{ 50 }$$

For massive disks, to the right of the vertical |A_b| = |A_d| line in Figure 4(a), one can set A ≈ A_d and relate it to B_d via Equations (6) and (8). As a result, Equation (50) becomes

$$\begin{eqnarray} &&e_c\approx \left|\frac{B_b+B_d\left(1+2\psi _1\psi _2^{-1} \right)}{A_d}\right|. \end{eqnarray} \tag{ 51 }$$

For the disk model considered here (p = 1, q = −1), one has ψ₁ = −0.5, ψ₂ = 1.5 and $1+2\psi _1\psi _2^{-1}=1/3$ so that e_c ≈ |A_d|⁻¹|B_b + B_d/3|. Also B_d > 0 while B_b is always negative; see Equations (7) and (8). Given that B_d∝e₀M_d it is then obvious that one can make e_c ≈ 0 by choosing e₀M_d such that |B_b| ≈ |B_d|/3. Thus, in the case of an aligned disk, a "valley" of low e_c is described by the relation $e_0\propto M_d^{-1}$ as long as |A_d| ≳ |A_b| (i.e., for massive disks).

From this discussion we see that e_c ≈ 0 for values of e₀ and M_d that are close to the curve

$$\begin{eqnarray} M_{d,|B_b|=|B_d|} &=& M_s\frac{15}{8(2-p)|\psi _2|} \frac{e_b}{e_0}\nonumber \\ \nonumber &\times & \left(\frac{a_p}{a_b}\right)^{2+p+q} \left(\frac{a_{\rm out}}{a_b}\right)^{2-p-q} \nonumber\\ & \approx & 1.2\times 10^{-3}\,M_\odot \frac{M_s}{0.4 \,M_\odot } \frac{e_b}{0.4}\frac{0.1}{e_0} \frac{a_{\rm out,5}^2}{a_{b,20}^4}a_{p,2}^2, \nonumber\\ \end{eqnarray} \tag{ 52 }$$

on which |B_b| = |B_d|; see Equations (7) and (8) in which we took p = 1, q = −1. This relation is shown by the blue line in Figure 4 and is quite close to the valley of low e_c.

Note that according to Equation (51), the value of e_c can be lowered globally in a massive disk if its structure is such that $1+2\psi _1\psi _2^{-1}=0$. However, this is not the case for the disk model used in this work.

The situation is different for the low mass, aligned disks to the left of the |A_d| = |A_b| (blue dashed) line in Figure 4(a). Here A ≈ A_b and B_b dominates over B_d for low enough M_d at a fixed e₀, which, in the terminology of SR15, corresponds to Case BB of planetesimal excitation. In this regime, Equation (50) shows that

$$\begin{eqnarray} &&e_c\rightarrow \left|e_g+\frac{B_b}{A_b}\right|= \left|e_g-\frac{5}{4}\frac{a_p}{a_b}e_b\right| \end{eqnarray} \tag{ 53 }$$

Our adopted radial scaling of e_g in the form Equation (1) with q = −1 results in a particular value of

$$\begin{eqnarray} &&e_0\left |_{e_c\rightarrow 0}=\frac{5}{4}\frac{a_{\rm out}}{a_b}e_b= 0.125\frac{a_{\rm out}/a_b}{0.25}\frac{e_b}{0.4}\right., \end{eqnarray} \tag{ 54 }$$

for which e_c → 0. This critical value of e₀ is independent of M_d explaining why the valley of low e_c starts going almost horizontally for M_d ≲ M_d|_{A = 0} in Figure 4(a).

Moreover, $e_0|_{e_c\rightarrow 0}$ is also independent of a_p, which means that e_c → 0 globally when $e_0\rightarrow e_0|_{e_c\rightarrow 0}$ in parts of the disk where |A_b| ≳ |A_d| and |B_b| ≳ |B_d|. This is the reason why in the upper left corner of Figure 4(c) e_c is considerably lower than in the same region of Figure 4(d), despite e₀ being an order of magnitude higher in the former case. Indeed, according to Equation (54), e₀ = 0.1 used in Figure 4(c) is very close to $e_0|_{e_c\rightarrow 0}$ for the adopted system parameters. As a result of this coincidence, e_c is strongly suppressed in the BB regime in a rather eccentric (e₀ = 0.1) disk.

A narrow region of low e_c stretching along the blue curve |B_b| = |B_d| in Figures 4(c) and (d) is the same valley of stability, but now revealing itself in M_d − a_p coordinates.¹ It may lie inside (for low e₀) as well as outside (for high e₀) of the secular resonance. Note that in Figure 4(c) the |A_d| = |A_b| and |B_d| = |B_b| curves fall almost on top of each other, which is a coincidence caused by our choice of e₀ = 0.1 in this case. Because of that, the valley of stability appears as a very narrow band of low e_c just to the left of the |B_d| = |B_b| curve in this panel.

If the disk is not aligned with the binary orbit and ϖ_d is not small, then both h_d and h_g are nonzero and contribute to e_c; see Equation (29). Moreover, for disks that are close to being anti-aligned with the binary, k_b and k_d have the same sign, eliminating the possibility of their mutual cancellation. As a result, the low-e_c valley at high e₀ and low M_d disappears as long as |ϖ_d − ϖ_b| ≳ 10°.

To summarize, the valley of stability creates favorable conditions for locally lowering planetesimal velocity in aligned disks around some particular locations even in low-mass disks with M_d ≲ 10⁻² M_☉.

7.3. Behavior of d_c

Next we discuss the behavior of the characteristic size d_c at which planetesimals of similar (but not equal) mass collide at the highest relative velocity ∼e_cv_K. Equation (31) shows that for a given value of e_c, critical size d_c scales inversely with the planetesimal precession rate|A|. If planetesimal precession is dominated by the potential of the secondary, then A = A_b, and one finds

$$\begin{eqnarray} d_c &=& \frac{C_D {\rm E}(\sqrt{3}/2)}{\pi \nu }\frac{r}{h} \frac{\Sigma _g}{\rho _p}\left(\frac{a_b}{a_p}\right)^3e_c \nonumber\\ \ms &\approx & 30 \; {\rm km} \frac{C_D}{\nu }\frac{0.1}{h/r} \frac{M_{d,-2}a_{b,20}^3}{a_{\rm out,5}}\frac{e_c}{0.1}a_{p,1}^{-4}, \end{eqnarray} \tag{ 55 }$$

where the numerical estimate is for a p = 1 disk and ρ_p = 3 g cm⁻³.

In the opposite case, when precession is dominated by the disk gravity and A = A_d, one obtains

$$\begin{eqnarray} d_c &=& \frac{3C_D {\rm E}(\sqrt{3}/2)}{8\pi ^2\psi _1}\frac{r}{h} \frac{M_p}{\rho _p a_p^2}e_c\nonumber \\ &\approx & 1 {\rm km} \frac{C_D}{\psi _1}\frac{0.1}{h/r} \frac{e_c}{0.1}M_{p,1}a_{p,1}^{-2}, \end{eqnarray} \tag{ 56 }$$

independent of the disk mass. It is obvious that in the disk-dominated case, d_c is much smaller than in the binary-dominated case for a_p ≲ 1AU, a fact predicted in R13.

This difference can be easily seen in Figure 8, where the situation depicted in panel (a) corresponds to the DD regime where Equation (56) applies. As a result, the planetesimal size for which the low-e_r "waist" in this figure is narrowest is around 1 km. In constrast, Figure 8(b) shows a situation in which the disk gravity has been turned off, so the dynamics are in the BB regime and Equation (55) applies. Not surprisingly, this pushes the characteristic d_p at the narrowest point of the waist to be about 30 km.

Using this reasoning, one might expect the critical "dangerous" size d_p at which e_r ∼ e_c for objects of comparable size to be smaller for more massive disks in which |A_d| ≫ |A_b|. However, this logic directly applies only if e_c were kept the same. In reality, changing A also directly affects the value of e_c; see Equation (29). Figure 5 shows that in practice, the behavior of d_c largely reflects that of e_c, with all the features of e_c maps (e.g., valleys of low d_c) present in d_c maps as well. In particular, the valley of stability shows up prominently in Figures 5(a) and (d).

Zoom In Zoom Out Reset image size

The only noticeable difference with Figure 4 is the increase of d_c with decreasing a_p in the high-M_d (DD) regime; see Figures 5(c) and (d), a behavior that is predicted by Equation (56). Also, in agreement with Equation (55), d_c decreases with increasing a_p in the outer disk for small M_d (upper left in Figures 5(c) and (d)) even though e_c varies weakly there. In this region, planetesimal dynamics are determined predominantly by the binary companion (BB regime of SR15) and Equation (55) applies.

Our next step is to study the behavior of the relative approach velocity v₁₂ between planetesimals with sizes d₁ and d₂. It is this velocity that determines the outcome of their collision.

We now provide a calculation of the distribution df₁₂/dv₁₂ of v₁₂ between the two planetesimal populations, one with the eccentricity vector e_p(d₁) and another with e_p(d₂). In previous sections, we have shown that after the initial transient period, when the free eccentricity damps out, the value of e_p becomes time-independent and is uniquely determined by the planetesimal size. Then the only additional orbital parameter that can give rise to the variation of the relative velocity v₁₂ is the difference in semi-major axes b₁₂ between approaching particles; see Equation (A7) of Appendix A. Using Equations (A4) and (A7), this relation can be written as

$$\begin{eqnarray} &&v_{12}=\Omega a_p\left[e_{12}^2-\frac{3}{4} \left(\frac{b_{12}}{a_p}\right)^2\right]^{1/2}, \end{eqnarray} \tag{ 57 }$$

where a is the mean semi-major axis of both planetesimals, and the condition of close approach x₁₂ = 0 is used. Note that in this expression we ignored the contribution of particle inclination to the velocity. This is a reasonable assumption since we expect eccentricity excitation in the binary plane to dominate over the out-of-plane excitation.

Ida et al. (1993) consider encounters between the two populations of objects with fixed eccentricity vectors e₁ = (k₁, h₁) and e₂ = (k₂, h₂). They derive the following expression for the flux of objects with an eccentricity e₂ approaching a given object with an eccentricity e₁ with random orbital phases, having a separation of their semi-major axes b₁₂ in the range (b₁₂, b₁₂ + db₁₂):

$$\begin{eqnarray} &&dF_{12}=\frac{1}{\pi ^2}\frac{\Sigma _2}{a_p m_2 i_{12}}\frac{v_{12}db_{12}}{\left[e_{12}^2 -(b_{12}/a_p)^2\right]^{1/2}}. \end{eqnarray} \tag{ 58 }$$

Here e₁₂ = [(h₁ − h₂)² + (k₁ − k₂)²]^1/2 is the relative eccentricity between the two particle populations, i₁₂ is their relative inclination, and Σ₂ is the surface density of objects with eccentricity e₂.

Using Equation (57), we can express db in Equation (58) via dv₁₂, resulting in a differential particle flux per unit v₁₂

$$\begin{eqnarray} &&\frac{dF_{12}}{dv_r}=\frac{4}{3\pi ^2}\frac{\Sigma _2 a_p}{m_2 i_{12}} \frac{e_{12}^2 -(3/4)(b_{12}/a_p)^2}{|b_{12}|\left[e_{12}^2 -(b_{12}/a_p)^2\right]^{1/2}}. \end{eqnarray} \tag{ 59 }$$

We now express b₁₂ via v₁₂ using Equation (57) and introduce

$$\begin{eqnarray} &&v_{\rm min}=\frac{1}{2}e_{12} n a_p,\quad v_{\rm max}=e_{12} n a_p. \end{eqnarray} \tag{ 60 }$$

Then it is clear that v_min < v₁₂ < v_max and we can re-write(59) as

$$\begin{eqnarray} &&\frac{dF_{12}}{dv_{12}}=\frac{1}{\pi ^2}\frac{\Sigma _2}{m_2 i_{12}} \frac{v_{12}^2}{\left[\left(v_{\rm max}^2-v_{12}^2\right) \left(v_{12}^2-v_{\rm min}^2\right)\right]^{1/2}}. \end{eqnarray} \tag{ 61 }$$

From this we find that the distribution of relative velocities df₁₂/dv₁₂ of different planetesimals normalized to unity is given by the following expression:

$$\begin{eqnarray} &&\frac{df_{12}}{dv_{12}}=\frac{v_{\rm max}^{-1}}{{\rm E}(\sqrt{3}/2)} \frac{v_{12}^2}{\left[\left(v_{\rm max}^2-v_{12}^2\right) \left(v_{12}^2-v_{\rm min}^2\right)\right]^{1/2}}. \end{eqnarray} \tag{ 62 }$$

Particle sizes enter into this expression only through e₁₂ via Equations (60).

This distribution of relative velocities is shown in Figure 6. It diverges at both v = v_min and v = v_max, but the total particle flux is finite and given by

$$\begin{eqnarray} &&F_{12}=\int _{v_{\rm min}}^{v_{\rm max}}\frac{dF_{12}}{dv_{12}}dv_{12}= \frac{{\rm E}(\sqrt{3}/2)}{\pi ^2} \frac{\Sigma _2 e_{\rm max}}{m_2 i_{12}}. \end{eqnarray} \tag{ 63 }$$

With the distribution function (61), one finds that the mean relative velocity 〈v₁₂〉 ≈ 0.81v_max = 0.81e₁₂n_pa_p, while the rms velocity is given by $v_{\rm rms}=\langle v_{12}^2\rangle ^{1/2}=0.828 e_{12} n_p a_p$.

Zoom In Zoom Out Reset image size

The results of the previous section clearly demonstrate that the relative velocity with which two planetesimals with sizes d₁ and d₂ approach each other prior to collision is determined by their relative eccentricity e₁₂ = |e_p(d₁) − e_p(d₂)|. Using solutions (22) and (23), it is trivial to show that

$$\begin{eqnarray} &&e_{12}=e_c\frac{\left|A\tau _{d,1}-A\tau _{d,2}\right|}{\sqrt{\left(1+A^2\tau _{d,1}^2\right) \left(1+A^2\tau _{d,2}^2\right)}}, \end{eqnarray} \tag{ 64 }$$

where τ_{d, i} ≡ τ_d(d_i), i = 1, 2. According to the results of Section 5, Aτ_d and, subsequently, e₁₂, are functions of (1) the sizes of the colliding planetesimals d_{1, 2} and (2) the binary parameters and local disk properties, which set the values of both e_c and d_c; see Equations (29) and (31). We already explored the latter in Section 7 and now we turn our attention to understanding e_{1, 2}(d₁, d₂).

In Figure 7 we map out e₁₂(d₁, d₂) (as well as the relative velocity v₁₂ = e_{1, 2}v_K) at the location of the planet a_p = 2 AU in the γ Cephei system for different characteristics of the disk, for which a model (1) with p = 1, q = −1 is adopted. We vary disk mass M_d, eccentricity at its outer edge e₀, and its orientation with respect to the binary orbit ϖ_d, one at a time, keeping other disk parameters fixed. All panels clearly show several key invariant features.

Zoom In Zoom Out Reset image size

First, there is a critical size of order d_c, around d₁ = d₂ ∼ (0.1–1) km, at which maps exhibit a "waist," in which e₁₂ is small for collisions of equal size bodies. Second, e₁₂ becomes small for encounters between both the small bodies, with d₁, d₂ ≲ d_c, and for large objects with d₁, d₂ ≳ d_c. Third, e₁₂ saturates at a value roughly independent of d₁ or d₂ for collisions of particles with very different sizes, i.e., when d₁ ≲ d_c ≲ d₂, and vice versa.

These gross features, as well as the variations of the overall velocity scale seen in these maps, are addressed below using the results of Section 5. Given that particles can be in different drag regimes—strong or weak—we will consider several possibilities separately.

9.1. Strong–Strong Encounters

When both planetesimals are in the strong drag regime, d₁, d₂ ≪ d_c, both |Aτ_{d, 1}| ≪ 1 and |Aτ_{d, 2}| ≪ 1. Then Equation (64) predicts that

$$\begin{eqnarray} e_{12}^{\rm ss} & \approx & e_c\left|\left|A\tau _{d,1}\right|- \left|A\tau _{d,2}\right|\right| \end{eqnarray} \tag{ 65 }$$

$$\begin{eqnarray} \hspace{15pt}& \approx & e_c\left|\left(\frac{d_1}{d_c}\right)^{1/2}- \left(\frac{d_2}{d_c}\right)^{1/2}\right|,\hspace{-15pt} \end{eqnarray} \tag{ 66 }$$

where we used Equation (34) to express |Aτ_d| in terms of planetesimal sizes. Since d_{1, 2} ≪ d_c in the strong drag limit, one finds that $e_r^{\rm ss}\lesssim e_c$, which explains the low values of e_r in the lower left corner in the maps in Figure 7.

Physically, in this regime, the relative velocity of two planetesimals is considerably lower than their individual velocities because of the apsidal alignment of their orbits by gas drag, see Marzari & Scholl 2000) and similar magnitudes of e_p.

9.2. Weak–Weak Encounters

When both planetesimals are in the weak drag regime |Aτ_{d, 1}| ≫ 1 and |Aτ_{d, 2}| ≫ 1, one finds using Equation (64) that

$$\begin{eqnarray} e_{12}^{\rm ww} & \approx & e_c\left||A\tau _{d,1}|^{-1}- |A\tau _{d,2}|^{-1}\right| \end{eqnarray} \tag{ 67 }$$

$$\begin{eqnarray} & \approx & e_c\left|\frac{d_c}{d_1}- \frac{d_c}{d_2}\right|, \end{eqnarray} \tag{ 68 }$$

where Equation (38) has been used. Since d_{1, 2} ≫ d_c in the weak drag limit, one again finds that $e_r^{\rm ww}\lesssim e_c$, explaining the low relative eccentricity in the upper right corner in the maps in Figure 7.

In this case, apsidal alignment is again at work, lowering e_r compared to e_p(d₁), e_p(d₂). However, now it is caused by the disk+binary gravity, which affects planetesimals in the same way when they are weakly coupled to gas. This is because the gas damps the free eccentricity, but is not strong enough to significantly change the forced eccentricity.

9.3. Weak–Strong Encounters

When one of the planetesimals (e.g., of size d₁) is in the strong drag regime, |Aτ_{d, 1}| ≪ 1, while the other is in the weak drag regime, |Aτ_{d, 2}| ≫ 1, Equation (64) shows that their relative eccentricity e₁₂ is just

$$\begin{eqnarray} &&e_{12}^{\rm sw}\approx e_c. \end{eqnarray} \tag{ 69 }$$

One can see that e₁₂ is roughly independent of the sizes of particles participating in an encounter.

9.4. Overall e₁₂ Scale as a Function of Disk Parameters

Our work extends and complements existing results on planetesimal dynamics in binaries in several important ways.

First, for the first time, our solutions for e_p in Section 5 simultaneously account for a number of key physical ingredients needed for a complete description of the secular dynamics of planetesimals in binaries: the gravity of both the eccentric disk and the eccentric companion as well as the gas drag, which causes orbital phasing of planetesimals and reduces their relative eccentricity in certain regimes.

Second, we provide a rigorous derivation of the equations of eccentricity evolution due to gas drag (15)–(18) in an eccentric disk. Previously, Adachi et al. (1976) derived analogous equations for the case of a circular disk, while Beaugé et al. (2010) proposed a set of empirical equations similar to Equations (15) and (16) but without proper calculation of the constant pre-factors.

Third, we derive an analytic expression (62) for the relative velocity distribution function df₁₂/dv₁₂ for locally homogeneous populations of objects with fixed eccentricity vectors, which is appropriate in the limit |e_p| ≪ 1 in the presence of gas drag. We also provide an in-depth analysis of e₁₂ behavior for objects of different sizes in systems with different parameters (Section 9). Previously the distribution of planetesimal encounter velocities has been explored only numerically by following a large number of trace particles in simulations of different kinds (Thébault et al. 2006, 2008, 2009; Paardekooper et al. 2008; Fragner et al. 2011). Thus, our derivation of df₁₂/dv_r represents an important analytical step in understanding planetesimal dynamics.

We now provide a more detailed comparison of our results with previous studies and discuss the limitations of this work.

10.1. Comparison of Different Dynamical Approximations

The main novelty of our study is the extension of the line of analytical investigation of disk gravity effects, which was started in R13 and SR15 for axisymmetric and non-axisymmetric disks, respectively, by including gas drag. Previous (semi-)analytical studies of planetesimal dynamics in binaries neglected the gravitational effect of the disk.

Our calculations account for both the precession of planetesimal orbits due to the axisymmetric part of the disk potential and the eccentricity excitation due to its non-axisymmetric component. Disk non-axisymmetry is modeled via its nonzero eccentricity, i.e., m = 1 distortion, which can be a function of radius. We expect this approximation to capture the key effect of the disk asymmetry, as higher-m distortions of the disk shape are relatively small (Marzari et al. 2012).

In Table 1 we summarize some existing (semi-)analytical treatments of planetesimal dynamics including this work (note that this list is not exhaustive), classified according to the physical ingredients that are taken into account. We primarily focus on studies of secular effects to put our work in proper context. Our current results cover all dynamical regimes listed in this table in appropriate limits. The majority of previous studies considered planetesimal dynamics in the presence of gas drag with only the direct binary gravitational perturbations taken into account (Marzari & Scholl 2000; Thébault et al. 2004, 2006, 2008, 2009; Paardekooper et al. 2008). As shown in SR15, this approximation is unwarranted as long as the disk mass M_d ≳ 10⁻² M_☉ since then the disk potential dominates gravitational perturbation.

Table 1. Different Approximations for Planetesimal Dynamics in Binaries

Gravitational Effects	W/O Gas Drag	With Gas Drag
Included
Binary companion only	2, 3	4, 5, 6, 8, 9, 10
Axisymmetric disk	7	1
and binary companion
Non-axisymmetric disk	8	1
and binary companion

Notes. (1) This work; (2) Giuppone et al. 2011; (3) Heppenheimer 1978; (4) Marzari & Scholl 2000; (5) Paardekooper et al. 2008; (6) Beaugé et al. 2010; (7) Rafikov 2013; (8) Silsbee & Rafikov 2015; (9) Thébault et al. 2006; (10) Xie & Zhou 2008.

Download table as:  ASCIITypeset image

We also provide full analytical solutions for test particle dynamics in a general precessing or non-precessing disk without companion perturbation; see Equations (44)–(47). Previously, Beaugé et al. (2010) studied this regime for a precessing disk but did not account for the gravitational effect of the disk (i.e., only gas drag was taken into account). In Figure 8 we illustrate the differences in various descriptions of planetesimal dynamics. It shows the relative eccentricity as a function of planetesimal sizes d₁ and d₂ at 1 AU in an aligned disk of M_d = 10⁻²M_p and e₀ = 0.1 around a primary of γ Cephei in four different limits. Panel (a) presents a full calculation with all physical ingredients (gas drag, gravity of both the eccentric disk and the binary companion) accounted for using the solutions obtained in Section 5.

Zoom In Zoom Out Reset image size

In panel (b) we show how things change if disk gravity is completely switched off by setting A_d = B_d = 0—an approximation common to a number of previous studies (Marzari & Scholl 2000; Thébault et al. 2004, 2006, 2008, 2009; Paardekooper et al. 2008; Beaugé et al. 2010; Leiva et al. 2013). One can see that without disk gravity, relative planetesimal velocities go up by a factor of several. Moreover, the "waist" between the two high-e₁₂ regions in panel (b) is narrowest at d₁ ∼ d₂ ∼ 10² km, which is considerably larger than in panel (a) where this happens for d ∼ 0.3 km objects. This difference is in complete agreement with Equations (55) and (56).

In panel (c) we account for the gravitational effect of a non-axisymmetric disk but neglect gas drag (τ_d → ∞), i.e., use Equations (19)–(21), as was done in SR15. In the absence of gas drag there is no apsidal alignment of planetesimal orbits and they approach each other at random phases. Also, e_r is independent of d₁ and d₂ (the size-dependent drag is absent) explaining uniform color in Figure 8(c). The absence of gas drag results in rather high relative velocities of planetesimals, making their survival in collisions problematic. Thus, apsidal alignment of planetesimal orbits and eccentricity suppression due to gas drag are very important for the proper description of their dynamics.

Finally, in panel (d) we retain only the axisymmetric component of the disk potential, neglecting the eccentricity excitation by the disk, i.e., B_d = 0 but A_d ≠ 0. In this limit, also neglecting gas drag (accounted for here), R13 predicted a dramatic lowering of e_p. A comparison with panel (a) clearly shows this is not the case when gas drag is included, which can be understood by noting that e_c in Equation (29) can significantly deviate from $|{\bf e}_p^{\rm n/drag}|$ because of the e_g contribution. This is why lowering $|{\bf e}_p^{\rm n/drag}|$ by setting B_d = 0 and increasing |A| does not necessarily result in smaller e_c, as expected in R13.

To summarize, simultaneously accounting for all the physical processes affecting planetesimals—gas drag, disk, and secondary gravity—is very important for understanding planetesimal growth. Omission of even a single physical ingredient can significantly affect the conclusions drawn from the dynamical calculations.

Previously, Kley & Nelson (2007) and Fragner et al. (2011) numerically explored planetesimal dynamics in gaseous disks that were evolved using direct hydrodynamical simulations. They accounted for the effect of disk gravity on planetesimal motion and at least some of their calculations assumed coplanarity of the disk and the binary. However, even though the setup of these studies is very similar to that of our present work, some subtle differences prevent direct comparison of their results. In particular, when estimating the relative velocities of planetesimals based on their orbit crossing, Kley & Nelson (2007) do not take into account the apsidal phasing of their orbits (Marzari & Scholl 2000), clearly obvious in their Figure 10. As a result they find very high relative speeds even between equal-size planetesimals, which we believe is an artifact of their neglect of apsidal phasing. Fragner et al. (2011) study the case of a circular binary in which apsidal phasing is naturally absent, resulting in high relative speeds of planetesimals. As a result, the applicability of calculations using circular binaries to understanding planetesimal dynamics in eccentric systems like γ Cep is not obvious.

10.2. Limitations of this Work

We studied secular dynamics of planetesimals and explored prospects for planet formation around one of the components of an eccentric binary. We believe that our study includes most, if not all, of the important physical ingredients relevant for this problem—perturbations due to the binary, gas drag, and gravitational effects of an eccentric disk. This is the first time planetesimal dynamics in binaries have been studied analytically in such generality. The analytical nature of our solutions for planetesimal dynamical variables allowed us to explore their dependence on system parameters in great detail.

Our main results can be summarized as follows.

• 1.  
  We find that under the action of gas drag as well as the gravitational effects of the binary companion and the eccentric disk, the planetesimal eccentricity vector e_p converges to a constant value depending on the planetesimal size and the disk and binary properties. We obtained complete analytical solutions for e_p in the case of a non-precessing disk and analyzed them in detail, extending the results of previous studies.
• 2.  
  We showed that relative particle–gas (Equation (33)) and particle–particle velocities can be expressed as simple functions of only two key parameters—the characteristic eccentricity e_c and planetesimal size d/d_c in units of characteristic size d_c, given by Equations (29) and (31). The behavior of these variables has been explored in detail in Section 7.
• 3.  
  We show that in massive disks containing enough gas to form giant planets (M_d ≳ 10⁻² M_☉), planetesimal dynamics are always in the regime where the apsidal precession of planetesimal orbits is dominated by disk gravity, i.e., in the DB or DD regimes in the classification of SR15. Significantly eccentric (e₀ ≳ 10⁻³) disks also dominate the eccentricity excitation of planetesimals by their gravity (DD regime). This emphasizes the key role of the disk gravity in relation to planet formation in binaries.
• 4.  
  We derive the explicit form of the relative velocity distribution between the populations of planetesimals with different sizes and show that it depends only on the relative eccentricity e₁₂ of the approaching objects.
• 5.  
  In disks aligned with the binary, planetesimals collide with lower velocities than in misaligned disks. Thus, planetesimal growth favors disk–binary apsidal alignment.
• 6.  
  We also present analytical results for the dynamics of planetesimals in precessing disks in certain limits.

Our results will be used in Paper II to understand planet formation in small separation binaries, such as γ Cep and α Cen. They can also be used to understand the circumbinary planet formation.

We are grateful to Jihad Touma for useful discussions.

Here we review local (or guiding center) approximation, which is often used in studies of planetesimal and galactic dynamics (Binney & Tremaine 2008) and forms the basis of the so-called Hill approximation (Hénon & Petit 1986; Hasegawa & Nakazawa 1990). In this approach, eccentric motion of a planetesimal is considered in a locally Cartesian frame (x_p, y_p), with e_x, e_y pointing in the radial and horizontal directions, correspondingly. The origin of this frame is in circular Keplerian motion at some characteristic semi-major axis a₀, which is close to the planetesimal semi-major axis a_p, so that b_p ≡ |a_p − a₀| ≪ a_p. Equations of motion for a particle of mass m_p subject to external force F = (F_x, F_y) can be reduced to

$$\begin{eqnarray} &&\ddot{x}_p-2n_p\dot{y}_p-3 n_p^2 x_p=F_x/m_p,\quad \ddot{y}_p+2n_p\dot{x}_p =F_y/m_p. \nonumber \\ \end{eqnarray} \tag{ A1 }$$

Provided that e_p ≪ 1, one can represent planetesimal motion unperturbed by external forces as

$$\begin{eqnarray} x_p&=&b_p-a_0(k_p\cos n_p t+h_p\sin n_p t),\nonumber \\ y_p &=&\psi _p-\frac{3}{2}n_p b_p t + 2a_0(k_p\sin n_p t-h_p\cos n_p t), \end{eqnarray} \tag{ A2 }$$

where ψ_p is a constant and e_p = (k_p, h_p). This is an exact solution of Equations (A1) with F = 0 and is a superposition of linear shear and epicyclic motion.

Assuming that fluid in a gaseous disk also moves on eccentric Keplerian orbits, motion of the gas can be represented by analogous equations

$$\begin{eqnarray} x_g&=&b_g-a_0(k_g\cos n_p t+h_g\sin n_p t), \nonumber \\ y_g &=&\psi _g-\frac{3}{2}n_p b_g t + 2a_0(k_g\sin n_p t-h_g\cos n_p t). \end{eqnarray} \tag{ A3 }$$

Relative motion of a particular fluid element and a particle is described using relative coordinates x_r = x_p − x_g, y_r = y_p − y_g. According to Equations (A2)

$$\begin{eqnarray} x_r&=&b_r-a_p(k_r\cos n_p t+h_r\sin n_p t), \nonumber \\ y_r &=&\psi _r-\frac{3}{2}n_p b_r t + 2a_p(k_r\sin n_p t-h_r\cos n_p t), \end{eqnarray} \tag{ A4 }$$

where k_r ≡ k_p − k_g, h_r ≡ h_p − h_g are the components of the relative eccentricity vector, b_r ≡ b_p − b_g is the semi-major axis separation between the particle and fluid element, and ψ_r ≡ ψ_p − ψ_g. We have also used the fact that a_g ≈ a₀ ≈ a_p and switched from a₀ to a_p.

Velocity of Keplerian motion in the local approximation is obtained by differentiating Equations (A4) with respect to time. In particular, relative particle–gas velocity is given by

$$\begin{eqnarray} v_{x,r}&=&n_p a_p(k_r\sin n_p t-h_r\cos n_p t), \nonumber \\ v_{y,r}&=&{-}\frac{3}{2}n_p b_r + 2n_p a_p(k_r\cos n_p t+h_r\sin n_p t). \end{eqnarray} \tag{ A5 }$$

Analogous formulae apply to the relative motion of two planetesimals with sizes d₁ and d₂, with the replacement e_r → e₁₂, b_r → b₁₂, (x_r, y_r) → (x₁₂, y₁₂), and so on. In particular, Equations (A4) shows that two objects with |b₁₂| < a_pe₁₂ can experience close approaches. When this happens, x₁₂ = y₁₂ = 0 and b₁₂ can be eliminated from Equations (A5) giving

$$\begin{eqnarray} &&v_{12,y}(x_{12}=0)=\frac{1}{2}n_p a_p(k_{12}\cos n_p t+h_{12}\sin n_p t), \end{eqnarray} \tag{ A6 }$$

(here e₁₂ = (k₁₂, h₁₂)) so that the relative approach velocity (i.e., the velocity unaffected by the mutual gravitational attraction of particles) is

$$\begin{eqnarray} &&v_{12}=n_p a_p\left[k_{12}^2+h_{12}^2-(3/4)(k_{12}\cos n_p t +h_{12}\sin n_p t)^2\right]^{1/2}.\nonumber \\ \end{eqnarray} \tag{ A7 }$$

Whenever a particle is affected by forces other than the stellar gravity, i.e., F ≠ 0, solutions (A2) are no longer strictly valid. However, one can still represent particle motion via these solutions, assuming that orbital elements osculate, i.e., evolve, in time. Hasegawa & Nakazawa (1990) derived equations for the orbital element evolution, in particular

$$\begin{eqnarray} \dot{a}_p&=&\dot{b}_p=\frac{2F_y}{n_p m_p},\quad\! \dot{k}_p=\frac{1}{n_p a_p m_p}(2F_y\cos n_p t+ F_x\sin n_p t),\nonumber \\ \dot{h}_p&=& \frac{1}{n_p a_p m_p}(2F_y\sin n_p t- F_x\cos n_p t). \end{eqnarray} \tag{ A8 }$$

For a given force expression F, these equations, after averaging over the orbital period, represent the extra terms entering the Equations (3) and (4).

Here we derive the full time-dependent solution for planetesimal eccentricity starting with arbitrary initial conditions and assuming that the gas drag is linear, i.e., τ_d in Equation (17) is a constant independent of e_p. We also include a possibility of uniform disk precession so that $\varpi _d(t)=\dot{\varpi }_d t+\varpi _{d0}$. Then Equations (3) and (4) represent a linear system of equations which can be trivially solved to give

$$\begin{eqnarray} &&\left\lbrace \begin{array}{@{}l@{}}k(t)\\ h(t) \end{array} \right\rbrace =e_{\rm free}e^{-t/\tau _d} \left\lbrace \begin{array}{@{}l@{}}\cos \left(At+\varpi _0\right)\\ \sin \left(At+\varpi _0\right) \end{array} \right\rbrace + \left\lbrace \begin{array}{@{}l@{}}k_{\rm f}\\ h_{\rm f} \end{array} \right\rbrace, \end{eqnarray} \tag{ B1 }$$

where the first term represents the free eccentricity, with e_free and ϖ₀ being constant, while the second is the forced eccentricity e_f = (k_f, h_f) = e_{f, b} + e_{f, d}, where e_{f, b} is given by Equation (25) and

$$\begin{eqnarray} {\bf e}_{f,d} &=& \left[\frac{e_g^2+\tau _d^2 B_d^2}{1+\tau _d^2 \left(A-\dot{\varpi }_d\right)^2}\right]^{1/2} \left\lbrace \begin{array}{l}\cos \left(\varpi _d(t)+\phi \right)\\ \sin \left(\varpi _d(t)+\phi \right) \end{array} \right\rbrace,\nonumber \\ \cos \phi &=&\frac{e_g-\tau _d^2 B_d\left(A-\dot{\varpi }_d\right)}{\left(e_g^2+\tau _d^2 B_d^2\right)^{1/2} \left[1+\tau _d^2\left(A-\dot{\varpi }_d\right)^2 \right]^{1/2}}. \end{eqnarray} \tag{ B2 }$$

In the limit of slow precession $|\dot{\varpi }_d|\ll |A|$ one finds that e_f is given by expressions (24)–(27). Generally, the relative planetesimal–gas eccentricity e_r = e_f − e_g is

$$\begin{eqnarray} {\bf e}_r &=& {\bf e}_{f,b} - \tau _d\frac{B_d+e_g\left(A-\dot{\varpi }_d\right)}{\left[1+\tau _d^2\left(A-\dot{\varpi }_d\right)^2\right]^{1/2}} \left\lbrace \begin{array}{@{}l@{}}\cos \left(\varpi _d(t)-\phi _r\right)\\ \sin \left(\varpi _d(t)-\phi _r\right) \end{array} \right\rbrace,\nonumber \\ \cos \phi _r&=&\frac{\tau _d\left(A-\dot{\varpi }_d\right)}{\left[1+\tau _d^2\left(A-\dot{\varpi }_d\right)^2 \right]^{1/2}}. \end{eqnarray} \tag{ B3 }$$

The first forced term e_{f, b} results from excitation by the binary companion. It is constant in time and is independent of $\dot{\varpi }_d$. The second term is induced by the disk via both its gravitational potential and gas drag. This contribution to e_f circulates at the disk precession frequency $\dot{\varpi }_d$ and its amplitude is sensitive to $\dot{\varpi }_d$.

Independent of the initial conditions (i.e., the values of e_free and ϖ₀), the free eccentricity contribution damps out on a characteristic timescale τ_d. As a result, in the long run, e_p inevitably converges to e_f.

In the limit of strong gas drag, τ_d → 0, one finds that e_f → e_g as expected, since drag is strong enough to align planetesimal orbits with fluid trajectories. In this limit, the relative eccentricity between planetesimals of different sizes having different damping times τ_{d, 1} and τ_{d, 2} is

$$\begin{eqnarray} && e_{12}\rightarrow |A-\dot{\varpi }_d| |\tau _{d,1}-\tau _{d,2}| \left[\left(e_c^{\rm pr}\sin \varpi _d\right)^2\right. \nonumber \\ &&\quad +\left.\left(\!e_c^{\rm pr}\cos \varpi _d +k_b\frac{A}{A-\dot{\varpi }_d}\right)^2\right]^{1/2}\!\!, \tau _{d,1},\tau _{d,2}\ll |A-\dot{\varpi }_d|^{-1},\nonumber \\ \end{eqnarray} \tag{ B4 }$$

where $e_c^{\rm pr}$ is given by Equation (45).

In the opposite extreme, τ_d → ∞ (weak drag), one finds ϕ → π and e_f reduces to the forced eccentricity value (with disk precession) obtained in SR15. The relative velocity becomes

$$\begin{eqnarray} && e_{12}\rightarrow \frac{1}{\left|A-\dot{\varpi }_d\right|} \left|\tau _{d,1}^{-1}-\tau _{d,2}^{-1}\right| \left[\left(e_c^{\rm pr}\sin \varpi _d\right)^2 \right.\nonumber \\ &&\quad \left. + \left(\!e_c^{\rm pr}\cos \varpi _d +k_b\frac{A-\dot{\varpi }_d}{A}\right)^2\right]^{1/2}\!, \tau _{d,1},\tau _{d,2}\gg |A-\dot{\varpi }_d|^{-1}.\nonumber \\ \end{eqnarray} \tag{ B5 }$$

Note that in this expression, k_b is multiplied by a factor different from that in Equation (B4). However, it is clear that in both limiting cases, e₁₂ ≪ e₁, e₂, i.e., the relative planetesimal eccentricity is much less than the individual eccentricities e₁ and e₂, a result that remains valid in a precessing disk.

In the case of quadratic drag (10), Figure 2(b) clearly shows the phenomenon of e_p convergence to a quasi-stationary limit cycle behavior, similar to the results of Appendix B. This behavior is further illustrated in Figure 9, where we show the dependence of the limit cycles on planetesimal size d_p and disk precession rate $\dot{\varpi }_d$. Because of the nonlinear drag law, the shapes of the limit cycles in general deviate from ellipses.

Zoom In Zoom Out Reset image size

Nevertheless, their gross features still can be understood in our linear solution (B3). In particular, limit cycles are not centered on (k_r, h_r) = 0 because of the binary companion perturbations, i.e., non-zero e_{f, b} varying as d_p (and τ_d) change. The amplitude of the limit cycles goes down for smaller d_p because τ_d is also smaller, which, according to Equation (B3), reduces the oscillating contribution to e_r. As we vary $\dot{\varpi }_d$ in Figure 9(b), the binary contribution stays unchanged and all limit cycles stay centered on the same point in h_r-k_r space.

Their sizes vary with $\dot{\varpi }_d$ as predicted by Equation (B3). They shrink at high $|\dot{\varpi }_d|\sim |A|$, in agreement with Equation (B3). For slow precession, $|\dot{\varpi }_d|\ll |A|$ limit cycles converge to the trajectory for the non-precessing disk solution (24) in which ϖ_d is set to vary as $\varpi _d(t)=\dot{\varpi }_d t+\varpi _{d0}$. Note that such a convergence to solution (24) is obvious only in the case of $|\dot{\varpi }_d \tau _d|\ll 1$, i.e., when gas drag allows e_p to quickly re-adjust to a new "quasi-static" solution as ϖ_d changes. This is the case shown in Figure 9(b). In the opposite case of $|\dot{\varpi }_d \tau _d|\gg 1$ (and $|\dot{\varpi }_d|\ll |A|$), this convergence is not obvious as the disk precession constantly drives free eccentricity, while the gas drag is not strong enough to quickly damp it. We leave detailed exploration of such topics to a future study.

Now, let us rewrite Equations (3) and (4) in terms of the relative particle–gas eccentricity components k_r and h_r:

$$\begin{eqnarray} &&\frac{dh_r}{dt}=Ak_r-\frac{h_r}{\tau _d}+B_b+ [\left(A-\dot{\varpi }_d\right)e_g+B_d]\cos \varpi _d(t),\nonumber \\ \end{eqnarray} \tag{ C1 }$$

$$\begin{eqnarray} &&\frac{dk_r}{dt}=-Ah_r-\frac{k_r}{\tau _d}- [\left(A-\dot{\varpi }_d\right)e_g+B_d]\sin \varpi _d(t),\nonumber \\ \end{eqnarray} \tag{ C2 }$$

with τ_d given by Equation (18) and dependent upon e_r.

The explicit time dependence of the last terms in these nonlinear equations precludes us from finding their general analytical solutions even in the case of the limit-cycle behavior. However, we can still obtain analytical results for planetesimal eccentricity in the two limiting cases, reviewed next.

First, one can assume that the binary companion dominates eccentricity forcing, which implies that the condition (42) is fulfilled. Then one can drop the last ϖ_d-dependent terms in Equations (C1) and (C2), removing the explicit time dependence from them. This is essentially equivalent to neglecting both the gravitational effect of the disk, i.e., |e_d| → 0, and the gas eccentricity, e_g, compared to |e_b|. As a result, we find a steady-state solution (43) for k_r ≈ k_p and h_r ≈ h_p, which is essentially the Equations (22) and (23) with |e_d|, |e_g| → 0. Then planetesimal dynamics are described by the analytical results of Section 5 with e_c ≈ |B_b/A|.

In the opposite extreme of weak eccentricity excitation by the binary companion we introduce new coordinates H ≡ k_gh_r − h_gk_r, K ≡ h_gh_r + k_gk_r (see Beaugé et al. 2010 for a similar treatment). Then the evolution of H and K is given by

$$\begin{eqnarray} &&\frac{dH}{dt}=\left(A-\dot{\varpi }_d\right)K-\frac{H}{\tau _d}+B_bk_g(t)+ e_g\left[\left(A-\dot{\varpi }_d\right)e_g+B_d\right],\nonumber \\ \end{eqnarray} \tag{ C3 }$$

$$\begin{eqnarray} &&\frac{dK}{dt}=-\left(A-\dot{\varpi }_d\right)H-\frac{K}{\tau _d}+B_bh_g(t). \end{eqnarray} \tag{ C4 }$$

When the eccentricity excitation by the companion is small, we can drop the B_b terms in these equations, removing the explicit time-dependence, which appears because of circulating k_g and h_g. As a result, we find the steady state solutions for H and K in the implicit form

$$\begin{eqnarray} &&K=e_c^{\rm pr}e_g\frac{\left(A-\dot{\varpi }_d\right)^2\tau _d^2}{1+\left(A-\dot{\varpi }_d\right)^2\tau _d^2},\quad H=-e_c^{\rm pr}e_g\frac{\left(A-\dot{\varpi }_d\right)\tau _d}{1+\left(A-\dot{\varpi }_d\right)^2\tau _d^2},\nonumber \\ \end{eqnarray} \tag{ C5 }$$

where τ_d is a function of the relative particle–gas eccentricity $e_r=e_g^{-1}(K^2+H^2)^{1/2}$ given by Equation (44). This solution corresponds to the eccentricity vector e_p fixed in a disk frame, which uniformly precesses at the rate $\dot{\varpi }_d$.

Using these solutions, it is trivial to show that e_p → e_{f, d} given by Equation (B2) with B_b set to zero. That in the weak binary perturbation regime we find the same expression for e_p as in the case of linear drag is not surprising: with B_b = 0 one finds that |e_r| is constant in time, so that τ_d is also constant. Then Equations (C3) and (C4) are the same as in the linear drag case and have the same steady state solutions (C5).