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 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 ≪ M c , particles perturbed by the disc gravity move on nearly-Keplerian orbits which evolve rather ⋆ E-mail: aas79@cam.ac.uk 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 ̟ p due to a co-planar point mass δm d orbiting with semi-major axis a, eccentricity e d , and apsidal angle ̟ d , upon smearing into elliptical rings, the secular disturbing function takes the form (Murray & Dermott 1999) valid for a > a p as well as a < a p , as long as particle orbits do not cross. Here b which obeys b (m) s (α). Explicit time independence of δ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 δR (Eq. 1) over the radial extent of the disc: where a in and a out 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 semimajor axes of the rings) is needed, significantly accelerating calculations 1 . Unfortunately, this straightforward procedure is illposed from the mathematical point of view. Indeed, it is well known that the Laplace coefficients b (m) 3/2 featured in Eq.
(1) diverge as b (m) 3/2 (α) → (1 − α) −2 when α → 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 in ≤ a p ≤ a 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 onedimensional 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 1998Tremaine , 2001Touma 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 → 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.

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 µ d (a), eccentricity e d (a), and apsidal angle ̟ 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 Σ d (a) instead of µ d (a); its relation to µ d for arbitrary profiles of e d and ̟ 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 µ d (a) does not change in time. The same statement is true for Σ d (a) to lowest order in (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 semimajor axis a p , eccentricity e p , and apsidal angle ̟ p : Here n p is the test-particle mean motion (n 2 p = GM c /a 3 p ), and we have introduced a two-component eccentricity vector for a test particle e p = e p (cos ̟ p , sin ̟ p ).
The coefficients A d and B d in Eq. (4) are related to the disc mass (or surface density) and eccentricity profiles in the following fashion: where e d = e d (a)(cos ̟ d (a), sin ̟ d (a)) is the eccentricity vector for an annular disc element 2 . Functions φ ij (α), 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 φ ij appearing in the literature are linear combi-2 We refer the reader to Heppenheimer (1980); Silsbee & Rafikov (2015); Davydenkova & Rafikov (2018) The softening parameter ǫ(α) appearing in this definition remains non-zero as α → 1, thus preventing the divergence of the softened Laplace coefficients B (m) s (α)). 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 relations for softened Laplace coefficients B (m) s (α, ǫ), 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, µ d (a) → m pl δ(a − a pl ) (where δ(z) is the Dirac delta-function), and set softening to zero (so that B (m) s (α)). Then one finds that R d reduces to the un-softened, orbit-averaged potential δR due to a planet with mass m pl and semi-major axis a pl , with the unsoftened coefficients φ ij 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 p .
We emphasize that the functional forms of φ ij are not s appearing in the unsoftened definition (8) -(9) by B (m) s . This can be seen in Table 1 where we summarize some of the expressions for φ ij (α) proposed in the literature and analyzed in this paper (see §2.1). Nevertheless, examination of these expressions shows that when ǫ 2 (α) → 0, the coefficients φ ij (α) do reduce to their unsoftened versions φ LL ij (α) 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 ≪ M c ) located at r 1 and r 2 , correspondingly, relative to the central mass, has a form 3 with i, j = 1, 2 and j i. Here 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 r 1 −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 standalone 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 F via ǫ 2 = [max(a 1 , a 2 )] −2 F (a 1 , a 2 ), where a 1,2 are the semi-major axes of the interacting particles (see Eq. A21). The most general expressions of φ ij entering the arbitrarily softened ring-ring disturbing function, (here i = 1, 2 and j i) is given by Eqs. (A22)-(A24) in terms of B (m) s (α, F ). In the above expression, we have defined a > = max(a 1 , a 2 ) and a < = min(a 1 , a 2 ) such that 4 α = 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 F the coefficients φ ij (α) do not obey certain symmetry properties when a/a p is replaced with a p /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 d and B d , we remind the reader that A d represents the precession rate of the free eccentricity vector of a test particle in the disc potential, while B d characterizes the torque exerted on the particle orbit by the non-axisymmetric component of the disc gravity. Corresponding forced eccentricity vector is e p, f = −B d /A d . In particular, test-particles initiated on circular orbits experience eccentricity oscillations of maximum amplitude e m p = 2 e p, f . As A d (a p ) and B d (a p ) uniquely determine R d 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 φ ij (α) are provided in Table 1.
Here β 2 c is the dimensionless softening parameter, treated as a constant, i.e. independent of distance. The physical interpretation of this manoeuvre is that β c , inhibiting the formal divergence of R d as a → 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 distancedependent radius b = β c max(a 1 , a 2 ). In Tremaine (1998) coefficients φ ij (α) were expressed as derivatives of B (m),Tr 1/2 with respect to α, see equations (26) of Tremaine (1998). These expressions, along with their versions modified using the recursive relations for Laplace coefficients (see Appendix C1), can be found in Table 1.

Formalism of Touma (2002) -T02
Touma (2002) However, in Touma (2002) the softening parameter ǫ 2 (α) = β 2 is no longer a constant but depends on the distance such that β = 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 .

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 ∼ h/a, which was assumed to be constant in that work (see also Ward 1989). Hahn (2003) demonstrated that the resultant φ ij (α) are functions of softened Laplace coefficients with constant H ≪ 1. In other words, the softening parameter is given by ǫ 2 (α) = H 2 (1 + α 2 ) in that work. The explicit expressions for φ ij (α) in terms of B Thus, their softening parameter is ǫ 2 (α) = S 2 α, 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 φ ij (α) are given by equations (7)-(9) of Teyssandier & Ogilvie (2016).
The aforementioned softening prescriptions have their softening parameters ǫ 2 (α) controlled by different constants β 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 ς.

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 particle 5 . 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 (α) ∝ log(1 − α). 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. 5 Note 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).

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̟ d /da = 0), truncated power-law (hereafter PL) disc as a simple example. We characterize surface density and eccentricity of such a disc by for a in ≤ a ≤ a out , where Σ 0 and e 0 are the pericentric surface density and eccentricity of the disc at some reference semimajor 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) where K = πGΣ 0 a p 0 a 1−p p and the dimensionless coefficients ψ 1 and ψ 2 are given by with α 1 = a in /a p and α 2 = a p /a out . The coefficients ψ 1 and ψ 2 are functions of the powerlaw indices (p and q), any softening parameter involved (through φ ij ), as well as the test-particle semi-major axis a p (through α 1,2 ). They are related to A d and B d via As shown in Appendix D, for certain ranges of powerlaw indices p and q both ψ 1 and ψ 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 α 1,2 → 0). For p and q in 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 Σ d (a) and e d (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 p (i.e. functions of p, q, and ς only), making them useful as simple metrics for judging the validity of different models of softening. , 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.

Behavior with respect to variation of softening
"softening" 6 ς 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 the un-softened Heppenheimer approach 7 . 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 8 find that ς 10 −3 guar-antees the convergence of ψ 1 and ψ 2 to within few per cent of the correct values for all p and q as long as a in ≪ a p ≪ a out (see Figure 4).
The same panels also indicate that ψ 1 (ς) and ψ 2 (ς) predicted by the softened formalism of Tremaine (1998) (green), while converging to finite values as ς = β c → 0, do not reproduce the SR15 results exactly in this limit. Indeed, one can see that even for the smallest adopted value of β c = 10 −3 , the softening prescription of Tremaine (1998) yields ψ 1 and ψ 2 different by tens of per cent from SR15. It is easy to demonstrate that these quantitative differences do not vanish by further decreasing β c . 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 β c → 0, the ψ 1 and ψ 2 curves computed using softening 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 Distance separating a test-particle and its neighboring disc rings 0 ∞ 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 p = 1 A. x -1 10 -8 10 -7 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 Distance separating a test-particle and its neighboring disc rings 0 ∞ x = 1m 10 -8 10 -7 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 -10 3 -10 2 -10 1 -10 0 -10 -1 -10 -2 -10 -3 -x -1 Figure 2. Behavior of the cumulative pre-factorsψ 1 (x) (panel A) andψ 2 (x) (panel B) of the disturbing function due to a power-law disc (p = 1, q = 0.5 and a in → 0, a out → ∞) 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)  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 predicts ψ 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 B (m),TO 3/2 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 p or α → 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 ψ 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 Σ d (a) is set to zero, are at a d,i = (1 − x)a p ≤ a p and a d,o = (1− x) −1 a p ≥ a p , respectively, with a single parameter x 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ψ 1 (x) andψ 2 (x), by appropriately changing the upper integration limits in the definitions (18)-(19), i.e. from 1 to α m ≡ 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ψ 1 (x) (Panel A) 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.
First, when the gap is wider than the characteristic softening length ςa p , i.e. ς x ≤ 1, the amplitudes of bothψ 1 (x) andψ 2 (x) 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ψ 1 (x) ∼ |ψ 2 (x)| ∼ 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, B (m) , and all φ ij effectively reduce to their classical counterparts φ LL ij given by Eqs. (8) -(9), which can be easily verified using the expressions listed in Table 1. The scaling ofψ 1 (x) and |ψ 2 (x)| with x is simply a result of asymptotic behavior of b Second, upon reaching their extrema at x ∼ ς, amplitudes ofψ 1 (x) andψ 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 ς, softening significantly modifies the behavior of B s (α). The modification is such that the softened interaction with the disc annuli ςa p away from the test-particle orbit starts to dynamically counteract the contribution of the more distant annuli (with x ≈ 1). As a result of this compensation,ψ 1 andψ 2 cross zero and change sign at some x = Cς 2 , where C ∼ 1 is a constant 9 .
At the same time,ψ TO16 1 andψ TO16 2 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 |ψ TO16 1 | and |ψ TO16 2 | converge keeps increasing as ς → 0. Moreover, both coefficients also never change sign, always predicting prograde precession (ψ TO16 1 > 0). The origin of this difference with other smoothing prescriptions will be addressed in §6.2.
Upon further decrease of x below ς 2 , bothψ 1 andψ 2 computed using models of Tr98, T02 and H03 ultimately converge to their corresponding values obtained for a con-9 For p = 1,ψ 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α in Eq. (18) -(19), we find C T r 98 = (π − 1)/2 and C H 03 = π; in agreement with Fig. 2. For other values of p and q, for which ψ 1 < 0 (c.f. Fig. 4), we numerically find that C varies by at most a factor of ten.
tinuous 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 testparticle orbit starts crossing the gap edge (i.e. for x e p ), with the disc annuli crossing the particle orbit giving rise to a negative contribution toψ 1 . Eventually, the shrinking of the gap bringsψ 1 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ψ 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,ψ 1 (x) ∼ x −1 scaling for ς x ≪ 1, decay ofψ 1 (x) for ς 2 x ς, and ultimate convergence to ψ 1 in a disc with no gap. The only obvious difference is the fact thatψ 1 does not cross zero 10 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 ς) 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 modelp 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 11 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, 10 This 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  Distance separating a test-particle and its neighboring disc rings 0 ∞ x -1 p = -0.5 H03 T02 Tr98 TO16 DR18 Softening 10 -3 10 -2 10 -1 x = 1m Figure 3. Same as Figure 2, but now for an axisymmetric powerlaw disc with p = −0.5. Note that for this disc model softenedψ 1 (x) does not cross zero and converges to a positive value as x → 0, in agreement with the results in Figure 1A.
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, e m p /e d (a) = |ψ 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 Σ d 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

C.
A.

DR18 H03 T02 Tr98
Quartic 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 a in = 0.1 AU to a out = 5 AU. Shown as a function of semi-major axis a p are the profiles of (A) the amplitude e m p of the disc-induced eccentricity oscillations, (B) the rate of disc-driven free precession A d , and (C) the coefficient B d 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. and linear eccentricity profile in the form for a in ≤ a ≤ a out (with a in = 0.1 AU, a out = 5 AU), wherẽ Σ 0 = 1153 g cm −2 andẽ 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 ς = 10 −3 . It shows the excellent agreement between  Figure 6. Same as Figure 5, but now for a Gaussian disc with Σ d (a) and e d (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. the radial profiles of A d , B d and e m p 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 m p = 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)  . The behavior of the free precession rate A d near the inner edge a in = 1 AU of a circular power-law disc with surface density Σ d (a) = 100 g cm −2 (10 AU/a) (Eq. 16). One can see that the expected divergent behavior of A d near the disc edge is reproduced by the softening prescription of Hahn (2003) in the 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 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 c = 1.5 AU with width w c = 0.18 and surface densityΣ 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 m p 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.

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 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 ς → 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 in of the disc for particles both inside (a p < a in ) and outside (a p > a in ) 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 a in = 1 AU to a out = 10 AU, where we have set a 0 = a 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 in |/a in ∼ 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, a in ≪ a p ≪ a out , see Fig. 1. For ς = 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 softening 12 . 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 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 Σ 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 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 Σ d (a) exhibits a discontinuity.
Finally, we note that the dynamics of particles orbiting outside the disc (where Σ 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.

DISCUSSION
Results of previous sections reveal a diversity of outcomes when different softening models are applied. Two modelsthose 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) 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 φ 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 orbitaveraging 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 d .
Using these results we show in Appendix B that the expressions for φ ij found by Touma (2002) and Hahn (2003) can be recovered exactly using this general framework if we set F (r 1 , r 2 ) = b 2 c and F (r 1 , r 2 ) = H 2 (r 2 1 + 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 B (m),Tr s defined by Eq. (12) suggests interaction potential (10) with F (r 1 , r 2 ) = β 2 c max(r 2 1 , r 2 2 ) for the softening model of Tremaine (1998). In Appendix B we show that propagating this form of F (r 1 , r 2 ) through our general framework results in the fol-lowing expressions for the coefficients φ ij : (1),Tr (2),Tr (1),Tr 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 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 φ ij involves second-order derivatives of F , 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 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 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 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 δ A d (a p ) and δB d (a p ) are proportional to the local disc surface density Σ d (a p ) and B (m),Tr 3/2 (α = 1) ∼ β −2 c , 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) Figure 1, but now using the expressions for φ i j given by Eqs. (26)(27) and Eqs. (32-33) obtained by propagating F(r 1 , r 2 ) = ς 2 max(r 2 1 , r 2 2 ) of Tremaine (1998) and F(r 1 , r 2 ) = ς 2 r 1 r 2 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.

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 B On the other hand, in Appendix B we show that softening parameter in the form ǫ 2 (α) = ς 2 α corresponds to softening function F (r 1 , r 2 ) = ς 2 r 1 r 2 in the two-point potential (10), see Eq. (A21). Propagating such a form of F (r 1 , r 2 ) through our general framework in Appendix A, we find the following expressions for the coefficients φ ij with ς = S (Ap-pendix B): 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 φ ij 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 φ ij given by Eqs. (32) -(33). One can see see that the correct implementation of the softening ǫ 2 (α) = S 2 α 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 13 as done, for instance, in Appendix A. In other words, a direct replacement of the classical Laplace coefficients b (m) 3/2 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(N 2 ), 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 13 An 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 (m) 1/2 with B (m) 1/2 in Eq. (7.1) of Murray & Dermott (1999) before applying the derivatives with respect to α. We note that this procedure could apply for all F(r 1 , r 2 ) with continuous first derivatives satisfying D 1 + D 2 = −1; see Appendix A.
the integral contribution to e.g. the precession rate. Indeed, purple curves in Figures 2 & 3 demonstrate this by showing the un-softenedψ 1 (x) andψ 2 (x) computed without accounting 14 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ψ 1 (x) andψ 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 ∆a separating a given test-particle orbit from nearest neighboring inner and outer disc rings should be quite small, ∆a/a p 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 p ), 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 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 ψ 1 for PL discs) by using the following global metric 15 Here f num (a i ) is the value of the metric basis (e.g. precession rate A d ) evaluated at the position a i of ith ring by summing up the contributions of all other rings in the disc, while f theor (a i ) 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 d , 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. M( 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)  . Scaling of number of softened annuli (rings) N with softening parameter ς to ensure convergence of disc-driven free precession A d (or ψ 1 ) in discretized discs to the expected results in continuous softened discs (Eqs. 5, 18). Calculations assume axisymmetric disc models extending from a in = 0.1 to a out = 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 M( 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. ous axisymmetric disc models as indicated in the legend 16 . Figure 9 shows that as ς → 0, the number of rings scales as N ∼ Cς −χ with 17 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.
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 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 ∼ 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 16 We exclude the softening method of Touma (2002) from this analysis as it introduces additional complexity due to the nature of softening parameter; ǫ 2 = b 2 /max(a 2 1 , a 2 2 ), see §2.1.2. 17 For example, the curve computed using the (corrected) model of Teyssandier & Ogilvie (2016) has C = 10.9 and χ = 1.91, while the one for Quartic disc has C = 7.2 and χ = 1.75. 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 p 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 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 4th 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 way 18 , 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 razorthin 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 18 Another way to calculate the softened disturbing function for arbitrarily high eccentricities is to numerically compute the ringring interaction potential, as was done by Touma et al. (2009). 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 interparticle 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 1. The coefficients φ i j (α) 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 < /a > where a > = max(a 1 , a 2 ) and a < = 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 φ LL i j 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 φ i j reported in Touma (2002) have been corrected in a subsequent paper of Touma & Sridhar (2012). (1),Tr (2),Tr (1),Tr 5/2

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 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 as 19 where F (r 1 , r 2 ) is an arbitrary softening function introduced to cushion the singularity which arises otherwise at null interparticle separations. In the above expression, f i is the true anomaly of the i th ring, ̟ 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 Ψ to second order in eccentricities e i valid for arbitrary F (r 1 , r 2 ).

A1 Expansion of the interaction function Ψ around small eccentricities
Following the classical techniques of celestial mechanics (see, Plummer 1918, Ch. XVI), we start by expanding Ψ around circular orbits. Using Taylor expansion we write where θ = M 1 − M 2 + ̟ 1 − ̟ 2 , M i represents the mean anomaly of the i th ring characterized with semi-major axis a i , and the linear operators D k are given by (Plummer 1918) 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 T defined by Eq. (A2) on the disturbing function of circular softened rings Ψ 0 , we make use of the elliptical expansions of r/a and f − M, to multiply individual terms appearing in 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 Ψ reduces to where the operators A, B and C acting on Ψ 0 are defined as We have used the fact that cos(M 1 − M 2 ) = cos θ cos(̟ 1 − ̟ 2 ) and sin(M 1 − M 2 ) = sin θ cos(̟ 1 − ̟ 2 ) in the secular regime (Plummer 1918).
Equipped with the expression (A7) for Ψ, we proceed to compute the action of operator T on Ψ 0 prior to orbit-averaging the resultant expression. With this in mind, we compute the action of several operators appearing in the definitions of A, B and C on Ψ 0 and list them below: D 1 D 4 Ψ 0 = a 1 a 2 sin θ Ψ 3 0 − 3a 1 a 2 sin θ a 2 1 − a 1 a 2 cos θ + where for conciseness we have written F instead of F (a 1 , a 2 ). Here, it is worthwhile to mention that, as far as the expansion technique is concerned, the terms ∂ i F (with i = 1, 2) appearing in the above expressions are the only difference brought upon by softening the Newtonian point-mass interaction (Eq. A1). Another set of operators useful in computing TΨ 0 is the following: which can be obtained by making use of the identity (D 1 + D 2 + 1)Ψ 0 = 1 2 (2F − a 1 ∂ 1 F − a 2 ∂ 2 F )Ψ 3 0 . Here, we note that for all softening functions F for which 2F − a 1 ∂ 1 F − a 2 ∂ 2 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 AΨ 0 = BΨ 0 (since D 2 3 = D 2 4 , see Eqs. (A8) and (A10)). As a result, the resultant orbit-averaged disturbing function (A7) 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 F (r 1 , r 2 ) = const 0.

A3 Orbit-averaging the interaction function Ψ
The expressions (A10)-(A17) allow the computation of Ψ = TΨ 0 , which needs to be time-averaged in order to recover the secular disturbing function. We do not show the cumbersome collated expression for TΨ 0 and proceed to the final step of orbit-averaging, which will conclude our derivation. In short, our goal is to compute which essentially reduces to computing the individual terms AΨ 0 , BΨ 0 and CΨ 0 . At the outset, it is important to note that each of the terms appearing in TΨ 0 (through AΨ 0 , BΨ 0 and CΨ 0 , or the operators they entail) are proportional to cos(mθ)Ψ 2s 0 . By making use of α = a < /a > , where a < = min(a 1 , a 2 ) and a > = max(a 1 , a 2 ), this combination can be reduced to For that reason, calculation of the orbit-averaged Ψ (by integrating over dθ) yields integrals of the form which is the generalization of the classical Laplace coefficients b (m) s (recovered when F (a 1 , a 2 ) = 0, see Eq. 2) with the dimensionless softening parameter see Eq. (7). Employing this notation, we present the simplified expressions of AΨ 0 , BΨ 0 and CΨ 0 obtained as a result of orbit-averaging: In equations (A22)-(A24), we have defined the dimensionless functions T i (α) such that where, as before, F ≡ F (a 1 , a 2 ), α = a < /a > and ∂ i ≡ ∂/∂a i . Note that the expressions for φ 11 and φ 22 swap definitions upon replacing a 1 by a 2 , whilst keeping α < 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 ⇌ a 2 while, at the same time, T 4 and T 7 (appearing in the second line of Eq. (A22)) translate to T 6 and T 8 (appearing in the second line of Eq. (A23)); and vice versa. These identities, when combined, yield the desired expression of Ψ = TΨ 0 ; see Eqs. (A7)-(A9). Subsequently, the softened ring-ring disturbing function in the form (11) is recovered, with the coefficients φ ij defined by Eqs. (A22) -(A24). 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 F (r 1 , r 2 ).
Note that in the absence of softening (i.e. F (r 1 , r 2 ) = 0) T i = 0 for all i and the classical expressions for φ LL 11 , φ LL 22 and φ LL 12 -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.
A4 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 = G ∫ S Σ(r d )Φ(r d , r p )dS , where Φ is the interaction potential given by equation (10), angle brackets indicate averaging over the orbit of the test particle given by r p and integration is carried out over the full surface of the disc S with 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 Φ(r d , 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 F on r d and 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: 2n p a 2 p ∫ aΣ(a)da a > 1 4 3a p F ′ F ′ + 4a p − 2 2F ′ + a p F ′′ a 2 p + a 2 + F − 12a p F a p B where prime denotes differentiation with respect to a p (e.g. F ′ = ∂F /∂a p ), a > = max(a p , a), α = 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. B d ) is somewhat more tedious but can nevertheless be done similar to Davydenkova & Rafikov (2018).

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 out /a in ≫ 1. In such situations, far from the disc edges one can take the limit of both α 1 = a in /a p and α 2 = a p /a 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 ψ 1 and ψ 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 φ ij (α), which appear in the integrands of each of ψ 1 and ψ 2 , in the limit of α ≈ 0. Using the Taylor expansions of softened Laplace coefficients B (m) s , we determined that both ψ 1 and ψ 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, β c ). This follows from the fact that for both Hahn (2003) and Tremaine (1998) we have φ 11 = φ 22 ∼ α 2 and φ 12 ∼ α 3 to lowest order in α. 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 φ 11 = φ 22 ≈ − 1 4 S 2 α + 3 8 (1+ 3 2 S 4 )α 2 and φ 12 ≈ 3 2 S 2 α 2 − 15 16 (1+5S 4 )α 3 which, in the limit of S → 0, translate to the same ranges for ψ 1 and ψ 2 convergence as Silsbee & Rafikov (2015). However, when S is relatively large, it is trivial to show that ψ 1 and ψ 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 ψ 1 and ψ 2 convergence are in line with the findings of Silsbee & Rafikov (2015) when the corresponding softening parameter b c → 0. However, when b c is non-zero, the ranges are narrowed down to −1 < p < 2 and −2 < p + q < 3 respectively. This paper has been typeset from a T E X/L A T E X file prepared by the author.