PLANET FORMATION IN BINARIES: DYNAMICS OF PLANETESIMALS PERTURBED BY THE ECCENTRIC PROTOPLANETARY DISK AND THE SECONDARY

Planet-hosting binary systems with separations of several tens of AU present an interesting testbed for planet formation theories. Strong gravitational perturbations induced by the companion excite high eccentricities of planetesimals out of which planets form. Agglomeration of these objects into bigger bodies in mutual collisions, most effective at low relative speeds because of gravitational focusing, may become very ineffective. In a strongly dynamically excited environment, planetesimals would destroy each other instead of growing. This fragmentation barrier presents a very serious problem for planetary growthin binaries.

This issue is particularly severe for binaries with small separation. At the moment, we know (Chauvin et al. 2011; Dumusque et al. 2012) of five planet-hosting systems with eccentric companions (eccentricities ≳ 0.4) and semimajor axes of about 20 AU. Three of them—HD196885, γ Cep, and HD 41004—harbor giant planets with masses above that of Jupiter at 1.6–2.6 AU. At these separations, the eccentricity of a free particle can easily reach 0.1 (Heppenheimer 1978), leading to collisions at speeds of several km s⁻¹ and resulting in destruction of even rather massive (several hundred kilometers in size) objects in collisions as well as smaller planetesimals. Two other systems—α Cen and Gl 86—harbor planets at ≲ 0.1 AU but even these objects have likely formed farther out and then migrated in.

Planetesimal agglomeration must proceed in gaseous protoplanetary disks. It has long been recognized that gas drag is an important agent of planetesimal dynamics (Marzari & Scholl 2000; Thébault et al. 2004, 2006, 2008, 2009; Paardekooper et al. 2008), helping lower relative speeds of planetesimals to some extent. Recently, it has also been realized that the gravitational field of a massive protoplanetary disk can have a strong effect on planetesimal dynamics. In particular, Rafikov (2013b, hereafter R13) has shown that an axisymmetric, massive gaseous disk drives fast precession of planetesimal orbits by its gravity, which effectively suppresses eccentricity excitation by the companion. This mechanism permits growth of even 10 km planetesimals at 2 AU as long as the disk is massive (∼0.1 M_☉) and axisymmetric.

At the same time, hydrodynamic simulations of protoplanetary disks in binaries always find that disks perturbed by the companion develop some degree of non-axisymmetry (Okazaki et al. 2002; Kley et al. 2008; Marzari et al. 2009; Paardekooper et al. 2008), which usually manifests itself as a non-zero disk eccentricity. Such a disk has a non-axisymmetric component of its gravitational field which affects planetesimals in a way similar to the binary companion. Thus, one expects an eccentric gaseous disk to drive planetesimal eccentricity excitation (in addition to that produced by the binary companion), an effect absent in the case of an axisymmetric disk studied in R13. Recent work of Marzari et al. (2013) supports this expectation by showing this effect to operate in circumbinary disks, which can also develop eccentric structure and drive eccentricity growth by their gravity.

The goal of this work is to analyze dynamics of planetesimals in the presence of gravitational perturbations due to both the binary companion and the eccentric disk. To focus on purely gravitational effects, we neglect gas drag in our calculations (it is taken into account in Rafikov & Silsbee 2015a, 2015b). We explore planetesimal dynamics in the secular approximation, neglecting short-period perturbations of planetesimal orbits that average out over the long time intervals. The majority of our results are derived for the case of a non-precessing disk, which is steady with respect to the orientation of the eccentric orbit of the secondary. However, we also explore planetesimal dynamics in the case of a precessing disk.

A significant part of this work is a derivation of the disturbing function due to an eccentric disk, which has been carried out for the first time. Because of the technical nature of this derivation, which we cover in Appendix A, it can be skipped at first reading. The main results are summarized in the main text.

The structure of the paper is as follows. We outline the problem set-up in Section 2 and present basic equations of planetesimal motion and their solutions for the case of non-precessing disk in Section 3. We analyze our solutions and describe four possible dynamical regimes for planetesimal eccentricity excitation in Section 4. Eccentricity behavior as a function of the distance from the primary is discussed in Section 5. In Section 6 we explore the case of a uniformly precessing disk. Our results are discussed in Section 7, where we cover the implications for planetesimal growth (Section 7.1), ways of lowering planetesimal eccentricity (Section 7.2), and comparison with existing numerical results (Section 7.3). Our findings are summarized in Section 8.

We consider a binary star in which the primary and secondary have masses M_p and M_s, and define ν ≡ M_s/M_p. The semimajor axis and eccentricity of the binary are a_b and e_b, and its orientation is specified by apsidal angle ϖ_b.

Coplanar with the binary and orbiting the primary star (this designation is arbitrary) is the eccentric gaseous disk with a non-axisymmetric surface density distribution Σ(r_d, ϕ_d). The disk is eccentric in the sense that trajectories of its fluid elements are confocal ellipses, which in general is not equivalent to Σ being constant along these ellipses (see the discussion of this approximation in Section 7). We define r_d to be the distance from the common focus of the elliptical fluid trajectories, and ϕ_d to be the polar angle with respect to the disk apsidal line; see Figure 1 for an illustration. For every such gaseous trajectory with semimajor axis a_d, we can define the disk surface density at the periastron Σ_p(a_d) and the eccentricity of the fluid trajectory e_d(a_d), which we will simply call disk eccentricity. In general, both Σ_p(a_d) and e_d(a_d) can be arbitrary functions of the fluid semimajor axis a_d, as long as e_d(a_d) varies slowly enough for the particle trajectories to be non-crossing (Ogilvie 2001).

Zoom In Zoom Out Reset image size

Statler (1999) has given the following expression for the surface density behavior in such a disk, assuming that the lines of apsides of all elliptical trajectories are aligned:

$$\begin{eqnarray} &&\Sigma (a_d,\phi _d)=\Sigma _p(a_d)\frac{1-e_d^2-\zeta e_d (1+e_d)}{1-e_d^2-\zeta e_d \left[e_d+\cos E(\phi _d)\right]}, \end{eqnarray} \tag{ 1 }$$

where Σ_p(a_d) is the surface density at the pericenter (ϕ_d = E = 0), as a function of the semimajor axis a_d, E(ϕ_d) is the eccentric anomaly (Murray & Dermott 1999) and ζ ≡ dln e_d(a_d)/dln a_d. Equation (1) has been generalized in Statler (2001) and Ogilvie (2001) to the case of the disk apsidal angle ϖ_d varying with a_d but we will not consider this additional complication here as it adds little new physics to our problem. Interestingly, Equation (1) predicts that surface density is constant along the elliptical fluid trajectory if e_d is not varying with a_d, i.e., ζ = 0.

Throughout this work we assume simple power-law scalings

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

for a_in < a_d < a_out, where a_in and a_out are the semimajor axes of the innermost and outermost fluid trajectories, and Σ₀ and e₀ are the pericenter surface density and eccentricity at the outer edge of the disk. If the semimajor axis of the innermost fluid trajectory a_in ≪ a_out, as expected for realistic disks, then Σ₀ can be directly related to the 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 as

$$\begin{eqnarray} &&\Sigma _0=\frac{2-p}{2\pi }\frac{M_d}{a_{\rm out}^2}, \end{eqnarray} \tag{ 3 }$$

where we neglected disk ellipticity (see below) and assumed p < 2, so that most of the disk mass is concentrated in its outer part.

We will neglect the precession of the binary apsidal line caused by the gravity of the circumprimary disk, as the corresponding precession period is considerably longer than other timescales of the problem. We will also focus predominantly on the case of a non-precessing disk. We cover the precessing disk case in Appendix C and Section 6.

Our focus is on the dynamics of planetesimals embedded in the gaseous disk. We characterize planetesimal orbits by semimajor axis a_p, eccentricity e_p, and apsidal angle ϖ_p.

Even though expression (1) does not assume e_d to be small, in the rest of the paper we will take both the disk and planetesimal eccentricities to be small, e_d(r) ≪ 1 and e_b ≪ 1.

We study planetesimal dynamics, taking into account gravitational perturbations from both the binary companion and the eccentric disk. We perform calculations in the secular approximation (Murray & Dermott 1999) by averaging the planetesimal disturbing function R over time, thus eliminating the short-period terms, and keeping only the slowly varying contributions up to second order in the planetesimal eccentricity e_p and to lowest order in disk eccentricity e_d (in all terms).

3.1. Disturbing Function Due to the Disk

In Appendix A we provide a detailed calculation of the planetesimal disturbing function R_d due to a non-axisymmetric disk with surface density and eccentricity distributions given by Equations (2). This calculation is very general and can be applied to an arbitrary eccentric disk, not necessarily around one of the components of the binary. In particular, it can be used to study planetesimal motion in a circumbinary disk. This calculation thus represents an important stand-alone result of this work.

We show in Appendix A that in the secular approximation and to lowest order in e_d and e_p, the disturbing function due to the eccentric disk with orientation ϖ_d (independent of the distance from the primary) has the form

$$\begin{equation} R_d = a_p^2n_p \left[\frac{1}{2}A_de_p^2 + B_de_p \cos {(\varpi _p - \varpi _d)}\right], \end{equation} \tag{ 4 }$$

where

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

where $n_p\equiv \sqrt{{{GM_p/a_p^3}}}$ is the planetesimal mean motion, and dimensionless constants ψ₁ and ψ₂ are given by Equations (A33) and (A34). In deriving this expression for R_d we used Equation (A31), in which we dropped the term independent of e_p.

Coefficients ψ₁ and ψ₂ are functions of the power-law indices p, q, characterizing the disk structure as well as the distance a_p with respect to the disk boundaries. Figure 2 shows the behavior of ψ₁ and ψ₂ for several values of p, q, and different α₁ ≡ a_in/a_p ⩽ 1, α₂ ≡ a_p/a_out ⩽ 1 computed according to Equations (A33) and (A34). One can see that for the selected values of p and q, both ψ₁ and ψ₂ converge to values depending only on p and q in the limit of α₁ → 0, α₂ → 0. Indeed, in Appendix A, we show that as long as

$$\begin{eqnarray} &&-1<p < 4\quad \;{\rm and}\quad -2 < p + q < 5 \end{eqnarray} \tag{ 7 }$$

the values of ψ₁ and ψ₂ are determined locally, by the surface density and e_d behavior in the vicinity of a_p (see Figure 3 for illustration). In this case, for a disk spanning more than about an order of magnitude in radius and a_in ≲ a_p ≲ a_out the gravitational effect of disk parts near the boundaries is not important. Then ψ₁ and ψ₂ only weakly depend on α_{1, 2} and can be well approximated by Equations (A37) and (A38). Their values in this limit are shown in Figure 4 as functions of p and p+q. This is how these coefficients often will be treated (i.e., as constants) in the following analysis, though as can be seen in Figure 2, this approximation breaks down near the boundaries of the disk.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We verified our analytical derivation of R_d given by Equations (4)–(6) in several different ways. In particular, in the case of an axisymmetric disk B_d = 0, we made sure that in this case R_d coincides with the expressions derived in R13 for a surface density profile with p = 1 and in Rafikov (2013b) for an arbitrary p, based on the results of Ward (1981). The accuracy of our results in the case of a non-axisymmetric disk is verified by direct integration of particle motion discussed in Section 3.5.

3.2. Disturbing Function Due to the Binary

Another perturbation to the planetesimal motion is provided by the companion star. For an external binary companion, this is given by (Murray & Dermott 1999)

$$\begin{equation} R_b = a_p^2n_p \left[\frac{1}{2}A_be_p^2 + B_b e_p \cos {(\varpi _p - \varpi _b)}\right], \end{equation} \tag{ 8 }$$

where

$$\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,\\ 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{ 9 }$$

Here α_b ≡ a_p/a_b and

$$\begin{eqnarray} &&b_{s}^{(j)}(\alpha) = \frac{1}{\pi } \int _0^{2\pi } \frac{\cos {(j\theta)} d\theta }{(1 - 2 \alpha \cos {\theta } + \alpha ^2)^s} \end{eqnarray} \tag{ 11 }$$

stands for the standard Laplace coefficient. The approximate expressions assume α_b ≪ 1, which is a reasonable assumption. Equations (9) and (10) are valid up to the leading order in e_b ≪ 1; more accurate expressions can be found in Heppenheimer (1978) or R13.

3.3. Full Planetesimal Disturbing Function

Given that the binary precession due to disk gravity is slow, the orientation of the orbital ellipse of the secondary can be approximated as fixed in time. Then, without loss of generality, we may choose the binary apsidal line as the reference direction, in which case ϖ_b = 0. The total (disk plus star) disturbing function R = R_d + R_b is then given by

$$\begin{eqnarray} &&R = a_p^2n_p \left [\frac{1}{2}Ae_p^2 + B_d e_p \cos (\varpi _p-\varpi _d) + B_b e_p \cos {\varpi _p}\right],\nonumber\\ \end{eqnarray} \tag{ 12 }$$

where

$$\begin{eqnarray} &&A = A_d + A_b. \end{eqnarray} \tag{ 13 }$$

We now introduce the planetesimal eccentricity vector e_p = (k_p, h_p), where

$$\begin{eqnarray} &&k_p = e_p \cos \varpi _p,\quad h_p = e_p \sin \varpi _p. \end{eqnarray} \tag{ 14 }$$

Then R can be written in terms of h_p and k_p as follows:

$$\begin{eqnarray} R &=& a_p^2n_p\left[\frac{1}{2}A\left(h_p^2 + k_p^2\right) + \left(B_b + B_d\cos \varpi _d\right)k_p\right. \nonumber \\ &&\left. + B_d \sin \varpi _d h_p\right]. \end{eqnarray} \tag{ 15 }$$

3.4. Evolution Equations and their Solution

In secular planar approximation, only the eccentricity e_p and apsidal angle ϖ_p of the planetesimal orbit vary in time. We study this process by following the evolution of k_p and h_p using Lagrange equations (Murray & Dermott 1999)

$$\begin{equation} \frac{dk_p}{dt} = -\frac{1}{n_pa_p^2} \frac{\partial R}{\partial h_p},\quad \frac{dh_p}{dt} = \frac{1}{n_pa_p^2} \frac{\partial R}{\partial k_p}. \end{equation} \tag{ 16 }$$

With R given by the expression (15), the evolution equations become

$$\begin{eqnarray} \frac{dk_p}{dt} & = & -Ah_p - B_d \sin \varpi _d.\\ \frac{dh_p}{dt} & = & Ak_p + B_b + B_d \cos \varpi _d, \end{eqnarray} \tag{ 17 }$$

This is the key system of equations for our work, valid as long as the orientation of elliptical fluid trajectories, determined by ϖ_d, is independent of radius.

Note that in deriving this system we did not make any assumptions regarding the time behavior of ϖ_d. Thus, ϖ_d in Equations (17) and (18) can be an arbitrary function of time, which makes this system of equations applicable to rigidly precessing disks as well as disks in which the common apsidal line librates around some equilibrium orientation.

However, for simplicity, we start with a case when ϖ_d = const, i.e., when the disk shape is fixed in the frame of the binary. The precessing disk case is covered in Section 6. We solve Equations (17) and (18) assuming an initially circular planetesimal orbit, i.e., k_p(0) = h_p(0) = 0. The solution

$$\begin{eqnarray} &&\left\lbrace\!\! \begin{array}{l}k_p(t)\\ h_p(t) \end{array} \!\!\right\rbrace & = & {\bf e}_p(t) = {\bf e}_{{\rm forced},b}+ {\bf e}_{{\rm forced},d}+{\bf e}_{\rm free},\\ {\bf e}_{{\rm forced},b} & = & -\frac{B_b}{A} \left\lbrace\!\! \begin{array}{l}1 \\ 0 \end{array} \!\!\right\rbrace,\\ {\bf e}_{{\rm forced},d} & = & -\frac{B_d}{A} \left\lbrace\!\! \begin{array}{l}\cos \varpi _d \\ \sin \varpi _d \end{array} \!\!\right\rbrace,\\ {\bf e}_{\rm free}(t) & = & A^{-1}\sqrt{B_d^2 + 2B_dB_b\cos {\varpi _d} + B_b^2} \nonumber \\ &&\times \left\lbrace\!\! \begin{array}{l}\cos (At + \phi) \\ \sin (At + \phi) \end{array} \!\!\right\rbrace, \end{eqnarray} \tag{ 19 }$$

is decomposed into three distinct contributions: e_{forced, b} is the forced eccentricity due to binary potential, e_{forced, d} is the forced eccentricity due to disk potential, and e_free(t) is the free eccentricity vector rotating at the precession rate A, with the phase ϕ given by equation

$$\begin{eqnarray} &&\sin \phi =\frac{B_d\sin \varpi _d}{\sqrt{B_d^2 + 2B_dB_b\cos {\varpi _d} + B_b^2}}. \end{eqnarray} \tag{ 23 }$$

Variation of eccentricity $e_p=(h_p^2+k_p^2)^{1/2}$ is given by a simple formula

$$\begin{eqnarray} &&e_p(t) = \frac{2}{A}\left|\sin \frac{At}{2}\right| \sqrt{B_d^2 + 2B_dB_b\cos {\varpi _d} + B_b^2}. \end{eqnarray} \tag{ 24 }$$

This result shows that the maximum eccentricity ranges between 2|(|B_d| − |B_b|)/A| and 2|(|B_d| + |B_b|)/A| depending on the value of ϖ_d.

For the subsequent discussion, we will be using a characteristic eccentricity of

$$\begin{equation} e_{\rm {char}} = 2 \frac{|B_b| + |B_d|}{|A|}, \end{equation} \tag{ 25 }$$

which is an upper bound on the e_p. This estimate ignores the dependence of e_p on ϖ_d and overlooks some interesting cases when e_p can be significantly lower than e_char, e.g., when

$$\begin{eqnarray} &&|B_b| \approx |B_d| \quad\;{\rm and}\quad \cos \varpi _d \approx -\;{\rm sgn}\;\left(B_dB_b\right), \end{eqnarray} \tag{ 26 }$$

(sgn(z) is a sign function) a possibility that is discussed in more detail in Section 7.2.5.

3.5. Comparison with Direct Orbit Integrations

To test our analytical prescription (4)–(6) for the disk disturbing function R_d, we compared our theory with the results of direct numerical integration of planetesimal motion in the gravitational field of an eccentric disk. We consider a disk extending from a_in = 0.1 AU to a_out = 5 AU and having Σ_p(1 AU) = 100 g cm⁻². To isolate effects of the disk gravity, we set the mass of the secondary to zero. The details of our numerical calculations are described in Appendix B.

Numerical results were then compared with analytical solutions obtained in the previous section, and the outcomes are shown in Figure 5 in the form of planetesimal eccentricity e_p and apsidal angle ϖ_p dependence on time. We tried different initial conditions for the planetesimal orbit, but in this figure, we concentrate on the case of zero initial planetesimal eccentricity, when the analytical solution is given by Equation (24) with B_b = 0.

Zoom In Zoom Out Reset image size

One can see that irrespective of the parameters of our integrations, the agreement between theory and numerical results is very good. The amplitude of e_p variation is always in excellent agreement with theory (the difference being less than a percent), even for a disk eccentricity at the outer edge as high as e₀ = 0.2; see panel (b). The period of secular oscillations is within several percent of our analytical prediction 2π/A_d given by Equation (5) in the high-eccentricity case e₀ = 0.2 for a disk with p = −q = 1. However, this discrepancy is considerably smaller in other cases shown.

Such deviations (infact much larger in amplitude) between orbit integrations and linear secular theory, predominantly in the period of the variation, have been previously documented in the case of perturbation by the eccentric binary companion alone (Thébault et al. 2006; Barnes & Greenberg 2006). Giuppone et al. (2011) find discrepancies in both the amplitude and period of e_p oscillations at the level of ∼50% when secular theory predicts e_p ≳ 0.1, but, as Figures 5(b) and (d) clearly demonstrate, the agreement between theory and simulations in the case of a disk is much better even when e_p is as high as 0.2–0.4. Most likely this is because the smooth mass distribution of the disk reduces the amplitude of its higher-order gravitational multipoles and allows secular theory, which goes only to octupole order, to better capture the main effects of the disk gravity.

The general conclusion one can draw from the comparisons shown in Figure 5 is that the secular theory for perturbations due to the disk developed in Appendix A works very well and our analytical results (A33) and (A34) for the behavior of coefficients ψ₁ and ψ₂ are correct.

We will now consider different regimes of planetesimal dynamics. We start by using Equation (3) to express the disk-related precession rate A_d and the eccentricity excitation coefficient B_d via the disk mass M_d:

$$\begin{eqnarray} A_d & = & (2-p)\psi _1 n_p \frac{M_d}{M_p}\left(\frac{a_p}{a_{\rm {out}}}\right)^{2-p},\\ B_d & = & \frac{2-p}{2}\psi _2 n_p \frac{M_d}{M_p}\left(\frac{a_p}{a_{\rm {out}}}\right)^{2-p} e_d(a_p) ;\end{eqnarray} \tag{ 27 }$$

see Equations (5) and (6). e_d(a_p) is given by Equation (2).

These expressions show that disk-driven planetesimal eccentricity is determined, in part, by the values of power-law indices p and q. Unfortunately, these parameters are rather poorly known for real protoplanetary disks. Based on standard accretion disk theory, R13 advocated the use of p ≈ 1 for the circumstellar disks in binaries. However, this choice is subject to uncertainty in our knowledge of the radial behavior of the viscous α-parameter, thermal structure of the disk, etc. Thus, in this work, we explore a range of values of p.

Equally uncertain is the choice of the disk eccentricity slope q. If one were to neglect the self-gravity, pressure, and viscous forces in the gaseous disk, then its fluid elements would behave as free particles perturbed by the binary companion and have their eccentricity scaling linearly with a_p (Heppenheimer 1978; also Equation (34)), e_d∝a_p, so that q = −1. This behavior is at least approximately supported by the numerical results of Okazaki et al. (2002) and semi-analytical calculations of Paardekooper et al. (2008) within a range of radii. Other authors find e_d to exhibit more complicated, non-power-law behavior (Kley et al. 2008; Marzari et al. 2009). Despite that, in this work, we will predominantly stick to using q = −1, but sometimes we will consider other values of q < 0.

All disk models considered in this paper have an eccentricity e_d increasing with radius and a surface density decreasing with radius. Under these natural assumptions, the disk should dominate the motion of planetesimals close to the primary star since A_b and B_b grow very rapidly with a_p (while A_d and B_d can even decay with a_p for certain values of p and q). Similarly, for less massive disks, in the outer parts of the disk, the binary dominates both the precession and eccentricity excitation of planetesimals. Then we may ignore the disk-driven perturbations unless the binary orbit is completely circular, in which case eccentricity excitation is solely due to the gravity of the elliptical disk.

Using Equations (9) and (27), we can quantify this logic by forming a ratio

$$\begin{eqnarray} \left|\frac{A_d}{A_b}\right| & = & \frac{4|(2-p)\psi _1|}{3} \frac{M_d}{\nu M_p} \left(\frac{a_b}{a_{\rm {out}}}\right)^{2-p} \left(\frac{a_p}{a_b}\right)^{-(1+p)}. \quad \end{eqnarray} \tag{ 29 }$$

Disk (binary) terms dominate the planetesimal precession rate when |A_d/A_b| ≳ 1 (|A_d/A_b| ≲ 1); see Figure 6.

Zoom In Zoom Out Reset image size

We do an analogous calculation for eccentricity excitation using Equations (10) and (28):

$$\begin{eqnarray} \left |\frac{B_d}{B_b}\right| & = & \frac{8|(2-p)\psi _2|}{15} \frac{e_0}{e_b}\frac{M_d}{\nu M_p}\left(\frac{a_b}{a_{\rm {out}}}\right)^{2-p-q} \left(\frac{a_p}{a_b}\right)^{-(2+p+q)}, \nonumber\\ \end{eqnarray} \tag{ 30 }$$

where e₀ is the disk eccentricity at its outer edge. Again, disk (binary) dominates planetesimal eccentricity excitation when |B_d/B_b| ≳ 1 (|B_d/B_b| ≲ 1); see Figure 6.

Conditions (29) and (30) define special locations in the disk, where the ratios |A_d/A_b|, |B_d/B_b| become equal to unity. We find that |A_d/A_b| = 1 at

$$\begin{eqnarray} a_A & = & a_b\left[\frac{4|\psi _1 (2 - p)|}{3} \frac{M_d}{\nu M_p} \left(\frac{a_b}{a_{\rm out}}\right)^{2-p}\right]^{1/(1+p)} \nonumber \\ & \approx & 0.16a_b\left[\frac{M_d/(\nu M_p)}{0.01} \frac{0.25}{a_{\rm out}/a_b}\right]^{0.5}, \end{eqnarray} \tag{ 31 }$$

while |B_d/B_b| = 1 at

$$\begin{eqnarray} a_B & = & a_b\left[\frac{8|\psi _2(2-p)|}{15}\frac{M_d}{\nu M_p} \frac{e_{0}}{e_b}\left(\frac{a_b}{a_{\rm {out}}}\right)^{2-p-q}\right]^{1/(2+p+q)} \nonumber \\ & \approx & 0.11a_b \frac{0.25}{a_{\rm out}/a_b} \left[\frac{M_d/(\nu M_p)}{0.01} \frac{e_0/e_b}{0.1}\right]^{0.5}, \end{eqnarray} \tag{ 32 }$$

where numerical estimates are for a disk model with p = 1, q = −1 (|ψ₁(1)| = 0.5, |ψ₂(0)| = 1.5).

For the parameters adopted in these estimates, both a_A and a_B lie within the disk, at separations of 2–3 AU for a_b = 20 AU (with a_out = 5 AU), which is outside the semimajor axes of the planets in binaries detected so far. The obvious implication is that these planets have formed in the part of the disk where secular effects were dominated by the disk gravity rather than by the secondary. This suggests that disk gravity plays a decisive role is determining planetesimal dynamics in the planet-building zone.

Using ratios (29) and (30), we now describe different possible regimes of the planetesimal eccentricity behavior, as illustrated in Figure 6. We identify each regime using a two-letter notation in which the first letter describes what dominates the planetesimal precession rate A, while the second refers to the dominance of eccentricity excitation (e.g., "Case DB" means that |A_d/A_b| ≳ 1 and |B_d/B_b| ≲ 1). In Figure 7 we map out these different dynamical regimes in the space of the scaled disk mass M_d/(νM_p) and planetesimal semimajor a_p/a_b for different disk models (combinations of p, q, e₀/e_b).

Zoom In Zoom Out Reset image size

4.1. Case DD: Disk Dominates Both Precession and Excitation

At small separations from the primary, a_p ≲ a_A, a_B, the disk dominates both the precession and the eccentricity excitation of planetesimals, so A ≈ A_d and |B_b| ≪ |B_d|. In this case, the characteristic planetesimal eccentricity (25) tends to

$$\begin{eqnarray} &&e_p^{\rm DD}(a_p)\rightarrow 2\left|\frac{B_d}{A_d}\right| = \left|\frac{\psi _2 }{\psi _1}\right|e_d(a_p). \end{eqnarray} \tag{ 33 }$$

In this regime, the maximum planetesimal eccentricity is of the order of the local disk eccentricity, since |ψ_{1, 2}| ∼ 1; for example, ignoring edge effects e_p(a_p) → 3e_d(a_p) for a p = 1, q = −1 disk. Thus, an elliptical disk is capable of exciting planetesimal eccentricity of the order of its own eccentricity e_d purely by its non-axisymmetric gravitational field. In this regime, planetesimal eccentricity should increase with a_p because e_d(a_p) is expected to be a growing function of a_p.

Figure 7 demonstrates that this dynamical regime is unavoidable for a_p ≲ 1 AU even for relatively small disk masses, down to M_d ∼ 10⁻³ M_☉.

4.2. Case BB: The Binary Dominates Both Precession and Excitation

In the opposite limit, far from the primary, as a_A, a_B ≲ a_p (which is, of course, possible only if a_A, a_B ≲ a_out), planetesimal dynamics are governed completely by the binary potential. The contribution from the disk is insignificant so that both A ≈ A_b and |B_d| ≪ |B_b|. This is the limit of planetesimal dynamics in a diskless binary, which has been investigated by Heppenheimer (1978).

In this case, planetesimal eccentricity is given by

$$\begin{eqnarray} &&e_p^{\rm BB}(a_p)\rightarrow 2\left|\frac{B_b}{A_b}\right| = \frac{5}{2} \frac{a_p}{a_b} e_b, \end{eqnarray} \tag{ 34 }$$

in agreement with Heppenheimer (1978).

Figure 7 shows that Case BB is important for a broad range of separations, down to 1 AU, when the disk mass is very small, ≲ 10⁻³ M_☉. However, for more massive disks with M_d ≳ 10⁻² M_☉, this regime never emerges for a_p < a_out. Thus, in compact binaries (a_b ∼ 20 AU) with massive disks the classical result of Heppenheimer (1978) may never actually apply.

4.3. Case BD: The Binary Dominates Precession, and the Disk Dominates Excitation

In between the two limiting cases covered in Sections 4.1 and 4.2, there are other dynamical regimes.

Provided that a_A < a_B, there exists a region in the disk with a_A ≲ a_p ≲ a_B, where planetesimal precession is dominated by the binary companion (A ≈ A_b), while eccentricity excitation is determined by the disk gravity (|B_d| ≫ |B_b|). In this limit, planetesimal eccentricity is given by

$$\begin{eqnarray} e_p^{\rm BD}(a_p) & \rightarrow & 2\left|\frac{B_d}{A_b}\right| = \frac{4|\psi _2(2-p)|}{3} e_d(a_p) \frac{M_d}{\nu M_p} \nonumber \\ && \times \left(\frac{a_b}{a_{\rm {out}}}\right)^{2-p}\left(\frac{a_p}{a_b}\right)^{-(1+p)}. \end{eqnarray} \tag{ 35 }$$

Using this expression and Equations (29) and (30), one can easily show that

$$\begin{eqnarray} e_p^{\rm BD}(a_p) & = & e_p^{\rm DD}(a_p)\left(\frac{a_p}{a_A}\right)^{-(1+p)} \nonumber \\ & = & e_p^{\rm BB}(a_p)\left(\frac{a_p}{a_B}\right)^{-(2+p+q)}. \end{eqnarray} \tag{ 36 }$$

Since a_A ≲ a_p ≲ a_B in Case BD, this result implies (for p > −1, p + q > −2) that $e_p^{\rm BB}(a_p)\lesssim e_p^{\rm BD}(a_p)\lesssim e_p^{\rm DD}(a_p)$. It is then clear that Case BD requires $e_p^{\rm DD}(a_p)\gtrsim e_p^{\rm BB}(a_p)$ locally, i.e., according to Equation (33), that the disk eccentricity e_d(a_p) be higher than planetesimal eccentricity $e_p^{\rm BB}$ in a diskless case for the values of p and q explored in this paper. This situation may not be easy to realize in practice since pressure and viscous forces may tend to reduce (and not increase) the eccentricity of fluid elements compared to that expected for test particles (i.e., $e_p^{\rm BB}$).

Figure 7 shows that indeed this dynamical regime requires rather special conditions to be realized, such as the relatively high value of the disk eccentricity e₀/e_b. Even then it typically occupies a narrow range of separations; see Figures 7(a) and (b). This is because disk models in these two panels have $e_d(a_p)\approx e_p^{\rm BB}(a_p)$, essentially eliminating the Case BD region. In Figure 7(c), we do display a model with $e_d(a_p)\gtrsim e_p^{\rm BB}(a_p)$ close to the primary (we take $e_d\propto a_p^{1/2}$, while $e_p^{\rm BB}\propto a_p$) so that Case BD emerges at small M_d and relatively small a_p. However, as we mentioned before, this may not be a typical situation.

4.4. Case DB: The Disk Dominates Precession, and the Binary Dominates Excitation

Now we look at the opposite case of a_B < a_A, which emerges when e₀/e_b is low. Within the range a_B ≲ a_p ≲ a_A, planetesimal precession is dominated by the disk gravity (A ≈ A_d), while eccentricity excitation is determined predominantly by the secondary star (|B_d| ≪ |B_b|). This is the approximation of a massive axisymmetric disk discussed in R13. In agreement with that work, we find the maximum eccentricity to follow

$$\begin{eqnarray} e_p^{\rm DB}(a_p)\rightarrow 2\left|\frac{B_b}{A_d}\right| & = & \frac{15}{8|\psi _1(2-p)|}e_b\frac{\nu M_p}{M_d} \nonumber \\ && \times \left(\frac{a_{\rm out}}{a_b}\right)^{2-p} \left(\frac{a_p}{a_b}\right)^{2+p}. \end{eqnarray} \tag{ 37 }$$

This expression and Equations (29) and (30) imply that

$$\begin{eqnarray} e_p^{\rm DB}(a_p) & = & e_p^{\rm DD}(a_p)\left(\frac{a_p}{a_B}\right)^{2+p+q} \nonumber \\ & = & e_p^{\rm BB}(a_p)\left(\frac{a_p}{a_A}\right)^{1+p}. \end{eqnarray} \tag{ 38 }$$

Because now a_B ≲ a_p ≲ a_A, we see that $e_p^{\rm DD}(a_p)\lesssim e_p^{\rm DB}(a_p)\lesssim e_p^{\rm BB}(a_p)$. Then it follows that the DB regime requires $e_d(a_p)\sim e_p^{\rm DD}(a_p)\lesssim e_p^{\rm BB}(a_p)$ in non-pathological cases.

According to Figure 7, this dynamical regime is rather common at low e₀/e_b, but is difficult to realize inside the disk for higher e₀/e_b. For some models (e.g., see Figures 7(d) and (e)), dynamics are in the DB regime over an extended region of the disk.

To illustrate the results of the previous section, in Figure 8 we show profiles of planetesimal eccentricity computed for different disk models. For reference, each of the panels displays planetesimal eccentricity for the diskless case ($e_p^{\rm BB}(a_p)$, dark blue, big-dotted) as well as e_p for the case with no secondary ($e_p^{\rm DD}(a_p)$, red dot-dashed). In the left panels, we have chosen a disk eccentricity e_d(a_p) very close to the eccentricity of a free particle in the binary potential, which explains why the curves of $e_p^{\rm BB}(a_p)$ and $e_p^{\rm DD}(a_p)$ almost overlap. In the right panels, e_d is reduced by an order of magnitude and the two curves are well separated.

Zoom In Zoom Out Reset image size

Note that the $e_p^{\rm DD}(a_p)$ curve does not follow the simple power law in a_p as one would have expected based on Equation (33) and the assumption of ψ₁, ψ₂ being constant—it clearly deviates from this simple form at the disk edges. This is because near the disk edge, boundary terms neglected in computing the Figure 4 start to affect the values of ψ₁ and ψ₂ in a non-trivial manner; see Figure 2.

We plot eccentricity profiles for different values of the disk mass. At the lowest disk mass, M_d = 10⁻³M_p, planetesimal eccentricity e_p starts out very high in the outer disk (in the BB regime; see Figure 7). A notable feature of this profile is the secular resonance located at ≈1.5 AU and causing e_p to diverge and stay above $e_p^{\rm BB}(a_p)$. Its existence was predicted in R13 and Rafikov (2013a) for the case of circumprimary and circumbinary disks respectively. Later, Meschiari (2014) confirmed the emergence of this resonance in massive circumbinary disks using numerical simulations of planetesimal dynamics.

The origin of this resonance lies in the fact that A_b is always positive, whereas for the disks that we are considering, A_d is negative; see Figure 4. This means that at a_A (see Equation (31)), where |A_d| = |A_b| one actually has A = 0 and our secular solution (3) diverges. Inward of the resonance, e_p rapidly goes down (in DB and DD regimes) and asymptotically approaches $e_p^{\rm DD}(a_p)$ for a_p ≲ 0.5 AU.

For a somewhat more massive disk M_d = 10⁻²M_p the e_p profile looks very different—it does not exhibit secular resonance (since a_A is now outside the outer disk edge a_out) and generally features lower values of e_p. This happens because with such a massive disk, planetesimal excitation is never in the BB regime. Disk gravity governs particle dynamics essentially through the whole disk.

This is even more so for the M_d = 10⁻¹M_p disk. At this high mass, planetesimal eccentricity curves closely follow $e_p^{\rm DD}(a_p)$ for all a_p. As a result, in low-e₀ disks, e_p can be appreciably lower than what it is if planetesimals are affected by the gravity of the binary companion alone, similar to the case studied in R13. Somewhat counerintuitively, adding an additional perturber—a massive disk—to the system does not heat it up dynamically but in fact reduces planetesimal random velocities.

In all cases, we see that e_p is above the smaller of the $e_p^{\rm BB}(a_p)$ and $e_p^{\rm DD}(a_p)$. Thus, nonzero disk eccentricity introduces a lower limit on the e_p value.

So far, we have been dealing with the case of a non-precessing disk that keeps its orientation fixed in the frame of the binary orbit. However, simulations often find that gas disks in binaries not only develop a non-zero eccentricity but also precess (Okazaki et al. 2002; Paardekooper et al. 2008; Marzari et al. 2009). Thus, it is important to discuss how planetesimal dynamics changes in the case of a precessing disk.

In Appendix C we present the extension of our solutions for the planetesimal eccentricity in Section 3.4 to the case of a disk that precesses as a solid body at a constant rate¹ $\dot{\varpi }_d$. We find that the eccentricity vector can again be separated into three distinct contributions; see Equation (C3): (1) a standard forced eccentricity vector due to the binary companion with an amplitude of |e_{forced, b}| = |B_b/A|, which is stationary in the binary frame, (2) a forced eccentricity vector due to the disk with an amplitude of $|{\bf e}_{{\rm forced},d}|=|B_d/(A-\dot{\varpi }_d)|$, rotating at a rate of $\dot{\varpi }_d$, and (3) the free eccentricity term with an amplitude of

$$\begin{eqnarray} |{\bf e}_{\rm free}| & = & \frac{1}{|A(A-\dot{\varpi }_d)|} [(AB_d)^2 + \left(B_b(A-\dot{\varpi }_d)\right)^2 \nonumber \\ && + 2AB_dB_b(A - \dot{\varpi }_d)\cos \varpi _{d0}]^{1/2} \end{eqnarray} \tag{ 39 }$$

rotating at the precession rate A (here ϖ_d0 is the value of ϖ_d at t = 0).

The expression for the characteristic eccentricity becomes more complicated and depends on the value of ϖ_d0. The maximum possible eccentricity (for planetesimals starting with h_p(0) = k_p(0) = 0) is reached when ϖ_d(0) = ϖ_d0 = 0 or π (disk and binary periapses aligned or anti-aligned initially), depending on the signs of A, B_d, and $A-\dot{\varpi }_d$. Then the maximum eccentricity is given by

$$\begin{equation} e_{\rm char}=2\left(\left|\frac{B_b}{A}\right|+ \left|\frac{B_d}{A-\dot{\varpi }_d}\right|\right). \end{equation} \tag{ 40 }$$

Comparing this expression with Equation (25), we conclude that disk precession does not affect planetesimal eccentricity behavior as long as $|\dot{\varpi }_d|\lesssim |A|$.

However, in the opposite case of $|\dot{\varpi }_d|\gtrsim |A|$ the disk-driven forced part of the eccentricity vector is suppressed compared to the case of no precession. This is because rapid precession of the disk (compared to the rate of planetesimal orbital precession) effectively averages out the non-axisymmetric part of the disk potential, considerably reducing the related eccentricity excitation. This has implications discussed in Section 7.2.2. At the same time, the forced eccentricity contribution due to the binary stays unchanged for planetesimals embedded in the precessing disk. We expect these asymptotic results to remain valid even in the case of non-uniform disk precession, both when it is much faster and much slower than |A|. However, all this discussion strictly applies only in the absence of gas drag.

We can put our findings in the context of existing results on the purely gravitational dynamics (i.e., not accounting for gas drag) of planetesimals in binaries. Heppenheimer (1978) explored planetesimal dynamics under the gravity of the companion alone. Our results reduce to his in the limit of a zero-mass disk, i.e., when planetesimal dynamics are in the BB regime; see Section 4.2.

It was first shown analytically in R13 that the gravity of a massive disk can significantly suppress planetesimal eccentricity excitation in binaries. The reason lies in the fast precession of planetesimal orbits caused by the disk gravity, which effectively averages out e_p forcing by the companion. This effect is present in our calculations as well and we reproduce the results of R13 in Case DB.

However, our work includes another important ingredient not considered previously in the framework of secular theory—gravitational forcing of planetesimal eccentricity by the disk itself, which should be present in addition to planetesimal precession if the disk is eccentric. While some numerical studies on this topic do exist (see Section 7.3), analytical understanding of their results has been hampered by the complexity of the problem.

In this work, we have provided the first (to the best of our knowledge) calculation of the eccentric disk potential in application to planetesimal dynamics. Using this prescription, we uncovered the existence of two entirely new regimes of planetesimal dynamics—Case BD (Section 4.3) and Case DD (Section 4.1)—in which eccentricity excitation by the disk exceeds that due to the secondary. The latter regime (DD) represents a very common situation in protoplanetary disks in binaries. As we have shown in Section 4, in many cases, planetesimal excitation is in the DD regime throughout the whole disk.

The significance of this dynamical regime also lies in the fact that the disk drives planetesimal eccentricities to a value of the order of the local disk eccentricity; see Equation (33). Although in the absence of any damping agents, the eccentricity of a particle starting on a circular orbit oscillates, see Equation (24), so that during some periods e_p ≪ e_char, most of the time e_p is of order e_char ∼ e_d in the regime DD. Thus, eccentricity of the disk gives rise to a lower limit on the characteristic planetesimal eccentricity (25), which is a very important finding.

In particular, it constrains the applicability of the axisymmetric disk approximation used in R13. Indeed, let us calculate e_p using Equation (37), which is identical to the result of R13, for a system with M_p = M_☉, ν = 0.3, e_b = 0.4, a_b = 20 AU harboring an axisymmetric disk with a_out = 5 AU, M_d = 10⁻² M_☉, p = 1. At a_p = 2 AU, we find e_p ≈ 10⁻², which is much less than it would be in a diskless case, $e_p^{\rm BB}\approx 0.1$; see Equation (34). However, for this result to hold in a non-axisymmetric disk, the disk eccentricity at 2 AU has to be less than 10⁻². Whether such a low e_d is realistic is not clear at the moment (see Section 7.3).

Our current results have been derived assuming that the disk affects planetesimals only via its gravitational field. In practice, planetesimals are also subject to gas drag, which has important consequences for their dynamics. First, gas drag lowers planetesimal velocities with respect to gas, which also lowers relative velocities between planetesimals, therefore positively affecting their survival in mutual collisions. Second, it has long been known that gas drag introduces apsidal alignment of planetesimal orbits (Marzari & Scholl 2000), which considerably reduces relative collision velocities of near-equal bodies. However, planetesimals of different sizes would still collide at high speeds suppressing growth (Thébault et al. 2006, 2008). Third, gas drag damps the free part of eccentricity; see Beaugé et al. (2010). This should affect the time dependence of planetesimal eccentricity, which in our case is given by Equation (24). We address the effects of gas drag on planetesimal dynamics in binaries in Rafikov & Silsbee (2015a).

Because of the neglect of gas drag, our current results are strictly valid only for relatively large objects, with sizes of several hundred kilometers. For such planetesimals, gas drag can be unimportant compared to purely gravitational forces during a rather long time span, and may thus be neglected. Inclusion of gas drag does not negate our finding that disk gravity from an eccentric disk leads to high encounter velocities between planetesimals, even of kilometer size. Our results also clearly show that purely gravitational effects alone, in the absence of dissipative forces, can give rise to non-trivial behavior of e_p (see, e.g., Sections 4.1 and 4.3) not captured in previous analyses of the problem.

We also note that our assumed surface density profile (1) and (2) may not fully capture the distribution of Σ in real disks. First, pressure forces drive differential precession in a hydrodynamical disk, which can be avoided only under rather special circumstances (Statler 2001). Second, these equations in their current form do not capture the possible presence of the density waves in the disk driven by the companion perturbation. They can be accounted for by assuming that the apsidal angle ϖ_d of the fluid trajectories varies with the distance in a particular fashion. For simplicity, we did not consider such a possibility in this work.

However, even if the expressions (1) and (2) are only approximate, this does not change our main conclusions about the key role of the disk gravity. Indeed, we find the values of A_d and B_d/e_g, which determine the disk effect on planetesimal dynamics, to not depend sensitively on the power-law indices p and q over a range of reasonable values; see Figure 4. Thus, we do not expect our results to change dramatically if the behavior of Σ(a_p) and e_d(a_p) were to deviate from the pure power laws in a_p.

Finally, short-term variability of the disk surface density can induce fast-changing torques on planetesimals. These effects cannot be captured by our secular (time-averaged) approach. However, we do not expect them to act coherently on long timescales; therefore they should be subdominant for the same reason that the short-period terms of the planetary perturbations play an insignificant role on long time intervals in classical celestial mechanics (Murray & Dermott 1999).

7.1. Implications for Planetesimal Growth

Planetesimal growth requires relative velocities of colliding bodies to be small, otherwise they get eroded or destroyed. We use our results to provide some insights on planetesimal accretion in binaries.

For bodies held together primarily by gravity, the threshold collision velocity at which planetesimals can still survive is about the escape speed. Guided by this logic, Moriwaki & Nakagawa (2004) and R13 use the simple criterion

$$\begin{eqnarray} e>e_{\rm {crit}} &=& \frac{2v_{\rm {esc}}}{v_k} \approx 6.1 \times 10^{-4} \frac{d}{10\;{\rm km}} \nonumber \\ & &\times \left(\frac{M_\odot }{M_p}\frac{a_p}{1\;{\rm AU}\;}\frac{\rho }{3\;{\rm g\; cm}^{-3}}\right)^{1/2} \end{eqnarray} \tag{ 41 }$$

as the condition for planetesimal destruction in collisions. Here v_esc is the escape speed from a planetesimal of a given radius d and ρ is the bulk density of planetesimal material.

In Figure 8, we display e_crit(a_p) with the dashed line and compare it with the characteristic e_p attained by planetesimals as a result of disk+secondary gravitational perturbation for different disk models. One can see that in all models where disk eccentricity e_d is high, comparable to the free-particle diskless eccentricity $e_p^{\rm BB}$ (Figures 8(a) and (b)), the disk does not help eliminate the fragmentation barrier. This is because e_p cannot drop below e_d and e_d is high. The situation is clearly more helpful for planetesimal growth in lower e_d cases; see Figures 8(c) and (d), even though it is still not as easy as in the case of the axisymmetric disk studied in R13. On the other hand, it has been noted in Rafikov (2013b) that the catastrophic destruction condition (41) is likely too conservative and underestimates the ability of planetesimals to survive in mutual collisions. This issue is addressed in more detail in Rafikov & Silsbee (2015b).

Another potential problem that may arise in low-mass disks with M_d ∼ 10⁻³ M_☉ is the presence of secular resonance in the disk; see Section 5 and Figure 8. There e_p becomes very large in a narrow range of a_p, making planetesimal collisions highly destructive. This phenomenon is non-local since high-e_p objects can penetrate other disk regions and destroy planetesimals there as well. However, this problem is unlikely to last for a long time as the small number of planetesimals from the vicinity of the secular resonance will be rapidly destroyed in collisions, leaving no more projectiles to destroy the remaining planetesimals in the rest of the disk.

7.2. Lowering Planetesimal Excitation

Motivated by our results and their implications for planetesimal accretion we next discuss different scenarios (in order of their likely significance) in which relative velocities of planetesimals affected only by the gravity of the gaseous disk and binary companion can be considerably lowered.

7.2.1. Intrinsically Low e_d

The major obstacle for planetesimal growth in high e_d disks has to do with our general result (Section 5) that e_p is always above the smaller of $e_p^{\rm BB}(a_p)$ and $e_p^{\rm DD}(a_p)\sim e_d$. Thus, one of the most straightforward ways of lowering collision speeds is for the disk to have low e_d either locally or globally for a long period of time. Our current understanding of eccentricity excitation in gaseous disks is based primarily on the results of numerical simulations, which are reviewed in Section 7.3. We describe possible ways of lowering e_d there.

7.2.2. Rapidly Precessing Disk

When discussing the possibility of disk precession in Section 6 we noted that the disk-induced contribution to the forced eccentricity can be effectively suppressed if the disk precesses faster than the planetesimals, i.e., if $|\dot{\varpi }_d|\gg |A|$. In this case, e_p can easily be below $e_p^{\rm DD}\sim e_d$. The remaining excitation due to the binary will keep e_p at the level of $e_p^{\rm DB}$, which is low because of the fast planetesimal precession driven by the massive disk; see Section 4.4. Thus, fast disk precession effectively brings planetesimal dynamics to the situation described in R13 and can serve as a mechanism for lowering planetesimal excitation as long as gas drag can be neglected (Rafikov & Silsbee 2015a). Whether the gaseous disk can precess at a rate exceeding |A| at separations of several AU, where the giant planets are detected in close binaries, should thus be explored in more detail.

7.2.3. Globally Suppressed Eccentricity Excitation

Another way of making $e_p^{\rm DD}$ low is hinted at by Equation (33), which shows that $e_p^{\rm DD}$ can be low even for high e_d if ψ₂ is very small. The same is true for $e_p^{\rm BD}$; see Equation (35). This is because the disk then produces zero contribution to the non-axisymmetric component of the disturbing function. Figure 4 shows that, ignoring the possible edge effects, this is possible, e.g., if p + q = −1. Our fiducial disk model, based on general ideas about the accretion disk physics and their eccentricity excitation, has p = 1 and q = −1, which is not compatible with this condition. However, our present understanding of protoplanetary disks in binaries does not allow us to exclude disk models with p + q = −1. Interestingly, when p = 0, then the disk also has zero contribution to the axisymmetric component; a disk with p = 0 (i.e., uniform disk), q = −1 affects neither planetesimal precession nor eccentricity excitation by its gravity in the absence of edge effects.

In Figure 9 we show an example of one such model having p = 0.5 and q = −1.5, with rather high disk eccentricity at the outer disk edge, e₀ = 0.25e_b = 0.05 (for e_b = 0.2). One can clearly see that, in this case, $e_p^{\rm DD}$ is low and comparable to e_crit at ∼AU separations. This should facilitate planetesimal growth on these scales. The disk-induced excitation for this model is not exactly zero due to edge effects. This is more of an issue near the outer edge of the disk because, as shown in Appendix A (and illustrated in Figure 3) p + q = −1 is closer to the line of convergence at the outer edge than at the inner edge.² This means that edge effects are more important in the outer disk for this set of power-law indices.

Zoom In Zoom Out Reset image size

It is also worth noting that some (though not all) of the lower planetesimal eccentricity in this figure as compared to Figure 8 is due to lower assumed disk eccentricity e_d in the inner part of the disk. However, the drop in e_p as one moves away from the inner edge of the disk reflects the drop in the non-axisymmetric part of the disk disturbing function as the inner edge effect becomes less important and we see that e_p drops well below the local value of e_d.

7.2.4. Locally Suppressed Eccentricity Excitation

Additionally, there are at least two ways in which e_p can be reduced locally, within a narrow range of semimajor axes. First, even if the disk does not have p + q = −1 globally, as we assumed in making Figure 9, there could be parts of the disk in which this condition is fulfilled for a range of a_p, for example, near the disk edges, where Σ_p should be petering out to zero, or near dead zones or opacity transitions where the material pileup is possible and a non-power-law scaling of Σ_p is likely. Our results do not directly apply to such situations since we assumed a purely power-law behavior of Σ_p(a_p) but based on them, we can expect that it might be possible to have ψ₂ close to zero at radii, near which locally computed p + q = −∂ln (e_dΣ)/∂ln a_p passes through −1 (edge effects mentioned in Section 7.2.3 may make the situation even more complicated). At this location, contributions of the inner and outer disks to ψ₂ should nearly cancel each other out, resulting in low $e_p^{\rm DD}$. Of course, e_p is lowered in this way only if the disk dominates eccentricity excitation, i.e., in Case DD.

7.2.5. Favorable Disk–Binary Orientation

Second, so far, we have always assumed planetesimal eccentricity to be given by the characteristic value e_char defined by the Equation (25). This approach ignores the dependence of the actual maximum planetesimal eccentricity upon the relative disk–binary orientation, obvious from Equation (24). In particular, in Section 3.4 we noted that whenever the conditions (26) are fulfilled, the maximum eccentricity is much lower than e_char. Because of the different dependence of B_d and B_b on a_p the first condition can be fulfilled only locally, within a narrow range of radii around a_p = a_B given by Equation (32). Since a_B lies within the disk only for relatively small M_d (see Equation (32)), we conclude that the first condition is fulfilled only for relatively light disks, M_d ≲ 10⁻² M_☉.

For most disk models considered in this work, one finds ψ₂ > 0 (see Figure 4) and B_d > 0 (Equation (6)), while B_b < 0 (Equation (10)). Then the second condition in (26) implies ϖ_d ≈ 0, i.e., that the binary and the disk apsidal lines need to be aligned for e_p to be suppressed at a_B. For the more atypical cases with ψ₂ < 0, one finds that the disk–secondary anti-alignment (ϖ_d ≈ π) is necessary to suppress e_p at a_B.

The actual value of ϖ_d for disks inside binaries is not well understood and Okazaki et al. (2002) find numerically that both alignment and anti-alignment are possible for the disks stationary in the binary frame. Needless to say, if the disk is precessing, it is no longer possible for it to be aligned or anti-aligned with the binary companion for a long time and the conditions (26) are no longer relevant.

7.3. Comparison with Numerical Studies

In this work, we explored the secular dynamics of planetesimals embedded in an eccentric gaseous disk, with implications for planet formation in binaries. We derived, for the first time, the analytical expression for the disturbing function of a body subject to gravity of a massive, eccentric, confocal, and coplanar disk in the limit when both the disk and planetesimal eccentricities are small (Appendix A). This expression has been used in Section 3.4 to understand secular excitation of e_p in the presence of both the non-axisymmetric disk and the binary companion. Assuming initially circular orbits and neglecting any dissipation (such as due to gas drag) in this work, we found the general analytical solution for the evolution of planetesimal eccentricity—Equation (24)—which shows that e_p oscillates from zero up to a maximum value.

Both the period and the amplitude of oscillations depend on properties of the disk and the secondary. Depending on which agent—disk or secondary—dominates planetesimal precession and eccentricity excitation, we find four distinct regimes for the e_p behavior. Two of them, in which the gravity of the eccentric disk dominates planetesimal eccentricity excitation, are novel results of this work. We have shown, in particular, that when the disk dominates both planetesimal precession and eccentricity excitation (the so-called Case DD; see Section 4.1), the characteristic planetesimal eccentricity e_p is of the order of the local disk eccentricity e_d. Thus, the value of e_d sets a lower limit on e_p and essentially determines the characteristic collision speeds of planetesimals. As a result, we generally find that eccentricity of the disk presents a serious obstacle for the growth of planetesimals with sizes of less than several tens of kilometers.

We then discuss possible ways of lowering e_p, which would be favorable for planetesimal growth (Section 7.2). One of them is for the disk to be massive, typically ≳ 10⁻² M_☉, so that (1) its own self-gravity reduces the disk eccentricity e_d as has been suggested by some simulations and (2) the disk gravity dominates planetesimal dynamics. Another possibility is for the disk to precess much faster than the precession rate of planetesimal orbits (Section 6). Some other ways of lowering e_p, both global and local (within a finite range of separations) are also described. These possibilities may represent pathways to planetesimal growth in at least a subset of protoplanetary disks in binary systems.

Despite the neglect of dissipative effects such as gas drag (accounted for in Rafikov & Silsbee 2015a, 2015b) the present study demonstrates the variety of planetesimal dynamical behaviors driven by the coupled gravitational perturbations of an eccentric disk and the binary. It thus represents an important step in building a complete picture of planetesimal dynamics in binaries.

The analytical description of the gravitational effects of the eccentric disk derived in this work (Appendix A) can be applied to a variety of other astrophysical problems: planetesimal dynamics in circumbinary disks (K. Silsbee & R. R. Rafikov, in preparation), dynamics of self-gravitating gaseous and stellar disks, and so on.

Here we present a calculation of the disturbing function due to an eccentric disk. We assume that the disk eccentricity and surface density are given by the power-law ansatz (2) and the apsidal angle is constant with radius. The latter assumption can be easily relaxed and analytical results can be obtained for ϖ_d varying as a power law of the semimajor axis of a fluid element.

There are different ways in which such a calculation can be approached. In particular, one can use an analogy with the Gauss averaging method (Murray & Dermott 1999), which treats the time-averaged potential of a point mass on an eccentric orbit as that produced by an elliptical wire along the orbit with the line density proportional to the time the planet spends at each point of its orbit. In the case of a gaseous disk, we can consider fluid in a narrow elliptical annulus between the two adjacent fluid trajectories. Because of the continuity equation, the line density of this fluid along the annulus is also proportional to the time fluid spends at a given location. Given that the density distributions are the same in two cases, one can simply employ the expression for the disturbing function given by the Gauss method. For example, secular contribution due to the outer disk becomes

$$\begin{eqnarray} R^{\rm Gauss}_{{\rm out}} &=& n_p^2 a_p^2 \int _{a_p}^{a_{\rm out}} \frac{2 \pi a \Sigma _p(a)}{M_p} \left[\frac{\alpha ^2}{8} b_{3/2}^{(1)}\left(\alpha \right) e_p^2\right. \nonumber\\ &&\left. - \frac{\alpha ^2}{4} b_{3/2}^{(2)}\left(\alpha \right) e_p e_d(a) \cos {(\varpi _p - \varpi _d)}\right]da, \end{eqnarray} \tag{ A1 }$$

where α = a_p/a. A similar expression can be written for the inner part of the disk as well. However, both $\int _1 b_{3/2}^{(1)}(\alpha) d \alpha$ and $\int _1 b_{3/2}^{(2)}(\alpha) d\alpha$ are non-convergent, as well as the sum of the inner and outer disk contributions in the vicinity of planetesimal orbit. This is a well-known problem with the Gauss expression for the secular disturbing function (Murray & Dermott 1999). For this reason, we are unable to use Gauss's method to calculate the disturbing function due to an eccentric disk for a planetesimal that is embedded in a disk with no gap.

Instead, we have resorted to a different approach previously used by Heppenheimer (1980) and Ward (1981) to compute the gravitational field of an axisymmetric disk with a power-law surface density profile. To use this approach for an elliptical disk, we had to come up with a number of important modifications. The idea behind this method is to compute the disturbing function directly as

$$\begin{eqnarray} &&R({\bf S})=G \left\langle \int _{\bf S} \frac{\Sigma (r_d, \phi _d)r_d dr_d d\phi _d}{\left(r_p^2+r_d^2 - 2r_pr_d\cos \theta \right)^{1/2}}\right\rangle, \end{eqnarray} \tag{ A2 }$$

where the integral is taken over the area of the disk S, angle brackets 〈...〉 represent the time averaging over planetesimal orbital motion, r_p is the (time-dependent) instantaneous radius of a planetesimal, and θ is the angle between vectors r_d and r_p; see Figure 1. According to this figure, ϕ_d is the polar angle counted from the disk periastron, ϕ_p is the angle of the planetesimal with respect to the planetesimal periastron, so that θ = ϕ_d + ϖ_d − ϕ_p − ϖ_p.

We divide the disk up into three regions as shown in Figure 10, so that S = S_c + S₀ − S_i. Here S_c is the annulus bounded by circles with radii equal to the periastron of the outer disk edge a_out[1 − e_d(a_out)] and the periastron of the inner disk edge a_in[1 − e_d(a_in)]; S_o is the outer crescent region bounded by the outer circle of S_c on the inside and the outermost elliptical trajectory on the outside; S_i is the inner crescent region bounded by the inner circle of S_c on the outside and the innermost elliptical trajectory on the inside. The full disturbing function of an eccentric disk is given by

$$\begin{equation} R({\bf S}) = R({\bf S}_c) + R({\bf S}_o) - R({\bf S}_i) . \end{equation} \tag{ A3 }$$

We now separately calculate the contributions due to different regions using an extension of the method employed by Heppenheimer (1980).

Zoom In Zoom Out Reset image size

A.1. Contribution from the Annular Region ${\boldsymbol S}_c$

We start by calculating the contribution from the annular region S_c. In the following, we define for brevity a_in = a₁, a_out = a₂, e_d(a_in) = e₁, e_d(a_out) = e₂, with r_{d, in} = a₁(1 − e₁) and r_{d, out} = a₂(1 − e₂) being the inner and outer radii of S_c. We can write

$$\begin{eqnarray} &&R({\bf S}_c)=\left \langle G \int ^{r_{d,{\rm out}}}_{r_{d,{\rm in}}} r_d dr_d \int _{0}^{2\pi } \frac{\Sigma (r_d, \phi _d)}{\sqrt{r_p^2+r_d^2 - 2r_pr_d\cos \theta }}d\theta \right \rangle ,\quad \end{eqnarray} \tag{ A4 }$$

where we have used the fact that dϕ_d = dθ.

As given in Statler (2001) and using Equation (1), to first order in e_d,

$$\begin{eqnarray} \Sigma (r_d, \phi _d) &=& \Sigma _p(r_d)+ e_d\left[\zeta (r_d)\Sigma _p(r_d)\left(\cos \phi _d - 1\right)\right.\nonumber\\ &&\left. + r_d \cos \phi _d \frac{d\Sigma _p(r_d)}{dr_d}\right] , \end{eqnarray} \tag{ A5 }$$

where ζ(r_d) was defined after Equation (1). Note that Σ_p is considered to be a function of the semimajor axis of a fluid element passing through a given point in the disk, as stated after Equation (1). For that reason Σ(r_d, 0) ≠ Σ_p(r_d) but Σ(r_d, 0) = Σ_p(r_d/(1 − e_d)) (to second order in e_d), i.e., at the semimajor axis r_d/(1 − e_d) for which r_d is the periastron distance.

Classical secular theory neglects terms in the disturbing function that are higher order than $e_p^2$ in planetesimal eccentricity; see Section 3.1. Thus, in our subsequent calculations we will retain only terms proportional to $e_p^2$ and e_pe_d; terms of higher order in e_d are neglected because, by assumption, e_d ≪ 1. As we will see below, terms with no ϕ_d dependence lead to corrections of the order of $e_p^2$. Therefore, we can drop the terms with no ϕ_d dependence, which are also proportional to e_d.

With this in mind, we write the contribution to the disturbing function from the annular component S_c as

$$\begin{eqnarray} R({\bf S}_c) & = & I_1 + I_2, \end{eqnarray} \tag{ A6 }$$

$$\begin{eqnarray} I_1 & = & \left\langle G \int _{r_{d,{\rm in}}}^{r_{d,{\rm out}}} dr_d \Sigma _p(r_d) r_d \int _{0}^{2\pi } \frac{d\theta }{\sqrt{r_p^2+r_d^2-2r_pr_d\cos {\theta }}} \right\rangle\qquad \end{eqnarray} \tag{ A7 }$$

$$\begin{eqnarray} I_2 & = & \left\langle G \int _{r_{d,{\rm in}}}^{r_{d,{\rm out}}} dr_d \left[\zeta (r_d)\Sigma _p(r_d) + r_d \frac{d\Sigma _p(r_d)}{dr_d}\right]e_d(r_d) r_d\right.\nonumber\\ &&\times\left. \int _{0}^{2\pi } \frac{\cos {(\theta + v)} d\theta }{\sqrt{r_p^2+r_d^2-2r_pr_d\cos {\theta }}} \right\rangle. \end{eqnarray} \tag{ A8 }$$

Here we expressed ϕ_d = θ + v, where v = ϖ_p − ϖ_d + ϕ_p; see Figure 1. We now evaluate these two contributions.

Evaluation of I₁. From the definition (11) of the Laplace coefficients, we can write the inner integral over θ in Equation (A7) as $(\pi /r_d)b_{1/2}^{(0)}(r_p/r_d)$ for r_p < r_d (outer disk) and $(\pi /r_p)b_{1/2}^{(0)} (r_d/r_p)$ for r_p > r_d (inner disk). Assuming a surface density prescription (2), we can write

$$\begin{eqnarray} I_1 & = & \pi G \Sigma _0\left\langle \int _{r_{d,{\rm in}}}^{r_p} \left(\frac{a_{\rm out}}{r_d}\right)^p \frac{r_d}{r_p} b_{1/2}^{(0)}\left(\frac{r_d}{r_p}\right)dr_d\right.\nonumber\\ &&\left.+ \int _{r_p}^{r_{d,{\rm out}}} \left(\frac{a_{\rm out}}{r_d}\right)^p b_{1/2}^{(0)}\left(\frac{r_p}{r_d}\right)dr_d \right\rangle. \end{eqnarray} \tag{ A9 }$$

We define an auxiliary function

$$\begin{equation} I(x, y, z) \equiv \int _x^1 \alpha ^y b_{1/2}^{(z)}(\alpha) d\alpha, \end{equation} \tag{ A10 }$$

and a new constant factor

$$\begin{equation} K = \pi G \Sigma _0 a_{\rm out}^p a_p^{1-p}. \end{equation} \tag{ A11 }$$

With these definitions we re-write expression (A9) as

$$\begin{equation} I_1 = K\left\langle \left(\frac{r_p}{a_p}\right)^{1-p}\left[I(a_1/r_p, 1-p, 0) + I(r_p/a_2, p-2, 0) \right] \right\rangle. \end{equation} \tag{ A12 }$$

We note that a₁ and a₂ in these expressions approximate r_{d, in} = a₁(1 − e₁) and r_{d, out} = a₂(1 − e₂), respectively. However, the difference is a correction linear in disk eccentricity e_d and should be ignored for I₁.

We now proceed to the last, time averaging, step. For illustration, we perform it first on the second integral in this expression by expanding it in a Taylor series in a small quantity r₂ − α₂, where r₂ = r_p/a₂, and α₂ = a_p/a₂. We have

$$\begin{eqnarray} I(r_p/a_2, p-2, 0)& =& I(\alpha _2, p-2, 0) - (r_2 - \alpha _2) \alpha _2^{p-2} b_{1/2}^{(0)}(\alpha _2)\nonumber\\ && - \frac{(r_2 - \alpha _2)^2}{2} \frac{d}{d\alpha _2}\left[\alpha _2^{p-2} b_{1/2}^{(0)}(\alpha _2)\right]. \end{eqnarray} \tag{ A13 }$$

We may relate r_p and a_p using the eccentric anomaly E as r_p = a_p(1 − e_pcos E). Then r₂ − α₂ = −α₂e_pcos E and

$$\begin{equation} \left(\frac{r_p}{a_p}\right)^{1-p} = 1 - (1-p)e_p \cos E - \frac{p(1-p)}{2} e_p^2\cos ^2 E.\quad \end{equation} \tag{ A14 }$$

Using these relations, the second integrand in (A12) becomes (retaining only terms up to $e_p^2$)

$$\begin{eqnarray} && K \left\langle 1 - (1-p)e_p \cos E - \frac{p(1-p)}{2} e_p^2\cos ^2 E \right\rangle I(\alpha _2, p-2, 0) \nonumber \\ &&\quad + K \left\langle [1 - (1-p)e_p \cos E ] e_p \cos E \alpha _2^{p - 1} b_{1/2}^{(0)}(\alpha _2)\right\rangle\nonumber\\ &&\quad -\frac{K}{2} \left\langle \alpha _2^2 e_p^2 \cos ^2 E \frac{\partial }{\partial \alpha _2} \left[\alpha _2^{p - 2} b_{1/2}^{(0)}(\alpha _2)\right] \right\rangle. \end{eqnarray} \tag{ A15 }$$

Using 〈cos E〉 = −e_p/2 and 〈cos ²E〉 = 1/2, Equation (A15) reduces to

$$\begin{eqnarray} &&K\left[1+\frac{e_p^2}{4}(1-p)(2-p) \right] I(\alpha _2, p-2, 0) + K\frac{e_p^2}{4}\nonumber\\ &&\quad\times\left[2(p-1)\alpha _2^{p-1}b_{1/2}^{(0)}(\alpha _2) - \frac{d}{d\alpha _2} \left[\alpha _2^p b_{1/2}^{(0)}(\alpha _2)\right] \right].\qquad \end{eqnarray} \tag{ A16 }$$

We can apply the identical procedure to the first integral of Equation (A12), resulting in

$$\begin{eqnarray} &&K\left[1+\frac{e_p^2}{4}(1-p)(2-p) \right] I(\alpha _1, 1-p, 0) +K\frac{e_p^2}{4}\nonumber\\ &&\quad\times \left[2(2-p) \alpha _1^{2-p} b_{1/2}^{(0)}(\alpha _1) - \frac{\partial }{\partial \alpha _1}\left[\alpha _1^{3-p}b_{1/2}^{(0)}(\alpha _1)\right]\right],\qquad \end{eqnarray} \tag{ A17 }$$

where α₁ = a₁/a_p. Note that the second order term in e_p must be included in $r_1 - \alpha _1 = \alpha _1 e_p \cos {E} + \alpha _1 e_p^2 \cos ^2{E}$, where r₁ = a₁/r_p. The sum of Equations (A16) and (A17) is equal to I₁ and represents the axisymmetric part of the disk disturbing function from the region S_c.

Calculation of I₂. In order to calculate I₂—the non-axisymmetric component of R(S_c)—we use the prescription (2) for Σ_p and e_d, and expand cos (θ + v):

$$\begin{eqnarray} I_2 & = & -G(p+q)\left\langle \int _{r_{d,{\rm in}}}^{r_{d,{\rm out}}} dr_d r_d \Sigma _0 e_0\left(\frac{a_{\rm out}}{r_d}\right)^{p+q}\right.\nonumber\\ &&\times\left. \int _{0}^{2\pi } \frac{\cos {\theta }\cos {v} - \sin {\theta } \sin {v}}{\sqrt{r_p^2+r_d^2-2r_pr_d\cos {\theta }}} d\theta \right\rangle. \end{eqnarray} \tag{ A18 }$$

In the inner integral over θ, the terms with sin θ in the numerator vanish upon integration, while the terms with cos θ result in Laplace coefficients $b_{1/2}^{(1)}$; see definition (11). Separately accounting for the contributions from the inner and outer disks when integrating over r_d, we obtain, analogous to Equation (A12)

$$\begin{eqnarray} I_2 &=& -Ke_d(a_p)(p + q)\left\langle \cos {v}\left(\frac{r_p}{a_p}\right)^{1-p-q}\right. \nonumber\\ &&\left. \times\left[ I(r_p/a_2, p + q - 2, 1)+I(a_1/r_p, 1\!-\!p\!-\!q, 1) \right] \!\right\rangle.\qquad \end{eqnarray} \tag{ A19 }$$

The final step of time-averaging is somewhat more challenging here because the cos v term introduces additional time dependence through ϕ_p. It can be taken care of using the definition v = (ϖ_p − ϖ_d) + ϕ_p and the relation ϕ_p = E + e_psin E accurate to linear order in e_p. As before, we also expand integrals in (A19) in a series in the small quantities r_{1, 2} − α_{1, 2}. Since we are not interested in the terms $O\left(e_d e_p^2\right)$ and higher order (small factor e_d is already present in Equation (A19)), we only expand to first order in e_p. As a result of tedious but straightforward calculation, we find

$$\begin{eqnarray} I_2&=& -Ke_d(a_p)e_p\cos {(\varpi _p - \varpi _d)}(p + q) \nonumber\\ &&\times \left [\frac{(p+q-3)}{2}I(\alpha _1, 1-p - q, 1) - \frac{1}{2} \alpha _1^{2-p-q}b_{1/2}^{(1)}(\alpha _1)\right. \nonumber \\ &&\left. + \frac{(p+q-3)}{2}I(\alpha _2, p + q \!-\! 2, 1) \!+\! \frac{1}{2} \alpha _2^{p+q-1}b_{1/2}^{(1)}(\alpha _2)\!\right ]. \qquad \end{eqnarray} \tag{ A20 }$$

This completes our calculation of R(S_c).

A.2. Contribution from the Inner Crescent ${\boldsymbol S}_i$

We now calculate the disturbing function R(S_i), given by Equation (A2) with integration carried out over the inner crescent S_i. The width of the crescent is O(e_d), meaning that we need to keep all variables only up to first order in e_p. The integrand, which led to an axisymmetric contribution in the case of R(S_c) now leads to a non-axisymmetric contribution when integrated over this non-axisymmetric region of the disk.

Consider an ellipse with a periastron distance a_{p, 1} = a₁(1 − e₁) and an apoastron distance a_{a, 1} = a₁(1 − e₁) bounding S_i on the outside. Define the angle ξ(r_d) as the angle between the periastron of this ellipse and the point of intersection of the ellipse and a circle of radius r_d, a_{p, 1} < r_d < a_{a, 1}. Linearizing the equation of an ellipse $r_d = a_1(1-e_1^2)/\left[1+e_1\cos \xi (r_d)\right]$ in e₁, we get r_d = a₁(1 − e₁cos ξ(r_d)). This yields

$$\begin{equation} \xi (r_d) = \arccos \frac{a_1-r_d}{e_1a_1}, \end{equation} \tag{ A21 }$$

where the arccos function is the inverse cosine function. We write explicitly

$$\begin{eqnarray} R({\bf S}_i) &=& \left\langle G \Sigma _0 a_{\rm out}^{p} \int _{a_{p,1}}^{a_{a,1}} \frac{r_d^{1-p}}{r_p} dr_d \int _{\xi -\Delta \varpi - \phi _p}^{2\pi -\xi -\Delta \varpi - \phi _p}\right.\nonumber\\ &&\times\left. \frac{d\theta }{\sqrt{1+ \alpha ^{\prime 2} - 2\alpha ^\prime \cos \theta }} \right\rangle, \end{eqnarray} \tag{ A22 }$$

where α' = r_d/r_p ≈ a₁/r_p, and Δϖ = ϖ_p − ϖ_d. Then using the relation

$$\begin{equation} \left(1 + \alpha ^{\prime 2} - 2\alpha ^\prime \cos {\theta }\right)^{-1/2} = \frac{1}{2} b_{1/2}^{(0)}(\alpha ^\prime) + \sum _{j = 1}^\infty b_{1/2}^{(j)}(\alpha ^\prime) \cos {(j \theta)} \end{equation} \tag{ A23 }$$

the inner integral over θ becomes

$$\begin{eqnarray} &&\left[\pi - \xi (r_d)\right] b_{1/2}^{(0)}(\alpha ^\prime) - \sum _{j = 1}^\infty \frac{2}{j} b_{1/2}^{(j)}(\alpha ^\prime) \nonumber\\ &&\quad\times \sin \left[j\xi (r_d)\right]\cos \left[j(\Delta \varpi + \phi _p)\right]. \end{eqnarray} \tag{ A24 }$$

Then we may write

$$\begin{eqnarray} R({\bf S}_i) & = & G \Sigma _0 a_{\rm out}^p \left \langle r_p^{-1}\int _{a_{p,1}}^{a_{a,1}} r_d^{1-p} dr_d \left [b_{1/2}^{(0)}\left(\frac{a_1}{r_p}\right)\left[\pi - \xi (r_d)\right] \right.\right.\nonumber \\ & &\left.\left.- \sum _{j = 1}^\infty \frac{2}{j} b_{1/2}^{(j)}\left(\frac{a_1}{r_p}\right) \sin \left[j\xi (r_d)\right] [\cos {(j\Delta \varpi)}\cos {(j\phi _p)}\right.\right.\nonumber\\ &&\left.\left. - \sin {(j\Delta \varpi)} \sin {(j\phi _p)}] \right ] \right \rangle. \end{eqnarray} \tag{ A25 }$$

We will use the following definite integrals

$$\begin{eqnarray} \int _{a_{p,1}}^{a_{a,1}} \left [\pi - \xi (r_d)\right ]dr_d &=& \pi a_1 e_1,\nonumber\\ \int _{a_{p,1}}^{a_{a,1}} \sin \left[\xi (r_d)\right] dr_d &=& \frac{\pi }{2} e_1 a_1,\nonumber\\ \int _{a_{p,1}}^{a_{a,1}} \sin \left[j\xi (r_d)\right] dr_d &=& 0 \end{eqnarray} \tag{ A26 }$$

for the integer j > 1. Then in (A25) we may ignore the terms in the sum with j > 1:

$$\begin{eqnarray} R({\bf S}_i) &=&\pi G \Sigma _0 \left(\frac{a_{\rm out}}{a_1}\right)^p e_1 a_1^2 \left\langle r_p^{-1} \left[b_{1/2}^{(0)}\left(\frac{a_1}{r_p}\right) - b_{1/2}^{(1)}\left(\frac{a_1}{r_p}\right)\right.\right.\nonumber\\ &&\left.\left.\times \left[\cos {\Delta \varpi } \cos {\phi _p} - \sin {\Delta \varpi } \sin {\phi _p}\right] \right] \right\rangle, \end{eqnarray} \tag{ A27 }$$

where e₁ and a₁ are the disk eccentricity and semimajor axis, respectively, evaluated at the inner edge. Using a₁/r_p = α₁ + e_pα₁cos E and the relation between ϕ_p and E, this becomes

$$\begin{eqnarray} R({\bf S}_i) &= & \pi G \Sigma _0 \left(\frac{a_{\rm out}}{a_1}\right)^p a_1 e_1 \left \langle (\alpha _1 + e_p \alpha _1 \cos {E}) \left \lbrace b_{1/2}^{(0)}(\alpha _1)\right.\right. \nonumber\\ &&\left.\left.+ e_p \alpha _1 \cos {E} \frac{ \partial b_{1/2}^{(0)}}{\partial \alpha _1} - \left[b_{1/2}^{(1)}(\alpha _1) + e_p \alpha _1 \cos {E} \frac{ \partial b_{1/2}^{(1)}}{\partial \alpha _1}\right]\right. \right. \nonumber \\ &&\left.\left. \times (\cos \Delta \varpi \cos E - \sin \Delta \varpi \sin E- e_p \cos {\Delta \varpi }\sin ^2{E}\right.\right.\nonumber\\ &&\left.\left. - e_p \sin {\Delta \varpi } \cos {E} \sin {E} )\right \rbrace \right \rangle. \end{eqnarray} \tag{ A28 }$$

Expanding all products in this expression, one gets a total of 20 terms. It is straightforward to angle-average them as before. Keeping only terms of order O(e₁e_p) and substituting e₁ = e₀(a_out/a₁)^q, we find that the disturbing function from the inner crescent is given by

$$\begin{eqnarray} R({\bf S}_i) &=&\frac{1}{2} K e_d(a_p)e_p \cos \left(\varpi _p - \varpi _d\right) \alpha _1^{2-p - q}\nonumber\\ &&\times\left[b_{1/2}^{(1)}(\alpha _1) - \alpha _1\frac{\partial b_{1/2}^{(1)}}{\partial \alpha _1}\right]. \end{eqnarray} \tag{ A29 }$$

A.3. Contribution from the Outer Crescent ${\boldsymbol S}_o$

The derivation of R(S_o) follows the same basic concept as that of R(S_i) except that now α' = r_p/r_d. As a result, one finds the contribution of the outer crescent S_o to be given by

$$\begin{eqnarray} R({\bf S}_o)&=&\frac{1}{2} K e_d(a_p)e_p\cos \left(\varpi _p - \varpi _d\right) \alpha _2^{p + q - 1}\nonumber\\ &&\times\left[2 b_{1/2}^{(1)}(\alpha _2) + \alpha _2 \frac{\partial b_{1/2}^{(1)}}{\partial \alpha _2}\right]. \end{eqnarray} \tag{ A30 }$$

A.4. Putting Everything Together

Plugging Equations (A16), (A17), (A20), (A29), (A30) into the expression (A3), we find that the total eccentric disk-induced disturbing function, including the eccentricity independent term and terms proportional to $e_p^2$ and e_de_p, is given by

$$\begin{equation} R = K\left[\psi _0 + \psi _1 e_p^2 + \psi _2 e_d(a_p)e_p \cos (\varpi _p - \varpi _d)\right] , \end{equation} \tag{ A31 }$$

with

$$\begin{eqnarray} \psi _0(\alpha _1,\alpha _2) & = & I(\alpha _1, 1-p, 0) + I(\alpha _2, p - 2, 0) \end{eqnarray} \tag{ A32 }$$

$$\begin{eqnarray} \psi _1(\alpha _1,\alpha _2) & = & \frac{1}{4}\left[(1-p)(2-p) \psi _0(\alpha _1,\alpha _2)\right.\nonumber\\ &&\left. + 2(p-1)\alpha _2^{p-1}b_{1/2}^{(0)}(\alpha _2) - \frac{d}{d\alpha _2} \left(\alpha _2^p b_{1/2}^{(0)}(\alpha _2)\right)\right. \nonumber \\ & &\left.+ 2(2-p) \alpha _1^{2-p} b_{1/2}^{(0)}(\alpha _1) - \frac{\partial }{\partial \alpha _1}\left(\alpha _1^{3-p}b_{1/2}^{(0)}(\alpha _1)\right)\right]\nonumber\\ \end{eqnarray} \tag{ A33 }$$

$$\begin{eqnarray} \psi _2(\alpha _1,\alpha _2) & = & -\frac{(p + q)(p+q-3)}{2} [I(\alpha _1, 1-p - q, 1)\nonumber \\ && + I(\alpha _2, p + q - 2, 1)] + \frac{\alpha _1^{2 - p - q}}{2} \nonumber\\ &&\times\left[\left(p+q-1\right)b_{1/2}^{(1)}(\alpha _1) + \alpha _1\frac{\partial b_{1/2}^{(1)}}{\partial \alpha _1}\right] + \frac{\alpha _2^{p + q - 1}}{2}\nonumber\\ &&\times \left[\left(2-p-q\right) b_{1/2}^{(1)}(\alpha _2) + \alpha _2 \frac{\partial b_{1/2}^{(1)}}{\partial \alpha }\right]. \end{eqnarray} \tag{ A34 }$$

This completes our calculation of the disturbing function due to an eccentric disk with properties given by Equation (2).

A.5. Asymptotic Behavior

Astrophysical disks typically span several orders of magnitude in radius. It is then plausible that far from the disk boundaries, we can ignore edge effects, i.e., the expression for the disturbing function does not depend on the a_in and a_out as a_in → 0 and a_out → ∞. This corresponds to the limit of α_{1, 2} → 0. Using a Taylor expansion,

$$\begin{eqnarray} &&b_{1/2}^{(0)}(\alpha) = 2 + \frac{\alpha ^2}{2},\quad b_{1/2}^{(1)}(\alpha) = \alpha + \frac{3}{8} \alpha ^3, \end{eqnarray} \tag{ A35 }$$

for small α in Equations (A33) and (A34), we determined that ψ₁ is convergent and independent of α_{1, 2} as α_{1, 2} → 0 goes to zero for −1 < p < 4. Similarly, ψ₂ is convergent as α_{1, 2} → 0 for −2 < p + q < 5. Convergence limits are illustrated in Figure 3.

Provided that the disturbing function is dominated by the local parts of the disk (i.e., the values of p and q fall within the white region in Figure 3) and the values of coefficients ψ_{1, 2} are independent of α_{1, 2} when the disk edges are well separated from the planetesimal semimajor axis (α_{1, 2} → 0), we can obtain simpler analytical expressions for these coefficients. Indeed, using the fact that $b_{1/2}^{(0)}(\alpha)=(4/\pi){\bf K}(\alpha)$, $b_{1/2}^{(1)}(\alpha)=(4/\pi \alpha)\left[{\bf K}(\alpha)-{\bf E}(\alpha)\right]$ (here E and K are complete elliptic integrals) and series expansions (Gradshteyn & Ryzhik 1994)

$$\begin{eqnarray} {\bf K}(\alpha)&=&\frac{\pi }{2}\left(\!1\!+\!\sum \limits _{n=1}^{\infty }A_n\alpha ^{2n}\!\right), {\bf E}(\alpha)=\frac{\pi }{2}\left(\!1\!-\!\sum \limits _{n=1}^{\infty }\frac{A_n}{2n-1}\alpha ^{2n}\!\right),\nonumber\\ A_n&=&\left[\frac{(2n)!}{2^{2n}(n!)^2}\right]^2, \nonumber\\ \end{eqnarray} \tag{ A36 }$$

we can provide asymptotic expressions for ψ_{1, 2} as follows:

$$\begin{eqnarray} \psi _1(0,0) & = & -\frac{1}{2}+\frac{(1-p)(2-p)}{2}\nonumber\\ &&\times \sum \limits _{n=1}^{\infty }\frac{(4n+1)A_n}{(2n+2-p)(2n+p-1)}, \end{eqnarray} \tag{ A37 }$$

$$\begin{eqnarray} \psi _2(0,0) & \rightarrow & \frac{3}{2}-(p+q)(p+q-3) \nonumber\\ &&\times \sum \limits _{n=2}^{\infty }\frac{2n(4n-1)A_n}{(2n-1)(2n+1-p-q)(2n-2+p+q)}. \nonumber\\ \end{eqnarray} \tag{ A38 }$$

The behavior of ψ₁(0, 0) and ψ₂(0, 0) as functions of p and p+q, respectively, are shown in Figure 4.

Here we describe the details of the numerical verification of our analytical results; see Section 3.5. We directly integrated orbits of planetesimals affected by the gravity of an eccentric disk using the MERCURY package (Chambers 1999). All our integrations employed the Bulirsch–Stoer algorithm (Press et al. 1992). Accelerations due to the gravity of an eccentric disk g_d (used as an input for our integrations) were computed at different positions and for different disk parameters via direct numerical integration as

$$\begin{eqnarray} &&{\bf g}_d({\bf r})=-G\int _{\bf S}\Sigma ({\bf r}_d) \frac{{\bf r}-{\bf r}_d}{|{\bf r}-{\bf r}_d|^3} d{\bf S}({\bf r}_d), \end{eqnarray} \tag{ B1 }$$

where dS(r_d) is a surface element centered on r_d. This two-dimensional integral was performed using standard integration by quadratures in SciPy. We used a small softening parameter in the integrand to better handle the singularity, and verified convergence to within a percent as we lowered the value of this parameter. The surface density of the eccentric disk was assumed to be given directly by Equation (1). In this calculation, we did not make an assumption of e_d ≪ 1 and thus were not expanding Equation (1) in powers of e_p (as opposed to Equation (A5)).

Here we explore the secular evolution of planetesimals in the case of a disk precessing according to a simple linear prescription $\varpi _d(t)=\varpi _{d0}+\dot{\varpi }_d t$. Plugging this into the Lagrange Equations (18) and (17), one finds the following solution for the components of the eccentricity vector (k_p, h_p) with initial conditions k_p(0) = 0, h_p(0) = 0:

$$\begin{eqnarray} k_p(t) & = & \frac{1}{A(A-\dot{\varpi }_d)} \lbrace AB_d\left[\cos (At+\varpi _{d0}) - \cos (\dot{\varpi }_d t + \varpi _{d0})\right]\nonumber\\ && + B_b(A-\dot{\varpi }_d)\left[\cos {(At)} - 1\right]\rbrace. \end{eqnarray} \tag{ C1 }$$

$$\begin{eqnarray} h_p(t) & = & \frac{1}{A(A-\dot{\varpi }_d)} \lbrace AB_d\left[\sin {(At+\varpi _{d0})} - \sin {(\dot{\varpi }_d t + \varpi _{d0})}\right]\nonumber\\ && + B_b(A-\dot{\varpi }_d)\sin {(At)}\rbrace, \end{eqnarray} \tag{ C2 }$$

which generalizes solution (19) to the case of non-zero precession. It can again be written as the sum of three distinct contributions as described in Section 6:

$$\begin{eqnarray} &&{\bf e}_p(t)=\left\lbrace\!\! \begin{array}{l}k_p(t)\\ h_p(t) \end{array} \!\!\right\rbrace = -\frac{B_b}{A} \left\lbrace\!\! \begin{array}{l}1 \\ 0 \end{array}\!\! \right\rbrace -\frac{B_d}{A-\dot{\varpi }_d} \left\lbrace\!\! \begin{array}{l}\cos \varpi _d(t) \\ \sin \varpi _d(t) \end{array}\!\! \right\rbrace \nonumber \\ &&\quad + \frac{\left[(AB_d)^2 + \left(B_b(A-\dot{\varpi }_d)\right)^2 + 2AB_dB_b(A - \dot{\varpi }_d)\cos \varpi _{d0}\right]^{1/2}}{A(A - \dot{\varpi }_d)}\nonumber\\ &&\quad \times \left\lbrace \!\!\begin{array}{l}\cos (At + \phi) \\ \sin (At + \phi) \end{array}\!\! \right\rbrace, \end{eqnarray} \tag{ C3 }$$

where ϕ is a phase defined analogous to Equation (23) and is a function of ϖ_d0, A, B_d and B_b. Note that in the case of a precessing disk, forced eccentricity due to the disk changes in time.