© L. Rodet et al. 2019

1 Introduction

In the context of the rapid increase of exoplanet discoveries, the need for efficient N-body simulations has become strong to model the evolution of complex systems and the interaction between planets, planets and debris disk, or within debris disks. Mixed variable symplectic integrators are widely used for dynamical simulations of planetary systems, as they present two major advantages with respect to other N-body integrators: first, they exhibit no long-term accumulation of energy error, which is essential to ensure orbital stability through the integration. On the other hand, they provide a gain of at least one order of magnitude in computation speed, for equivalent accuracy, because they allow one to adopt a much larger time-step than other integrators for the same result. In 1991, Wisdom and Holman devise the first symplectic map specifically designed for N-body problems with a central dominant mass (Wisdom & Holman 1991). Since then, numerous codes implemented this structure that are still widely used today (e.g., Swift, Levison & Duncan 1994, Mercury, Chambers 1999).

Yet, symplectic integrators can model the interactions between multiple stars, moon, or simply planets whose mass are non negligible with respect to the central mass as well. They are versatile tools well suited to characterize the great diversity of extrasolar system architectures, well beyond the framework of our Solar System. Efforts were made to extend the scheme to binary stars in two modified versions of Mercury (Chambers et al. 2002), but it could not be generalized to multiple systems with other hierarchies. In this context, Beust (2003) designed a symplectic scheme valid for any type of hierarchical architecture, and implemented it with Swift HJS. This generalized the theoretical frame of Wisdom and Holman to any hierarchical system.

However, in Swift HJS, the hierarchical structure of the system is given at the beginning of the run and must be preserved along the integration. This is a severe limitation as it prevents the efficient modeling of non stable hierarchies with e.g. orbital captures (planets, dust), whereas such situations may be numerous among young systems. With Swift HJS, handling accurately such configurations is only possible adopting a very small time-step, which is of course not optimal. This motivated us to build a new version of Swift HJS, ODEA, that tackles this issue. The code is available online¹.

In the following, we describe the new code in detail, and present a full application to the complex system of HD 106906. Before that, we present this system and our motivations for modeling it and using it as a benchmark for our new code.

The system HD 106906 (HIP 59960) is located at a distance of 103.3 ± 0.5 pc (Brown et al. 2018) and belongs to the Lower Centaurux Crux (LCC) group, which is a subgroup of the Scorpius-Centaurus (Sco-Cen) OB association (De Zeeuw et al. 1999). The LCC group has a mean age of 15 ± 3 Myr, with an age spread of 6 Myr (Pecaut & Mamajek 2016). HD 106906 is a 2.58 ± 0.04 M_⊙ spectroscopic binary star, on an eccentric (0.66) and tight (0.6 au) orbit (Lagrange et al. 2019). Moreover, high contrast imaging has revealed an asymmetric debris disk (Kalas et al. 2015; Lagrange et al. 2016) and a giant planet on a wide orbit (projected separation from the binary: 735 ± 5 au, Bailey et al. 2013). At such a separation, the planet relative motion can not be detected with present imaging instruments on a reasonable time basis. The orbital inclination with respect to the plane of the disk is probably significant (20°), but a coplanar configuration cannot be excluded. The planet mass has been estimated at 11 ± 9 M_J mass from hot-start models by Daemgen et al. (2017).

Two major scenarios compete for the formation of giant planets (e.g., Baruteau et al. 2016). In the core accretion scenario, planets begin their formation with the growth of dust grains and the formation of planetesimals, that will slowly accrete each other to form terrestrial planets or planetary cores. On the other hand, the gravitational instability scenario is a faster process that is able to form giant planets at large separation from an instability in the protoplanetary disk. In both cases, planet formation takes place in the primordial gaseous disk. Forming a giant planet at 700 au ormore from any central star appears very unlikely in any of those scenarios, first due to the lack of circumstellar gas at that distance, and second because the corresponding formation timescale would exceed the lifetime of the gaseous disk. This led Rodet et al. (2017) to propose a dynamical scenario to account for the planet’s current separation. The scenario involves a traditional planetary formation within the gaseous disk, an inward migration and a subsequent scattering by the binary. However, for the planet to remain bound, an external perturbation such as a fly-by is necessary in order to reduce its eccentricity and stabilize its orbit in a bound configuration.

Recently, De Rosa & Kalas (2019) investigated the stellar neighborhood of system HD 106906 in Gaia DR2 (Brown et al. 2018), and discovered two stars that have recently come within 1 pc of the central binary HD 106906 AB. Given the uncertainty on the perturbers distances and radial velocities, De Rosa & Kalas concluded that there was a possibility that the fly-by was dynamically significant for the planet evolution history. This motivates us to reinvestigate the Rodet et al. (2017) scenario, using ODEA, to check this possibility.

2 Algorithm

2.1 Structure of the code: Swift HJS

Let us consider the gravitational N-body problem, with masses , positions and impulsions . The Hamiltonian is (1)

where G is the constant of gravitation and r_ij = ∥r_j −r_i∥ is the distance between bodies i and j.

In the current version of Swift HJS, as in the other similar codes, the integrator do not solve H exactly, but a surrogate Hamiltonian. The latter ischosen to be close to the real one, and exactly solvable. In that case, the algorithm is symplectic: it exactly preserves the areas in phase space and exhibit no long-term drift of the energy.

In order to design a proper in orbital mechanics, the key idea is to split the Hamiltonian into two integrable parts: (2)

Several splitting have been suggested (e.g., Wisdom & Holman 1991; Saha & Tremaine 1994; Chambers 1999), most of them consisting on a Keplerian part and a perturbation part. Both parts are then integrable within computer round-off errors. corresponds to the successive integration of these parts separately. For a second order symplectic integrator, a so-called leap-frog method can be used. It consists in integrating H_B for Δt∕2 (kick), then H_A for Δt (drift), then again H_B for Δt∕2 (kick), where Δt is the time step.

Swift HJS is based on the Hierarchical Jacobi Symplectic method introduced by Beust (2003), where the description is based on orbits instead of on bodies. An orbit consists in a collection of two non-empty sets of bodies, the set of centers and the set of satellites, that have empty intersection. In all problems in orbital mechanics, a hierarchy can then be defined as a collection of orbits comprising all bodies satisfying the following rule: for all couples of orbit k and l≠ k, one of the three subsequent propositions apply

• orbits k and l have no common bodies (orbits k and l are foreign);

• orbit k is comprised in the centers or satellites of orbit l (orbit k is inner to orbit l);

• orbit l is comprised in the centers or satellites of orbit k (orbit k is outer to orbit l).

A so-defined hierarchy is made of exactly N − 1 orbits. In Swift HJS, the orbits are numbered from 2 to N. Finally, we define μ_k and η_k as the total mass of the satellites and centers, respectively, in orbit k. The total dynamical mass in orbit k is then M_k = μ_k + η_k and the reduced mass .

In this formalism, a new set of N coordinates are designed with a Jacobi-like approach: is the relative position of the center of mass of orbit k’s satellites with respect to that of its centers, and is the relative conjugate momentum. The first coordinates and are the position and impulsion of the center of mass. These positions and conjugate momenta derive from a canonical transformation that let the Hamiltonian invariant. They can be expressed with the bodies coordinates as

The Hamiltonian can then be split as follows

When the hierarchy is sufficiently clear (that is if the orbits are almost Keplerian), H_B ≪ H_A. As H_A is a Keplerian Hamiltoniandescribing N − 1 independent orbits, the drift consists of evolving each Keplerian orbits. On the other hand, as H_B depends exclusively on the positions, the kick consists of a linear raise of the velocities, with accelerations .

Fig. 1 Example of hierarchy change in the case of a capture. At first the red body orbits the yellow–blue pair. After a strong interaction, it captures the small blue body.

2.2 Building a new hierarchy

The above scheme is well adapted to lightly perturbed Keplerian orbits in a fixed hierarchy, but becomes strongly unsuitable if the initial hierarchy evolves, whether temporarily or definitively (see example Fig. 1).

Thus, when the hierarchy is not relevant anymore (that is the splitting in the initial H_A and H_B does not optimize the error), a module of the algorithm will design a new hierarchy from the current positions of the bodies. For this, the algorithm computes a two-dimensional symmetric array that compiles the Keplerian acceleration between two bodies , where M_k is the sum of the masses. The strongest acceleration gives the first orbit, then the two bodies are replaced by their center of mass and the array is updated, and again until the last orbit comprises all bodies. We first checked that this algorithm always returns the existing hierarchy when no change is expected. Then, if the computed hierarchy is different than the current one, the hierarchy must be changed.

If the hierarchy needs to be changed, so is the time-step Δt. We choose a Keplerian-like time , where (7)

if orbit k is bound or if its smallest approach has not yet occurred, or (8)

otherwise. The choice to adapt or not the time step is given to the user.

Strictly speaking, when changing the hierarchy, the symplectic nature of the algorithm does not hold anymore, as the splitting of the Hamiltonian is entirely based on the hierarchy. This is also true for any change of the time step. A new approximate Hamitonian is integrated from an already approximated scheme, which means that the error budget raises potentially at each hierarchy change. However, the algorithm is designed for orbital dynamics, where systems are not subject to frequent reorganization of their architecture. Designing a new Hamiltionan when the initial hierarchy is not suited anymore allows to limitate the error on each orbit, which will otherwise become out of control. This is basically the same problem as the one raised by close encouters in planetary dynamics. When handling close encounters (Levison & Duncan 1994; in Swift RMVS) and (Chambers 1999; in Mercury) temporarily change the way of splitting the Hamiltonian when transferring to H_A the part of H_B that concerns the close encounter, even sometimes changing the hierarchy to planetocentric (in the latest version of Mercury the use of a smooth criterion that weights the different perturbing terms allows the map to remain symplectic with a continuous Hamiltonian while handling close encounters; Rein et al. 2019). Conceptually, a close encounter within a planetary system can be viewed as a temporary change of hierarchy that eventually returns to the initial hierarchy. Here we are concerned by changes that can be permanent.

2.3 Checking the relevance of the hierarchy

Performing a hierarchy change is quite costly: all the acceleration couples have to be computed and must be compared and updated for the definition of each of the N − 1 orbits (multiple operations that scale as O(N³)). Checking for a possible change at each time-step, with the result that most of the time the current hierarchy would be left unchanged, would thus amount to a considerable loss of efficiency. Prior to launching the entire hierarchy re-building process, an efficient algorithm with a simpler criterion must be applied to check whether it is appropriate or not. The most exact criterion would be the theoretical energy error associated to the symplectic mapping, as it gives us an objectiveestimate of the relevance of the numerical scheme. However, its computation is tedious (grows as N⁴, see Appendix B). The criterion must be fast to compute (maximum as N³, like the accelerations) and correlated to the error.

In Mercury (Chambers 1999), the criterion to spot close encounters is the ratio between the relative distances and the Hill radii, assuming the latter roughly constant. This is a legit criterion for the study of the Solar system, but it is not relevant to our case. Indeed, the Hill radius is not easy to compute for eccentric orbit, it depends strongly on the orbital parameters (which is subject to variation in the general case) and it is not satisfyingly correlated to the errors in a complex architecture.

We choose to compute at each step the ratio for each orbit k, where and are the accelerations respectively induced by H_A and H_B (Eqs. (5) and (6)). We declare the hierarchy questionable if it is higher than 0.2 for at least one orbit. The computation of that criterion also scales as O(N³) in theory, but it uses the acceleration a^B that is already computed in any step of the integration, so that the extra cost remains limited. Thus, with this criterion, the problem can keep a non-optimal hierarchy if the associated error remains small. This can be adjusted by changing the value of the threshold, which is a free parameter of the code. This might be useful in situations when two similar hierarchies become alternately optimal, to prevent the algorithm to perform numerous changes that will have a negative effect on long term conservation properties.

2.4 The case of test particles

The study of planetary systems often involved the study of debris belts. In N-body simulations, the dust is modeled at first order by massless bodies (or test particless) that interact with the massive bodies but not with each other. Test particles must be specifically considered in ODEA as the handling of their hierarchy is slightly different. Indeed, they are the only satellites of their orbit and their orbit is invisible to the bodies and other test particles evolution. When looking for a new hierarchy, ODEA will not considerthe test particles, for it searches foremost to optimize the energy error budget related to the massive bodies.

When the hierarchy of the massive bodies changes, each test particle must find its natural orbit given its relative position. A similar procedure to the hierarchy building of massive bodies is then performed. For a consistent hierarchy, the test particles have 2N − 1 possibilities for their orbit: around one massive bodies (N) or around one orbit (N − 1). Thus, for each test particle, a 2N − 1 array is computed, compiling the Keplerian accelerations. The maximal element will correspond to the new particle configuration.

Finally, a test particle may also be subject to a hierarchy change, independently of the massive bodies architecture evolution. Thus, the acceleration ratio criterion is computed at each time step to check the suitability of the particle orbit, and a new orbital configuration is investigated if necessary following the previous procedure.

Fig. 2 Maximum relative energy error over a 1 Myr evolution of HD 106906 including the two fly-bys, for ODEA and Swift HJS, as a function of the time-step assumed.

2.5 Comparison with other codes

Several algorithms have been introduced since the formalization of the first mixed-variable symplectic map for orbital mechanics, including the widely used Mercury (Chambers 1999). Most of them are designed to work in Solar-System-like hierarchy. Chambers et al. (2002) introduced two algorithms, derived from Mercury, to model planetary motions in binary systems. However, to our knowledge, no symplectic integrator are able to integrate indifferently any types of hierarchy, or a more complex hierarchy, except from Swift HJS.

Moreover, no mixed-variable integrator that we know of are designed to handle long or definitive hierarchy change. Such situations can be encountered in case of a stellar fly-by, or of a capture of debris disk dust by a stellar or planetary companion. The subsequent study of system HD 106906 is a perfect example of situations that can not be tackled by ordinary symplectic algorithms: binary fly-by and dust capture. Figure 2 illustrates the gain of energy precision that ODEA allows on the HD 106906’s fly-bys test case (the parameters of the corresponding simulation are presented in Appendix A). The relative energy error is here entirely dominated by the close encounters. By changing the hierarchy, ODEA decreases the error of one or two orders of magnitude compared to Swift HJS. Moreover, it allows the energy error to decrease after the encounter in case of definitive hierarchy change.

On an other hand, Rein & Spiegel (2014) argue that a high-order classical integrator is quicker and more accurate than symplectic integrators. This may be true for some complex cases, or if we aim for a very high precision. However, symplectic integrators have encoded the exact resolution of the Keplerian motion, while a classical integrator makes no hypothesis for the form of the motion, and has to solve from scratch the differential equations of motion. Thus, for lightly perturbed Keplerian motion, symplectic algorithms are certainly more practical than classical integrators. The time steps can be large without endangering the stability of the orbits.

For example, in the case of HD 106906, the simulations involved very different scales, from the planet periastron to the wide hyperbolic orbit of the perturbers. A classical integrator would have to adapt its time step to the smallest distance, while a symplectic integrator can adopt a larger timescale without compromising the stability of the planet orbit. This is illustrated in Table 1, where the two symplectic integrators ODEA and Swift HJS can achieve a reasonable precision with a large time step (same simulation than for Fig. 2, but for the entire 15 Myr evolution). Outside the fly-bys, they reach a precision similar to that of the classical integrator IAS 15, that has to decrease regularly its time step toresolve the periastron passage, increasing the computation time. We also ran the simulation with the Bulirsch-Stoer implementation of Press et al. (1989), with a precision constraint on the trajectory of order 10⁻⁷ (similar to the value reached by ODEA). The energy error grows very rapidly, and the computation time is already significantly larger than the other codes.

We also point out that Swift HJS never makes the assumption that the orbits we are considering are actually bound. The only requirement is that the sum of the Keplerian interactions associated with the hierarchy (i.e. H_A) must represent most of the full Hamiltonian. Some of the orbits we are considering can thus be hyperbolic, and this will be the case in a fly-by configuration. The Kepler solver used to integrate H_A handles bound or unbound orbits as well.

3 Application to system HD 106906

3.1 Characterizing the perturbers

Searching for potential stellar perturbers in Sco-Cen during the previous 15 Myr, De Rosa & Kalas (2019) identified two perturbers in LCC (Pecaut et al. 2012): HIP 59716 and HIP 59721. Located around 11 pc (projected 0.5°) from HD 106906 and 0.5 pc (projected 30″) from each other, their relative velocities suggest an encounter with HD 106906 a few million years ago. The coordinates and velocities of the three systems are summarized in Table 1 of De Rosa & Kalas (2019). As can be seen on Fig. 3, the relative separation and velocity between HD 106906 and its perturbers lie essentially on the direction to Earth. Unfortunately, the quantities projected in this direction (distance and radial velocity) have the larger observational uncertainties, which creates a high dispersion on the closest encounters, in particular for the most promising candidate HIP 59716 (Fig. 4).

We note that the relative velocities between each systems (~4 km s⁻¹) are four times higher than the velocity dispersion reported for LCC (1.13 ± 0.07 km s⁻¹; Madsen et al. 2002), that was used in Rodet et al. (2017). We will see in Sect. 3.3 that the effect of a fly-by is inversely proportionalto the velocity of the passing star.

The masses of HIP 59716 and HIP 59721 have been estimated, respectively, 1.37 M_⊙ for HIP 59716 and 1.22 M_⊙ for HIP 59721 from the spectral types. HD 106906 binary mass has been estimated to 2.58 ± 0.04 M_⊙ from radial velocity and interferometric measurements by Lagrange et al. (2019).

Fig. 3 Representation of HD 106906, HIP 59716 and HIP 59721 current positions and velocities in HD 106906 rest frame (disk lies in the Y –Z plane, observed extension in the Y -direction).

3.2 Simulating the encounters

N-body simulations performed by De Rosa & Kalas (2019) indicate that the galactic gravitational potential has a negligible influence on the characteristics of the encounters. Moreover, the binarity of HD 106906 does not affect the encounters, because of the very high ratio between the closest approaches and the binary separation (>1000). In order to efficiently determine the parameters of the encounters, we first performed 10 000 simulations with ODEA, including three bodies: HD 106906 ABb (2.58 + 0.01 M_⊙), HIP 59716 and HIP 59721. The mass of HD 106906 and the algorithm that we present here are the only differences with De Rosa & Kalas (2019) study at that point.

The initialization of the simulations is designed with a Monte-Carlo approach, following De Rosa & Kalas (2019). The 3 × 6 parameters and their respective precision arethe right ascension α (0.05 mas), the declination δ (0.002 mas), the parallax π (0.05 mas), the proper motion of the right ascension μ_α cosδ (0.05 mas yr⁻¹), the proper motion of the declination μ_δ (0.05 mas yr⁻¹) and the radial velocity γ (up to 1.7 km s⁻¹). The parameters are drawn from a normal distribution centered on their measured values, with a dispersion equal to the observations uncertainties, taking into account the correlations given by Gaia catalog. Then, we trace back the stars trajectory to observe the encounters.

Most of the simulations follow the same hierarchy evolution, represented on Fig. 5: the first fly-by involves HIP 59716 and the second HIP 59721, before the two perturbers get very close at each other as can be seen today. The hierarchy will thus naturally evolves to take into account the successive encounters. Computing the eccentricity of several sets of configurations, we evaluated that the two perturbers have currently a 2.1 ± 0.1% chance of being gravitationally bound to each other. However, De Rosa & Kalas (2019) point out that the probability of them having such similar angular positions and proper motions without being bound are extremely low.

We launched 10 000 simulations for 15 Myr, corresponding to a backward evolution from our days to the formations of the stars. The different timescales of the simulations are summarized in Table 2. At first sight, 10 000 simulations may not seem enough to correctly sample the 18 parameters confidence intervals. However, most of the parameters are strongly constrained, the only strong uncertainties being the perturbers relative radial velocities and distances, that is 4 parameters. Thus, these are the critical parameters that must be correctly sampled, and 10 000 is then a sufficient number. The initial time-step was set to 1000 yr, with outputs every 1000 yr. To account for the possibility of the two perturbers being bound, we performed an additional 10 000 simulations with only bound configurations. It comes down essentially to selecting only the configurations where the perturbers have similar radial velocities.

The distances at closest approach were computed for each simulation (Fig. 6). Most of the encounters occur with a closest approach between 0.3 and 2 pc, with a maximal probability around 0.6 pc, consistent with the results of De Rosa & Kalas (2019). We then reviewed the simulations for which a close (<0.1 pc) fly-by occurred, from any one or both of the two perturbers. 359 configurations were selected, that is around 4% of the total number of studied configurations. In most cases (≳90%), HIP 59716 encounters HIP 106906 at the shortest distance. For the bound configurations, the peak is around 0.4 pc but the number of close fly-bys is roughly the same. HIP 59716 coordinates distributions are presented on Fig. 8. Most of the parameters of the configurations with close fly-bys are drawn randomly within the configurations, except for the radial velocity, where we see that the configurations leading to a close fly-by correspond to the higher radial velocities (closer to the radial velocity of HIP 59721). The distributions for the two other bodies are presented on Figs. E.1 and E.2.

The distributions of the time and velocities of the perturber at closest approach are represented in Fig. 7 (only the cases where the distance was less than 0.1 pc). Most of the encounters occur between 4 and 2 Myr ago, with a velocity between 2 and 6 km s⁻¹.

Fig. 4 Two dimension histograms of the coordinates of the intersection points between the perturbers trajectories and the X–Y plane, assuming linear trajectories.

Fig. 5 Representation of a typical evolution of the hierarchy in the three-body simulations of HD 106906 fly-bys with ODEA. All orbits here are hyperbolic.

Fig. 6 Distribution of the distances at closest approach. The following study will focus on the red part, that corresponds to fly-by closer than 0.1 pc (3.6% of the configurations).

Fig. 7 Distribution of the times and velocities at closest approach, for the cases where the distance at closest approach is less than 0.1 pc.

3.3 Effect on the planet

3.3.1 Setup

Once the configurations for which a close fly-by occur within the 15 Myr of the system life have been identified, we launch a new set of simulations, this time including the planet. The bodies are initialized at their position at the end of the first simulation, that is at their position 15 Myr ago. HD 106906 is separated into two bodies, namely the binary HD 106906 AB (2.58 M_⊙), and the planet HD 106906 ABb (0.01 M_⊙). The simulations are launched from 15 Myr ago to the present epoch, so that the final outcome represents the current positions of the bodies. The time-step was set to 100 yr, with outputs every 1000 yr.

In the study of Rodet et al. (2017), the destabilization of the planet takes place after a violent encounter with the central binary, in the beginning of the system’s life. The outcome was either a definitive ejection on a hyperbolic trajectory, or a transitional state where the eccentricity raised dramatically without passing 1. The probability of the different outcomes depends on the characteristics of the encounter, which is highly underconstrained. In the case of a hyperbolic trajectory, a subsequent stabilization by a fly-by must be precisely synchronized, and is thus difficult to achieve. Thus, we study here the case of a highly eccentric transitional bound orbit. The periastron should roughly correspond to the separation of the planet when the perturbation occurred, around 1 au. On the other hand, the apoastron will remain mostly unchanged after a fly-by. The current projected separation implies a minimal value of 730 au. Moreover, the probability is higher to observe the planet near apoastron: it spends 2/3 of its time at a separation greater than 700 au for an apoastron of 1000 au, and 95% for an apoastron of 3000 au. All in all, two sets of simulations are performed, where the planet is initialized with a periastron of 1 au and an apoastron of 1000 (a = 500.5 au, e = 0.998) or 3000 au (a = 1500.5 au, e = 0.9993).

The necessary energy to completely eject the planet is , where a_p is the initial semi-major axis of the planet and M_HD106906 the mass of the host binary. From its current position close to the central binary, a definitive ejection requires around 1 M_⊙au² yr⁻². A proportion of 2 × 10⁻³ less corresponds to an elliptic trajectory with apoastron 1000 au, and 2 × 10⁻⁴ M_⊙au² yr⁻² less corresponds to 10 000 au. Thus, from an energetic point of view, reaching a high apoastron on a still bound orbit in the ejection process is nearly as costly as being definitely ejected.

For a fly-by to have a meaningful role in the dynamical history of the planet, it has to decrease the planet eccentricity by increasing the periastron to a safer value (an increase of the order of the astronomical unit at least). The time-scale of the fly-by is much larger than the orbital period of the planet, so that the initial position of the planet on its orbit is not a relevant parameter in the simulations. Moreover, in our scenario, the planet formed within the disk, so that its orbit was initially coplanar with the disk mid-plane. We assume that the planet apoastron is aligned with the observed extension of the disk. A close encounter with the central binary will retain this coplanarity if the inclination of the binary orbit is similar to that of the disk plane, which seems likely from the first estimates of its orbital parameters (Lagrange et al. 2019). As the fly-by is likely to keep the apoastron roughly unchanged and the eccentricity high (consistent with the observed patterns of the disk according to Jílková & Zwart 2015; Nesvold et al. 2017; Rodet et al. 2017), this is consistent with the current position of the planet.

Fig. 8 Initial distribution (today) of HIP 59716 coordinates and velocities for the 10 000 simulated cases (green), and for the 359 cases where a fly-by closer than 0.1 pc occurred (red).

3.3.2 Results

The conclusion of the study depends essentially on the possibility for the fly-by to increase significantly the periastron. This effect is stongly correlated to the distance at closest approach. We thus represented the periastron change with respect to the distance at closest approach for the outputs of the two sets of simulations on Figs. 9 and 10.

Whether for a 1000 or 3000 au apoastron, a 0.1 pc encounter is not enough to significantly raise the periastron: a closer fly-by is required. For the 1000 au apoastron case, the distance at closest approach must be less than 0.01 pc, that is 2000 au. For the 3000 au apoastron case, the destabilization is certainly easier, but the distance at closest approach must still be less than 0.05 pc, that is 10 000 au. For suchdistances, the results are essentially identical for the bound cases, as the separation between the two perturbers is greater of similar than the distance at closest approach with HD 106906. On our initial 10 000 draws, respectively, 2 and 20 resulted in a periastron increase superior to 1 au for the 1000 and 3000 au apoastron cases, and 1 and 2 lead to the ejection of the planet (for distance at closest approach similar or less than the planet semi-major axis).

Moreover, coplanarity of the planet orbit with the disk plan is expected if the planet formed within the disk. The current projected planet misalignment with the disk plane is currently estimated at 23°, although a lower angle (and even coplanarity) would be possible if the planet true separation is greater than its projected separation (≳ 3000 au for coplanarity). A 23° misalignment corresponds to a minimalaltitude of ~280 au above the disk plane, and such gain of altitude is rarely seen in the simulations, even in the most favorable case of a high initial apoastron. This would suggest that the misalignment (or part of it at least) is an illusion due to projection effects.

Fig. 9 Periastron increase with respect to the distance at closest approach, from N-body simulations (dots) and theoretical approaches (lines), for the closer fly-bys, and for an initial planetary apoastron of 1000 au. The grey part corresponds to a periastron change inferior to + 1 au, which will not secure the planet stability.

Fig. 10 Periastron increase with respect to the distance at closest approach, from N-body simulations (dots) and theoretical approaches (lines), for the closer fly-bys, and for an initial planetary apoastron of 3000 au. The grey part corresponds to a periastron change inferior to + 1 au, which will not secure the planet stability.

3.3.3 Theory

We first study the periastron increase as a function of the distance at closest approach, and compare it to the theoretical predictions. The computation of the following theoretical formula is explained in the appendix. The simplest approach is the impulse approximation, where the fly-by is assumed to be instantaneous and trigger a sudden velocity change on the planet. Although this cannot be considered as representative for the reality if we compare the fly-by time-scale with the orbital period of the planet, this approximation often provides a good estimate. In this framework, Brunini & Fernandez (1996) show that the fly-by increases the planet velocity by: (9)

where v_p is the planet velocity, M_* is the perturber’s mass, V its velocity at closest approach, D its distance at closest approach, and a_p the planet semi-major axis. This formula nevertheless applies to circular orbits only (Brunini & Fernandez 1996). By supposing that the new orbit intersects the old one at apoastron, the planet eccentricity e_p takes part, and we have a change of semi-major axis Δa_p = −a_pΔe_p, which gives a change of periastron Δperi = −2a_pΔe_p. Finally, one gets (see Appendix C): (10)

It can be adapted to an eccentric orbit, as was done in Rodet et al. (2017), by supposing that the perturbations occur only at apoastron. Then, stating that the apoastron is preserved, one gets Δa_p = −a_pΔe_p∕(1 + e_p) and Δperi = −2a_pΔe_p∕(1 + e_p). Finally, using Eq. (9) to quantify the velocity increase at apoastron, one gets (see Appendix C): (11)

On the other hand, a more rigorous approach is to compute the secular evolution of the orbital elements of the planet during the passage of the perturber. Heggie & Rasio (1996) used that method to determine the eccentricity increase of a companion, and found a complex formula depending on all 6 orbital elements of the perturber’s orbit. In this framework, the semi-major axis is invariant throughout the fly-by. Considering a coplanar orbit and a perturber’s eccentricity significantly higher than 1 (strongly unbound orbit), the maximum is: (12)

The three theoretical predictions are represented on Figs. 9 and 10: circular impulse, apoastron impulse and secular approximation. They all correspond to maximum values, as the true periastron evolution depends on the angular characteristics of the encounter. The velocity V is set to its mean value over all closest approaches, around 4 km s⁻¹. M_* was set to 1.3 M_⊙, but the increase depends weakly on the perturber’s exact mass. The eccentricity e_p is set to its initial value, an approximation that becomes less relevant when Δe_p ≳ 1 − e_p = 2 × 10⁻³ (for closest approach less or around 0.01 pc).

We see on Fig. 9 that the periastron change is best modeled by the secular approximation, but is also correctly approached by the impulse approximation at apoastron. It suggests that the effect of both perturbers on the planet can be estimated by the effect of the perturber that had the closest approach. This is also true for the cases where the two perturbers are bound (see Appendix E).

Furthermore, we seek to estimate if the fly-by could account for the possible misalignment of the planet with the debris disk plane. Depending on the exact value of the argument of periastron ω_p, a very eccentric orbit does not necessarily have a large elevation above the disk plane, even if it is highly inclined. To have a meaningful plan misalignment, the planet should have an inclination change combined with a shift of the argument of its periastron that results in a significant elevation above the disk plane. For any Keplerian orbit, the maximum elevation z_max above the reference plane is given by: (13)

Obviously, with e_p ~ 1 and ω_p ~ 0 or π, z_max remains small irrespective of the value of i_p.

We thus computed the change in z_max, inspiring from Heggie & Rasio (1996). The details are explained in Appendix D. The resulting maximal altitude is represented in Figs. 11 and 12.

Fig. 11 Maximal altitude with respect to the distance at closest approach, from N-body simulations (dots) and secular theoretical approach (line), for the closer fly-bys, and for an initial planetary apoastron of 1000 au. The grey line indicates the projected elevation of the planet.

Fig. 12 Maximal altitude with respect to the distance at closest approach, from N-body simulations (dots) and secular theoretical approach (line), for the closer fly-bys, and for an initial planetary apoastron of 3000 au. The grey line indicates the projected elevation of the planet.

3.3.4 Discussion

From both approaches, theoretical and numerical, in the most favorable case, it appears that a fly-by has a significant impact on the planet (periastron increase above 1 au) only if its closest approach is less than 0.05 pc, that is 10 000 au. This corresponds to a small subset among the initial draws, not because of an incompatibility with the observations, but because of the high dispersion of closest approaches, underconstrained by the observations.

We checked that the distance at closest approach is not correlated to the time at closest approach, nor to the velocity at closest approach. Considering the compatibility between our results and the dynamical scenario proposed in Rodet et al. (2017), the time of the fly-by must be considered. Given our simulations, the closest approach occurred likely 2 to 4 Myr ago (3 ± 1 Myr). However, our scenario account for the ejection of the planet only in the beginning of the system life, when protoplanetary disk is still present and can effectively trigger planetary migration. Given the disk lifetime for massive stars (~ 3 Myr, Ribas et al. 2015)and the system assumed age (15 Myr), 2 to 4 Myr ago is significantly too late for the fly-by to have a decisive role. However, a younger age for the system (10 Myr, compatible with LCC age spread of 6 Myr) could still account for this discrepancy.

3.4 Effect on the disk

The effects of a fly-by on a disk may be significant, depending on the parameters of the encounter. The case of a dynamically efficient fly-by can be observed in system HD 141569, where the ongoing encounter has been deeply studied in Reche et al. (2009). In this system, the fly-by could be responsible for truncation, spiral formation, collisional evolution, eccentricity and inclination raise. In our study, the effect of the fly-by on test-particles will be essentially similar to that on the planet. Since the test particles in a debris disk have a nearly circular orbit, the fly-by will increase the eccentricity, significantly or not depending on the distance of closest approach. Moreover, all fly-by characteristics being equal, particles inclination will be excited differently depending on their distance to the host star. The disk might then be warped. The sensitivity of the scattered-light images of the disk are not sufficient to reveal a weak warp, but the warp can induce further instabilities and asymmetries in the disk that could account for its non-standard shape.

We chose among the previous cases a situation with a very short distance at closest approach (1000 au), with a medium relative inclination (~45°) and ran a simulation with the three massive bodies (HD 106906 ABb and the perturbers) and 1000 test particles. The particles have initially semi-major axes evenly shared between 10 and 600 au, eccentricity below 0.05, and an inclination spread of 2 degrees. The simulation was launched for 100 000 yr around the fly-by epoch, with a time step of 1 yr. The resulting disk is represented in Fig. 13.

On the other hand, the repeating passing of the planet within the disk would have stronger consequences. If a very small percentage is ejected over one period (≲0.01%), the mean eccentricity of the particles raise from 0.02 at each passage. For the disk to remain long-lived in its current shape, Jílková & Zwart (2015; non collisional simulations) and Nesvold et al. (2017; collisional simulations) estimated that the planet orbit should not cross the disk. Thus, the planet periastron should be outer to the observed ~100 au outer disk radius. Within our scenario, it means that this enlargement of the periastron occurred rather quickly, whether or not it was caused entirely by the fly-by. In any case, the planet interactions would have cover the track of the fly-by-induced perturbations

The new structure of the code allows to estimate the percentage of dust capture by the planet. It turns out that temporary (less than 10 yr) capture is experienced by about 5% of the dust at each passage, but no permanent captures were produced.

Fig. 13 Orbital elements of the test particles after a close fly-by. The lightly blue zones represent the initial configuration.

4 Conclusion

In this paper, we present the N-body mixed-variable code ODEA, that is able to study multiple systems in evolving architectures. We use it tostudy the rare planetary system HD 106906. We confirm that the two stars identified by De Rosa & Kalas (2019) could have helped stabilizing the planet after a destabilization by its host binary star. This scenario could account for the wide separation of the planet, its possible elevation with respect to the disk plane, as well as the structures evidenced within the disk.

However, the significance of the encounter strongly depends on the distances at closest approach. With the current precision on the three systems configuration (especially the relative radial velocities and distances), it is not possible to establish the role of the flybys. To circularize the planet orbit if it was previously ejected on a wide trajectory, a fly-by closer than 0.05 pc is needed (assuming apoastron ≤3000 au), which is oneorder of magnitude below the uncertainty on the closest approach. The simulations show that the angular

configuration is favorable when this condition is met.

Any indication of HD 106906 b relative motion would be helpful to constrain its orbit, and thus its dynamical history. More precise parallaxes and radial velocities for HIP 59716 and HIP 59721 are necessary to constrain the distances at closest approach, and conclude on the effect of the fly-bys on the system dynamical evolution.

ODEA handles hierarchy changes in systems with non-Solar-system-type architectures. It can model efficient captures and fly-bys. Through a criterion based on accelerations ratios, a new hierarchy is defined when the current is perturbed. ODEA’s natural upgrade is the implementation of a Mercury-like approach to handle close encounters, that is transitional states of non-Keplerian movements.

Acknowledgements

We thank the anonymous referee for reviewing our work and for insightful and constructive comments which improved the manuscript. The project is supported by CNRS, by the Agence Nationale de la Recherche (ANR-14-CE33-0018, GIPSE), the OSUG@2020 labex and the Programme National de Planétologie (PNP, INSU) and Programme National de Physique Stellaire (PNPS, INSU). Most of the computations presented in this paper were performed using the Froggy platform of the CIMENT infrastructure (https://ciment.ujf-grenoble.fr), which is supported by the Rhône-Alpes region (GRANT CPER07_13 CIRA), the OSUG@2020 labex (reference ANR10 LABX56) and the Equip@Meso project (reference ANR-10-EQPX-29-01) of the programme Investissements d’Avenir, supervised by the Agence Nationale pour la Recherche. P.K. and R.J.D.R. thank support from NSF AST-1518332, NASA NNX15AC89G and NNX15AD95G/NEXSS. This work benefited from NASA’s Nexus for Exoplanet System Science (NExSS) research coordination network sponsored by NASA’s Science Mission Directorate.

Appendix A Theoretical error associated with the symplectic mapping

Splitting the Hamiltonian with a kick-drift-kick approach, as described in Sect. 2.1, the energy error that we get is (Saha & Tremaine 1994) (A.1)

where the Poisson brackets are defined as follows: (A.2)

In Swift HJS, H_A and H_B are given by Eqs. (5) and (6), which can be computed respectively with O(N) and O(N²) operations. Thus, computing {{H_A, H_B}, H_A} already requires O(N⁴) operations.

Appendix B Test simulation for ODEA

The simulation that was used for the study of the performance of ODEA belongs to the 10 000 simulations performed in Sect. 3.2, in the case where the two perturbers were bound. We chose a simulation with a small distance of closest approach (~2000 au), so that there is an effect on the planet orbital elements. The fly-bys occurred between 4 and 2 Myr ago, so that we restricted ourselves to this time interval when studying the energy error with respect to the time-scale (to limit the floating point round-off error associated with the large distances).

Appendix C Derivation of the changes of planet periastron due to the fly-by in the impulse approximation

C.1 Circular impulse

The expression of the change of the planet velocity is given in Eq. (9). Supposing that the new orbit intersects the old one at apoastron or periastron, we have Δa_p = a_pΔe_p. Moreover, the velocity of the planet if on a circular orbit is . Thus, the eccentricity is

and the periastron is then given by Δperi = Δa_p − aΔe_p = −2a_pΔe_p.

C.2 Apoastron impulse

Stating instead that the apoastron is preserved, one gets Δa_p = −a_pΔe_p∕(1 + e_p). Within the impulse framework, the change of velocity involves the velocity at apoastron, so that the velocity writes . Thus,

which gives

and the periastron is then given by Δperi = Δa_p(1 − e_p) − aΔe_p = −2a_pΔe_p∕(1 + e_p).

Appendix D Derivation of the changes of planet orbital characteristics due to the fly-by in the secular approximation

D.1 Perturbative potential

We inspire from Heggie & Rasio (1996) to derive the first-order perturbation of the planet orbital elements in the secular approximation.

Following Heggie & Rasi, we number respectively 1, 2 and 3 HD 106906 central star, HD 106906 b and one of the stellar perturber. The position of the planet relative to its host star is denoted by r, and the position of the third body relative to HD 106906 center of mass is denoted by R. In this framework, the evolution of the planet orbit verifies:

where U is the perturbative potential.

In the secular approximation, U is averaged over the orbit of HD 106906 planetary orbit. The implicit assumptions is that all orbital elements but the anomaly have a longer evolution timescale than the orbital period. As we are interested in the first order evolution, we only integrate the dominant part in a_p∕a (quadripole order). Then, we use Lagrange equations to retrieve the evolution of the eccentricity, the inclination and the longitude of periastron.

D.2 Eccentricity and periastron change

After we first averaged over the planet orbital motion, the secular evolution of the eccentricity obtained at the quadrupole level writes:

where the x–y plane is the initial plane of the planet (plane of the disk), and the x direction is given by the planet initial periastron. To compute the first order of the change of e after the fly-by, we integrate de∕dt along time from −∞ to + ∞ by fixing all variables to their initial values but the angular evolution of the stellar perturber.

Heggie & Rasio computed in their Eq. (7) the change in eccentricity as a function of the angular parameters of the encounter, and we exactly retrieve their expression. The maximum efficiency is obtained for a coplanar encounter, where all the transferred angular momentum apply only on the eccentricity. Stating that the eccentricity of the perturber’s orbit is significantly more than 1 (V = 3 km s⁻¹ and D = 1 pc gives e ~ 500, D = 0.1 pc gives e ~ 50), we obtain

where Ω is the longitude of the ascending node and ω the argument of the periastron of the perturber hyperbolic orbit. The maximum is obtained for Ω + ω = π∕4. Moreover, the eccentricity e depends on D, V and GM_tot as so that . Thus, the eccentricity change satisfies: (D.1)

On the other hand, the semi-major axis is constant in the secular approximation. The periastron is then given by Δperi = −aΔe_p.

D.3 Inclination change

The secular evolution of the inclination obtained at the quadrupole level writes:

We then integrate as before to compute the change of inclination Δi_p.

The maximum is reached for i = π∕4, Ω = π∕2 and ω = π∕2. Thus, we obtain

D.4 Longitude of the periastron change

The secularevolution of the total longitude of the periastron obtained at the quadrupole level writes:

We then integrate as before to compute the change of inclination .

The maximum is reached for i = π∕2, Ω = 0 and ω = 0. Thus, we obtain

D.5 Maximal altitude

The maximum altitude z_max reached by the planet on its orbit is given as a function of its orbital elements: (D.2)

It thus depends on the evolution of a_p, e_p, i_p and ω_p.

Due to the term sin(i_p), the same approach than above leads to neglecting all evolution but that of the inclination. It is consistent with the fact that in the previous expressions, Δ_ip ≫ Δe_p, Δi_p when the eccentricity tends to 1. We get:

However, this estimate is not valid anymore when Δ_ip approaches π∕2, that is when sin(i_p) approaches 1. At this point, the estimates of Δe_p and must be taken into account. In order to comprise all the different evolution scales, we thus simply estimate the maximal altitude by replacing directly the computed evolution in the definition formula: (D.5)

where , and .

Appendix E Additional materials for HD 106906 fly-by simulations

Figures E.1 and E.2 represents the distribution of the coordinates of the bodies in the simulations. Figures E.3 and E.4 describe the case where the two perturbers are bound. The coordinates of the bodies are drawn from the observational constraints with the same process that for the non-bound case, but we discarded the configurations where the eccentricity of the relative orbit is greater than 1. The resulting semi-major axis and eccentricity distributions are presented here, along with the effect of the fly-bys on the planet periastron, which is very similar to the non-bound case.

Fig. E.1 Initial distribution (today) of HD 106906 coordinates and velocities for the 10 000 simulated cases (green), and for the 359 cases where a fly-by closer than 0.1 pc occurred (red).

Fig. E.2 Initial distribution (today) of HIP 59721 coordinates and velocities for the 10 000 simulated cases (green), and for the 359 cases where a fly-by closer than 0.1 pc occurred (red).

Fig. E.3 Semi-major axis and eccentricity distributions for the relative orbit of the two perturbers HIP 59716 and HIP 59721, assuming they are bound.

Fig. E.4 Periastron increase with respect to the distance at closest approach, from N-body simulations (dots) and theoretical approaches (lines), for the closer fly-bys, for an initial planetary apoastron of 1000 au, in the case where the two perturbers are bound.

References

All Tables

All Figures

	Fig. 1 Example of hierarchy change in the case of a capture. At first the red body orbits the yellow–blue pair. After a strong interaction, it captures the small blue body.
In the text

	Fig. 2 Maximum relative energy error over a 1 Myr evolution of HD 106906 including the two fly-bys, for ODEA and Swift HJS, as a function of the time-step assumed.
In the text

	Fig. 3 Representation of HD 106906, HIP 59716 and HIP 59721 current positions and velocities in HD 106906 rest frame (disk lies in the Y –Z plane, observed extension in the Y -direction).
In the text

	Fig. 4 Two dimension histograms of the coordinates of the intersection points between the perturbers trajectories and the X–Y plane, assuming linear trajectories.
In the text

	Fig. 5 Representation of a typical evolution of the hierarchy in the three-body simulations of HD 106906 fly-bys with ODEA. All orbits here are hyperbolic.
In the text

	Fig. 6 Distribution of the distances at closest approach. The following study will focus on the red part, that corresponds to fly-by closer than 0.1 pc (3.6% of the configurations).
In the text

	Fig. 7 Distribution of the times and velocities at closest approach, for the cases where the distance at closest approach is less than 0.1 pc.
In the text

	Fig. 8 Initial distribution (today) of HIP 59716 coordinates and velocities for the 10 000 simulated cases (green), and for the 359 cases where a fly-by closer than 0.1 pc occurred (red).
In the text

	Fig. 9 Periastron increase with respect to the distance at closest approach, from N-body simulations (dots) and theoretical approaches (lines), for the closer fly-bys, and for an initial planetary apoastron of 1000 au. The grey part corresponds to a periastron change inferior to + 1 au, which will not secure the planet stability.
In the text

	Fig. 10 Periastron increase with respect to the distance at closest approach, from N-body simulations (dots) and theoretical approaches (lines), for the closer fly-bys, and for an initial planetary apoastron of 3000 au. The grey part corresponds to a periastron change inferior to + 1 au, which will not secure the planet stability.
In the text

	Fig. 11 Maximal altitude with respect to the distance at closest approach, from N-body simulations (dots) and secular theoretical approach (line), for the closer fly-bys, and for an initial planetary apoastron of 1000 au. The grey line indicates the projected elevation of the planet.
In the text

	Fig. 12 Maximal altitude with respect to the distance at closest approach, from N-body simulations (dots) and secular theoretical approach (line), for the closer fly-bys, and for an initial planetary apoastron of 3000 au. The grey line indicates the projected elevation of the planet.
In the text

	Fig. 13 Orbital elements of the test particles after a close fly-by. The lightly blue zones represent the initial configuration.
In the text

	Fig. E.1 Initial distribution (today) of HD 106906 coordinates and velocities for the 10 000 simulated cases (green), and for the 359 cases where a fly-by closer than 0.1 pc occurred (red).
In the text

	Fig. E.2 Initial distribution (today) of HIP 59721 coordinates and velocities for the 10 000 simulated cases (green), and for the 359 cases where a fly-by closer than 0.1 pc occurred (red).
In the text

	Fig. E.3 Semi-major axis and eccentricity distributions for the relative orbit of the two perturbers HIP 59716 and HIP 59721, assuming they are bound.
In the text

	Fig. E.4 Periastron increase with respect to the distance at closest approach, from N-body simulations (dots) and theoretical approaches (lines), for the closer fly-bys, for an initial planetary apoastron of 1000 au, in the case where the two perturbers are bound.
In the text