Dynamical friction-driven orbital circularisation in rotating discs: a semi-analytical description

We present and validate a novel semi-analytical approach to study the effect of dynamical friction on the orbits of massive perturbers in rotating stellar discs. We find that dynamical friction efficiently circularises the orbit of co-rotating perturbers, while it constantly increases the eccentricity of counter-rotating ones until their angular momenta reverse, then once again promoting circularisation. Such"drag toward circular corotation"could shape the distribution of orientations of kinematically decoupled cores in disc galaxies, naturally leading to the observed larger fraction of co-rotating cores.


INTRODUCTION
The orbital evolution of massive perturbers (MPs, i.e. objects considerably heavier than the single bodies forming the galactic structure within which the perturber is moving) under the effect of dynamical friction (DF, see Binney & Tremaine 2008) has been studied in many different contexts and over a broad range of spacial scales: from the pairing and mergers of galaxies at kpc scales (see De Rosa et al. 2020, for a recent review), through the sinking and tidal disruption of globular clusters in the central kpc of galaxies (Binney & Tremaine 2008), down to the formation of massive black hole (MBH) binaries at pc-scales (MBHBs, see e.g.Dotti et al. 2012, for a review).
The interaction between a MP and a rotating system embedding it (either stellar or gaseous) has been found to reverse the direction of the MP angular momentum if it is initially anti-aligned with that of the background, or, if the two are aligned, to circularise the MP orbit.Such behaviour has first been observed by Dotti et al. (2006) in simulations of MBHs in circumnuclear gaseous discs, common in ultraluminous infrared galaxies (Sanders & Mirabel 1996;Downes & Solomon 1998), and observed to form in gas rich galaxy mergers (Barnes 2002;Mayer et al. 2007;Capelo & Dotti 2017).This behaviour, dubbed here "drag toward circular corotation" 1 (DTCC hereafter), has been demonstrated to be independent of the nature of the background, working in nuclear stellar discs as well (Dotti et al. 2007), and then observed to be present in larger discs during galaxy minor mergers (Callegari et al. 2011;Fiacconi et al. 2013).Dotti et al. (2006) proposed a qualitative model in which DF is responsible for the DTCC, by accelerating the MP in the direction opposite to its instantaneous motion relative to the local environment.DF always decelerates counter-rotating MPs, increasing their eccentricity without changing the orbital plane, until the DF effect integrated over one orbit suffices to reverse the MP motion; DF instead increases the angular momentum of co-rotating MPs at the apocentre (where the MP is slower than the background) and decreases it at the pericentre, circularising the orbit.
Contrary to the cases of DF exerted by spherical structures with isotropic velocity fields (e.g.Chandrasekhar 1 The angular momentum reversal of initially retrograde orbits was originally referred to as "angular momentum flip" (e.g.Dotti et al. 2006).This convention could be misleading, since the orbital plane of the orbit does not evolve gradually passing through a perpendicular configuration.For this reason we decided to adopt the new conventional name.
1943), an analytical description of DF in composite systems with at least one rotating disc-like structure (such as a disc galaxy, or a nuclear disc embedded in the larger bulge) has not been developed yet.As a consequence, the model put forward by Dotti et al. (2006) and adopted as physical explanation for the DTCC by the all the above-mentioned studies has not been proved to date.Such proof is long due, because of the numerous and important consequences of the DTCC, among which: the slow-down of the MP orbital shrinking (see e.g. the discussion in Mayer 2013), the enhanced accretion onto orbiting MBHs once they circularise (Dotti et al. 2009;Callegari et al. 2011), the alignment of the spins of MBHs in circumnuclear discs, well before they bind in a binary (Dotti et al. 2010) and the consequent low recoils expected at coalescence (but for a small tail of fast recoiling remnants, see Lousto et al. 2012).
In this paper we describe a simple analytical form for the dynamical friction on MPs in rotating systems (Section 2), and integrate it numerically in composite systems including discs.We compare the outcome of this new semianalytical description with the results of high resolution Nbody simulations of the same galactic models (Section 3).Finally, as discussed in Section 4, the agreement between the two numerical tests demonstrates that DF is indeed the only responsible for the DTCC in rotating systems.In the same section we comment on the possible implications of the DTCC for kinematically decoupled cores, and discuss some future expansions and applications of our model.

METHOD
Employing a C++ implementation of the Bulirsch-Stoer (Bulirsch & Stoer 1966) integration scheme,2 we numerically integrate the differential equations describing the motion of a MP when subjected to the global conservative gravitational potential φ generated by the galaxy mass distribution, as well as because of the action of DF.
In order to obtain analytical equations of motion, we model the host galaxy through appropriate density-potential pairs (ρ/φ), in particular we choose: • a Navarro-Frenk-White profile (NFW, Navarro et al. 1997), with mass M h and scale radius a h , to describe the dark matter (DM) halo (1) • an Hernquist profile (Hernquist 1990), with mass M b and scale radius a b , to model a compact stellar bulge • an exponential disc profile (see e.g.Binney & Tremaine 2008), with mass M d , scale radius R d and scale height z d , to represent the galactic stellar disc where R = x 2 + y 2 , w = R/(2R d ), while I, K are modified Bessel functions of the first and second kind respectively.The disc potential is known only in the disc equatorial plane, where the orbits of the perturbers are set.We note that the first term in equation ( 6) is strictly valid for a razor-thin disc, while the second term is an approximate correction that takes into account the non-zero thickness of the exponential disc (see e.g.Kuijken & Gilmore 1989).
Under the above assumptions, the total mass density and total gravitational potential (in the z = 0 plane) are simply given by and the particle m p at position r p = (x p , y p , z p ) with velocity v p = (v p x , v p y , v p z ) is therefore subjected to a conservative acceleration In addition to the conservative gravitational force, the motion of the MP is also affected by the "dissipative" DF force. 3In the case of spherically symmetric non-dissipative systems, the effect of DF is often modeled via the Chandrasekhar (1943) formula: where ) is a dimensionless parameter relating the perturber velocity to the isotropic velocity dispersion σ, while ln(1 + Λ 2 ) (or more frequently ln Λ) is commonly known as the Coulomb logarithm, with Λ = p max /p min being the ratio between the maximum and minimum impact parameters.The above expression for DF is derived in the simplistic assumption of a MP moving in an infinite and homogeneous density background whose velocity distribution is strictly Maxwellian.In spite of this, equation (10) has been proven to reproduce the inspiral of a MP surprisingly well for a vast range of spherically symmetric and isotropic density profiles, once the radius-dependent local velocity dispersion σ(r) is plugged in the equation (see, e.g.Tremaine & Weinberg 1984).
Recently, a series of investigations (Hashimoto et al. 2003;Just & Peñarrubia 2005;Just et al. 2011;Petts et al. 2015Petts et al. , 2016) ) showed that the adoption of position dependent minimum and maximum impact parameters further improves the agreement with MP orbital decays obtained in full N-body simulations.Motivated by this, we implement equation ( 10) separately for the halo and bulge components, adopting position-dependent velocity dispersion σ(r) x and mass density ρ(r) x (where x refers either to the halo or to the bulge) and the following expressions for the maximum and minimum impact parameters: where γ x represents the logarithmic slope of the density profiles of the two spherical components, r is the radial coordinate and D p is the physical radius of the MP, which is zero for the case of an MBH.Note that, in the case of multiple component systems, as the one we present in this paper, the value of σ(r) x should be formally derived self-consistently accounting for all components.However, the comparison with the results of a full N-body simulation discussed in the next section demonstrates that, for the parameters we consider, our simple approximation is sufficient to recover the main features of the perturber orbital evolution.Contrary to the isotropic cases, the DF acceleration due to a rotating disc component is not necessarily opposed to the perturber direction of motion with respect to the centre of mass of the whole system, but, as discussed in the introduction, primarily depends on the direction of the relative velocity vector between the perturber and the surrounding medium.In fact, as long as the orbit is counter-rotating or, if co-rotating and eccentric, DF acts instantaneously in the direction opposite to the velocity of the perturber with respect to the surrounding matter in its proximity, making the dependence on the velocity direction more relevant than that on the specific shape of the density profile. 4Being interested in confirming the DTCC effect, we limit our investigations to non-circular co-rotating cases, and therefore we make the simplifying assumption that the DF in the disc is primarily determined by the local distribution of matter and velocity around the MP.In particular, in the expression: we assume that the distribution function f (v a ) is a delta function around the circular velocity at the position of the MP, i.e. yielding where Due to the paucity of detailed studies on dynamical friction in rotating systems, a set of equations like (11) for discs has not been discussed in literature so far.Nevertheless, in order to obtain a simple form for the acceleration, we introduce a physically motivated prescription, based on the Coulomb logarithm, similar in spirit to the spherical case.In particular, we define Λ = p max,d /p min,d , with the maximum impact parameter taken equal to the scale-height of the disc, as the underlying assumption of quasi-homogeneous three-dimensional density fails on scales larger than z d .The minimum impact parameter is defined instead as an effective influence radius based on the relative velocity The additional term 0.01v 2 circ at the denominator is added only for purely operational purposes.It represents a numerical floor to p min,d that avoids its divergence and therefore a spurious suppression of the DF force due to the disc once the circular velocity is met. 5  In addition to the above prescription and trying to account for the other possible dependencies of DF that cannot be accounted for in this simple framework (e.g. the faster decrease of density in the vertical direction w.r.t. the radial one), we introduce an additional tunable constant A disc in equation ( 14), yielding to The value of the constant A disc is determined by comparing the semi-analytical orbital evolution of the perturber with that obtained in an N-body simulation, as discussed in Section 3.
Finally, including all terms presented above, the total acceleration acting on m p reads which allows us to numerically integrate the orbital motion in the galactic disc.

RESULTS: COMPARISON WITH N-BODY SIMULATIONS
To validate our analytic prescription we perform highresolution N-body simulations of an isolated galaxy with the code gizmo (Hopkins 2015), descendent of gadget2 (Springel 2005), where gravity is solved via a multipole hierarchical method based on a Barnes-Hut Oct-Tree structure, that guarantees an optimal balance between accuracy and performance.
Our galaxy model (used both in the semi-analytical calculation and in the N-body runs) includes a dark matter halo, described in the semi-analytical framework by a NFW profile and modeled as a Hernquist (1990) model in the Nbody run.In the latter case, the scale radius is set according to the halo concentration parameter, assumed to be c = 10.A stellar exponential disc and a central spherical bulge (also 5 Actually, if p min,d diverges then Λ 2 → 0 making ln(1+Λ 2 ) vanish.We checked that our setup is practically insensitive to any change up to an order of magnitude, both up and downwards, of the factor suppressing v circ .modelled as a spherical Hernquist profile) are embedded in the DM component.The model parameters are summarized in Tab. 1.
The choice of a massive bulge component was made mainly to enforce the Toomre-stability criterion over the whole stellar disc (as can be seen from the Toomre Q parameter profile in Fig. 1), thus suppressing the formation of strong local over-densities.Such features would add a stochastic forcing onto the perturber any time a close encounter takes place, hence contaminating the net DTCC effect always present whenever the background density is rotating (see Fiacconi et al. 2013;del Valle et al. 2015;Roškar et al. 2015;Lupi et al. 2016;Tamburello et al. 2017;Souza Lima et al. 2017).Our choice minimizes such stochastic perturbations and therefore allows for a clean test of the hypothesis that the DTCC effect is only due to DF.
We create the initial conditions for the N-body runs using GalIC (Yurin & Springel 2014), an iterative tool to create galaxy models in collisionless equilibrium, in which the particle orbits are integrated at each iteration for about ten orbital times.
We then set a MP with mass M MP = 10 8 M at a separation r MP = 5 kpc from the galaxy centre, moving along an elliptic orbit in the disc equatorial plane with eccentricity e MP = 0.7.We consider two orbital configurations, with the MP either co-rotating or counter-rotating relative to the stellar disc.The mass and spatial resolution of the N-body simulations are chosen to guarantee that the dynamical friction onto the MP is properly resolved (see the discussion in Pfister et al. 2017), as shown in Tab. 1.
In order to tune the free parameter A disc , we performed several tests to minimize the differences between the semianalytical integration and the full N-body simulations, finding that the best value for A disc is ≈ 0.5.This best fit value, of the order of unity, suggests that our modeling catches all the major physical properties of disc DF.Hence, the final form of the disc DF acceleration used to obtain the results presented in this section is: with Λ = p max,d /p min,d . 6Fig. 2 and Fig. 3 compare the orbital evolution obtained via the full N-body runs (dashed orange lines) with our semianalytical prescription (solid blue lines).In both figures, the left-hand panels show the perturber radial separation 7 (upper panel), its angular momentum normalised to its start-6 During the perturber evolution we found that the average value of the factor containing the Coulomb logarithm (i.e.A disc 2π ln(1 + Λ 2 )) lies around 20-22 for the prograde case and around 28-32 for the retrograde case.The difference can be ascribed to the higher relative velocity attained in the retrograde case.Additionally, in order to test the robustness of the employed prescription, we also tried to repeat the fitting procedure for A disc with a slightly different galaxy with z d /R d = 0.15 rather than 0.2 (see Tab. 1).Again we found that the best value for A disc is ≈ 0.5.Moreover, we also confirmed that Λ remains practically unchanged with respect to the z d /R d = 0.2 disc. 7The high resolution of the N-body run ensures that the perturber orbit remains close to the disc equator.Evolution of an initially prograde perturber.The left-hand panels, from top to bottom, refer respectively to the perturber separation, angular momentum (normalized to its initial magnitude) and eccentricity as a function of time.The right-hand panel shows the perturber orbit in the galactic equatorial plane.In all panels, our semi-analytical approach is shown as a blue solid line and the N-body run as an orange dashed line.ing value (middle panel) and its orbital eccentricity 8 (lower 8 As long as the perturber lies in the equatorial plane, the potential depends on the radial coordinate only.This allows us to employ well known results valid for central potentials.At each time step we evaluate the total energy E and angular momentum h (per unit mass), and assuming them as instantaneously conserved, we obtain an expression for the radial velocity, i.e.
Setting v r = 0, we can then numerically solve for r to obtain the two roots that characterise bound orbits, i.e. the values of panel) as a function of time.The single right-hand panel instead reports the perturber motion in the galactic equatorial plane.Fig. 2 shows the orbital evolution for a prograde perturber, while in Fig. 3 we consider an initially retrograde case.In spite of our simplistic assumptions, the semianalytical evolution closely matches that of the full N-body simulation.This is particularly evident for the perturber radial separation, orbital eccentricity and motion in the disc for the pericentre and apocentre, from which the eccentricity is ultimately evaluated.
plane.The angular momentum evolution is initially well reproduced, but, especially for the prograde case, it starts to deviate at later times.This behaviour is due to the assumption of locality made in our framework, which allows us to safely employ the semi-analytical prescription as long as the MP orbit is eccentric or, if circular, counter-rotating, while it becomes unreliable when the perturber orbit approaches a co-rotating circular one, i.e. when the relative velocity between the perturber and the disc at the perturber location goes to zero and equation ( 19) diverges.This limitation of our semi-analytical prescription is the reason why we limit ourselves to the study of the DTCC, and we do not focus on reproducing the whole evolution down to the centre of the composite system.Nevertheless, we can note that the orbital evolution is very well reproduced up to the point when the eccentricity approaches zero.Remarkably, as clear from the bottom left panels of Fig. 2 and Fig. 3, both for the prograde and retrograde cases, the time for which e → 0 is practically the same of the corresponding N-body simulation.
Circularisation of initially retrograde orbits takes longer (by about a factor of two) with respect to the prograde case (as shown by the comparison between Fig. 2 and Fig. 3).The difference in the timescales is due to the fact that, in order to circularise the MP orbit, DF has first to reverse the MP angular momentum, forcing the retrograde configuration into a prograde one.Still, despite the simplistic assumptions made, the time at which the sign reversal happens turns out to be quite similar to that of the N-body run and even more remarkably the eccentricity evolution gets closely reproduced up to the time of circularisation.We stress again that our goal is to model the motion of eccentric perturbers (or circular but counter-rotating), where velocities relative to the stellar background are non-negligible.In this regime our approach closely follows the outcomes of N-body simulations confirming that the angular momentum reversal, as the circularisation, is determined by DF.

DISCUSSION AND CONCLUSIONS
In this work we considered a simple semi-analytical approach to describe the motion of MPs inside rotating axisymmetric structures.We validated our prescription against full Nbody simulations, finding a remarkable good agreement.We verified that DF drives the MP motion towards circular orbits, even when initially the MP counter-rotates.In this last case, first the MP angular momentum gets reversed causing the MP to co-rotate with the disc, then, once on a prograde orbit, DF circularises it.
Here, we focused on rotating flat structures, but this feature does not specifically require a disc-like geometry.The DTCC effect is indeed common in any rotating structure and it could play an influential role even in fast rotating ellipticals (Emsellem et al. 2011) or galactic nuclei (see e.g. the discussion in Sesana et al. 2011;Rasskazov & Merritt 2017).Still, since a proper analytical description for such cases is not as easy as for the current case, we do not explore that further, limiting our discussion to the disc case only.
In addition to the consequences of the DTCC effect for massive black hole binaries (discussed in the Introduction and in the references therein), the DTCC also has a direct impact on the properties of kinematically decoupled stellar cores (KDCs).KDCs are central stellar components with a rotation axis that is misaligned from the main stellar body of the galaxy.They are characterised by an abrupt change of more than 30 degrees in the kinematic position angle and are observed in a significant fraction of early-type galaxies (∼ 7% Krajnović et al. 2011).Their presence is a clear signature of an external accretion event, with smaller KDCs (diameter < 1 kpc) associated with more recent accretion events, either caused by gas accretion followed by subsequent in-situ star formation or caused by a minor merger (McDermid et al. 2006;Raimundo et al. 2013).
In general, one could expect that if a minor merger occurs with a random orientation, the probability of observing decoupled cores that co-rotate or counter-rotate with respect to the main stellar body of the galaxy would be approximately similar.However, this does not seem to be the case, at least for the distribution of externally accreted gas.Davis et al. (2011) find a larger fraction of co-rotating than counter-rotating gas in early type galaxies, specifically in a ∼ 65/35 ratio between co-rotating gas and gas with large misalignment (> 30 deg).It is still not clear why such trend is observed, but some suggested reasons are that minor merger accretion may be anisotropic, occurring along a preferential direction such as the major axes of the dark matter halo (Deason et al. 2011;Davis et al. 2011), or that there is a drag towards co-rotation due to the interaction with already present gas (e.g. in the hot halo of the galaxy, Davis et al. 2011).If a KDC forms from the gas that has been externally accreted, then the statistics for the stellar core misalignment should be the same as for the gas since they would share similar dynamics.
The work presented here suggests an alternative mechanism to promote the higher relative occurrence of co-rotating KDSc.Our work predicts that for galaxies with discs, such as S0s and spiral galaxies, DF can result in the inversion of the angular momentum of a massive perturber, from a counter-rotating orbit to a co-rotating orbit.This is especially relevant for minor mergers, where we expect the accreted satellite galaxy to be a massive perturber and at a later stage to form a KDC.If one considers that the original direction of the minor merger is randomly oriented, this angular momentum inversion would result in a larger fraction of co-rotating KDCs than counter-rotating KDCs formed after minor merger events.The consequence of this alternative mechanism is an overall tendency for co-rotating KDCs in galaxies with discs, and its relevance can be observationally tested with the study of larger samples of KDCs and with the characterization of the host IGM properties.The fact that the DF mechanism studied here operates even in extremely gas poor systems, makes it a an important mechanism to consider even in early-type galaxies with low native gas content.In addition, after circularisation, DF removes efficiently angular momentum leading the perturber to the host nucleus without requiring a finely tuned initial angular momentum to justify the fall of externally accreted material down to the inner kpc of its host.
In analogy with KDCs, DTCC could have a similar impact on the kinematic properties of nuclear star clusters.Such systems have been proposed to form either via in-situ star formation (e.g.Loose et al. 1982), and/or via the accretion of multiple stellar clusters reaching the centre of the system via DF (e.g.Tremaine et al. 1975).The in-situ chan-nel has been always regarded as more promising if one wants to explain the kinematic properties of nuclear star clusters, and in particular their significant rotation.On the other hand, recently, Tsatsi et al. (2017) showed that even the formation via clusters accretion can bring to a net significant rotation, even if the orbits of the inspiralling clusters are isotropic.Our result strengthens their findings as, if we account for DTCC, clusters coming from the galactic disc should be accreted with the same direction in angular momentum (if their circularisation occurs before the cluster gets tidally disrupted).
Finally, in a future study we plan to expand our semianalytical treatment to the case of circular perturbers (corresponding to the final configuration of any initial configuration as a consequence of the DTCC effect) and to outof-plane geometries.This prescription will be instrumental in order to fully characterize the DF time distribution of infalling perturbers as a function of their host galaxy properties.

Figure 1 .
Figure 1.Toomre Q parameter for the stellar disc component in our galaxy model ICs.The black dashed line corresponds to the minimum Q to avoid strong non-axisymmetric perturbations developing in the disc.
Figure2.Evolution of an initially prograde perturber.The left-hand panels, from top to bottom, refer respectively to the perturber separation, angular momentum (normalized to its initial magnitude) and eccentricity as a function of time.The right-hand panel shows the perturber orbit in the galactic equatorial plane.In all panels, our semi-analytical approach is shown as a blue solid line and the N-body run as an orange dashed line.

Figure 3 .
Figure 3. Same as Fig. 2 but for an initially retrograde perturber.

Table 1 .
Yurin & Springel 2014cal parameters of the three galaxy components: halo, bulge and disc.For the DM halo, GalIC assumes an Hernquist profile with the scale radius set by the concentration parameter of the halo (assumed c = 10; seeYurin & Springel 2014, for full details).N part and ε correspond to the number of particles and the Plummer-equivalent gravitational softening employed for each component respectively.