Potential softening and eccentricity dynamics in razor-thin,   nearly-Keplerian discs

Antranik A. Sefilian; Roman R. Rafikov

arXiv:1904.07592·astro-ph.EP·September 25, 2019

Potential softening and eccentricity dynamics in razor-thin, nearly-Keplerian discs

Antranik A. Sefilian, Roman R. Rafikov

PDF

TL;DR

This paper evaluates various softening methods in secular disc dynamics, identifying models that accurately reproduce eccentricity behavior and highlighting the computational demands for precise simulations, especially near disc edges.

Contribution

It develops a general framework for computing the secular disturbing function with arbitrary softening, and assesses the convergence and accuracy of existing softening models in disc eccentricity dynamics.

Findings

01

Some softening models converge to correct behavior as softening approaches zero.

02

Accurate eccentricity dynamics require a large number of interacting elements, especially for small softening.

03

Very small softening is necessary near disc edges to capture boundary dynamics accurately.

Abstract

In many astrophysical problems involving discs (gaseous or particulate) orbiting a dominant central mass, gravitational potential of the disc plays an important dynamical role. Its impact on the motion of external objects, as well as on the dynamics of the disc itself, can usually be studied using secular approximation. This is often done using softened gravity to avoid singularities arising in calculation of the orbit-averaged potential --- disturbing function --- of a razor-thin disc using classical Laplace-Lagrange theory. We explore the performance of several softening formalisms proposed in the literature in reproducing the correct eccentricity dynamics in the disc potential. We identify softening models that, in the limit of zero softening, give results converging to the expected behavior exactly, approximately or not converging at all. We also develop a general framework for…

Tables2

Table 1. Table 1: The coefficients ϕ i j ( α ) subscript italic-ϕ 𝑖 𝑗 𝛼 \phi_{ij}(\alpha) of the secular disturbing function with softened gravity featured in Eqs. ( 5 )-( 6 ), which govern the individual secular ring-ring interaction (Eq. 11 ), adopted from the literature (listed in the first column). Here α 𝛼 \alpha is defined such that α = a < / a > 𝛼 subscript 𝑎 subscript 𝑎 \alpha=a_{<}/a_{>} where a > = max ( a 1 , a 2 ) subscript 𝑎 max subscript 𝑎 1 subscript 𝑎 2 a_{>}=\text{max}(a_{1},a_{2}) and a < = min ( a 1 , a 2 ) subscript 𝑎 min subscript 𝑎 1 subscript 𝑎 2 a_{<}=\text{min}(a_{1},a_{2}) . The softened interactions under consideration are those of Tremaine ( 1998 ) , Touma ( 2002 ) , Hahn ( 2003 ) and Teyssandier & Ogilvie ( 2016 ) – see § 2.1 for further details. For reference, the expressions of ϕ i j LL superscript subscript italic-ϕ 𝑖 𝑗 LL \phi_{ij}^{\rm LL} corresponding to the (unsoftened) Newtonian ring-ring interaction (i.e. classical Laplace-Lagrange formalism) are also shown in the top row. The Laplace coefficients which are softened by the introduction of a softening parameter ϵ 2 ( α ) superscript italic-ϵ 2 𝛼 \epsilon^{2}(\alpha) are defined in Eq. ( 7 ). Note that the expressions of ϕ i j subscript italic-ϕ 𝑖 𝑗 \phi_{ij} reported in Touma ( 2002 ) have been corrected in a subsequent paper of Touma & Sridhar ( 2012 ) .

Formalism	$ϵ^{2} (α)$	$ϕ_{11}$	$ϕ_{12}$	$ϕ_{22}$
Laplace-Lagrange	–	$\frac{1}{8} α b_{3 / 2}^{(1)}$	$- \frac{1}{4} α b_{3 / 2}^{(2)}$	$ϕ_{11}$
Tremaine (1998) (Tr98)	$β_{c}^{2}$	$\frac{1}{8} (2 α \frac{d}{d α} + α^{2} \frac{d^{2}}{d α^{2}}) ℬ_{1 / 2}^{(0), Tr}$	$\frac{1}{4} (2 - 2 α \frac{d}{d α} - α^{2} \frac{d^{2}}{d α^{2}}) ℬ_{1 / 2}^{(1), Tr}$	$ϕ_{11}$
		$= \frac{1}{8} α [ℬ_{3 / 2}^{(1), Tr} - 3 α β_{c}^{2} ℬ_{5 / 2}^{(0), Tr}]$	$= - \frac{1}{4} α [ℬ_{3 / 2}^{(2), Tr} - 3 α β_{c}^{2} ℬ_{5 / 2}^{(1), Tr}]$
Touma (2002) (T02)	$β^{2} = b_{c}^{2} / a_{>}^{2}$	$- \frac{5}{8} α ℬ_{3 / 2}^{(1), T} + \frac{3}{16} α^{2} ℬ_{5 / 2}^{(0), T} + \frac{3}{8} α (1 + α^{2}) ℬ_{5 / 2}^{(1), T}$	$\frac{9}{8} α ℬ_{3 / 2}^{(0), T} + \frac{1}{8} α ℬ_{3 / 2}^{(2), T} - \frac{9}{8} α (1 + α^{2}) ℬ_{5 / 2}^{(0), T}$	$- \frac{5}{8} α ℬ_{3 / 2}^{(1), T} + \frac{3}{16} α^{2} ℬ_{5 / 2}^{(0), T} + \frac{3}{8} α (1 + α^{2}) ℬ_{5 / 2}^{(1), T}$
		$- \frac{15}{16} α^{2} ℬ_{5 / 2}^{(2), T} - \frac{3}{8} α β^{2} (α ℬ_{5 / 2}^{(0), T} - ℬ_{5 / 2}^{(1), T})$	$+ \frac{21}{16} α^{2} ℬ_{5 / 2}^{(1), T} + \frac{3}{8} α (1 + α^{2}) ℬ_{5 / 2}^{(2), T} + \frac{3}{16} α^{2} ℬ_{5 / 2}^{(3), T}$	$- \frac{15}{16} α^{2} ℬ_{5 / 2}^{(2), T} - \frac{3}{8} β^{2} (ℬ_{5 / 2}^{(0), T} - α ℬ_{5 / 2}^{(1), T})$
Hahn (2003) (H03)	$H^{2} (1 + α^{2})$	$\frac{1}{8} α [ℬ_{3 / 2}^{(1), H} - 3 α H^{2} (2 + H^{2}) ℬ_{5 / 2}^{(0), H}]$	$- \frac{1}{4} α [ℬ_{3 / 2}^{(2), H} - 3 α H^{2} (2 + H^{2}) ℬ_{5 / 2}^{(1), H}]$	$ϕ_{11}$
Teyssandier & Ogilvie (2016) (TO16)	$S^{2} α$	$\frac{1}{8} α ℬ_{3 / 2}^{(1), TO}$	$- \frac{1}{4} α ℬ_{3 / 2}^{(2), TO}$	$ϕ_{11}$

Table 2. Table 2: The functional forms of the coefficients T i ( α ) subscript 𝑇 𝑖 𝛼 T_{i}(\alpha) given by Eqs. ( 59 )-( 61 ) appearing in the orbit-averaged disturbing function due to two coplanar (arbitrarily) softened rings (Eq. 56 - 58 ) such that α ≡ a < / a > ≤ 1 𝛼 subscript 𝑎 subscript 𝑎 1 \alpha\equiv a_{<}/a_{>}\leq 1 . The first column lists the softening prescriptions analyzed in this work (see § 2.1 ), while the second column shows the specific forms of the softening function ℱ ( r 1 , r 2 ) ℱ subscript 𝑟 1 subscript 𝑟 2 \mathcal{F}(r_{1},r_{2}) in Eq. ( 35 ). The corresponding expressions for the dimensionless softening parameters ϵ 2 ( α ) = a > − 2 ℱ ( a 1 , a 2 ) superscript italic-ϵ 2 𝛼 superscript subscript 𝑎 2 ℱ subscript 𝑎 1 subscript 𝑎 2 \epsilon^{2}(\alpha)=a_{>}^{-2}\mathcal{F}(a_{1},a_{2}) (Eq. 55 ) entering in the definition of softened Laplace coefficients (Eq. 54 ) are also shown. Here, Θ ( x ) Θ 𝑥 \Theta(x) represents the Heaviside step function and δ ( x ) = d Θ ( x ) / d x 𝛿 𝑥 𝑑 Θ 𝑥 𝑑 𝑥 \delta(x)=d\Theta(x)/dx stands for Dirac delta-function.

Method	$ℱ (r_{1}, r_{2})$	$ϵ^{2} (α)$	$T_{1} (α)$	$T_{2} (α)$	$T_{3} (α)$	$T_{4} (α)$	$T_{5} (α)$	$T_{6} (α)$	$T_{7} (α)$	$T_{8} (α)$
H03	$H^{2} (r_{1}^{2} + r_{2}^{2})$	$H^{2} (1 + α^{2})$	$0$	$2 H^{2} (1 + α^{2})$	$0$	$0$	$4 α H^{2} (2 + H^{2})$	$0$	$2 H^{2} \frac{a_{1}}{a_{2}}$	$2 H^{2} \frac{a_{2}}{a_{1}}$
T02	$b_{c}^{2}$	$β^{2} = {(b_{c} / a_{>})}^{2}$	$2 β^{2}$	$0$	$0$	$0$	$0$	$0$	$0$	$0$
Tr98	$β_{c}^{2} \max (r_{1}^{2}, r_{2}^{2})$	$β_{c}^{2}$	$0$	$2 β_{c}^{2}$	$- 2 β_{c}^{2} δ (α - 1)$	$0$	$4 α β_{c}^{2}$	$0$	$\frac{2 β_{c}^{2}}{α} Θ (a_{1} - a_{2})$	$\frac{2 β_{c}^{2}}{α} Θ (a_{2} - a_{1})$
TO16	$S^{2} r_{1} r_{2}$	$S^{2} α$	$0$	$2 α S^{2}$	$S^{2}$	$0$	$S^{2} (S^{2} α + 2 α^{2} + 2)$	$0$	$S^{2}$	$S^{2}$

Equations131

\delta R=\frac{G\delta m_{d}a_{p}}{a^{2}}\bigg{[}\frac{1}{8}b_{3/2}^{(1)}\bigg{(}\frac{a_{p}}{a}\bigg{)}e_{p}^{2}-\frac{1}{4}b_{3/2}^{(2)}\bigg{(}\frac{a_{p}}{a}\bigg{)}e_{p}e_{d}\cos(\varpi_{p}-\varpi_{d})\bigg{]},

\delta R=\frac{G\delta m_{d}a_{p}}{a^{2}}\bigg{[}\frac{1}{8}b_{3/2}^{(1)}\bigg{(}\frac{a_{p}}{a}\bigg{)}e_{p}^{2}-\frac{1}{4}b_{3/2}^{(2)}\bigg{(}\frac{a_{p}}{a}\bigg{)}e_{p}e_{d}\cos(\varpi_{p}-\varpi_{d})\bigg{]},

b_{s}^{(m)}(\alpha)=\frac{2}{\pi}\int\limits_{0}^{\pi}\cos(m\theta)\bigg{[}1+\alpha^{2}-2\alpha\cos\theta\bigg{]}^{-s}d\theta,

b_{s}^{(m)}(\alpha)=\frac{2}{\pi}\int\limits_{0}^{\pi}\cos(m\theta)\bigg{[}1+\alpha^{2}-2\alpha\cos\theta\bigg{]}^{-s}d\theta,

R_{d} = \int_{a_{in}}^{a_{out}} δ R,

R_{d} = \int_{a_{in}}^{a_{out}} δ R,

R_{d}={n_{p}a_{p}^{2}}\bigg{[}\frac{1}{2}\mathsf{A}_{d}(a_{p})\mathbf{e}_{p}^{2}+\bm{\mathsf{B}}_{d}(a_{p})\cdot\mathbf{e}_{p}\bigg{]}.

R_{d}={n_{p}a_{p}^{2}}\bigg{[}\frac{1}{2}\mathsf{A}_{d}(a_{p})\mathbf{e}_{p}^{2}+\bm{\mathsf{B}}_{d}(a_{p})\cdot\mathbf{e}_{p}\bigg{]}.

A_{d} (a_{p}) = \frac{2 G}{n _{p} a _{p}^{3}}

A_{d} (a_{p}) = \frac{2 G}{n _{p} a _{p}^{3}}

B_{d} (a_{p}) = \frac{G}{n _{p} a _{p}^{3}}

B_{d} (a_{p}) = \frac{G}{n _{p} a _{p}^{3}}

\mathcal{B}_{s}^{(m)}(\alpha,\epsilon)=\frac{2}{\pi}\int\limits_{0}^{\pi}\cos(m\theta)\bigg{[}1+\alpha^{2}-2\alpha\cos\theta+\epsilon^{2}(\alpha)\bigg{]}^{-s}d\theta.

\mathcal{B}_{s}^{(m)}(\alpha,\epsilon)=\frac{2}{\pi}\int\limits_{0}^{\pi}\cos(m\theta)\bigg{[}1+\alpha^{2}-2\alpha\cos\theta+\epsilon^{2}(\alpha)\bigg{]}^{-s}d\theta.

ϕ_{11}^{LL} (α) = ϕ_{22}^{LL} (α) = \frac{1}{8} α b_{3/2}^{(1)} (α),

ϕ_{11}^{LL} (α) = ϕ_{22}^{LL} (α) = \frac{1}{8} α b_{3/2}^{(1)} (α),

ϕ_{12}^{LL} (α) = - \frac{1}{4} α b_{3/2}^{(2)} (α),

\Phi_{i}(\mathbf{r}_{1},\mathbf{r}_{2})=-{G}m_{j}\big{[}(\mathbf{r}_{1}-\mathbf{r}_{2})^{2}+\mathcal{F}(r_{1},r_{2})\big{]}^{-1/2},

\Phi_{i}(\mathbf{r}_{1},\mathbf{r}_{2})=-{G}m_{j}\big{[}(\mathbf{r}_{1}-\mathbf{r}_{2})^{2}+\mathcal{F}(r_{1},r_{2})\big{]}^{-1/2},

R_{i}=\frac{{G}m_{j}}{a_{>}}\bigg{[}\phi_{11}(\alpha)~{}e_{1}^{2}+\phi_{22}(\alpha)~{}e_{2}^{2}+\phi_{12}(\alpha)~{}e_{1}e_{2}\cos(\varpi_{1}-\varpi_{2})\bigg{]},

R_{i}=\frac{{G}m_{j}}{a_{>}}\bigg{[}\phi_{11}(\alpha)~{}e_{1}^{2}+\phi_{22}(\alpha)~{}e_{2}^{2}+\phi_{12}(\alpha)~{}e_{1}e_{2}\cos(\varpi_{1}-\varpi_{2})\bigg{]},

\mathcal{B}_{s}^{(m),\mathrm{Tr}}=\frac{2}{\pi}\int\limits_{0}^{\pi}\cos(m\theta)\bigg{[}1+\alpha^{2}-2\alpha\cos\theta+\beta_{c}^{2}\bigg{]}^{-s}d\theta.

\mathcal{B}_{s}^{(m),\mathrm{Tr}}=\frac{2}{\pi}\int\limits_{0}^{\pi}\cos(m\theta)\bigg{[}1+\alpha^{2}-2\alpha\cos\theta+\beta_{c}^{2}\bigg{]}^{-s}d\theta.

\mathcal{B}_{s}^{(m),\mathrm{T}}=\frac{2}{\pi}\int\limits_{0}^{\pi}\cos(m\theta)\bigg{[}1+\alpha^{2}-2\alpha\cos\theta+\beta^{2}\bigg{]}^{-s}d\theta.

\mathcal{B}_{s}^{(m),\mathrm{T}}=\frac{2}{\pi}\int\limits_{0}^{\pi}\cos(m\theta)\bigg{[}1+\alpha^{2}-2\alpha\cos\theta+\beta^{2}\bigg{]}^{-s}d\theta.

\mathcal{B}_{s}^{(m),\mathrm{H}}=\frac{2}{\pi}\int\limits_{0}^{\pi}\cos(m\theta)\bigg{[}1+\alpha^{2}-2\alpha\cos\theta+H^{2}(1+\alpha^{2})\bigg{]}^{-s}d\theta

\mathcal{B}_{s}^{(m),\mathrm{H}}=\frac{2}{\pi}\int\limits_{0}^{\pi}\cos(m\theta)\bigg{[}1+\alpha^{2}-2\alpha\cos\theta+H^{2}(1+\alpha^{2})\bigg{]}^{-s}d\theta

\mathcal{B}_{s}^{(m),\mathrm{TO}}=\frac{2}{\pi}\int\limits_{0}^{\pi}\cos(m\theta)\bigg{[}1+\alpha^{2}-2\alpha\cos\theta+S^{2}\alpha\bigg{]}^{-s}d\theta.

\mathcal{B}_{s}^{(m),\mathrm{TO}}=\frac{2}{\pi}\int\limits_{0}^{\pi}\cos(m\theta)\bigg{[}1+\alpha^{2}-2\alpha\cos\theta+S^{2}\alpha\bigg{]}^{-s}d\theta.

\displaystyle\Sigma_{d}(a)=\Sigma_{0}\bigg{(}\frac{a_{0}}{a}\bigg{)}^{p},~{}~{}~{}~{}~{}~{}e_{d}(a)=e_{0}\bigg{(}\frac{a_{0}}{a}\bigg{)}^{q}

\displaystyle\Sigma_{d}(a)=\Sigma_{0}\bigg{(}\frac{a_{0}}{a}\bigg{)}^{p},~{}~{}~{}~{}~{}~{}e_{d}(a)=e_{0}\bigg{(}\frac{a_{0}}{a}\bigg{)}^{q}

R_{d}=K\big{[}\psi_{1}e_{p}^{2}+\psi_{2}e_{p}e_{d}(a_{p})\cos(\varpi_{p}-\varpi_{d})\big{]},

R_{d}=K\big{[}\psi_{1}e_{p}^{2}+\psi_{2}e_{p}e_{d}(a_{p})\cos(\varpi_{p}-\varpi_{d})\big{]},

ψ_{1}

ψ_{1}

ψ_{2}

A_{d} (a_{p}) = \frac{2 K}{n _{p} a _{p}^{2}} ψ_{1}, B_{d} (a_{p}) = \frac{K}{n _{p} a _{p}^{2}} e_{d} (a_{p}) ψ_{2} .

A_{d} (a_{p}) = \frac{2 K}{n _{p} a _{p}^{2}} ψ_{1}, B_{d} (a_{p}) = \frac{K}{n _{p} a _{p}^{2}} e_{d} (a_{p}) ψ_{2} .

\psi_{1}^{\mathrm{Tr98}}=-\frac{1}{2\sqrt{\beta_{c}^{2}+1}}+\frac{\mathbf{E}\bigg{(}2/\sqrt{\beta_{c}^{2}+4}\bigg{)}}{\pi\sqrt{\beta_{c}^{2}+4}}=-\frac{1}{2}+\frac{1}{2\pi}+\mathcal{O}(\beta_{c}^{2})

\psi_{1}^{\mathrm{Tr98}}=-\frac{1}{2\sqrt{\beta_{c}^{2}+1}}+\frac{\mathbf{E}\bigg{(}2/\sqrt{\beta_{c}^{2}+4}\bigg{)}}{\pi\sqrt{\beta_{c}^{2}+4}}=-\frac{1}{2}+\frac{1}{2\pi}+\mathcal{O}(\beta_{c}^{2})

ψ_{1}^{TO16} \approx \frac{1}{2 S} + O (S), ψ_{2}^{TO16} \approx - \frac{1}{S} + O (S)

ψ_{1}^{TO16} \approx \frac{1}{2 S} + O (S), ψ_{2}^{TO16} \approx - \frac{1}{S} + O (S)

Σ_{d} (a) = \tilde{Σ}_{0} \frac{( a _{out} - a ) ^{2} ( a _{in} - a ) ^{2}}{( a _{out} - a _{in} ) ^{4}},

Σ_{d} (a) = \tilde{Σ}_{0} \frac{( a _{out} - a ) ^{2} ( a _{in} - a ) ^{2}}{( a _{out} - a _{in} ) ^{4}},

e_{d}(a)=\tilde{e}_{0}\bigg{(}1+\frac{a_{\rm out}-a}{a_{\rm out}-a_{\rm in}}\bigg{)}

e_{d}(a)=\tilde{e}_{0}\bigg{(}1+\frac{a_{\rm out}-a}{a_{\rm out}-a_{\rm in}}\bigg{)}

\Sigma_{d}(a)=\tilde{\Sigma}_{0}\exp\bigg{\{}\frac{4-[(a/a_{c})+(a_{c}/a)]^{2}}{w_{c}}\bigg{\}}

\Sigma_{d}(a)=\tilde{\Sigma}_{0}\exp\bigg{\{}\frac{4-[(a/a_{c})+(a_{c}/a)]^{2}}{w_{c}}\bigg{\}}

ϕ_{11} = ϕ_{22}

ϕ_{11} = ϕ_{22}

ϕ_{12}

δ A_{d} (a_{p})

δ A_{d} (a_{p})

δ B_{d} (a_{p})

δ ψ_{1}

δ ψ_{1}

δ ψ_{2}

ϕ_{11}

ϕ_{11}

\displaystyle=\frac{\alpha}{8}\bigg{[}\mathcal{B}_{3/2}^{(1),\mathrm{TO}}+\frac{1}{2}S^{2}\mathcal{B}_{3/2}^{(0),\mathrm{TO}}-\frac{3}{4}S^{2}(2+2\alpha^{2}+S^{2}\alpha)\mathcal{B}_{5/2}^{(0),\mathrm{TO}}\bigg{]},

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Potential softening and eccentricity dynamics in razor-thin, nearly-Keplerian discs

Antranik A. Sefilian1

and Roman R. Rafikov1,2

1 Department of Applied Mathematics and Theoretical Physics, CMS, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK

2 Institute of Advanced Study, Einstein Drive, Princeton, NJ 08540, USA E-mail: [email protected]

(Accepted XXX. Received YYY; in original form ZZZ)

Abstract

In many astrophysical problems involving discs (gaseous or particulate) orbiting a dominant central mass, gravitational potential of the disc plays an important dynamical role. Its impact on the motion of external objects, as well as on the dynamics of the disc itself, can usually be studied using secular approximation. This is often done using softened gravity to avoid singularities arising in calculation of the orbit-averaged potential — disturbing function — of a razor-thin disc using classical Laplace-Lagrange theory. We explore the performance of several softening formalisms proposed in the literature in reproducing the correct eccentricity dynamics in the disc potential. We identify softening models that, in the limit of zero softening, give results converging to the expected behavior exactly, approximately or not converging at all. We also develop a general framework for computing secular disturbing function given an arbitrary softening prescription for a rather general form of the interaction potential. Our results demonstrate that numerical treatments of the secular disc dynamics, representing the disc as a collection of $N$ gravitationally interacting annuli, are rather demanding: for a given value of the (dimensionless) softening parameter, $\varsigma\ll 1$ , accurate representation of eccentricity dynamics requires $N\sim C\varsigma^{-\chi}\gg 1$ , with $C\sim O(10)$ , $1.5\lesssim\chi\lesssim 2$ . In discs with sharp edges a very small value of the softening parameter $\varsigma$ ( $\lesssim 10^{-3}$ ) is required to correctly reproduce eccentricity dynamics near the disc boundaries; this finding is relevant for modelling planetary rings.

keywords:

celestial mechanics — methods: analytical — planet-disc interactions — planets and satellites: rings

††pubyear: 2019††pagerange: Potential softening and eccentricity dynamics in razor-thin, nearly-Keplerian discs –D

1 Introduction

Astrophysical discs orbiting a central mass $M_{c}$ are ubiquitous in a variety of contexts – galactic, stellar, and planetary (Latter et al., 2017). In many instances, masses of such discs $M_{d}$ are much less than the central object mass. Despite this fact, gravity of such discs can still play an important dynamical role in the orbital evolution of their constituent particles as well as the dynamics of external objects (e.g. Goldreich & Tremaine, 1979; Heppenheimer, 1980; Ward, 1981; Kocsis & Tremaine, 2011; Kazandjian & Touma, 2013; Teyssandier et al., 2013; Meschiari, 2014; Silsbee & Rafikov, 2015; Petrovich et al., 2019; Sefilian & Touma, 2019). Consequently, characterizing dynamical effects of disc gravity is important.

Whenever $M_{d}\ll M_{c}$ , particles perturbed by the disc gravity move on nearly-Keplerian orbits which evolve rather slowly. This justifies the use of the so-called secular approximation which implies averaging of the fast-evolving dynamical variables over the orbits of particles under consideration (Murray & Dermott, 1999). The orbit-averaging procedure, also known as Gauss’ method, is equivalent to calculating the time-averaged potential due to orbiting point masses by smearing them into massive elliptical "wires" (having shape of their eccentric orbits) with non-uniform linear density proportional to the time spent by an object at a particular phase of its orbit. Such orbit-averaged potential, also known as secular disturbing function $R_{d}$ , fully determines the secular dynamics of the system.

For a test particle with semi-major axis $a_{p}$ , eccentricity $e_{p}$ , and apsidal angle $\varpi_{p}$ due to a co-planar point mass $\delta m_{d}$ orbiting with semi-major axis $a$ , eccentricity $e_{d}$ , and apsidal angle $\varpi_{d}$ , upon smearing into elliptical rings, the secular disturbing function takes the form (Murray & Dermott, 1999)

[TABLE]

valid for $a>a_{p}$ as well as $a<a_{p}$ , as long as particle orbits do not cross. Here $b_{s}^{(m)}(\alpha)$ is the Laplace coefficient defined by

[TABLE]

which obeys $b_{s}^{(m)}(\alpha^{-1})=\alpha^{2s}b_{s}^{(m)}(\alpha)$ . Explicit time independence of $\delta R$ guarantees that the semi-major axes of the secularly interacting objects stay fixed.

When considering gravitational effects of a razor-thin continuous disc with smooth distribution of surface density, a straightforward way to compute the secular disturbing function would be to orbit-average the disc potential (obtained by direct integration over its full surface) along the particle orbit. However, this procedure involves a triple integration (two-dimensional integral over the disc surface and orbit averaging) and is numerically challenging.

A more efficient approach lies in representing the disc as a collection of massive, nested, confocal elliptical "wires" (also referred to as "annuli" or "rings" in this work) with fixed semi-major axes (e.g. Touma et al., 2009; Batygin, 2012). Due to the additive nature of gravity, the disturbing function due to a disc can be represented as a sum of individual contributions in the form (1) produced by all wires, which amounts to integration of $\delta R$ (Eq. 1) over the radial extent of the disc:

[TABLE]

where $a_{\rm in}$ and $a_{\rm out}$ are the semi-major axes of the inner and outer disc edges. In this case, provided that $\delta R$ is known as a function of $a$ , only a single integration (over the semi-major axes of the rings) is needed, significantly accelerating calculations111The Laplace coefficients entering in $\delta R$ can be easily evaluated, without relying on integration over $\theta$ in Eq. (2), by expressing them through elliptic integrals, see Appendix C.3..

Unfortunately, this straightforward procedure is ill-posed from the mathematical point of view. Indeed, it is well known that the Laplace coefficients $b_{3/2}^{(m)}$ featured in Eq. (1) diverge as $b_{3/2}^{(m)}(\alpha)\rightarrow(1-\alpha)^{-2}$ when $\alpha\to 1$ . This implies that the radial integration in Eq. (3) encounters an essential singularity at $a=a_{p}$ . As a result, for a co-planar particle orbiting inside a razor-thin disc, $a_{\rm in}\leq a_{p}\leq a_{\rm out}$ , this direct way of computing $R_{d}$ does not converge to a finite value.

This divergence, as well as the pressing need for having an efficient way of computing $R_{d}$ (via a one-dimensional integration over $a$ only), have motivated the development of alternative analytic approaches for calculating $R_{d}$ . These approaches can be generally grouped into two classes. Calculations of one kind are rooted in the derivation of the potential of an axisymmetric disc with power law surface density profile presented in Heppenheimer (1980), which does not suffer from the singularity of Laplace-Lagrange secular theory. A number of subsequent studies used this approach (Ward, 1981) and extended it to the case of eccentric discs, both apsidally aligned (Silsbee & Rafikov, 2015; Davydenkova & Rafikov, 2018) and misaligned (Davydenkova & Rafikov, in prep.). Higher order (in eccentricity) extensions of this approach have also been developed (Sefilian & Touma, 2019). This framework for treating secular dynamics has been extensively verified using direct orbit integrations under different conditions (Silsbee & Rafikov, 2015; Fontana & Marzari, 2016; Davydenkova & Rafikov, 2018). In this work, we refer to this type of calculation as the unsoftened Heppenheimer’s method.

Unfortunately, by construction Heppenheimer’s method is inapplicable in situations where the disc eccentricity rapidly varies with semi-major axis, potentially resulting in orbit crossings (Davydenkova & Rafikov, 2018). An alternative approach, which avoids this problem, while at the same time alleviating the aforementioned singularity, is to use softened gravity by spatially smoothing the Newtonian point-mass potential in various ways – both analytically (e.g. Tremaine, 1998, 2001; Touma, 2002; Hahn, 2003; Touma & Sridhar, 2012; Teyssandier & Ogilvie, 2016) and numerically (e.g. Touma et al., 2009). In these models, the classical Laplace-Lagrange disturbing function (Eq. 1) is modified by softening the interaction potential in some way to circumvent the divergence of $R_{d}$ as $a\to a_{p}$ . In this method orbit crossing does not lead to problems as long as the softening scale is finite. However, a physical justification for a specific form of softening (absent in the Heppenheimer (1980) approach) often remains unclear, making the introduction of softening rather arbitrary.

The primary goal of our present work is to assess how well the different calculations relying on potential softening reproduce secular dynamics driven by the gravity of a razor-thin disc. The main metric we use in this exercise is the convergence of the results of such calculations to the true secular evolution (represented by the un-softened Heppenheimer method) in the limit of vanishing softening, when the limit of Newtonian gravity is recovered. Complementary to this, we develop a general framework for computing the well-behaved secular disturbing function for a broad range of softened gravitational potentials.

Our work is organized as follows. We describe the general analytical expressions governing the orbit-averaged potential due to a coplanar disc of arbitrary structure and arbitrary softening prescription in §2. Having provided a brief account of the different softened potentials under our probe and the un-softened approach of Heppenheimer in §2.1 and §2.2, respectively, we analyze their performance in reproducing the correct secular dynamics for various disc models in §3, §4 and §5. We discuss and briefly summarize our results in §6 and §7 respectively. Technical details of our calculations can be found in Appendices.

2 Disturbing function due to a disk

Prior to providing the details of different softening prescriptions examined in this work in §2.1, we briefly summarize some of their common features. The ultimate goal of all these prescriptions is the calculation of the disturbing function $R_{d}$ due to gravity of a (generally eccentric) disc comprised of massive objects (stars, planetesimals, ring particles) or fluid elements (in gaseous discs) moving on Keplerian orbits.

We consider the disc to be razor-thin and coplanar. Mass distribution of such a disc can be uniquely characterized by the mass density per unit semi-major axis $\mu_{d}(a)$ , eccentricity $e_{d}(a)$ , and apsidal angle $\varpi_{d}(a)$ of the trajectories of its constituent elements, as functions of the semi-major axis $a$ . In practice, it is often convenient to use the surface density at periastron $\Sigma_{d}(a)$ instead of $\mu_{d}(a)$ ; its relation to $\mu_{d}$ for arbitrary profiles of $e_{d}$ and $\varpi_{d}$ has been established in Statler (2001), Davydenkova & Rafikov (2018) and Davydenkova & Rafikov (in prep.). Constancy of semi-major axis in secular theory implies that $\mu_{d}(a)$ does not change in time. The same statement is true for $\Sigma_{d}(a)$ to lowest order in $e_{d}$ since $\mu_{d}(a)\approx 2\pi a\Sigma_{d}(a)+O(e_{d})$ (Davydenkova & Rafikov, 2018).

Close inspection of the various softening methods for computing secular disc potential (§2.1) reveals that all of them arrive at the following general form of the disturbing function for a test particle moving on an orbit with the semi-major axis $a_{p}$ , eccentricity $e_{p}$ , and apsidal angle $\varpi_{p}$ :

[TABLE]

Here $n_{p}$ is the test-particle mean motion ( $n_{p}^{2}=GM_{c}/a_{p}^{3}$ ), and we have introduced a two-component eccentricity vector for a test particle $\mathbf{e}_{p}=e_{p}(\cos\varpi_{p},\sin\varpi_{p})$ .

The coefficients $\mathsf{A}_{d}$ and $\bm{\mathsf{B}}_{d}$ in Eq. (4) are related to the disc mass (or surface density) and eccentricity profiles in the following fashion:

[TABLE]

where $\mathbf{e}_{d}=e_{d}(a)(\cos\varpi_{d}(a),\sin\varpi_{d}(a))$ is the eccentricity vector for an annular disc element222We refer the reader to Heppenheimer (1980); Silsbee & Rafikov (2015); Davydenkova & Rafikov (2018) for the expressions of $A_{d}$ and $B_{d}$ computed using the un-softened Heppenheimer method for different disc models..

Functions $\phi_{ij}(\alpha)$ , $i,j=1,2$ entering these expressions fully characterize the softened ring-ring secular interaction, see Eq. (11). They are unique for each potential softening prescription, with explicit forms for the models that we explore in this work specified in Table 1. This Table shows that coefficients $\phi_{ij}$ appearing in the literature are linear combinations of softened Laplace coefficients $\mathcal{B}_{s}^{(m)}$ defined by

[TABLE]

The softening parameter $\epsilon(\alpha)$ appearing in this definition remains non-zero as $\alpha\to 1$ , thus preventing the divergence of the softened Laplace coefficients $\mathcal{B}_{s}^{(m)}(\alpha,\epsilon)$ at $\alpha=1$ (unlike the classical $b_{s}^{(m)}(\alpha)$ ). The explicit form of $\epsilon(\alpha)$ is different for every softening method considered in this work, see §2.1 and Table 1. Appendix C collates some useful relations for softened Laplace coefficients $\mathcal{B}_{s}^{(m)}(\alpha,\epsilon)$ , as well as their approximate asymptotic behavior and relationships to complete elliptic integrals.

The mathematical structure of $R_{d}$ given by Eq. (4) is similar to that of the classical Laplace-Lagrange planetary theory (Murray & Dermott, 1999), see Eq. (1). Indeed, let us consider mass distribution of a point mass smeared along an elliptical orbit, $\mu_{d}(a)\rightarrow m_{\rm pl}\delta(a-a_{\rm pl})$ (where $\delta(z)$ is the Dirac delta-function), and set softening to zero (so that $\mathcal{B}_{s}^{(m)}(\alpha,\epsilon\to 0)\to b_{s}^{(m)}(\alpha)$ ). Then one finds that $R_{d}$ reduces to the un-softened, orbit-averaged potential $\delta R$ due to a planet with mass $m_{\rm pl}$ and semi-major axis $a_{\rm pl}$ , with the unsoftened coefficients $\phi_{ij}$ in the form (Murray & Dermott, 1999)

[TABLE]

see Eq. (1).

Accordingly, it is intuitive to think of Eqs. (4)-(6) as the continuous version of classical Laplace-Lagrange planetary theory, modified by the introduction of non-zero softening parameter $\epsilon$ to avoid the mathematical divergence of the classical disturbing function as $a\to a_{p}$ .

We emphasize that the functional forms of $\phi_{ij}$ are not simple replacements of $b_{s}^{(m)}$ appearing in the unsoftened definition (8) - (9) by $\mathcal{B}_{s}^{(m)}$ . This can be seen in Table 1 where we summarize some of the expressions for $\phi_{ij}(\alpha)$ proposed in the literature and analyzed in this paper (see §2.1). Nevertheless, examination of these expressions shows that when $\epsilon^{2}(\alpha)\to 0$ , the coefficients $\phi_{ij}(\alpha)$ do reduce to their unsoftened versions $\phi_{ij}^{\rm LL}(\alpha)$ given by Eqs. (8) - (9).

In Appendix A we show that the form of the disturbing function given by Eqs. (4)-(6) is generic for a wide class of softening models (and not just the ones covered in §2.1), for which the interaction potential between the two masses $m_{1}$ and $m_{2}$ ( $m_{i}\ll M_{c}$ ) located at ${\bf r}_{1}$ and ${\bf r}_{2}$ , correspondingly, relative to the central mass, has a form333Note that the inter-particle force resulting from such potential does not, in general, obey Newton’s third law (as long as $\mathcal{F}(r_{1},r_{2})\neq$ const).

[TABLE]

with $i,j=1,2$ and $j\neq i$ . Here $\mathcal{F}(r_{1},r_{2})$ represents an arbitrary softening function introduced to cushion the singularity which arises otherwise at null inter-particle separations. Note that in general this potential may depend not only on the relative distance between the two masses $\mathbf{r}_{1}-\mathbf{r}_{2}$ , but also on their distances to the dominant central mass $r_{1}$ , $r_{2}$ .

Explicit demonstration of the connection between the potential (10) and $R_{d}$ given by Eq. (4) represents a stand-alone result of this work. In particular, our calculations in Appendix A, which can be skipped at first reading, show that the softening parameter $\epsilon$ featured in the definition (7) is related to $\mathcal{F}$ via $\epsilon^{2}=[{\max(a_{1},a_{2})}]^{-2}\mathcal{F}(a_{1},a_{2})$ , where $a_{1,2}$ are the semi-major axes of the interacting particles (see Eq. 55). The most general expressions of $\phi_{ij}$ entering the arbitrarily softened ring-ring disturbing function,

[TABLE]

(here $i={1,2}$ and $j\neq i$ ) is given by Eqs. (56)-(58) in terms of $\mathcal{B}_{s}^{(m)}(\alpha,\mathcal{F})$ . In the above expression, we have defined $a_{>}={\rm max}(a_{1},a_{2})$ and $a_{<}={\rm min}(a_{1},a_{2})$ such that 444Here we clarify that the definitions of $\phi_{11}(\alpha)$ and $\phi_{22}(\alpha)$ , even when different (see Table 1 and Appendix A), are swapped upon interchanging $a_{1}$ with $a_{2}$ but keeping, by construction, $\alpha=a_{<}/a_{>}<1$ – see Eqs. (56), (57) for details. $\alpha=a_{<}/a_{>}$ .

Note that in equations (5) and (6) we split integration over $a$ in two parts: over the part of the disc interior to $a_{p}$ , and exterior to it. We do this because for some softening functions $\mathcal{F}$ the coefficients $\phi_{ij}(\alpha)$ do not obey certain symmetry properties when $a/a_{p}$ is replaced with $a_{p}/a$ , see Eq. (69). Moreover, in general $\phi_{11}$ and $\phi_{22}$ are not necessarily identical as in classical Laplace-Lagrange theory (i.e. Eq. 8); see Table 1 and Appendix A for further details.

As to the physical meaning of $\mathsf{A}_{d}$ and $\bm{\mathsf{B}}_{d}$ , we remind the reader that $\mathsf{A}_{d}$ represents the precession rate of the free eccentricity vector of a test particle in the disc potential, while $\bm{\mathsf{B}}_{d}$ characterizes the torque exerted on the particle orbit by the non-axisymmetric component of the disc gravity. Corresponding forced eccentricity vector is $\mathbf{e}_{p,f}=-\bm{\mathsf{B}}_{d}/A_{d}$ . In particular, test-particles initiated on circular orbits experience eccentricity oscillations of maximum amplitude $e_{p}^{m}=2\left|\mathbf{e}_{p,f}\right|$ .

As $\mathsf{A}_{d}(a_{p})$ and $\mathsf{B}_{d}(a_{p})$ uniquely determine $R_{d}$ for different forms of softening, comparison of their behavior in the limit of $\epsilon\to 0$ with that found in the unsoftened Heppenheimer (1980) approach (validated in Silsbee & Rafikov 2015; Fontana & Marzari 2016; Davydenkova & Rafikov 2018) is sufficient to assess the validity of a particular softening model, see §3.

2.1 Summary of existing softening models

Here we provide a brief description of the four different softening prescriptions that have been previously proposed in the literature. Corresponding expressions for their softening parameters $\epsilon^{2}(\alpha)$ and coefficients $\phi_{ij}(\alpha)$ are provided in Table 1.

2.1.1 Formalism of Tremaine (1998) – Tr98

Tremaine (1998) suggested an expression for the secular disturbing function due to a continuous disc, which uses modified Laplace coefficients in the form

[TABLE]

Here $\beta_{c}^{2}$ is the dimensionless softening parameter, treated as a constant, i.e. independent of distance. The physical interpretation of this manoeuvre is that $\beta_{c}$ , inhibiting the formal divergence of $R_{d}$ as $a\to a_{p}$ , can be viewed as the disc aspect ratio. Within this prescription, it is intuitive to think of the eccentric "wires" that comprise the disc as having a distance-dependent radius $b=\beta_{c}\max(a_{1},a_{2})$ . In Tremaine (1998) coefficients $\phi_{ij}(\alpha)$ were expressed as derivatives of $\mathcal{B}_{1/2}^{(m),\mathrm{Tr}}$ with respect to $\alpha$ , see equations (26) of Tremaine (1998). These expressions, along with their versions modified using the recursive relations for Laplace coefficients (see Appendix C.1), can be found in Table 1.

2.1.2 Formalism of Touma (2002) – T02

Touma (2002) derived the orbit-averaged potential of a disc by assuming individual particles comprising the disc to interact via Plummer potential with a fixed length scale $b_{c}$ (Binney & Tremaine, 2008). Smearing particles into gravitating eccentric wires, Touma (2002) (see also Touma & Sridhar, 2012) derived the expressions (equations (6) of Touma (2002)) for $\phi_{ij}(\alpha)$ in the form of linear combinations of softened Laplace coefficients $\mathcal{B}_{s}^{(m),\mathrm{T}}$ , similar to those of Tremaine (1998):

[TABLE]

However, in Touma (2002) the softening parameter $\epsilon^{2}(\alpha)=\beta^{2}$ is no longer a constant but depends on the distance such that $\beta=b_{c}/\max(a_{1},a_{2})$ . Within this formalism, one can think of a disc as comprised of nested annuli with a constant thickness $b_{c}$ .

2.1.3 Formalism of Hahn (2003) – H03

Hahn (2003) computed the orbit-averaged interaction between two eccentric wires by accounting for their vertical thickness. The vertical extent $h$ of a ring effectively softens its gravitational potential over a dimensionless scale $H\sim h/a$ , which was assumed to be constant in that work (see also Ward, 1989). Hahn (2003) demonstrated that the resultant $\phi_{ij}(\alpha)$ are functions of softened Laplace coefficients

[TABLE]

with constant $H\ll 1$ . In other words, the softening parameter is given by $\epsilon^{2}(\alpha)=H^{2}(1+\alpha^{2})$ in that work. The explicit expressions for $\phi_{ij}(\alpha)$ in terms of $\mathcal{B}_{s}^{(m),\mathrm{H}}$ are given by equations (17) of Hahn (2003).

2.1.4 Formalism of Teyssandier &

Ogilvie (2016) – TO16

Teyssandier & Ogilvie (2016) modified the unsoftened expressions (8), (9) for $\phi_{ij}^{\rm LL}(\alpha)$ by simply replacing the usual Laplace coefficients $b_{s}^{(m)}$ with softened versions defined such that

[TABLE]

Thus, their softening parameter is $\epsilon^{2}(\alpha)=S^{2}\alpha$ , where $S$ is a dimensionless constant. According to the authors, this substitution approximates the process of vertical averaging over the disc with constant aspect ratio $S$ , and alleviates the classical singularity. The corresponding expressions for $\phi_{ij}(\alpha)$ are given by equations (7)-(9) of Teyssandier & Ogilvie (2016).

The aforementioned softening prescriptions have their softening parameters $\epsilon^{2}(\alpha)$ controlled by different constants — $\beta_{c},b_{c},H,$ and $S$ . For this reason, in what follows – with some abuse of notation – we will collectively refer to these constants as “softening parameters" and denote them by $\varsigma$ .

2.2 The unsoftened Heppenheimer method

A different approach to computing the disturbing function of a razor-thin disc has been developed by Heppenheimer (1980) without resorting to any form of softened gravity (see also Ward, 1981). The essence of this method is in computing the potential by direct integration over the disc surface before expanding the integral limits (which involve instantaneous particle position $r$ ) in terms of small eccentricity of a test particle555Note that the order of these procedures is opposite to what is usual in the Laplace-Lagrange treatment (e.g. Murray & Dermott, 1999). For further details, see e.g. Heppenheimer (1980).. This expansion is followed by time-averaging over the orbit of a test particle.

The outcome of this procedure is a set of expressions, akin to Eq. (4)-(6), which are convergent throughout the disc, in contrast to the classical Laplace-Lagrange theory. Mathematically, this convergent behavior is due to the fact that the emergent expressions contain Laplace coefficients $b_{1/2}^{(m)}(\alpha)$ – and not $b_{3/2}^{(m)}$ – which diverge only weakly (logarithmically) as $\alpha\to 1$ : $b_{1/2}^{(m)}(\alpha)\propto\log(1-\alpha)$ . As a result, upon integrating these expressions over the radial extent of the disc, one obtains a convergent and finite result for $R_{d}$ . Physically, convergent expression is only natural since the calculation of the disk potential by direct two-dimensional integration over its surface is fully convergent at every point in the disc. The Heppenheimer’s method simply allows one to properly capture this property, unlike the standard Laplace-Lagrange procedure (when applied to continuous discs).

In his pioneering calculation, Heppenheimer (1980) applied this method to axisymmetric power-law discs to recover the orbit-averaged disc potential to second order in eccentricities. This calculation has been subsequently extended to more general disc structures (Silsbee & Rafikov, 2015; Davydenkova & Rafikov, 2018) (hereafter, SR15 and DR18 respectively), as well as to higher order in eccentricities (Sefilian & Touma, 2019). This framework has been extensively verified for eccentric discs using direct integrations of test particle orbits in actual disc potentials (e.g. SR15, Fontana & Marzari, 2016, DR18), validating this approach.

3 Comparison: Power-Law Discs

Our goal is to examine the performance of different softening prescriptions outlined in §2.1 in comparison with the results obtained using the un-softened Heppenheimer method (§2.2).

We start this exercise using a model of apse-aligned (i.e. $d\varpi_{d}/da=0$ ), truncated power-law (hereafter PL) disc as a simple example. We characterize surface density and eccentricity of such a disc by

[TABLE]

for $a_{\rm in}\leq a\leq a_{\rm out}$ , where $\Sigma_{0}$ and $e_{0}$ are the pericentric surface density and eccentricity of the disc at some reference semi-major axis $a_{0}$ .

Plugging this anzatz into Eqs. (4) – (6), the secular disturbing function $R_{d}$ due to PL discs can be simplified to (Silsbee & Rafikov, 2015)

[TABLE]

where $K=\pi G\Sigma_{0}a_{0}^{p}a_{p}^{1-p}$ and the dimensionless coefficients $\psi_{1}$ and $\psi_{2}$ are given by

[TABLE]

with $\alpha_{1}=a_{\rm in}/a_{p}$ and $\alpha_{2}=a_{p}/a_{\rm out}$ .

The coefficients $\psi_{1}$ and $\psi_{2}$ are functions of the power-law indices ( $p$ and $q$ ), any softening parameter involved (through $\phi_{ij}$ ), as well as the test-particle semi-major axis $a_{p}$ (through $\alpha_{1,2}$ ). They are related to $\mathsf{A}_{d}$ and $\mathsf{B}_{d}$ via

[TABLE]

As shown in Appendix D, for certain ranges of power-law indices $p$ and $q$ both $\psi_{1}$ and $\psi_{2}$ converge to values depending only on $p$ and $q$ and a softening parameter used, provided that the test-particle orbit is well-separated from the disc boundaries (i.e. in the limit $\alpha_{1,2}\to 0$ ). For $p$ and $q$ in these ranges (determined in Appendix D for each of the considered softened formalisms, similar to SR15), the coefficients $\psi_{1}$ and $\psi_{2}$ are determined by the local behavior of $\Sigma_{d}(a)$ and $e_{d}(a)$ in the vicinity of test-particle semi-major axis.

Given this, we first focus on infinitely extended ( $\alpha_{1,2}\to 0$ ) PL discs with $p$ and $q$ within these ranges (we defer discussion of secular dynamics near the disc edges to §5). Then, $\psi_{1}$ and $\psi_{2}$ become independent of $a_{p}$ (i.e. functions of $p$ , $q$ , and $\varsigma$ only), making them useful as simple metrics for judging the validity of different models of softening.

3.1 Behavior with respect to variation of softening

Figure 1 illustrates the behavior of $\psi_{1}$ and $\psi_{2}$ predicted by each of the softening formalisms described in §2.1 for an infinite PL disc, shown as a function of the corresponding “softening"666The softening length $b_{c}$ present in the formulation of Touma (2002) is scaled by the test-particle semi-major axis $a_{p}$ in all the figures where we present results for infinite PL discs. We do this to properly collate the results computed by different softening formalisms in one figure. $\varsigma$ for two different sets of $p,q$ (indicated in panel B). For reference, black horizontal lines show the values of $\psi_{1}$ and $\psi_{2}$ expected from the calculations of SR15 using the un-softened Heppenheimer approach777Equations (A37) and (A38) in Silsbee & Rafikov (2015) provide analytic expressions for $\psi_{1}$ and $\psi_{2}$ , respectively, for infinite PL discs..

The left panels of Figure 1 illustrate the behavior of the softening models of Tremaine (1998), Touma (2002) and Hahn (2003). They demonstrate that the latter two formalisms predict $\psi_{1}$ and $\psi_{2}$ in quantitative agreement with the unsoftened calculations of SR15: results of both Touma (2002) (blue) and Hahn (2003) (red) converge to the SR15 results as their corresponding softening $\varsigma$ approaches zero; both the amplitude and sign of $\psi_{1}$ and $\psi_{2}$ are reproduced. It is also evident that, depending on disc model, $\psi_{1}$ and $\psi_{2}$ converge to values given by SR15 at different values of softening. Nevertheless, we generally888For particles with orbits near sharp disc edges, we find that smaller values of $\varsigma$ is required to recover the expected dynamics, see §5. find that $\varsigma\lesssim 10^{-3}$ guarantees the convergence of $\psi_{1}$ and $\psi_{2}$ to within few per cent of the correct values for all $p$ and $q$ as long as $a_{\rm in}\ll a_{p}\ll a_{\rm out}$ (see Figure 4).

The same panels also indicate that $\psi_{1}(\varsigma)$ and $\psi_{2}(\varsigma)$ predicted by the softened formalism of Tremaine (1998) (green), while converging to finite values as $\varsigma=\beta_{c}\to 0$ , do not reproduce the SR15 results exactly in this limit. Indeed, one can see that even for the smallest adopted value of $\beta_{c}=10^{-3}$ , the softening prescription of Tremaine (1998) yields $\psi_{1}$ and $\psi_{2}$ different by tens of per cent from SR15. It is easy to demonstrate that these quantitative differences do not vanish by further decreasing $\beta_{c}$ . For instance, when $p=1$ , the coefficient $\psi_{1}$ can be evaluated analytically as

[TABLE]

in agreement with Panel A ( $\mathbf{E}(k)$ is the complete elliptic integral of a second kind). At the same time, the unsoftened approach of SR15 predicts $\psi_{1}=-1/2$ for $p=1$ disc. Moreover, close inspection of Fig. 1A,B shows that, in the limit of $\beta_{c}\to 0$ , the $\psi_{1}$ and $\psi_{2}$ curves computed using softening model of Tremaine (1998) are offset vertically from the unsoftened calculations by $1/2\pi$ and $-1/\pi$ , respectively, for any $(p,q)$ – see also Fig. 4. We will analyze reasons for this quantitative discrepancy in §6.1.

Right panels of Fig. 1 show the behavior of $\psi_{1}$ (Panel C) and $\psi_{2}$ (Panel D) as a function of “softening", $\varsigma=S$ , resulting from the approach of Teyssandier & Ogilvie (2016). There are several features to note here. First, this model predicts $\psi_{1}>0$ for all values of softening $S$ and disc models (i.e. $p$ and $q$ ), implying prograde free precession. This is in contrast with the other softening prescriptions, as well as SR15, which correctly capture retrograde free precession for $p=1$ and prograde for $p=-0.5$ (see Panel A). Similarly, $\psi_{2}$ is always negative, contrary to the expectations (see Panel B). Second, in the limit of $S\to 0$ , both $\psi_{1}$ and $\psi_{2}$ attain values independent of the disc model, which is clearly inconsistent with the dependence on $(p,q)$ seen in Figure 1A, B. Third, and most importantly, both $\psi_{1}$ and $\psi_{2}$ diverge as the softening $S\to 0$ . Indeed, it suffices to employ the asymptotic expansion of the Laplace coefficients $\mathcal{B}_{3/2}^{(m),\mathrm{TO}}$ in the limit of $\alpha\rightarrow 1$ (Eq. 72) to demonstrate that both $\psi_{1}$ and $\psi_{2}$ (Eqs. 18 - 19) behave as

[TABLE]

as $S\to 0$ for all values of $p$ and $q$ . The behavior shown in Fig. 1C, D agrees with these asymptotic expressions.

3.2 Details of convergence of different softening prescriptions

Different softening prescriptions explored in this work are designed to modify the behavior of the integrand in equations (5)-(6) primarily in the vicinity of the test particle orbit, i.e. as $a\to a_{p}$ or $\alpha\to 1$ . For this reason, it is interesting to look in more detail on how this modification actually allows each softening model to achieve (or not) the expected results. This exercise also illustrates the contribution of different parts of the disc to secular dynamics.

To this goal we compute the values of $\psi_{1}$ and $\psi_{2}$ in an infinitely extended PL disc, like in §3.1, but now with a narrow clean gap (in semi-major axis) just around the test particle orbit, and explore the effect of varying the width of this gap (Ward, 1981). The inner and outer edges of the gap, in which $\Sigma_{d}(a)$ is set to zero, are at $a_{d,i}=(1-x)a_{p}\leq a_{p}$ and $a_{d,o}=(1-x)^{-1}a_{p}\geq a_{p}$ , respectively, with a single parameter $x$ controlling the gap width. As $x\to 0$ , the width of the gap goes to zero. We compute secular coefficients in such a gapped disc denoted $\tilde{\psi}_{1}(x)$ and $\tilde{\psi}_{2}(x)$ , by appropriately changing the upper integration limits in the definitions (18)-(19), i.e. from 1 to $\alpha_{m}\equiv 1-x$ . This eliminates gravitational effect of the disc annuli with $a_{d,i}(x)<a<a_{d,o}(x)$ .

In Figure 2 we display the behavior of $\tilde{\psi}_{1}(x)$ (Panel A) and $\tilde{\psi}_{2}(x)$ (Panel B) as a function of $x=1-\sqrt{a_{d,i}/a_{d,o}}$ for various values of softening $\varsigma$ to highlight the effects of different softening prescriptions. The calculations assume a base PL disc model with $p=1$ and $q=0.5$ (recall that $\psi_{1}$ depends on $p$ , while $\psi_{2}$ depends on $p+q$ ; Eqs. 18, 19). There are several notable features in this figure.

First, when the gap is wider than the characteristic softening length $\varsigma a_{p}$ , i.e. $\varsigma\lesssim x\leq 1$ , the amplitudes of both $\tilde{\psi}_{1}(x)$ and $\tilde{\psi}_{2}(x)$ increase from zero at $x=1$ (infinitely wide gap) to their maximum values reached at $x\sim\varsigma$ . In all cases $\psi_{1}$ is positive, meaning prograde precession of a test particle orbit in a wide gap, in agreement with the unsoftened results of Ward (1981) and Davydenkova & Rafikov (2018) — secular effect of a collection of distant disc "wires" conforms to expectations of the classical Laplace-Largange theory (i.e. prograde precession).

In the range $\varsigma\lesssim x\ll 1$ we find that $\tilde{\psi}_{1}(x)\sim|\tilde{\psi}_{2}(x)|\sim x^{-1}$ , irrespective of the softening model used; their maximum values are always $\sim\varsigma^{-1}$ . This convergent behavior is easy to understand since for $\varsigma\lesssim x$ the role of softening is negligible, $\mathcal{B}_{s}^{(m)}(\alpha,\varsigma)\approx b_{s}^{(m)}(\alpha)$ , and all $\phi_{ij}$ effectively reduce to their classical counterparts $\phi_{ij}^{\rm LL}$ given by Eqs. (8) - (9), which can be easily verified using the expressions listed in Table 1. The scaling of $\tilde{\psi}_{1}(x)$ and $|\tilde{\psi}_{2}(x)|$ with $x$ is simply a result of asymptotic behavior of $b_{3/2}^{(m)}(\alpha)\rightarrow(1-\alpha)^{-2}$ as $\alpha\to 1$ , upon radial integration in Eqs. (18) – (19).

Second, upon reaching their extrema at $x\sim\varsigma$ , amplitudes of $\tilde{\psi}_{1}(x)$ and $\tilde{\psi}_{2}(x)$ computed using softening prescriptions of Tr98, T02 and H03 start decreasing as $x$ decreases. In the range of semi-major axes corresponding to $x\lesssim\varsigma$ , softening significantly modifies the behavior of $\mathcal{B}_{s}^{(m)}(\alpha,\varsigma)$ away from the divergent behavior of $b_{s}^{(m)}(\alpha)$ . The modification is such that the softened interaction with the disc annuli $\lesssim\varsigma a_{p}$ away from the test-particle orbit starts to dynamically counteract the contribution of the more distant annuli (with $x\approx 1$ ). As a result of this compensation, $\tilde{\psi}_{1}$ and $\tilde{\psi}_{2}$ cross zero and change sign at some $x=C\varsigma^{2}$ , where $C\sim 1$ is a constant999For $p=1$ , $\tilde{\psi}_{1}$ becomes analytic for the softened formalisms of both H03 and Tr98 allowing us to quantify the value of $C$ . Performing the integral over $d\alpha$ in Eq. (18) - (19), we find $C_{Tr98}=(\pi-1)/2$ and $C_{H03}=\pi$ ; in agreement with Fig. 2. For other values of $p$ and $q$ , for which $\psi_{1}<0$ (c.f. Fig. 4), we numerically find that $C$ varies by at most a factor of ten..

At the same time, $\tilde{\psi}_{1}^{\mathrm{TO16}}$ and $\tilde{\psi}_{2}^{\mathrm{TO16}}$ calculated according to Teyssandier & Ogilvie (2016) clearly show different behavior. Instead of decreasing in amplitude as $x\lesssim\varsigma$ , they remain essentially constant, having reached their saturated values $\sim\varsigma^{-1}$ at $x\sim\varsigma$ . This explains the lack of convergence with $S$ obvious in Figure 1C, D, since the values to which $|\tilde{\psi}_{1}^{\mathrm{TO16}}|$ and $|\tilde{\psi}_{2}^{\mathrm{TO16}}|$ converge keeps increasing as $\varsigma\to 0$ . Moreover, both coefficients also never change sign, always predicting prograde precession ( $\tilde{\psi}_{1}^{\mathrm{TO16}}>0$ ). The origin of this difference with other smoothing prescriptions will be addressed in §6.2.

Upon further decrease of $x$ below $\varsigma^{2}$ , both $\tilde{\psi}_{1}$ and $\tilde{\psi}_{2}$ computed using models of Tr98, T02 and H03 ultimately converge to their corresponding values obtained for a continuous disc (i.e. for $x=0$ , see Fig. 1) independent of the assumed value of $\varsigma$ .

We note that the opposite contributions to e.g. $\psi_{1}$ produced by the distant ( $x\gtrsim\varsigma$ , positive) and nearby (i.e. with $x\lesssim\varsigma$ , negative) disc annuli is not unique to softened gravity. Indeed, both Ward (1981) and Davydenkova & Rafikov (2018), using the un-softened Heppenheimer method, found that a particle orbit fully embedded in a $p=1$ disc has negative precession rate, whereas a particle orbiting fully in the gap precesses in the positive sense (and at high rate if the gap is narrow). As the gap width is reduced, a smooth transition between the two regimes must occur as the test-particle orbit starts crossing the gap edge (i.e. for $x\lesssim e_{p}$ ), with the disc annuli crossing the particle orbit giving rise to a negative contribution to $\tilde{\psi}_{1}$ . Eventually, the shrinking of the gap brings $\tilde{\psi}_{1}$ to a finite negative value (for $p=1$ disc) as $x\to 0$ . This sequence is very similar to the behavior we find with softened gravity for $x\lesssim\varsigma$ .

In Figure 3 we show calculations for $\tilde{\psi}_{1}(x)$ similar to those in Fig. 2A but for a different disc model — axisymmetric PL disc with $p=-0.5$ . In this case unsoftened calculations (e.g. SR15) predict that disc gravity should drive prograde precession of a test particle in a smooth disc. One can clearly see that many of the features present in Fig. 2 are reproduced for this model as well: discrepancy between the TO16 model and others, $\tilde{\psi}_{1}(x)\sim x^{-1}$ scaling for $\varsigma\lesssim x\ll 1$ , decay of $\tilde{\psi}_{1}(x)$ for $\varsigma^{2}\lesssim x\lesssim\varsigma$ , and ultimate convergence to $\psi_{1}$ in a disc with no gap. The only obvious difference is the fact that $\tilde{\psi}_{1}$ does not cross zero101010This is the case for all power-law disc models with $p<0$ or $p>3$ for which the expected free precession rate is positive, see Fig. 4. for this disc model with $p=-0.5$ .

To summarize, Figs. 2, 3 indicate that secular dynamics in softened power-law discs is dictated by the delicate balance of the opposing contributions due to nearby (i.e. with $x\lesssim\varsigma$ ) and distant disc annuli (i.e. with $x\gtrsim\varsigma$ ), in qualitative agreement with the unsoftened results of Ward (1981). These figures also demonstrate that the softening prescription of TO16 yields inaccurate results due to its inability to capture the dynamical effects of disc annuli adjacent to the test-particle orbit (those with $x\lesssim\varsigma$ ), see §6.2. We will discuss additional implications of these calculations in §6.3.

3.3 Variation of disc model — $p$ and $q$

We now examine the dependence of $\psi_{1}$ and $\psi_{2}$ on the specifics of the disc model reflected in power-law indices $p$ and $q$ . Fig. 4A,B illustrates the results based on different softening prescriptions111111We do not present results obtained by the method of Teyssandier & Ogilvie (2016). assuming a softening value of $\varsigma=10^{-3}$ (for which Fig. 1A, B suggests good convergence of $\psi_{1}$ and $\psi_{2}$ ). For reference, black open circles show the expected behavior of $\psi_{1}$ and $\psi_{2}$ computed by Silsbee & Rafikov (2015) using the un-softened Heppenheimer approach.

It is clear that the softened formalisms of both Touma (2002) and Hahn (2003) perfectly reproduce the expected behavior of the pre-factors $\psi_{1}$ and $\psi_{2}$ as a function of $p$ and $q$ (i.e. for various PL disc models). On the other hand, the prescription of Tremaine (1998) predicts a behavior of $\psi_{1}$ and $\psi_{2}$ only in qualitative agreement with the expected results: the computed values of secular coefficients deviate by tens of per cent from that of SR15. For all values of $p$ and $q$ , the formalism of Tremaine (1998) yields an additional positive contribution to $\psi_{1}$ equal to $1/2\pi$ and a negative contribution to $\psi_{2}$ equal to $-1/\pi$ (these offsets are highlighted in Fig. 4A,B by scale bars). Although these differences are not very significant, they lead to (1) predicting a wrong sign for the test-particle free-precession rate for $p\approx 0$ or $p\approx 3$ (for which SR15 yields $\psi_{1}\approx 0$ ), and (2) a mismatch of tens of per cent between the disc-driven forced eccentricity oscillations, $e_{p}^{m}/e_{d}(a)=|\psi_{2}/\psi_{1}|$ , and the expectations based on SR15. The latter point is illustrated in Figure 4C.

4 Comparison: non-Power-Law Discs

We now turn our attention to the performance of the different softening prescriptions for more general discs. Namely, we focus on two apse-aligned, non-PL disc models previously studied by Davydenkova & Rafikov (2018) based on the unsoftened Heppenheimer method. The dynamics in such non-PL discs, according to DR18, differ from the PL discs in a very important way: the free-precession of test-particles can naturally change from retrograde to prograde (and vice versa) within such discs. Furthermore, an important feature of the models considered below is that $\Sigma_{d}$ smoothly goes to zero at finite radii in a manner that does not give rise to the edge effects, see DR18 and §5.

4.1 Quartic Disc Model

We start by looking at the secular dynamics in the potential of a Quartic disc characterized by the surface density

[TABLE]

and linear eccentricity profile in the form

[TABLE]

for $a_{\rm in}\leq a\leq a_{\rm out}$ (with $a_{\rm in}=0.1$ AU, $a_{\rm out}=5$ AU), where $\tilde{\Sigma}_{0}=1153$ g cm*-2* and $\tilde{e}_{0}=0.01$ are normalization constants (one of the models in DR18).

Figure 5 summarizes the salient features of secular dynamics in the potential of such a disc adopting a softening value of $\varsigma=10^{-3}$ . It shows the excellent agreement between the radial profiles of $A_{d}$ , $B_{d}$ and $e_{p}^{m}$ computed using the un-softened calculations of Davydenkova & Rafikov (2018) and those computed using softening prescriptions of Touma (2002) and Hahn (2003). Similar to the case of PL discs, we find that the softening prescription of Tremaine (1998) yields results which agree qualitatively with the expected results but differ quantitatively. Deviations of $A_{d}$ and $B_{d}$ computed using this model from Davydenkova & Rafikov (2018), in particular, modify the locations at which $A_{d}$ and $B_{d}$ become zero. This explains the slight shift in the semi-major axes at which $e_{p}^{m}=2B_{d}/A_{d}$ goes through zero or diverges, see Figure 5.

The difference between the Tremaine (1998) and Touma (2002) calculations illustrated here could be relevant for understanding the quantitative differences between the studies of Tremaine (2001) and Gulati et al. (2012) who analyzed the slow ( $m=1$ ) modes supported by softened Kuzmin discs with softening prescriptions $b\propto r$ and $b=\text{const}$ respectively.

4.2 Gaussian Rings

Next we investigate secular dynamics in the potential of another disc model from DR18 — a Gaussian ring with the surface density profile

[TABLE]

centered around $a_{c}=1.5$ AU with width $w_{c}=0.18$ and surface density $\tilde{\Sigma}_{0}=100$ g cm*-2* at $a_{c}$ . The eccentricity profile is still given by Eq. (24).

In Figure 6 we plot the behavior of the corresponding $A_{d}$ , $B_{d}$ and $e_{p}^{m}$ for the three (convergent) softened formalisms with $\varsigma=10^{-3}$ , together with those of unsoftened Heppenheimer method (DR18, in black). Once again, the results obtained using the formalisms of Touma (2002) and Hahn (2003) fall on top of the expectations. However, for this disc model the formalism of Tremaine (1998) reproduces the un-softened calculations of Davydenkova & Rafikov (2018) quite well: the relative deviations are always less than $10\%$ . This improvement will be discussed further in §6.1.

5 Effects of proximity to the disc edge

So far the disc models that we explored were either infinitely extended (§3) or had surface density smoothly petering out to zero at finite radii (§4). This allowed us to not worry about the effects of sharp disc edges — discontinuous drops of the surface density — on secular dynamics, which are known to be important (Silsbee & Rafikov, 2015; Davydenkova & Rafikov, 2018).

We now relax this assumption and examine the performance of different softening models in the vicinity of a sharp edge of the disc, where surface density drops discontinuously from a finite value to zero at a finite semi-major axis $a=a_{\rm edge}$ . To that effect we analyze the behavior of secular coefficient $A_{d}$ computed using the formalism of Hahn (2003) (we verified that softening prescriptions of Touma (2002) and Tremaine (1998) give very similar results in the limit $\varsigma\rightarrow 0$ ) for different values of softening (results for $B_{d}$ are very similar) near the disc edge. Figure 7 shows the run of $A_{d}$ near the inner edge $a_{\rm in}$ of the disc for particles both inside ( $a_{p}<a_{\rm in}$ ) and outside ( $a_{p}>a_{\rm in}$ ) the disc as predicted by the formalism of Hahn (2003). The calculation assumes circular PL disc with $p=1$ and $\Sigma_{0}=100$ g cm*-2* extending between $a_{\rm in}=1$ AU to $a_{\rm out}=10$ AU, where we have set $a_{0}=a_{\rm out}$ (Eq. 16).

The unsoftened calculations based on Heppenheimer (1980) invariably predict that the free eccentricity precession rate $A_{d}$ , as well as $B_{d}$ , should diverge as the sharp edge of the disc is approached (e.g. SR15, DR18). Tremaine (2001) also found precession rate to diverge near the edge of a Jacobs-Sellwod ring (Jacobs & Sellwood, 2001). This is indeed the case as shown by the dashed curve computed using SR15.

The softened calculation using Hahn (2003) does largely reproduce this behavior. However, we find that very close to the ring edge (at $|a-a_{\rm in}|/a_{\rm in}\sim 10^{-3}$ ) the agreement is achieved only for $\varsigma\leq 10^{-4}$ , which is considerably smaller than the values ( $\varsigma\sim 10^{-2}$ ) required to reproduce the dynamics of particles far from the disc edges, $a_{\rm in}\ll a_{p}\ll a_{\rm out}$ , see Fig. 1. For $\varsigma=10^{-2}$ the softened calculation predicts $A_{d}$ different from the SR15 results near the disc edge by more than an order of magnitude. Thus, accurately capturing secular dynamics near the sharp edges of discs/rings requires using very small values of softening121212On the other hand, this condition is relaxed when the edge is not exactly sharp but rather has a finite width $\Delta r$ over which the disc surface density smoothly peters out to zero; in this case $\varsigma$ only needs to be $\lesssim\Delta r/r$ .. This finding could be problematic, for instance, for numerical modeling of planetary rings, often found to have very sharp edges (Graps et al., 1995; Tiscareno, 2013).

Note that in Fig. 7 softened $A_{d}$ passes through zero exactly at $a_{\rm in}$ , showing two sharp peaks of opposite signs just around this radius. Similar behavior was found by Davydenkova & Rafikov (2018) for zero-thickness discs with $\Sigma_{d}$ dropping sharply but continuously near the edge, demonstrating that variation of the sharpness of the edge is akin to softening gravity. In the case of truly zero-thickness disc and no softening (e.g. SR15) the segment of $A_{d}$ curve connecting the two peaks turns into a vertical line at $a_{\rm in}$ .

Similar divergent behavior of $A_{d}$ (and $B_{d}$ ) arises also at the outer edge of the disc considered in Fig. 7 and, in general, at any radius within a disc where $\Sigma_{d}(a)$ exhibits a discontinuity.

Finally, we note that the dynamics of particles orbiting outside the disc (where $\Sigma_{d}(a)=0$ ) is successfully reproduced by the classical Laplace-Lagrange theory without adopting any softening prescription (e.g. see Petrovich et al., 2019). Indeed, outside the radial extent of the disc semi-major axis overlap (i.e. $a_{p}=a$ ) is naturally excluded thus avoiding the classical singularity. Outside the disc the unsoftened calculations based on the Heppenheimer method (e.g. SR15, DR18) reduce to the Laplace-Lagrange theory exactly.

6 Discussion

Results of previous sections reveal a diversity of outcomes when different softening models are applied. Two models — those of Hahn (2003) and Touma (2002) — successfully reproduce the un-softened calculations based on the Heppenheimer method in the limit of zero softening. In the same limit, the formalism of Tremaine (1998) yields convergent results which are, however, different from the un-softened calculations, typically by tens of per cent. Finally, the softening method of Teyssandier & Ogilvie (2016) does not lead to convergent results in the limit of vanishing softening parameter. Interestingly, the two successful models (Hahn, 2003; Touma, 2002) have been derived using rather different underlying assumptions (see §2.1.2 & 2.1.3), producing different mathematical expressions for $\phi_{ij}$ (see Table 1), and yet their results are consistent with the un-softened calculations as $\varsigma\to 0$ .

To understand this variation of outcomes, we developed a general framework for computing secular coefficients $\phi_{ij}$ (thus fully determining the softened secular model via Eqs. (4)-(6)) given an arbitrary softened two-point interaction potential in the form (10). This procedure involves orbit-averaging the softened potential along the particle trajectories; its details are presented in Appendix A. There is also an alternative approach, sketched in Appendix A.4, which assumes the disc to be a continuous entity from the start. Both of them arrive at the same expressions for $R_{d}$ .

Using these results we show in Appendix B that the expressions for $\phi_{ij}$ found by Touma (2002) and Hahn (2003) can be recovered exactly using this general framework if we set $\mathcal{F}(r_{1},r_{2})=b_{c}^{2}$ and $\mathcal{F}(r_{1},r_{2})=H^{2}(r_{1}^{2}+r_{2}^{2})$ , respectively, in the expression (10) for the two-point potential. This approach also allows us to address some of the questions raised above, which we do in §6.1 & §6.2 below.

6.1 On the softening prescription of Tremaine (1998)

Results of §3 & §4 indicate that the softening prescription of Tremaine (1998) – unlike that of Touma (2002) and Hahn (2003) – leads to quantitative differences compared to the un-softened calculations. We now demonstrate where these differences come from.

The form of the softened Laplace coefficient $\mathcal{B}_{s}^{(m),\mathrm{Tr}}$ defined by Eq. (12) suggests interaction potential (10) with $\mathcal{F}(r_{1},r_{2})=\beta_{c}^{2}{\rm max}(r_{1}^{2},r_{2}^{2})$ for the softening model of Tremaine (1998). In Appendix B we show that propagating this form of $\mathcal{F}(r_{1},r_{2})$ through our general framework results in the following expressions for the coefficients $\phi_{ij}$ :

[TABLE]

These expressions are different from the entries in the Table 1 for Tremaine (1998) in a single but very important way — presence of terms involving Dirac $\delta$ -function. Such terms arise because the form of $\mathcal{F}(r_{1},r_{2})$ adopted in Tremaine (1998) is not sufficiently smooth — its first derivative is discontinuous at $r_{1}=r_{2}$ , while the calculation of $\phi_{ij}$ involves second-order derivatives of $\mathcal{F}$ , see Eqs. (59)-(61), as well as Eq. (62). Such singular terms do not arise in other types of softening prescriptions examined in our work since they all use infinitely differentiable versions of $\mathcal{F}(r_{1},r_{2})$ . Thus, these terms should not be interpreted as representing some kind of “self-interaction" within the disc, they merely reflect the mathematical smoothness properties of $\mathcal{F}$ used in Tremaine (1998).

Presence of these terms in Eqs. (26)-(27) introduces corrections to coefficients $A_{d}$ and $B_{d}$ (Eqs. 5, 6) in apse-aligned discs in the form

[TABLE]

Accounting for these corrections, we confirmed that the correct (un-softened) behavior of the coefficients of $R_{d}$ can be reproduced for the non-PL discs – Quartic and Gaussian models, see §4. Note that $\delta A_{d}(a_{p})$ and $\delta B_{d}(a_{p})$ are proportional to the local disc surface density $\Sigma_{d}(a_{p})$ and $\mathcal{B}_{3/2}^{(m),\mathrm{Tr}}(\alpha=1)\sim\beta_{c}^{-2}$ , see Eq. (72). This likely explains the improved agreement between the calculations of Tremaine (1998) and Davydenkova & Rafikov (2018) for Gaussian rings (see Fig. 6), which feature mass concentration in a narrow range of radii (in contrast to the Quartic model, see Fig. 5).

For PL discs the terms proportional to $\delta$ -function in Eqs. (26)-(27) give rise to corresponding modifications of the coefficients $\psi_{1}$ and $\psi_{2}$ defined by Eqs. (18)-(19):

[TABLE]

see Eqs. (20). These corrections exactly match the offsets seen in Fig. 4 between the calculations of Tremaine (1998) and the un-softened calculations, thus explaining the origin of these uniform shifts. We also confirmed this explanation in Fig. 8, where we show the convergence of modified Tremaine (1998) coefficients to the correct un-softened values as softening is varied for 2 values of $p$ and $q$ .

To summarize, Eqs. (26)-(27) should replace the expressions given by Eq. (26) of Tremaine (1998) in applications to continuous discs. However, when considering the interaction of two individual annuli with different semi-major axes (like in the classical Laplace-Largange theory), one has $\alpha\neq 1$ and terms in Eqs. (26)-(27) containing $\delta$ -function naturally vanish, reducing $\psi_{1}$ and $\psi_{2}$ back to the expressions quoted in Tremaine (1998).

6.2 On the softening prescription of Teyssandier &

Ogilvie (2016)

We now turn our attention to the model of Teyssandier & Ogilvie (2016) trying to understand its distinct (divergent) behavior. From the expression for $\mathcal{B}_{s}^{(m),\mathrm{TO}}$ in Eq. (15) one infers that this model features softening parameter in the form $\epsilon^{2}(\alpha)=S^{2}\alpha$ . To soften secular interaction Teyssandier & Ogilvie (2016) directly substituted $b_{3/2}^{(m)}$ in the classical expressions (8, 9) for $\phi_{ij}^{\rm LL}$ with $\mathcal{B}_{3/2}^{(m),\mathrm{TO}}$ , see §2.1.4; this simple swap of Laplace coefficients has not been justified rigorously.

On the other hand, in Appendix B we show that softening parameter in the form $\epsilon^{2}(\alpha)=\varsigma^{2}\alpha$ corresponds to softening function $\mathcal{F}(r_{1},r_{2})=\varsigma^{2}r_{1}r_{2}$ in the two-point potential (10), see Eq. (55). Propagating such a form of $\mathcal{F}(r_{1},r_{2})$ through our general framework in Appendix A, we find the following expressions for the coefficients $\phi_{ij}$ with $\varsigma=S$ (Appendix B):

[TABLE]

Approach of Teyssandier & Ogilvie (2016) accounts for only the first terms in Eqs. (32), (33), with coefficients which are $O(S^{0})$ , see Table 1. However, as we show below, the correct behavior of $\phi_{ij}$ as $S\to 0$ is guaranteed only when all the terms present in the above expressions are taken into account.

To demonstrate this, in Figure 8 we repeat the same convergence study as in §3.1 but with the modified $\phi_{ij}$ given by Eqs. (32) – (33). One can see see that the correct implementation of the softening $\epsilon^{2}(\alpha)=S^{2}\alpha$ proposed by TO16 leads to the recovery of the expected test-particle dynamics in infinite PL discs; this is very different from the divergent behavior obvious in Fig. 1C, D. Similar to Hahn (2003) and Touma (2002), both $\psi_{1}$ and $\psi_{2}$ smoothly converge to their expected unsoftened values in the limit of $S\to 0$ for various PL disc models (i.e. $p$ and $q$ ). Further tests using other disc models, looking at the edge effects, etc. reinforce this conclusion.

This discussion strongly suggests that for any adopted form of softening, the expansion of the secular disturbing function must be performed following a certain rigorous procedure 131313An analogous method is to modify the literal expansion of disturbing function (see Murray & Dermott, 1999, Ch. 6) to account for softened interactions (e.g. Tr98, Lee et al., 2019, H03). This could be done by replacing $b_{1/2}^{(m)}$ with $\mathcal{B}_{1/2}^{(m)}$ in Eq. (7.1) of Murray & Dermott (1999) before applying the derivatives with respect to $\alpha$ . We note that this procedure could apply for all $\mathcal{F}(r_{1},r_{2})$ with continuous first derivatives satisfying $D_{1}+D_{2}=-1$ ; see Appendix A. as done, for instance, in Appendix A. In other words, a direct replacement of the classical Laplace coefficients $b_{3/2}^{(m)}$ in Eq. (1) with their softened analogues is, evidently, not sufficient for obtaining a well-behaved softened version of Laplace-Lagrange theory for co-planar discs.

6.3 Implications for numerical applications

In numerical studies of secular dynamics, self-gravitating discs are often treated as a collection of $N$ eccentric annuli (rings), with prescribed spacing (justified by the constancy of the semi-major axis), interacting gravitationally with each other (e.g. Touma et al., 2009; Batygin, 2012). This representation approximates a continuous particulate or fluid disc in the limit of $N\rightarrow\infty$ .

Computational cost associated with the evaluation of mutual ring-ring interactions in this setup, going as $\mathcal{O}(N^{2})$ , imposes limitations on the number of rings that can be used in practice. This is typically not a problem for the un-softened calculations, which converge to the expected full disc result even with a relatively coarse radial sampling of the integral contribution to e.g. the precession rate. Indeed, purple curves in Figures 2 & 3 demonstrate this by showing the un-softened $\tilde{\psi}_{1}(x)$ and $\tilde{\psi}_{2}(x)$ computed without accounting141414Note that, technically, in the un-softened case this mathematical procedure is not equivalent to introducing an actual physical gap in the disc, as the latter would result in additional boundary terms. for the contributions from $a_{d,i}<a_{p}<a_{d,o}$ (see §3.2) to the integral terms in the un-softened expressions of Davydenkova & Rafikov (2018). These curves converge to the correct full disc result without exhibiting large variations in $\tilde{\psi}_{1}(x)$ and $\tilde{\psi}_{2}(x)$ , typical for softened cases.

On the contrary, the results for the softened gravity presented in §3.2 do elicit concern about the number of rings $N$ that is needed to accuratly capture the eccentricity dynamics of continuous razor-thin discs. Indeed, Figs. 2 and 3 reveal that the expected secular dynamics can be recovered using various softened gravity prescriptions only when one properly accounts for the gravitational effects of all disc annuli, including those very close to the orbit of particle under consideration. Indeed, we demonstrated that to reproduce both the magnitude and the sign of e.g. the precession rate, the distance $\Delta a$ separating a given test-particle orbit from nearest neighboring inner and outer disc rings should be quite small, $\Delta a/a_{p}\lesssim 0.1\varsigma^{2}$ . Only then does the delicate cancellation of large (in magnitude) contributions produced by different parts of the disc recovers the expected result. Thus, the separation between the modeled disc rings has to be substantially lower than the softening length itself ( $\varsigma a_{p}$ ), meaning that $N$ has to be very large, $N\gtrsim 10\varsigma^{-2}$ . This could easily make numerical studies of the eccentricity dynamics in discs very challenging.

We further confirmed this expectation by studying the convergence of disc-driven free precession rate in numerically discretized softened discs to the precession rate $A_{d}$ computed exactly for continuous softened discs (Eqs. 5, 18). To this end, we represented a given disc model as a collection of $N$ logarithmically-spaced rings, and measured the agreement between the radial profiles of theoretical and numerical results for $A_{d}$ (or $\psi_{1}$ for PL discs) by using the following global metric151515For PL discs, we neglect rings within $10\%$ of disc edges when computing $\mathcal{M}(\psi_{1})$ .

[TABLE]

Here $f_{\rm num}(a_{i})$ is the value of the metric basis (e.g. precession rate $A_{d}$ ) evaluated at the position $a_{i}$ of $i$ th ring by summing up the contributions of all other rings in the disc, while $f_{\rm theor}(a_{i})$ is the analogous quantity computed in the limit of a continuous disc, i.e. as $N\to\infty$ (it is given by the non-discretized version of Eq. (5) if $f=A_{d}$ , or Eq. (18) if $f=\psi_{1}$ ). Repeating this calculation for various combinations of $(N,\varsigma)$ , we can determine the smallest number of rings $N(\varsigma)$ that ensures the desired convergence to within, e.g. $\sim 10\%$ (i.e. $\mathcal{M}(f)\sim 0.1$ ), for a given value of softening $\varsigma$ .

Figure 9 depicts a sample of the results obtained using the softening methods of Hahn (2003), Tremaine (1998) and (rectified) Teyssandier & Ogilvie (2016) (see §6.2) for various axisymmetric disc models as indicated in the legend161616We exclude the softening method of Touma (2002) from this analysis as it introduces additional complexity due to the nature of softening parameter; $\epsilon^{2}=b^{2}/\textrm{max}(a_{1}^{2},a_{2}^{2})$ , see §2.1.2.. Figure 9 shows that as $\varsigma\to 0$ , the number of rings scales as $N\sim C\varsigma^{-\chi}$ with171717For example, the curve computed using the (corrected) model of Teyssandier & Ogilvie (2016) has $C=10.9$ and $\chi=1.91$ , while the one for Quartic disc has $C=7.2$ and $\chi=1.75$ . $C\sim 10$ and $\chi\approx(1.8-1.9)$ . The only notable exception is the Gaussian ring, for which convergence is faster (i.e. $N\propto\varsigma^{-1.5}$ ), probably because of mass concentration in a narrow range of radii.

We note that the proportionality constant $C$ in the $N(\varsigma)$ relation is not perfectly defined in the sense that it depends on the (i) desired accuracy (roughly inversely proportional to $\mathcal{M}(f)$ ), (ii) adopted metric of accuracy (mild dependence), and (iii) softening prescription used — Fig. 9 shows that discretized calculations using softening model of Hahn (2003) require substantially lower (by $\sim 2$ ) number of annuli than those using the models of Teyssandier & Ogilvie (2016) and Tremaine (1998). Nevertheless, these results further reinforce the requirement of large number of rings, with $N\sim\varsigma^{-2}$ , to capture the expected secular eccentricity dynamics in nearly-Keplerian discs.

Qualitatively similar results were stated in Hahn (2003) who showed that the secular effects of a continuous disc can be recovered only when the disc rings are sufficiently numerous that their radial separation is below the softening length. Although, interestingly, Hahn (2003) and Lee et al. (2019) claimed good convergence of the precession rate to the expected value already for $N\sim\mathcal{O}(\varsigma^{-1})$ (however, note that Lee et al. (2019) also included effects of gas pressure in their calculations, in addition to disc gravity). In our case, the condition on the separation between disc rings motivated by Figs. 2 & 3 (i.e. $\Delta a/a_{p}\lesssim 0.1\varsigma^{2}$ ), along with the results presented in Fig. 9, indicate that accurate representation of eccentricity dynamics in a cold, razor-thin disc requires a very large number of rings $N$ whenever small values of the softening parameter are used.

As we have shown in §5, very small values of softening $\varsigma\lesssim 10^{-3}$ are, in fact, necessary to accurately capture eccentricity dynamics near the sharp edges of thin discs. This suggests that $N$ has to be prohibitively large when softened gravity is applied e.g. to study the dynamics of planetary ring (Goldreich & Tremaine, 1979; Chiang & Goldreich, 2000; Pan & Wu, 2016), which are known to have sharp edges.

6.4 Further generalizations and extensions

All calculations in this work are based on the expansion of the secular disturbing function $R_{d}$ due to a coplanar disc — softened and unsoftened — to second order in eccentricities. This approximation may yield inaccurate results when the disc or particle eccentricities are high, e.g. in the vicinity of secular resonances where $A_{d}(a_{p})=0$ (Davydenkova & Rafikov, 2018), see Figs. 5, 6. Such situations may necessitate a higher-order extension of the disc potential.

Such an exercise was pursued recently by Sefilian & Touma (2019) who presented a calculation of $R_{d}$ to $4$ th order in eccentricities based on the un-softened method of Heppenheimer (1980). The general framework for calculating $R_{d}$ with arbitrary softening prescriptions presented in Appendix A can also be extended to higher order in eccentricities in similar way181818Another way to calculate the softened disturbing function for arbitrarily high eccentricities is to numerically compute the ring-ring interaction potential, as was done by Touma et al. (2009)., see e.g. Touma & Sridhar (2012). We expect that conclusions similar to those drawn from our analysis in §3-5 will also apply to the higher-order expansions.

Additionally, although we only analyzed coplanar configurations in this work, the general framework presented in Appendix A may be extended to account for non-coplanar configurations and study the inclination dynamics.

7 Summary

In this work we investigated the applicability of softened gravity for computing the orbit-averaged potential of razor-thin eccentric discs. We compared disc-driven secular dynamics of coplanar test-particles computed using softening prescriptions available in the literature with the calculations based on the unsoftened method of Heppenheimer (1980). Our findings are summarized below.

•

We confirmed that the softening methods of both Touma (2002) and Hahn (2003) correctly reproduce eccentricity dynamics of razor-thin discs in the limit of vanishing softening parameter $\varsigma$ for all disc models.

•

The softening prescription proposed in Tremaine (1998) yields convergent results as $\varsigma\rightarrow 0$ . However, quantitative differences of up to $\sim(20-30)\%$ from the unsoftened calculations are observed. We demonstrate that these differences arise because of the insufficient smoothness of the inter-particle interaction assumed in Tremaine (1998).

•

The softening formalism suggested in Teyssandier & Ogilvie (2016) does not result in convergent results in the limit of zero softening.

•

Very small values of the (dimensionless) softening parameter are required for correctly reproducing secular eccentricity dynamics near sharp edges of disks/rings.

•

We developed a general analytical framework for computing the secular disturbing function between two co-planar rings with arbitrary interaction potential of rather general form (Eq. 10). This framework accurately reproduces the orbit-averaged razor-thin disc potential as $\varsigma\to 0$ for a wide class of softened gravity models.

•

Using this general framework, we demonstrated that an accurate implementation of the softened potentials suggested in both Tremaine (1998) and Teyssandier & Ogilvie (2016) leads to the recovery of the expected dynamical behavior in the limit of small softening.

•

Our results suggest that the numerical treatments of the secular eccentricity dynamics in softened, nearly-Keplerian discs must obey important constraints. Namely, a fine numerical sampling (i.e. large number $N$ of discrete annuli representing the disc, with $N\sim C\varsigma^{-\chi}$ , $C\sim O(10)$ , $1.5\lesssim\chi\lesssim 2$ ) is required to ensure that the correct secular behavior is properly captured by such calculations when $\varsigma$ is small. This finding has important ramifications for numerical treatments of planetary rings with sharp edges.

In the future our results for the disc-driven eccentricity dynamics may be extended to higher order in eccentricity, as well as generalized for treating inclination dynamics.

Acknowledgements

We express our gratitude to Scott Tremaine and Jihad Touma for a number of stimulating discussions, which have led to substantial improvements of the manuscript. We are also grateful to Gordon Ogilvie, Jean Teyssandier, Yoram Lithwick, and Cristobal Petrovich for useful discussions, and an anonymous referee for constructive comments. A.A.S. acknowledges a scholarship by the Gates Cambridge Trust (OPP1144), while R.R.R. was supported by NASA via grant 15-XRP15-2-0139. Open Access for this article was funded by the Bill & Melinda Gates Foundation.

Appendix A Calculation of the secular ring-ring interaction

Here we present a calculation of the secular disturbing function due to two co-planar rings interacting with each other via softened gravity in the form (10). We do not assume any specific form for the softening function $\mathcal{F}$ apart from requiring it to be a function of the instantaneous positions of interacting particles with respect to the centre of the system. We first write the ring-ring interaction function as191919Note that we do not deal with the indirect part of the potential – which is left unsoftened – as it contains only periodic terms and does not affect the secular dynamics (Murray & Dermott, 1999).

[TABLE]

where $\mathcal{F}(r_{1},r_{2})$ is an arbitrary softening function introduced to cushion the singularity which arises otherwise at null inter-particle separations. In the above expression, $f_{i}$ is the true anomaly of the $i^{th}$ ring, $\varpi_{i}$ is its longitude of periapse and $r_{i}$ is its instantaneous position, $i=1,2$ . Our goal is to obtain the orbit-averaged expansion of $\Psi$ to second order in eccentricities $e_{i}$ valid for arbitrary $\mathcal{F}(r_{1},r_{2})$ .

A.1 Expansion of the interaction function $\Psi$ around small eccentricities

Following the classical techniques of celestial mechanics (see, Plummer, 1918, Ch. XVI), we start by expanding $\Psi$ around circular orbits. Using Taylor expansion we write

[TABLE]

with

[TABLE]

where $\theta=M_{1}-M_{2}+\varpi_{1}-\varpi_{2}$ , $M_{i}$ represents the mean anomaly of the $i^{\rm th}$ ring characterized with semi-major axis $a_{i}$ , and the linear operators $D_{k}$ are given by (Plummer, 1918)

[TABLE]

Note that this expansion, as well as subsequent steps, is completely symmetric with respect to interchanging the particle indices.

Next, in order to calculate the action of the operator $\mathbb{T}$ defined by Eq. (36) on the disturbing function of circular softened rings $\Psi_{0}$ , we make use of the elliptical expansions of $r/a$ and $f-M$ ,

[TABLE]

to multiply individual terms appearing in $\mathbb{T}$ , keep the ones up to second order in eccentricities, and drop all terms which do not contain the difference of mean anomalies, $k(M_{1}-M_{2})$ , as they are evidently periodic and vanish upon orbit-averaging. Performing this procedure and dropping an irrelevant constant term, one can demonstrate that $\Psi$ reduces to

[TABLE]

where the operators $\mathbb{A}$ , $\mathbb{B}$ and $\mathbb{C}$ acting on $\Psi_{0}$ are defined as

[TABLE]

We have used the fact that $\cos(M_{1}-M_{2})=\cos\theta\cos(\varpi_{1}-\varpi_{2})$ and $\sin(M_{1}-M_{2})=\sin\theta\cos(\varpi_{1}-\varpi_{2})$ in the secular regime (Plummer, 1918).

A.2 Computation of the action of relevant operators

Equipped with the expression (41) for $\Psi$ , we proceed to compute the action of operator $\mathbb{T}$ on $\Psi_{0}$ prior to orbit-averaging the resultant expression. With this in mind, we compute the action of several operators appearing in the definitions of $\mathbb{A}$ , $\mathbb{B}$ and $\mathbb{C}$ on $\Psi_{0}$ and list them below:

[TABLE]

where for conciseness we have written $\mathcal{F}$ instead of $\mathcal{F}(a_{1},a_{2})$ . Here, it is worthwhile to mention that, as far as the expansion technique is concerned, the terms $\partial_{i}\mathcal{F}~{}(\text{with}~{}i=1,~{}2)$ appearing in the above expressions are the only difference brought upon by softening the Newtonian point-mass interaction (Eq. 35). Another set of operators useful in computing $\mathbb{T}\Psi_{0}$ is the following:

[TABLE]

which can be obtained by making use of the identity $(D_{1}+D_{2}+1)\Psi_{0}=\frac{1}{2}(2\mathcal{F}-a_{1}\partial_{1}\mathcal{F}-a_{2}\partial_{2}\mathcal{F})\Psi_{0}^{3}$ . Here, we note that for all softening functions $\mathcal{F}$ for which $2\mathcal{F}-a_{1}\partial_{1}\mathcal{F}-a_{2}\partial_{2}\mathcal{F}=0$ , one finds $D_{1}+D_{2}=-1$ . Consequently, in such cases, the operators $D_{1}(D_{1}+1)$ and $D_{2}(D_{2}+1)$ become identical rendering $\mathbb{A}\Psi_{0}=\mathbb{B}\Psi_{0}$ (since $D_{3}^{2}=D_{4}^{2}$ , see Eqs. (42) and (44)). As a result, the resultant orbit-averaged disturbing function (41) is symmetric in $e_{1}$ and $e_{2}$ , similar to the case of classical Laplace-Lagrange theory. This is not true in general, for instance, when $\mathcal{F}(r_{1},r_{2})=$ const $\neq 0$ .

A.3 Orbit-averaging the interaction function $\Psi$

The expressions (44)-(51) allow the computation of $\Psi=\mathbb{T}\Psi_{0}$ , which needs to be time-averaged in order to recover the secular disturbing function. We do not show the cumbersome collated expression for $\mathbb{T}\Psi_{0}$ and proceed to the final step of orbit-averaging, which will conclude our derivation. In short, our goal is to compute

[TABLE]

which essentially reduces to computing the individual terms $\langle\mathbb{A}\Psi_{0}\rangle$ , $\langle\mathbb{B}\Psi_{0}\rangle$ and $\langle\mathbb{C}\Psi_{0}\rangle$ . At the outset, it is important to note that each of the terms appearing in $\mathbb{T}\Psi_{0}$ (through $\mathbb{A}\Psi_{0}$ , $\mathbb{B}\Psi_{0}$ and $\mathbb{C}\Psi_{0}$ , or the operators they entail) are proportional to $\cos(m\theta)\Psi_{0}^{2s}$ . By making use of $\alpha={a_{<}}/{a_{>}}$ , where $a_{<}={\rm min}(a_{1},a_{2})$ and $a_{>}={\rm max}(a_{1},a_{2})$ , this combination can be reduced to

[TABLE]

For that reason, calculation of the orbit-averaged $\Psi$ (by integrating over $d\theta$ ) yields integrals of the form

[TABLE]

which is the generalization of the classical Laplace coefficients $b_{s}^{(m)}$ (recovered when $\mathcal{F}(a_{1},a_{2})=0$ , see Eq. 2) with the dimensionless softening parameter

[TABLE]

see Eq. (7). Employing this notation, we present the simplified expressions of $\langle\mathbb{A}\Psi_{0}\rangle$ , $\langle\mathbb{B}\Psi_{0}\rangle$ and $\langle\mathbb{C}\Psi_{0}\rangle$ obtained as a result of orbit-averaging:

[TABLE]

In equations (56)-(58), we have defined the dimensionless functions $T_{i}(\alpha)$ such that

[TABLE]

where, as before, $\mathcal{F}\equiv\mathcal{F}(a_{1},a_{2})$ , $\alpha=a_{<}/a_{>}$ and $\partial_{i}\equiv\partial/\partial a_{i}$ . Note that the expressions for $\phi_{11}$ and $\phi_{22}$ swap definitions upon replacing $a_{1}$ by $a_{2}$ , whilst keeping $\alpha<1$ by construction. This can be understood by first noting that functions $T_{i}$ with $i=1,~{}2,~{}3,$ and $5$ are invariant under $a_{1}\leftrightharpoons a_{2}$ while, at the same time, $T_{4}$ and $T_{7}$ (appearing in the second line of Eq. (56)) translate to $T_{6}$ and $T_{8}$ (appearing in the second line of Eq. (57)); and vice versa.

These identities, when combined, yield the desired expression of $\langle\Psi\rangle=\langle\mathbb{T}\Psi_{0}\rangle$ ; see Eqs. (41)-(43). Subsequently, the softened ring-ring disturbing function in the form (11) is recovered, with the coefficients $\phi_{ij}$ defined by Eqs. (56) – (58). This completes our calculation of the secular ring-ring interaction between two softened coplanar rings, up to second order in eccentricity and valid for arbitrary softening functions $\mathcal{F}(r_{1},r_{2})$ .

Note that in the absence of softening (i.e. $\mathcal{F}(r_{1},r_{2})=0$ ) $T_{i}=0$ for all $i$ and the classical expressions for $\phi_{11}^{\rm LL}$ , $\phi_{22}^{\rm LL}$ and $\phi_{12}^{\rm LL}$ — Eqs. (8)-(9) — are recovered. Finally, we mention that the expansion technique exploited here can be used to recover the orbit-averaged disturbing function valid to arbitrary order in eccentricity, as well as inclinations.

A.4 Alternative calculation: secular disc-particle interaction

Calculations presented above describe the orbit-averaged coupling between the two individual annuli, which subsequently need to be integrated over the semi-major axes of the disc elements to represent the effect of a continuous disc. In principle, one can also arrive at the expressions (4) by assuming a continuous mass distribution in the disc from the start and performing a calculation similar to that in Davydenkova & Rafikov (2018). Namely, one would need to compute $R_{d}=\langle G\int_{S}\Sigma({\bf r}_{d})\Phi({\bf r}_{d},{\bf r}_{p})dS\rangle$ , where $\Phi$ is the interaction potential given by equation (10), angle brackets indicate averaging over the orbit of the test particle given by ${\bf r}_{p}$ and integration is carried out over the full surface of the disc $S$ with ${\bf r}_{d}$ denoting the location of a disc element. To obtain the expression for $R_{d}$ accurate to second order in eccentricities one would need to expand $\Phi({\bf r}_{d},{\bf r}_{p})$ to second order in particle and disc eccentricities by e.g. writing $r_{p}=a_{p}(1-e_{p}\cos E_{p})$ , where $E_{p}$ is the eccentric anomaly of the particle orbit. This expansion should explicitly account for the dependence of $\mathcal{F}$ on ${\bf r}_{d}$ and ${\bf r}_{p}$ . Averaging the resulting expressions over $E_{p}$ , one would arrive at the proper expression for $R_{d}$ in the form (4).

In particular, after a lengthy but straightforward calculation this method gives the following expression for the disc-driven precession rate:

[TABLE]

where prime denotes differentiation with respect to $a_{p}$ (e.g. $\mathcal{F}^{\prime}=\partial\mathcal{F}/\partial a_{p}$ ), $a_{>}=\max(a_{p},a)$ , $\alpha=\min(a_{p},a)/\max(a_{p},a)$ and integration is done over the semi-major axis $a$ of the disk elements. Calculation of the non-axisymmetric part of $R_{d}$ resulting from non-zero disk eccentricity (i.e. $\bm{\mathsf{B}}_{d}$ ) is somewhat more tedious but can nevertheless be done similar to Davydenkova & Rafikov (2018).

Appendix B Specific Cases of $\mathcal{F}(r_{1},r_{2})$

The general framework developed in Appendix A allows us to recover the expressions of $\phi_{ij}$ arrived at by Touma (2002) and Hahn (2003) upon specifying certain functional forms of $\mathcal{F}(r_{1},r_{2})$ . Indeed, Touma (2002) performed the same calculations as presented in Appendix A for the case of Plummer potential – $\mathcal{F}(r_{1},r_{2})=b_{c}^{2}$ – to second order in eccentricities, and later to fourth order in eccentricities (Touma & Sridhar, 2012). Furthermore, we find that the results obtained by Hahn (2003) can be recovered from our general framework by setting $\mathcal{F}(r_{1},r_{2})=H^{2}(r_{1}^{2}+r_{2}^{2})$ . For reference, the functional forms of $T_{i}$ for these forms of $\mathcal{F}(r_{1},r_{2})$ , along with their softening parameters $\epsilon^{2}(\alpha)$ , are summarized in Table 2, which can be used to show that Eqs. (56)-(58) reduce to those in Table 1 after some algebra with the aid of the recursive relationships for $\mathcal{B}_{s}^{(m)}$ presented in Appendix C.

As to the formalism of Teyssandier & Ogilvie (2016), we find, using their softening prescription of $\mathcal{F}(r_{1},r_{2})=S^{2}r_{1}r_{2}$ , that our general framework yields $\phi_{ij}$ expressions different from those reported by Teyssandier & Ogilvie (2016). Indeed, we first note that in this case, $T_{1}=T_{4}=T_{6}=0$ (Table 2) rendering the expressions of $\phi_{11}$ and $\phi_{22}$ identical such that

[TABLE]

Using the recursive relationships listed in Appendix C.1, the above expression can be simplified further. Indeed, Eq. (67) with $m=1$ and $s=5/2$ and Eq. (66) with $m=1$ and $s=3/2$ read

[TABLE]

respectively. Inserting the above two identities in Eq. (63) one arrives at Eq. (32). Similarly, the expression of $\phi_{12}$ (Eq. 58) can be simplified with the aid of Eq. (68) (with $m=0,~{}s=3/2$ ), Eq. (67) (with $m=2,~{}s=5/2$ ) and Eq. (66) (with $m=2,~{}s=3/2$ ) resulting in Eq. (33) after some algebra. As discussed in §6.2, the terms in Eqs. (32)-(33) explicitly proportional to $S^{2}$ are absent in the original formulation of Teyssandier & Ogilvie (2016) (see Table 1).

Similarly, for the formalism of Tremaine (1998), propagating their functional form of $\mathcal{F}(r_{1},r_{2})=\beta_{c}^{2}{\rm max}(r_{1}^{2},r_{2}^{2})$ through our general framework, we arrive at the expressions for $\phi_{ij}(\alpha)$ differing from those reported in Tremaine (1998) in a very special way: we find $\phi_{ij}$ to contain additionl terms proportional to $T_{3}(\alpha)\sim\delta(\alpha-1)$ , where $\delta(x)$ is the Dirac delta-function. Such terms are absent in the original formulation of Tremaine (1998) (see Tables 1, 2). Emergence of these terms can be easily demonstrated by first noting that in this case $\phi_{11}=\phi_{22}$ (as $T_{1}=T_{4}=T_{6}=0$ ), employing the recursive relationships for Laplace coefficients (in a similar order as done above for TO16) to simplify the general expressions of $\phi_{11}(=\phi_{22})$ and $\phi_{12}$ , and finally arriving at Eqs. (26), (27). The ramifications of this finding is discussed in Section 6.1.

Appendix C Generalized Laplace coefficients

As demonstrated in Appendix A, softening the Newtonian point-mass potential by an arbitrary function $\mathcal{F}(r_{1},r_{2})$ modifies the definition of the Laplace coefficients as shown by Eqs. (7), (54) by the introduction of a softening parameter $\epsilon^{2}(\alpha)=a_{>}^{-2}\mathcal{F}(a_{1},a_{2})$ (Eq. 55), $0\leq\alpha=a_{<}/a_{>}\leq 1$ . Here, we present some useful recursive relationships amongst different generalized Laplace coefficients $\mathcal{B}_{s}^{(m)}(\alpha)$ , along with their asymptotic behavior in the limits of $\alpha\rightarrow 0,~{}1$ as well as their relationship to complete elliptic integrals.

C.1 Recursive Relations

Generalizing the results for the usual (unsoftened) Laplace coefficients $b_{s}^{(m)}$ (e.g. Plummer, 1918, p. 159), the following relationships can be easily obtained for the generalized Laplace coefficients defined by Eq. (7),(54):

[TABLE]

The difference with the classical recursive relations for $b_{s}^{(m)}$ amounts to substituting the combination $1+\alpha^{2}$ appearing in the case of ordinary Laplace coefficients with $1+\alpha^{2}+\epsilon^{2}(\alpha)$ .

Another useful expression relating the generalized Laplace coefficients of arguments $\alpha$ and $\alpha^{-1}$ is

[TABLE]

Note that the above relationship is valid only as long as the softening parameter satisfies $\alpha^{2}\epsilon^{2}(1/\alpha)=\epsilon^{2}(\alpha)$ . For instance, this condition is violated when the softening parameter $\epsilon$ has no dependence on $\alpha$ (e.g. that of Tremaine (1998), see Table 1).

C.2 Asymptotic Behavior

Here we derive approximate expressions for $\mathcal{B}_{s}^{(m)}$ in the asymptotic limits; for $\alpha\rightarrow 0$ and $\alpha\rightarrow 1$ .

Case 1: In the limit of $\alpha\approx 0$ , one can factor out the term $1+\alpha^{2}+\epsilon^{2}(\alpha)$ from the integrand of $\mathcal{B}_{s}^{(m)}$ to expand the denominator around $\gamma^{-1}\approx 0$ , where $\gamma=(2\alpha)^{-1}[1+\alpha^{2}+\epsilon^{2}(\alpha)]$ . This allows us to approximate $\mathcal{B}_{s}^{(m)}$ as

[TABLE]

Using the orthogonality of the cosine functions, it is straightforward to show that

[TABLE]

Case 2: In the opposite limit of $x=1-\alpha\approx 0$ , the dominant contribution to $\mathcal{B}_{s}^{(m)}$ comes from $\theta\ll 1$ (Goldreich & Tremaine, 1980). Thus one can set $\cos(m\theta)\to 1$ in the numerator, approximate $\cos\theta\approx 1-\theta^{2}/2$ in the denominator and extend the integration limit to infinity. Furthermore, setting $\alpha=1$ (i.e. $x=0$ ) everywhere except when it appears in the combination $1-\alpha$ , the generalized Laplace coefficient can be approximated as

[TABLE]

where $\epsilon^{2}_{\alpha=1}$ is the softening parameter evaluated at $\alpha=1$ .

C.3 Relationship to elliptic integrals

Here we express the generalized Laplace coefficients $\mathcal{B}_{s}^{(m)}$ in terms of complete elliptic integrals. These expressions can be used for rapid numerical evaluation of the generalized Laplace coefficients without relying on numerical integration of Eq. (54) (or Eq. (7)). Let us write, as before, $2\alpha\gamma=1+\alpha^{2}+\epsilon^{2}(\alpha)$ and define $\chi=\sqrt{2/(\gamma+1)}$ such that, for any general softening parameter $\epsilon^{2}(\alpha)$ , we have $0\leq\chi\leq 1$ and $\gamma\geq 1$ . Now let us express $\mathcal{B}_{s}^{(m)}$ in terms of $\gamma$ to write

[TABLE]

Introducing complete elliptic integrals $\mathbf{K}(\chi)=\int_{0}^{\pi/2}\left(1-\chi^{2}\sin^{2}\phi\right)^{-1/2}d\phi$ and $\mathbf{E}(\chi)=\int_{0}^{\pi/2}\left(1-\chi^{2}\sin^{2}\phi\right)^{1/2}d\phi$ , we find that

[TABLE]

These expressions permit efficient numerical evaluation of arbitrarily softened Laplace coefficients as functions of $\alpha$ , since effective algorithms for computing $\mathbf{K}$ and $\mathbf{E}$ exist (Press et al., 2002).

Appendix D Convergence Criterion for the pre-factors of power-law discs

Astrophysical discs often extend over a few orders of magnitude in radius so that $a_{\rm out}/a_{\rm in}\gg 1$ . In such situations, far from the disc edges one can take the limit of both $\alpha_{1}=a_{\rm in}/a_{p}$ and $\alpha_{2}=a_{p}/a_{\rm out}$ going to zero, provided that the gravitational potential of a power-law disc is insensitive to the locations of the disc boundaries (see Eqs. 18, 19). Then the pre-factors $\psi_{1}$ and $\psi_{2}$ of the disturbing function converge to values depending only on the power-law indices $p$ and $p+q$ respectively, as well as on the adopted softening prescription.

The conditions on the values of $p$ and $q$ which guarantee this convergence can be determined by expanding the coefficients $\phi_{ij}(\alpha)$ , which appear in the integrands of each of $\psi_{1}$ and $\psi_{2}$ , in the limit of $\alpha\approx 0$ . Using the Taylor expansions of softened Laplace coefficients $\mathcal{B}_{s}^{(m)}$ , we determined that both $\psi_{1}$ and $\psi_{2}$ calculated using the softening methods of Hahn (2003) and Tremaine (1998) (as well as its rectified version) are convergent as long as $-1<p<4$ and $-2<p+q<5$ , respectively, for all values of softening (i.e. $H,\beta_{c}$ ). This follows from the fact that for both Hahn (2003) and Tremaine (1998) we have $\phi_{11}=\phi_{22}\sim\alpha^{2}$ and $\phi_{12}\sim\alpha^{3}$ to lowest order in $\alpha$ . These ranges of $p$ and $p+q$ are in line with the findings of Silsbee & Rafikov (2015).

As to the (rectified) softening model of Teyssandier & Ogilvie (2016), a similar exercise yields that $\phi_{11}=\phi_{22}\approx-\frac{1}{4}S^{2}\alpha+\frac{3}{8}(1+\frac{3}{2}S^{4})\alpha^{2}$ and $\phi_{12}\approx\frac{3}{2}S^{2}\alpha^{2}-\frac{15}{16}(1+5S^{4})\alpha^{3}$ which, in the limit of $S\to 0$ , translate to the same ranges for $\psi_{1}$ and $\psi_{2}$ convergence as Silsbee & Rafikov (2015). However, when $S$ is relatively large, it is trivial to show that $\psi_{1}$ and $\psi_{2}$ are convergent over limited ranges of $0<p<3$ and $-1<p+q<4$ , respectively. A similar analysis for the softening method of Touma (2002) reveals that the ranges for $\psi_{1}$ and $\psi_{2}$ convergence are in line with the findings of Silsbee & Rafikov (2015) when the corresponding softening parameter $b_{c}\to 0$ . However, when $b_{c}$ is non-zero, the ranges are narrowed down to $-1<p<2$ and $-2<p+q<3$ respectively.

Bibliography35

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1Batygin (2012) Batygin K., 2012, Nature , 491, 418 · doi ↗
2Binney & Tremaine (2008) Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition. Princeton University Press
3Chiang & Goldreich (2000) Chiang E. I., Goldreich P., 2000, Ap J , 540, 1084 · doi ↗
4Davydenkova & Rafikov (2018) Davydenkova I., Rafikov R. R., 2018, Ap J , 864, 74 · doi ↗
5Fontana & Marzari (2016) Fontana A., Marzari F., 2016, A&A , 589, A 133 · doi ↗
6Goldreich & Tremaine (1979) Goldreich P., Tremaine S., 1979, AJ , 84, 1638 · doi ↗
7Goldreich & Tremaine (1980) Goldreich P., Tremaine S., 1980, Ap J , 241, 425 · doi ↗
8Graps et al. (1995) Graps A. L., Showalter M. R., Lissauer J. J., Kary D. M., 1995, AJ , 109, 2262 · doi ↗

TL;DR

Contribution

Findings

Abstract

Peer Reviews

Videos

Potential softening and eccentricity dynamics in razor-thin, nearly-Keplerian discs

Abstract

keywords:

1 Introduction

2 Disturbing function due to a disk

2.1 Summary of existing softening models

2.1.1 Formalism of Tremaine (1998) – Tr98

2.1.2 Formalism of Touma (2002) – T02

2.1.3 Formalism of Hahn (2003) – H03

2.1.4 Formalism of Teyssandier &

2.2 The unsoftened Heppenheimer method

3 Comparison: Power-Law Discs

3.1 Behavior with respect to variation of softening

3.2 Details of convergence of different softening prescriptions

3.3 Variation of disc model — ppp and qqq

4 Comparison: non-Power-Law Discs

4.1 Quartic Disc Model

4.2 Gaussian Rings

5 Effects of proximity to the disc edge

6 Discussion

6.1 On the softening prescription of Tremaine (1998)

6.2 On the softening prescription of Teyssandier &

6.3 Implications for numerical applications

6.4 Further generalizations and extensions

7 Summary

Acknowledgements

Appendix A Calculation of the secular ring-ring interaction

A.1 Expansion of the interaction function Ψ\PsiΨ around small eccentricities

A.2 Computation of the action of relevant operators

A.3 Orbit-averaging the interaction function Ψ\PsiΨ

A.4 Alternative calculation: secular disc-particle interaction

Appendix B Specific Cases of F(r1,r2)\mathcal{F}(r_{1},r_{2})F(r1​,r2​)

Appendix C Generalized Laplace coefficients

C.1 Recursive Relations

C.2 Asymptotic Behavior

C.3 Relationship to elliptic integrals

Appendix D Convergence Criterion for the pre-factors of power-law discs

3.3 Variation of disc model — $p$ and $q$

A.1 Expansion of the interaction function $\Psi$ around small eccentricities

A.3 Orbit-averaging the interaction function $\Psi$

Appendix B Specific Cases of $\mathcal{F}(r_{1},r_{2})$