Secular Evolution Driven by Massive Eccentric Disks/Rings: An Apsidally Aligned Case

Our calculation of the orbit-averaged disturbing function R_d due to an eccentric disk uses the general mathematical procedure outlined in a seminal paper of Heppenheimer (1980) for the calculation of the disturbing function due to an axisymmetric disk. In this approach, the expansion of the disturbing function in terms of a small parameter—test particle eccentricity—proceeds differently from the classical Laplace–Lagrange theory (Murray & Dermott 1999). The resultant expression for R_d does not contain non-integrable singularities at the particle semimajor axis, resulting in a convergent expression for the disturbing function. In other words, this calculation does not require introduction of an ad hoc softening of the potential. This method was later used by Ward (1981) to study the stability of the early solar system perturbed by an axisymmetric protoplanetary disk.

Silsbee & Rafikov (2015b) extended the method of Heppenheimer (1980) to the case of nonaxisymmetric, apsidally aligned, eccentric disks with Σ given by Equation (2). However, their work was restricted to disks with power-law profiles of Σ_d(a) and e_d(a). Here we generalize this calculation even further to cover the disks with arbitrary behaviors of Σ_d(a) and e_d(a).

As a result of a rather lengthy derivation, the mathematical details of which are presented in Appendix A, we arrive at the following expression for the disturbing function due to an apsidally aligned disk:

$$\begin{eqnarray}&&{R}_{d}={a}_{p}^{2}{n}_{p}\left[\displaystyle \frac{1}{2}{A}_{d}{e}_{p}^{2}+{B}_{d}{e}_{p}\cos ({\varpi }_{p}-{\varpi }_{d})\right],\end{eqnarray} \tag{ 3 }$$

where n_p = (GM_c/a_p³)^1/2 is the particle mean motion, with the coefficients A_d and B_d (having dimensions of [s⁻¹] and discussed in more detail in Section 2.2) given by Equations (45)–(49). The key underlying simplifications making this calculation possible are that e_d ≪ 1, as well as d e_d/d ln a ≪ 1; see Equation (22).

Introducing a two-component eccentricity vector e_p =(k_p, h_p) = e_p(cosϖ_p, sin ϖ_p) for a test particle, as well as the auxiliary vector

$$\begin{eqnarray}&&{{\boldsymbol{B}}}_{d}={B}_{d}(\cos {\varpi }_{d},\sin {\varpi }_{d}),\end{eqnarray} \tag{ 4 }$$

Equation (3) can be rewritten as

$$\begin{eqnarray}{R}_{d} & = & {a}_{p}^{2}{n}_{p}\left[\displaystyle \frac{1}{2}{A}_{d}{{\boldsymbol{e}}}_{p}^{2}+{{\boldsymbol{e}}}_{p}\cdot {{\boldsymbol{B}}}_{d}\right]\\ & = & {a}_{p}^{2}{n}_{p}\left[\displaystyle \frac{{A}_{d}}{2}({k}_{p}^{2}+{h}_{p}^{2})+{B}_{d}({k}_{p}\cos {\varpi }_{d}+{h}_{p}\sin {\varpi }_{d})\right].\end{eqnarray} \tag{ 5 }$$

Note that B_d is different from the disk eccentricity vector e_d = e_d(cosϖ_d, sin ϖ_d); in fact, B_d involves a convolution of e_d(a) with a complicated kernel (see Equation (48)) over the radial extent of the disk.

The mathematical structure of R_d in Equation (3) is identical to that of the conventional Laplace–Lagrange secular disturbing function (Murray & Dermott 1999). The difference between them lies in the explicit dependence of the coefficients A_d and B_d on a_p. The behavior of A_d(a_p) and B_d(a_p) is determined, eventually, by the radial profiles of the disk surface density Σ_d and eccentricity e_d. The exploration of this behavior is the major focus of our study.

Equation (45) for A_d consists of two parts. One of them, ${A}_{d}^{\mathrm{bulk}}$, involves integral convolution over the disk semimajor axis a of a (prescribed) disk surface density ${{\rm{\Sigma }}}_{d}(a)$, as well as its radial derivatives up to second order, all of which enter linearly; see Equation (46). The convolution kernel involves a Laplace coefficient ${b}_{1/2}^{(0)}(\alpha )$, which features a weak (logarithmic) singularity arising as α = a/a_p → 1 (i.e., for disk annuli close to the particle orbit). This singularity is, however, fully integrable and vanishes upon integration over the radial extent of the disk (see Appendix A). As always, the Laplace coefficients are defined as

$$\begin{eqnarray}&&{b}_{s}^{(j)}(\alpha )={\pi }^{-1}{\int }_{0}^{2\pi }\displaystyle \frac{\cos (j\psi )d\psi }{{(1-2\alpha \cos \psi +{\alpha }^{2})}^{s}}.\end{eqnarray} \tag{ 6 }$$

In addition, whenever the disk has sharp edges (or discontinuous transitions in surface density), Equation (45) for A_d features boundary terms ${A}_{d}^{\mathrm{edge}}$, given by Equation (47). These terms involve ${b}_{1/2}^{(0)}(\alpha )$ and its derivative ${b}_{1/2}^{(0)^{\prime} }(\alpha )$, both evaluated at α = a_e/a_p, where a_e is the semimajor axis of the disk edge. As the semimajor axis of a test particle a_p approaches the disk edge, both ${b}_{1/2}^{(0)}(\alpha )$ and ${b}_{1/2}^{(0)^{\prime} }(\alpha )$ diverge: ${b}_{1/2}^{(0)}(\alpha )\sim \mathrm{ln}| 1-\alpha |$, while ${b}_{1/2}^{(0)^{\prime} }(\alpha )\sim | 1-\alpha {| }^{-1}$, where in this case $\alpha ={a}_{e}/{a}_{p}\to 1$. As a result, boundary terms, as well as A_d, diverge at the sharp disk edge. This singularity of the secular disturbing function is further explored in Section 4.

The divergence of A_d at the disk edge does not arise if ${{\rm{\Sigma }}}_{d}(a)$ goes to zero at the boundary in a sufficiently smooth fashion. Indeed, ${b}_{1/2}^{(0)^{\prime} }(\alpha )$ in the boundary terms in Equation (47) is multiplied by Σ_d at the edge; ${b}_{1/2}^{(0)}(\alpha )$ is multiplied by both Σ_d and its radial derivative ${{\rm{\Sigma }}}_{d}^{{\prime} }$. As a result, boundary terms vanish whenever ${{\rm{\Sigma }}}_{d}\propto | a-{a}_{e}{| }^{\kappa }$ and $\kappa \gt 1$ near the edge at a_e. Thus, in disks with Σ_d smoothly (faster than linearly in $| a-{a}_{e}|$) turning to zero at a finite semimajor axis a_e, the coefficient A_d has no boundary contributions and, hence, does not diverge at the boundary. The same is also true for disks without edges, which have their surface density smoothly declining to zero as a → 0 and $a\to \infty$. Both possibilities will be explored further in Section 4.

Similarly, Equation (45) for B_d involves radial integrals over the product of Σ_d(a) and the disk eccentricity e_d(a), as well as their derivatives up to second order (contribution ${B}_{d}^{\mathrm{bulk}}$ given by Equation (48)), in addition to the boundary terms ${B}_{d}^{\mathrm{edge}}$ represented by Equation (48). Once again, ${B}_{d}^{\mathrm{edge}}$ vanishes whenever Σ_d goes to zero at the disk edges sufficiently rapidly, e.g., as ${{\rm{\Sigma }}}_{d}\sim | a-{a}_{e}{| }^{\kappa }$ with κ > 1.

Other features of A_d(a_p) and B_d(a_p) behavior and their effect on secular dynamics will be discussed in Section 4–7.

One approach to calculating orbital evolution in the potential of an eccentric disk uses Lagrange equations for the evolution of orbital elements, ${\dot{k}}_{p}=-{({a}_{p}^{2}{n}_{p})}^{-1}\partial {R}_{d}/\partial {h}_{p}$, ${\dot{h}}_{p}={({a}_{p}^{2}{n}_{p})}^{-1}\partial {R}_{d}/\partial {k}_{p}$, where we use the secular expression of Equation (3) for the disturbing function R_d. As a result, one finds

$$\begin{eqnarray}&&{\dot{k}}_{p}=-{A}_{d}{h}_{p}-{B}_{d}\sin {\varpi }_{d},\,\,{\dot{h}}_{p}={A}_{d}{k}_{p}+{B}_{d}\cos {\varpi }_{d},\end{eqnarray} \tag{ 7 }$$

with the general solution given by the superposition of the free and forced eccentricity vectors (Murray & Dermott 1999):

$$\begin{eqnarray}&&{{\boldsymbol{e}}}_{p}={{\boldsymbol{e}}}_{\mathrm{free}}(t)+{{\boldsymbol{e}}}_{\mathrm{forced}},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{{\boldsymbol{e}}}_{\mathrm{free}}(t)={e}_{\mathrm{free}}(\cos ({A}_{d}t+{\varpi }_{0}),\sin ({A}_{d}t+{\varpi }_{0})),\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{{\boldsymbol{e}}}_{\mathrm{forced}}={e}_{\mathrm{forced}}(\cos {\varpi }_{d},\sin {\varpi }_{d}).\end{eqnarray} \tag{ 10 }$$

Here the forced eccentricity is

$$\begin{eqnarray}&&{e}_{\mathrm{forced}}(a)=-\displaystyle \frac{{B}_{d}(a)}{{A}_{d}(a)},\end{eqnarray} \tag{ 11 }$$

and constants e_free > 0 and ϖ₀ are such that at t = 0, Equation (8) satisfies the initial condition ${{\boldsymbol{e}}}_{p}(0)=(k(0),h(0))={e}_{p}(0)(\cos {\varpi }_{p}(0),\sin {\varpi }_{p}(0))$, where ${e}_{p}(0)$ and ${\varpi }_{p}(0)$ are the initial eccentricity and periastron angle of an object.

In particular, an object starting on a circular orbit (e_p(0) = 0) has e_free = e_forced, and its motion is described by

$$\begin{eqnarray}{e}_{p}(t) & = & \left|\displaystyle \frac{2{B}_{d}}{{A}_{d}}\sin \displaystyle \frac{{A}_{d}t}{2}\right|,\\ \tan {\varpi }_{p}(t) & = & \tan \left(\displaystyle \frac{{A}_{d}t}{2}+{\varpi }_{d}+\displaystyle \frac{\pi }{2}\right),\end{eqnarray} \tag{ 12 }$$

where ϖ_p stays in the range (0, 2π). We will use this solution for numerical verification of our semi-analytical results, although any alternative solution corresponding to different initial conditions would work just as well.

Full determination of e_p requires calculating A_d and B_d for the given profiles of the disk surface density Σ_d(a) and eccentricity ϖ_d(a). We do this with the help of Equations (45)–(49) by numerically evaluating the corresponding integrals. Despite the weak (logarithmic) singularity of the integrands, the integrals themselves are convergent and calculated directly without resorting to any type of softening of the integrand near its singular point (which is often done in studies of disk dynamics; see, e.g., Tremaine 2001; Hahn 2003; Touma et al. 2009). The lack of an ad hoc parameter (softening length) in our calculation is its important distinctive feature.

Once A_d and B_d are calculated, Equations (12) and (13) (or, more generally, Equations (8)–(11)) provide a full semi-analytical description of the particle motion in the disk potential.

An alternative method to compute particle motion uses the direct orbit integrator MERCURY (Chambers 1999), which employs the Bulirsch–Stoer algorithm (Press et al. 1992). This integration uses as its input the gravitational acceleration g due to the disk potential, which is calculated on a grid of (245, 100) points⁷ in (r, φ). Values of g on the grid are then interpolated to provide accurate accelerations everywhere in the disk.

At each grid point, g is computed by direct summation of ∇ Φ_d (see Equation (21)) over the disk surface, with the surface density given by Equation (2) in its exact form, i.e., not performing small-e_d expansion. This two-dimensional integral is convergent in the Cauchy principal value sense, even though the integrand diverges at the point where acceleration is calculated. To avoid this mathematical singularity, our numerical evaluation is performed at a small height (10⁻³ au or smaller) above the disk; we make sure that the result is convergent with respect to the value of this height. Note that this procedure is not equivalent to the introduction of softening in our theoretical calculation.

We integrate particle orbits starting with e_p = 0. This allows us to directly compare the results with Equations (12) and (13) following from our semi-analytical calculation, thus directly verifying the accuracy of our framework for the secular evolution in the potential of an eccentric disk.

In our comparison effort, we use the following model distributions of the disk eccentricity e_d(a):

$$\begin{eqnarray}&&{e}_{p}(r)={e}_{0}\left(1+\displaystyle \frac{{a}_{2}-a}{{a}_{2}-{a}_{1}}\right),\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{e}_{p}(r)={e}_{0}\displaystyle \frac{{(a-{a}_{1})}^{2}{(a-{a}_{2})}^{2}}{{({a}_{2}-{a}_{1})}^{4}},\end{eqnarray} \tag{ 14 }$$

and surface density at periastron Σ_d(a),

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{d}(a)={{\rm{\Sigma }}}_{0}\displaystyle \frac{{(a-{a}_{1})}^{2}{(a-{a}_{2})}^{2}}{{({a}_{2}-{a}_{1})}^{4}},\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{d}(a)={{\rm{\Sigma }}}_{0}\exp \left[\displaystyle \frac{4-{[(a/{a}_{c})+({a}_{c}/a)]}^{2}}{0.18}\right],\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{d}(a)={{\rm{\Sigma }}}_{0}\left[1+4\sin \left(\pi \displaystyle \frac{a-{a}_{1}}{{a}_{2}-{a}_{1}}\right)\right].\end{eqnarray} \tag{ 17 }$$

with a₁ = 0.1 au, a₂ = 5 au, and a_c = 1.5 au and e₀ and Σ₀ being the normalization factors that we vary in our models. In our calculations, we always assume the mass of the central star to be 1 M_⊙.

Figure 2 illustrates the behavior of these profiles of Σ_d and e_d. In the majority of our calculations, we use the linear eccentricity profile given by Equation (13), corresponding to a disk that is everywhere eccentric. Disks with the eccentricity profile given by Equation (14) studied in Section 4.3.2 have circular inner and outer edges but are eccentric in between.

Zoom In Zoom Out Reset image size

The surface density profile given by Equation (15) corresponds to a disk in which Σ_d smoothly goes to zero at the edges at a₁ and a₂. Based on the discussion in Section 2.2, we expect A_d and B_d not to diverge at the edges of such a disk, which is verified in Section 4.1.1. Equation (16) describes a disk without edges (Section 4.1.2), which has its surface density exponentially decreasing for both a ≪ a_c and a ≫ a_c. This disk also does not feature boundary terms (as it has no boundaries). Finally, the Σ_d profile given by Equation (17) describes a disk with discontinuous drops of the surface density at the edges. In this disk, we expect A_d and B_d to diverge at the boundaries, which is verified in Section 4.1.3.

Figure 3 illustrates the 2D distributions of the surface density obtained via Equation (2) with some of the Σ_d and e_d profiles listed above. One can see that Σ(r, φ) can have a rather complicated structure, depending on the particular disk model used. The contours of constant Σ(r, φ) (left) are very different from the elliptical trajectories (right) of the mass elements giving rise to the surface density distribution in the disk.

Zoom In Zoom Out Reset image size

We now test the accuracy of our secular theory in disks with different Σ_d(a) profiles. This allows us to examine the different possible behaviors of the secular coefficients A_d and B_d, as well as to see how sensitive the agreement with the numerical results is to the different features of the Σ_d(a) distributions.

4.1.1. Σ_d Smoothly Vanishing at the Edges

We start by looking at the secular effect of a disk, in which Σ_d smoothly goes to zero at the boundaries, i.e., the one with Σ_d given by Equation (15). In Figure 5, we display the theoretical radial profiles of A_d, B_d, and 2B_d/A_d for such a disk. The curves of A_d(a_p) are compared with the numerically determined ${\dot{\varpi }}_{\sec }$ (blue dots), while the theoretical $2{B}_{d}({a}_{p})/{A}_{d}({a}_{p})$ is compared against the numerical ${e}_{p}^{{\rm{m}}}$. This particular calculation uses the linear eccentricity profile e_d(a) given by Equation (13) with relatively low e₀ = 0.01, so that the maximum value of the disk eccentricity (reached at its inner edge) is 0.02. Thus, the requirement e_d ≪ 1 necessary for the secular results of Equations (3) and (45)–(49) to apply is fulfilled.

Zoom In Zoom Out Reset image size

One can see an almost perfect match between the numerical results and theory in terms of both the amplitude of the disk-induced particle eccentricity and the period of the associated orbital precession. Theoretical calculation easily reproduces even the very subtle features of the ${e}_{p}^{{\rm{m}}}(a)$ behavior, including the variations happening on very short radial scales manifesting themselves as sharp features in Figure 5. There are several features to note in this figure.

First, A_d has a different sign in different parts of the disk: precession of the free eccentricity vector is prograde near the boundaries of the disk, while it is retrograde away from the edges. This change of sign of ${\dot{\varpi }}_{\sec }$ was clear already in Figure 4 (which is drawn for the same disk model as Figure 5, at the locations highlighted with circles in the latter), but Figure 5 provides a much more detailed representation of the different characteristics of the secular effect of the disk.

Second, B_d also changes sign as a_p varies. As a result of sign variations of both A_d and B_d, the forced eccentricity vector ${{\boldsymbol{e}}}_{\mathrm{forced}}$ can be aligned with the apsidal line of the disk in some intervals of a_p and anti-aligned with it in others.

Third, both numerical and analytical ${e}_{p}^{{\rm{m}}}$ exhibit formal singularity at two distinct locations in this disk, a ≈ 1.54 and ≈4.39 au. The origin of these singularities can be traced to Equation (11) for the forced eccentricity (recall that our theory predicts ${e}_{p}^{{\rm{m}}}=2{e}_{\mathrm{forced}}$), which has A_d in the denominator, and the fact that A_d crosses zero (i.e., the direction of e_free precession changes) at these locations; see Figure 5(b). We will discuss these singularities in more detail in Section 4.2.

Fourth, at certain values of the semimajor axis (at ≈0.92 and ≈4.39 au), the disk-induced forced eccentricity vanishes. This happens because B_d changes sign at these locations so that B_d → 0; see Figure 5(c) and Equation (11). Interestingly, one of the semimajor axes where ${e}_{p}^{{\rm{m}}}\to 0$ (${a}_{p}\approx 4.39$ au) is located in the immediate proximity of the singularity of ${e}_{p}^{{\rm{m}}}$. This occurs because in this part of the disk, A_d and B_d go through zero at almost the same (but still slightly different) values of a_p.

Fifth, for the surface density profile given by Equation (15), we find that A_d and B_d remain finite everywhere in the disk, including the boundaries. This is in line with the expectations outlined in Section 2.2 for the disks with Σ_d smoothly going to zero at the edges, which do not result in divergent boundary terms in A_d and B_d.

4.1.2. Disk without Boundaries

Next we consider secular dynamics in the potential of a Gaussian ring with Σ_d given by Equation (16), which does not feature well-defined edges. We plot the behavior of the corresponding ${e}_{p}^{{\rm{m}}}$, A_d, and B_d in Figure 6.

Zoom In Zoom Out Reset image size

One can see that many features of secular dynamics present in the case studied in Section 4.1.1 are present here as well: both A_d and B_d are finite, vary with a_p, and change sign, while ${e}_{p}^{{\rm{m}}}$ exhibits both singularities and nulls. This is likely related to the fact that these two Σ_d profiles are morphologically similar: they exhibit no discontinuities, are single-peaked, and smoothly decay to zero away from the peak. The only obvious difference is that in the Gaussian case in Figure 6, the nulls of ${e}_{p}^{{\rm{m}}}$ are located very close to both singularities of ${e}_{p}^{{\rm{m}}}$. However, the significance of this difference is unclear.

4.1.3. Σ_d Sharply Truncated at the Edges

We also examined secular behavior for the disk model in which surface density displays a discontinuous drop to zero at the inner and outer edges, namely, the one represented by Equation (17). Sharply truncated Σ distributions are rather typical in planetary rings and other astrophysical disks, and it is important to study the effect of their gravity on the dynamics of embedded objects.

Figure 7 illustrates secular dynamics in the potential of such a disk. It again shows the variation of sign of A_d, which results in the emergence of two singularities of ${e}_{p}^{{\rm{m}}}$—one very close to the inner edge of the disk at a_p = 0.11 au and another at ${a}_{p}\approx 1.28\,\mathrm{au}$. However, for this disk model, B_d does not change sign—it always stays positive. As a result, there are no nulls of ${{\boldsymbol{e}}}_{p}^{{\rm{m}}}$ for orbits in the potential of such a disk.

Zoom In Zoom Out Reset image size

An even more dramatic difference with the cases explored in Sections 4.1.1–4.1.2 is the behavior of A_d and B_d near the disk edges. Unlike the two previous cases, in which both coefficients remained finite everywhere at all radii, in a sharply truncated disk, A_d and B_d exhibit singularity as the disk edge is approached. This is very clearly seen at the outer⁸ boundary of the disk in Figure 7, where both theoretical curves and the results of orbit integrations exhibit divergent behavior. Unfortunately, we cannot probe this divergence in great detail numerically, as the orbits of test particles start crossing the edge of the disk.

At the same time, the test particle eccentricity ${e}_{p}^{{\rm{m}}}$ remains finite at the disk edges, even though both A_d and B_d are singular there. This is because both coefficients of the disturbing function given by Equation (3) diverge in similar fashion at the edge, so that ${e}_{\mathrm{forced}}$ remains finite there.

A common feature for all ${{\rm{\Sigma }}}_{d}$ distributions examined in Section 4.1 is the emergence of multiple singularities of ${e}_{p}^{{\rm{m}}}$. At these locations, ${e}_{p}^{{\rm{m}}}$ formally diverges, and the assumption ${e}_{p}\,\ll 1$ used in deriving our secular disturbing function (Equations (3)–(5)) breaks down. The way in which this happens is illustrated in Figure 8, where we compare the numerical and analytical results in the vicinity of one of the singularities (near a_p = 1.54 au) in a disk with ${{\rm{\Sigma }}}_{d}(a)$ given by Equation (15); see Figure 5.

Zoom In Zoom Out Reset image size

One can see that as the theoretical singularity is approached, the behavior of ${\dot{\varpi }}_{\sec }$ starts to deviate from the prediction of Equations (45)–(47). As a result, ${\dot{\varpi }}_{\sec }$ goes through zero at a location slightly different from the one where A_d = 0. Note that even at the point where ${\dot{\varpi }}_{\sec }=0$, particle eccentricity remains finite (even though it reaches values close to 1). This means that Equation (11) is no longer valid when ${e}_{p}\sim 1$ and that additional, higher-order terms become important in addition to the lowest-order secular potential contribution (Equation (3)). This discrepancy could be at least partly ameliorated by including higher-order (in e_p) terms in the calculation of R_d, as was done recently in Sefilian & Touma (2018) for the case of power-law disks.

As evidenced by Figures 5 and 6, ${e}_{p}^{{\rm{m}}}$ singularities often occur in the immediate vicinity of nulls of B_d. This results in a characteristic shape of these singularities, with ${e}_{p}^{{\rm{m}}}$ sharply dropping to zero in close proximity to the singularity. This leads to a dramatic difference in the eccentricities of particles with almost identical semimajor axes, resulting in their orbits crossing. Such locations thus provide a natural environment for particle collisions.

It is not clear why in some cases conditions A_d = 0 and B_d = 0 get realized at almost the same value of a. Inspection of the integrands in Equations (46) and (48) does not reveal an obvious reason for that to be the case. Interestingly, we find such "null-singularity" pairs only in the two disks without sharp edges (Section 4.1.1–4.1.2), for which the boundary terms ${A}_{d}^{\mathrm{edge}}$ and ${B}_{d}^{\mathrm{edge}}$ in Equation (45) vanish. The disk with sharply truncated ${{\rm{\Sigma }}}_{d}$ (see Section 4.1.3 and Figure 7) has singularities without neighboring nulls of B_d (in fact, B_d does not change sign in this disk). Whether this outcome is due to the nontrivial boundary terms in this disk model is not clear. These curious properties of ${e}_{p}^{{\rm{m}}}$ singularities deserve further investigation.

Next we examine how our secular theory fares against changes of the disk eccentricity. First, we explore the effect of uniformly varying just the amplitude of eccentricity (Section 4.3.1), keeping the radial profile of e_d(a) the same. We then look at the effect of a different radial profile of e_d(a) on the agreement between the theory and numerical calculations (Section 4.3.2).

4.3.1. Variation of the Disk Eccentricity Amplitude

In Figure 9, we plot ${e}_{p}^{{\rm{m}}}(a)$ for the disk model with ${{\rm{\Sigma }}}_{d}$ and e_d given by Equations (15) and (13), where we set ${{\rm{\Sigma }}}_{0}\,\approx 1100$ g cm⁻² but vary the eccentricity normalization e₀ as indicated in the figure (which is very similar to Figure 5(a)).

Zoom In Zoom Out Reset image size

One can see that our theory works surprisingly well in predicting ${e}_{p}^{{\rm{m}}}(a)$, even when particle eccentricity reaches values in excess of e_p = 0.5, when one would naively expect the description based on the lowest-order secular disturbing function given by Equation (3) to break down. For the e_d profile given by Equation (13), the maximum value of e_d, reached at the inner boundary of the disk, is 2e₀. Thus, even for disks with inner eccentricities reaching ${e}_{d}({a}_{1})=0.4$, our theory performs quite well.

The amplitude of eccentricity oscillations ${e}_{p}^{{\rm{m}}}(a)$ is just one metric by which the performance of our theory can be judged. Another obvious one is the free precession rate ${\dot{\varpi }}_{\sec }$. In Figure 10, we plot the ratio of ${\dot{\varpi }}_{\sec }$ to its analytical counterpart A_d as a function of the disk eccentricity normalization e₀. In the framework of our secular calculation, the precession rate should not depend on e₀ and must be equal to A_d.

Zoom In Zoom Out Reset image size

Figure 10 shows that this is not really the case and ${\dot{\varpi }}_{\sec }$ does deviate from A_d when e₀ is nonnegligible. The agreement between these two frequencies is somewhat less impressive than for ${e}_{p}^{{\rm{m}}}(a)$, with ${\dot{\varpi }}_{\sec }$ deviating from A_d by tens of percent already for e₀ = 0.1. Note that the discrepancy between the numerical and analytical secular frequencies is a strong function of the semimajor axis a_p, with the largest deviation occurring in the vicinity of the ${\dot{\varpi }}_{\sec }$ singularity (just interior of it) at 1.3 au. Thus, one should be careful when applying our lowest-order secular calculation to characterize the orbital evolution of test particles at certain locations in the disk with ${e}_{d}\gtrsim 0.1$.

4.3.2. Variation of the Disk Eccentricity Profile

Next we study the effect of changing the radial profile of the disk eccentricity e_d(a). Figure 11 is the analog of Figure 9 but made for a disk with an eccentricity profile given by Equation (14).

Zoom In Zoom Out Reset image size

The first thing to note is that the radial profile of the theoretical ${e}_{p}^{{\rm{m}}}(a)$ still shows two clear singularities at the same radial locations as in Figure 9. This is easy to understand, since ${e}_{p}^{{\rm{m}}}$ diverges at radii where ${A}_{d}\to 0$, and A_d is independent of the disk eccentricity profile (it depends only on the ${{\rm{\Sigma }}}_{d}(a)$ profile). As a result, singularities of ${e}_{p}^{{\rm{m}}}$ stay at fixed locations even though e_d(a) varies.

Second, the detailed shape of ${e}_{p}^{{\rm{m}}}(a)$ is notably different from that shown in Figures 5 and 9, despite the fact that the ${{\rm{\Sigma }}}_{d}(a)$ profile is the same in both cases. While previously, ${e}_{p}^{{\rm{m}}}$ was dropping to zero right next to the outer singularity at ${a}_{p}\approx 4.39\,\mathrm{au}$, in Figure 11, this happens near the inner singularity at $\approx 1.5\,\mathrm{au}$. Another null of ${e}_{p}^{{\rm{m}}}$, more distant from the singularity, has also swapped its location and now lies in the outer part of the disk in Figure 11.

Third, the agreement between the theory and direct orbit integrations is somewhat worse for the disk with an e_d(a) profile (Equation (14)). This is especially noticeable to the right of the inner singularity. For example, the ${e}_{p}^{{\rm{m}}}$ differs by $\approx 50 \%$ from the theoretical expectation at ${a}_{p}=2\,\mathrm{au}$ for e₀ = 1.75 (black curve, which corresponds to the maximum eccentricity in the disk of $\approx 0.1$). Away from the region $1.5\mbox{--}2\,\mathrm{au}$, the agreement is generally better, even though the particle eccentricity excited by the disk potential can be quite high, ${e}_{p}\gtrsim 0.1\mbox{--}0.2$.

We hypothesize that the reduced accuracy with which our secular theory predicts secular dynamics in the case of a disk with an e_d(a) profile given by Equation (14) is caused by the fact that this disk features a rather nonaxisymmetric surface density distribution, ${\rm{\Sigma }}(r,\varphi )$. Indeed, Figures 3(c) and (d) show that the disk with e_d(a) given by Equation (14) features a well-defined concentration of mass in the outer region near the apastra of the constituent trajectories. As a result, for certain values of a, surface density ${\rm{\Sigma }}({r}_{d},{\varphi }_{d})$ exhibits substantial variations along the eccentric trajectories of the mass elements comprising the disk; see Figure 3(d). According to Equation (2), this is only possible if $\zeta {e}_{d}$ is not small for these values of a, meaning that the assumption $| \zeta {e}_{d}| \ll 1$ underlying our expansion given by Equation (22) is not fulfilled, which explains the deviations seen for this particular disk model.

By contrast, the disk with the eccentricity profile given by Equation (13) shows a more axisymmetric distribution of ${\rm{\Sigma }}(r,\varphi )$ (see Figures 3(a) and (b)), meaning smaller $| \zeta {e}_{d}|$ and higher accuracy of the expansion given by Equation (22). Thus, we expect that our secular theory should perform better for disks without highly nonuniform azimuthal features in ${\rm{\Sigma }}({r}_{d},{\varphi }_{d})$.

Our derivation of the secular disturbing function in Appendix A explicitly assumes that the particle orbits fully within the disk. In other words, in the case of a disk with ${{\rm{\Sigma }}}_{d}$ dropping to zero at some finite a₁ and a₂, the particle semimajor axis a_p satisfies ${a}_{1}\lt {a}_{p}\lt {a}_{2}$, and its eccentric orbit does not cross the boundaries of the disk. Obviously, if the disk has no edges, as in, e.g., the case of a Gaussian ring given by Equation (16), then the particle orbit is always fully enclosed within the disk.

However, close examination of our derivation of R_d in Appendix A demonstrates that it should also apply equally well outside the disk with sharp edges, as long as the particle orbit does not cross the disk boundaries. For example, in the case of a particle orbiting outside the outer edge of the disk (${a}_{p}\gt {a}_{2}$), first, one drops the contribution of the outer disk (i.e., $a\gt {a}_{p}$, as there is no disk material there) when computing R_d, and, second, the integration in the inner disk runs not up to $\alpha =1$ but only up to $\alpha ={a}_{2}/{a}_{p}\lt 1$. As a result, the resultant Equations (45)–(49) for the coefficients A_d and B_d apply without modification, even when ${a}_{p}\lt {a}_{1}\lt {a}_{2}$ or ${a}_{1}\lt {a}_{2}\lt {a}_{p}$.

To verify this claim, in Figures 12 and 13, we show radial profiles of ${e}_{p}^{{\rm{m}}}$ and ${\dot{\varpi }}_{\sec }$, as well as those of $2{e}_{\mathrm{forced}}$, A_d, and B_d for particles orbiting inside and outside (correspondingly) the radial extent of the disk. Disk parameters are the same as in Figure 7; in particular, ${{\rm{\Sigma }}}_{d}$ is given by Equation (17) with sharp edges at ${a}_{1}=0.1\,\mathrm{au}$ and ${a}_{2}=5\,\mathrm{au}$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

One can see that, as in the previous cases illustrated in Figures 5–7, there is excellent agreement between our theory and direct orbit integrations, as long as the disk and particle eccentricities are low. Both the theoretical and the numerical values of the particle eccentricity extrapolated to the disk edges match the corresponding values extrapolated from inside the disk; see Figure 7. At the same time, both A_d and B_d diverge as ${a}_{p}\to {a}_{1}-0$ and ${a}_{p}\to {a}_{2}+0$, mirroring the singularity of these coefficients identified previously (Section 2.2).

These results extend the applicability of our calculation of R_d to (coplanar) particles having arbitrary semimajor axes relative to the disk, as long as their orbits do not cross the edge of the disk where ${{\rm{\Sigma }}}_{d}$ discontinuously drops to zero. However, it can be shown that even this latter constraint can be removed, extending the applicability of our results even further. We do not dwell on this point here,⁹ deferring it to a future study.

Some astrophysical disks are known to have very sharp edges. For example, using Voyager 2 occultation data, Graps et al. (1995) demonstrated that the ring of Uranus has an ${{\rm{\Sigma }}}_{d}$ steeply dropping to zero at the ring boundaries. The edges of the Saturn rings are also known to be very sharp (Tiscareno 2013). Outside the solar system, eclipse data reveal the sharpness of the edge of the circumbinary ring around the young star KH 15D (Winn et al. 2006).

Our calculations predict that sharp edges result in the divergence of the secular potential of a razor-thin disk, leading to a divergence in the precession rate near these locations (Section 4.1.3). This outcome, resulting from the nonvanishing boundary terms, was previously pointed out in Silsbee & Rafikov (2015a) for truncated power-law disks, and now we generalize it for other models of eccentric disks with edges. This prediction is nicely confirmed by the direct integration of particle orbits in a particular disk model with sharp edges; see Figure 7 and Section 4.1.3.

The divergence of ${\dot{\varpi }}_{\sec }$ (or A_d) near the sharp edge of a zero-thickness disk can be traced to the fact that the (in-plane) gravitational acceleration in this region behaves as ${g}_{d}\propto \mathrm{ln}| {\rm{\Delta }}r|$, where ${\rm{\Delta }}r$ is the separation from the edge. Specializing to the case of an axisymmetric disk, one then finds the free precession rate (Fontana & Marzari 2016)

$$\begin{eqnarray}&&{\dot{\varpi }}_{\sec }=-\displaystyle \frac{n}{2{{rg}}_{{\rm{K}}}}\displaystyle \frac{{{dg}}_{d}}{{dr}}\propto {({\rm{\Delta }}r)}^{-1}\end{eqnarray} \tag{ 18 }$$

(it should be ${g}_{{\rm{K}}}={{GM}}_{\star }/{r}^{2}$) near the edge at the leading order. The divergent behavior ${\dot{\varpi }}_{\sec }\propto {({\rm{\Delta }}r)}^{-1}$ coincides with the scaling of the boundary terms in the expression (Equation (47)) for ${A}_{d}^{\mathrm{edge}}$; see Section 2.2. Analogous singularities should arise at any radius in the disk where ${{\rm{\Sigma }}}_{d}(a)$ exhibits a discontinuity.

A disk with ${{\rm{\Sigma }}}_{d}$ dropping to zero smoothly over a narrow but finite range of a would not have A_d diverging there, as the boundary terms (${A}_{d}^{\mathrm{edge}}$ and ${B}_{d}^{\mathrm{edge}}$) vanish for smooth ${{\rm{\Sigma }}}_{d}$ profiles. Nevertheless, A_d still exhibits a nontrivial behavior in this region. This is illustrated in Figure 14, where we plot ${\dot{\varpi }}_{\sec }={A}_{d}(a)$ (top panel) for rings with several ${{\rm{\Sigma }}}_{d}(a)$ profiles¹⁰ of different degree of steepness near the boundaries (bottom panel). One can see that near the edge, ${\dot{\varpi }}_{\sec }$ exhibits very rapid variation, changing from large negative values just inside the edge to large positive values just outside the edge.

Zoom In Zoom Out Reset image size

This behavior can be understood by noticing that Equation (46) for ${A}_{d}^{\mathrm{bulk}}={A}_{d}$ contains a second derivative of the surface density ${{\rm{\Sigma }}}_{d}$ inside the integral. Very close to the boundary, where ${{\rm{\Sigma }}}_{d}$ is still high, ${{\rm{\Sigma }}}_{d}^{{\prime\prime} }(a)$ is large and negative. Due to the logarithmic singularity of ${b}_{1/2}^{(0)}(\alpha )$ as $\alpha \to 1$, the radial convolution in Equation (46) enhances the contribution of this region (with large ${{\rm{\Sigma }}}_{d}^{{\prime\prime} }\lt 0$) to A_d, resulting in rapid retrograde precession in this part of the disk. On the other hand, just slightly away from the boundary, where ${{\rm{\Sigma }}}_{d}$ is low, ${{\rm{\Sigma }}}_{d}^{{\prime\prime} }(a)$ is large and positive, driving fast prograde precession there.

As the sharpness of the ${{\rm{\Sigma }}}_{d}$ profile near the boundaries increases, so does the magnitude of ${{\rm{\Sigma }}}_{d}^{{\prime\prime} }(a)$. As a result, the amplitude of ${\dot{\varpi }}_{\sec }={A}_{d}$ on both sides of the edge grows. In the limit of the infinitely sharp ${{\rm{\Sigma }}}_{d}$ transition at the edge, the behavior of A_d becomes singular, with the change of sign at the edge, in agreement with the expectation of Equation (18). In our formalism, this is accounted for mathematically by the appearance of the boundary terms of Equation (47), whereas ${{\rm{\Sigma }}}_{d}^{{\prime\prime} }$ in the integral of Equation (46) remains finite.

In a disk with small but finite vertical thickness $h\ll a$, the behavior of the coefficients of R_d would be slightly different. In such a disk, the rise of A_d and B_d would saturate at a finite value ∝h⁻¹ as the edge is approached. This transition is easy to understand by setting ${\rm{\Delta }}r\sim h$ in Equation (18).

The divergence of A_d near the sharp edges can have important implications for, e.g., the dynamics of planetary rings. Even though ${{\rm{\Sigma }}}_{d}$ remains finite at the edge, our results demonstrate that particle orbits should precess very rapidly and at a rapidly changing rate as their semimajor axes get closer to the edge. This naturally leads to particle orbit crossing resulting in their collisions, helping redistribute angular momentum near the disk edge (Chiang & Goldreich 2000; Chiang & Culter 2003).

Disk gaps present another example of sharp surface density gradients. Gaps could form naturally as a result of gas clearing by the gravitational torques due to massive planets orbiting within the disk.

Ward (1981) looked at the effect of gaps on the free precession rate of test particles in axisymmetric disks with a power-law profile of ${{\rm{\Sigma }}}_{d}$, finding ${\dot{\varpi }}_{\sec }$ to be negative within the disk but positive within the gap. Ward (1981) modeled the gap by setting ${{\rm{\Sigma }}}_{d}$ to zero within a range of semimajor axes, which makes his result quite natural. Indeed, a disk with such a gap can be viewed as a combination of two disjoint disks with sharp edges. An object orbiting within a gap is thus moving exterior to an inner truncated disk and interior to an outer truncated disk. Our results in Sections 4.1.3, 4.4, and 5.1 demonstrate that interaction with each of the disks drives prograde free precession of such an object, with the combined effect being simply a linear superposition (also prograde) of the two contributions.

Using our theoretical formalism, we can explore how the results of Ward (1981) change for more realistic (i.e., less sharp) gap profiles. We consider a disk, which is a combination of a wide ring with flat ${{\rm{\Sigma }}}_{d}$ given by Equation (50), and a gap of width w and relative depth d (changing from d = 0 for no gap to d = 1 for ${{\rm{\Sigma }}}_{d}=0$ in the gap center), with the profile given by Equation (51). In Figure 15, we show how ${A}_{d}={\dot{\varpi }}_{\sec }$ varies as we change the gap depth. One can see that A_d inside the gap region, which is negative in the absence of a ${{\rm{\Sigma }}}_{d}$ depression, gradually decreases in magnitude, crosses zero, and becomes large and positive as the gap depth is increased. At the edges of the gap, A_d exhibits a nontrivial structure reminiscent of that seen in Figure 14. If our gap had a more abrupt drop of ${{\rm{\Sigma }}}_{d}$ at its edges, we would have converged to the case explored by Ward (1981).

Zoom In Zoom Out Reset image size

Figure 16 looks at the effect of variation of the gap width, keeping its depth constant. One can see that narrower gaps have a more pronounced effect on free precession in the center of the gap. This is because the edges of such gaps are sharper and also closer to the orbit of a precessing object. In this case, the intuition developed in Figure 14 again suggests that A_d should be large and positive, as is observed in Figure 16. To summarize, the effect of a gap on free precession seems to be determined more by the sharpness of the ${{\rm{\Sigma }}}_{d}$ gradient than by the gap width or depth separately.

Zoom In Zoom Out Reset image size

The nontrivial behavior of ${\dot{\varpi }}_{\sec }$ discussed in Sections 5.1–5.2 is caused by the localized sharp features in ${{\rm{\Sigma }}}_{d}(a)$. However, our results in Sections 4.1.1 and 4.1.3, confirmed by direct orbit integration, clearly demonstrate that A_d can also easily change sign within a disk with a smooth distribution of ${{\rm{\Sigma }}}_{d}$. Since, according to Equation (9), ${A}_{d}={\dot{\varpi }}_{\sec }$ is the rate of precession of the free eccentricity vector, this implies that free precession of a test particle can be both prograde and retrograde, depending on its location within the disk. The possibility of a change of sign of free precession may appear rather surprising, given that the simple truncated power-law models usually employed to study the gravitational effect of thin disks tend to predict ${\dot{\varpi }}_{\sec }\lt 0$, i.e., retrograde free precession (Heppenheimer 1980; Ward 1981; Batygin et al. 2011). Prograde precession is (locally) possible in disks featuring gaps, i.e., sharp drops of ${{\rm{\Sigma }}}_{d}$ (Ward 1981; see also Section 5.2), but the disks explored in Sections 4.1.1 and 4.1.2 have rather smooth profiles of ${{\rm{\Sigma }}}_{d}(a)$.

On the other hand, in Silsbee & Rafikov (2015a), it was shown that even the truncated power-law disks can exhibit prograde free precession far from the disk edges for certain values of the density slope $p=-d\mathrm{ln}{{\rm{\Sigma }}}_{d}/d\mathrm{ln}a$, e.g., when $p\lt 0$, so that ${{\rm{\Sigma }}}_{d}$ grows with a. Given that near the edges (but still within the disk), one expects ${A}_{d}\lt 0$ (see Section 5.1), the direction of free precession must change sign at some location inside such a disk. In light of this observation, our current results simply show that the change of sign of free precession is a rather common phenomenon for arbitrary profiles of ${{\rm{\Sigma }}}_{d}$.

Unfortunately, predicting the direction of the free precession (i.e., the sign of A_d) at a specific semimajor axis a in a disk with a given profile ${{\rm{\Sigma }}}_{d}(a)$ is not an easy task. Even in a disk without sharp boundaries, when the boundary term ${A}_{d}^{\mathrm{edge}}$ vanishes and ${A}_{d}={A}_{d}^{\mathrm{bulk}}$ is represented by Equation (46), it is generally not straightforward to predict a priori the sign of this integral term. Indeed, a smooth, continuous ${{\rm{\Sigma }}}_{d}(a)$ gradually decaying to zero at finite (e.g., given by Equation (15)) or infinite (e.g., Equation (16)) boundaries would necessarily have ${{\rm{\Sigma }}}_{d}^{{\prime\prime} }(a)$ changing sign within the disk. Integration over a provides a nontrivial, nonlocal mapping between the global behavior of ${{\rm{\Sigma }}}_{d}^{{\prime\prime} }(a)$ and the value (and sign) of A_d.

As discussed in Sections 4.1 and 4.2, a change of sign of A_d inside the disk is also important because it gives rise to very high (formally divergent) values of the test particle eccentricity at radii where A_d = 0, as long as the disk eccentricity e_d is nonzero. Previously, a similar effect—a localized singularity of e_p—was identified in studies of planet formation within stellar binaries, both analytically (Rafikov 2013a, 2013b; Rafikov & Silsbee 2015a, 2015b; Silsbee & Rafikov 2015b, 2015a) and numerically (Meschiari 2014). Its origin could be traced to a secular resonance caused by the cancellation of prograde precession due to the binary companion and retrograde precession driven by the disk gravity in the presence of the nonzero binary torque (and disk torque, if the disk is eccentric). The importance of these singularities for planet formation in binaries lies in the fact that high values of e_p lead to very energetic collisions between planetesimals, resulting in their destruction and hampering planetary accretion. In this work, we show that the same mechanism works even without a binary companion—disk torque is always present when ${e}_{d}\ne 0$ and results in divergent e_p whenever ${A}_{d}\to 0$, which naturally happens in our disks.

As shown in Rafikov & Silsbee (2015b) and Silsbee & Rafikov (2015b), an opposite effect is also possible in disks in binaries—at certain locations, the eccentricities of test particles affected by the combined potential of a companion and an eccentric disk become very small, as a result of the cancellation of the corresponding torques.¹¹ Again, in our case, this happens even without a binary companion—at the locations where B_d = 0, one naturally finds ${e}_{p}\to 0$, as a result of the cancellation of torques arising from different parts of the same disk. This can be seen, e.g., at $\approx 0.9\,\mathrm{au}$ in Figures 5 and 9 and at $\approx 4\,\mathrm{au}$ in Figure 11. At these locations, the relative velocities of colliding objects naturally become very small, promoting their agglomeration (rather than fragmentation) and growth.

We start by computing the integral over the circular annulus ${{\mathbb{D}}}_{c}$. When using Equation (22), it is natural to divide $R({{\mathbb{D}}}_{c})$ into two parts: one containing $\cos {\varphi }_{d}$ (we call it I₁) and one independent of ${\varphi }_{d}$ (I₂). We evaluate the ${\varphi }_{d}$-independent contribution first:

$$\begin{eqnarray}{I}_{1} & = & \langle G{\displaystyle \int }_{{r}_{d,\mathrm{in}}}^{{r}_{d,\mathrm{out}}}{{dr}}_{d}\,{{\rm{\Sigma }}}_{d}({r}_{d}){r}_{d}{\displaystyle \int }_{0}^{2\pi }\displaystyle \frac{d\theta }{\sqrt{{r}_{p}^{2}+{r}_{d}^{2}-2{r}_{p}{r}_{d}\cos \theta }}\rangle .\end{eqnarray} \tag{ 24 }$$

This integral has a singularity at r_p = r_d, so we break it into two parts, ${r}_{p}\lt {r}_{d}$ and ${r}_{p}\gt {r}_{d}$, while also rewriting the inner integral as a Laplace coefficient ${b}_{1/2}^{(0)}$; see Equation (6):

$$\begin{eqnarray}{I}_{1} & = & \pi G\left\langle {\displaystyle \int }_{{r}_{d,\mathrm{in}}}^{{r}_{p}}{{\rm{\Sigma }}}_{d}({r}_{d})\displaystyle \frac{{r}_{d}}{{r}_{p}}{b}_{1/2}^{(0)}\left(\displaystyle \frac{{r}_{d}}{{r}_{p}}\right){{dr}}_{d}\right.\left.\ +\ {\displaystyle \int }_{{r}_{p}}^{{r}_{d,\mathrm{out}}}{{\rm{\Sigma }}}_{d}({r}_{d}){b}_{1/2}^{(0)}\left(\displaystyle \frac{{r}_{p}}{{r}_{d}}\right){{dr}}_{d}\right\rangle .\end{eqnarray} \tag{ 25 }$$

Now let us make a change of variables in both integrals so as to get rid of the time dependence in the Laplace coefficients, entering via r_p. In the first integral, we set $\alpha ={r}_{d}/{r}_{p}\lt 1$, while in the second, we use $\alpha ={r}_{p}/{r}_{d}\lt 1$. We will also approximate ${r}_{d,\mathrm{in}}\approx {a}_{1}$ and ${r}_{d,\mathrm{out}}\approx {a}_{2}$, which, as we will see, introduces error of negligible order. Now we have

$$\begin{eqnarray}{I}_{1} & = & {I}_{1}^{(0)}+{I}_{1}^{(1)}=\pi G\left\langle {\displaystyle \int }_{{a}_{1}/{r}_{p}}^{1}{r}_{p}{{\rm{\Sigma }}}_{d}(\alpha {r}_{p})\alpha {b}_{1/2}^{(0)}(\alpha )d\alpha \right.\left.\ +\ {\displaystyle \int }_{{r}_{p}/{a}_{2}}^{1}{r}_{p}{{\rm{\Sigma }}}_{d}({r}_{p}/\alpha ){\alpha }^{-2}{b}_{1/2}^{(0)}(\alpha )d\alpha \right\rangle .\end{eqnarray} \tag{ 26 }$$

Our goal is to get rid of triangular brackets, that is, to time average the integrals. The time dependence is hidden in r_p and has the form ${r}_{p}(t)={a}_{p}(1-{e}_{p}\cos E(t))$, where ${e}_{p}\ll 1$ and E is the eccentric anomaly of the particle orbit. Since we are only interested in terms that are no more than quadratic in eccentricity (of both particle and disk), we next expand the integrals in a Taylor series over ${r}_{p}-{a}_{p}=-{a}_{p}{e}_{p}\cos E(t)$ up to second order and average over the planetesimal orbit using the relations $\langle \cos E\rangle =-{e}_{p}/2$, $\langle {\cos }^{2}E\rangle =1/2$. Note that the dependence on r_p is also present in the limits of integration, which, upon this Taylor expansion, give rise to the nonintegral boundary terms. We omit the tedious process of writing out the Taylor expansion and subsequent averaging and present just the result,

$$\begin{eqnarray}{I}_{1}^{(0)} & \approx & \pi {{Ga}}_{p}\displaystyle \frac{{e}_{p}^{2}}{2}{\displaystyle \int }_{{a}_{1}/{a}_{p}}^{1}\left[\displaystyle \frac{d}{{{dr}}_{p}}(\alpha {r}_{p}{{\rm{\Sigma }}}_{d}(\alpha {r}_{p}))\right.{\left.\left.+\displaystyle \frac{{a}_{p}}{2}\displaystyle \frac{{d}^{2}}{{{dr}}_{p}}(\alpha {r}_{p}{{\rm{\Sigma }}}_{d}(\alpha {r}_{p}))\right]\right|}_{{r}_{p}={a}_{p}}{b}_{1/2}^{(0)}(\alpha )d\alpha +\pi G\displaystyle \frac{{e}_{p}^{2}}{4}\\ & & \times \,\left[-\displaystyle \frac{{a}_{1}^{3}}{{a}_{p}^{2}}{{\rm{\Sigma }}}_{d}({a}_{1}){b}_{1/2}^{(0)^{\prime} }\left(\displaystyle \frac{{a}_{1}}{{a}_{p}}\right)\right.+\displaystyle \frac{{a}_{1}^{3}}{{a}_{p}}{{\rm{\Sigma }}}_{d}^{{\prime} }({a}_{1}){b}_{1/2}^{(0)}\left(\displaystyle \frac{{a}_{1}}{{a}_{p}}\right)+\left.\displaystyle \frac{{a}_{1}^{2}}{{a}_{p}}{{\rm{\Sigma }}}_{d}({a}_{1}){b}_{1/2}^{(0)}\left(\displaystyle \frac{{a}_{1}}{{a}_{p}}\right)\right],\end{eqnarray} \tag{ 27 }$$

where we dropped the terms independent of e_p. Such terms vanish when substituted into the Lagrange evolution equations. For ${I}_{1}^{(1)}$, we similarly have

$$\begin{eqnarray}&&{I}_{1}^{(1)}\approx \pi {{Ga}}_{p}^{2}\displaystyle \frac{{e}_{p}^{2}}{2}{\displaystyle \int }_{{a}_{p}/{a}_{2}}^{1}\left[\displaystyle \frac{1}{{r}_{p}}\displaystyle \frac{d}{{{dr}}_{p}}({r}_{p}{{\rm{\Sigma }}}_{d}({\alpha }^{-1}{r}_{p}))\right.{\left.\left.+\displaystyle \frac{1}{2}\displaystyle \frac{{d}^{2}}{{{dr}}_{p}^{2}}({r}_{p}{{\rm{\Sigma }}}_{d}({\alpha }^{-1}{r}_{p}))\right]\right|}_{{r}_{p}={a}_{p}}{\alpha }^{-2}{b}_{1/2}^{(0)}(\alpha )d\alpha \,-\,\pi G\displaystyle \frac{{e}_{p}^{2}}{4}\\ &&\,\times \,\left[2{a}_{2}{{\rm{\Sigma }}}_{d}({a}_{2}){b}_{1/2}^{(0)}\left(\displaystyle \frac{{a}_{p}}{{a}_{2}}\right)\right.\left.\,+\,{a}_{p}{{\rm{\Sigma }}}_{d}({a}_{2}){b}_{1/2}^{(0)^{\prime} }\left(\displaystyle \frac{{a}_{p}}{{a}_{2}}\right)+{a}_{2}^{2}{{\rm{\Sigma }}}_{d}^{{\prime} }({a}_{2}){b}_{1/2}^{(0)}\left(\displaystyle \frac{{a}_{p}}{{a}_{2}}\right)\right].\end{eqnarray} \tag{ 28 }$$

Now we turn to calculation of the integral contribution I₂ containing ${\varphi }_{d}=\theta +v$, where we defined $v={\varphi }_{p}+{\varpi }_{p}-{\varpi }_{d}$. Let us introduce a new auxiliary function, $Q({r}_{d})=d[e({r}_{d}){{\rm{\Sigma }}}_{d}({r}_{d})]/{{dr}}_{d}.$ Then,

$$\begin{eqnarray}{I}_{2} & = & G\langle {\displaystyle \int }_{{r}_{d,\mathrm{in}}}^{{r}_{d,\mathrm{out}}}Q({r}_{d}){r}_{d}^{2}{\displaystyle \int }_{0}^{2\pi }\displaystyle \frac{\cos \theta \cos {vd}\theta }{\sqrt{{r}_{p}^{2}+{r}_{d}^{2}-2{r}_{p}{r}_{d}\cos \theta }}{{dr}}_{d}\rangle ,\end{eqnarray} \tag{ 29 }$$

where we expanded $\cos (\theta +v)$ and took into account that the term proportional to $\sin \theta$ evaluates to zero.

Next, we divide the integral over dr_d into two parts, one for ${r}_{d}\lt {r}_{p}$ and another for ${r}_{d}\gt {r}_{p}$, just like it was done for I₁. We again approximate ${r}_{d,\mathrm{in}}\approx {a}_{1}$ and ${r}_{d,\mathrm{out}}\approx {a}_{2}$, as the linear corrections to these expressions would lead to a contribution, which is third order in eccentricity. We then again make a change of variables in both resulting integrals, setting $\alpha ={r}_{d}/{r}_{p}\lt 1$ in the first one and $\alpha ={r}_{p}/{r}_{d}\lt 1$ in the second. As a result, using the definition of the Laplace coefficient ${b}_{1/2}^{(1)}$, we obtain the following expression:

$$\begin{eqnarray}{I}_{2} & = & {I}_{2}^{(0)}+{I}_{2}^{(1)}=\pi G\left\langle {\displaystyle \int }_{{a}_{1}/{r}_{p}}^{1}\cos v\,{r}_{p}^{2}Q(\alpha {r}_{p}){\alpha }^{2}{b}_{1/2}^{(1)}(\alpha )d\alpha \right.\left.\,+{\displaystyle \int }_{{r}_{p}/{a}_{2}}^{1}\cos v\,{r}_{p}^{2}Q({\alpha }^{-1}{r}_{p}){\alpha }^{-3}{b}_{1/2}^{(1)}(\alpha )d\alpha \right\rangle .\end{eqnarray} \tag{ 30 }$$

There are two main differences in our subsequent expansion over r_p − a_p compared to the previous case. First, $Q({r}_{d})\sim O({e}_{d})$; thus, in the expansion, we only need to retain terms up to first order in e_p. Second, now we have the time dependence of the integrand, also through $\cos v(t)=\cos ({\rm{\Delta }}\varpi +{\varphi }_{p}(t))$ $=\,\cos ({\rm{\Delta }}\varpi +E(t)+{e}_{p}\sin E(t))$, where we defined ${\rm{\Delta }}\varpi ={\varpi }_{p}-{\varpi }_{d}$ and used the relation ${\varphi }_{p}=E+{e}_{p}\sin E$. The easiest way to deal with the calculation is to consider r_p as a function of e_p, ${r}_{p}={a}_{p}(1-{e}_{p}\cos E),$ and to expand the integrals directly over e_p. We start with ${I}_{2}^{(0)}$:

$$\begin{eqnarray}{I}_{2}^{(0)} & = & \pi G\langle {\displaystyle \int }_{{a}_{1}/({a}_{p}(1-{e}_{p}\cos E))}^{1}\cos ({\rm{\Delta }}\varpi +E+{e}_{p}\sin E){r}_{p}^{2}({e}_{p},E)Q(\alpha {r}_{p}({e}_{p},E)){\alpha }^{2}{b}_{1/2}^{(1)}(\alpha )d\alpha \rangle .\end{eqnarray} \tag{ 31 }$$

We need the zeroth and first-order (in ${e}_{p}\ll 1$) terms of the Taylor expansion of this expression at e_p = 0. We will omit the technicalities of writing out the expansion and present the result:

$$\begin{eqnarray}&&{I}_{2}^{(0)}\approx \pi G{\displaystyle \int }_{{a}_{1}/{a}_{p}}^{1}{a}_{p}^{2}Q(\alpha {a}_{p}){\alpha }^{2}{b}_{1/2}^{(1)}(\alpha )d\alpha \langle \cos ({\rm{\Delta }}\varpi +E)\rangle \,-\,\pi {{Ge}}_{p}\displaystyle \frac{{a}_{1}^{3}}{{a}_{p}}{b}_{1/2}^{(1)}\displaystyle \frac{{a}_{1}}{{a}_{p}}Q({a}_{1})\langle \cos ({\rm{\Delta }}\varpi +E)\cos (E)\rangle \\ &&\,-\,\pi {{Ge}}_{p}{\displaystyle \int }_{{a}_{1}/{a}_{p}}^{1}({a}_{p}^{2}Q(\alpha {a}_{p})\langle \sin ({\rm{\Delta }}\varpi +E)\sin (E)\rangle \,+\,\left.{a}_{p}{\left.\displaystyle \frac{d}{{{dr}}_{p}}({r}_{p}^{2}Q(\alpha {r}_{p}))\right|}_{{r}_{p}={a}_{p}}\langle \cos ({\rm{\Delta }}\varpi +E)\cos (E)\rangle \right)\,{\alpha }^{2}{b}_{1/2}^{(1)}(\alpha )d\alpha .\end{eqnarray} \tag{ 32 }$$

Using the averages of the trigonometric functions of E,

$$\begin{eqnarray}\langle \cos ({\rm{\Delta }}\varpi +E)\rangle & = & -\displaystyle \frac{{e}_{p}}{2}\cos {\rm{\Delta }}\varpi ,\,\,\langle \sin ({\rm{\Delta }}\varpi +E)\sin (E)\rangle \,= & \langle \cos ({\rm{\Delta }}\varpi +E)\cos (E)\rangle =\displaystyle \frac{1}{2}\cos {\rm{\Delta }}\varpi ,\end{eqnarray} \tag{ 33 }$$

based on $\langle \sin E\rangle =0$, $\langle \sin E\cos E\rangle =0$, and $\langle {\sin }^{2}E\rangle =1/2$, the final result for ${I}_{2}^{(0)}$ becomes

$$\begin{eqnarray}&&{I}_{2}^{(0)}\approx -\pi G\displaystyle \frac{{e}_{p}}{2}\cos {\rm{\Delta }}\varpi \left[{\int }_{{a}_{1}/{a}_{p}}^{1}(4Q(\alpha {a}_{p})\right.\left.\,+\,{a}_{p}{\left.\displaystyle \frac{d}{{{dr}}_{p}}(Q(\alpha {r}_{p}))\right|}_{{r}_{p}={a}_{p}}\right){a}_{p}^{2}{\alpha }^{2}{b}_{1/2}^{(1)}(\alpha )d\alpha \left.\,+\,{a}_{1}^{2}{b}_{1/2}^{(1)}\left(\displaystyle \frac{{a}_{p}}{{a}_{1}}\right)Q({a}_{1})\right].\end{eqnarray} \tag{ 34 }$$

Following a similar procedure, for ${I}_{2}^{(1)}$, we have

$$\begin{eqnarray}{I}_{2}^{(1)} & \cong & -\pi G\displaystyle \frac{{e}_{p}}{2}\cos {\rm{\Delta }}\varpi [{\int }_{{a}_{p}/{a}_{2}}^{1}(4Q\,\displaystyle (\frac{{a}_{p}}{\alpha })\,+{a}_{p}{\left.\displaystyle \frac{d}{{{dr}}_{p}}\left(Q\left(\displaystyle \frac{{r}_{p}}{\alpha }\right)\right)\right|}_{{r}_{p}={a}_{p}})\displaystyle \frac{{a}_{p}^{2}}{{\alpha }^{3}}{b}_{1/2}^{(1)}(\alpha )d\alpha \,-{a}_{2}^{2}{b}_{1/2}^{(1)}\displaystyle (\frac{{a}_{p}}{{a}_{2}})Q({a}_{2})].\end{eqnarray} \tag{ 35 }$$

Equations (27), (28), (34), and (35) provide the desired contribution to the disturbing function from the circular annulus ${{\mathbb{D}}}_{c}$, since $R({{\mathbb{D}}}_{c})={I}_{1}^{(0)}+{I}_{1}^{(1)}+{I}_{2}^{(0)}+{I}_{2}^{(1)}$.

Now let us turn to the inner crescent ${{\mathbb{D}}}_{i}$. As the eccentricity of our disk is small, the width of this crescent is already $O({e}_{d})$; thus, we only need to keep the terms up to first order in e_p in our expansion. Following the recipe in Silsbee & Rafikov (2015b), we write

$$\begin{eqnarray}&&R({{\mathbb{D}}}_{i})=-\left\langle G{\displaystyle \int }_{{a}_{1}(1-{e}_{1})}^{{a}_{1}(1+{e}_{1})}\displaystyle \frac{{r}_{d}}{{r}_{p}}{\rm{\Sigma }}({r}_{d}){{dr}}_{d}\right.\left.{\displaystyle \int }_{\xi -{\rm{\Delta }}\varpi -{\varphi }_{p}}^{2\pi -\xi -{\rm{\Delta }}\varpi -{\varphi }_{p}}\displaystyle \frac{d\theta }{\sqrt{1+{\alpha }^{2}-2\alpha \cos \theta }}\right\rangle ,\end{eqnarray} \tag{ 36 }$$

where $\xi ({r}_{d})=\arccos (({a}_{1}-{r}_{d})/{e}_{1}{a}_{1})$ determines the azimuthal extent of the crescent for a given r_d (Silsbee & Rafikov 2015b), and $\alpha ={r}_{d}/{r}_{p}$.

For the function in the inner integral, we can use the Fourier expansion,

$$\begin{eqnarray}&&{(1+{\alpha }^{2}-2\alpha \cos \theta )}^{-1/2}=\displaystyle \frac{1}{2}{b}_{1/2}^{(0)}(\alpha )\,+\,\displaystyle \sum _{j=1}^{\infty }{b}_{1/2}^{(j)}(\alpha )\cos (j\theta ),\end{eqnarray} \tag{ 37 }$$

letting us calculate the inner integral as

$$\begin{eqnarray}&&{\displaystyle \int }_{\xi -{\rm{\Delta }}\varpi -{\varphi }_{p}}^{2\pi -\xi -{\rm{\Delta }}\varpi -{\varphi }_{p}}\left(\displaystyle \frac{1}{2}{b}_{1/2}^{(0)}(\alpha )+\displaystyle \sum _{j=1}^{\infty }{b}_{1/2}^{(j)}(\alpha )\cos (j\theta )\right)d\theta \,=\,(\pi -\xi ({r}_{d})){b}_{1/2}^{(0)}(\alpha )+{\left.\displaystyle \sum _{j=1}^{\infty }\displaystyle \frac{1}{j}{b}_{1/2}^{(j)}(\alpha )\sin (j\theta )\right|}_{\xi -{\rm{\Delta }}\varpi -{\varphi }_{p}}^{2\pi -\xi -{\rm{\Delta }}\varpi -{\varphi }_{p}}\\ &&=\,(\pi -\xi ({r}_{d})){b}_{1/2}^{(0)}(\alpha )-\displaystyle \sum _{j=1}^{\infty }\displaystyle \frac{2}{j}{b}_{1/2}^{(j)}(\alpha )\sin \ (j\xi ({r}_{d}))\cos (j({\rm{\Delta }}\varpi +{\varphi }_{p})).\end{eqnarray} \tag{ 38 }$$

The next step is to remember that we only need to keep terms up to $O({e}_{d})$ in Equation (36) and to notice that the interval over which radial integration is performed is already $O({e}_{d})$. As a result, in Equation (36), we can set ${r}_{d}\approx {a}_{1}$ and ${{\rm{\Sigma }}}_{d}({r}_{d})\approx {{\rm{\Sigma }}}_{d}({a}_{1})$ and take them out of the integral. We can also set $\alpha \approx {a}_{1}/{r}_{p}$ in Equation (38) for the same reason. As a result, we find

$$\begin{eqnarray}&&R({{\mathbb{D}}}_{i})\approx -G{\rm{\Sigma }}({a}_{1})\left\langle \displaystyle \frac{{a}_{1}}{{r}_{p}}{\displaystyle \int }_{{a}_{1}(1-{e}_{1})}^{{a}_{1}(1+{e}_{1})}\left[(\pi -\xi ({r}_{d})){b}_{1/2}^{(0)}\left(\displaystyle \frac{{a}_{1}}{{r}_{p}}\right)\right.\right.\left.\left.\,-\,\displaystyle \sum _{j=1}^{\infty }\displaystyle \frac{2}{j}{b}_{1/2}^{(j)}\left(\displaystyle \frac{{a}_{1}}{{r}_{p}}\right)\sin (j\xi ({r}_{d}))\cos (j({\rm{\Delta }}\varpi +{\varphi }_{p}))\right]{{dr}}_{d}\right\rangle \\ &&\,=\,-\pi {{Ge}}_{1}{a}_{1}^{2}{\rm{\Sigma }}({a}_{1})\left\langle {r}_{p}^{-1}{b}_{1/2}^{(0)}\left(\displaystyle \frac{{a}_{1}}{{r}_{p}}\right)-{r}_{p}^{-1}{b}_{1/2}^{(1)}\right.\left.\left(\displaystyle \frac{{a}_{1}}{{r}_{p}}\right)\cos ({\rm{\Delta }}\varpi +{\varphi }_{p})\right\rangle ,\end{eqnarray} \tag{ 39 }$$

where we have used the fact that

$$\begin{eqnarray}&&{\displaystyle \int }_{{a}_{1}(1-{e}_{1})}^{{a}_{1}(1+{e}_{1})}\sin (j\xi ({r}_{d})){{dr}}_{d}\,=\,{e}_{1}{a}_{1}{\displaystyle \int }_{0}^{\pi }\sin (j\xi )\sin (\xi )d\xi ={\delta }_{j1}\displaystyle \frac{\pi }{2}{e}_{1}{a}_{1},\end{eqnarray} \tag{ 40 }$$

where ${\delta }_{{ij}}$ is the Kronecker delta.

Now we just need to expand Equation (39) up to first order in e_p and then orbit average the resultant expression, which is straightforward using Equation (33)

$$\begin{eqnarray}&&\left\langle {r}_{p}^{-1}\left[{b}_{1/2}^{(0)}\left(\displaystyle \frac{{a}_{1}}{{r}_{p}}\right)-{b}_{1/2}^{(1)}\left(\displaystyle \frac{{a}_{1}}{{r}_{p}}\right)\cos ({\rm{\Delta }}\varpi +{\varphi }_{p})\right]\right\rangle \approx \,{a}_{p}^{-1}\left\langle {b}_{1/2}^{(0)}\left(\displaystyle \frac{{a}_{1}}{{a}_{p}}\right)-{b}_{1/2}^{(1)}\left(\displaystyle \frac{{a}_{1}}{{a}_{p}}\right)\cos ({\rm{\Delta }}\varpi +E)\right\rangle \\ &&+\,{a}_{p}{e}_{p}\left\langle \cos E\,\left[{a}_{p}^{-2}{b}_{1/2}^{(0)}\left(\displaystyle \frac{{a}_{1}}{{a}_{p}}\right)+{a}_{p}^{-1}{a}_{1}{a}_{p}^{-2}{b}_{1/2}^{(0)^{\prime} }\left(\displaystyle \frac{{a}_{1}}{{a}_{p}}\right)\right]\right\rangle -\,\left\langle {a}_{p}{e}_{p}\cos E\,{a}_{p}^{-2}{b}_{1/2}^{(1)}\left(\displaystyle \frac{{a}_{1}}{{a}_{p}}\right)\cos ({\rm{\Delta }}\varpi +E)\right.\\ &&+\,{a}_{p}{e}_{p}\cos E\,{a}_{p}^{-1}{a}_{1}{a}_{p}^{-2}{b}_{1/2}^{(1)^{\prime} }\left(\displaystyle \frac{{a}_{1}}{{a}_{p}}\right)\cos ({\rm{\Delta }}\varpi +E)\left.\,-\,{e}_{p}{a}_{p}^{-1}{b}_{1/2}^{(1)}\left(\displaystyle \frac{{a}_{1}}{{a}_{p}}\right)\sin E\,\sin ({\rm{\Delta }}\varpi +E)\right\rangle \\ &&=\,{a}_{p}^{-1}{b}_{1/2}^{(0)}\left(\displaystyle \frac{{a}_{1}}{{a}_{p}}\right)+\displaystyle \frac{{e}_{p}}{2}{a}_{p}^{-1}\left[{b}_{1/2}^{(1)}\left(\displaystyle \frac{{a}_{1}}{{a}_{p}}\right)\right.\left.\,-\,\displaystyle \frac{{a}_{1}}{{a}_{p}}{b}_{1/2}^{(1)^{\prime} }\left(\displaystyle \frac{{a}_{1}}{{a}_{p}}\right)\right]\cos {\rm{\Delta }}\varpi .\end{eqnarray} \tag{ 41 }$$

As a result, by dropping the first term not containing e_p, one finally finds

$$\begin{eqnarray}R({{\mathbb{D}}}_{i}) & \approx & -\pi G{\rm{\Sigma }}({a}_{1}){e}_{1}\displaystyle \frac{{e}_{p}}{2}\displaystyle \frac{{a}_{1}^{2}}{{a}_{p}}\left[{b}_{1/2}^{(1)}\left(\displaystyle \frac{{a}_{1}}{{a}_{p}}\right)\right.\left.\,-\displaystyle \frac{{a}_{1}}{{a}_{p}}{b}_{1/2}^{(1)^{\prime} }\left(\displaystyle \frac{{a}_{1}}{{a}_{p}}\right)\right]\cos {\rm{\Delta }}\varpi .\end{eqnarray} \tag{ 42 }$$

This case has only minor differences from the previous one, so we will omit the details:

$$\begin{eqnarray}&&R({{\mathbb{D}}}_{o})=\left\langle G{\displaystyle \int }_{{a}_{2}(1-{e}_{2})}^{{a}_{2}(1+{e}_{2})}{\rm{\Sigma }}({r}_{d}){{dr}}_{d}{\displaystyle \int }_{\xi -{\rm{\Delta }}\varpi -{\varphi }_{p}}^{2\pi -\xi -{\rm{\Delta }}\varpi -{\varphi }_{p}}\right.\left.\displaystyle \frac{d\theta }{\sqrt{1+{\alpha }^{2}-2\alpha \cos \theta }}\right\rangle \,\approx \,\pi {{Ge}}_{2}{a}_{2}{\rm{\Sigma }}({a}_{2})\left\langle {b}_{1/2}^{(0)}\left(\displaystyle \frac{{r}_{p}}{{a}_{2}}\right)\right.\left.\,-\,{b}_{1/2}^{(1)}\left(\displaystyle \frac{{r}_{p}}{{a}_{2}}\right)\cos ({\rm{\Delta }}\varpi +{\varphi }_{p})\right\rangle \\ &&\,\approx \,\pi {{Ge}}_{2}{a}_{2}{\rm{\Sigma }}({a}_{2})\left[{b}_{1/2}^{(0)}\left(\displaystyle \frac{{a}_{p}}{{a}_{2}}\right)+{e}_{p}{b}_{1/2}^{(1)}\left(\displaystyle \frac{{a}_{p}}{{a}_{2}}\right)\right.\,\left.\cos {\rm{\Delta }}\varpi +\displaystyle \frac{{e}_{p}}{2}{a}_{p}{a}_{2}^{-1}{b}_{1/2}^{(1)^{\prime} }\left(\displaystyle \frac{{a}_{p}}{{a}_{2}}\right)\cos {\rm{\Delta }}\varpi \right].\end{eqnarray} \tag{ 43 }$$

Thus, again dropping the first, e_p-independent term, we eventually arrive at

$$\begin{eqnarray}R({{\mathbb{D}}}_{o}) & \approx & \pi G{\rm{\Sigma }}({a}_{2}){e}_{2}\displaystyle \frac{{e}_{p}}{2}{a}_{2}\left[2{b}_{1/2}^{(1)}\left(\displaystyle \frac{{a}_{p}}{{a}_{2}}\right)\right.\left.+\,\displaystyle \frac{{a}_{p}}{{a}_{2}}{b}_{1/2}^{(1)^{\prime} }\left(\displaystyle \frac{{a}_{p}}{{a}_{2}}\right)\right]\cos {\rm{\Delta }}\varpi .\end{eqnarray} \tag{ 44 }$$

The full expression resulting from evaluating the integral in Equation (21) is given by ${R}_{d}={I}_{1}^{(0)}+{I}_{1}^{(1)}+{I}_{2}^{(0)}+{I}_{2}^{(1)}+R({{\mathbb{D}}}_{i})+R({{\mathbb{D}}}_{o})$ and can be written out explicitly using Equations (27), (28), (34), (35), (42), and (44). When doing this, we introduce two important modifications.

First, it is possible to combine the different terms corresponding to the outer and inner disks into a single expression using certain properties of the Laplace coefficients. For example, we can combine ${I}_{1}^{(0)}$ and ${I}_{1}^{(1)}$ by changing the variable $\alpha \to {\alpha }^{-1}$ in Equation (28) for ${I}_{1}^{(1)}$ and then manipulate the result using the fact that ${b}_{s}^{(j)}({\alpha }^{-1})={\alpha }^{2s}{b}_{s}^{(j)}(\alpha )$. This procedure turns ${I}_{1}^{(1)}$ in Equation (28) into the expression analogous to Equation (27) but with the limits of integration running from 1 to ${a}_{2}/{a}_{p}$. This allows us to combine it with ${I}_{1}^{(0)}$ to get an expression analogous to Equation (27) with integration running from ${a}_{1}/{a}_{p}\lt 1$ to ${a}_{2}/{a}_{p}\gt 1$ (for orbits within the disk, ${a}_{1}\lt {a}_{p}\lt {a}_{2}$). The integrand of this final expression has a weak (logarithmic) singularity at $\alpha =1$, but the integral itself is fully convergent. This procedure is then repeated to combine ${I}_{2}^{(0)}$ and ${I}_{2}^{(1)}$. In addition, differentiating the identity ${b}_{s}^{(j)}({\alpha }^{-1})={\alpha }^{2s}{b}_{s}^{(j)}(\alpha )$, we obtain the relation for ${b}_{s}^{(j)^{\prime} }({\alpha }^{-1})$, which allows us to also merge $R({{\mathbb{D}}}_{i})$ and $R({{\mathbb{D}}}_{o})$ into a single expression.

Second, in all these expressions, we switch from d/dr_p to derivatives with respect to the argument (denoted by prime), e.g., ${dQ}(\alpha {r}_{p})/{{dr}}_{p}=\alpha {Q}^{{\prime} }(a){| }_{a=\alpha {r}_{p}}$.

As a result of performing these steps, one finds that R_d can be written in the form of Equation (3) with the coefficients A_d and B_d given by the following expressions:

$$\begin{eqnarray}&&{A}_{d}={A}_{d}^{\mathrm{bulk}}+{A}_{d}^{\mathrm{edge}},\,\,\,\,{B}_{d}={B}_{d}^{\mathrm{bulk}}+{B}_{d}^{\mathrm{edge}},\end{eqnarray} \tag{ 45 }$$

$$\begin{eqnarray}&&{A}_{d}^{\mathrm{bulk}}=\displaystyle \frac{\pi G}{2{a}_{p}{n}_{p}}{\int }_{{a}_{1}/{a}_{p}}^{{a}_{2}/{a}_{p}}\alpha {b}_{1/2}^{(0)}(\alpha ){({a}^{2}{{\rm{\Sigma }}}_{d}(a))}^{{\prime\prime} }{| }_{a=\alpha {a}_{p}}d\alpha ,\end{eqnarray} \tag{ 46 }$$

$$\begin{eqnarray}{A}_{d}^{\mathrm{edge}} & = & -\displaystyle \frac{\pi G}{2{a}_{p}{n}_{p}}\left[{\left(\displaystyle \frac{a}{{a}_{p}}\right)}^{2}{b}_{1/2}^{(0)}\left(\displaystyle \frac{a}{{a}_{p}}\right){(a{{\rm{\Sigma }}}_{d}(a))}^{{\prime} }\right.-{\left.\left.{\left(\displaystyle \frac{a}{{a}_{p}}\right)}^{3}{b}_{1/2}^{(0)^{\prime} }\left(\displaystyle \frac{a}{{a}_{p}}\right){{\rm{\Sigma }}}_{d}(a)\right]\right|}_{a={a}_{1}}^{a={a}_{2}},\end{eqnarray} \tag{ 47 }$$

$$\begin{eqnarray}{B}_{d}^{\mathrm{bulk}} & = & -\displaystyle \frac{\pi G}{2{n}_{p}}{\displaystyle \int }_{{a}_{1}/{a}_{p}}^{{a}_{2}/{a}_{p}}{\alpha }^{2}{b}_{1/2}^{(1)}(\alpha )[4{({{\rm{\Sigma }}}_{d}(a){e}_{d}(a))}^{{\prime} }+a{({{\rm{\Sigma }}}_{d}(a){e}_{d}(a))}^{{\prime\prime} }]{| }_{a=\alpha {a}_{p}}d\alpha ,\end{eqnarray} \tag{ 48 }$$

and

$$\begin{eqnarray}&&{B}_{d}^{\mathrm{edge}}=\displaystyle \frac{\pi G}{2{a}_{p}^{3}{n}_{p}}\left[{a}^{3}{b}_{1/2}^{(1)}\left(\displaystyle \frac{a}{{a}_{p}}\right){({{\rm{\Sigma }}}_{d}(a){e}_{d}(a))}^{{\prime} }\right.\left.\,+\,{a}^{2}\left({b}_{1/2}^{(1)}\left(\displaystyle \frac{a}{{a}_{p}}\right)-\displaystyle \frac{a}{{a}_{p}}{b}_{1/2}^{(1)^{\prime} }\left(\displaystyle \frac{a}{{a}_{p}}\right)\right){{\rm{\Sigma }}}_{d}(a){e}_{d}(a)\right]{| }_{a={a}_{1}}^{a={a}_{2}},\end{eqnarray} \tag{ 49 }$$

where $F(a){| }_{a={a}_{1}}^{a={a}_{2}}\equiv F(a={a}_{2})-F(a={a}_{1})$. If we consider a disk model with power-law profiles of ${{\rm{\Sigma }}}_{d}(a)\propto {a}^{-p}$ and ${e}_{d}(a)\propto {a}^{-q}$, Equations (45)–(49) reduce to the corresponding formulae in Silsbee & Rafikov (2015b).