==== Front Mon Not R Astron Soc Mon Not R Astron Soc MNRAS Monthly Notices of the Royal Astronomical Society 0035-8711 1365-2966 Oxford University Press MN-489-03-4176 10.1093/mnras/stz2412 Article Potential softening and eccentricity dynamics in razor-thin, nearly-Keplerian discs Softened potentials of discsSefilian & RafikovSefilian Antranik A. 1★ Rafikov Roman R. 12 1 Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK 2 Institute of Advanced Study, Einstein Drive, Princeton, NJ 08540, USA ★ E-mail: aas79@cam.ac.uk 02 9 2019 02 9 2019 489 3 4176 4195 © 2019 The Authors2019This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.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. celestial mechanicsmethods: analyticalplanet-disc interactionsplanets and satellites: rings ==== Body 1 INTRODUCTION Astrophysical discs orbiting a central mass Mc are ubiquitous in a variety of contexts – galactic, stellar, and planetary (Latter et al. 2017). In many instances, masses of such discs Md 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 Md ≪ Mc, 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 Rd, fully determines the secular dynamics of the system. For a test particle with semi-major axis ap, eccentricity ep, and apsidal angle ω¯p due to a co-planar point mass δmd orbiting with semi-major axis a, eccentricity ed, and apsidal angle ω¯d, upon smearing into elliptical rings, the secular dis-turbing function takes the form (Murray & Dermott 1999) δR=Gδmdapa218b3/21apaep2−14b3/22apaepedcosω¯p−ω¯d,ϖϖ(1) valid for a > ap as well as a < ap, as long as particle orbits do not cross. Here bsmα is the Laplace coefficient defined by bsmα=2π∫0πcosmθ1+α2−2αcosθ−sdθ,(2) which obeys bsmα−1=α2sbsmα. 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: Rd=∫ainaoutδR,(3) 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 calculations1. Unfortunately, this straightforward procedure is illposed from the mathematical point of view. Indeed, it is well known that the Laplace coefficients b3/2m featured in Eq.(1) diverge as b3/2mα→1−α−2 when α → 1. This implies that the radial integration in Eq. (3) encounters an essential singularity at a = ap . As a result, for a co-planar particle orbiting inside a razor-thin disc, ain ≤ ap ≤ aout, this direct way of computing Rd does not converge to a finite value. This divergence, as well as the pressing need for having an efficient way of computing Rd (via a one-dimensional integration over a only), have motivated the development of alternative analytic approaches for calculating Rd . 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 Rd as a → ap . In this method orbit crossing does not lead to problems as long as the softening scale is finite. However, a physical justification for a specific form of softening (absent in the Heppenheimer (1980) approach) often remains unclear, making the introduction of softening rather arbitrary. The primary goal of our present work is to assess how well the different calculations relying on potential softening reproduce secular dynamics driven by the gravity of a razor-thin disc. The main metric we use in this exercise is the convergence of the results of such calculations to the true secular evolution (represented by the un-softened Heppenheimer method) in the limit of vanishing softening, when the limit of Newtonian gravity is recovered. Complementary to this, we develop a general framework for computing the well-behaved secular disturbing function for a broad range of softened gravitational potentials. Our work is organized as follows. We describe the general analytical expressions governing the orbit-averaged potential due to a coplanar disc of arbitrary structure and arbitrary softening prescription in §2. Having provided a brief account of the different softened potentials under our probe and the un-softened approach of Heppenheimer in §2.1 and §2.2, respectively, we analyze their performance in reproducing the correct secular dynamics for various disc models in §3, §4 and §5. We discuss and briefly summarize our results in §6 and §7 respectively. Technical details of our calculations can be found in Appendices. 2 DISTURBING FUNCTION DUE TO A DISK Prior to providing the details of different softening prescriptions examined in this work in §2.1, we briefly summarize some of their common features. The ultimate goal of all these prescriptions is the calculation of the disturbing function Rddue 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 ed(a), and apsidal angle ω¯da 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 ed 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 ed since µd(a) ≈ 2πaΣd(a) + O(ed) (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 ap, eccentricity ep, and apsidal angle ω¯p: Rd=npap212Αdapep2+Βdap⋅ep.(4) Here np is the test-particle mean motion (np2=GMc/ap3), and we have introduced a two-component eccentricity vector for a test particle ep=epcos ω¯p,sin ω¯p. The coefficients Ad and Bd in Eq. (4) are related to the disc mass (or surface density) and eccentricity profiles in the following fashion: Αdap=2Gnpap3× ∫ainapμdaϕ22aapda +∫apaoutμdaapaϕ11apada,(5) Βdap=Gnpap3× ∫ainapμdaedaϕ12aapda +∫apaoutμdaedaapaϕ12apada(6) where ed=edacos ω¯da,sin ω¯da is the eccentricity vec-tor for an annular disc element2. 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 combinations of softened Laplace coefficients Bsm defined by Bsmα,∈=2π∫0πcosmθ1+α2−2α cosθ+∈2α−sdθ(7) Table 1 The coefficients ϕij(α) 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 Bsmα,∈ 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 ϕij 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α+α2d2dα2B1/20,Tr=18αB3/21,Tr−3αβc2B5/20,Tr 142−2αddα−α2d2dα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αH22+H2B5/20,H −14αB3/22,H−3αH22+H2B5/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 Bsmα,∈ at α = 1 (unlike the classical bsmα. 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 Bsmα,∈, as well as their approximate asymptotic behavior and relationships to complete elliptic integrals. The mathematical structure of Rd 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) → mplδ(a − apl) (where δ(z) is the Dirac delta-function), and set softening to zero (so that Bsmα,∈→0→bsmα). Then one finds that Rd 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 ϕij in the form (Murray & Dermott 1999) ϕ11LLα=ϕ11LLα=18αb3/21α(8) ϕ11LLα=−14αb3/22α(9) 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 → ap. We emphasize that the functional forms of ϕij are not simple replacements of bsm appearing in the unsoftened definition (8) - (9) by Bsm. 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). Nev-ertheless, examination of these expressions shows that when ∈2(α) → 0, the coefficients ϕij(α) do reduce to their unsoft-ened versions ϕijLLα 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 (mi ≪ Mc) located at r1 and r2, correspondingly, relative to the central mass, has a form3 Φir1,r2=−Gmjr1,r22+Fr1,r2−1/2(10) 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 Rd 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 ϕij entering the arbitrarily softened ring-ring disturbing function, Ri=Gmja>ϕ11α e12+ϕ22α e22+ϕ12α e1e2cosω¯1−ω¯2,(11) (here i = 1, 2 and j ≠ i) is given by Eqs. (A22)-(A24) in terms of Bsmα,f. In the above expression, we have defined a> = max(a1, a2) and a< = min(a1, a2) such that 4 α = a. Note that in equations (5) and (6) we split integration over a in two parts: over the part of the disc interior to ap, and exterior to it. We do this because for some softening functions ℱ the coefficients ϕij(α) do not obey certain symmetry properties when a/ap is replaced with ap/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 Ad and Bd, we remind the reader that Ad represents the precession rate of the free ec-centricity vector of a test particle in the disc potential, while Bd characterizes the torque exerted on the particle orbit by the non-axisymmetric component of the disc gravity. Corresponding forced eccentricity vector is ep, f = −Bd/Ad. In par-ticular, test-particles initiated on circular orbits experience eccentricity oscillations of maximum amplitude epm=2|ep,f|. As Ad(ap) and Bd(ap) uniquely determine Rd 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. 2.1 Summary of existing softening models Here we provide a brief description of the four different softening prescriptions that have been previously proposed in the literature. Corresponding expressions for their softening parameters ∈2(α) and coefficients ϕij(α) are provided in Table 1. 2.1.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 BSm,Tr=2π∫0πcosmθ1+α2−2α cosθ+βc2−Sdθ.(12) Here βc2 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 Rd as a → ap, 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 ϕij (α) were expressed as derivatives of B1/2m,Tr 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. 2.1.2 Formalism of Touma (2002) – T02 Touma (2002) derived the orbit-averaged potential of a disc by assuming individual particles comprising the disc to interact via Plummer potential with a fixed length scale bc (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 ϕij(α) in the form of linear combinations of softened Laplace coefficients Bsm,T, similar to those of Tremaine (1998): BSm,T=2π∫0πcosmθ1+α2−2α cosθ+β2−Sdθ.(13) However, in Touma (2002) the softening parameter ∈2(α) = β2 is no longer a constant but depends on the distance such that β = bc/max(a1, a2). Within this formalism, one can think of a disc as comprised of nested annuli with a con-stant thickness bc. 2.1.3 Formalism of Hahn (2003) – H03 Hahn (2003) computed the orbit-averaged interaction between two eccentric wires by accounting for their vertical thickness. The vertical extent h of a ring effectively softens its gravitational potential over a dimensionless scaleH ∼ 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 BSm,H=2π∫0πcosmθ1+α2−2α cosθ+H21+α2−Sdθ(14) with constant H ≪ 1. In other words, the softening parameter is given by ∈2(α) = H2(1 + α2) in that work. The explicit expressions for ϕij(α) in terms of Bsm,H are given by equations (17) of Hahn (2003). 2.1.4 Formalism of Teyssandier & Ogilvie (2016) – TO16 Teyssandier & Ogilvie (2016) modified the unsoftened expressions (8), (9) for ϕijLLα by simply replacing the usual Laplace coefficients bsm with softened versions defined such that BSm,TO=2π∫0πcosmθ1+α2−2α cosθ+S2α−Sdθ.(15) 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 ϕ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, bc, 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 ς. 2.2 The unsoftened Heppenheimer method A different approach to computing the disturbing function of a razor-thin disc has been developed by Heppenheimer (1980) without resorting to any form of softened gravity (see also Ward 1981). The essence of this method is in computing the potential by direct integration over the disc surface before expanding the integral limits (which involve instanta-neous particle position r) in terms of small eccentricity of a test particle5. 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 b1/2mα – and not b3/2m – which diverge only weakly (logarithmically) as α→1:b1/2mα∝log1−α. As a result, upon in-tegrating these expressions over the radial extent of the disc,one obtains a convergent and finite result for Rd . 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. 3 COMPARISON: POWER-LAW DISCS Our goal is to examine the performance of different softening prescriptions outlined in §2.1 in comparison with the results obtained using the un-softened Heppenheimer method (§2.2). We start this exercise using a model of apse-aligned (i.e. dω¯d/da=0), truncated power-law (hereafter PL) disc as a simple example. We characterize surface density and eccentricity of such a disc by ∑da=∑0a0ap, eda=e0a0aq(16) 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 Rd due to PL discs can be simplified to (Silsbee & Rafikov 2015) Rd=Kψ1ep2+ψ2epedapcosω¯p−ω¯d,(17) where K=πG∑0a0pap1−p and the dimensionless coefficients ψ1 and ψ2 are given by ψ1=2∫α11α1−pϕ22αdα+2∫α21αp−2ϕ11αdα,(18) ψ2=2∫α11α1−p−qϕ12αdα+2∫α21αp+q−2ϕ12αdα,(19) with α1 = ain/ap and α2 = ap/aout. The coefficients ψ1 and ψ2 are functions of the power-law indices (p and q), any softening parameter involved (through ϕij), as well as the test-particle semi-major axis ap (through α1,2). They are related to Ad and Bd via Adap=2Knpap2ψ1, Bdap=Knpap2edapψ2.(20) 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 Σd(a) and ed(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 ap (i.e. functions of p, q, and ς only), making them useful as simple metrics for judg-ing the validity of different models of softening. 3.1 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”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 theun-softened Heppenheimer approach7. 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 generally8 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 ≪ ap ≪ aout (see Figure 4). 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 ς = β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 todemonstrate that these quantitative differences do not vanish by further decreasing βc. For instance, when p = 1, the coefficient ψ1 can be evaluated analytically as ψ1Tr98=12βc2+1+E2/βc2+4πβc2+4=−12+12π+Oβc2(21) 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 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 B3/2m,TO in the limit of α → 1 (Eq. C7) to demonstrate that both ψ1 and ψ2 (Eqs. 18 - 19) behave as ψ1TO16≈12S+OS, ψ2TO16≈1S+OS(22) as S → 0 for all values of p and q. The behavior shown in Fig. 1C, D agrees with these asymptotic expressions. 3.2 Details of convergence of different softening prescriptions Different softening prescriptions explored in this work are designed to modify the behavior of the integrand in equations (5)-(6) primarily in the vicinity of the test particle orbit, i.e. as a → ap 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 Σd(a) is set to zero, are at ad,i = (1 − x)ap ≤ ap andad,o = (1− x)−1ap ≥ ap, 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 ψ˜1x and ψ˜2x, 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 ad,i(x) < a < ad,o(x). In Figure 2 we display the behavior of ψ˜1x (Panel A) and ψ˜2x (Panel B) as a function of x=1−ad,i/ad,0 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 ψ˜1x (panel A) and ψ˜2x (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 ςap, i.e. ς ≲ x ≤ 1, the amplitudes of both ψ˜1x and ψ˜2x 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 ψ˜1x∼|ψ˜2x|∼x−1 ∼ 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, Bsmα,ς≈bsmα, and all ϕij effectively reduce to their classical counterparts ϕijLL given by Eqs. (8) - (9), which can be easily verified using the expressions listed in Table 1. The scaling of ψ˜1x and ψ˜2x with x is simply a result of asymptotic behavior of b3/2mα→1−α−2 → (1 − α)−2 as α → 1, upon radial integration in Eqs. (18) – (19). Second, upon reaching their extrema at x ∼ ς, amplitudes of ψ˜1x and ψ˜2x 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 Bsmα,ς away from the divergent behavior of bsmα. The modification is such that the softened interaction with the disc annuli ≲ ςap 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, ψ˜1 and ψ˜2 cross zero and change sign at some x = Cς2, where C ∼ 1 is a constant9. At the same time, ψ˜1TO16 and ψ˜2TO16 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 ψ˜1TO16 and ψ˜2TO16 converge keeps increasing as ς → 0. Moreover, both coefficients also never change sign, always predicting prograde precession (ψ˜1TO16>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 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 ≲ ep), 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 ψ˜1x 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, ψ˜1x∼x−1 ∼ x−1 scaling for ς ≲ x ≪ 1, decay of ψ˜1x 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 zero10 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 ψ˜1x 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. 3.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 prescriptions11 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, epm/eda=ψ2/ψ1 = |ψ2/ψ1|, and the expectations based on SR15. The latter point is illustrated in Figure 4C. 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 epm of eccentricity oscillations (normalized by disc eccentricity ed) 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 epm 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 epm curves (panel C). 4 COMPARISON: NON-POWER-LAW DISCS We now turn our attention to the performance of the different softening prescriptions for more general discs. Namely, we focus on two apse-aligned, non-PL disc models previously studied by Davydenkova & Rafikov (2018) based on the unsoftened Heppenheimer method. The dynamics in such non-PL discs, according to DR18, differ from the PL discs in a very important way: the free-precession of test-particles can naturally change from retrograde to prograde (and vice versa) within such discs. Furthermore, an important feature of the models considered below is that Σd smoothly goes to zero at finite radii in a manner that does not give rise to the edge effects, see DR18 and §5. 4.1 Quartic Disc Model We start by looking at the secular dynamics in the potential of a Quartic disc characterized by the surface density ∑da=∑0aout−a2ain−a2aout−ain4,(23) and linear eccentricity profile in the form eda=e˜01+aout−aaout−ain(24) for ain ≤ a ≤ aout (with ain = 0.1 AU, aout = 5 AU), where ∑˜0=1153 g cm−2 = 1153 g cm−2 and e˜0=0.01 = 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 Ad, Bd and epm 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 Ad and Bd computed using this model from Davydenkova & Rafikov (2018), in particular, modify the locations at which Ad and Bd become zero. This explains the slight shift in the semi-major axes at which epm=2Bd/Ad = 2Bd/Ad 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 ap are the profiles of (A) the amplitude epm of the disc-induced eccentricity oscillations, (B) the rate of disc-driven free precession Ad, and (C) the coefficient Bd 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. 4.2 Gaussian Rings Next we investigate secular dynamics in the potential of another disc model from DR18 — a Gaussian ring with the surface density profile ∑da=∑˜0exp4−a/ac+ac/a2wc(25) centered around ac = 1.5 AU with width wc = 0.18 and surface density ∑˜0= = 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 Ad, Bd and epm 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 Σd(a) and ed(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. 5 EFFECTS OF PROXIMITY TO THE DISC EDGE So far the disc models that we explored were either infinitely extended (§3) or had surface density smoothly petering out to zero at finite radii (§4). This allowed us to not worry about the effects of sharp disc edges — discontinuous drops of the surface density — on secular dynamics, which are known to be important (Silsbee & Rafikov 2015; Davydenkova & Rafikov 2018). We now relax this assumption and examine the performance of different softening models in the vicinity of a sharp edge of the disc, where surface density drops discontinuously from a finite value to zero at a finite semi-major axis a = aedge. To that effect we analyze the behavior of secular coefficient Ad 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 Bd are very similar) near the disc edge. Figure 7 shows the run of Ad near the inner edge ain of the disc for particles both inside (ap < ain) and outside (ap > 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 Ad near the inner edge ain = 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 Ad near the disc edge is reproduced by the softening prescription of Hahn (2003) in the Ad 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 Ad, as well as Bd, 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 ≪ ap ≪ aout, see Fig. 1. For ς = 10−2 the softened calculation predicts Ad 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 softening12. 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 Ad 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 Σ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 Ad curve connecting the two peaks turns into a vertical line at ain. Similar divergent behavior of Ad (and Bd) 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. ap = a) is naturally excluded thus avoiding the classical singularity. Outside the disc the unsoftened calculations based on the Heppenheimer method (e.g. SR15, DR18) reduce to the Laplace-Lagrange theory exactly. 6 DISCUSSION Results of previous sections reveal a diversity of outcomes when different softening models are applied. Two models — those of Hahn (2003) and Touma (2002) — successfully reproduce the un-softened calculations based on the Heppenheimer method in the limit of zero softening. In the same limit, the formalism of Tremaine (1998) yields convergent results which are, however, different from the un-softened calculations, typically by tens of per cent. Finally, the softening method of Teyssandier & Ogilvie (2016) does not lead to convergent results in the limit of vanishing softening parameter. Interestingly, the two successful models (Hahn 2003; Touma 2002) have been derived using rather different underlying assumptions (see §2.1.2 & 2.1.3), producing different mathematical expressions for ϕij (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 ϕij (thus fully determining the softened secular model via Eqs. (4)-(6)) given an arbitrary softened two-point interaction potential in the form (10). This procedure involves orbit-averaging the softened potential along the particle trajectories; its details are presented in Appendix A. There is also an alternative approach, sketched in Appendix A4, which assumes the disc to be a continuous entity from the start. Both of them arrive at the same expressions for Rd. 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 ℱ (r1, r2) = bc2 and ℱ (r1, r2) = H2r12+r22, 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 Bsm,Tr defined by Eq. (12) suggests interaction potential (10) with ℱ (r1, r2) = βc2maxr12,r22 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 ϕij: ϕ11=ϕ22=α8B3/21,Tr−3αβc2B5/20,Tr−δα−1βc2B3/20,Tr,(26) ϕ12=α4B3/22,Tr−3αβc2B5/21,Tr−δα−1βc2B3/21,Tr.(27) 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 ϕij 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 Ad and Bd (Eqs. 5, 6) in apse-aligned discs in the form δAdap=−πG2npapβc2∑dapB3/20,Tr|α=1,(28) δBdap=+πG2npapβc2∑dapedapB3/21,Tr|α=1,(29) Accounting for these corrections, we confirmed that the correct (un-softened) behavior of the coefficients of Rd can be reproduced for the non-PL discs – Quartic and Gaussian models, see §4. Note that δ Ad(ap) and δBd(ap) are proportional to the local disc surface density £d(ap) and B3/2m,Tr(α = 1) ∼ βc−2, 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): δψ1=−14βc2B3/20,Tr|α=1=−12π+Oβc2,(30) δψ2=12βc2B3/21,Tr|α=1=1π+Oβc2,(31) see Eqs. (20). These corrections exactly match the offsets seen in Fig. 4 between the calculations of Tremaine (1998) and the un-softened calculations, thus explaining the origin of these uniform shifts. We also confirmed this explanation in Fig. 8, where we show the convergence of modified Tremaine (1998) coefficients to the correct un-softened values as softening is varied for 2 values of p and q. To summarize, Eqs. (26)-(27) should replace the expressions given by Eq. (26) of Tremaine (1998) in applications to continuous discs. However, when considering the interaction of two individual annuli with different semi-major axes (like in the classical Laplace-Largange theory), one has α ≠ 1 and terms in Eqs. (26)-(27) containing δ-function naturally vanish, reducing ψ1 and ψ2 back to the expressions quoted in Tremaine (1998). 6.2 On the softening prescription of Teyssandier & Ogilvie (2016) We now turn our attention to the model of Teyssandier & Ogilvie (2016) trying to understand its distinct (divergent) behavior. From the expression for Bsm,TO 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 b3/2m in the classical expressions (8, 9) for ϕijLL with B3/2m,TO, see §2.1.4; this simple swap of Laplace coefficients has not been justified rigorously. On the other hand, in Appendix B we show that softening parameter in the form ∈2(α) = ς2α corresponds to softening function ℱ(r1, r2) = ς2r1r2 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 ϕij with ς = S (Appendix B): ϕ11=ϕ22=α8B3/21,TO+12S2B3/20,TO−34S2B3/20,TO−34S22+2α2+S2αB5/20,TO,(32) ϕ12=−α4B3/22,TO+12S2B3/21,TO−34S22+2α2+S2αB5/21,TO.(33) 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 ϕ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(α) = 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. Figure 8 Similar to Figure 1, but now using the expressions for ϕij given by Eqs. (26-27) and Eqs. (32-33) obtained by propagating fr1,r2=ς2maxr12r22 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. 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 b3/2m in Eq. (1) with their softened analogues is, evidently, not sufficient for obtaining a well-behaved softened version of Laplace-Lagrange theory for co-planar discs. 6.3 Implications for numerical applications In numerical studies of secular dynamics, self-gravitating discs are often treated as a collection of N eccentric annuli (rings), with prescribed spacing (justified by the constancy of the semi-major axis), interacting gravitationally with each other (e.g. Touma et al. 2009; Batygin 2012). This representation approximates a continuous particulate or fluid disc in the limit of N → ∞. 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 ψ˜1x and ψ˜2x computed without accounting14 for the contributions from ad,i < ap < ad,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 ψ˜1x and ψ˜2x, 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/ap ≲ 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 (ςap), 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 Ad (or ψ1 for PL discs) by using the following global metric15 Μf=∫ainaoutftheora−fnuma2da∫ainaoutftheor2ada.(34) Here fnum(ai) is the value of the metric basis (e.g. precession rate Ad) evaluated at the position ai of ith ring by summing up the contributions of all other rings in the disc, while ftheor(ai) 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 = Ad, 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 legend16. Figure 9 shows that as ς → 0, the number of rings scales as N ∼ Cς−χ with17 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 Ad (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/ap ≲ 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. 6.4 Further generalizations and extensions All calculations in this work are based on the expansion of the secular disturbing function Rd 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 Ad(ap) = 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 Rd to 4th order in eccentricities based on the un-softened method of Heppenheimer (1980). The general framework for calculating Rd with arbitrary softening prescriptions presented in Appendix A can also be extended to higher order in eccentricities in similar way18, see e.g. Touma & Sridhar (2012). We expect that conclusions similar to those drawn from our analysis in §3-5 will also apply to the higher-order expansions. Additionally, although we only analyzed coplanar configurations in this work, the general framework presented in Appendix A may be extended to account for non-coplanar configurations and study the inclination dynamics. 7 SUMMARY In this work we investigated the applicability of softened gravity for computing the orbit-averaged potential of razor-thin eccentric discs. We compared disc-driven secular dynamics of coplanar test-particles computed using softening prescriptions available in the literature with the calculations based on the unsoftened method of Heppenheimer (1980). Our findings are summarized below. We confirmed that the softening methods of both Touma (2002) and Hahn (2003) correctly reproduce eccentricity dynamics of razor-thin discs in the limit of vanishing softening parameter ς 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. ACKNOWLEDGEMENTS We express our gratitude to Scott Tremaine and Jihad Touma for a number of stimulating discussions, which have led to substantial improvements of the manuscript. We are also grateful to Gordon Ogilvie, Jean Teyssandier, Yoram Lithwick, and Cristobal Petrovich for useful discussions, and an anonymous referee for constructive comments. A.A.S. acknowledges a scholarship by the Gates Cambridge Trust (OPP1144), while R.R.R. was supported by NASA via grant 15-XRP15-2-0139. Open Access for this article was funded by the Bill & Melinda Gates Foundation. 1 The Laplace coefficients entering in δR can be easily evaluated, without relying on integration over θ in Eq. (2), by expressing them through elliptic integrals, see Appendix C3. 2 We refer the reader to Heppenheimer (1980); Silsbee & Rafikov (2015); Davydenkova & Rafikov (2018) for the expressions of Ad and Bd computed using the un-softened Heppenheimer method for different disc models. 3 Note that the inter-particle force resulting from such potential does not, in general, obey Newton’s third law (as long as ℱ(r1, r2) ≠ const). 4 Here we clarify that the definitions of ϕ11(α) and ϕ22(α), even when different (see Table 1 and Appendix A), are swapped upon interchanging a1 with a2 but keeping, by construction, α = a < 1 – see Eqs. (A22), (A23) for details. 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). 6 The softening length bc present in the formulation of Touma (2002) is scaled by the test-particle semi-major axis ap in all the figures where we present results for infinite PL discs. We do this to properly collate the results computed by different softening formalisms in one figure. 7 Equations (A37) and (A38) in Silsbee & Rafikov (2015) provide analytic expressions for ψ1 and ψ2, respectively, for infinite PL discs. 8 For particles with orbits near sharp disc edges, we find that smaller values of ς is required to recover the expected dynamics, see §5. 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 CTr98 = (π − 1)/2 and CH03 = π; 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. 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 Fig. 4. 11 We do not present results obtained by the method of Teyssandier & Ogilvie (2016). 12 On the other hand, this condition is relaxed when the edge is not exactly sharp but rather has a finite width ∆r over which the disc surface density smoothly peters out to zero; in this case ς only needs to be ≲ ∆r/r. 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 b1/2m with B1/2m in Eq. (7.1) of Murray & Dermott (1999) before applying the derivatives with respect to α. We note that this procedure could apply for all ℱ(r1, r2) with continuous first derivatives satisfying D1 + D2 = −1; see Appendix A. 14 Note that, technically, in the un-softened case this mathematical procedure is not equivalent to introducing an actual physical gap in the disc, as the latter would result in additional boundary terms. 15 For PL discs, we neglect rings within 10% of disc edges when computing ℳ(ψ1). 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 = b2/max(a12,a22), 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. 18 Another way to calculate the softened disturbing function for arbitrarily high eccentricities is to numerically compute the ring-ring interaction potential, as was done by Touma et al. (2009). 19 Note that we do not deal with the indirect part of the potential – which is left unsoftened – as it contains only periodic terms and does not affect the secular dynamics (Murray & Dermott 1999). 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 ℱ 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 as19 ψ=r1−r22+Fr1,r2−1/2 =r12+r22−2r1r2cosf1−f2+ω¯1−ω¯2+Fr1,r2−1/2,(A1) where ℱ(r1, r2) is an arbitrary softening function introduced to cushion the singularity which arises otherwise at null inter-particle separations. In the above expression, fi is the true anomaly of the ith ring, ∈i is its longitude of periapse and ri is its instantaneous position, i = 1, 2. Our goal is to obtain the orbit-averaged expansion of Ψ to second order in eccentricities ei valid for arbitrary ℱ(r1, r2). 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 ψ=explogr1a1D1+logr2a2D2+f1−M1D3+f2−M2D4ψ0 ≡ Τψ0(A2) with ψ0=a12+a12−2a1a2cosθ+Fa1,a2−1/2,(A3) where θ = M1 − M2 + ∈1 − ∈2, Mi represents the mean anomaly of the ith ring characterized with semi-major axis ai, and the linear operators Dk are given by (Plummer 1918) D1=a1∂∂a1≡a1∂1 , D2=a2∂∂a2≡a2∂2 , and D3=−D4=∂∂θ . (A4) 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 ℱ − M, a−1rD=1−ecosM⋅D+12e21−cos2M⋅D+14e21+cos2M⋅DD−1+Oe3,(A5) expf−MD = 1+2esinM⋅D+54e2sin2M⋅D+e21−cos2M⋅D2+Oe3(A6) 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(M1 − M2), 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 ψ=Τψ0≡Αψ0e12+Βψ0e22+ℂψ0e1e2cosω¯1−ω¯2,(A7) where the operators A, B and C acting on Ψ0 are defined as Α ≡ D32+14D1D1+1, Β≡D42+14D2D2+1,(A8) ℂ ≡ cosθ2D3D4+12D1D2−sinθD2D3−D1D4.(A9) We have used the fact that cos(M1 − M2) = cos θ cos(∈1 − ∈2) and sin(M1 − M2) = sin θ cos(∈1 − ∈2) in the secular regime (Plummer 1918). A2 Computation of the action of relevant operators 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: D32ψ0 = D42 = 3a12a22sin2θψ05−a1a2cosθψ03,(A10) D1D2ψ0=a1a2cosθ−12∂1∂2Fψ03+3a22+a1a2cosθ+a22∂2Fa12−a1a2cosθ+a12∂1Fψ05,(A11) D2D3ψ0 = −a1a2sinθψ03+3a1a2sinθa22−a1a2cosθ+a22∂2Fψ05,(A12) D1D4ψ0 = a1a2sinθψ03−3a1a2sinθa12−a1a2cosθ+a12∂1Fψ05,(A13) D1ψ03 = −3a12−a1a2cosθ+a12∂1Fψ05,(A14) D2ψ03 = −3a22−a1a2cosθ+a22∂2Fψ05,(A15) where for conciseness we have written ℱ instead of ℱ(a1, a2). Here, it is worthwhile to mention that, as far as the expansion technique is concerned, the terms ∂i ℱ (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: D1D1+1ψ0 = −D1D2ψ0+12D12F−a1∂1F−a2∂2Fψ03,(A16) D2D2+1ψ0 = −D1D2ψ0+12D22F−a1∂1F−a2∂2Fψ03,(A17) which can be obtained by making use of the identity D1+D2+1Ψ0=122f−a1∂1f−a2∂2fΨ03. Here, we note that for all softening functions ℱ for which 2ℱ − a1 ∂1 ℱ − a2 ∂2 ℱ = 0, one finds D1 + D2 = −1. Consequently, in such cases, the operators D1(D1 + 1) and D2(D2 + 1) become identical rendering AΨ0 = BΨ0 (since D32=D42, see Eqs. (A8) and (A10)). As a result, the resultant orbit-averaged disturbing function (A7) is symmetric in e1 and e2, similar to the case of classical Laplace-Lagrange theory. This is not true in general, for instance, when ℱ(r1, r2) = 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 ψ=Τψ0=12π∫02πΤψ0dθ,(A18) which essentially reduces to computing the individual terms ΑΨ0, ΒΨ0 and ℂΨ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 cosmθΨ02s. By making use of α = a< /a>, where a< = min(a1, a2) and a> = max(a1, a2), this combination can be reduced to cosmθψ02s=a>−2scosmθ1+α2−2αcosθ+a>−2Fa1,a2−s.(A19) For that reason, calculation of the orbit-averaged Ψ (by integrating over dθ) yields integrals of the form Bsmα≡2π∫0πcosmθ1+α2−2αcosθ+∈2α−sdθ,(A20) which is the generalization of the classical Laplace coefficients bsm (recovered when ℱ (a1, a2) = 0, see Eq. 2) with the dimensionless softening parameter ∈2α≡a>−2Fa1,a2,(A21) see Eq. (7). Employing this notation, we present the simplified expressions of ΑΨ0, ΒΨ0 and a12,a22 obtained as a result of orbit-averaging: a>Αψ0α≡ϕ11α = α2−54Β3/21+38αΒ5/20+341+α2Β5/21−158αΒ5/22+38T2Β5/21−316T5B5/20+18T3+α−1T4B3/20−38T1a1a2B5/20−B5/21+12T7B5/20,(A22) a>Bψ0α≡ϕ22α = α2−54Β3/21+38αΒ5/20+341+α2Β5/21−158αΒ5/22+38T2Β5/21−316T5B5/20+18T3+α−1T6B3/20−38T1a1a2B5/20−B5/21+12T8B5/20,(A23) a>ℂψ0α≡ϕ12α = α2−94Β3/20+14Β3/22+38αΒ5/23−218αΒ5/21+341+α2Β5/22−941+α2B5/20−14T3B3/21−98T2B5/20+38T5B5/21+38T2B5/22.(A24) In equations (A22)-(A24), we have defined the dimensionless functions Ti(α) such that T1 = a>−22F−a1∂1F−a2∂2F, T2 =α∂1Fa2+∂2Fa1, T3 = ∂1∂2F,(A25) T4 = a1a>2∂12F−a1∂1F−a2∂2F, T5 = α2∂1Fa1+2∂2Fa2+∂1Fa1∂2Fa2,(A26) T6 = a2a>2∂22F−a1∂1F−a2∂2F, T7 = α2−1∂1F, T8=a1−1∂2F,(A27) where, as before, ℱ ≡ ℱ (a1, a2), α = a< /a> and ∂i ≡ ∂/∂ai . Note that the expressions for φ11 and φ22 swap definitions upon replacing a1 by a2, whilst keeping α < 1 by construction. This can be understood by first noting that functions Ti with i = 1, 2, 3, and 5 are invariant under a1 ⇌ a2 while, at the same time, T4 and T7 (appearing in the second line of Eq. (A22)) translate to T6 and T8 (appearing in the second line of Eq. (A23)); and vice versa. These identities, when combined, yield the desired expression of D1+D2+1Ψ0=122f−a1∂1f−a2∂2fΨ03; 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 ℱ (r1, r2). Note that in the absence of softening (i.e. ℱ (r1, r2) = 0) Ti = 0 for all i and the classical expressions for D32=D42, ΑΨ0 and ΒΨ0 — 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 ℂΨ0, where Φ is the interaction potential given by equation (10), angle brackets indicate averaging over the orbit of the test particle given by rp and integration is carried out over the full surface of the disc S with rd denoting the location of a disc element. To obtain the expression for Rd accurate to second order in eccentricities one would need to expand Φ(rd, rp) to second order in particle and disc eccentricities by e.g. writing rp = ap (1 − ep cos Ep), where Ep is the eccentric anomaly of the particle orbit. This expansion should explicitly account for the dependence of ℱ on rd and rp. Averaging the resulting expressions over Ep, one would arrive at the proper expression for Rd in the form (4). In particular, after a lengthy but straightforward calculation this method gives the following expression for the disc-driven precession rate: Ad=πG2npap2∫a∑adaa> 143apF′F′+4ap−22F′+apF″ap2+a2+F−12apFapB5/20αa>4+αB3/21α−F′−apF″apαB5/21αa>2,(A28) where prime denotes differentiation with respect to ap (e.g. ℱ′ = ∂ℱ /∂ap), a> = max(ap, a), α = min(ap, a)/max(ap, a) and integration is done over the semi-major axis a of the disk elements. Calculation of the non-axisymmetric part of Rd resulting from non-zero disk eccentricity (i.e. Bd) is somewhat more tedious but can nevertheless be done similar to Davydenkova & Rafikov (2018). APPENDIX B SPECIFIC CASES OF ℱ(R1, R2) The general framework developed in Appendix A allows us to recover the expressions of φij arrived at by Touma (2002) and Hahn (2003) upon specifying certain functional forms of ℱ (r1, r2). Indeed, Touma (2002) performed the same calculations as presented in Appendix A for the case of Plummer potential – cosmθΨ02s – to second order in eccentricities, and later to fourth order in eccentricities (Touma & Sridhar 2012). Furthermore, we find that the results obtained by Hahn (2003) can be recovered from our general framework by setting bsm. For reference, the functional forms of Ti for these forms of ℱ(r1, r2), along with their softening parameters ∈2(α), are summarized in Table B1, which can be used to show that Eqs. (A22)-(A24) reduce to those in Table 1 after some algebra with the aid of the recursive relationships for ΑΨ0 presented in Appendix C. Table B1 The functional forms of the coefficients Ti (α) 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 ∈2a=a>−2Fa1,a2 (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 H2r12+r22 H2(1 + α2) 0 2H2(1 + α2) 0 0 4αH2(2 + H2) 0 2H2a1a2 2H2a2a1 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αΘa1−a2 2βc2αΘa2−a1 T016 S2r1r2 S2α 0 2αS2 S2 0 S2(S2α + 2α2 + 2) 0 S2 S2 As to the formalism of Teyssandier & Ogilvie (2016), we find, using their softening prescription of ℱ (r1, r2) = S2r1r2, that our general framework yields φij expressions different from those reported by Teyssandier & Ogilvie (2016). Indeed, we first note that in this case, T1 = T4 = T6 = 0 (Table B1) rendering the expressions of φ11 and φ22 identical such that ϕ11=ϕ22=α8−5B3/21,TO+32αB5/20,TO+31+α2+S2αB5/21,TO−152αB5/22,TO−34S2S2α+2α2+2B5/20,TO+12S2B3/20,TO(B1) Using the recursive relationships listed in Appendix C1, the above expression can be simplified further. Indeed, Eq. (C2) with m = 1 and s = 5/2 and Eq. (C1) with m = 1 and s = 3/2 read 31+α2+S2αB5/21,TO=−3α2B5/22,TO+152αB5/20,TO,(B2) −6B3/21,TO=9αB5/22,|TO−B5/20,TO,(B3) respectively. Inserting the above two identities in Eq. (B1) one arrives at Eq. (32). Similarly, the expression of φ12 (Eq. A24) can be simplified with the aid of Eq. (C3) (with m = 0, s = 3/2), Eq. (C2) (with m = 2, s = 5/2) and Eq. (C1) (with m = 2, s = 3/2) resulting in Eq. (33) after some algebra. As discussed in §6.2, the terms in Eqs. (32)-(33) explicitly proportional to S2 are absent in the original formulation of Teyssandier & Ogilvie (2016) (see Table 1). Similarly, for the formalism of Tremaine (1998), propagating their functional form of ΒΨ0 through our general framework, we arrive at the expressions for φij(α) differing from those reported in Tremaine (1998) in a very special way: we find φij to contain additionl terms proportional to T3(α) ∼ δ(α − 1), where δ(x) is the Dirac delta-function. Such terms are absent in the original formulation of Tremaine (1998) (see Tables 1, B1). Emergence of these terms can be easily demonstrated by first noting that in this case φ11 = φ22 (as T1 = T4 = T6 = 0), employing the recursive relationships for Laplace coefficients (in a similar order as done above for TO16) to simplify the general expressions of φ11(= φ22) and φ12, and finally arriving at Eqs. (26), (27). The ramifications of this finding is discussed in Section 6.1. APPENDIX C GENERALIZED LAPLACE COEFFICIENTS As demonstrated in Appendix A, softening the Newtonian point-mass potential by an arbitrary function ℱ (r1, r2) modifies the definition of the Laplace coefficients as shown by Eqs. (7), (A20) by the introduction of a softening parameter fr1,r2=βc2maxr12,r22 (Eq. A21), 0 ≤ α = a< /a> ≤ 1. Here, we present some useful recursive relationships amongst different generalized Laplace coefficients ∈2α=a>−2fa1,a2, along with their asymptotic behavior in the limits of α → 0, 1 as well as their relationship to complete elliptic integrals. C1 Recursive Relations Generalizing the results for the usual (unsoftened) Laplace coefficients Bsmα (e.g. Plummer 1918, p. 159), the following relationships can be easily obtained for the generalized Laplace coefficients defined by Eq. (7), (A20): mBsm=sαBs+1m−1−sαBs+1m+1,(C1) m1+α2+∈2Bsm=αm+1−sBsm+1+αm+s−1Bsm−1,(C2) m+sBsm=s1+α2+∈2Bs+1m−2sαBs+1m+1,(C3) The difference with the classical recursive relations for bsm amounts to substituting the combination 1 + α2 appearing in the case of ordinary Laplace coefficients with 1 + α2 + ∈2(α). Another useful expression relating the generalized Laplace coefficients of arguments α and α−1 is Bsmα−1=α2sBsmα.(C4) Note that the above relationship is valid only as long as the softening parameter satisfies α2∈2(1/α) = ∈2(α). For instance, this condition is violated when the softening parameter ∈ has no dependence on α (e.g. that of Tremaine (1998), see Table 1). C2 Asymptotic Behavior Here we derive approximate expressions for bsm in the asymptotic limits; for α → 0 and α → 1. Case 1: In the limit of α ≈ 0, one can factor out the term 1 + α2 + ∈2(α) from the integrand of Bsm to expand the denominator around γ−1 ≈ 0, where γ = (2α)−1[1 + α2 + ∈2(α)]. This allows us to approximate Bsm as Bsmα≈2π(2αγ)s∫0πcosmθ×1+sγcosθ+ss+12γ2cos2θ+ss+1s+26γ3cos3θdθ.(C5) Using the orthogonality of the cosine functions, it is straightforward to show that Bsm≈αmFm2αγs+m' as α→0, where Fm=2 if m=02s if m=1ss+1 if m=213ss+1s+2 if m=3(C6) Case 2: In the opposite limit of x = 1−α ≈ 0, the dominant contribution to Bsm comes from θ ≪ 1 (Goldreich & Tremaine 1980). Thus one can set cos(mθ) → 1 in the numerator, approximate cosθ ≈ 1 − θ2/2 in the denominator and extend the integration limit to infinity. Furthermore, setting α = 1 (i.e. x = 0) everywhere except when it appears in the combination 1− α, the generalized Laplace coefficient can be approximated as Bsm≈2π∫0∞dθx2+θ2+∈α=12s=2πx2+∈α=12−1 if s=3/22/3x2+∈α=12−2 if s=5/2(C7) where Bsm is the softening parameter evaluated at α = 1. C3 Relationship to elliptic integrals Here we express the generalized Laplace coefficients ∈α=12 in terms of complete elliptic integrals. These expressions can be used for rapid numerical evaluation of the generalized Laplace coefficients without relying on numerical integration of Eq. (A20) (or Eq. (7)). Let us write, as before, 2αγ = 1 + α2 + ǫ2(α) and define χ=2/γ+1 such that, for any general softening parameter ∈2(α), we have 0 ≤ χ ≤ 1 and γ ≥ 1. Now let us express Bsm in terms of γ to write Bsm=21−sπαs∫0πcosmθγ−cosθsdθ(C8) Introducing complete elliptic integrals Κχ=∫0π/21−χ2sin2ϕ−1/2dϕ and Εχ=∫0π/21−χ2sin2ϕ−1/2dϕ we find that B3/20=2Exπαγ−12αγ+1 B3/21=2−(γ−1)K(x)+γExπαγ−12αγ+1(C9) B3/22=2−4γγ−1Kx+4γ2−3Exπαγ−12αγ+1 B3/23=23−γ−1(32γ2−5)K(x)+γ32y2−29Exπαγ−12αγ+1(C10) B5/20=4−γ−1Kx+4γEx3π2α5/2γ+13/2γ−12 B5/21=4−γγ−1Kx+γ2+3Ex3π2α5/2γ+13/2γ−12(C11) B5/22=4γ−1(4γ2−5)Kx−4γγ2−2Ex3π2α5/2γ+13/2γ−12 B5/23=4γγ−132γ2−33Kx−32γ4−57γ2+21Ex3π2α5/2γ+13/2γ−12(C12) These expressions permit efficient numerical evaluation of arbitrarily softened Laplace coefficients as functions of α, since effective algorithms for computing K and E exist (Press et al. 2002). APPENDIX D CONVERGENCE CRITERION FOR THE PRE-FACTORS OF POWER-LAW DISCS Astrophysical discs often extend over a few orders of magnitude in radius so that aout/ain ≫ 1. In such situations, far from the disc edges one can take the limit of both α1 = ain/ap and α2 = ap /aout 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 Bsm, 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≈−14S2α+381+32|S4α2  and  ϕ12≈32S2α2−15161+5S4α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 bc → 0. However, when bc 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 TEX/LATEX file prepared by the author. ==== Refs REFERENCES Batygin K. , 2012 , Nature , 491 , 418 Binney J. , Tremaine S. , 2008 , Galactic Dynamics: Second Edition . Princeton University Press Chiang E. I. , Goldreich P. , 2000 , ApJ , 540 , 1084 Davydenkova I. , Rafikov R. R. , 2018 , ApJ , 864 , 74 Fontana A. , Marzari F. , 2016 , A&A , 589 , A133 Goldreich P. , Tremaine S. , 1979 , AJ , 84 , 1638 Goldreich P. , Tremaine S. , 1980 , ApJ , 241 , 425 Graps A. L. , Showalter M. R. , Lissauer J. J. , Kary D. M. , 1995 , AJ , 109 , 2262 Gulati M. , Saini T. D. , Sridhar S. , 2012 , MNRAS , 424 , 348 Hahn J. M. , 2003 , ApJ , 595 , 531 Heppenheimer T. A. , 1980 , Icarus , 41 , 76 Jacobs V. , Sellwood J. A. , 2001 , ApJ , 555 , L25 Kazandjian M. V. , Touma J. R. , 2013 , MNRAS , 430 , 2732 Kocsis B. , Tremaine S. , 2011 , MNRAS , 412 , 187 Latter H. N. , Ogilvie G. I. , Rein H. , 2017 , preprint , (arXiv:1701.04312) Lee W.-K. , Dempsey A. M. , Lithwick Y. , 2019 , ApJ , 872 , 184 Meschiari S. , 2014 , ApJ , 790 , 41 Murray C. D. , Dermott S. F. , 1999 , Solar system dynamics . Cambridge University Press Pan M. , Wu Y. , 2016 , ApJ , 821 , 18 Petrovich C. , Wu Y. , Ali-Dib M. , 2019 , AJ , 157 , 5 Plummer H. C. K. , 1918 , An introductory treatise on dynamical astronomy . Cambridge University Press Press W. H. , Teukolsky S. A. , Vetterling W. T. , Flannery B. P. , 2002 , Numerical recipes in C++ : the art of scientific computing . Cambridge University Press Sefilian A. A. , Touma J. R. , 2019 , AJ , 157 , 59 Silsbee K. , Rafikov R. R. , 2015 , ApJ , 798 , 71 Statler T. S. , 2001 , AJ , 122 , 2257 Teyssandier J. , Ogilvie G. I. , 2016 , MNRAS , 458 , 3221 Teyssandier J. , Terquem C. , Papaloizou J. C. B. , 2013 , MNRAS , 428 , 658 Tiscareno M. S. , 2013 , Planetary Rings . Springer , p. 309 , doi:10.1007/978-94-007-5606-97 Touma J. R. , 2002 , MNRAS , 333 , 583 Touma J. R. , Sridhar S. , 2012 , MNRAS , 423 , 2083 Touma J. R. , Tremaine S. , Kazandjian M. V. , 2009 , MNRAS , 394 , 1085 Tremaine S. , 1998 , AJ , 116 , 2015 Tremaine S. , 2001 , AJ , 121 , 1776 Ward W. R. , 1981 , Icarus , 47 , 234 Ward W. R. , 1989 , ApJ , 336 , 526