Dynamical evolution of massive perturbers in realistic multicomponent galaxy models I: implementation and validation

Bonetti, Matteo; Bortolas, Elisa; Lupi, Alessandro; Dotti, Massimo

doi:10.1093/mnras/stab222

ABSTRACT

Galaxies are self-gravitating structures composed by several components encompassing spherical, axial, and triaxial symmetry. Although real systems feature heterogeneous components whose properties are intimately connected, semi-analytical approaches often exploit the linearity of the Poisson’s equation to represent the potential and mass distribution of a multicomponent galaxy as the sum of the individual components. In this work, we expand the semi-analytical framework developed in Bonetti et al. (2020) by including both a detailed implementation of the gravitational potential of exponential disc (modelled with a sech² and an exponential vertical profile) and an accurate prescription for the dynamical friction experienced by massive perturbers (MP) in composite galaxy models featuring rotating disc structures. Such improvements allow us to evolve arbitrary orbits either within or outside the galactic disc plane. We validate the results obtained by our numerical model against public semi-analytical codes as well as full N-body simulations, finding that our model is in excellent agreement to the codes it is compared with. The ability to reproduce the relevant physical processes responsible for the evolution of MP orbits and its computational efficiency make our framework perfectly suited for large parameter-space exploration studies.

gravitation, methods: numerical, stars: kinematics and dynamics, galaxies: nuclei, galaxies: structure

1 INTRODUCTION

In the last decades, more and more interest has grown towards the modelling of orbital motions within the potential of a galaxy (e.g. Bovy 2015; Boubert, Erkal & Gualandris 2020; Granados et al. 2021). This has to be attributed to the advent of missions like Gaia (Gaia Collaboration et al. 2016, 2018), which is tracking the orbits of stars and stellar clusters sailing though the Galactic potential, and will release an 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 the rates of such mergers, one needs to model the inspiral phase of an intruder MBH within the potential of the host galaxy, down to the scale at which it will form a binary with another black hole that may already sit in the centre (Begelman, Blandford & Rees 1980).

A detailed modelling of galactic orbits can only be achieved with a sufficiently good description of the associated galactic potential wells; in addition, if one wants to accurately follow the orbit of a massive perturber (MP, i.e. an object whose mass is significantly larger than the characteristic stellar mass), it is crucial to further account for the effect of dynamical friction (DF; Chandrasekhar 1943; Binney & Tremaine 2008).

At the present time, state-of-the art numerical simulations with sub-grid recipes for unresolved physical processes promise to be the best tool to properly model the orbital evolution of stars and MPs within galaxies (e.g. Tremmel et al. 2015; Pfister et al. 2017, 2019). However, they cannot be adopted when one needs to perform a throughout exploration of the parameter space as they are typically very computationally expensive. For this, different tools, such as semi-analytical approaches, may be preferred owing to their much lower computational cost that comes at the expense of some precision.

In semi-analytical models, the motion of a test particle subject to the galactic potential can be followed more and more accurately as the description of the potential associated to each galactic component (the dark matter (DM) halo, disc, bulge, etc.) gets closer to the one of real galaxies. While the spherically symmetric galactic components can be approximated reasonably well by well-known potential–density pairs featuring a simple analytical form (e.g. Plummer 1911; Dehnen 1993; Navarro, Frenk & White 1997), the same cannot be said for less idealized galactic structures such as discs. In fact, their light distribution is often well described by an exponential disc density profile that unfortunately does not admit a close analytical form for the related potential, meaning that its computation has to be performed numerically and the associated computational cost may be relatively large. For this, less physically motivated disc potentials featuring a simpler, close analytical form are often adopted in the literature, even if their approximation of the potential of a physical galactic disc may be relatively poor.

Even neglecting the aforementioned issues, a proper description of the galactic potential in semi-analytical models is not enough when the target particle is a MP, as DF may significantly affect the evolution of its orbit. DF arises as a response of the background to the perturbing mass, and it typically results in a deceleration of the latter that gradually sinks towards the centre of its host system. Its standard semi-analytical description grounds its basis on the seminal work of Chandrasekhar (1943), who first analytically modelled the motion of an MP in an infinite and homogeneous stellar background.¹ In spite of its simplicity, Chandrasekhar’s approach has been proven to work remarkably well in many scenarios, sometimes far beyond the expected regime of applicability (e.g. White 1983), and despite missing the physics of resonant interactions between the MP and the host system (Tremaine & Weinberg 1984; Weinberg 1986). Owing to its success, the semi-analytical treatment for DF has seen a substantial development throughout the last decades, as attested by the vast piece of recent literature focusing on improving Chandrasekhar’s prescriptions and their agreement with N-body integrations (e.g. Hashimoto, Funato & Makino 2003; Just & Peñarrubia 2005; Read et al. 2006; Just et al. 2011; Arca-Sedda & Capuzzo-Dolcetta 2014; Petts, Gualandris & Read 2015; Petts, Read & Gualandris 2016). However, the work in this direction has been limited to very idealized, spherically symmetric and isotropic systems that can only describe a narrow range of galaxy environments; in fact, astrophysical systems feature much more complexity than this: an extension of the semi-analytical DF treatment to composite systems, possibly featuring disky components with rotational support, would thus greatly enhance the realm of applicability of semi-analytical DF.

N-body experiments already addressed the effect of DF in a composite, rotationally supported environment; they highlighted that DF acting in rotationally supported systems induces the circularization of a MP initially on a prograde, eccentric orbit, and reverses the angular momentum of counter-rotating MPs, then again promoting circularization (Dotti, Colpi & Haardt 2006; Callegari et al. 2011a; Fiacconi et al. 2013a); this effect appears to be independent of the nature of the background (Dotti et al. 2007). However, such past studies never attempted a quantitative modelling of the phenomenon. Only very recently, Bonetti et al. (2020) successfully modelled for the first time the DF-induced ‘drag towards circular co-rotation’ in a semi-analytical fashion in a multicomponent galaxy that featured a DM halo, a stellar disc, and a bulge; their prescription showed remarkable agreement with N-body simulations down to the circularization phase. However, their study was limited to on-plane orbits, and their approach can be adopted only before circularization occurs.

Here, we expand the work presented in Bonetti et al. (2020) by improving the DF treatment in rotating discs with the development of a more physical DF prescription. Such prescription features no arbitrary tunable parameters and through the employment of a more realistic (but still isotropic) disc velocity distribution function, we extend its validity also to circular orbits. Further to DF treatment, we also provide a completely revisited implementation for the acceleration evaluation due to an exponential disc mass distribution, now allowing arbitrary orbits and not only equatorial motion.

On a longer time-span, this work is intended to be the first of a series aiming at providing a semi-analytical computational framework that encodes the most realistic as possible dynamics of MPs in multicomponent galaxies, but at the same time still featuring a high computational efficiency, able to guarantee vast and systematic parameter space explorations in terms of diverse galactic profiles and MP initial properties/trajectories.

The paper is organized as follows: in Section 2, we describe the improvements brought to the framework developed in Bonetti et al. (2020); in Section 3, we validate our updated set-up against existing semi-analytical codes and full N-body simulations featuring the same initial conditions. Finally, in Section 4, we draw our conclusions.

2 METHOD

As in Bonetti et al. (2020), we consider a multicomponent galaxy model comprising a Navarro–Frenk–White (NFW; Navarro et al. 1997) profile to model the DM halo, an Hernquist (Hernquist 1990) profile to describe a compact stellar bulge and a thick exponential disc profile to characterize the galactic disc. On top of that, we model the effect of DF on the motion of MPs through the employment of dissipative forces.

In this section, we detail the procedure we followed to (i) evaluate the potential and accelerations generated by a thick exponential disc in the general case, (ii) discuss possible improvement for the DF force in spherically symmetric systems, and (iii) extend the validity of our prescription for DF in rotating discs to nearly circular orbits.

2.1 Exponential disc conservative dynamics

We consider a finite-thickness exponential disc whose mass density is described by

$$\begin{eqnarray*} \rho _d(R,z) = \frac{M_d}{4\pi R_\mathrm{ d}^2 z_d} {\rm e}^{-R/R_d} {\rm sech}^2\left(\frac{z}{z_d}\right), \end{eqnarray*}$$

(1)

where M_d and R_d are the total mass and scale radius of the disc, while z_d is a parameter shaping the profile along the z direction. Such profile is our choice for the modelling of general galactic discs, since it can well reproduce the light profiles of real galaxies and it is physically motivated, being the vertical distribution obtained for an isothermal self-gravitating disc (Spitzer 1942). For these reasons, it is usually employed to generate the initial conditions of N-body simulations (see e.g. Yurin & Springel 2014), allowing for a validation of the results of our semi-analytical prescriptions against N-body runs. Nevertheless, the double exponential (DE) disc, i.e. a disc with vertical density profile decaying as |${\rm e}^{-|z|/z_d}$|⁠, is widely used and well reproduces the vertical profiles of edge-on disc galaxies (e.g. de Grijs, Peletier & van der Kruit 1997). In order to have a proper comparison, we also implemented this second profile in our framework.²

In general, to obtain the gravitational potential from a given density distribution, one has to solve the Poisson’s equation. Practically, when deviating from spherical symmetry, a complete analytical expression can be found only in a handful of cases and unfortunately the thick exponential case does not belong to those. Therefore, to obtain the gravitational potential (and consequently the accelerations) of such mass distribution at arbitrary (R, z) pairs, we necessarily need to recur to numerical methods. Unfortunately, the general solution of the Poisson’s equation that can be obtained via Green’s function results in multidimensional integrals, implying that obtaining the solution can be pretty expensive in computational terms. Thus, if possible, it is advisable to reduce the expressions as much as possible before recurring to numerical integration.

It turns out that, when considering an axisymmetric density distribution, the Poisson’s equation can be solved in terms of a Hankel’s transform (see e.g. Kuijken & Gilmore 1989). This approach allows us to write the gravitational potential as a one-dimensional integral (see Appendix A for a complete derivation), given by

$$\begin{eqnarray*} \phi (R,z) = - \frac{G M_d}{2 R_\mathrm{ d}^3 z_d} \int _0^\infty ~{\rm d}k J_0(k R) \frac{I_z(k)}{\left(R_\mathrm{ d}^{-2}+k^2\right)^{3/2}}, \end{eqnarray*}$$

(2)

where J₀ is the Bessel function of the first kind, and I_z(k) is a function depending on both k and z through the Gauss hypergeometric function ₂F₁ (see Appendixes A and B).

The radial and vertical accelerations can be readily obtained by taking the (negative) derivative of equation (2) with respect to R and z and then performing the integration over the k variable

$$\begin{eqnarray*} a_R = -\frac{\partial \phi }{\partial R} = - \frac{G M_d}{2 R_\mathrm{ d}^3 z_d} \int _0^\infty ~{\rm d}k k J_1(k R) \frac{I_z(k)}{\left(R_\mathrm{ d}^{-2}+k^2\right)^{3/2}}, \end{eqnarray*}$$

(3)

$$\begin{eqnarray*} a_z = -\frac{\partial \phi }{\partial z} = - \frac{G M_d}{2 R_\mathrm{ d}^3 z_d} \int _0^\infty ~{\rm d}k J_0(k R) \frac{-\partial _z I_z(k)}{\left(R_\mathrm{ d}^{-2}+k^2\right)^{3/2}}. \end{eqnarray*}$$

(4)

Note that the integration of Bessel functions could be problematic given their highly oscillatory behaviour. Nevertheless, if appropriate quadrature techniques are employed, the integrals can be computed quite efficiently using a few hundreds of points, still guaranteeing a good level of accuracy (see Appendix A).

2.2 DF in spherical profiles

DF, though not dissipative in nature, in the context of semi-analytical calculations is usually modelled as a drag force acting opposite to the velocity of a MP. Despite the effort to improve the original derivation by Chandrasekhar (1943; see e.g. Mulder 1983; Tremaine & Weinberg 1984; Weinberg 1986; Colpi, Mayer & Governato 1999, and reference therein) the rather simple original expression, i.e.

$$\begin{eqnarray*} \mathbf {a}_{\rm df} = -2\pi G^2 \ln (1+\Lambda ^2) m_p \rho \left({\rm erf}(X) - \frac{2 X{\rm e}^{-X^2}}{\sqrt{\pi }}\right) \frac{\mathbf {v}_{p}}{|\mathbf {v}_{p}|^3}, \end{eqnarray*}$$

(5)

is still widely used, mostly because of its simplicity. In the above formula, m_p and |$\mathbf {v}_p$| are the MP mass and velocity, ρ is the background density, Λ = p_max/p_min is the ratio between the maximum and minimum impact parameter, while |$X=v_p/(\sqrt{2}\sigma)$| is the ratio of the MP velocity over the velocity dispersion. Though systematically employed, equation (5) is subjected to a number of assumptions that are not always justified in realistic galactic contexts, i.e.:

The motion of scattering particles with large impact parameters is rectilinear;
The background density is homogeneous and isotropic;
The velocity distribution function is isotropic and Maxwellian with constant dispersion σ;
The maximum impact parameter is chosen equal to a maximum scale-length of the system, while p_min is kept fixed for all possible encounters and equal to a characteristic scale length (or, in the case of point-like perturbers, the radius below which a 90° deflection occurs). Specifically, this last assumption implies that only encounters with velocity lower than the MP velocity contribute to the deceleration.

Despite the above assumptions are clearly not fully compatible with realistic galaxies, the introduction of some a posteriori modifications (mostly affecting points 2 and 3 above) to equation (5) allows to reproduce the results from N-body simulations with reasonable accuracy. In particular, when limiting the DF treatment to spherically symmetric profiles (which model to a fairly good approximation elliptical galaxies or galaxy bulges), the adoption of the position-dependent mass density, velocity dispersion, and minimum and maximum impact parameters substantially improves the agreement between semi-analytical models and full N-body simulations (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-dependent velocity dispersion σ(r) and mass density ρ(r), and the following expressions for the maximum and minimum impact parameters:

$$\begin{eqnarray*} p_{\rm max} &=& r/\gamma , \nonumber \\ p_{\rm min} &=& \max \left(\frac{G m_p}{v_p^2+\sigma (r)^2}, D_\mathrm{ p}\right), \nonumber \\ \gamma &=& - \frac{~{\rm d}\ln \rho }{~{\rm d}\ln r}, \end{eqnarray*}$$

(6)

with γ denoting the logarithmic slope of the density profile, r the radial coordinate, and D_p the physical radius of the MP (which can be set to zero for the case of an MBH).

Even with the above modifications, the acceleration computed with the functional form of equation (5) considers only the contribution of scattering particles moving slower than the MP. Although this does not represent a serious issue for the majority of stellar systems, the deceleration described by equation (5) could severely underestimate the real DF when the number of background objects moving faster than the MP is large. Specifically, Antonini & Merritt (2012) investigated the DF experienced by a MP in galactic cores dominated by a central MBH and they found that, when the slope of the density profile is γ ≲ 1, the contribution of the so-called ‘fast-moving stars’ can become comparable to or even larger than that of slow moving ones.

In order to take into account such fast population, one needs to drop some assumptions and consider the general Chandrasekhar formula (see equations 25 and 26 of Chandrasekhar 1943) that, despite being more general, usually does not allow for an analytical form, i.e.³

$$\begin{eqnarray*} \mathbf {a}_{\rm df} &=& -4\pi G^2 m_p \rho (r) \frac{\mathbf {v}_p}{|\mathbf {v}_p|^3} \nonumber \\ &&\times \int _0^{v_{\rm esc}} \frac{J(v_p,v_\star ,p_{\rm max})}{8\pi v_\star } 4\pi v_\star ^2 f(v_\star) ~{\rm d}v_\star , \end{eqnarray*}$$

(7)

$$\begin{eqnarray*} J(v_p,v_\star ,p_{\rm max}) &=& \int _{|v_p-v_\star |}^{v_p+v_\star } ~{\rm d}V \left(1+\frac{v_p^2-v_\star ^2}{V^2}\right) \nonumber\\ &&\times \,\ln \left(1+\frac{p_{\rm max}^2 V^4}{G^2 m_p^2}\right). \end{eqnarray*}$$

(8)

Equation (7) inevitably implies a higher computational burden due to the unavoidable numerical computation of the integral over the distribution function. Given that, in the vast majority of astrophysical situations, the inclusion of fast moving stars represents only a minor correction (but see Section 4 for some caveats of our choice), in this work, we adopt the simplified expression considering only slow-moving stars, together with the modifications introduced by different authors that improve its range of applicability.

2.3 DF in rotating discs

Bonetti et al. (2020) derived a prescription to describe the DF acceleration acting on a MP determined by a rotating disc,

$$\begin{eqnarray*} \mathbf {a}_{\rm df,disc} = -A_{\rm disc} 2\pi G^2 \ln (1+\Lambda ^2) m_p \rho _d(R,0) \frac{\mathbf {v}_p-\mathbf {v}_c(R)}{|\mathbf {v}_p-\mathbf {v}_c(R)|^3}, \end{eqnarray*}$$

(9)

where |$\mathbf {v}_p,m_p$| are the velocity and mass of the MP, |$\mathbf {v}_c(R)$| is the circular velocity at radius R determined by the total gravitational potential, while Λ = p_max,d/p_min,d, which enters an effective Coulomb logarithm, is chosen as the ratio of the maximum impact parameter, taken equal to the z_d parameter of the disc, and a minimum one given by

$$\begin{eqnarray*} p_{\rm min, d} &=& \frac{G m_p}{v_{\rm rel}^2 + 0.01 v_c^2},\nonumber \\ v_{\rm rel} &=& |\mathbf {v}_p-\mathbf {v}_c|. \end{eqnarray*}$$

(10)

Finally, the factor A_disc is a tunable constant set by comparing the semi-analytical orbital evolution of the MP with that obtained in an N-body simulation. Bonetti et al. (2020) found that this factor is of the order of unity, enforcing that their modelling catches the main physics of DF in discs.

The main limitation of this prescription lies on the dependence of the DF acceleration on the relative velocity with respect to circular velocity. When the MP velocity approaches the circular one, equation (9) diverges, thus the orbital decay cannot be followed further.

To overcome this limitation, here we consider a less idealized distribution function for the stars in the disc. Specifically, we replace the delta function with an isotropic Gaussian centred around the local rotational velocity. This choice allows us to employ the very same functional form of the standard DF Chandrasekhar formulation, but with the velocity dependence on the perturber velocity in the galactic frame replaced by the relative velocity with respect to the local medium (see e.g. Kashlinsky 1986, for a similar derivation). In addition, with this new prescription, we also discard the tunable A_disc parameter previously introduced in Bonetti et al. (2020). In particular, we obtain

$$\begin{eqnarray*} \mathbf {a}_{\rm df,disc} &=& -2\pi G^2 \ln (1+\Lambda ^2) m_p \rho _d(R,z) \nonumber \\ &&\times \left({\rm erf}(X) - \frac{2 X{\rm e}^{-X^2}}{\sqrt{\pi }}\right) \frac{\mathbf {v}_{\rm rel}}{|\mathbf {v}_{\rm rel}|^3}, \end{eqnarray*}$$

(11)

where |$\mathbf {v}_{\rm rel} = \mathbf {v}_p - \mathbf {v}_{\rm rot}(R)$|⁠, while |$X = |\mathbf {v}_{\rm rel}|/(\sqrt{2}\sigma _R)$|⁠, with σ_R denoting the radial velocity dispersion. For the Coulomb logarithm, we replace the expression of p_min,d appearing in equation (10) by

$$\begin{eqnarray*} p_{\rm min, d} = \frac{G m_p}{v_{\rm rel}^2 + \sigma _R^2}, \end{eqnarray*}$$

(12)

which takes into account both the relative velocity with surrounding medium, which can be different from the circular velocity, and an intrinsic dispersion that sets the minimum effective distance characterizing the encounters between the MP and the objects in the disc.

The presence of a velocity dispersion σ_R implies that the disc is not fully rotationally supported. We can characterize the rotational velocity by following Hernquist (1993),

$$\begin{eqnarray*} |\mathbf {v}_{\rm rot}|^2 = v_\mathrm{ c}^2(R) + \sigma _R^2(R) \left(1-\frac{\kappa ^2}{4\Omega ^2} - 2\frac{R}{R_\mathrm{ d}}\right), \end{eqnarray*}$$

(13)

in which v_c denotes the circular velocity at R, while κ and Ω (both depending on R) are the epicyclic and angular frequency, respectively. In evaluating such quantities for multicomponent galaxy models, we consider the total potential rather than that generated by the disc only. As in Hernquist (1993), we also assume that the velocity dispersion changes as

$$\begin{eqnarray*} \sigma _R^2(R) \propto {\rm e}^{-R/R_\mathrm{ d}}, \end{eqnarray*}$$

(14)

where the normalization is evaluated assuming equality between σ_R and the critical radial velocity dispersion (augmented by a factor Q)⁴ at a specific radius, here 2R_d, i.e.

$$\begin{eqnarray*} \sigma _{R, \rm crit} = Q \times \frac{3.36 G \ \Sigma (2 R_\mathrm{ d})}{\kappa (2 R_\mathrm{ d})}. \end{eqnarray*}$$

(15)

We can note that since in equation (13) the term in parenthesis on the right-hand side is usually negative, then the rotational velocity results (often only slightly) lower than the local circular velocity. However, as already pointed out in Hernquist (1993), at small radii, the approximations yielding equation (13) can break down and provide nonphysical v_rot (e.g. an imaginary velocity). To avoid unnecessary complications, we decided to fix the rotational velocity equal to a fraction of the local circular velocity (specifically 0.95v_c(R)) anytime equation (13) gives imaginary values. We expect our choice not to dramatically impact the evolution, since at small radii the dynamics is likely dominated by the stellar bulge (and the associated DF), while the contribution of the disc should be subdominant. We anyway try to address this issue by also directly extracting the rotational profile from an N-body realization, as we will further discuss in Section 3.

3 RESULTS

In this section, we compare our semi-analytical setup with other similar implementations for the disc potential and the DF treatment, and against full N-body simulations. We start by considering the conservative dynamics in exponential discs, highlighting how a different vertical disc profile can affect the orbits. We then analyse our DF prescriptions in multicomponent galaxy models, and we compare them against the results from N-body simulations featuring the same physical system and MP orbital initial conditions. We employ the public code gizmo (Hopkins 2015) to run the N-body simulations, using the same N-body setup of Bonetti et al. (2020), whose parameters are reported in Table 1 for completeness.⁵

Table 1.

Open in new tab

Structural parameters of the employed galaxy model. For each component, N_part and ε correspond to the number of particles and the Plummer-equivalent gravitational softening employed in the N-body runs.

	Halo	Bulge	Disc
Mass	\|$1.1 \times 10^{12}~\rm M_{\odot }$\|	\|$2.2\times 10^{9}~\rm M_{\odot }$\|	\|$4.4\times 10^{10}~\rm M_{\odot }$\|
Scale radius	37 (21) kpc	0.96 kpc	4.25 kpc
Vertical scale	–	–	0.85 kpc
Profile	Hernquist (NFW)	Hernquist	Exponential
N_part	10⁶	5 × 10⁵	2 × 10⁶
ε	40 pc	10 pc	10 pc

	Halo	Bulge	Disc
Mass	\|$1.1 \times 10^{12}~\rm M_{\odot }$\|	\|$2.2\times 10^{9}~\rm M_{\odot }$\|	\|$4.4\times 10^{10}~\rm M_{\odot }$\|
Scale radius	37 (21) kpc	0.96 kpc	4.25 kpc
Vertical scale	–	–	0.85 kpc
Profile	Hernquist (NFW)	Hernquist	Exponential
N_part	10⁶	5 × 10⁵	2 × 10⁶
ε	40 pc	10 pc	10 pc

Table 1.

Open in new tab

Structural parameters of the employed galaxy model. For each component, N_part and ε correspond to the number of particles and the Plummer-equivalent gravitational softening employed in the N-body runs.

	Halo	Bulge	Disc
Mass	\|$1.1 \times 10^{12}~\rm M_{\odot }$\|	\|$2.2\times 10^{9}~\rm M_{\odot }$\|	\|$4.4\times 10^{10}~\rm M_{\odot }$\|
Scale radius	37 (21) kpc	0.96 kpc	4.25 kpc
Vertical scale	–	–	0.85 kpc
Profile	Hernquist (NFW)	Hernquist	Exponential
N_part	10⁶	5 × 10⁵	2 × 10⁶
ε	40 pc	10 pc	10 pc

	Halo	Bulge	Disc
Mass	\|$1.1 \times 10^{12}~\rm M_{\odot }$\|	\|$2.2\times 10^{9}~\rm M_{\odot }$\|	\|$4.4\times 10^{10}~\rm M_{\odot }$\|
Scale radius	37 (21) kpc	0.96 kpc	4.25 kpc
Vertical scale	–	–	0.85 kpc
Profile	Hernquist (NFW)	Hernquist	Exponential
N_part	10⁶	5 × 10⁵	2 × 10⁶
ε	40 pc	10 pc	10 pc

3.1 Orbits in isolated exponential discs

Here, we investigate the conservative dynamics in single-component galaxies modelled with exponential discs. In our framework, we consider both the exponential disc with z-profile decaying as sech²(z/z_d; see Section 2.1) as well as the so-called DE disc, i.e. a profile also decaying as an exponential along z.

3.1.1 Consistency checks

In Fig. 1, we report the relative error in the radial (left column) and vertical (right column) accelerations of an exponential disc with a sech²(z/z_d) profile, computed by our semi-analytical framework through a modified version of the DE quadrature technique (see e.g. Michalski & Mosig 2016). In the upper panels, we compare our calculation to the numerical evaluation of equations (3) and (4) obtained via the trapezoidal rule with ∼ 2 × 10⁷ points and assuming a k_max ≈ 10⁵. We can appreciate that the relative errors are always very small (10⁻¹²–10⁻¹¹) over several decades in both R and z. Only well above 10R/R_d the relative error starts to grow significantly, but still without exceeding 10⁻⁶. This trend is associated with the behaviour of the Bessel function in the integrand function, which oscillates significantly for large R, mildly reducing the accuracy of the numerical integration.

Figure 1.

Relative error for the radial (left column) and vertical (right column) accelerations obtained with our code and considering an exponential disc with vertical sech² profile. Upper panels: relative errors are evaluated against a numerical integration of equations (3) and (4) with 2 × 10⁷ points using trapezoidal rule in the range k ∈ [0, 10⁵]. Lower panels: relative errors are obtained comparing our acceleration with those computed by the software agama employing an exponential disc with identical parameters.

Open in new tab Download slide

In the lower panels, instead, we evaluate the relative errors against the accelerations obtained from the software agama (Vasiliev 2019), in which the disc density profile is expanded in multipoles⁶ and the potential (and the corresponding force) is computed via 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⁻⁶ for a vast portion of the (R − z) plane, while larger deviations arise at the borders of the plane. Given that the computation in agama relies on the multipole expansion and makes systematic use of spline interpolation, the quantities that agama evaluates are necessarily approximated and subjected to some numerical noise. Still, errors are quite small in the relevant portion of the (R − z) plane and reach 10⁻² only for R and z that are much larger or much smaller than R_d and z_d. This does not represent a real problem as in realistic galaxy models other galactic components are expected to dominate the dynamics at large/small distances, for instance dark matter haloes (at large distances) or stellar bulges (close to the center), hence we can consider our implementation quite robust.

3.1.2 Orbit integration in exponential discs

We next compare the orbital evolution that we obtain⁷ with those of the software galpy (Bovy 2015). Specifically, we analyse the trajectories obtained in an exponential disc with mass |$M_d = 4.4 \times 10^{10} ~\rm M_{\odot }$|⁠, scale radius R_d = 4.25 kpc, and z_d = 0.85 kpc and with a vertical profile given by either sech²(z/z_d) or |${\rm e}^{-|z|/z_d}$| (see also Table 1). In the top panel of Fig. 2, we compare the orbit for the sech² case (blue solid for our implementation versus orange dashed for galpy)⁸ by showing the trajectory projections in the x − y, x − z, R − z planes. The orbit is evolved for approximately 2 Gyr starting from an initial radial separation of R = 5 kpc and an initial velocity along the y and z directions equal to 0.5 times the local circular velocity. From the top panel, we can infer that the orbit in our implementation closely matches that of galpy and to better quantify the agreement, in the bottom panel of the figure, we evaluate the relative error on the spherical radius between the two runs, finding that is at the level of 10⁻² at most. This small difference is probably due to the different approaches followed to compute the disc potential: specifically through the employment of the Hankel’s transform in our case, which reduces the problem to the numerical computation of an integral (see equation 2), or via a multipole expansion approximation in the galpy framework. We perform the same comparison by also considering the DE disc profile. As for the sech², we found a good agreement therefore enforcing the robustness of our framework. Another important comparison with galpy concerns the computational cost, which for a semi-analytical framework has to necessarily be modest in order to fully exploit its power. For each run, we find a total run-time of the order of few tens of seconds (on a single core of a standard laptop), similar to that of galpy (but only once an adequate accuracy, similar to that in our framework, is also adopted in galpy). Moreover, to further speed-up the orbital integration, we have also implemented the possibility of storing in advance the accelerations computed on a (R − z) adaptive grid and then employ a bi-cubic interpolation to obtain a_R and a_z on arbitrary points.

$Top: orbit projections in the x − y, x − z, and R − z planes when considering a sech2 exponential disc with mass $M_d = 4.4 \times 10^{10} ~\rm M_{\odot }$, scale radius Rd = 4.25 kpc, and zd = 0.85 kpc. We show the comparison between the trajectories obtained within our framework (solid blue lines) and those obtained with the software galpy (dashed orange lines). The initial velocity is set equal to 0.5vc in both y and z directions, with vc the circular velocity at a cylindrical radius R = 5 kpc. Bottom: relative error on the spherical radius between our evolution and that obtained with galpy. Note how the relative error remains the most at the level of 10−2 after 2 Gyr of evolution.$

Figure 2.

Top: orbit projections in the x − y, x − z, and R − z planes when considering a sech² exponential disc with mass |$M_d = 4.4 \times 10^{10} ~\rm M_{\odot }$|⁠, scale radius R_d = 4.25 kpc, and z_d = 0.85 kpc. We show the comparison between the trajectories obtained within our framework (solid blue lines) and those obtained with the software galpy (dashed orange lines). The initial velocity is set equal to 0.5v_c in both y and z directions, with v_c the circular velocity at a cylindrical radius R = 5 kpc. Bottom: relative error on the spherical radius between our evolution and that obtained with galpy. Note how the relative error remains the most at the level of 10⁻² after 2 Gyr of evolution.

Open in new tab Download slide

Finally, we compare the orbital shape of a point mass when subject to the potential of a DE disc and a sech² one. We consider the same disc and initial conditions as described above, but here we set the velocity in the z direction equal to 0.1 times the local circular velocity. The results are shown in Fig. 3. In both the left-hand and right-hand panels, we show the three projections over the x − y, x − z, and y − z planes as sub-panels (from top left to bottom right). We can appreciate that the orbits resemble each other, at least qualitatively, in the sense that they cover more or less the same volume of the disc. Still, from a quantitative point of view, e.g. checking the position reached at a specific time, the orbits are well distinct and evolve on slightly different timescales. Hence, a precise reconstruction of orbits in a realistic galaxy requires the employment of the density/potential pair that better describes such structure that, given the variety of profiles observed, would require an ad hoc model for each individual case. Such exercise can be easily preformed using our procedure, which provides both DE and sech² exponential disc implementations, but is beyond the scope of our current validation study.

$Comparison between orbits in a DE disc (solid blue line) and in an exponential disc with a z-profile described by $\rm sech^2$ (dashed orange lines). In the top (bottom) panel, the in-plane initial velocity is set equal to 0.5 (1.5) the local circular velocity, while vertical velocity is set to 0.1 this value. Each sub-panel shows the orbit projections in the x − y, x − z, and R − z planes.$

Figure 3.

Comparison between orbits in a DE disc (solid blue line) and in an exponential disc with a z-profile described by |$\rm sech^2$| (dashed orange lines). In the top (bottom) panel, the in-plane initial velocity is set equal to 0.5 (1.5) the local circular velocity, while vertical velocity is set to 0.1 this value. Each sub-panel shows the orbit projections in the x − y, x − z, and R − z planes.

Open in new tab Download slide

3.2 DF in multicomponents models

Now, we consider the MP evolution under the DF combined effect of all galactic components. Before comparing our semi-analytical evolution with that of an N-body simulation, we first assess the relative strength of DF from different galactic components, specifically halo, bulge, and disc. In the upper panel of Fig. 4, we show a generic MP radial decay achievable within our framework, while in the bottom panel, we show the corresponding instantaneous relative contribution of DF generated by each of the considered galactic component (see labels) normalized to the total acceleration suffered by the MP. From the bottom panel, we can clearly appreciate that, for a galaxy with structural parameters chosen as in Tab. 1, our framework predicts an early evolution where the dominant DF is actually generated by the disc component, and we therefore have to expect specific phenomenology connected to the net rotation of the disc structure, for instance the ‘drag towards circular co-rotation’. Although sub-dominant contributions by the halo and bulge components are also present, only the halo has a non-negligible effect, being the bulge too compact to properly affect the MP orbit quite outside the scale radius. However, as the evolution proceeds, the above picture reverses. As soon as the MP reach a radial separation of ∼2 kpc (comparable with the bulge length-scale), the DF contribution from disc weakens and DF from spherical components takes over, especially the contribution from the stellar bulge that dominates the MP orbital decay at small scales. From the above findings, we expect, in a generic disc galaxy modelled with at least three components, the evolution of a MP to be initially driven by the interaction with the disc, therefore making its modelization quite relevant.

Figure 4.

Upper panel: MP radial evolution. Lower panel: relative contributions of the DF accelerations from the three galactic components (each acceleration modulus is normalized to the total acceleration, i.e. both dissipative and conservative).

Open in new tab Download slide

We now proceed to validate our semi-analytical implementation against full N-body simulations. In Fig. 5, we show, from top to bottom, the time evolution of radial separation, the z component of the angular momentum (normalized to its initial value), and the eccentricity of a MP orbiting in the galactic mid-plane (i.e. co-planar with the galactic disc). ⁹ We here compare the evolution obtained with our semi-analytical setup, as described in Section 2 (solid lines), and with N-body simulations (dashed lines). From the figure, we can infer that our semi-analytical setup can qualitatively catch the relevant phenomenology of the evolution displayed by the N-body simulations, such as the trends shown by the angular momentum and eccentricity. For instance, in the prograde case, we clearly observe that DF is responsible for a quite rapid circularization of the MP orbit well before the MP enters the bulge-dominated region (≲1 kpc), while in the initially retrograde case, first DF greatly increases the orbital eccentricity until the angular momentum reverses, then when the orbit becomes prograde it again promotes the circularization. Despite we can capture the behaviour of the N-body runs, still the timescales are not exactly reproduced. Indeed, although the evolution of a co-planar prograde orbit (left-hand panel) seems to be well recovered by our framework, the same does not hold when looking at an initially co-planar retrograde trajectory (right-hand panel), in which the semi-analytical decay looks faster. Such differences, however, can be interpreted by reminding that the N-body simulations also have some limitations and a fair comparison with our semi-analytical framework has to necessarily take into account the limited spatial and mass resolution of N-body runs. Specifically, we need to consider that in N-body codes all interactions below the softening length get considerably weakened and, as a consequence, the resulting DF force acting on a MP is smaller with respect to an ideal and arbitrarily resolved physical system (see e.g. Tremmel et al. 2015; Pfister et al. 2017).

Figure 5.

Time evolution of the radial separation, the normalized 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-hand panel: co-planar prograde orbit. Right-hand panel: co-planar initially retrograde orbit.

Open in new tab Download slide

To account for the limited N-body resolution, we can act on the parameters that enter the expression of the DF force. According to equation (12), we can observe that the minimum impact parameter p_min depends on the MP velocity and specifically for the disc DF it depends on the relative velocity between the MP and the surrounding medium. A small p_min translates into larger DF accelerations through the term ln (1 + Λ²). For the retrograde case, v_rel can be quite large, therefore p_min can decrease significantly, becoming smaller than the softening length of the N-body simulation. To assess the influence of the softening in the N-body framework, we thus force the minimum impact parameter p_min to be larger than 2.8 times the Plummer-equivalent softening length of the corresponding particle type, i.e. halo, bulge, or disc. The evolution resulting from softened semi-analytical setup is shown in Fig. 6 with blue lines. Both the prograde and retrograde orbits (left-hand and right-hand panels, respectively) now show a slower decay, and in particular the retrograde case appears to be much more similar to the N-body evolution. The prograde inspiral however seems now too slow.

Figure 6.

Same as Fig. 5, 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 p_min for each component (disc, halo, and 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 velocity directly extracted from the N-body simulations.

Open in new tab Download slide

Another important effect that needs to be considered is the possible difference between the rotational velocity profile employed in our semi-analytic approach, based on the approximated equation (13; Hernquist 1993), and that of the simulation. As shown in Fig. 7, a direct comparison of the two velocity profiles reveal that, as we move to smaller radii, the difference in the rotational velocity profile between the one used in the semi-analytical approach and that in the N-body simulation gets enhanced. If the rotational velocity is significantly different from the local circular one when the orbit becomes nearly circular (i.e. above 400 Myr in the prograde case), then the disc DF force does not vanish (see equations 11 and 13), resulting in a slightly faster decay. In order to account for this effect, we implemented the rotational velocity profile extracted from the simulation and fitted with a rational function of the radius (see the dashed black line in Fig. 7)

$$\begin{eqnarray*} |\mathbf {v}_{\rm rot}| = \frac{a (R^2 + b R)}{R^2 + c R + d} + K \quad {\rm [km\, s^{-1}]} \end{eqnarray*}$$

(16)

with a = 1.25 × 10², b = 1.24 × 10², c = 6.03 × 10¹, d = 1.91 × 10², and K = −4.41. This modification is enough to improve the agreement between the two approaches, as shown by the green lines in Fig. 6. Moreover, the apparent good agreement between the N-body run and the semi-analytical approach in the left-hand panel of Fig. 5 actually results from a fortunate coincidence related to the different weighting of the DF acceleration generated by each galactic components, specifically a slightly stronger (because not softened) DF from the spherical components and a weaker DF from the disc, due to a velocity profile closer to the circular one.

Figure 7.

Circular velocity (thick blue line), rotational velocity as predicted by equation (13; thin orange line) and extracted from N-body simulation (red crossed line) as a function of radius. At small R, the large oscillations of N-body v_rot are due to the limited particle coverage in the central regions. To avoid spurious numerical noise, we employed an analytical fit of the N-body v_rot (black dashed line).

Open in new tab Download slide

Although the orbits we have considered so far were co-planar with the disc, our semi-analytical framework allows us to explore orbits with arbitrary inclinations (see Section 2.1 for details). Here, we check the orbital decay of MPs starting from off-plane configurations, as reported in Fig. 8, ι = 45° (upper left-hand panels), ι = 135° (upper right-hand panels), and ι = 90° (bottom panels); see additional off-plane cases shown in Appendix D. For each case, instead of the eccentricity, we show here the evolution of the inclination. For the ι = 90° case, we also report the total angular momentum evolution, being L_z ≈ 0. As in Fig. 6, for each inclination, we compare the evolution derived from N-body simulations (orange dashed lines) with that computed in our semi-analytical framework, again taking care of including the softening of the corresponding N-body run (solid blue lines) and the N-body rotational profile (solid green lines). From the figure, we can clearly see that the semi-analytical evolution matches remarkably well the N-body one, correctly reproducing both the observed phenomenology and the decay timescale. Minor differences are of course present, but these are unavoidable given the discrete nature and the limited mass resolution of the N-body approach. Similarly to the planar counterparts, also in these inclined runs a better agreement is met when both the softening and the rotational curve of the N-body are considered.

Figure 8.

Same as Fig. 6, but considering three different initial inclinations, 45° (upper left), 135° (upper right), and 90° (bottom). The bottom panel of each case shows the inclination evolution, while for the ι = 90° run the total angular momentum is shown instead of its z component.

Open in new tab Download slide

We stress again that the very good agreement that we observe in both planar and non-planar cases is achieved only when some input from the N-body runs is provided. Indeed, a fair comparison has to necessarily take into account the limited resolution of N-body simulations. Still, for realistic galaxies, we expect that a faster decay, similar to the one we find when considering our standard framework (see Fig. 5), is actually more physical and should provide a better timescale estimate, valid at least down to radial separations not much smaller than the bulge scale radius. Indeed, well inside this scale radius, we expect the correction from fast moving star to start being relevant, hence a more general DF formulation should be used. Despite this, our framework has proved to excellently capture the physics driving MP motion in realistic multicomponent galaxy models with rotationally supported structures. Moreover, given its quite modest computational cost compared to full N-body simulations, our framework represents a particularly attractive solution for large parameters space explorations in terms of galaxy structures as well as MP masses and orbital configurations.

4 DISCUSSION AND CONCLUSIONS

In this paper, we developed a semi-analytical framework to integrate the motion of particles in realistic galactic potentials composed of multiple components. Among several implemented potential–density pairs (spherical, axisymmetric, and triaxial), we included the cases of exponential discs with vertical density decaying either as |${\rm e}^{-|z|/z_d}$| or sech²(z/z_d), for which we developed an efficient setup to numerically compute its potential and, consequently, the acceleration field. Though more complex to tackle, these exponential profiles have to be preferred to other analytical profiles that can only partially reproduce the exponential decay of the surface brightness of observed disc galaxies. For example, the Miyamoto–Nagai (MN; Miyamoto & Nagai 1975) disc is one of the most common analytical profiles employed to model the potential of disc galaxies. Still, this profile can significantly differ than the one of realistic galaxies, especially at large radii where it falls off as a power law rather than exponentially (see discussion in chapter 2 of Binney & Tremaine 2008). To enhance the density decay, a combination of three MN discs is usually employed, although this requires one of the discs to have a negative mass, often leading to the appearance of non-physical negative density regions in the (R − z) plane (see e.g. Barros, Lépine & Dias 2016; Bonetti et al. 2019). This strongly supports our choice about the direct employment of the exponential disc profile.

In addition to a more realistic description of the density and potential of disc galaxies, our semi-analytical framework can accurately reproduce the orbital decay of MPs due to DF. For an arbitrary mass distribution, an analytical description of the MP orbital decay from first principles is generally unfeasible, since such complex phenomenon depends both on local and non-local effects (see e.g. the recent discussion in Tamfal et al. 2020). Any semi-analytical formulation must therefore introduce a number of approximations. We leverage on the fact that generally, at a first level of approximation, local effects are more relevant in the vast majority of common astrophysical situations. This allows us to employ the formalism derived by Chandrasekhar (1943) based on two-body interactions. We described the DF for both the spherical profiles and the disc. For the latter, which features a net rotation, we describe DF as a function of the relative velocity between the MP and the local surrounding medium. We compared the evolution in the semi-analytical framework against full N-body simulations, and we showed that, when additional information about the limited resolution of the latter is included in our setup, the agreement between the two is excellent.

We conclude by discussing some possible short and mid-term improvements of the newly discussed model. First, even considering non-evolving hosts with smooth density profiles, the DF description could be improved in multiple ways. For instance, when computing the DF expression for the disc, we assumed an isotropic Gaussian as distribution function. A natural extension is to consider an anisotropic version of such distribution function encoding different velocity dispersion along different direction, such as σ_R = σ_ϕ > σ_z (see e.g. Binney 1977). Despite this choice inevitably requires the numerical computation of some integrals, the task could likely be optimized in order to maintain good computational efficiency. A further (indirect) improvement for disc DF consists in the development of more realistic expressions for the rotational pattern |$\mathbf {v}_{\rm rot}$| and σ_R profile, such that a well-behaved form is maintained also in the central regions.

A more challenging aspect concerns how the distribution functions of the individual galactic components are evaluated. In principle, such functions should be computed self-consistently considering the whole galaxy potential, while in the current analysis the velocity distributions of the bulge and halo are computed as if they were in isolation. In addition, we follow the approximation discussed in Chandrasekhar (1943), in which only stars moving slower than the MP contribute to DF. Such approximation is less and less valid for cored spherical systems (Antonini & Merritt 2012). While such approximation has a limited impact for Hernquist-like profiles (justifying the agreement we find with numerical simulations), it could severely affect the orbital decay of MPs in bulges with close to constant density cores (Antonini & Merritt 2012).¹⁰

We also limited our analysis to MPs with fixed masses, while in general a substantial mass evolution is expected in the case of galaxy mergers. For instance, during the early stages of the inspiral, the decaying MBH is expected to be embedded in a sizable fraction of its original stellar core. The resulting MP is therefore characterized by both an increased effective mass (determining a more efficient DF) and a larger size. However, as the decay proceeds, tidal effects exerted by the host galaxy on the MP tend to gradually erode the residual stellar envelope, until we are left with a ‘naked’ MBH (see e.g. Van Wassenhove et al. 2014). On the other hand, if large reservoirs of gas are available and promptly funnelled toward the secondary nucleus, the decaying MBH could accrete a considerable amount 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 assumption of smooth/axisymmetric density distributions. Specifically, when considering the galactic gaseous medium, it has been found that in young galaxies, gas can be quite turbulent and even clumpy. This introduces stochastic patterns in the motion of MPs that, depending on the specific gas conditions, can significantly alter the otherwise smooth decay (Fiacconi et al. 2013b; Roškar et al. 2015; Souza Lima et al. 2017; Tamburello et al. 2017). Stocasticity is also typical in marginally Toomre unstable galactic discs that may lead to the formation of bars and/or spirals. Such structures strongly deviate from axi-symmetry and the torques that they exert on inspiralling MPs can significantly disturb their orbits, in the most extreme situations even scattering them away (Bortolas et al. 2020).

Finally, a more realistic semi-analytical model should also consider the possible time evolution of the host galaxy during the orbital evolution of the MP. If the decay from large separation lasts for a quite long time, of the order of several Gyr, the host galaxy can significantly evolve, in particular at high redshift, where accretion of pristine gas from the cosmic filaments can substantially change the total mass and possibly affect the main galaxy geometry (see e.g. the discussion in Rosas-Guevara et al. 2020). Major galaxy mergers can also strongly perturb the galaxy structure and determining radical changes for any MP evolution. Unfortunately, given the very complex physics involved in such mergers, a semi-analytical description is challenging and full numerical simulation should be considered.

We plan to address the above-mentioned caveats by including them in our semi-analytical framework, to provide a fast and as accurate as possible description of the orbital decay of satellites on to evolving disc galaxies.

ACKNOWLEDGEMENTS

We thank Jo Bovy for insightful discussion and for his precious help in setting up an exponential disc with sech² vertical profile within galpy. Numerical calculations have been made possible through a CINECA-INFN agreement, providing access to resources on GALILEO and MARCONI at CINECA. MD, MB, and AL acknowledge funding from MIUR under the grant PRIN 2017-MB8AEZ. AL acknowledges support from the European Research Council No. 740120 ‘INTERSTELLAR’. EB acknowledges support from the Swiss National Science Foundation under the Grant 200020_178949 and from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program ERC-2018-COG under grant agreement No. 818691 (B Massive).

DATA AVAILABILITY STATEMENT

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

1

Under similar assumptions, Ostriker (1999) describes the analogous phenomenon within a gaseous medium.

2

The DE disc also does not admit any closed form for the gravitational potential, but requires numerical integration as given in equation (2), where the sole difference is in the function I(k) given by equation A9 of Kuijken & Gilmore (1989).

3

The integral over the relative velocity J(v_p, v_⋆, p_max) actually has an analytical form. It is usually quoted in an approximate form that is well behaved everywhere except when v_p and v_⋆ are close (see e.g. equation 7 of Antonini & Merritt 2012). For completeness, we report the full expression in Appendix C.

4

Here, we assume Q = 1.5.

5

In the framework of N-body simulations of isolated galaxies, the Hernquist profile is generally used to model the dark matter halo instead of the well known NFW profile. This choice is barely related to the finite mass of the Hernquist profile that instead diverges for the NFW profile. We note that well within their scale radii, the NFW and Hernquist profiles have the very same shape and therefore when focusing on the dynamics inside a galaxy, the two profiles are practically the same (Springel, Di Matteo & Hernquist 2005).

6

For the comparison in Fig. 1, we adopt within agama a multipole expansion ranging between 10⁻⁶ and 500 R_d; we used a logarithmically spaced grid in R with 1000 grid points and maximum order of the expansion l_max = 100, while the properties of the z grid are kept as default in agama (see Vasiliev 2019, for more details).

7

See Bonetti et al. (2016) for details about the orbit integration method.

8

The sech² exponential disc within the galpy framework is accessible through the SCF basis-function-expansion (see galpy documentation).

9

As in Bonetti et al. (2020), the eccentricity is found by assuming instantaneous conservation of energy and angular momentum and solving for the two radial roots (i.e. peri- and apocentre) such that the radial velocity vanishes. The orbital inclination is instead computed by finding the angle between the angular momentum and the z axis.

10

We stress that the results presented in Antonini & Merritt (2012) refer to the 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 of the host central MBH. While the same distribution function has been used for MPs at much larger separations (e.g. Li, Bogdanović & Ballantyne 2020), a self-consistent implementation would require the use of the correct distribution function at every scale.

11

We refer the interested reader to Ogata (2005) and Michalski & Mosig (2016) for a complete description of the method.

12

A notable exception is the mpmath python package that unfortunately cannot be efficiently used in our c++ implementation.

13

Alternatively, the interested reader can express ₂F₁ in terms of elementary functions employing known mathematical relations. Moreover, following the steps outlined in Chandrasekhar (1941), a more complicated version of the J integral, without neglecting some non-dominant terms, can be derived to obtain a more general, though cumbersome, expression.

REFERENCES

Abramowitz

M.

,

Stegun

I. A.

,

1972

,

Handbook of Mathematical Functions

.

Dover Press

,

New York

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Amaro-Seoane

P.

et al. ,

2017

,

preprint (arXiv:1702.00786)

Antonini

F.

,

Merritt

D.

,

2012

,

ApJ

,

745

,

83

10.1088/0004-637X/745/1/83

Crossref

Arca-Sedda

M.

,

Capuzzo-Dolcetta

R.

,

2014

,

ApJ

,

785

,

51

10.1088/0004-637X/785/1/51

Crossref

Barros

D. A.

,

Lépine

J. R. D.

,

Dias

W. S.

,

2016

,

A&A

,

593

,

A108

10.1051/0004-6361/201527535

Crossref

Begelman

M. C.

,

Blandford

R. D.

,

Rees

M. J.

,

1980

,

Nature

,

287

,

307

10.1038/287307a0

Crossref

Binney

J.

,

1977

,

MNRAS

,

181

,

735

10.1093/mnras/181.4.735

Crossref

Binney

J.

,

Tremaine

S.

,

2008

,

Galactic Dynamics

, 2nd edn.

Princeton Univ. Press

,

Princeton, NJ

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Bonetti

M.

,

Haardt

F.

,

Sesana

A.

,

Barausse

E.

,

2016

,

MNRAS

,

461

,

4419

10.1093/mnras/stw1590

Crossref

Bonetti

M.

,

Perego

A.

,

Dotti

M.

,

Cescutti

G.

,

2019

,

MNRAS

,

490

,

296

10.1093/mnras/stz2554

Crossref

Bonetti

M.

,

Bortolas

E.

,

Lupi

A.

,

Dotti

M.

,

Raimundo

S. I.

,

2020

,

MNRAS

,

494

,

3053

10.1093/mnras/staa964

Crossref

Bortolas

E.

,

Capelo

P. R.

,

Zana

T.

,

Mayer

L.

,

Bonetti

M.

,

Dotti

M.

,

Davies

M. B.

,

Madau

P.

,

2020

,

MNRAS

,

498

,

3601

10.1093/mnras/staa2628

Crossref

Boubert

D.

,

Erkal

D.

,

Gualandris

A.

,

2020

,

MNRAS

,

497

,

2930

10.1093/mnras/staa2211

Crossref

Bovy

J.

,

2015

,

ApJS

,

216

,

29

10.1088/0067-0049/216/2/29

Crossref

Callegari

S.

,

Kazantzidis

S.

,

Mayer

L.

,

Colpi

M.

,

Bellovary

J. M.

,

Quinn

T.

,

Wadsley

J.

,

2011a

,

ApJ

,

729

,

85

10.1088/0004-637X/729/2/85

Crossref

Callegari

S.

,

Kazantzidis

S.

,

Mayer

L.

,

Colpi

M.

,

Bellovary

J. M.

,

Quinn

T.

,

Wadsley

J.

,

2011b

,

ApJ

,

729

,

85

10.1088/0004-637X/729/2/85

Crossref

Capelo

P. R.

,

Dotti

M.

,

2017

,

MNRAS

,

465

,

2643

10.1093/mnras/stw2872

Crossref

Capelo

P. R.

,

Volonteri

M.

,

Dotti

M.

,

Bellovary

J. M.

,

Mayer

L.

,

Governato

F.

,

2015

,

MNRAS

,

447

,

2123

10.1093/mnras/stu2500

Crossref

Casertano

S.

,

1983

,

MNRAS

,

203

,

735

10.1093/mnras/203.3.735

Crossref

Chandrasekhar

I. S.

,

1941

,

ApJ

,

93

,

285

10.1086/144265

Crossref

Chandrasekhar

S.

,

1943

,

ApJ

,

97

,

255

10.1086/144517

Crossref

Colpi

M.

,

Mayer

L.

,

Governato

F.

,

1999

,

ApJ

,

525

,

720

10.1086/307952

Crossref

de Grijs

R.

,

Peletier

R. F.

,

van der Kruit

P. C.

,

1997

,

A&A

,

327

,

966

Dehnen

W.

,

1993

,

MNRAS

,

265

,

250

10.1093/mnras/265.1.250

Crossref

Dotti

M.

,

Colpi

M.

,

Haardt

F.

,

2006

,

MNRAS

,

367

,

103

10.1111/j.1365-2966.2005.09956.x

Crossref

Dotti

M.

,

Colpi

M.

,

Haardt

F.

,

Mayer

L.

,

2007

,

MNRAS

,

379

,

956

10.1111/j.1365-2966.2007.12010.x

Crossref

Fiacconi

D.

,

Mayer

L.

,

Roškar

R.

,

Colpi

M.

,

2013a

,

ApJ

,

777

,

L14

10.1088/2041-8205/777/1/L14

Crossref

Fiacconi

D.

,

Mayer

L.

,

Roškar

R.

,

Colpi

M.

,

2013b

,

ApJ

,

777

,

L14

10.1088/2041-8205/777/1/L14

Crossref

Gaia Collaboration

et al. .,

2016

,

A&A

,

595

,

A1

10.1051/0004-6361/201629272

Crossref

Gaia Collaboration

et al. .,

2018

,

A&A

,

616

,

A1

10.1051/0004-6361/201833051

Crossref

Gradshteyn

I. S.

,

Ryzhik

I. M.

,

2007

,

Table of Integrals, Series, and Products

, 2th edn.

Elsevier/Academic Press

,

Amsterdam

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Granados

A.

,

Torres

D.

,

Castañeda

L.

,

Henao

L.

,

Vanegas

S.

,

2021

,

New Astron.

,

82

,

101456

Hashimoto

Y.

,

Funato

Y.

,

Makino

J.

,

2003

,

ApJ

,

582

,

196

10.1086/344260

Crossref

Hernquist

L.

,

1990

,

ApJ

,

356

,

359

10.1086/168845

Crossref

Hernquist

L.

,

1993

,

ApJS

,

86

,

389

10.1086/191784

Crossref

Hopkins

P. F.

,

2015

,

MNRAS

,

450

,

53

10.1093/mnras/stv195

Crossref

Just

A.

,

Peñarrubia

J.

,

2005

,

A&A

,

431

,

861

10.1051/0004-6361:20041175

Crossref

Just

A.

,

Khan

F. M.

,

Berczik

P.

,

Ernst

A.

,

Spurzem

R.

,

2011

,

MNRAS

,

411

,

653

10.1111/j.1365-2966.2010.17711.x

Crossref

Kashlinsky

A.

,

1986

,

ApJ

,

306

,

374

10.1086/164350

Crossref

Kuijken

K.

,

Gilmore

G.

,

1989

,

MNRAS

,

239

,

571

10.1093/mnras/239.2.571

Crossref

Li

K.

,

Bogdanović

T.

,

Ballantyne

D. R.

,

2020

,

ApJ

,

896

,

113

10.3847/1538-4357/ab93c6

Crossref

Michalski

K. A.

,

Mosig

J. R.

,

2016

,

J. Electromagn. Waves Appl.

,

30

,

281

Miyamoto

M.

,

Nagai

R.

,

1975

,

PASJ

,

27

,

533

Mulder

W. A.

,

1983

,

A&A

,

117

,

9

Navarro

J. F.

,

Frenk

C. S.

,

White

S. D. M.

,

1997

,

ApJ

,

490

,

493

10.1086/304888

Crossref

Ogata

H.

,

2005

,

Publications of the Research Institute for Mathematical Sciences

,

41

:

Ostriker

E. C.

,

1999

,

ApJ

,

513

,

252

10.1086/306858

Crossref

Petts

J. A.

,

Gualandris

A.

,

Read

J. I.

,

2015

,

MNRAS

,

454

,

3778

10.1093/mnras/stv2235

Crossref

Petts

J. A.

,

Read

J. I.

,

Gualandris

A.

,

2016

,

MNRAS

,

463

,

858

10.1093/mnras/stw2011

Crossref

Pfister

H.

,

Lupi

A.

,

Capelo

P. R.

,

Volonteri

M.

,

Bellovary

J. M.

,

Dotti

M.

,

2017

,

MNRAS

,

471

,

3646

10.1093/mnras/stx1853

Crossref

Pfister

H.

,

Volonteri

M.

,

Dubois

Y.

,

Dotti

M.

,

Colpi

M.

,

2019

,

MNRAS

,

486

,

101

10.1093/mnras/stz822

Crossref

Plummer

H. C.

,

1911

,

MNRAS

,

71

,

460

10.1093/mnras/71.5.460

Crossref

Read

J. I.

,

Goerdt

T.

,

Moore

B.

,

Pontzen

A. P.

,

Stadel

J.

,

Lake

G.

,

2006

,

MNRAS

,

373

,

1451

10.1111/j.1365-2966.2006.11022.x

Crossref

Rosas-Guevara

Y.

et al. ,

2020

,

MNRAS

,

491

,

2547

10.1093/mnras/stz3180

Crossref

Roškar

R.

,

Fiacconi

D.

,

Mayer

L.

,

Kazantzidis

S.

,

Quinn

T. R.

,

Wadsley

J.

,

2015

,

MNRAS

,

449

,

494

10.1093/mnras/stv312

Crossref

Souza Lima

R.

,

Mayer

L.

,

Capelo

P. R.

,

Bellovary

J. M.

,

2017

,

ApJ

,

838

,

13

10.3847/1538-4357/aa5d19

Crossref

Spitzer

,

Lyman

J.

,

1942

,

ApJ

,

95

,

329

10.1086/144407

Crossref

Springel

V.

,

Di Matteo

T.

,

Hernquist

L.

,

2005

,

MNRAS

,

361

,

776

10.1111/j.1365-2966.2005.09238.x

Crossref

Tamburello

V.

,

Capelo

P. R.

,

Mayer

L.

,

Bellovary

J. M.

,

Wadsley

J. W.

,

2017

,

MNRAS

,

464

,

2952

10.1093/mnras/stw2561

Crossref

Tamfal

T.

,

Mayer

L.

,

Quinn

T. R.

,

Capelo

P. R.

,

Kazantzidis

S.

,

Babul

A.

,

Potter

D.

,

2020

,

preprint (arXiv:2007.13763)

Tremaine

S.

,

Weinberg

M. D.

,

1984

,

MNRAS

,

209

,

729

10.1093/mnras/209.4.729

Crossref

Tremmel

M.

,

Governato

F.

,

Volonteri

M.

,

Quinn

T. R.

,

2015

,

MNRAS

,

451

,

1868

10.1093/mnras/stv1060

Crossref

Van Wassenhove

S.

,

Capelo

P. R.

,

Volonteri

M.

,

Dotti

M.

,

Bellovary

J. M.

,

Mayer

L.

,

Governato

F.

,

2014

,

MNRAS

,

439

,

474

10.1093/mnras/stu024

Crossref

Vasiliev

E.

,

2019

,

MNRAS

,

482

,

1525

10.1093/mnras/sty2672

Crossref

Verbiest

J. P. W.

et al. ,

2016

,

MNRAS

,

458

,

1267

10.1093/mnras/stw347

Crossref

Weinberg

M. D.

,

1986

,

ApJ

,

300

,

93

10.1086/163785

Crossref

White

S. D. M.

,

1983

,

ApJ

,

274

,

53

10.1086/161425

Crossref

Yurin

D.

,

Springel

V.

,

2014

,

MNRAS

,

444

,

62

10.1093/mnras/stu1421

Crossref

APPENDIX A: POTENTIAL AND ACCELERATION FOR THICK EXPONENTIAL DISCS

In this appendix, we detail the derivation leading to the form of potential and accelerations employed in this work. Our derivation follows those presented in Casertano (1983) and Kuijken & Gilmore (1989), with the exception that we express our results in terms of the Gauss hypergeometric function ₂F₁(a, b; c; z; see e.g. Abramowitz & Stegun 1972, for full details; see also Appendix B for the employed-specific implementation).

We consider a thick exponential disc with density profile given by:

$$\begin{eqnarray*} \rho (R,z) = \rho _{d,0} \ \rho _R(R) \ \rho _z(z), \end{eqnarray*}$$

(A1)

where

$$\begin{eqnarray*} \rho _R(R) = {\rm e}^{-\alpha R}, \end{eqnarray*}$$

(A2)

$$\begin{eqnarray*} \rho _z(z) = {\rm sech}^2\left(\frac{\beta z}{2}\right) = \cosh ^{-2}\left(\frac{\beta z}{2}\right), \end{eqnarray*}$$

(A3)

where ρ_d,0 represents the normalization constant (see equation 1), α = 1/R_d is the inverse of the scale radius, while β = 2/z_d is the inverse of the parameter shaping the vertical profile.

In order to obtain the gravitational potential ϕ generated by the above mass distribution, we need to solve the Poisson’s equation, i.e.

$$\begin{eqnarray*} \nabla ^{2} \phi (R,z) = 4\pi G \rho (R,z). \end{eqnarray*}$$

(A4)

When axial symmetry holds, equation (A4) can be solved in terms of Hankel’s transform, defined as

$$\begin{eqnarray*} \tilde{f}(k) = \int _0^\infty ~{\rm d}R \ R \ J_{0}(k R) f(R), \end{eqnarray*}$$

(A5)

$$\begin{eqnarray*} f(R) = \int _0^\infty ~{\rm d}k \ k \ J_{0}(k R) \tilde{f}(k), \end{eqnarray*}$$

(A6)

where J₀ is the 0-th order Bessel function of the first kind, and f(R) is a generic function.

By taking the Hankel’s transform of both sides of equation (A4), we obtain the linear non-homogeneous second order differential equation

$$\begin{eqnarray*} -k^2 \tilde{\phi }(k,z) + \frac{\partial ^2\tilde{\phi }}{\partial z}(k,z) = 4\pi G \ \tilde{\rho }(k,z), \end{eqnarray*}$$

(A7)

that, once solved through standard techniques, yields

$$\begin{eqnarray*} \tilde{\phi }(k,z) = - \frac{2\pi G}{k} \int _{-\infty }^{+\infty } ~{\rm d}\zeta \ {\rm e}^{-k |z-\zeta |} \tilde{\rho }(k,\zeta), \end{eqnarray*}$$

(A8)

provided that the Hankel transform of the density profile, |$\tilde{\rho }(k,z)$|⁠, vanishes when |z| → ±∞ (actually our case). Through equation (A6), the potential reads

$$\begin{eqnarray*} \phi (R,z) = - 2\pi G \int _0^\infty ~{\rm d}k J_0(k R) \int _{-\infty }^{+\infty } ~{\rm d}\zeta \ {\rm e}^{-k |z-\zeta |} \tilde{\rho }(k,\zeta). \end{eqnarray*}$$

(A9)

Since the density profile in equation (A1) is factorized in the radial and vertical part, the Hankel transform only affects the former, such that equation (A9) reads

$$\begin{eqnarray*} \phi (R,z) &=& - 2\pi G \rho _{d,0} \int _0^\infty ~{\rm d}k J_0(k R) \nonumber \\ &&\times \, \int _0^\infty ~{\rm d}u \ u J_0(uk)\rho _R(u) \int _{-\infty }^{+\infty } ~{\rm d}\zeta \ {\rm e}^{-k |z-\zeta |} \rho _z(\zeta). \end{eqnarray*}$$

(A10)

The u-integral can be evaluated analytically (see e.g. Kuijken & Gilmore 1989) and gives

$$\begin{eqnarray*} \int _0^\infty ~{\rm d}u \ u J_0(uk) {\rm e}^{-\alpha u} = \frac{\alpha }{(\alpha ^2+k^2)^{3/2}}, \end{eqnarray*}$$

(A11)

while the ζ-integral needs some additional manipulations. In particular, we can first split the integration domain and change the sign of the integration variable of the [−∞, z] part, obtaining

$$\begin{eqnarray*} I_z(k) = \int _{-z}^{+\infty } ~{\rm d}\zeta \ \frac{{\rm e}^{-k (z+\zeta)}}{\cosh ^2\left(\frac{\beta \zeta }{2}\right)} + \int _{z}^{+\infty } ~{\rm d}\zeta \ \frac{{\rm e}^{-k (\zeta -z)}}{\cosh ^2\left(\frac{\beta \zeta }{2}\right)}, \end{eqnarray*}$$

(A12)

that after the variable change ζ ± z = (2/β)y we reduce to

$$\begin{eqnarray*} I_z(k) &=& \frac{2}{\beta }\left[ \int _{0}^{+\infty } ~{\rm d}y \ \frac{{\rm e}^{-2ky/\beta }}{\cosh ^2\left(y-\frac{\beta z}{2}\right)}\right. \nonumber\\ &&\left.+\, \int _{0}^{+\infty } ~{\rm d}y \ \frac{{\rm e}^{-2ky/\beta }}{\cosh ^2\left(y+\frac{\beta z}{2}\right)} \right]. \end{eqnarray*}$$

(A13)

Finally, we use the symbolic software Mathematica to express it in a closed form in terms of the Gauss hypergeometric function ₂F₁, i.e.

$$\begin{eqnarray*} I_z(k) &=& \frac{4}{\beta }\Biggl \lbrace 1 - \frac{k}{k+\beta } \biggl [{\rm e}^{-z \beta} {}_{2}{F}_{1}\left(1,1+\frac{k}{\beta };2+\frac{k}{\beta };-{\rm e}^{-z \beta }\right) \nonumber \\ &&+\, {\rm e}^{z \beta} {}_{2}{F}_1\left(1,1+\frac{k}{\beta };2+\frac{k}{\beta };-{\rm e}^{z \beta }\right) \biggr ] \Biggr \rbrace . \end{eqnarray*}$$

(A14)

Going back to equation (A10), the final form of the gravitational potential can be written as

$$\begin{eqnarray*} \phi (R,z) = - 2\pi G \alpha \rho _{d,0} \int _0^\infty ~{\rm d}k J_0(k R) \frac{I_z(k)}{(\alpha ^2+k^2)^{3/2}}. \end{eqnarray*}$$

(A15)

To obtain the radial and vertical accelerations necessary to integrate the equations of motion, we take the (negative) gradient of equation (A15), obtaining

$$\begin{eqnarray*} a_R = -\frac{\partial \phi }{\partial R} &= - 2\pi G \alpha \rho _{d,0} \int _0^\infty ~{\rm d}k k J_1(k R) \frac{I_z(k)}{(\alpha ^2+k^2)^{3/2}}, \end{eqnarray*}$$

(A16)

$$\begin{eqnarray*} a_z = -\frac{\partial \phi }{\partial z} &= - 2\pi G \alpha \rho _{d,0} \int _0^\infty ~{\rm d}k J_0(k R) \frac{-\partial _z I_z(k)}{(\alpha ^2+k^2)^{3/2}}, \end{eqnarray*}$$

(A17)

where ∂_zI_z(k) is given by

$$\begin{eqnarray*} \partial _z I_z(k) &=& \frac{4 k^2 \text{sgn}(z)}{\beta (k+\beta)} \Biggl [ -{\rm e}^{-z \beta } {}_{2}{F}_1\left(1,1+\frac{k}{\beta };2+\frac{k}{\beta };-{\rm e}^{-z \beta }\right) \nonumber \\ &&+\, {\rm e}^{z \beta } {}_{2}{F}_1\left(1,1+\frac{k}{\beta };2+\frac{k}{\beta };-{\rm e}^{z \beta }\right) \nonumber \\ &&-\, \frac{k+\beta }{k} \tanh \left(\frac{z \beta }{2}\right)\Biggr ]. \end{eqnarray*}$$

(A18)

Equations (A15)–(A17) express the gravitational potential, the radial and vertical accelerations in form of one-dimensional integrals. Unfortunately, no further simplifications are possible and no closed form expressions can be found, therefore those integrals have to be evaluated numerically. Despite the apparent simplicity of dealing with one-dimensional integrals, the presence of the Bessel functions J₀ and J₁ can make the task quite challenging given their highly oscillatory behaviour, requiring therefore specific integration schemes to obtain sensible outcomes. Nevertheless, once this task is properly accomplished, it allows us to integrate the orbits with arbitrary initial conditions. Specifically, to numerically integrate such kind of oscillatory functions, we employ a special variant of the DE rule.¹¹

In the special case of equatorial motion (i.e. z = 0), equation (A14) can be simplified, reducing equation (A12) to

$$\begin{eqnarray*} I_0(k) = 2\int _{0}^{+\infty } ~{\rm d}\zeta \ \frac{{\rm e}^{-k \zeta }}{\cosh ^2\left(\frac{\beta \zeta }{2}\right)}, \end{eqnarray*}$$

(A19)

that with some manipulations evaluates to (see e.g. Gradshteyn & Ryzhik 2007, p. 383)

$$\begin{eqnarray*} I_0(k) = \frac{4}{\beta } \Biggl \lbrace \frac{k}{\beta }\biggl [ \psi \left(\frac{k}{2\beta }+\frac{1}{2}\right) - \psi \left(\frac{k}{2\beta }\right) \biggr ] - 1 \Biggr \rbrace , \end{eqnarray*}$$

(A20)

where the function ψ is the digamma function (see e.g. Abramowitz & Stegun 1972).

Finally, since the functions I_z(k), I₀(k), and ∂_zI_z(k) contain addition and subtractions of special functions, round-off error can represent a serious issue for the numerical evaluation of the integrals when k → +∞. To circumvent this possible problem, we can replace the expression of those functions with their asymptotic expansions for k → +∞, i.e.

$$\begin{eqnarray*} I_z(k) = \frac{8 \ {\rm e}^{-z \beta }}{(1+{\rm e}^{-z \beta })^2 k} + \mathcal {O}(1/k^3), \end{eqnarray*}$$

(A21)

$$\begin{eqnarray*} I_0(k) = \frac{2}{k} - \frac{\beta ^2}{k^3} + \mathcal {O}(1/k^5), \end{eqnarray*}$$

(A22)

$$\begin{eqnarray*} \partial _z I_z(k) = \frac{8 \beta \ {\rm e}^{-z \beta } (1-{\rm e}^{-z \beta })}{(1+{\rm e}^{-z \beta })^3 \ (k+\beta)} + \mathcal {O}(1/k^3). \end{eqnarray*}$$

(A23)

APPENDIX B: COMPUTATION OF THE GAUSS HYPERGEOMETRIC FUNCTION

The Gauss hypergeometric function ₂F₁(a, b; c; z) is a special function defined for |z| < 1 by the power series

$$\begin{eqnarray*} _2F_1(a,b;c;z) = \frac{\Gamma (c)}{\Gamma (a)\Gamma (b)} \sum _{n=0}^\infty \frac{\Gamma (a+n)\Gamma (b+n)}{\Gamma (c+n)} \frac{z^n}{n!}, \end{eqnarray*}$$

(B1)

and by analytic continuation elsewhere on the complex plane. In the above definition, Γ(·) represents the gamma function. The Gauss hypergeometric function is very general and includes as limiting cases several other mathematical functions. A detailed discussion about the properties and relations involving the hypergeometric function is definitely beyond the scope of this work, for which we are only interested in the specific case when the hypergeometric function assumes the form

$$\begin{eqnarray*} _2F_1(1,1+y;2+y;z), \end{eqnarray*}$$

(B2)

with |$z \in \mathbb {R}$| always negative and |$y \in \mathbb {R}$| and positive.

Several software implementations of the ₂F₁ function exist, but, at least to our knowledge, none of them is able to provide a sensible evaluation for all possible combinations of (y, z), for which reliable results are usually provided only when |z| < 1.¹² We therefore implemented our own algorithm to evaluate ₂F₁ for |z| > 1, relying instead on the GNU GSL mathematical library for |z| ≤ 1.

In order to efficiently compute ₂F₁, we first note that two limiting cases exist, i.e. when y → 0 and y → +∞. In these specific situations, ₂F₁ assumes the following special values

$$\begin{eqnarray*} _2F_1(1,1;2;z) &= -\frac{\ln (1-z)}{z}, \end{eqnarray*}$$

(B3)

$$\begin{eqnarray*} _2F_1(1,b;b;z) &= (1-z)^{-1}, \end{eqnarray*}$$

(B4)

allowing us to express ₂F₁ as one of the limiting value plus small corrections.

Specializing equation (B1) with parameters in equation (B2), we get

$$\begin{eqnarray*} _2F_1(1,1+y;2+y;z) = (1+y) \sum _{n=0}^\infty \frac{z^n}{1+y+n}. \end{eqnarray*}$$

(B5)

In the case of y → 0, we can expand the function as

$$\begin{eqnarray*} _2F_1(1,1+y;2+y;z) = (1+y) \sum _{k=0}^\infty (-1)^k y^k \sum _{n=0}^\infty \frac{z^n}{(1+n)^{k+1}}, \end{eqnarray*}$$

(B6)

that, changing the summation index to m = n + 1, yields

$$\begin{eqnarray*} _2F_1(1,1+y;2+y;z) = \frac{1+y}{z} \sum _{k=0}^\infty (-1)^k y^k \sum _{m=1}^\infty \frac{z^m}{m^{k+1}}. \end{eqnarray*}$$

(B7)

where the series over m defines the polylogarithm function (Abramowitz & Stegun 1972) of order (k + 1), denoted Li_{k + 1}(z). Starting from the first order of the function Li₁(z) = −ln (1 − z)/z, all lower (i.e. the 0-th and negative) orders can be expressed in closed form exploiting the recurrence relation

$$\begin{eqnarray*} z \frac{\partial Li_s}{\partial z}(z) = Li_{s-1}(z). \end{eqnarray*}$$

(B8)

On the other hand, above the second positive order (included), a closed form cannot be found. Exploiting the polylogarithm function, equation (B7) can be written as

$$\begin{eqnarray*} _2F_1(1,1+y;2+y;z) &=& -\frac{1+y}{z}\ln (1-z) \nonumber\\ &&+\, \frac{1+y}{z} \sum _{k=1}^\infty (-1)^k y^k Li_{k+1}(z). \end{eqnarray*}$$

(B9)

Depending on the value of y, we can then truncate the series to match a desired accuracy. Finally, by setting y = 0 in the above equation, we can check that it correctly simplifies into equation (B3).

On the opposite limit, i.e. when y → +∞, we follow the very same procedure outlined above, this time expanding the term containing y as

$$\begin{eqnarray*} \frac{1+y}{1+y+n} = 1 + \sum _{l=1}^\infty (-1)^l \frac{n (1+n)^{l-1}}{y^l}, \end{eqnarray*}$$

(B10)

that, inserted into equation (B5), gives

$$\begin{eqnarray*} _2F_1(1,1+y;2+y;z) = \sum _{n=0}^\infty z^n + \sum _{n=0}^\infty z^n \sum _{l=1}^\infty (-1)^l \ \frac{n(1+n)^{l-1}}{y^l}. \nonumber\\ \end{eqnarray*}$$

(B11)

The first term involving only powers of z is the series definition for 1/(1 − z). The second term, instead, can be further expanded exploiting the binomial theorem, leading to

$$\begin{eqnarray*} &&{{}_2F_1(1,1+y;2+y;z)}\nonumber \\ &&{\quad= (1-z)^{-1} + \sum _{l=1}^\infty \frac{(-1)^l}{y^l} \sum _{k=0}^{l-1}\binom{l-1}{k} \sum _{n=1}^\infty n^{k+1} z^n,}\nonumber \\ &&{\quad= (1-z)^{-1} + \sum _{l=1}^\infty \frac{(-1)^l}{y^l} \sum _{k=0}^{l-1}\binom{l-1}{k} \ Li_{-(k+1)}(z),} \end{eqnarray*}$$

(B12)

where, in the third line, we recognized the definition of the polylogarithm function. Since the polylogarithm function appears with negative orders, each of them can be easily evaluated through the recurrence relation in equation (B8).

Finally, in all cases for which the limiting forms for y → 0 or y → +∞ are not suitable, we evaluate the hypergeometric function using its series definition in equation (B1). Formally, this provides a sensible result only when |z| < 1, but exploiting the linear transformation property of ₂F₁ (Abramowitz & Stegun 1972) and performing the transformation z → 1/z, an alternative expression of ₂F₁ as a function of 1/z can be obtained. Once reduced and specialized to the case of interest, it reads

$$\begin{eqnarray*} _2F_1(1,1+y;2+y;z) = \frac{1+y}{z}\sum _{n=0}^\infty \frac{1}{(n-y)z^n} - \frac{\pi (1+y)}{\sin (\pi y) \ z^{1+y}}. \nonumber\\ \end{eqnarray*}$$

(B13)

APPENDIX C: USEFUL EXPRESSIONS FOR DF WITH FAST MOVING STARS

Here, we provide the analytical form of the integral over the variable V that appears in equation (8), i.e.

$$\begin{eqnarray*} J(v_p,v_\star ,p_{\rm max}) = \int _{|v_p-v_\star |}^{v_p+v_\star } ~{\rm d}V \left(1+\frac{v_p^2-v_\star ^2}{V^2}\right) \ln \left(1+\frac{p_{\rm max}^2 V^4}{G^2 m_p^2}\right). \nonumber\\ \end{eqnarray*}$$

(C1)

where v_p and v_⋆ are the MP and star velocities, whereas p_max is the maximum impact parameter. To compute the integral, we need to consider the two different cases in which v_p > v_⋆ and v_p < v_⋆. For convenience, we express the result in terms of the Gauss hypergeometric function F(•) = ₂F₁(1, 5/4, 9/4, •).¹³

When v_p > v_⋆, the integration yields

$$\begin{eqnarray*} &&{\frac{1}{8 v_\star } \int _{v_p-v_\star }^{v_p+v_\star } ~{\rm d}V \left(1+\frac{v_p^2-v_\star ^2}{V^2}\right) \ln \left(1+\frac{p_{\rm max}^2 V^4}{G^2 m_p^2}\right)} \nonumber \\ &&{\quad = \frac{1}{20 v_\star }\Biggl [ -2\left(\frac{G m_p}{p_{\rm max}}\right)^2 \frac{v_p+v_\star }{(v_p-v_\star)^4} F(\alpha) +2\left(\frac{p_{\rm max}}{G m_p}\right)^2 }\nonumber \\ &&{\qquad\times \,(v_p-v_\star)^5 F(\beta){+}2\left(\frac{G m_p}{p_{\rm max}}\right)^2 \frac{v_p{-}v_\star }{(v_p{+}v_\star)^4} F(\gamma) {-}2\left(\frac{p_{\rm max}}{G m_p}\right)^2 } \nonumber \\ &&{\qquad(v_p+v_\star)^5 F(\delta) +5v_\star \Biggl (4 + \ln \Bigl (1 + \frac{p_{\rm max}^4}{G^4 m_p^4}(v_p^2-v_\star ^2)^4} \nonumber\\ &&{\qquad+ \frac{2 p_{\rm max}^2}{G^2 m_p^2} (v^4+6v^2v_\star ^2+v_\star ^4) \Bigr) \Biggr) \Biggr ],}\nonumber\\ \end{eqnarray*}$$

(C2)

while, when v_p < v_⋆, we obtain

$$\begin{eqnarray*} &&{\frac{1}{8 v_\star } \int _{-v_p+v_\star }^{v_p+v_\star } ~{\rm d}V \left(1+\frac{v_p^2-v_\star ^2}{V^2}\right) \ln \left(1+\frac{p_{\rm max}^2 V^4}{G^2 m_p^2}\right)} \nonumber \\ &&{\quad = \frac{1}{10 v_\star }\Biggl [ -10 v_p +\left(\frac{G m_p}{p_{\rm max}}\right)^2 \frac{v_p+v_\star }{(v_p-v_\star)^4} F(\alpha) -\left(\frac{p_{\rm max}}{G m_p}\right)^2 }\nonumber \\ &&{\qquad (v_p-v_\star)^5 F(\beta)+\left(\frac{G m_p}{p_{\rm max}}\right)^2 \frac{v_p-v_\star }{(v_p+v_\star)^4} F(\gamma) -\left(\frac{p_{\rm max}}{G m_p}\right)^2 } \nonumber \\ &&{\qquad (v_p+v_\star)^5 F(\delta) +\frac{5}{2}v_\star \ln \left(\frac{1 + \frac{p_{\rm max}^2}{G^2 m_p^2}(v_p+v_\star)^4}{1 + \frac{p_{\rm max}^2}{G^2 m_p^2}(v_p-v_\star)^4} \right) \Biggr ],} \end{eqnarray*}$$

(C3)

where in both expression we have defined

$$\begin{eqnarray*} \alpha = -\frac{G^2 m_p^2}{p_{\rm max}^2 (v_p-v_\star)^4}, \end{eqnarray*}$$

(C4)

$$\begin{eqnarray*} \beta = -\frac{p_{\rm max}^2 (v_p-v_\star)^4}{G^2 m_p^2}, \end{eqnarray*}$$

(C5)

$$\begin{eqnarray*} \gamma = -\frac{G^2 m_p^2}{p_{\rm max}^2 (v_p+v_\star)^4}, \end{eqnarray*}$$

(C6)

$$\begin{eqnarray*} \delta = -\frac{p_{\rm max}^2 (v_p+v_\star)^4}{G^2 m_p^2}. \end{eqnarray*}$$

(C7)

We note that, if we assume the logarithmic term in equation (C3) to be independent of V, then the resulting integral vanishes, and the equation simplifies to the well known result in which only particles moving slower than the MP contribute to DF.

APPENDIX D: ADDITIONAL TEST CASES

Here, we show some additional test cases of the evolution computed within our framework and compared to a corresponding N-body simulation. In particular, in Fig. D1, we consider low inclination orbits with the initial z coordinates ranging from 0.5 kpc (top left) to 2.0 kpc (bottom right). For each case, we show the time evolution of the radial separation that of the angular momentum normalized to its staring value and that of the inclination. As for the cases commented in Section 3, we obtain a fairly good agreement when the limited spatial resolution of the N-body simulation is taken into account and also when the tangential velocity profiles is used. As previously noted, the semi-analytical evolution reproduces quite well the eccentric inspiral of the MP, with a slightly larger dephasing only when the orbit tends to become more circular. This could hint to the fact that, despite the reasonably good recovered evolution, the assumption of locality made in our framework could be less reasonable as the orbit circularizes and that second order effects, connected with the complexity and non-linearity of DF, could start to play a role.

Figure D1.

Time evolution of radial separation, z component of the angular momentum (normalized to its staring value) and orbital inclination for four cases featuring an increasing initial z_max ranging from 0.5 kpc (top left) to 2 kpc (bottom right). In all plots, the semi-analytical evolution is shown with solid blue lines, while the comparison with the corresponding N-body simulation is given by the dashed orange lines.

Open in new tab Download slide

Finally, note that the large variations in inclination that can be observed at later times are due to spurious noise determined by the sparser particle number in the central regions.

This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://dbpia.nl.go.kr/journals/pages/open_access/funder_policies/chorus/standard_publication_model)

Download all slides

Month:	Total Views:
January 2021	1
February 2021	10
March 2021	13
April 2021	6
May 2021	8
June 2021	5
July 2021	2
August 2021	4
September 2021	7
October 2021	9
December 2021	2
January 2022	3
March 2022	1
April 2022	2
June 2022	5
July 2022	3
August 2022	3
September 2022	1
October 2022	6
February 2023	3
March 2023	1
April 2023	2
May 2023	2
June 2023	3
July 2023	1
September 2023	2
November 2023	2
December 2023	3
January 2024	3
February 2024	7
March 2024	11
April 2024	5
May 2024	18
June 2024	8
July 2024	9
August 2024	5
September 2024	12
October 2024	8
November 2024	5
December 2024	31
January 2025	6
February 2025	2
March 2025	10
April 2025	1
May 2025	6

Article Contents

Dynamical evolution of massive perturbers in realistic multicomponent galaxy models I: implementation and validation

ABSTRACT

1 INTRODUCTION

2 METHOD

2.1 Exponential disc conservative dynamics

2.2 DF in spherical profiles

2.3 DF in rotating discs

3 RESULTS

3.1 Orbits in isolated exponential discs

3.1.1 Consistency checks

3.1.2 Orbit integration in exponential discs

3.2 DF in multicomponents models

4 DISCUSSION AND CONCLUSIONS

ACKNOWLEDGEMENTS

DATA AVAILABILITY STATEMENT

Footnotes

REFERENCES

APPENDIX A: POTENTIAL AND ACCELERATION FOR THICK EXPONENTIAL DISCS

APPENDIX B: COMPUTATION OF THE GAUSS HYPERGEOMETRIC FUNCTION

APPENDIX C: USEFUL EXPRESSIONS FOR DF WITH FAST MOVING STARS

APPENDIX D: ADDITIONAL TEST CASES

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

Article Contents

Dynamical evolution of massive perturbers in realistic multicomponent galaxy models I: implementation and validation Free

ABSTRACT

1 INTRODUCTION

2 METHOD

2.1 Exponential disc conservative dynamics

2.2 DF in spherical profiles

2.3 DF in rotating discs

3 RESULTS

3.1 Orbits in isolated exponential discs

3.1.1 Consistency checks

3.1.2 Orbit integration in exponential discs

3.2 DF in multicomponents models

4 DISCUSSION AND CONCLUSIONS

ACKNOWLEDGEMENTS

DATA AVAILABILITY STATEMENT

Footnotes

REFERENCES

APPENDIX A: POTENTIAL AND ACCELERATION FOR THICK EXPONENTIAL DISCS

APPENDIX B: COMPUTATION OF THE GAUSS HYPERGEOMETRIC FUNCTION

APPENDIX C: USEFUL EXPRESSIONS FOR DF WITH FAST MOVING STARS

APPENDIX D: ADDITIONAL TEST CASES

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only

Dynamical evolution of massive perturbers in realistic multicomponent galaxy models I: implementation and validation