Secular encounters

Analytic computation of the secular effects of encounters on a binary: features arising from second-order perturbation theory

Adrian S. Hamers and Johan Samsing
Institute for Advanced Study, School of Natural Sciences, Einstein Drive, Princeton, NJ 08540, USA
Department of Astrophysical Sciences, Princeton University, Peyton Hall, 4 Ivy Lane, Princeton, NJ 08544, USA
E-mail: hamers@ias.eduE-mail: jsamsing@gmail.com
Abstract

Binary-single interactions play a crucial role in the evolution of dense stellar systems such as globular clusters. In addition, they are believed to drive black hole (BH) binary mergers in these systems. A subset of binary-single interactions are secular encounters, for which the third body approaches the binary on a relatively wide orbit, and such that it is justified to average the equations of motion over the binary’s orbital phase. Previous works used first-order perturbation theory to compute the effects of such secular encounters on the binary. However, this approach can break down for highly eccentric binaries, which are important for BH binary mergers and gravitational wave sources. Here, we present an analytic computation using second-order perturbation techniques, valid to the quadrupole-order approximation. In our calculation, we take into account the instantaneous back-reaction of the binary to the third body, and compute corrections to previous first-order results. Using singly-averaged and direct 3-body integrations, we demonstrate the validity of our expressions. In particular, we show that the eccentricity change for highly eccentric binaries can reach a plateau, associated with a large inclination change, and can even reverse sign. These effects are not captured by previous first-order results. We provide a simple script to conveniently evaluate our analytic expressions, including routines for numerical integration and verification.

keywords:
gravitation – celestial mechanics – stars: kinematics and dynamics – globular clusters: general – stars: black holes
pubyear: 2019pagerange: Analytic computation of the secular effects of encounters on a binary: features arising from second-order perturbation theoryA.2

1 Introduction

There exists a large variety of astrophysical settings in which a binary system is perturbed by a passing object in a wider orbit. In the planetary context, stellar flybys can drive wide binaries to high eccentricities, leading to strong tidal interactions and/or stellar collisions (Kaib & Raymond, 2014), or destabilizing planetary systems (Kaib et al., 2013). In the Solar system, stellar encounters play an important role in disturbing the Oort cloud, thus transporting comets into the inner Solar system (Oort 1950, and, e.g, Heisler et al. 1987; Duncan et al. 1987; Dybczyński 2002; Serafin & Grothues 2002; Fouchard et al. 2011; Higuchi & Kokubo 2015). Such perturbations may also be responsible for triggering white dwarf pollution in extra-solar systems (e.g., Veras et al. 2014).

In dense stellar systems, flybys in which the binary remains bound after the encounter without exchanges can be considered a sub-type of more general binary-single interactions (e.g., Hut & Bahcall 1983; Hut 1983; Heggie & Sweatman 1991). Such interactions are key to driving the late-stage evolution of dense stellar systems such as globular clusters (e.g., Spitzer 1987; Binney & Tremaine 2008).

If the perturber passes the binary with a periapsis distance that is significantly larger than the binary semimajor axis, then the encounter is typically ‘secular’, i.e., the motion of the perturber is much slower than the orbital motion of the binary. In this case, it is appropriate to expand the Hamiltonian in terms of the ratio of the binary separation to the separation of the perturber to the binary center of mass, and to average the Hamiltonian over the binary motion. These procedures result in the simpler ‘singly-averaged’ (SA) equations of motion, which are computationally less intensive to solve compared to direct 3-body integrations. However, it is still necessary to numerically solve a set of first-order ordinary differential equations (ODEs) in order to obtain the new binary properties after the perturber’s passage.

It is also possible to apply first-order perturbation techniques, i.e., to analytically integrate the SA equations of motion assuming the perturber moves on a hyperbolic or parabolic trajectory, and ignoring changes of the binary’s orbital elements during the perturber’s passage. This approach was adopted in previous works (e.g., Heggie 1975; Heggie & Rasio 1996; Hamers 2018). In particular, the equations derived by Heggie & Rasio (1996) are commonly used to take into account the effects of (distant) encounters on a binary in Monte Carlo-style computations (e.g., Rasio & Heggie 1995; Spurzem et al. 2009; Geller et al. 2019).

However, there are situations in which the first-order perturbation approach can break down. In Heggie & Rasio (1996), it was discussed that their method no longer applies if the eccentricity change is of the order of the initial eccentricity, i.e., . However, as we will show here, when the binary is already highly eccentric, then the first-order approach can break down even if is small, in particular, if . The latter can occur in highly eccentricity binaries perturbed by distant (i.e., weak) encounters, which can drive relatively large changes in the binary’s orbital angular momentum. A binary can be driven to high eccentricities as a result of strong (non-secular) encounters in dense stellar systems. The breakdown of the first-order method originates from the fact that the binary’s elements can change significantly during the passage of the perturber. Therefore, it is not justified to assume that these elements are constant when integrating the equations of motion.

In this paper, we present an analytic computation of the change of the binary’s properties in the secular regime using second-order perturbation theory. In short, a Fourier expansion is used to describe the binary’s instantaneous response to the perturber, and this response is substituted back into the equations of motion in order to obtain correction terms to the first-order result. This approach is similar to that by Luo et al. (2016) applied to hierarchical triples (i.e., with a bound third object instead of an unbound one). In triples, similar interactions on the timescale of the third body can modulate the long-term secular interactions (e.g., Brown 1936; Ćuk & Burns 2004; Ivanov et al. 2005; Katz & Dong 2012; Antonini & Perets 2012; Seto 2013; Antognini et al. 2014; Bode & Wegg 2014; Lei et al. 2018; Grishin et al. 2018b). The correction terms give rise to behavior such as orbital flips and sign reversals in the eccentricity change that are not captured by first-order theory.

In an accompanying paper (Samsing et al., in prep.), we investigate the astrophysical implications of secular encounters on highly eccentric binaries, for which the first-order perturbation approach can break down. In particular, the properties of black hole (BH) binary mergers and their subsequent gravitational wave (GW) signals are discussed.

This paper is structured as follows. We give an overview of previous results in Section 2. In Section 3, we apply second-order perturbation techniques, and derive analytic expressions for the changes of the binary’s properties in response to secular encounters. These results, which are summarized (for the case of parabolic encounters) in equation (3.2.2), can be interpreted as additions to the work of Heggie & Rasio (1996). In Section 4, we provide a simple Python script to conveniently evaluate our analytic expressions, illustrate the regime in which the first-order approximation breaks down, and demonstrate the correctness of our second-order results using numerically-obtained SA solutions, and direct 3-body integrations. We discuss our results in Section 5, and conclude in Section 6.

2 Previous results

2.1 Preliminaries

Figure 1: Sketch of the configuration. A binary (center of mass ‘CM’) is perturbed by a passing object (mass ) on a parabolic or hyperbolic orbit (eccentricity ) with periapsis distance to the binary center of mass (, where is the binary’s semimajor axis).

Consider a binary with semimajor axis and eccentricity ; the component masses are denoted as and , with the total binary mass . The instantaneous eccentricity or Laplace-Runge-Lenz vector is given by , where is the relative separation between the components, dots denote derivatives with respect to time, and hats denote unit vectors. The normalized angular-momentum vector is , with magnitude . Apart from the orbital phase, the binary state is completely described in terms of and (note that the six degrees of freedom in these two vectors reduce to effectively four, since , and – together with the orbital phase, the binary’s state is described by five quantities).

We assume that the binary is perturbed by a third body with mass , passing by on a parabolic or hyperbolic orbit with periapsis distance , and separation relative to the binary center of mass. In the case of a parabolic (hyperbolic) orbit, the perturber orbit’s eccentricity is (). Without loss of generality, we set the angular-momentum vector of the perturber’s orbit to be aligned with the -axis, and the eccentricity vector along the -axis (i.e., the periapsis of the perturber’s orbit points along the -axis). A sketch of the configuration is shown in Fig. 1.

We also assume that the perturber moves sufficiently slowly such that the ‘secular’ approximation applies, i.e., the mean motion of the binary is much faster than the perturber’s angular speed at periapsis. This condition can be written as (e.g., Hamers 2018)

(1)

This condition typically implies that the binary is ‘hard’, i.e., its orbital energy is much larger than the typical kinetic energy of perturbers (Heggie, 1975). More quantitatively, the ‘hardness’ of a binary in a cluster with velocity dispersion can be defined as , i.e., the ratio between the (absolute value of the) specific binding energy, and the typical perturber kinetic energy. Writing , this implies that and are related according to

(2)

Therefore, if the binary is hard, , and

(3)

Unless and/or , this implies that , i.e., a secular encounter.

Some of the results below are presented in terms of the binary’s orbital elements , in addition to the orbital vectors and . We define the inclination , argument of periapsis , and longitude of the ascending node in the usual way, i.e., the relation between the orbital angles and orbital vectors is

(4a)
(4b)

We emphasize that our definition of the orbital elements differs from that of Heggie & Rasio (1996). The latter authors chose the binary orientation to be fixed, and used orbital elements to describe the orientation of the third body’s orbit. In contrast, we fix the orientation of the third body’s orbit, and use orbital vectors/elements to describe the binary’s orientation. In practice, this implies , , and , where elements with the subscript HR96 refer to Heggie & Rasio (1996), and elements without subscripts refer to the elements adopted here.

2.2 Equations of motion

We expand the Hamiltonian of the 3-body system in terms of the ratio of the separation of the binary, , to the separation of the third body relative to the binary center of mass, . To lowest order in , i.e., (quadrupole order), and after averaging over the binary mean motion assuming a Kepler orbit, the equations of motion for the binary orbital elements read (e.g., Hamers 2018)

(5a)
(5b)

Here, is the true anomaly of the perturber’s (parabolic or hyperbolic) orbit, and the small parameter

(6)

Generally, , where

(7)

Assuming that lies in the -plane, we can write whereas for the binary, and . Therefore, equations (5) can be written in the explicit form

(8a)
(8b)

We will refer to equations (5) and (8) as the SA equations of motion: averaged over the binary, but still explicitly a function of the perturber’s orbital phase .

2.3 First-order equations

The SA equations of motion constitute a system of ODEs and can be solved numerically. In addition, as in previous works (Heggie & Rasio 1996; Hamers 2018), the SA equations can be integrated over time from minus to plus infinity, or equivalently, over with , to obtain the changes in and . In the first-order (FO) approximation, the orbital vectors are assumed to be constant during the integration, i.e.,

(9a)
(9b)

In equations (9), the integrands are evaluated at the initial values of the eccentricity and angular-momentum vectors, and the latter are assumed to be constant throughout the perturber’s passage (the components of and in the explicit expressions in equations 9 should be interpreted as the components of the initial vectors; for brevity, we omitted the subscript ‘0’). In other words, in the FO approximation, the binary is not allowed to respond to the instantaneous perturbation of the third body.

In the limit when the change in the eccentricity vector is small, equation (9) yields the following more manageable expression for the scalar eccentricity change,

(10)

It can be shown that equation (10) is equivalent to the quadrupole-order result of Heggie & Rasio (1996; see also Appendix A2.2 of Hamers 2018).

3 Analytic second-order computation

The FO approximation works well in the limit when the changes in the binary’s state are small. However, in some cases, the change in and can be significant compared to the initial values (for example, when the initial eccentricity is already very high). In this case, it is no longer justified to assume that and are constant as in equations (9), i.e., FO perturbation theory is insufficient (see also Hamers 2018).

To derive improved expressions, we adopt a similar approach as Luo et al. (2016). In short, we develop Fourier expansions of the SA equations of motion. The orbital vectors, , are also written in the form of Fourier expansions. By plugging the Fourier expansions of and into the equations of motion, we obtain explicit relations for the Fourier coefficients associated with and . Subsequently, we integrate the equations of motion, now taking into account the leading order Fourier coefficients. The result consists of the original FO term, plus new additional terms to second order (SO). In other words, using Fourier expansions, we compute the instantaneous response of the binary to the perturber, and use the updated expressions to more accurately calculate the net effects on the binary elements.

3.1 Fourier expansions of the SA equations of motion

3.1.1 Calculation

First, we develop the Fourier series expansions of equations (8), i.e.,

(11a)
(11b)

with the Fourier coefficients given by

(12)
(13)

Here, the SA equations of motion in the integrands are to be evaluated at constant and (i.e., the initial values). Explicit expressions for these Fourier coefficients are given in Appendix A.1. Note that, trivially,

(14a)
(14b)

The summations in equations (11) run from to . In the case of a bound third body (Luo et al., 2016), the Fourier coefficients vanish identically for . It turns out that, in our case of an unbound third body, the Fourier coefficients do not vanish for . In principle, an infinite number of terms should be included. However, in practice, it suffices to restrict to a finite . In our explicit expressions, we set (see Section 3.1.2 below for the justification).

Next, we assume that the binary orbital vectors can be written in terms of similar Fourier series but with the addition of a linear term, i.e.,

(15a)
(15b)

Here, , , and are constant vectors, and the Fourier coefficients , , and associated with and are understood to be functions of the initial values of and . We substitute equations (15) into equations (11) by differentiating the former with respect to . Subsequently, we equate the coefficients of all terms (constant terms, and the sine and cosine terms) on both sides of the resulting equality. This procedure gives the following relations between the coefficients of the orbital vectors, and , and those of the equations of motion,

(16a)
(16b)

The necessity of the linear terms in equations (15) becomes clear by noting that without them, the -derivatives of equations (15) would not contain a term independent of , whereas such a term is present in equations (11). We remark that the linear term was omitted in Luo et al. (2016).

Equations (15) describe the instantaneous response of the binary to the perturber in terms of an infinite Fourier series. They can be written in the more compact form

(17a)
(17b)

Here, and are the initial state vectors, and and are (known) functions of , defined by

(18a)
(18b)

3.1.2 Testing the Fourier expansions

Figure 2: Eccentricity (top panels) and inclination (bottom panels) of a binary perturbed by a third body as a function of the third body’s orbital phase . Fixed parameters are , , , , , and . In the left (right)-hand panels, (), corresponding to (). Black solid lines show and by numerically integrating the SA equations of motion, i.e., equations (8); red lines correspond to the Fourier series in equations (17), where we take to be either 1, 2 or 3 and shown with the dotted, dashed and solid lines, respectively.

To demonstrate the correctness and usefulness of the Fourier expansions, we show in Fig. 2 the eccentricity and inclination of a binary perturbed by a third body as a function of the third body’s orbital phase . We set , , , , , and . In the left (right)-hand panels, (), corresponding to (). The black solid lines show and obtained by numerically integrating the SA equations of motion, i.e., equations (8). Within the limit of the SA and quadrupole-order approximation, these solutions can be considered to be exact. Red lines correspond to the Fourier series in equations (17), where we take to be either 1, 2 or 3.

In the case of , the Fourier series accurately describe and , provided that . There are no noticeable differences between and . For (i.e., larger ), the Fourier series still give a reasonable approximation of the true and , although there are significant differences, in particular for . The discrepancy in the case between the numerical solutions and the Fourier series can be traced to the fact that we included terms of order in equations (17), and ignored terms of order . Evidently, higher-order terms become increasingly important as increases. Note that convergence with respect to is distinct from convergence with respect to . In the case of , the series in still converges for .

3.2 Calculating corrections to second order in

3.2.1 General case: hyperbolic orbits

The next step is to substitute and up to and including first order in into the equations of motion, equations (11), where the Fourier coefficients , , , , , and are now all evaluated using equations (17), and subsequently integrating over from to .

To first order in , this procedure simply yields

(19a)
(19b)

i.e., we recover the FO result (note that, to zeroth order in , , and similarly for the sine terms and the other Fourier coefficients). For future convenience, we defined the two functions

(20)

To second order in , the general result is (much) more complicated. Specifically,

(21a)
(21b)

Here, and , which are composed of and respectively, are known functions of , , and . The functions and are easily obtained through

(22)

Here, and are the -independent terms of equations (17), i.e.,

(23)

In other words, and can be obtained by taking the FO result (i.e., and for the eccentricity and angular-momentum vector, respectively), and replacing the initial state vectors with and given by equations (23). Note that this result is independent of the value of . The explicit substitution relations in equations (23) are given in Appendix A.2. The functions and do depend on the assumed value of , and are generally very complicated. We therefore do not give the explicit expressions for and here. Instead, we implemented all required functions in an easy-to-use script that is freely available and described below (see Section 4).

In compact notation, the SO result can be written as

(24a)
(24b)

In the limit of small changes, the scalar eccentricity change to second order in is given by

(25)

Generally, the inclination change is given by

(26)

3.2.2 The parabolic limit

In the parabolic limit (), the SO expressions become much less complicated. Specifically,

(27a)
(27b)

The above expressions imply that, in the limit of small changes and parabolic orbits, the scalar eccentricity change is given by

(28)

In equation (3.2.2), we also wrote the result in terms of the orbital elements. Using that and , together with evaluated at (cf. equation 6), it is straightforward to show that equation (3.2.2) is consistent to order with equation (8) of Heggie & Rasio (1996). The novel result of this paper is the SO term in equation (3.2.2).

4 Exemplifying the correctness and relevance of the second-order expressions

In this section, we illustrate interesting behavior in the limit of high binary eccentricities, in which case the analytic FO expressions for the scalar eccentricity change can break down and a SO calculation is warranted. We demonstrate the correctness of our SO expressions by comparing to results obtained by numerically integrating the SA equations, and by directly integrating the 3-body equations of motion (i.e., without averaging over the binary’s motion).

We freely provide111https://github.com/hamers/flybys a simple Python script to compute the FO and SO expressions. In this script, which requires only the numpy and scipy libraries, we also include functionality to integrate the SA equations of motion numerically, and to carry out direct 3-body integration. The script can be used to quickly generate all the figures presented in this paper.

4.1 Phase evolution

Figure 3: Binary eccentricity (top panel) and inclination (bottom panel) as a function of the perturber’s true anomaly obtained by numerically integrating the SA equations of motion, assuming , , , , and . In the left (right)-hand panels, (. In each panel, three values of are chosen indicated in the legends. The red horizontal lines show according to FO perturbation theory (equation 10), with the line styles corresponding to the legends. The FO clearly breaks down in some cases

In Fig. 3, we show the binary eccentricity and inclination as a function of the perturber’s true anomaly assuming , , , , and . In the left (right)-hand panels, (. In each panel, three values of are chosen, ranging from relatively large, (), to small (). These results were obtained by numerically integrating the SA equations of motion. The red horizontal lines show according to FO perturbation theory (equation 10).

We first consider the case when the initial binary eccentricity is (left-hand panel in Fig. 3). For relatively large , , the changes in eccentricity and inclination are small. The FO equation accurately describes the scalar eccentricity change. However, as decreases, the agreement worsens, and for , the FO equation even yields , which is clearly not physical (note that, according to both SA and direct integration, the binary remains bound with a high eccentricity).

In the right-hand panel of Fig. 3, the initial binary eccentricity is even higher, (all other parameters are the same as in the left-hand panel). As shown with the SA integrations, as decreases, the maximum eccentricity reached during the perturber’s passage increases. This continues until, at a critical value of , the maximum eccentricity reached during the passage approaches unity, and the binary subsequently ‘bounces back’ to lower eccentricity. Consequently, the net scalar eccentricity change becomes negative. This phenomenon is also associated with a large change in the inclination, and is similar to the flipping of orbits associated with very high eccentricities in hierarchical systems such as triples (e.g., Lithwick & Naoz 2011), and higher-order systems (e.g., Hamers & Lai 2017; Grishin et al. 2018a). The FO equation fails to describe this ‘bounce’ effect, and incorrectly predicts a net increase of the eccentricity, and which is .

4.2 Series of

Figure 4: Scalar binary eccentricity (top panels) and inclination (bottom panels) change as a function of . Fixed parameters are , , , , and . In the left (right)-hand panels, the binary’s inclination is (). Results are shown from different techniques: analytic FO (black dashed lines), analytic SO (colored solid lines), numerical SA (colored open circles), and 3-body integration (colored stars). Blue (red) colors correspond to positive (negative) changes. The black horizontal dotted lines show , the maximum allowed (positive) change. The green dotted vertical line shows corresponding to (see equation 1); encounters with less than approximately this value are within the non-secular regime. The two vertical black dotted lines show the values of for which we can analytically calculated the locations of the sign change and plateau in using the SO expressions (see Section 4.2.2).

4.2.1 General behavior

As described above, for high initial binary eccentricities and assuming that for weak perturbations the binary’s eccentricity is increased, for sufficiently strong perturbations (small ), the binary eccentricity shows a ‘bounce’ phenomenon, resulting in a net negative scalar eccentricity change. Here, we describe this effect more quantitatively and show that our SO results can be used to accurately describe the binary’s eccentricity and angular-momentum changes in the associated parameter space.

In Fig. 4, we show the scalar eccentricity change and inclination in the top and bottom panels, respectively, as a function of . The fixed parameters are , , , , and . In the left (right)-hand panels, the binary’s inclination is (). We show analytic FO results (black dashed lines), analytic SO results (colored solid lines), numerical SA results (colored open circles), and 3-body integration results (colored stars). Blue (red) colors correspond to positive (negative) changes.

The 3-body integrations were carried out by formulating the 3-body equations of motion into ODE form and solving this set of equations using the odeint routine from the Python package scipy. Typical relative energy errors were (for small , corresponding to few binary orbits), to (for large , corresponding to many binary orbits). These energy errors are somewhat large, but sufficient for our purposes.

We first describe the behavior shown in Fig. 4 according to the SA and 3-body integrations. For large , the SA and 3-body integrations agree well. As decreases, starts to flatten, and reaches a local maximum, at . This plateau in is associated with a large inclination change (see the bottom panels in Fig. 4). In the case when , the plateau corresponds to , the maximum allowed positive change (shown with the black horizontal dotted lines).

For smaller , decreases with increasing , and at a critical value of , , changes sign. Subsequently, increases, with a different slope compared to before the sign change (in fact, this slope is , as shown below in Section 4.2.2). At small , , the SA and 3-body results start to differ significantly. This is the result of the breakdown of the SA approximation: for very small ( approaching unity), the binary’s motion is no longer much faster than the perturber’s motion, therefore it is no longer justified to average over the binary motion. This is also supported by the green vertical dotted line, which shows the value of for which (see equation 1). The scatter in the 3-body results originates from the binary phase-dependence in this regime (we assumed a fixed initial binary phase of ).

The FO expression gives a good description at large , but clearly fails as . The SO results, shown with the solid colored lines, agree well with the SA results for any and, therefore, with the 3-body results, unless .

Figure 5: Scalar eccentricity changes as a function of for different values of and (indicated in the top of each panel). Refer to the caption of Fig. 4 for a description of the meaning of the various lines and symbols.

In Fig. 5, we show similar results for different parameter choices in four panels with (parabolic encounter) or , and or . The SO expression generally gives a good description. Some deviation can be seen in the case when and .

4.2.2 Useful properties of the SO solution

Evidently, the SO solution removes the need for numerical integration of the equations of motion. The latter can be achieved with the SA approximation, or with the most general but also most computationally expensive method of direct 3-body integration. The implied speed-up can be advantageous when a large number of encounters need to be evaluated, e.g., in Monte Carlo approaches.

In addition, the SO equations can be used to gain more insight into the effect of secular encounters. In the case of the series of discussed in the previous section, the value of (corresponding to a particular ) for the sign change of can be simply expressed using equation (25), i.e.,