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

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, ς ≪ 1, accurate representation of eccentricity dynamics requires N ∼ Cς-χ ≫ 1, with C ∼ O(10), 1.5 ≲ χ ≳. In discs with sharp edges a very small value of the softening parameter ς (≲ 10-3) is required to correctly reproduce eccentricity dynamics near the disc boundaries; this finding is relevant for modelling planetary rings.

INTRODUCTION

Astrophysical discs orbiting a central mass M are ubiquitous in a variety of contexts – galactic, stellar, and planetary (Latter et al. 2017). In many instances, masses of such discs M 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 ≪ M, 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, fully determines the secular dynamics of the system.

For a test particle with semi-major axis a, eccentricity e, and apsidal angle due to a co-planar point mass δm orbiting with semi-major axis a, eccentricity e, and apsidal angle , upon smearing into elliptical rings, the secular dis-turbing function takes the form (Murray & Dermott 1999)

valid for a > a as well as a < a, as long as particle orbits do not cross. Here is the Laplace coefficient defined by

which obeys . Explicit time indepen-dence of δR guarantees that the semi-major axes of the sec-ularly 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 contribu-tions in the form (1) produced by all wires, which amounts to integration of δR (Eq. 1) over the radial extent of the disc:

where ain and aout are the semi-major axes of the inner and outer disc edges. In this case, provided that δR is known as a function of a, only a single integration (over the semi-major axes of the rings) is needed, significantly accelerating calculations.

Unfortunately, this straightforward procedure is illposed from the mathematical point of view. Indeed, it is well known that the Laplace coefficients featured in Eq.(1) diverge as when α → 1. This implies that the radial integration in Eq. (3) encounters an essential singularity at a = a . As a result, for a co-planar particle orbiting inside a razor-thin disc, ain ≤ a ≤ aout, this direct way of computing R does not converge to a finite value.

This divergence, as well as the pressing need for having an efficient way of computing R (via a one-dimensional integration over a only), have motivated the development of alternative analytic approaches for calculating R . 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 or-der (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 con-ditions (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 as a → a . 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.

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 Rdue 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 µ(a), eccentricity e(a), and apsidal angle 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 Σ(a) instead of µ(a); its relation to µ for arbitrary profiles of e and has been established in Statler (2001), Davydenkova & Rafikov (2018) and Davydenkova & Rafikov (in prep.). Constancy of semi-major axis in secular theory implies that µ(a) does not change in time. The same statement is true for Σ(a) to lowest order in e since µ(a) ≈ 2πaΣ(a) + O(e) (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, eccentricity e, and apsidal angle :

Here n is the test-particle mean motion , and we have introduced a two-component eccentricity vector for a test particle .

The coefficients A and B in Eq. (4) are related to the disc mass (or surface density) and eccentricity profiles in the following fashion:

where is the eccentricity vec-tor for an annular disc element.

Functions φ(α), 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 φ appearing in the literature are linear combinations of softened Laplace coefficients defined by

Table 1

The coefficients ϕ(α) 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 α is defined such that α = a where a> = max(a1, a2) and a< = min(a1, a2). 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 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(α) are defined in Eq. (7). Note that the expressions of ϕ reported in Touma (2002) have been corrected in a subsequent paper of Touma & Sridhar (2012).

Formalism	∈2(α)	ϕ11	ϕ12	ϕ22
Laplace-Lagrange	–	18αb3/21	−14αb3/22	ϕ11
Tremaine (1998) (Tr98)	βc2	182αddα+α^2d^2dα^2B1/20,Tr=18αB3/21,Tr−3αβc2B5/20,Tr	142−2αddα−α^2d^2dα^2B1/21,Tr=−14αB3/22,Tr−3αβc2B5/21,Tr	ϕ11
Touma (2002) (T02)	β^2=bc2/a>2	−58αB3/21,T+316α^2B5/20,T+38α1+α^2B5/21,T−1516α^2B5/22,T−38αβ^2αB5/20,T−B5/21,T	98αB3/20,T+18αB3/22,T−98α1+α^2B5/20,T+2116α^2B5/21,T+38α1+α^2B5/22,T+316α^2B5/23,T	−58αB3/21,T+316α^2B5/20,T+38α1+α^2B5/21,T−1516α^2B5/22,T−38αβ^2αB5/20,T−B5/21,T
Hahn (2003) (H03)	H2(1 + α2)	18αB3/21,H−3αH^22+H^2B5/20,H	−14αB3/22,H−3αH^22+H^2B5/21,H	ϕ11
Teyssandier & Ogilvie (2016) (TO16)	S2α	18αB3/21,TO	−14αB3/22,TO	ϕ11

The softening parameter ∈(α) appearing in this definition remains non-zero as α → 1, thus preventing the divergence of the softened Laplace coefficients at α = 1 (unlike the classical . The explicit form of ǫ(α) is different for every softening method considered in this work, see §2.1 and Table 1. Appendix C collates some useful rela-tions for softened Laplace coefficients , as well as their approximate asymptotic behavior and relationships to complete elliptic integrals.

The mathematical structure of R 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, µ(a) → mplδ(a − apl) (where δ(z) is the Dirac delta-function), and set softening to zero (so that ). Then one finds that R reduces to the un-softened, orbit-averaged potential δR due to a planet with mass mpl and semi-major axis apl, with the unsoftened coefficients ϕ in the form (Murray & Dermott 1999)

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 ∈ to avoid the mathematical divergence of the classical disturbing function as a → a.

We emphasize that the functional forms of ϕ are not simple replacements of appearing in the unsoftened definition (8) - (9) by . This can be seen in Table 1 where we summarize some of the expressions for ϕ(α) proposed in the literature and analyzed in this paper (see §2.1). Nev-ertheless, examination of these expressions shows that when ∈2(α) → 0, the coefficients ϕ(α) do reduce to their unsoft-ened versions 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 m1 and m2 (m ≪ M) located at r1 and r2, correspondingly, relative to the central mass, has a form

with i, j = 1, 2 and j ≠ i. Here ℱ(r1, r2) 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 r1 −r2, but also on their distances to the dominant central mass r1, r2.

Explicit demonstration of the connection between the potential (10) and R 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 ∈ featured in the definition (7) is related to ℱ via ∈2 = [max(a1, a2)]−2ℱ(a1, a2), where a1,2 are the semi-major axes of the interacting particles (see Eq. A21). The most general expressions of ϕ entering the arbitrarily softened ring-ring disturbing function,

(here i = 1, 2 and j ≠ i) is given by Eqs. (A22)-(A24) in terms of . In the above expression, we have defined a> = max(a1, a2) and a< = min(a1, a2) such that α = 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, and exterior to it. We do this because for some softening functions ℱ the coefficients ϕ(α) do not obey certain symmetry properties when a/a is replaced with a/a, see Eq. (C4). Moreover, in general ϕ11 and ϕ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 A and B, we remind the reader that A represents the precession rate of the free ec-centricity vector of a test particle in the disc potential, while B characterizes the torque exerted on the particle orbit by the non-axisymmetric component of the disc gravity. Corresponding forced eccentricity vector is e = −B/A. In par-ticular, test-particles initiated on circular orbits experience eccentricity oscillations of maximum amplitude .

As A(a) and B(a) uniquely determine R for different forms of softening, comparison of their behavior in the limit of ∈ → 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.

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 ∈2(α) and coefficients ϕ(α) are provided in Table 1.

Formalism of In 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

Here is the dimensionless softening parameter, treated as a constant, i.e. independent of distance. The physical interpretation of this manoeuvre is that β, inhibiting the formal divergence of R as a → a, 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 = βcmax(a1, a2). In Tremaine (1998) coefficients ϕ (α) were expressed as derivatives of with respect to α, see equations (26) of Tremaine (1998). These expressions, along with their versions modified using the re-cursive relations for Laplace coefficients (see Appendix C1),can be found in Table 1.

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 (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 ϕ(α) in the form of linear combinations of softened Laplace coefficients , similar to those of Tremaine (1998):

However, in Touma (2002) the softening parameter ∈2(α) = β2 is no longer a constant but depends on the distance such that β = b/max(a1, a2). Within this formalism, one can think of a disc as comprised of nested annuli with a con-stant thickness b.

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 scaleH ∼ h/a, which was assumed to be constant in that work (see also Ward 1989). Hahn (2003) demonstrated that the resultant ϕ(α) are functions of softened Laplace coefficients

with constant H ≪ 1. In other words, the softening parameter is given by ∈2(α) = H2(1 + α2) in that work. The explicit expressions for ϕ(α) in terms of are given by equations (17) of Hahn (2003).

Formalism of Teyssandier & Ogilvie (2016) – TO16

Teyssandier & Ogilvie (2016) modified the unsoftened expressions (8), (9) for by simply replacing the usual Laplace coefficients with softened versions defined such that

Thus, their softening parameter is ∈2(α) = S2α, 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 ϕ(α) are given by equations (7)-(9) of Teyssandier & Ogilvie (2016).

The aforementioned softening prescriptions have their softening parameters ∈2(α) controlled by different constants — β, b, 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 ς.

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 instanta-neous particle position r) in terms of small eccentricity of a test particle. 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 – and not – which diverge only weakly (logarithmically) as . As a result, upon in-tegrating these expressions over the radial extent of the disc,one obtains a convergent and finite result for R . 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 pot entials (e.g. SR15, Fontana & Marzari 2016, DR18), validating this approach.

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. ), truncated power-law (hereafter PL) disc as a simple example. We characterize surface density and eccentricity of such a disc by

for ain ≤ a ≤ aout, where Σ0 and e0 are the pericentric surface density and eccentricity of the disc at some reference semi-major axis a0.

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

where and the dimensionless coefficients ψ1 and ψ2 are given by

with α1 = ain/a and α2 = a/aout.

The coefficients ψ1 and ψ2 are functions of the power-law indices (p and q), any softening parameter involved (through ϕ), as well as the test-particle semi-major axis a (through α1,2). They are related to A and B via

As shown in Appendix D, for certain ranges of power-law indices p and q both ψ1 and ψ2 converge to values de-pending 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 α1,2 → 0). For p and qin these ranges (determined in Appendix D for each of the considered softened formalisms, similar to SR15), the coefficients ψ1 and ψ2 are determined by the local behavior of Σ(a) and e(a) in the vicinity of test-particle semi-major axis.

Given this, we first focus on infinitely extended (α1,2 → 0) PL discs with p and q within these ranges (we defer discussion of secular dynamics near the disc edges to §5). Then, ψ1 and ψ2 become independent of a (i.e. functions of p, q, and ς only), making them useful as simple metrics for judg-ing the validity of different models of softening.

Behavior with respect to variation of softening

Figure 1 illustrates the behavior of ψ1 and ψ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” ς for two different sets of p, q (indicated in panel B). For reference, black horizontal lines show the values of ψ1 and ψ2 expected from the calculations of SR15 using theun-softened Heppenheimer approach.

Figure 1

Behavior of the axisymmetric (ψ1, Eq. (18), top panels) and non-axisymmetric (ψ2, Eq. (19), bottom panels) components of the softened gravitational potential due to an infinite power-law disc as a function of softening ς. The calculations assume two different disc structures specified by the values of p and q shown by different line types as explained in legend. For clarity, the results obtained by the softened formalisms of Tremaine (1998), Touma (2002) and Hahn (2003) are collated in the left panels and those obtained by the softening method of Teyssandier & Ogilvie (2016) are shown in the right panels. The left panels also show the ψ1 and ψ2 obtained by SR15 not assuming any softening (black horizontal lines). See text (§3.1) for details.

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 ψ1 and ψ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 ς approaches zero; both the amplitude and sign of ψ1 and ψ2 are reproduced. It is also evident that, depending on disc model, ψ1 and ψ2 converge to values given by SR15 at different values of softening. Nevertheless, we generally find that ς ≲ 10−3 guarantees the convergence of ψ1 and ψ2 to within few per cent of the correct values for all p and q as long as ain ≪ a ≪ aout (see Figure 4).

Figure 4

Dependence of the coefficients ψ1 (panel A) and ψ2 (panel B) on the power-law disc model represented by the indices p and p + q, respectively. Panel C shows the amplitude of eccentricity oscillations (normalized by disc eccentricity e) induced by disc gravity. Results for softened formalisms of Hahn (2003) (in red), Touma (2002) (in blue) and Tremaine (1998) (in green) are computed using softening ς = 10−3. Calculations assume infinitely extended disc (i.e. no edge effects). For reference, open black circles show the profiles of ψ1, ψ2 and as computed by SR15: curves for Hahn (2003) and Touma (2002) fall on top of them, while those for Tremaine (1998) show constant offset in terms of both ψ1 and ψ2 (illustrated by scale bars in panels A,B) resulting in deviation between curves (panel C).

The same panels also indicate that ψ1(ς) and ψ2(ς) pre-dicted by the softened formalism of Tremaine (1998) (green), while converging to finite values as ς = β → 0, do not reproduce the SR15 results exactly in this limit. Indeed, one can see that even for the smallest adopted value of β = 10−3, the softening prescription of Tremaine (1998) yields ψ1 and ψ2 different by tens of per cent from SR15. It is easy todemonstrate that these quantitative differences do not vanish by further decreasing β. For instance, when p = 1, the coefficient ψ1 can be evaluated analytically as

in agreement with Panel A (E(k) is the complete elliptic integral of a second kind). At the same time, the unsoftened approach of SR15 predicts ψ1 = −1/2 for p = 1 disc. Moreover, close inspection of Fig. 1A,B shows that, in the limit of β → 0, the ψ1 and ψ2 curves computed using soften-ing model of Tremaine (1998) are offset vertically from the unsoftened calculations by 1/2π and −1/π, 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 ψ1 (Panel C) and ψ2 (Panel D) as a function of “softening”, ς = S, resulting from the approach of Teyssandier & Ogilvie (2016). There are several features to note here. First, this model pre-dicts ψ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, ψ2 is always negative, contrary to the expectations (see Panel B). Second, in the limit of S → 0, both ψ1 and ψ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 ψ1 and ψ2 diverge as the softening S → 0. Indeed, it suffices to employ the asymptotic expansion of the Laplace coefficients in the limit of α → 1 (Eq. C7) to demonstrate that both ψ1 and ψ2 (Eqs. 18 - 19) behave as

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

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 → a or α → 1. For this reason, it is interest-ing 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 ψ1 and ψ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 Σ(a) is set to zero, are at a = (1 − x)a ≤ a anda = (1− x)−1a ≥ a, respectively, with a single parameterx controlling the gap width. As x → 0, the width of the gap goes to zero. We compute secular coefficients in such a gapped disc denoted and , by appropriately changing the upper integration limits in the definitions (18)- (19), i.e. from 1 to α ≡ 1 − x. This eliminates gravitational effect of the disc annuli with a(x) < a < a(x).

In Figure 2 we display the behavior of (Panel A) and (Panel B) as a function of for various values of softening ς 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 ψ1 depends on p, while ψ2 depends on p + q; Eqs. 18, 19). There are several notable features in this figure.

Figure 2

Behavior of the cumulative pre-factors (panel A) and (panel B) of the disturbing function due to a power-law disc (p = 1, q = 0.5 and ain → 0, aout → ∞) with softened gravity, shown as a function of x — relative separation between a given test-particle orbit and the nearest neighboring disc rings. Formalisms of Hahn (2003), Touma (2002), Tremaine (1998) and Teyssandier & Ogilvie (2016) are shown by different colors as indicated in panel (A), for different values of softening (shown by different line types). The purple lines represent results obtained by the unsoftened expressions of Davydenkova & Rafikov (2018) (DR18) based on the Heppenheimer method (see §6.3). Insets illustrate the behavior as x → 0 for the three convergent softened formalisms — see text (§3.2) for more details.

First, when the gap is wider than the characteristic softening length ςa, i.e. ς ≲ x ≤ 1, the amplitudes of both and increase from zero at x = 1 (infinitely wide gap) to their maximum values reached at x ∼ ς. In all cases ψ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 ς ≲ x ≪ 1 we find that ∼ x−1, irrespective of the softening model used; their maximum values are always ∼ ς−1. This convergent behavior is easy to understand since for ς ≲ x the role of softening is negligible, , and all ϕ effectively reduce to their classical counterparts given by Eqs. (8) - (9), which can be easily verified using the expressions listed in Table 1. The scaling of and with x is simply a result of asymptotic behavior of → (1 − α)−2 as α → 1, upon radial integration in Eqs. (18) – (19).

Second, upon reaching their extrema at x ∼ ς, amplitudes of and computed using softening prescrip-tions of Tr98, T02 and H03 start decreasing as x decreases. In the range of semi-major axes corresponding to x ≲ ς, softening significantly modifies the behavior of away from the divergent behavior of . The modification is such that the softened interaction with the disc annuli ≲ ςa away from the test-particle orbit starts to dy-namically counteract the contribution of the more distant annuli (with x ≈ 1). As a result of this compensation, and cross zero and change sign at some x = Cς2, where C ∼ 1 is a constant.

At the same time, and calculated according to Teyssandier & Ogilvie (2016) clearly show different behavior. Instead of decreasing in amplitude as x ≲ ς, they remain essentially constant, having reached their saturated values ∼ ς−1 at x ∼ ς. This explains the lack of convergence with S obvious in Figure 1C, D, since the values to which and converge keeps increasing as ς → 0. Moreover, both coefficients also never change sign, always predicting prograde precession (). The origin of this difference with other smoothing prescriptions will be addressed in §6.2.

Upon further decrease of x below ς2, both and 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 ς.

We note that the opposite contributions to e.g. ψ1 produced by the distant (x ≳ ς, positive) and nearby (i.e. with x ≲ ς, 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 ≲ e), with the disc annuli crossing the particle orbit giving rise to a negative contribution to . Eventually, the shrinking of the gap brings to a finite negative value (for p = 1 disc) as x → 0. This sequence is very similar to the behavior we find with softened gravity for x ≲ ς.

In Figure 3 we show calculations for 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, ∼ x−1 scaling for ς ≲ x ≪ 1, decay of for ς2 ≲ x ≲ ς, and ultimate convergence to ψ1 in a disc with no gap. The only obvious difference is the fact that does not cross zero for this disc model with p = −0.5.

Figure 3

Same as Figure 2, but now for an axisymmetric power-law disc with p = −0.5. Note that for this disc model softened does not cross zero and converges to a positive value as x → 0, in agreement with the results in Figure 1A.

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 ≲ ς) and distant disc annuli (i.e. with x ≳ ς), 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 ≲ ς), see §6.2. We will discuss additional implications of these calculations in §6.3.

Variation of disc model — p and q

We now examine the dependence of ψ1 and ψ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 prescriptions assuming a softening value of ς = 10−3 (for which Fig. 1A, B suggests good convergence of ψ1 and ψ2). For reference, black open circles show the expected behavior of ψ1 and ψ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 ψ1 and ψ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 ψ1 and ψ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 ψ1 equal to 1/2π and a negative contribution to ψ2 equal to −1/π (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 ≈ 0 or p ≈ 3 (for which SR15 yields ψ1 ≈ 0), and (2) a mismatch of tens of per cent between the disc-driven forced eccentricity oscillations, = |ψ2/ψ1|, and the expectations based on SR15. The latter point is illustrated in Figure 4C.

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 Σ smoothly goes to zero at finite radii in a manner that does not give rise to the edge effects, see DR18 and §5.

Quartic Disc Model

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

and linear eccentricity profile in the form

for ain ≤ a ≤ aout (with ain = 0.1 AU, aout = 5 AU), where = 1153 g cm−2 and = 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 ς = 10−3. It shows the excellent agreement between the radial profiles of A, B and 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 and B computed using this model from Davydenkova & Rafikov (2018), in particular, modify the locations at which A and B become zero. This explains the slight shift in the semi-major axes at which = 2B/A goes through zero or diverges, see Figure 5.

Figure 5

Performance of different softening formalisms (different colors) with softening parameter ς = 10−3 in the potential of a Quartic disc, see Eq. (23), with the eccentricity profile (24). The disc extends from ain = 0.1 AU to aout = 5 AU. Shown as a function of semi-major axis a are the profiles of (A) the amplitude of the disc-induced eccentricity oscillations, (B) the rate of disc-driven free precession A, and (C) the coefficient B appearing in the non-axisymmetric part of the disturbing function (4). The black lines represent the expected unsoftened results as computed by Davydenkova & Rafikov (2018). Curves for Hahn (2003) and Touma (2002) fall on top of the unsoftened results, while the softening method of Tremaine (1998) shows only qualitative agreement.

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 ∝ r and b = const respectively.

Gaussian Rings

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

centered around a = 1.5 AU with width w = 0.18 and surface density = 100 g cm−2 at ac . The eccentricity profile is still given by Eq. (24).

In Figure 6 we plot the behavior of the corresponding A, B and for the three (convergent) softened formalisms with ς = 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.

Figure 6

Same as Figure 5, but now for a Gaussian disc with Σ(a) and e(a) given by Eq. (25) and (24) respectively. Note that for this disc model the formalism of Tremaine (1998) (green) shows quite good agreement with the unsoftened results, even at the quantitative level. See text (§4.2) for details.

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 = aedge. To that effect we analyze the behavior of secular coefficient A 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 ς → 0) for different values of softening (results for B are very similar) near the disc edge. Figure 7 shows the run of A near the inner edge ain of the disc for particles both inside (a < ain) and outside (a > ain) the disc as predicted by the formalism of Hahn (2003). The calculation assumes circular PL disc with p = 1 and Σ0 = 100 g cm−2 extending between ain = 1 AU to aout = 10 AU, where we have set a0 = aout (Eq. 16).

Figure 7

The behavior of the free precession rate A near the inner edge ain = 1 AU of a circular power-law disc with surface density Σ(a) = 100 g cm−2 (10 AU/a) (Eq. 16). One can see that the expected divergent behavior of A near the disc edge is reproduced by the softening prescription of Hahn (2003) in the A limit ς → 0. However, very near the sharp edge of the disc ς has to be very small for quantitative accuracy to be attained. Similar results can be obtained by the softened formalisms of both Touma (2002) and Tremaine (1998).

The unsoftened calculations based on Heppenheimer (1980) invariably predict that the free eccentricity precession rate A, as well as B, 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 − ain|/ain ∼ 10−3) the agreement is achieved only for ς ≤ 10−4, which is considerably smaller than the values (ς ∼ 10−2) required to reproduce the dynamics of particles far from the disc edges, ain ≪ a ≪ aout, see Fig. 1. For ς = 10−2 the softened calculation predicts A 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 softening. 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 passes through zero exactly at ain, showing two sharp peaks of opposite signs just around this radius. Similar behavior was found by Davydenkova & Rafikov (2018) for zero-thickness discs with Σ 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 curve connecting the two peaks turns into a vertical line at ain.

Similar divergent behavior of A (and B) arises also at the outer edge of the disc considered in Fig. 7 and, in general, at any radius within a disc where Σ(a) exhibits a discontinuity.

Finally, we note that the dynamics of particles orbiting outside the disc (where Σ(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 = 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.

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 ϕ (see Table 1), and yet their results are consistent with the un-softened calculations as ς → 0.

To understand this variation of outcomes, we developed a general framework for computing secular coefficients ϕ (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 A4, which assumes the disc to be a continuous entity from the start. Both of them arrive at the same expressions for R.

Using these results we show in Appendix B that the expressions for ϕ found by Touma (2002) and Hahn (2003) can be recovered exactly using this general framework if we set ℱ (r1, r2) = and ℱ (r1, r2) = H2, 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.

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 defined by Eq. (12) suggests interaction potential (10) with ℱ (r1, r2) = for the softening model of Tremaine (1998). In Appendix B we show that propagating this form of ℱ(r1, r2) through our general framework results in the following expressions for the coefficients ϕ:

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 δ-function. Such terms arise because the form of ℱ(r1, r2) adopted in Tremaine (1998) is not sufficiently smooth — its first derivative is discontinuous at r1 = r2, while the calculation of ϕ involves second-order derivatives of ℱ, see Eqs. (A25)-(A27), as well as Eq. (A28). Such singular terms do not arise in other types of softening prescriptions examined in our work since they all use infinitely differentiable versions of ℱ(r1, r2). 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 ℱ used in Tremaine (1998).

Presence of these terms in Eqs. (26)-(27) introduces corrections to coefficients A and B (Eqs. 5, 6) in apse-aligned discs in the form

Accounting for these corrections, we confirmed that the correct (un-softened) behavior of the coefficients of R can be reproduced for the non-PL discs – Quartic and Gaussian models, see §4. Note that δ A(a) and δB(a) are proportional to the local disc surface density £d(a) and (α = 1) ∼ , see Eq. (C7). 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 δ-function in Eqs. (26)-(27) give rise to corresponding modifications of the coefficients ψ1 and ψ2 defined by Eqs. (18)-(19):

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.

Figure 8

Similar to Figure 1, but now using the expressions for ϕ given by Eqs. (26-27) and Eqs. (32-33) obtained by propagating of Tremaine (1998) and ℱ(r1, r2) = ς r1r2 of Teyssandier & Ogilvie (2016), respectively, through the general framework outlined in Appendix A. Shown as a function of softening ς are ψ1 (panel A) and ψ2 (panel B) for two PL disc models specified by p and q indicated in panel A. Black lines represent the expectations based on Silsbee & Rafikov (2015), to which the new expressions for ψ1 and ψ2 successfully converge as ς → 0.

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 α ≠ 1 and terms in Eqs. (26)-(27) containing δ-function naturally vanish, reducing ψ1 and ψ2 back to the expressions quoted in Tremaine (1998).

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 in Eq. (15) one infers that this model features softening parameter in the form ∈2(α) = S2α. To soften secular interaction Teyssandier & Ogilvie (2016) directly disc result even with a relatively coarse radial sampling of substituted in the classical expressions (8, 9) for with , 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 ∈2(α) = ς2α corresponds to softening function ℱ(r1, r2) = ς2r in the two-point potential (10), see Eq. (A21). Propagating such a form of ℱ(r1, r2) through our general framework in Appendix A, we find the following expressions for the coefficients ϕ with ς = S (Appendix B):

Approach of Teyssandier & Ogilvie (2016) accounts for only the first terms in Eqs. (32), (33), with coefficients which are O(S0), see Table 1. However, as we show below, the correct behavior of ϕ as S → 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 ϕ given by Eqs. (32) – (33). One can see see that the correct implementation of the softening ∈2(α) = S2α 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 ψ1 and ψ2 smoothly converge to their expected unsoftened values in the limit of S → 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 as done, for instance, in Appendix A. In other words, a direct replacement of the classical Laplace coefficients 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.

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 → ∞.

Computational cost associated with the evaluation of mutual ring-ring interactions in this setup, going as O(N2), imposes limitations on the number of rings that can be used in practice. This is typically not a problem for the unsoftened 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 and computed without accounting for the contributions from a < a < a (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 and , 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 ∆a separating a given test-particle orbit

from nearest neighboring inner and outer disc rings should be quite small, ∆a/a ≲ 0.1ς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 (ςa), meaning that N has to be very large, N ≳ 10ς−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 Ad 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 (or ψ1 for PL discs) by using the following global metric

Here fnum(a) is the value of the metric basis (e.g. precession rate A) evaluated at the position a of ith ring by summing up the contributions of all other rings in the disc, while ftheor(a) is the analogous quantity computed in the limit of a continuous disc, i.e. as N → ∞ (it is given by the non-discretized version of Eq. (5) if f = A, or Eq. (18) if f = ψ1). Repeating this calculation for various combinations of (N, ς), we can determine the smallest number of rings N(ς) that ensures the desired convergence to within, e.g. ∼ 10% (i.e. ℳ(f) ∼ 0.1), for a given value of softening ς.

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 legend. Figure 9 shows that as ς → 0, the number of rings scales as N ∼ Cς−χ with C ∼ 10 and χ ≈ (1.8 – 1.9). The only notable exception is the Gaussian ring, for which convergence is faster (i.e. N ∝ ς−1.5), probably because of mass concentration in a narrow range of radii.

Figure 9

Scaling of number of softened annuli (rings) N with softening parameter ς to ensure convergence of disc-driven free precession A (or ψ1) in discretized discs to the expected results in continuous softened discs (Eqs. 5, 18). Calculations assume axisymmetric disc models extending from ain = 0.1 to aout = 5 AU: two PL discs (specified by p), a Quartic disc (same as Fig. 5) and a Gaussian ring (same as Fig. 6). We have used the softening methods of Hahn (2003), Tremaine (1998) and (corrected) Teyssandier & Ogilvie (2016), as specified in the panel. Convergence is measured using the metric ℳ(f) defined by Eq. (34). One can see that, when ς ≲ 0.1, N ∼ Cς−β, with C ∼ 10 and 1.5 ≲ χ ≲ 2. Similar results can be obtained for eccentric discs, and other softening prescriptions. See text (§6.3) for details.

We note that the proportionality constant C in the N(ς) relation is not perfectly defined in the sense that it depends on the (i) desired accuracy (roughly inversely proportional to ℳ(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 ∼ 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 ∼ ς−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 ∼ O(ς−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. ∆a/a ≲ 0.1ς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 ς ≲ 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.

Further generalizations and extensions

All calculations in this work are based on the expansion of the secular disturbing function R 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(a) = 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 to 4th order in eccentricities based on the un-softened method of Heppenheimer (1980). The general framework for calculating R with arbitrary softening prescriptions presented in Appendix A can also be extended to higher order in eccentricities in similar way, 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.

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 ς for all disc models.

The softening prescription proposed in Tremaine (1998) yields convergent results as ς → 0. However, quantitative differences of up to ∼ (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 ς → 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 ∼ Cς−χ, C ∼ O(10), 1.5 ≲ χ ≲ 2) is required to ensure that the correct secular behavior is properly captured by such calculations when ς 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.

Table B1

The functional forms of the coefficients T (α) given by Eqs. (A25)-(A27) appearing in the orbit-averaged disturbing function due to two coplanar (arbitrarily) softened rings (Eq. A22-A24) such that α ≡ a< /a> ≤ 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 ℱ(r1, r2) in Eq. (A1). The corresponding expressions for the dimensionless softening parameters (Eq. A21) entering in the definition of softened Laplace coefficients (Eq. A20) are also shown. Here, Θ(x) represents the Heaviside step function and δ(x) = dΘ(x)/dx stands for Dirac delta-function.

Methods	ℱ(r1, r2)	∈2(α)	T1(α)	T2(α)	T3(α)	T4(α)	T5(α)	T6(α)	T7(α)	T8(α)
H03	H^2r12+r22	H2(1 + α2)	0	2H2(1 + α2)	0	0	4αH2(2 + H2)	0	2H^2a_1a_2	2H^2a_2a_1
T02	bc2	β2 = (bc/a>)2	2β2	0	0	0	0	0	0	0
Tr98	βc2maxr12,r22	βc2	0	2βc2	−2βc2δα−1	0	4αβc2	0	2βc2αΘa_1−a_2	2βc2αΘa_2−a_1
T016	S2r1r2	S2α	0	2αS2	S2	0	S2(S2α + 2α2 + 2)	0	S2	S2