Dynamical evolution of massive perturbers in realistic multi-component galaxy models I: implementation and validation
Matteo Bonetti, Elisa Bortolas, Alessandro Lupi, Massimo Dotti
MMNRAS , 1–15 (2020) Preprint 30 October 2020 Compiled using MNRAS L A TEX style file v3.0
Dynamical evolution of massive perturbers in realisticmulti-component galaxy models I: implementation and validation
Matteo Bonetti, , , Elisa Bortolas, Alessandro Lupi and Massimo Dotti , , Dipartimento di Fisica “G. Occhialini”, Università degli Studi di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy INFN, Sezione di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy INAF, Osservatorio Astronomico di Brera, Via E. Bianchi 46, I-23807, Merate, Italy Center for Theoretical Astrophysics and Cosmology, Institute for Computational Science, University of Zürich, Winterthurerstrasse190, CH-8057 Zürich, Switzerland Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
Galaxies are self-gravitating structures composed by several components encompassing spher-ical, axial and triaxial symmetry. Although real systems feature heterogeneous componentswhose properties are intimately connected, semi-analytical approaches often exploit the lin-earity of the Poisson’s equation to represent the potential and mass distribution of a multi-component galaxy as the sum of the individual components. In this work, we expand thesemi-analytical framework developed in Bonetti et al. (2020) by including both a detailedimplementation of the gravitational potential of exponential disc (modelled with a sech andan exponential vertical profile) and an accurate prescription for the dynamical friction expe-rienced by massive perturbers in composite galaxy models featuring rotating disc structures.Such improvements allow us to evolve arbitrary orbits either within or outside the galactic discplane. We validate the results obtained by our numerical model against public semi-analyticalcodes as well as full N-body simulations, finding that our model is in excellent agreementto the codes it is compared with. The ability to reproduce the relevant physical processesresponsible for the evolution of massive perturber orbits and its computational efficiency makeour framework perfectly suited for large parameter-space exploration studies. Key words: galaxies: nuclei – stars: kinematics and dynamics – gravitation – galaxies:structure – methods: numerical
In the last decades, more and more interest has grown towards themodelling of orbital motions within the potential of a galaxy (e.g.Bovy 2015; Boubert et al. 2020; Granados et al. 2021). This hasto be attributed to the advent of missions like Gaia (Gaia Collabo-ration et al. 2016, 2018), which is tracking the orbits of stars andstellar clusters sailing though the Galactic potential, and will releasean unprecedented catalogue of Galactic stellar orbits. In addition,gravitational wave observatories such as LISA (Amaro-Seoane et al.2017) and PTA (Verbiest et al. 2016) will probe massive black hole(MBH) mergers along the cosmic time; in order to interpret therates of such mergers, one needs to model the inspiral phase of anintruder MBH within the potential of the host galaxy, down to thescale at which it will form a binary with another black hole that mayalready sit in the centre (Begelman et al. 1980).A detailed modelling of galactic orbits can only be achievedwith a sufficiently good description of the associated galactic poten-tial wells; in addition, if one wants to accurately follow the orbit ofa massive perturber (MP, i.e. an object whose mass is significantly larger than the characteristic stellar mass), it is crucial to further ac-count for the effect of dynamical friction (DF, Chandrasekhar 1943;Binney & Tremaine 2008).At the present time, state-of-the art numerical simulations withsub-grid recipes for unresolved physical processes promise to be thebest tool to properly model the orbital evolution of stars and MPswithin galaxies (e.g. Tremmel et al. 2015; Pfister et al. 2017, 2019).However, they cannot be adopted when one needs to perform athroughout exploration of the parameter space as they are typicallyvery computationally expensive. For this, different tools, such assemi-analytical approaches, may be preferred owing to their muchlower computational cost, that comes at the expense of some preci-sion.In semi-analytical models, the motion of a test particlesubject to the galactic potential can be followed more and moreaccurately as the description of the potential associated to eachgalactic component (the dark matter halo, disc, bulge etc.) getscloser to the one of real galaxies. While the spherically symmetricgalactic components can be approximated reasonably well by © a r X i v : . [ a s t r o - ph . GA ] O c t M. Bonetti et al. well-known potential-density pairs featuring a simple analyticalform (e.g. Navarro et al. 1997; Plummer 1911; Dehnen 1993), thesame cannot be said for less idealized galactic structures such asdiscs. In fact, their light distribution is often well described byan exponential disc density profile, which unfortunately does notadmit a close analytical form for the related potential, meaning thatits computation has to be performed numerically and the associatedcomputational cost may be relatively large. For this, less physicallymotivated disc potentials featuring a simpler, close analyticalform are often adopted in the literature, even if their approxima-tion of the potential of a physical galactic disc may be relatively poor.Even neglecting the aforementioned issues, a proper descrip-tion of the galactic potential in semi-analytical models is not enoughwhen the target particle is a MP, as DF may significantly affect theevolution of its orbit. DF arises as a response of the background tothe perturbing mass, and it typically results in a deceleration of thelatter, which gradually sinks towards the centre of its host system.Its standard semi-analytical description grounds its basis on theseminal work of Chandrasekhar (1943), who first analytically mod-elled the motion of an MP in an infinite and homogeneous stellarbackground. In spite of its simplicity, Chandrasekhar’s approachhas been proven to work remarkably well in many scenarios, some-times far beyond the expected regime of applicability (e.g. White1983), and despite missing the physics of resonant interactions be-tween the MP and the host system (Tremaine & Weinberg 1984;Weinberg 1986). Owing to its success, the semi-analytical treat-ment for DF has seen a substantial development throughout the lastdecades, as attested by the vast piece of recent literature focusing onimproving Chandrasekhar’s prescriptions and their agreement withN-body integrations (e.g. Hashimoto et al. 2003; Just & Peñarru-bia 2005; Just et al. 2011; Arca-Sedda & Capuzzo-Dolcetta 2014;Read et al. 2006; Petts et al. 2015, 2016). However, the work in thisdirection has been limited to very idealized, spherically symmetricand isotropic systems, which can only describe a narrow range ofgalaxy environments; in fact, astrophysical systems feature muchmore complexity than this: an extension of the semi-analytical DFtreatment to composite systems, possibly featuring disky compo-nents with rotational support, would thus greatly enhance the realmof applicability of semi-analytical DF.N-body experiments already addressed the effect of DF ina composite, rotationally supported environment; they highlightedthat DF acting in rotationally supported systems induces the circular-ization of a MP initially on a prograde, eccentric orbit, and reversesthe angular momentum of counter-rotating MPs, then again promot-ing circularization (Dotti et al. 2006; Callegari et al. 2011a; Fiacconiet al. 2013a); this effect appears to be independent of the nature ofthe background (Dotti et al. 2007). However, such past studies neverattempted a quantitative modelling of the phenomenon. Only veryrecently, Bonetti et al. (2020) successfully modelled for the firsttime the DF-induced “drag towards circular co-rotation” in a semi-analytical fashion in a multi-component galaxy that featured a darkmatter halo, a stellar disk and a bulge; their prescription showedremarkable agreement with N-body simulations down to the cir-cularization phase. However, their study was limited to on-planeorbits, and their approach can be adopted only before circulariza-tion occurs.Here, we expand the work presented in Bonetti et al. (2020) by Under similar assumptions, Ostriker (1999) describel the analogous phe-nomenon within a gaseous medium. improving the DF treatment in rotating discs in order to extend itsvalidity to circular orbits, but also providing a completely revisitedimplementation for the acceleration evaluation due to an exponentialdisc mass distribution, now allowing arbitrary orbits and not onlyequatorial motion. The paper is organised as follows: in Section 2we describe the improvements brought to the framework developedin Bonetti et al. (2020); in Section 3 we validate our updated setupagainst existing semi-analytical codes and full N-body simulationsfeaturing the same initial conditions. Finally, in Section 4 we drawour conclusions.
As in Bonetti et al. (2020), we consider a multi-component galaxymodel comprising a Navarro-Frenk-White (NFW, Navarro et al.1997) profile to model the dark matter (DM) halo, an Hernquist(Hernquist 1990) profile to describe a compact stellar bulge and athick exponential disc profile to characterise the galactic disc. Ontop of that we model the effect of DF on the motion of massiveperturbers through the employment of dissipative forces.In this section we detail the procedure we followed to i) evaluatethe potential and accelerations generated by a thick exponential discin the general case, ii) discuss possible improvement for the DF forcein spherically symmetric systems and iii) extend the validity of ourprescription for DF in rotating discs to nearly circular orbits.
We consider a finite-thickness exponential disc whose mass densityis described by 𝜌 𝑑 ( 𝑅, 𝑧 ) = 𝑀 𝑑 𝜋𝑅 𝑑 𝑧 𝑑 e − 𝑅 / 𝑅 𝑑 sech (cid:18) 𝑧𝑧 𝑑 (cid:19) , (1)where 𝑀 𝑑 , 𝑅 𝑑 , 𝑧 𝑑 are the total mass, scale radius and scale height ofthe disc, respectively. Such profile is our choice for the modelling ofgeneral galactic discs, since it can well reproduce the light profilesof real galaxies and it is physically motivated, being the verticaldistribution obtained for an isothermal self-gravitating disc (Spitzer1942). For these reasons, it is usually employed to generate the ini-tial conditions of N-body simulations (see e.g. Yurin & Springel2014), allowing for a validation of the results of our semi-analyticalprescriptions against N-body runs. Nevertheless, the double expo-nential disc, i.e. a disc with vertical density profile decaying ase −| 𝑧 |/ 𝑧 𝑑 , is widely used and well reproduces the vertical profiles ofedge-on disc galaxies (e.g. de Grijs et al. 1997). In order to have aproper comparison, we also implemented this second profile in ourframework. In general, to obtain the gravitational potential from a givendensity distribution, one has to solve the Poisson’s equation. Prac-tically, when deviating from spherical symmetry, a complete an-alytical expression can be found only in a handful of cases andunfortunately the thick exponential case does not belong to those.Therefore, to obtain the gravitational potential (and consequentlythe accelerations) of such mass distribution at arbitrary ( 𝑅, 𝑧 ) pairs, The double exponential disc also does not admit any closed form forthe gravitational potential, but requires numerical integration as given inequation (2), where the sole difference is in the function 𝐼 ( 𝑘 ) given byequation A9 of Kuijken & Gilmore (1989). MNRAS , 1–15 (2020) P evolution in realistic galaxy models we necessarily need to recur to numerical methods. Unfortunately,the general solution of the Poisson’s equation that can be obtainedvia Green’s function results in multidimensional integrals, implyingthat obtaining the solution can be pretty expensive in computationalterms. Thus, if possible, it is advisable to reduce the expressions asmuch as possible before recurring to numerical integration.It turns out that, when considering an axi-symmetric densitydistribution, the Poisson’s equation can be solved in terms of a Han-kel’s transform (see e.g. Kuijken & Gilmore 1989). This approachallows us to write the gravitational potential as a 1-dimensionalintegral (see Appendix A for a complete derivation), given by 𝜙 ( 𝑅, 𝑧 ) = − 𝐺 𝑀 𝑑 𝑅 𝑑 𝑧 𝑑 ∫ ∞ d 𝑘𝐽 ( 𝑘 𝑅 ) 𝐼 𝑧 ( 𝑘 ) (cid:16) 𝑅 − 𝑑 + 𝑘 (cid:17) / , (2)where 𝐽 is the Bessel function of the first kind, and 𝐼 𝑧 ( 𝑘 ) is a func-tion depending on both 𝑘 and 𝑧 through the Gauss hypergeometricfunction 𝐹 (see Appendix A and B).The radial and vertical accelerations can be readily obtainedby taking the (negative) derivative of equation (2) with respect to 𝑅 and 𝑧 and then performing the integration over the 𝑘 variable 𝑎 𝑅 = − 𝜕𝜙𝜕𝑅 = − 𝐺 𝑀 𝑑 𝑅 𝑑 𝑧 𝑑 ∫ ∞ d 𝑘 𝑘𝐽 ( 𝑘 𝑅 ) 𝐼 𝑧 ( 𝑘 ) (cid:16) 𝑅 − 𝑑 + 𝑘 (cid:17) / , (3) 𝑎 𝑧 = − 𝜕𝜙𝜕𝑧 = − 𝐺 𝑀 𝑑 𝑅 𝑑 𝑧 𝑑 ∫ ∞ d 𝑘𝐽 ( 𝑘 𝑅 ) − 𝜕 𝑧 𝐼 𝑧 ( 𝑘 ) (cid:16) 𝑅 − 𝑑 + 𝑘 (cid:17) / . (4)Note that the integration of Bessel functions could be problematicgiven their highly oscillatory behaviour. Nevertheless, if appropriatequadrature techniques are employed, the integrals can be computedquite efficiently using a few hundreds of points, still guaranteeing agood level of accuracy (see Appendix A). DF, though not dissipative in nature, in the context of semi-analyticalcalculations is usually modelled as a drag force acting opposite tothe velocity of a MP. Despite the effort to improve the originalderivation by Chandrasekhar (1943) (see e.g. Tremaine & Weinberg1984; Weinberg 1986; Mulder 1983; Colpi et al. 1999, and referencetherein) the rather simple original expression, i.e. a df = − 𝜋𝐺 ln ( + Λ ) 𝑚 𝑝 𝜌 (cid:32) erf ( 𝑋 ) − 𝑋 e − 𝑋 √ 𝜋 (cid:33) v 𝑝 | v 𝑝 | , (5)is still widely used, mostly because of its simplicity. In the aboveformula, 𝑚 𝑝 and v 𝑝 are the MP mass and velocity, 𝜌 is the back-ground density, Λ = 𝑝 max / 𝑝 min is the ratio between the maximumand minimum impact parameter, while 𝑋 = 𝑣 𝑝 /(√ 𝜎 ) is the ratio ofthe MP velocity over the velocity dispersion. Though systematicallyemployed, equation (5) is subjected to a number of assumptions thatare not always justified in realistic galactic contexts, i.e.:1. the motion of scattering particles with large impact parametersis rectilinear;2. the background density is homogeneous and isotropic;3. the velocity distribution function is isotropic and Maxwellianwith constant dispersion 𝜎 ;4. the maximum impact parameter is chosen equal to a maxi-mum scale-length of the system, while 𝑝 min is kept fixed for all possible encounters and equal to a characteristic scale length (or, inthe case of point-like perturbers, the radius below which a 90 ◦ de-flection occurs). Specifically, this last assumption implies that onlyencounters with velocity lower than the MP velocity contribute tothe deceleration.Despite the above assumptions are clearly not fully compatiblewith realistic galaxies, the introduction of some a-posteriori mod-ifications (mostly affecting points 2 and 3 above) to equation (5)allows to reproduce the results from N-body simulations with rea-sonable accuracy. In particular, when limiting the DF treatment tospherically symmetric profiles (which model to a fairly good ap-proximation elliptical galaxies or galaxy bulges), the adoption ofthe position-dependent mass density, velocity dispersion and min-imum and maximum impact parameters substantially improves theagreement between semi-analytical models and full 𝑁 -body sim-ulations (see e.g. Hashimoto et al. 2003; Just & Peñarrubia 2005;Just et al. 2011; Petts et al. 2015, 2016). Following Bonetti et al.(2020), we implement equation (5) adopting the position-dependentvelocity dispersion 𝜎 ( 𝑟 ) and mass density 𝜌 ( 𝑟 ) , and the followingexpressions for the maximum and minimum impact parameters: 𝑝 max = 𝑟 / 𝛾,𝑝 min = max (cid:32) 𝐺𝑚 𝑝 𝑣 𝑝 + 𝜎 ( 𝑟 ) , 𝐷 𝑝 (cid:33) ,𝛾 = − d ln 𝜌 d ln 𝑟 , (6)with 𝛾 denoting the logarithmic slope of the density profile, 𝑟 theradial coordinate and 𝐷 𝑝 the physical radius of the MP (which canbe set to zero for the case of an MBH).Even with the above modifications, the acceleration computedwith the functional form of equation (5) considers only the contri-bution of scattering particles moving slower than the MP. Althoughthis does not represent a serious issue for the majority of stellarsystems, the deceleration described by equation (5) could severelyunderestimate the real DF when the number of background objectsmoving faster than the MP is large. Specifically, Antonini & Merritt(2012) investigated the DF experienced by a MP in galactic coresdominated by a central MBH and they found that, when the slopeof the density profile is 𝛾 (cid:46)
1, the contribution of the so-called“fast-moving stars” can become comparable to or even larger thanthat of slow moving ones.In order to take into account such fast population, one needsto drop some assumptions and consider the general Chandrasekharformula (see equations 25 and 26 of Chandrasekhar 1943), thatdespite being more general, usually does not allow for an analyticalform, i.e. a df = − 𝜋𝐺 𝑚 𝑝 𝜌 ( 𝑟 ) v 𝑝 | v 𝑝 | ×× ∫ 𝑣 esc 𝐽 ( 𝑣 𝑝 , 𝑣 ★ , 𝑝 max ) 𝜋𝑣 ★ 𝜋𝑣 ★ 𝑓 ( 𝑣 ★ ) d 𝑣 ★ , (7) The integral over the relative velocity 𝐽 ( 𝑣 𝑝 , 𝑣 ★ , 𝑝 max ) actually has ananalytical form. It is usually quoted in an approximate form that is wellbehaved everywhere except when 𝑣 𝑝 and 𝑣 ★ are close (see e.g. equation 7of Antonini & Merritt 2012). For completeness we report the full expressionin Appendix C.MNRAS000
1, the contribution of the so-called“fast-moving stars” can become comparable to or even larger thanthat of slow moving ones.In order to take into account such fast population, one needsto drop some assumptions and consider the general Chandrasekharformula (see equations 25 and 26 of Chandrasekhar 1943), thatdespite being more general, usually does not allow for an analyticalform, i.e. a df = − 𝜋𝐺 𝑚 𝑝 𝜌 ( 𝑟 ) v 𝑝 | v 𝑝 | ×× ∫ 𝑣 esc 𝐽 ( 𝑣 𝑝 , 𝑣 ★ , 𝑝 max ) 𝜋𝑣 ★ 𝜋𝑣 ★ 𝑓 ( 𝑣 ★ ) d 𝑣 ★ , (7) The integral over the relative velocity 𝐽 ( 𝑣 𝑝 , 𝑣 ★ , 𝑝 max ) actually has ananalytical form. It is usually quoted in an approximate form that is wellbehaved everywhere except when 𝑣 𝑝 and 𝑣 ★ are close (see e.g. equation 7of Antonini & Merritt 2012). For completeness we report the full expressionin Appendix C.MNRAS000 , 1–15 (2020) M. Bonetti et al. 𝐽 ( 𝑣 𝑝 , 𝑣 ★ , 𝑝 max ) = ∫ 𝑣 𝑝 + 𝑣 ★ | 𝑣 𝑝 − 𝑣 ★ | d 𝑉 (cid:32) + 𝑣 𝑝 − 𝑣 ★ 𝑉 (cid:33) ln (cid:32) + 𝑝 𝑉 𝐺 𝑚 𝑝 (cid:33) . (8)Equation (7) inevitably implies a higher computational bur-den due to the unavoidable numerical computation of the integralover the distribution function. Given that in the vast majority ofastrophysical situations the inclusion of fast moving stars representsonly a minor correction (but see Section 4 for some caveats of ourchoice), in this work we adopt the simplified expression consideringonly slow moving stars, together with the modifications introducedby different authors that improve its range of applicability. Bonetti et al. (2020) derived a prescription to describe the DFacceleration acting on a MP determined by a rotating disc, a df , disc = − 𝐴 disc 𝜋𝐺 ln ( + Λ ) 𝑚 𝑝 𝜌 𝑑 ( 𝑅, ) v 𝑝 − v 𝑐 ( 𝑅 )| v 𝑝 − v 𝑐 ( 𝑅 )| , (9)where v 𝑝 , 𝑚 𝑝 are the velocity and mass of the MP, v 𝑐 ( 𝑅 ) is thecircular velocity at radius 𝑅 determined by the total gravitational po-tential, while Λ = 𝑝 max , d / 𝑝 min , d , that enters an effective Coulomblogarithm, is chosen as the ratio of the maximum impact parameter,taken equal to the scale-height of the disc, and a minimum one givenby 𝑝 min , d = 𝐺𝑚 𝑝 𝑣 + . 𝑣 𝑐 ,𝑣 rel = | v 𝑝 − v 𝑐 | . (10)Finally, the factor 𝐴 disc is a tunable constant set by comparing thesemi-analytical orbital evolution of the MP with that obtained in anN-body simulation. Bonetti et al. (2020) found that this factor is ofthe order of unity, enforcing that their modeling catches the mainphysics of DF in discs.The main limitation of this prescription lies on the dependenceof the DF acceleration on the relative velocity with respect to cir-cular velocity. When the MP velocity approaches the circular one,equation (9) diverges, thus the orbital decay cannot be followedfurther.To overcome this limitation, here we consider a less idealizeddistribution function for the stars in the disc. Specifically, we replacethe delta function with an isotropic Gaussian centered around therelative velocity between the MP and the local rotational motion.This choice allows us to employ the very same functional form ofthe standard DF Chandrasekhar formulation, but with the velocitydependence on the perturber velocity in the galactic frame replacedby the relative velocity with respect to the local medium (see e.g.Kashlinsky 1986, for a similar derivation). In particular, we obtain a df , disc = − 𝜋𝐺 ln ( + Λ ) 𝑚 𝑝 𝜌 𝑑 ( 𝑅, 𝑧 ) ×× (cid:32) erf ( 𝑋 ) − 𝑋 e − 𝑋 √ 𝜋 (cid:33) v rel | v rel | , (11)where v rel = v 𝑝 − v rot ( 𝑅 ) , while 𝑋 = | v rel |/(√ 𝜎 𝑅 ) , with 𝜎 𝑅 denoting the radial velocity dispersion. For the Coulomb logarithm,we replace the expression of 𝑝 min , d appearing in equation (10) by 𝑝 min , d = 𝐺𝑚 𝑝 𝑣 + 𝜎 𝑅 , (12) Halo Bulge DiscMass 1 . × M (cid:12) . × M (cid:12) . × M (cid:12) Scale Radius 37 ( ) kpc 0 .
96 kpc 4 .
25 kpcScale Height - - 0 .
85 kpcProfile Hernquist (NFW) Hernquist Exponential 𝑁 part × × 𝜀
40 pc 10 pc 10 pc
Table 1.
Structural parameters of the employed galaxy model. For eachcomponent, 𝑁 part and 𝜀 correspond to the number of particles and thePlummer-equivalent gravitational softening employed in the N-body runs. which takes into account both the relative velocity with surround-ing medium, that can be different from the circular velocity, andan intrinsic dispersion, which sets the minimum effective distancecharacterising the encounters between the MP and the objects in thedisc. The presence of a velocity dispersion 𝜎 𝑅 implies that the discis not fully rotationally supported. We can characterise the rotationalvelocity by following Hernquist (1993), | v rot | = 𝑣 𝑐 ( 𝑅 ) + 𝜎 𝑅 ( 𝑅 ) (cid:18) − 𝜅 Ω − 𝑅𝑅 𝑑 (cid:19) , (13)in which 𝑣 𝑐 denotes the circular velocity at 𝑅 , while 𝜅 and Ω (bothdepending on 𝑅 ) are the epicyclic and angular frequency, respec-tively. In evaluating such quantities for multi-component galaxymodels, we consider the total potential rather than that generatedby the disc only. As in Hernquist (1993), we also assume that thevelocity dispersion changes as 𝜎 𝑅 ( 𝑅 ) ∝ e − 𝑅 / 𝑅 𝑑 , (14)where the normalisation is evaluated assuming equality between 𝜎 𝑅 and the critical radial velocity dispersion (augmented by a factor 𝑄 ) at a specific radius, here 2 𝑅 𝑑 , i.e. 𝜎 𝑅, crit = 𝑄 × . 𝐺 Σ ( 𝑅 𝑑 ) 𝜅 ( 𝑅 𝑑 ) . (15)We can note that since in equation (13) the term in parenthe-sis on the right hand side is usually negative, then the rotationalvelocity results (often only slightly) lower than the local circularvelocity. However, as already pointed out in Hernquist (1993), atsmall radii the approximations yielding equation (13) can breakdown and provide nonphysical 𝑣 rot (e.g. an imaginary velocity). Toavoid unnecessary complications we decided to fix the rotationalvelocity equal to a fraction of the local circular velocity (specifi-cally 0 . 𝑣 𝑐 ( 𝑅 ) ) anytime equation (13) gives nonphysical results.We expect our choice not to dramatically impact the evolution, sinceat small radii the dynamics is likely dominated by the stellar bulge(and the associated DF), while the contribution of the disc should besubdominant. We anyway try to address this issues by also directlyextracting the rotational profile from an N-body realisation, as wewill further discuss in Section 3. In this section, we compare our semi-analytical setup with othersimilar implementations for the disc potential and the DF treatment, Here we assume 𝑄 = .
5. MNRAS , 1–15 (2020)
P evolution in realistic galaxy models − − − R/R d − − − − − − l og z / z d − − − − − relative error a R − − − − − r e l a t i v ee rr o r − − − R/R d − − − − − − l og z / z d − − − − − − relative error a z − − − − − r e l a t i v ee rr o r − − − R/R d − − − − − − l og z / z d relative error a R − − − − − − r e l a t i v ee rr o r − − − R/R d − − − − − − l og z / z d relative error a z − − − − − − r e l a t i v ee rr o r Figure 1.
Relative error for the radial (left column) and vertical (right column) accelerations obtained with our code. Upper panels: relative errors are evaluatedagainst a numerical integration of equations (3) and (4) with 2 × points using trapezoidal rule in the range 𝑘 ∈ [ , ] . Lower panels: relative errors areobtained comparing our acceleration with those computed by the software AGAMA employing an exponential disc with identical parameters. and against full N-body simulations. We start by considering theconservative dynamics in exponential discs, highlighting how adifferent vertical disc profile can affect the orbits. We then analyseour DF prescriptions in multi-component galaxy models and wecompare them against the results from N-body simulations featuringthe same physical system and MP orbital initial conditions. Weemploy the public code GIZMO (Hopkins 2015) to run the N-bodysimulations, using the same N-body setup of Bonetti et al. (2020),whose parameters are reported in Table 1 for completeness. Here we investigate the conservative dynamics in exponential discs.In our framework we consider both the exponential disc with 𝑧 -profile decaying as sech ( 𝑧 / 𝑧 𝑑 ) (see Section 2.1) as well as theso-called double exponential disc, i.e. a profile also decaying as anexponential along 𝑧 . In Fig. 1, we report the relative error in the radial (left column) andvertical (right column) accelerations of an exponential disc with asech ( 𝑧 / 𝑧 𝑑 ) profile, computed by our semi-analytical frameworkthrough a modified version of the double exponential quadraturetechnique (see e.g. Michalski & Mosig 2016). In the upper panels,we compare our calculation to the numerical evaluation of equa-tions (3) and (4) obtained via the trapezoidal rule with ∼ × points and assuming a 𝑘 max ≈ . We can appreciate that therelative errors are always very small (10 − –10 − ) over severaldecades in both 𝑅 and 𝑧 . Only well above 10 𝑅 / 𝑅 𝑑 the relative errorstarts to grow significantly, but still remaining below 10 − . Thistrend is associated with the behaviour of the Bessel function in theintegrand function, that oscillates significantly for large 𝑅 , mildlyreducing the accuracy of the numerical integration.In the lower panels, instead, we evaluate the relative errorsagainst the accelerations obtained from the software AGAMA(Vasiliev 2019), in which the disc density profile is expanded in mul- MNRAS000
Relative error for the radial (left column) and vertical (right column) accelerations obtained with our code. Upper panels: relative errors are evaluatedagainst a numerical integration of equations (3) and (4) with 2 × points using trapezoidal rule in the range 𝑘 ∈ [ , ] . Lower panels: relative errors areobtained comparing our acceleration with those computed by the software AGAMA employing an exponential disc with identical parameters. and against full N-body simulations. We start by considering theconservative dynamics in exponential discs, highlighting how adifferent vertical disc profile can affect the orbits. We then analyseour DF prescriptions in multi-component galaxy models and wecompare them against the results from N-body simulations featuringthe same physical system and MP orbital initial conditions. Weemploy the public code GIZMO (Hopkins 2015) to run the N-bodysimulations, using the same N-body setup of Bonetti et al. (2020),whose parameters are reported in Table 1 for completeness. Here we investigate the conservative dynamics in exponential discs.In our framework we consider both the exponential disc with 𝑧 -profile decaying as sech ( 𝑧 / 𝑧 𝑑 ) (see Section 2.1) as well as theso-called double exponential disc, i.e. a profile also decaying as anexponential along 𝑧 . In Fig. 1, we report the relative error in the radial (left column) andvertical (right column) accelerations of an exponential disc with asech ( 𝑧 / 𝑧 𝑑 ) profile, computed by our semi-analytical frameworkthrough a modified version of the double exponential quadraturetechnique (see e.g. Michalski & Mosig 2016). In the upper panels,we compare our calculation to the numerical evaluation of equa-tions (3) and (4) obtained via the trapezoidal rule with ∼ × points and assuming a 𝑘 max ≈ . We can appreciate that therelative errors are always very small (10 − –10 − ) over severaldecades in both 𝑅 and 𝑧 . Only well above 10 𝑅 / 𝑅 𝑑 the relative errorstarts to grow significantly, but still remaining below 10 − . Thistrend is associated with the behaviour of the Bessel function in theintegrand function, that oscillates significantly for large 𝑅 , mildlyreducing the accuracy of the numerical integration.In the lower panels, instead, we evaluate the relative errorsagainst the accelerations obtained from the software AGAMA(Vasiliev 2019), in which the disc density profile is expanded in mul- MNRAS000 , 1–15 (2020)
M. Bonetti et al. − y [ k p c ] − x [kpc] − z [ k p c ] − y [kpc] this workgalpy . . . . . × − − − − − − ∆ r / r Figure 2.
Left panel: orbit projections in the 𝑥 − 𝑦 , 𝑥 − 𝑧 and 𝑦 − 𝑧 planes when considering a sech exponential disc with mass 𝑀 𝑑 = . × M (cid:12) , scaleradius 𝑅 𝑑 = .
25 kpc and scale height 𝑧 𝑑 = .
85 kpc. We show the comparison between the trajectories obtained within our framework (solid blue lines) andthose obtained with the software galpy (dashed orange lines). The initial velocity is set equal to 0 . 𝑣 𝑐 in the both 𝑦 and 𝑧 directions, with 𝑣 𝑐 the circularvelocity at a cylindrical radius 𝑅 = − after 2 Gyr of evolution. tipoles and the potential (and the corresponding force) is computedvia the Poisson’s equation for each term of the multipole expansion.In this case, we find a larger error, but still at the level of 10 − fora vast portion of the ( 𝑅 − 𝑧 ) plane, while larger deviations arise atthe borders of the plane. Given that the computation in AGAMA re-lies on the multipole expansion and makes systematic use of splineinterpolation, the quantities that AGAMA evaluates are necessarilyapproximated and subjected to some numerical noise. Still, errorsare quite small in the relevant portion of the ( 𝑅 − 𝑧 ) plane and reach10 − only for 𝑅 and 𝑧 that are much larger or much smaller than 𝑅 𝑑 and 𝑧 𝑑 . This does not represent a real problem as in realisticgalaxy models other galactic components are expected to dominatethe dynamics at large/small distances, for instance dark matter halos(at large distances) or stellar bulges (close to the center), hence wecan consider our implementation quite robust. We next compare the orbital evolution that we obtain with thoseof the software galpy (Bovy 2015). Specifically, we analyse thetrajectories obtained in an exponential disc with mass 𝑀 𝑑 = . × M (cid:12) , scale radius 𝑅 𝑑 = .
25 kpc and scale height 𝑧 𝑑 = .
85 kpc and with a vertical profile given by either sech ( 𝑧 / 𝑧 𝑑 ) or e −| 𝑧 |/ 𝑧 𝑑 . In the left panel of Fig. 2, we compare the orbit for thesech case (blue solid for our implementation vs orange dashed forgalpy) by showing the trajectory projections in the 𝑥 − 𝑦 , 𝑥 − 𝑧 , For the comparison in Fig. 1 we adopt within AGAMA a multipole ex-pansion ranging between 10 − 𝑅 𝑑 and 500 𝑅 𝑑 ; we used a logarithmicallyspaced grid in 𝑅 with 1000 grid points and maximum order of the expan-sion 𝑙 max = 𝑧 grid are kept as default inAGAMA (see Vasiliev 2019, for more details). The sech exponential disc within the galpy framework is accessiblethrough the SCF basis-function-expansion (see galpy documentation). 𝑦 − 𝑧 planes. The orbit is evolved for approximately 2 Gyr start-ing from an initial radial separation of 𝑅 = 𝑦 and 𝑧 directions equal to 0.5 times the local cir-cular velocity. From the left panel, we can infer that the orbit inour implementation closely matches that of galpy and to betterquantify the agreement, in the right panel of the figure, we evalu-ate the relative error on the spherical radius between the two runs,finding that is at the level of 10 − at most. This small difference isprobably due to the different approaches followed to compute thedisc potential: specifically through the employment of the Hankel’stransform in our case, which reduces the problem to the numericalcomputation of an integral (see equation 2), and instead through amultipole expansion approximation in the galpy framework. Weperform the same comparison by also considering the double expo-nential disc profile. As for the sech we found a good agreement,therefore enforcing the robustness of our framework. Another im-portant comparison with galpy concerns the computational cost,that for a semi-analytical framework has to necessarily be modest inorder to fully exploit its power. For each run, we find a total run-timeof the order of few tens of seconds (on a single core of a standardlaptop), similar to that of galpy (but only once an adequate accu-racy, similar to that in our framework, is also adopted in galpy).Moreover, to further speed-up the orbital integration, we have alsoimplemented the possibility of storing in advance the accelerationscomputed on a ( 𝑅 − 𝑧 ) adaptive grid and then employ a bi-cubicinterpolation to obtain 𝑎 𝑅 and 𝑎 𝑧 on arbitrary points.Finally, we compare the orbital evolution of a point mass whensubject to the potential of a double exponential disc and a sech one.We consider the same disc and initial conditions as described above,but we also add a velocity component in the 𝑧 direction equal to 0.1times the local circular velocity. The results are shown in Fig. 3.In both the left and right panels we show the three projectionsover the 𝑥 − 𝑦 , 𝑥 − 𝑧 and 𝑦 − 𝑧 planes as sub-panels (from top leftto bottom right). We can appreciate that the orbits resemble eachother, at least qualitatively, in the sense that they cover more or MNRAS , 1–15 (2020)
P evolution in realistic galaxy models − y [ k p c ] − x [kpc] − . . . z [ k p c ] z profile: e − z/z d z profile: sech ( z/z d ) − y [kpc] − y [ k p c ] −
10 0 10 x [kpc] − . . . z [ k p c ] z profile: e − z/z d z profile: sech ( z/z d ) −
10 0 10 y [kpc] Figure 3.
Comparison between orbits in a double exponential disc (solid blue line) and in an exponential disc with a 𝑧 -profile described by the hyperbolicsecant (dashed orange lines). In the left (right) panel the in-plane initial velocity is set equal to 0.5 (1.5) the local circular velocity, while vertical velocity is setto 0.1 this value. r [ k p c ] SA N-body . . . L z / | L z , | × e r [ k p c ] SA N-body − L z / | L z , | × e Figure 4.
Time evolution of the radial separation, the normalised angular momentum, and the eccentricity obtained with the semi-analytical framework (SA,blue solid lines) and an N-body simulation (orange dashed lines), employing the same initial conditions and physical parameters. Left panel: prograde orbit.Right panel: initially retrograde orbit. less the same volume of the disc. Still, from a quantitative pointof view, e.g. checking the position reached at a specific time, theorbits are well distinct and evolve on slightly different timescales.Hence, a precise reconstruction of orbits in a realistic galaxy requiresthe employment of the density/potential pair that better describessuch structure, that, given the variety of profiles observed, wouldrequire an ad-hoc model for each individual case. Such exercisecan be easily preformed using our procedure, that provides bothdouble exponential and sech exponential disc implementations,but is beyond the scope of our current validation study. Now, we consider the MP evolution under the combined DF effect ofall galactic components. In Fig. 4, we show, from top to bottom, howthe radial separation, the 𝑧 component of the angular momentum(normalised to its initial value), and the eccentricity of a MP orbitingin the galactic mid-plane evolve. We here compare the evolutionobtained with our semi-analytical setup, as described in Section 2(solid lines), and with N-body simulations (dashed lines). Fromthe figure, we can infer that our setup can qualitatively catch therelevant phenomenology of the evolution, such as the trends shownby the angular momentum and eccentricity, but the timescales are MNRAS000
Time evolution of the radial separation, the normalised angular momentum, and the eccentricity obtained with the semi-analytical framework (SA,blue solid lines) and an N-body simulation (orange dashed lines), employing the same initial conditions and physical parameters. Left panel: prograde orbit.Right panel: initially retrograde orbit. less the same volume of the disc. Still, from a quantitative pointof view, e.g. checking the position reached at a specific time, theorbits are well distinct and evolve on slightly different timescales.Hence, a precise reconstruction of orbits in a realistic galaxy requiresthe employment of the density/potential pair that better describessuch structure, that, given the variety of profiles observed, wouldrequire an ad-hoc model for each individual case. Such exercisecan be easily preformed using our procedure, that provides bothdouble exponential and sech exponential disc implementations,but is beyond the scope of our current validation study. Now, we consider the MP evolution under the combined DF effect ofall galactic components. In Fig. 4, we show, from top to bottom, howthe radial separation, the 𝑧 component of the angular momentum(normalised to its initial value), and the eccentricity of a MP orbitingin the galactic mid-plane evolve. We here compare the evolutionobtained with our semi-analytical setup, as described in Section 2(solid lines), and with N-body simulations (dashed lines). Fromthe figure, we can infer that our setup can qualitatively catch therelevant phenomenology of the evolution, such as the trends shownby the angular momentum and eccentricity, but the timescales are MNRAS000 , 1–15 (2020)
M. Bonetti et al. r [ k p c ] SA, soft N-body SA soft, N-body v rot . . . L z / | L z , | × e r [ k p c ] SA, soft N-body SA soft, N-body v rot − L z / | L z , | × e Figure 5.
Same as Fig. 4, but modifying the semi-analytical setup to take into account the finite spatial resolution of the N-body simulations (solid blue lines)as well as the rotational velocity in the N-body run (solid green lines). Specifically, the minimum impact parameter 𝑝 min for each component (disc, halo, bulge)is fixed at 2.8 times the respective Plummer-equivalent softening length, while the evolution described by the green lines also employs the rotational velocitydirectly extracted from the N-body simulations. not exactly reproduced. Indeed, although the evolution of a progradeorbit (left panel) seems to be well recovered by our framework,the same does not hold when looking at an initially retrogradetrajectory (right panel), in which the semi-analytical decay looksfaster. Such differences, however, can be interpreted by remindingthat also the N-body simulations have some limitations and a faircomparison with our semi-analytical framework has to necessarilytake into account the limited spatial and mass resolution of N-bodyruns. Specifically, we need to consider that in N-body codes allinteractions below the softening length get considerably weakenedand, as a consequence, the resulting DF force acting on a MP issmaller with respect to an ideal and arbitrarily resolved physicalsystem (see, e.g. Tremmel et al. 2015; Pfister et al. 2017).To account for the limited N-body resolution we can act on theparameters that enter the expression of the DF force. According toequation (12), we can observe that the minimum impact parameter 𝑝 min depends on the MP velocity and specifically for the disc DF itdepends on the relative velocity between the MP and the surround-ing medium. A small 𝑝 min translates into larger DF accelerationsthrough the term ln ( + Λ ) . For the retrograde case, 𝑣 rel can be quitelarge, therefore 𝑝 min can decrease significantly, becoming smallerthan the softening length of the N-body simulation. To assess theinfluence of the softening in the N-body framework we thus forcethe minimum impact parameter 𝑝 min to be larger than 2.8 times thePlummer-equivalent softening length of the corresponding particletype, i.e. halo, bulge or disc. The evolution resulting from softenedsemi-analytical setup is shown in Fig. 5 with blue lines. Both the pro-grade and retrograde orbits (left and right panel respectively) nowshow a slower decay, and in particular the retrograde case appearsto be much more similar to the N-body evolution. The progradeinspiral however seems now too slow.Another important effect that needs to be considered is the pos-sible difference between the rotational velocity profile employedin our semi-analytic approach, based on the approximated equa-tion (13) (Hernquist 1993), and that of the simulation. As shown inFig. 6, a direct comparison of the two velocity profiles reveal that, R [kpc]050100150200 v [ k m / s ] v c v rot , Hernquist v rot , N − body analytical fit Figure 6.
Circular velocity (thick blue line), rotational velocity as predictedby equation (13) (thin orange line) and extracted from N-body simulation(red crossed line) as a function of radius. At small 𝑅 , the large oscillations ofN-body 𝑣 rot are due to the limited particle coverage in the central regions. Toavoid spurious numerical noise we employed an analytical fit of the N-body 𝑣 rot (black dashed line). as we move to smaller radii, the difference in the rotational velocityprofile between the one used in the semi-analaytical approach andthat in the N-body simulation gets enhanced. If the rotational ve-locity is significantly different from the local circular one when theorbit becomes nearly circular (i.e. above 400 Myr in the progradecase), then the disc DF force does not vanish (see equations 11 and13), resulting in a slightly faster decay. In order to account for this MNRAS , 1–15 (2020)
P evolution in realistic galaxy models r [ k p c ] SA, soft N-body SA soft, N-body v rot L z / | L z , | .
00 0 .
25 0 .
50 0 .
75 1 . × ι [ d e g ] r [ k p c ] SA, soft N-body SA soft, N-body v rot − L z / | L z , | .
00 0 .
25 0 .
50 0 .
75 1 . × ι [ d e g ] r [ k p c ] SA, soft N-body SA soft, N-body v rot L t o t / | L t o t , | .
00 0 .
25 0 .
50 0 .
75 1 . × ι [ d e g ] Figure 7.
Same as Fig. 5, but considering three different initial inclinations, 45 ◦ (upper left), 135 ◦ (upper right), 90 ◦ (bottom). The bottom panel of each caseshows the inclination evolution, while for the 𝜄 = ◦ run the total angular momentum is shown instead of its 𝑧 component. effect, we implemented the rotational velocity profile extracted fromthe simulation and fitted with a rational function of the radius (seethe dashed black line in Fig. 6) in our semi-analytical framework.This modification is enough to improve the agreement between thethe two approaches, as shown by the green lines in Fig. 5. Moreover,the apparent good agreement between the N-body run and the semi-analytical approach in the left panel of Fig. 4 actually results froma fortunate coincidence related to the different weighting of the DFacceleration generated by each galactic components, specifically aslightly stronger (because not softened) DF from the spherical com-ponents and a weaker DF from the disc, due to a velocity profilecloser to the circular one.Although the orbits we have considered so far were co-planarwith the disc, our semi-analytical framework allows us to exploreorbits with arbitrary inclinations (see Section 2.1 for details). Herewe check the orbital decay of MPs starting from off-plane con-figurations, as reported in in Fig. 7, 𝜄 = ◦ (upper left panels), 𝜄 = ◦ (upper right panels) and 𝜄 = ◦ (bottom panels). For each case, instead of the eccentricity, we show here the evolution of theinclination. For the 𝜄 = ◦ case we also report the total angularmomentum evolution, being 𝐿 𝑧 ≈
0. As in Fig. 5, for each incli-nation we compare the evolution derived from N-body simulations(orange dashed lines) with that computed in our semi-analyticalframework, again taking care of including the softening of the cor-responding N-body run (solid blue lines) and the N-body rotationalprofile (solid green lines). From the figure, we can clearly see thatthe semi-analytical evolution matches remarkably well the N-bodyone, correctly reproducing both the observed phenomenology andthe decay timescale. Minor differences are of course present, butthese are unavoidable given the discrete nature and the limited massresolution of the N-body approach. Similarly to the planar counter-parts, also in these inclined runs a better agreement is met when boththe softening and the rotational curve of the N-body are considered.We stress again that the very good agreement that we observe inboth planar and non-planar cases is achieved only when some inputfrom the N-body runs is provided. Indeed, a fair comparison has to
MNRAS000
MNRAS000 , 1–15 (2020) M. Bonetti et al. necessarily take into into account the limited resolution of N-bodysimulations. Still, for realistic galaxies we expect that a faster decay,similar to the one we find when considering our standard frame-work (see Fig. 4), is actually more physical and should provide abetter timescale estimate, valid at least down to radial separationsnot much smaller than the bulge scale radius. Indeed, well insidethis scale radius we expect the correction from fast moving star tostart being relevant, hence a more general DF formulation should beused. Despite this, our framework has proved to excellently capturethe physics driving MP motion in realistic multi-component galaxymodels with rotationally supported structures. Moreover, given itsquite modest computational cost compared to full N-body simula-tions, our framework represents a particularly attractive solution forlarge parameters space explorations in terms of galaxy structures aswell as MP masses and orbital configurations.
In this paper, we developed a semi-analytical framework to integratethe motion of particles in realistic galactic potentials composedof multiple components. Among several implemented potential-density pairs (spherical, axis-symmetric and triaxial), we includedthe cases of exponential discs with vertical density decaying eitheras e −| 𝑧 |/ 𝑧 𝑑 or sech ( 𝑧 / 𝑧 𝑑 ) , for which we developed an efficientsetup to numerically compute its potential and, consequently, theacceleration field. Though more complex to tackle, these exponen-tial profiles have to be preferred to other analytical profiles, that canonly partially reproduce the exponential decay of the surface bright-ness of observed disc galaxies. For example, the Miyamoto-Nagai(MN, Miyamoto & Nagai 1975) disc is one of the most commonanalytical profiles employed to model the potential of disc galaxies.Still, this profile can significantly differ than the one of realisticgalaxies, especially at large radii where it falls off as a power lawrather than exponentially (see discussion in chapter 2 of Binney &Tremaine 2008). To enhance the density decay, a combination ofthree MN discs is usually employed, although this requires one ofthe discs to have a negative mass, often leading to the appearance ofnon-physical negative density regions in the ( 𝑅 − 𝑧 ) plane (see e.g.Barros et al. 2016; Bonetti et al. 2019). This strongly supports ourchoice about the direct employment of the exponential disc profile.In addition to a more realistic description of the density andpotential of disc galaxies, our semi-analytical framework can accu-rately reproduce the orbital decay of MPs due to DF. For an arbitrarymass distribution, an analytical description of the MP orbital de-cay from first principles is generally unfeasible, since such complexphenomenon depends both on local and non-local effects (see, e.g.,the recent discussion in Tamfal et al. 2020). Any semi-analyticalformulation must therefore introduce a number of approximations.We leverage on the fact that generally, at a first level of approxima-tion, local effects are more relevant in the vast majority of commonastrophysical situations. This allows us to employ the formalism de-rived by Chandrasekhar (1943) based on two-body interactions. Wedescribed the DF for both the spherical profiles and the disc. For thelatter, which features a net rotation, we describe DF as a functionof the relative velocity between the MP and the local surroundingmedium. We compared the evolution in the semi-analytical frame-work against full N-body simulations, and we showed that, whenadditional information about the limited resolution of the latter isincluded in our setup, the agreement between the two is excellent.We conclude by discussing some possible improvements ofthe newly discussed model. First, even considering non-evolving hosts with smooth density profiles, the DF description could be im-proved in multiple ways. In particular, the distribution function ofthe individual components of the host galaxy should be computedself-consistently considering the whole galaxy potential, while inthe current analysis the velocity distributions of the bulge and haloare computed as if they were in isolation. In addition, we follow theapproximation discussed in Chandrasekhar (1943), in which onlystars moving slower than the MP contribute to DF. Such approxi-mation is less and less valid for cored spherical systems (Antonini& Merritt 2012). While such approximation has a limited impactfor Hernquist-like profiles (justifying the agreement we find withnumerical simulations), it could severely affect the orbital decayof MPs in bulges with close to constant density cores (Antonini &Merritt 2012). We also limited our analysis to MPs with fixed masses, while ingeneral a substantial mass evolution is expected in the case of galaxymergers. For instance, during the early stages of the inspiral, thedecaying MBH is expected to be embedded in a sizable fraction ofits original stellar core. The resulting MP is therefore characterisedby both an increased effective mass (determining a more efficientDF) and a larger size. However, as the decay proceeds, tidal effectsexerted by the host galaxy on the MP tend to gradually erode theresidual stellar envelope, until we are left with a “naked” MBH(see e.g. Van Wassenhove et al. 2014). On the other hand, if largereservoirs of gas are available and promptly funnelled toward thesecondary nucleus, the decaying MBH could accrete a considerableamount of it and thus gradually increasing its mass as time passes(Callegari et al. 2011b; Capelo et al. 2015; Capelo & Dotti 2017).Additional improvements can be considered relaxing the as-sumption of smooth/axi-symmetric density distributions. Specifi-cally, when considering the galactic gaseous medium, it has beenfound that in young galaxies gas can be quite turbulent and evenclumpy. This introduces stochastic patterns in the motion of MPs,that depending on the specific gas conditions can significantly alterthe otherwise smooth decay (Fiacconi et al. 2013b; Roškar et al.2015; Tamburello et al. 2017; Souza Lima et al. 2017). Stocasticityis also typical in marginally Toomre unstable galactic discs, whichmay lead to the formation of bars and/or spirals. Such structuresstrongly deviates from axi-symmetry and the torques that they exerton inspiralling MPs can significantly disturb their orbits, in the mostextreme situations even scattering them away (Bortolas et al. 2020).Finally, a more realistic semi-analytical model should also con-sider the possible time evolution of the host galaxy during the orbitalevolution of the MP. If the decay from large separation lasts for aquite long time, of the order of several Gyr, the host galaxy can sig-nificantly evolve, in particular at high redshift, where accretion ofpristine gas from the cosmic filaments can substantially change thetotal mass and possibly affect the main galaxy geometry (see e.g. thediscussion in Rosas-Guevara et al. 2020). Major galaxy mergers canalso strongly perturb the galaxy structure and determining radicalchanges for any MP evolution. Unfortunately, given the very com-plex physics involved in such mergers, a semi-analytical descriptionis challenging and full numerical simulation should be considered.We plan to address the above-mentioned caveats by including We stress that the results presented in Antonini & Merritt (2012) refer tothe specific case of spherical stellar systems in quasi-Keplerian potentials,i.e. it is valid only for the final stages of DF within the influence sphere ofthe host central MBH. While the same distribution function has been usedfor MPs at much larger separations (e.g. Li et al. 2020), a self-consistentimplementation would require the use of the correct distribution function atevery scale. MNRAS , 1–15 (2020)
P evolution in realistic galaxy models them in our semi-analytical framework, to provide a fast and asaccurate as possible description of the orbital decay of satellitesonto evolving disc galaxies. ACKNOWLEDGEMENTS
We thank Jo Bovy for insightful discussion and for his precioushelp in setting up an exponential disc with sech vertical profilewithin galpy. Numerical calculations have been made possiblethrough a CINECA-INFN agreement, providing access to resourceson GALILEO and MARCONI at CINECA. AL acknowledges sup-port from the European Research Council No. 740120 ‘INTER-STELLAR’. EB acknowledges support from the Swiss NationalScience Foundation under the Grant 200020_178949.
DATA AVAILABILITY STATEMENT
The data underlying this article will be shared on reasonable requestto the corresponding author.
REFERENCES
Abramowitz M., Stegun I. A., 1972, Handbook of Mathematical Functions.Dover PublicationsAmaro-Seoane P., et al., 2017, arXiv e-prints, p. arXiv:1702.00786Antonini F., Merritt D., 2012, ApJ, 745, 83Arca-Sedda M., Capuzzo-Dolcetta R., 2014, ApJ, 785, 51Barros D. A., Lépine J. R. D., Dias W. S., 2016, A&A, 593, A108Begelman M. C., Blandford R. D., Rees M. J., 1980, Nature, 287, 307Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition. PrincetonUniveristy PressBonetti M., Perego A., Dotti M., Cescutti G., 2019, MNRAS, 490, 296Bonetti M., Bortolas E., Lupi A., Dotti M., Raimundo S. I., 2020, MNRAS,494, 3053Bortolas E., Capelo P. R., Zana T., Mayer L., Bonetti M., Dotti M., DaviesM. B., Madau P., 2020, MNRAS, 498, 3601Boubert D., Erkal D., Gualandris A., 2020, MNRAS, 497, 2930Bovy J., 2015, ApJS, 216, 29Callegari S., Kazantzidis S., Mayer L., Colpi M., Bellovary J. M., Quinn T.,Wadsley J., 2011a, ApJ, 729, 85Callegari S., Kazantzidis S., Mayer L., Colpi M., Bellovary J. M., Quinn T.,Wadsley J., 2011b, ApJ, 729, 85Capelo P. R., Dotti M., 2017, MNRAS, 465, 2643Capelo P. R., Volonteri M., Dotti M., Bellovary J. M., Mayer L., GovernatoF., 2015, MNRAS, 447, 2123Casertano S., 1983, MNRAS, 203, 735Chandrasekhar I. S., 1941, ApJ, 93, 285Chandrasekhar S., 1943, ApJ, 97, 255Colpi M., Mayer L., Governato F., 1999, ApJ, 525, 720Dehnen W., 1993, MNRAS, 265, 250Dotti M., Colpi M., Haardt F., 2006, MNRAS, 367, 103Dotti M., Colpi M., Haardt F., Mayer L., 2007, MNRAS, 379, 956Fiacconi D., Mayer L., Roškar R., Colpi M., 2013a, ApJ, 777, L14Fiacconi D., Mayer L., Roškar R., Colpi M., 2013b, ApJ, 777, L14Gaia Collaboration et al., 2016, A&A, 595, A1Gaia Collaboration et al., 2018, A&A, 616, A1Gradshteyn I. S., Ryzhik I. M., 2007, Table of integrals, series, and products,seventh edn. Elsevier/Academic Press, AmsterdamGranados A., Torres D., Castañeda L., Henao L., Vanegas S., 2021, New A,82, 101456Hashimoto Y., Funato Y., Makino J., 2003, ApJ, 582, 196Hernquist L., 1990, ApJ, 356, 359Hernquist L., 1993, ApJS, 86, 389 Hopkins P. F., 2015, MNRAS, 450, 53Just A., Peñarrubia J., 2005, A&A, 431, 861Just A., Khan F. M., Berczik P., Ernst A., Spurzem R., 2011, MNRAS, 411,653Kashlinsky A., 1986, ApJ, 306, 374Kuijken K., Gilmore G., 1989, MNRAS, 239, 571Li K., Bogdanović T., Ballantyne D. R., 2020, ApJ, 896, 113Michalski K. A., Mosig J. R., 2016, Journal of Electromagnetic Waves andApplications, 30, 281Miyamoto M., Nagai R., 1975, PASJ, 27, 533Mulder W. A., 1983, A&A, 117, 9Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493Ogata H., 2005, Publications of the Research Institute for MathematicalSciences, 41Ostriker E. C., 1999, ApJ, 513, 252Petts J. A., Gualandris A., Read J. I., 2015, MNRAS, 454, 3778Petts J. A., Read J. I., Gualandris A., 2016, MNRAS, 463, 858Pfister H., Lupi A., Capelo P. R., Volonteri M., Bellovary J. M., Dotti M.,2017, MNRAS, 471, 3646Pfister H., Volonteri M., Dubois Y., Dotti M., Colpi M., 2019, MNRAS,486, 101Plummer H. C., 1911, MNRAS, 71, 460Read J. I., Goerdt T., Moore B., Pontzen A. P., Stadel J., Lake G., 2006,MNRAS, 373, 1451Rosas-Guevara Y., et al., 2020, MNRAS, 491, 2547Roškar R., Fiacconi D., Mayer L., Kazantzidis S., Quinn T. R., Wadsley J.,2015, MNRAS, 449, 494Souza Lima R., Mayer L., Capelo P. R., Bellovary J. M., 2017, ApJ, 838, 13Spitzer Lyman J., 1942, ApJ, 95, 329Tamburello V., Capelo P. R., Mayer L., Bellovary J. M., Wadsley J. W., 2017,MNRAS, 464, 2952Tamfal T., Mayer L., Quinn T. R., Capelo P. R., Kazantzidis S., Babul A.,Potter D., 2020, arXiv e-prints, p. arXiv:2007.13763Tremaine S., Weinberg M. D., 1984, MNRAS, 209, 729Tremmel M., Governato F., Volonteri M., Quinn T. R., 2015, MNRAS, 451,1868Van Wassenhove S., Capelo P. R., Volonteri M., Dotti M., Bellovary J. M.,Mayer L., Governato F., 2014, MNRAS, 439, 474Vasiliev E., 2019, MNRAS, 482, 1525Verbiest J. P. W., et al., 2016, MNRAS, 458, 1267Weinberg M. D., 1986, ApJ, 300, 93White S. D. M., 1983, ApJ, 274, 53Yurin D., Springel V., 2014, MNRAS, 444, 62de Grijs R., Peletier R. F., van der Kruit P. C., 1997, A&A, 327, 966MNRAS000
Abramowitz M., Stegun I. A., 1972, Handbook of Mathematical Functions.Dover PublicationsAmaro-Seoane P., et al., 2017, arXiv e-prints, p. arXiv:1702.00786Antonini F., Merritt D., 2012, ApJ, 745, 83Arca-Sedda M., Capuzzo-Dolcetta R., 2014, ApJ, 785, 51Barros D. A., Lépine J. R. D., Dias W. S., 2016, A&A, 593, A108Begelman M. C., Blandford R. D., Rees M. J., 1980, Nature, 287, 307Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition. PrincetonUniveristy PressBonetti M., Perego A., Dotti M., Cescutti G., 2019, MNRAS, 490, 296Bonetti M., Bortolas E., Lupi A., Dotti M., Raimundo S. I., 2020, MNRAS,494, 3053Bortolas E., Capelo P. R., Zana T., Mayer L., Bonetti M., Dotti M., DaviesM. B., Madau P., 2020, MNRAS, 498, 3601Boubert D., Erkal D., Gualandris A., 2020, MNRAS, 497, 2930Bovy J., 2015, ApJS, 216, 29Callegari S., Kazantzidis S., Mayer L., Colpi M., Bellovary J. M., Quinn T.,Wadsley J., 2011a, ApJ, 729, 85Callegari S., Kazantzidis S., Mayer L., Colpi M., Bellovary J. M., Quinn T.,Wadsley J., 2011b, ApJ, 729, 85Capelo P. R., Dotti M., 2017, MNRAS, 465, 2643Capelo P. R., Volonteri M., Dotti M., Bellovary J. M., Mayer L., GovernatoF., 2015, MNRAS, 447, 2123Casertano S., 1983, MNRAS, 203, 735Chandrasekhar I. S., 1941, ApJ, 93, 285Chandrasekhar S., 1943, ApJ, 97, 255Colpi M., Mayer L., Governato F., 1999, ApJ, 525, 720Dehnen W., 1993, MNRAS, 265, 250Dotti M., Colpi M., Haardt F., 2006, MNRAS, 367, 103Dotti M., Colpi M., Haardt F., Mayer L., 2007, MNRAS, 379, 956Fiacconi D., Mayer L., Roškar R., Colpi M., 2013a, ApJ, 777, L14Fiacconi D., Mayer L., Roškar R., Colpi M., 2013b, ApJ, 777, L14Gaia Collaboration et al., 2016, A&A, 595, A1Gaia Collaboration et al., 2018, A&A, 616, A1Gradshteyn I. S., Ryzhik I. M., 2007, Table of integrals, series, and products,seventh edn. Elsevier/Academic Press, AmsterdamGranados A., Torres D., Castañeda L., Henao L., Vanegas S., 2021, New A,82, 101456Hashimoto Y., Funato Y., Makino J., 2003, ApJ, 582, 196Hernquist L., 1990, ApJ, 356, 359Hernquist L., 1993, ApJS, 86, 389 Hopkins P. F., 2015, MNRAS, 450, 53Just A., Peñarrubia J., 2005, A&A, 431, 861Just A., Khan F. M., Berczik P., Ernst A., Spurzem R., 2011, MNRAS, 411,653Kashlinsky A., 1986, ApJ, 306, 374Kuijken K., Gilmore G., 1989, MNRAS, 239, 571Li K., Bogdanović T., Ballantyne D. R., 2020, ApJ, 896, 113Michalski K. A., Mosig J. R., 2016, Journal of Electromagnetic Waves andApplications, 30, 281Miyamoto M., Nagai R., 1975, PASJ, 27, 533Mulder W. A., 1983, A&A, 117, 9Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493Ogata H., 2005, Publications of the Research Institute for MathematicalSciences, 41Ostriker E. C., 1999, ApJ, 513, 252Petts J. A., Gualandris A., Read J. I., 2015, MNRAS, 454, 3778Petts J. A., Read J. I., Gualandris A., 2016, MNRAS, 463, 858Pfister H., Lupi A., Capelo P. R., Volonteri M., Bellovary J. M., Dotti M.,2017, MNRAS, 471, 3646Pfister H., Volonteri M., Dubois Y., Dotti M., Colpi M., 2019, MNRAS,486, 101Plummer H. C., 1911, MNRAS, 71, 460Read J. I., Goerdt T., Moore B., Pontzen A. P., Stadel J., Lake G., 2006,MNRAS, 373, 1451Rosas-Guevara Y., et al., 2020, MNRAS, 491, 2547Roškar R., Fiacconi D., Mayer L., Kazantzidis S., Quinn T. R., Wadsley J.,2015, MNRAS, 449, 494Souza Lima R., Mayer L., Capelo P. R., Bellovary J. M., 2017, ApJ, 838, 13Spitzer Lyman J., 1942, ApJ, 95, 329Tamburello V., Capelo P. R., Mayer L., Bellovary J. M., Wadsley J. W., 2017,MNRAS, 464, 2952Tamfal T., Mayer L., Quinn T. R., Capelo P. R., Kazantzidis S., Babul A.,Potter D., 2020, arXiv e-prints, p. arXiv:2007.13763Tremaine S., Weinberg M. D., 1984, MNRAS, 209, 729Tremmel M., Governato F., Volonteri M., Quinn T. R., 2015, MNRAS, 451,1868Van Wassenhove S., Capelo P. R., Volonteri M., Dotti M., Bellovary J. M.,Mayer L., Governato F., 2014, MNRAS, 439, 474Vasiliev E., 2019, MNRAS, 482, 1525Verbiest J. P. W., et al., 2016, MNRAS, 458, 1267Weinberg M. D., 1986, ApJ, 300, 93White S. D. M., 1983, ApJ, 274, 53Yurin D., Springel V., 2014, MNRAS, 444, 62de Grijs R., Peletier R. F., van der Kruit P. C., 1997, A&A, 327, 966MNRAS000 , 1–15 (2020) M. Bonetti et al.
APPENDIX A: POTENTIAL AND ACCELERATION FORTHICK EXPONENTIAL DISCS
In this appendix, we detail the derivation leading to the form ofpotential and accelerations employed in this work. Our derivationfollows those presented in Casertano (1983) and Kuijken & Gilmore(1989), with the exception that we express our results in terms of theGauss hypergeometric function 𝐹 ( 𝑎, 𝑏 ; 𝑐 ; 𝑧 ) (see e.g. Abramowitz& Stegun 1972, for full details; see also Appendix B for the em-ployed specific implementation).We consider a thick exponential disc with density profile givenby: 𝜌 ( 𝑅, 𝑧 ) = 𝜌 𝑑, 𝜌 𝑅 ( 𝑅 ) 𝜌 𝑧 ( 𝑧 ) , (A1)where 𝜌 𝑅 ( 𝑅 ) = e − 𝛼𝑅 , (A2) 𝜌 𝑧 ( 𝑧 ) = sech (cid:18) 𝛽𝑧 (cid:19) = cosh − (cid:18) 𝛽𝑧 (cid:19) , (A3)where 𝜌 𝑑, represents the normalisation constant (see equation 1), 𝛼 = / 𝑅 𝑑 is the inverse of the scale radius, while 𝛽 = / 𝑧 𝑑 is twicethe inverse of the scale-height.In order to obtain the gravitational potential 𝜙 generated by theabove mass distribution, we need to solve the Poisson equation, i.e. ∇ 𝜙 ( 𝑅, 𝑧 ) = 𝜋𝐺 𝜌 ( 𝑅, 𝑧 ) . (A4)When axial symmetry holds, equation (A4) can be solved in termsof Hankel’s transform, defined as˜ 𝑓 ( 𝑘 ) = ∫ ∞ d 𝑅 𝑅 𝐽 ( 𝑘 𝑅 ) 𝑓 ( 𝑅 ) , (A5) 𝑓 ( 𝑅 ) = ∫ ∞ d 𝑘 𝑘 𝐽 ( 𝑘 𝑅 ) ˜ 𝑓 ( 𝑘 ) , (A6)where 𝐽 is the 0-th order Bessel function of the first kind, and 𝑓 ( 𝑅 ) is a generic function.By taking the Hankel’s transform of both sides of equa-tion (A4), we obtain the linear non-homogeneous second orderdifferential equation − 𝑘 ˜ 𝜙 ( 𝑘, 𝑧 ) + 𝜕 ˜ 𝜙𝜕𝑧 ( 𝑘, 𝑧 ) = 𝜋𝐺 ˜ 𝜌 ( 𝑘, 𝑧 ) , (A7)that, once solved through standard techniques, yields˜ 𝜙 ( 𝑘, 𝑧 ) = − 𝜋𝐺𝑘 ∫ +∞−∞ d 𝜁 e − 𝑘 | 𝑧 − 𝜁 | ˜ 𝜌 ( 𝑘, 𝜁 ) , (A8)provided that the Hankel transform of the density profile, ˜ 𝜌 ( 𝑘, 𝑧 ) ,vanishes when | 𝑧 | → ±∞ (actually our case). Through equa-tion (A6), the potential reads 𝜙 ( 𝑅, 𝑧 ) = − 𝜋𝐺 ∫ ∞ d 𝑘𝐽 ( 𝑘 𝑅 ) ∫ +∞−∞ d 𝜁 e − 𝑘 | 𝑧 − 𝜁 | ˜ 𝜌 ( 𝑘, 𝜁 ) . (A9)Since the density profile in equation (A1) is factorized in theradial and vertical part, the Hankel transform only affects the former,such that equation (A9) reads 𝜙 ( 𝑅, 𝑧 ) = − 𝜋𝐺 𝜌 𝑑, ∫ ∞ d 𝑘𝐽 ( 𝑘 𝑅 ) ×× ∫ ∞ d 𝑢 𝑢𝐽 ( 𝑢𝑘 ) 𝜌 𝑅 ( 𝑢 ) ∫ +∞−∞ d 𝜁 e − 𝑘 | 𝑧 − 𝜁 | 𝜌 𝑧 ( 𝜁 ) . (A10)The 𝑢 -integral can be evaluated analytically (see e.g. Kuijken& Gilmore 1989) and gives ∫ ∞ d 𝑢 𝑢𝐽 ( 𝑢𝑘 ) e − 𝛼𝑢 = 𝛼 ( 𝛼 + 𝑘 ) / , (A11)while the 𝜁 -integral needs some additional manipulations. In par-ticular, we can first split the integration domain and change the signof the integration variable of the [−∞ , 𝑧 ] part, obtaining 𝐼 𝑧 ( 𝑘 ) = ∫ +∞− 𝑧 d 𝜁 e − 𝑘 ( 𝑧 + 𝜁 ) cosh (cid:16) 𝛽𝜁 (cid:17) + ∫ +∞ 𝑧 d 𝜁 e − 𝑘 ( 𝜁 − 𝑧 ) cosh (cid:16) 𝛽𝜁 (cid:17) , (A12)that after the variable change 𝜁 ± 𝑧 = ( / 𝛽 ) 𝑦 we reduce to 𝐼 𝑧 ( 𝑘 ) = 𝛽 ∫ +∞ d 𝑦 e − 𝑘𝑦 / 𝛽 cosh (cid:16) 𝑦 − 𝛽𝑧 (cid:17) + ∫ +∞ d 𝑦 e − 𝑘𝑦 / 𝛽 cosh (cid:16) 𝑦 + 𝛽𝑧 (cid:17) . (A13)Finally, we use the symbolic software Mathematica to express it ina closed form in terms of the Gauss hypergeometric function 𝐹 ,i.e. 𝐼 𝑧 ( 𝑘 ) = 𝛽 (cid:40) − 𝑘𝑘 + 𝛽 (cid:20) e − 𝑧𝛽 𝐹 (cid:18) , + 𝑘𝛽 ; 2 + 𝑘𝛽 ; − e − 𝑧𝛽 (cid:19) + e 𝑧𝛽 𝐹 (cid:18) , + 𝑘𝛽 ; 2 + 𝑘𝛽 ; − e 𝑧𝛽 (cid:19)(cid:21)(cid:41) . (A14)Going back to equation (A10), the final form of the gravita-tional potential can be written as 𝜙 ( 𝑅, 𝑧 ) = − 𝜋𝐺𝛼𝜌 𝑑, ∫ ∞ d 𝑘𝐽 ( 𝑘 𝑅 ) 𝐼 𝑧 ( 𝑘 )( 𝛼 + 𝑘 ) / . (A15)To obtain the radial and vertical accelerations necessary tointegrate the equations of motion, we take the (negative) gradientof equation (A15), obtaining 𝑎 𝑅 = − 𝜕𝜙𝜕𝑅 = − 𝜋𝐺𝛼𝜌 𝑑, ∫ ∞ d 𝑘 𝑘𝐽 ( 𝑘 𝑅 ) 𝐼 𝑧 ( 𝑘 )( 𝛼 + 𝑘 ) / , (A16) 𝑎 𝑧 = − 𝜕𝜙𝜕𝑧 = − 𝜋𝐺𝛼𝜌 𝑑, ∫ ∞ d 𝑘𝐽 ( 𝑘 𝑅 ) − 𝜕 𝑧 𝐼 𝑧 ( 𝑘 )( 𝛼 + 𝑘 ) / , (A17)where 𝜕 𝑧 𝐼 𝑧 ( 𝑘 ) is given by MNRAS , 1–15 (2020)
P evolution in realistic galaxy models 𝜕 𝑧 𝐼 𝑧 ( 𝑘 ) = 𝑘 sgn ( 𝑧 ) 𝛽 ( 𝑘 + 𝛽 ) (cid:34) − e − 𝑧𝛽 𝐹 (cid:18) , + 𝑘𝛽 ; 2 + 𝑘𝛽 ; − e − 𝑧𝛽 (cid:19) + e 𝑧𝛽 𝐹 (cid:18) , + 𝑘𝛽 ; 2 + 𝑘𝛽 ; − e 𝑧𝛽 (cid:19) − 𝑘 + 𝛽𝑘 tanh (cid:18) 𝑧𝛽 (cid:19)(cid:35) . (A18)Equations (A15–A17) express the gravitational potential, theradial and vertical accelerations in form of 1-dimensional integrals.Unfortunately, no further simplifications are possible and no closedform expressions can be found, therefore those integrals have to beevaluated numerically. Despite the apparent simplicity of dealingwith 1-dimensional integrals, the presence of the Bessel functions 𝐽 and 𝐽 can make the task quite challenging given their highly os-cillatory behaviour, requiring therefore specific integration schemesto obtain sensible outcomes. Nevertheless, once this task is prop-erly accomplished, it allows us to integrate the orbits with arbitraryinitial conditions. Specifically, to numerically integrate such kindof oscillatory functions we employ a special variant of the DoubleExponential (DE) rule. In the special case of equatorial motion (i.e. 𝑧 = 𝐼 ( 𝑘 ) = ∫ +∞ d 𝜁 e − 𝑘𝜁 cosh (cid:16) 𝛽𝜁 (cid:17) , (A19)that with some manipulations evaluates to (see e.g. Gradshteyn &Ryzhik 2007, p. 383) 𝐼 ( 𝑘 ) = 𝛽 (cid:40) 𝑘𝛽 (cid:20) 𝜓 (cid:18) 𝑘 𝛽 + (cid:19) − 𝜓 (cid:18) 𝑘 𝛽 (cid:19)(cid:21) − (cid:41) , (A20)where the function 𝜓 is the digamma function (see e.g. Abramowitz& Stegun 1972).Finally, since the functions 𝐼 𝑧 ( 𝑘 ) , 𝐼 ( 𝑘 ) , and 𝜕 𝑧 𝐼 𝑧 ( 𝑘 ) containaddition and subtractions of special functions, roundoff error canrepresent a serious issue for the numerical evaluation of the integralswhen 𝑘 → +∞ . To circumvent this possible problem, we can replacethe expression of those functions with their asymptotic expansionsfor 𝑘 → +∞ , i.e. 𝐼 𝑧 ( 𝑘 ) = − 𝑧𝛽 ( + e − 𝑧𝛽 ) 𝑘 + O( / 𝑘 ) , (A21) 𝐼 ( 𝑘 ) = 𝑘 − 𝛽 𝑘 + O( / 𝑘 ) , (A22) 𝜕 𝑧 𝐼 𝑧 ( 𝑘 ) = 𝛽 e − 𝑧𝛽 ( − e − 𝑧𝛽 )( + e − 𝑧𝛽 ) ( 𝑘 + 𝛽 ) + O( / 𝑘 ) . (A23) We refer the interested reader to Ogata (2005) and Michalski & Mosig(2016) for a complete description of the method.
APPENDIX B: COMPUTATION OF THE GAUSSHYPERGEOMETRIC FUNCTION
The Gauss hypergeometric function 𝐹 ( 𝑎, 𝑏 ; 𝑐 ; 𝑧 ) is a special func-tion defined for | 𝑧 | < 𝐹 ( 𝑎, 𝑏 ; 𝑐 ; 𝑧 ) = Γ ( 𝑐 ) Γ ( 𝑎 ) Γ ( 𝑏 ) ∞ ∑︁ 𝑛 = Γ ( 𝑎 + 𝑛 ) Γ ( 𝑏 + 𝑛 ) Γ ( 𝑐 + 𝑛 ) 𝑧 𝑛 𝑛 ! , (B1)and by analytic continuation elsewhere on the complex plane. Inthe above definition Γ (·) represents the gamma function. The Gausshypergeometric function is very general and contains as specialcases several other mathematical functions. A detailed discussionabout the properties and relations involving the hypergeometricfunction is definitely beyond the scope of this work, for which weare only interested in the specific case when the hypergeometricfunction assumes the form 𝐹 ( , + 𝑦 ; 2 + 𝑦 ; 𝑧 ) , (B2)with 𝑧 ∈ R always negative and 𝑦 ∈ R and positive.Several software implementations of the 𝐹 function exist,but, at least to our knowledge, none of them is able to providea sensible evaluation for all possible combinations of ( 𝑦, 𝑧 ) , forwhich reliable results are usually provided only when | 𝑧 | < We therefore implemented our own algorithm to evaluate 𝐹 for | 𝑧 | >
1, relying instead on the GNU GSL mathematical library for | 𝑧 | ≤ 𝐹 , we first note that twolimiting cases exist, i.e. when 𝑦 → 𝑦 → +∞ . In these specificsituations, 𝐹 assumes the following special values 𝐹 ( ,
1; 2; 𝑧 ) = − ln ( − 𝑧 ) 𝑧 , (B3) 𝐹 ( , 𝑏 ; 𝑏 ; 𝑧 ) = ( − 𝑧 ) − , (B4)allowing us to express 𝐹 as one of the limiting value plus smallcorrections.Specialising equation (B1) with parameters in equation (B2),we get 𝐹 ( , + 𝑦 ; 2 + 𝑦 ; 𝑧 ) = ( + 𝑦 ) ∞ ∑︁ 𝑛 = 𝑧 𝑛 + 𝑦 + 𝑛 . (B5)In the case of 𝑦 →
0, we can expand the function as 𝐹 ( , + 𝑦 ; 2 + 𝑦 ; 𝑧 ) = ( + 𝑦 ) ∞ ∑︁ 𝑘 = (− ) 𝑘 𝑦 𝑘 ∞ ∑︁ 𝑛 = 𝑧 𝑛 ( + 𝑛 ) 𝑘 + , (B6)that, changing the summation index to 𝑚 = 𝑛 +
1, yields 𝐹 ( , + 𝑦 ; 2 + 𝑦 ; 𝑧 ) = + 𝑦𝑧 ∞ ∑︁ 𝑘 = (− ) 𝑘 𝑦 𝑘 ∞ ∑︁ 𝑚 = 𝑧 𝑚 𝑚 𝑘 + . (B7)where the series over 𝑚 defines the polylogarithm function(Abramowitz & Stegun 1972) of order ( 𝑘 + ) , denoted 𝐿𝑖 𝑘 + ( 𝑧 ) .Starting from the 1st order of the function 𝐿𝑖 ( 𝑧 ) = − ln ( − 𝑧 )/ 𝑧 , all A notable exception is the mpmath
Pyhton package, which unfortunatelycannot be efficiently used in our C++ implementation.MNRAS000
Pyhton package, which unfortunatelycannot be efficiently used in our C++ implementation.MNRAS000 , 1–15 (2020) M. Bonetti et al. lower (i.e. the 0-th and negative) orders can be expressed in closedform exploiting the recurrence relation 𝑧 𝜕𝐿𝑖 𝑠 𝜕𝑧 ( 𝑧 ) = 𝐿𝑖 𝑠 − ( 𝑧 ) . (B8)On the other hand, above the 2nd positive order (included) a closedform cannot be found. Exploiting the polylogarithm function, equa-tion (B7) can be written as 𝐹 ( , + 𝑦 ; 2 + 𝑦 ; 𝑧 ) = − + 𝑦𝑧 ln ( − 𝑧 )+ + 𝑦𝑧 ∞ ∑︁ 𝑘 = (− ) 𝑘 𝑦 𝑘 𝐿𝑖 𝑘 + ( 𝑧 ) . (B9)Depending on the value of 𝑦 we can then truncate the series to matcha desired accuracy. Finally, by setting 𝑦 = 𝑦 → +∞ , we follow the verysame procedure outlined above, this time expanding the term con-taining 𝑦 as 1 + 𝑦 + 𝑦 + 𝑛 = + ∞ ∑︁ 𝑙 = (− ) 𝑙 𝑛 ( + 𝑛 ) 𝑙 − 𝑦 𝑙 , (B10)that, inserted into equation (B5), gives 𝐹 ( , + 𝑦 ; 2 + 𝑦 ; 𝑧 ) = ∞ ∑︁ 𝑛 = 𝑧 𝑛 + ∞ ∑︁ 𝑛 = 𝑧 𝑛 ∞ ∑︁ 𝑙 = (− ) 𝑙 𝑛 ( + 𝑛 ) 𝑙 − 𝑦 𝑙 . (B11)The first term involving only powers of 𝑧 is the series definitionfor 1 /( − 𝑧 ) . The second term, instead, can be further expandedexploiting the binomial theorem, leading to 𝐹 ( , + 𝑦 ; 2 + 𝑦 ; 𝑧 ) == ( − 𝑧 ) − + ∞ ∑︁ 𝑙 = (− ) 𝑙 𝑦 𝑙 𝑙 − ∑︁ 𝑘 = (cid:18) 𝑙 − 𝑘 (cid:19) ∞ ∑︁ 𝑛 = 𝑛 𝑘 + 𝑧 𝑛 , = ( − 𝑧 ) − + ∞ ∑︁ 𝑙 = (− ) 𝑙 𝑦 𝑙 𝑙 − ∑︁ 𝑘 = (cid:18) 𝑙 − 𝑘 (cid:19) 𝐿𝑖 −( 𝑘 + ) ( 𝑧 ) , (B12)where, in the third line, we recognised the definition of the poly-logarithm function. Since the polylogarithm function appears withnegative orders, each of them can be easily evaluated through therecurrence relation in equation (B8).Finally, in all cases for which the limiting forms for 𝑦 → 𝑦 → +∞ are not suitable, we evaluate the hypergeometricfunction using its series definition in equation (B1). Formally, thisprovides a sensible result only when | 𝑧 | <
1, but exploiting the lineartransformation property of 𝐹 (Abramowitz & Stegun 1972) andperforming the transformation 𝑧 → / 𝑧 , an alternative expressionof 𝐹 as a function of 1 / 𝑧 can be obtained. Once reduced andspecialised to the case of interest, it reads 𝐹 ( , + 𝑦 ; 2 + 𝑦 ; 𝑧 ) = + 𝑦𝑧 ∞ ∑︁ 𝑛 = ( 𝑛 − 𝑦 ) 𝑧 𝑛 − 𝜋 ( + 𝑦 ) sin ( 𝜋𝑦 ) 𝑧 + 𝑦 . (B13) APPENDIX C: USEFUL EXPRESSIONS FOR DF WITHFAST MOVING STARS
Here, we provide the analytical form of the integral over the variable 𝑉 that appears in equation (8), i.e. 𝐽 ( 𝑣 𝑝 , 𝑣 ★ , 𝑝 max ) = ∫ 𝑣 𝑝 + 𝑣 ★ | 𝑣 𝑝 − 𝑣 ★ | d 𝑉 (cid:32) + 𝑣 𝑝 − 𝑣 ★ 𝑉 (cid:33) ln (cid:32) + 𝑝 𝑉 𝐺 𝑚 𝑝 (cid:33) . (C1)where 𝑣 𝑝 and 𝑣 ★ are the MP and star velocities, whereas 𝑝 max is the maximum impact parameter. To compute the integral, weneed to consider the two different cases in which 𝑣 𝑝 > 𝑣 ★ and 𝑣 𝑝 < 𝑣 ★ . For convenience, we express the result in terms of theGauss hypergeometric function 𝐹 (•) = 𝐹 ( , / , / , •) . When 𝑣 𝑝 > 𝑣 ★ , the integration yields18 𝑣 ★ ∫ 𝑣 𝑝 + 𝑣 ★ 𝑣 𝑝 − 𝑣 ★ d 𝑉 (cid:32) + 𝑣 𝑝 − 𝑣 ★ 𝑉 (cid:33) ln (cid:32) + 𝑝 𝑉 𝐺 𝑚 𝑝 (cid:33) = 𝑣 ★ (cid:34) − (cid:18) 𝐺𝑚 𝑝 𝑝 max (cid:19) 𝑣 𝑝 + 𝑣 ★ ( 𝑣 𝑝 − 𝑣 ★ ) 𝐹 ( 𝛼 ) + (cid:18) 𝑝 max 𝐺𝑚 𝑝 (cid:19) ( 𝑣 𝑝 − 𝑣 ★ ) 𝐹 ( 𝛽 )+ (cid:18) 𝐺𝑚 𝑝 𝑝 max (cid:19) 𝑣 𝑝 − 𝑣 ★ ( 𝑣 𝑝 + 𝑣 ★ ) 𝐹 ( 𝛾 ) − (cid:18) 𝑝 max 𝐺𝑚 𝑝 (cid:19) ( 𝑣 𝑝 + 𝑣 ★ ) 𝐹 ( 𝛿 )+ 𝑣 ★ (cid:32) + ln (cid:16) + 𝑝 𝐺 𝑚 𝑝 ( 𝑣 𝑝 − 𝑣 ★ ) + 𝑝 𝐺 𝑚 𝑝 ( 𝑣 + 𝑣 𝑣 ★ + 𝑣 ★ ) (cid:17)(cid:33)(cid:35) , (C2)while, when 𝑣 𝑝 < 𝑣 ★ , we obtain18 𝑣 ★ ∫ 𝑣 𝑝 + 𝑣 ★ − 𝑣 𝑝 + 𝑣 ★ d 𝑉 (cid:32) + 𝑣 𝑝 − 𝑣 ★ 𝑉 (cid:33) ln (cid:32) + 𝑝 𝑉 𝐺 𝑚 𝑝 (cid:33) = 𝑣 ★ (cid:34) − 𝑣 𝑝 + (cid:18) 𝐺𝑚 𝑝 𝑝 max (cid:19) 𝑣 𝑝 + 𝑣 ★ ( 𝑣 𝑝 − 𝑣 ★ ) 𝐹 ( 𝛼 ) − (cid:18) 𝑝 max 𝐺𝑚 𝑝 (cid:19) ( 𝑣 𝑝 − 𝑣 ★ ) 𝐹 ( 𝛽 )+ (cid:18) 𝐺𝑚 𝑝 𝑝 max (cid:19) 𝑣 𝑝 − 𝑣 ★ ( 𝑣 𝑝 + 𝑣 ★ ) 𝐹 ( 𝛾 ) − (cid:18) 𝑝 max 𝐺𝑚 𝑝 (cid:19) ( 𝑣 𝑝 + 𝑣 ★ ) 𝐹 ( 𝛿 )+ 𝑣 ★ ln (cid:169)(cid:173)(cid:173)(cid:171) + 𝑝 𝐺 𝑚 𝑝 ( 𝑣 𝑝 + 𝑣 ★ ) + 𝑝 𝐺 𝑚 𝑝 ( 𝑣 𝑝 − 𝑣 ★ ) (cid:170)(cid:174)(cid:174)(cid:172)(cid:35) , (C3) Alternatively, the interested reader can follow the steps outlined in Chan-drasekhar (1941) to obtain an expression in terms of elementary functions.MNRAS , 1–15 (2020)
P evolution in realistic galaxy models where in both expression we have defined 𝛼 = − 𝐺 𝑚 𝑝 𝑝 ( 𝑣 𝑝 − 𝑣 ★ ) , (C4) 𝛽 = − 𝑝 ( 𝑣 𝑝 − 𝑣 ★ ) 𝐺 𝑚 𝑝 , (C5) 𝛾 = − 𝐺 𝑚 𝑝 𝑝 ( 𝑣 𝑝 + 𝑣 ★ ) , (C6) 𝛿 = − 𝑝 ( 𝑣 𝑝 + 𝑣 ★ ) 𝐺 𝑚 𝑝 . (C7)We note that, if we assume the logarithmic term in equa-tion (C3) to be independent of 𝑉 , then the resulting integral vanishes,and the equation simplifies to the well known result in which onlyparticles moving slower than the MP contribute to DF. This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000