Issue 
A&A
Volume 631, November 2019



Article Number  A139  
Number of page(s)  14  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201935728  
Published online  08 November 2019 
ODEA: Orbital Dynamics in a complex Evolving Architecture
Application to the planetary system HD 106906
^{1}
Université Grenoble Alpes, CNRS, IPAG,
38000
Grenoble,
France
email: laetitia.rodet@univgrenoblealpes.fr
^{2}
Kavli Institute for Particle Astrophysics and Cosmology, Stanford University,
Stanford,
CA
94305,
USA
^{3}
Department of Astronomy, University of California,
Berkeley,
CA
94720,
USA
Received:
18
April
2019
Accepted:
10
September
2019
Context. Mixedvariable symplectic integrators are widely used in orbital dynamics. However, they have been developed for Solar systemtype architectures, and can not handle evolving hierarchy, in particular in systems with two or more stellar components. Such configuration may have occurred in the history of HD 106906, a tight pair of Ftype stars surrounded by a debris disk and a planetarymass companion on a wide orbit.
Aims. We present the new algorithm ODEA, based on the symplectic algorithm Swift HJS, that can model any system (binary,...) with unstable architecture. We study the peculiar system HD 106906 as a testcase for the code.
Methods. We define and compute a criterion based on acceleration ratios to indicate when the initial hierarchy is not relevant anymore. A new hierarchy is then computed. The code is applied to study the two recently evidenced flybys that occurred on system HD 106906, to determine if they could account for the wide orbit of the planet. Thousands of simulations have been performed to account for the uncertainty on the perturbers coordinates and velocities.
Results. The algorithm is able to handle any change of hierarchy, temporary or not. We used it to fully model HD 106906 encounters. The simulations confirm that the flybys could have stabilized the planet orbit, and show that it can account for the planet probable misalignment with respect to the disk plane as well as the disk morphology. However, that requires a small distance at closest approach (≲0.05 pc), and this configuration is not guaranteed.
Conclusions. ODEA is a very good choice for the study of nonSolar type architecture. It can now adapt to an evolving hierarchy, and is thus suitable to study capture of planets and dust. Further observations of the perturbers, in particular their radial velocity, are required to conclude on the effects of the flyby on system HD 106906.
Key words: methods: numerical / celestial mechanics / planets and satellites: dynamical evolution and stability / planets and satellites: individual: HD 106906 / planet–star interactions / stars: kinematics and dynamics
© L. Rodet et al. 2019
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
In the context of the rapid increase of exoplanet discoveries, the need for efficient Nbody 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 Nbody integrators: first, they exhibit no longterm 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 timestep than other integrators for the same result. In 1991, Wisdom and Holman devise the first symplectic map specifically designed for Nbody 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 timestep, 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^{1}.
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 ScorpiusCentaurus (ScoCen) 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 hotstart 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 flyby 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 flyby 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 Nbody problem, with masses ${({m}_{i})}_{i=1,\mathrm{..},N}$, positions ${({r}_{i})}_{i=1,\mathrm{..},N}$ and impulsions ${({p}_{i})}_{i=1,\mathrm{..},N}$. The Hamiltonian is $$H={\sum}_{i=1}^{N}\frac{{p}_{i}{}^{2}}{2{m}_{i}}{\sum}_{1\le i<j\le N}\frac{G{m}_{i}{m}_{j}}{{r}_{ij}},$$(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$\tilde{H}$. 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 longterm drift of the energy.
In order to design a proper $\tilde{H}$ in orbital mechanics, the key idea is to split the Hamiltonian into two integrable parts: $$H={H}_{\text{A}}+{H}_{\text{B}}.$$(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 roundoff errors. $\tilde{H}$ corresponds to the successive integration of these parts separately. For a second order symplectic integrator, a socalled leapfrog 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 nonempty 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 sodefined 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 ${{m}^{\prime}}_{k}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}{\mu}_{k}{\eta}_{k}/{M}_{k}$.
In this formalism, a new set of N coordinates ${({{r}^{\prime}}_{k},{{p}^{\prime}}_{k})}_{i=1,\mathrm{..},N}$ are designed with a Jacobilike approach: ${{r}^{\prime}}_{k}$ is the relative position of the center of mass of orbit k’s satellites with respect to that of its centers, and ${{p}^{\prime}}_{k}$ is the relative conjugate momentum. The first coordinates ${{r}^{\prime}}_{1}$ and ${{p}^{\prime}}_{1}$ 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 $$\begin{array}{ll}{{r}^{\prime}}_{k}\hfill & ={\sum}_{i\text{,satellitesofk}}\frac{{m}_{i}{r}_{i}}{{\mu}_{k}}{\sum}_{i\text{,centersofk}}\frac{{m}_{i}{r}_{i}}{{\eta}_{k}},\hfill \\ {{p}^{\prime}}_{k}\hfill & ={{m}^{\prime}}_{k}\left({\sum}_{i\text{,satellitesofk}}\frac{{p}_{i}}{{\mu}_{k}}{\sum}_{i\text{,centersofk}}\frac{{p}_{i}}{{\eta}_{k}}\right).\hfill \end{array}$$
The Hamiltonian can then be split as follows $$\begin{array}{ll}{H}_{\text{A}}\hfill & ={\sum}_{k=2}^{N}\frac{{{p}^{\prime}}_{k}{}^{2}}{2{{m}^{\prime}}_{k}}\frac{G{\mu}_{k}{\eta}_{k}}{{{r}^{\prime}}_{k}};\hfill \\ {H}_{\text{B}}\hfill & ={\sum}_{k=2}^{N}\frac{G{\mu}_{k}{\eta}_{k}}{{{r}^{\prime}}_{k}}{\sum}_{1\le i<j<\le N}\frac{G{m}_{i}{m}_{j}}{{r}_{ij}}.\hfill \end{array}$$
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 ${a}_{k}^{B}\equiv 1/{{m}^{\prime}}_{k}\partial {H}_{\text{B}}/\partial {{r}^{\prime}}_{k}$.
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 twodimensional symmetric array that compiles the Keplerian acceleration between two bodies ${a}_{k}^{\text{Kep}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}G{M}_{k}/{r}_{ij}^{2}$, 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 timestep Δt. We choose a Keplerianlike time ${\mathrm{min}}_{k}{T}_{k}/20$, where $${T}_{k}=\sqrt{\frac{4{\pi}^{2}{a}_{k}^{3}\mid 1{e}_{k}{\mid}^{3}}{G{M}_{k}}}$$(7)
if orbit k is bound or if its smallest approach has not yet occurred, or $${T}_{k}=\sqrt{\frac{4{\pi}^{2}{{r}^{\prime}}_{k}{}^{3}}{G{M}_{k}}}$$(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^{3})). Checking for a possible change at each timestep, 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 rebuilding 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^{4}, see Appendix B). The criterion must be fast to compute (maximum as N^{3}, 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 ${a}_{k}^{B}/{a}_{k}^{\text{Kep}}$ for each orbit k, where ${a}_{k}^{\text{Kep}}$ and ${a}_{k}^{B}$ are the accelerations $\stackrel{\mathrm{..}}{{{r}^{\prime}}_{k}}$ 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^{3}) 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 nonoptimal 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 Nbody 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 flybys, for ODEA and Swift HJS, as a function of the timestep assumed. 
2.5 Comparison with other codes
Several algorithms have been introduced since the formalization of the first mixedvariable symplectic map for orbital mechanics, including the widely used Mercury (Chambers 1999). Most of them are designed to work in SolarSystemlike 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 mixedvariable integrator that we know of are designed to handle long or definitive hierarchy change. Such situations can be encountered in case of a stellar flyby, 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 flyby and dust capture. Figure 2 illustrates the gain of energy precision that ODEA allows on the HD 106906’s flybys 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 highorder 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 flybys, 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 BulirschStoer implementation of Press et al. (1989), with a precision constraint on the trajectory of order 10^{−7} (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 flyby configuration. The Kepler solver used to integrate H_{A} handles bound or unbound orbits as well.
Nbody simulations of the 15 Myr past evolution of HD 106906 including the two flybys.
3 Application to system HD 106906
3.1 Characterizing the perturbers
Searching for potential stellar perturbers in ScoCen 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^{−1}) are four times higher than the velocity dispersion reported for LCC (1.13 ± 0.07 km s^{−1}; Madsen et al. 2002), that was used in Rodet et al. (2017). We will see in Sect. 3.3 that the effect of a flyby 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
Nbody 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 MonteCarlo 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^{−1}), the proper motion of the declination μ_{δ} (0.05 mas yr^{−1}) and the radial velocity γ (up to 1.7 km s^{−1}). 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 flyby 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 timestep 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) flyby 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 flybys is roughly the same. HIP 59716 coordinates distributions are presented on Fig. 8. Most of the parameters of the configurations with close flybys are drawn randomly within the configurations, except for the radial velocity, where we see that the configurations leading to a close flyby 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^{−1}.
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 threebody simulations of HD 106906 flybys with ODEA. All orbits here are hyperbolic. 
Timescales of the HD 106906 simulations.
Fig. 6 Distribution of the distances at closest approach. The following study will focus on the red part, that corresponds to flyby 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 flyby 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 timestep 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 flyby 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 flyby. 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 $\frac{1}{2}G{M}_{\text{HD106906}}/{a}_{p}$, where a_{p} is the initial semimajor 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^{2} yr^{−2}. A proportion of 2 × 10^{−3} less corresponds to an elliptic trajectory with apoastron 1000 au, and 2 × 10^{−4} M_{⊙}au^{2} yr^{−2} 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 flyby 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 timescale of the flyby 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 midplane. 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 flyby 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 flyby closer than 0.1 pc occurred (red). 
3.3.2 Results
The conclusion of the study depends essentially on the possibility for the flyby 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 flyby 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 semimajor 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 Nbody simulations (dots) and theoretical approaches (lines), for the closer flybys, 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 Nbody simulations (dots) and theoretical approaches (lines), for the closer flybys, 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 flyby 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 flyby timescale with the orbital period of the planet, this approximation often provides a good estimate. In this framework, Brunini & Fernandez (1996) show that the flyby increases the planet velocity by: $$\Delta {v}_{p}\lesssim \frac{2G{M}_{*}}{V{D}^{2}}{a}_{p},$$(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 semimajor 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 semimajor axis Δa_{p} = −a_{p}Δe_{p}, which gives a change of periastron Δperi = −2a_{p}Δe_{p}. Finally, one gets (see Appendix C): $$\Delta \text{peri}\lesssim 8\frac{G{M}_{*}}{\sqrt{G{M}_{\text{HD106906}}}}\frac{{a}_{p}^{\frac{5}{2}}}{V{D}^{2}}.$$(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): $$\Delta \text{peri}\lesssim 8\frac{G{M}_{*}}{\sqrt{G{M}_{\text{HD106906}}}}\frac{{a}_{p}^{\frac{5}{2}}}{V{D}^{2}}\frac{\sqrt{(1{e}_{p})(1+{e}_{p})}}{3{e}_{p}}.$$(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 semimajor axis is invariant throughout the flyby. Considering a coplanar orbit and a perturber’s eccentricity significantly higher than 1 (strongly unbound orbit), the maximum is: $$\Delta \text{peri}\lesssim \frac{5}{2}\frac{G{M}_{*}}{\sqrt{G{M}_{\text{HD106906}}}}\frac{{a}_{p}^{\frac{5}{2}}}{V{D}^{2}}{e}_{p}\sqrt{1e{p}^{2}}.$$(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^{−1}. 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^{−3} (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 flyby 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: $${z}_{\text{max}}={a}_{p}\mathrm{sin}({i}_{p})\left(\sqrt{1{e}_{p}^{2}{\mathrm{cos}}^{2}({\omega}_{p})}+{e}_{p}\mathrm{sin}({\omega}_{p})\right).$$(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 Nbody simulations (dots) and secular theoretical approach (line), for the closer flybys, 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 Nbody simulations (dots) and secular theoretical approach (line), for the closer flybys, 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 flyby 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 flyby 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 flyby 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 flyby on a disk may be significant, depending on the parameters of the encounter. The case of a dynamically efficient flyby can be observed in system HD 141569, where the ongoing encounter has been deeply studied in Reche et al. (2009). In this system, the flyby could be responsible for truncation, spiral formation, collisional evolution, eccentricity and inclination raise. In our study, the effect of the flyby on testparticles will be essentially similar to that on the planet. Since the test particles in a debris disk have a nearly circular orbit, the flyby will increase the eccentricity, significantly or not depending on the distance of closest approach. Moreover, all flyby 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 scatteredlight 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 nonstandard 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 semimajor 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 flyby 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 longlived 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 flyby. In any case, the planet interactions would have cover the track of the flybyinduced 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 flyby. The lightly blue zones represent the initial configuration. 
4 Conclusion
In this paper, we present the Nbody mixedvariable 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 flyby 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 flybys on the system dynamical evolution.
ODEA handles hierarchy changes in systems with nonSolarsystemtype architectures. It can model efficient captures and flybys. 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 Mercurylike approach to handle close encounters, that is transitional states of nonKeplerian 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 (ANR14CE330018, 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.ujfgrenoble.fr), which is supported by the RhôneAlpes region (GRANT CPER07_13 CIRA), the OSUG@2020 labex (reference ANR10 LABX56) and the Equip@Meso project (reference ANR10EQPX2901) 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 AST1518332, 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 kickdriftkick approach, as described in Sect. 2.1, the energy error that we get is (Saha & Tremaine 1994) $$\tilde{H}=H\frac{\Delta {t}^{2}}{12}\left\{\{{H}_{\text{A}},{H}_{\text{B}}\},{H}_{\text{A}}+\frac{1}{2}{H}_{\text{B}}\right\}+O(\Delta {t}^{4}),$$(A.1)
where the Poisson brackets are defined as follows: $$\{f,g\}={\sum}_{i}\frac{\partial f}{\partial {r}_{i}}\frac{\partial g}{\partial {p}_{i}}\frac{\partial f}{\partial {p}_{i}}\frac{\partial g}{\partial {r}_{i}}.$$(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^{2}) operations. Thus, computing {{H_{A}, H_{B}}, H_{A}} already requires O(N^{4}) 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 flybys 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 timescale (to limit the floating point roundoff error associated with the large distances).
Appendix C Derivation of the changes of planet periastron due to the flyby 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 ${v}_{p}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\sqrt{G{M}_{\text{HD106906}}/{a}_{p}}$. Thus, the eccentricity is $$\Delta {e}_{p}=\frac{\Delta {a}_{p}}{{a}_{p}}=2\frac{\Delta {v}_{p}}{{v}_{p}}\lesssim 4\frac{G{M}_{*}}{\sqrt{G{M}_{\text{HD106906}}}}\frac{{a}_{p}^{\frac{3}{2}}}{V{D}^{2}}$$
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 ${v}_{p}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\sqrt{G{M}_{\text{tot}}/{a}_{p}}\sqrt{(1e)/(1+e)}$. Thus, $$\frac{\Delta {v}_{p}}{{v}_{p}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{\Delta {a}_{p}}{2{a}_{p}}\frac{\Delta {e}_{p}}{1e{p}^{2}}=\frac{\Delta {e}_{p}}{2(1+{e}_{p})}\frac{\Delta {e}_{p}}{1e{p}^{2}}$$
which gives $$\begin{array}{ll}\Delta {e}_{p}\hfill & =2\frac{\Delta {v}_{p}}{{v}_{p}}\frac{1{e}_{p}^{2}}{3{e}_{p}}\hfill \\ \hfill & \lesssim 4\frac{G{M}_{*}}{\sqrt{G{M}_{\text{HD106906}}}}\frac{{a}_{p}^{\frac{3}{2}}}{V{D}^{2}}\frac{{(1+ep)}^{\frac{3}{2}}\sqrt{1{e}_{p}}}{3{e}_{p}}\hfill \end{array}$$
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 flyby in the secular approximation
D.1 Perturbative potential
We inspire from Heggie & Rasio (1996) to derive the firstorder 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: $$\begin{array}{l}\ddot{r}=\frac{G{M}_{12}}{{r}^{3}}r+\nabla U\\ U=\frac{G{m}_{3}{M}_{12}}{{m}_{1}{m}_{2}}\left(\frac{{m}_{2}}{R\frac{{m}_{1}}{{M}_{12}}r}\frac{{m}_{1}}{R+\frac{{m}_{2}}{{M}_{12}}r}\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}=\frac{{G}_{m}3{r}^{2}}{2{R}^{3}}\left(3{\left(\frac{r.R}{rR}\right)}^{2}1\right)+O\left({\left(\frac{r}{R}\right)}^{3}\right)\end{array}$$
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: $$\frac{\text{d}{e}_{p}}{\text{d}t}=\frac{15G{m}_{3}{R}_{x}{R}_{y}{a}_{p}^{\frac{3}{2}}{e}_{p}\sqrt{1e{p}^{2}}}{2{R}^{5}\sqrt{G{M}_{12}}}$$
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 flyby, 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^{−1} and D = 1 pc gives e ~ 500, D = 0.1 pc gives e ~ 50), we obtain $$\Delta {e}_{p}=\frac{5}{2}\frac{{M}_{*}}{\sqrt{{M}_{\text{HD106906}}{M}_{\text{tot}}}}\frac{{a}_{p}^{\frac{3}{2}}}{{D}^{\frac{3}{2}}}\frac{{e}_{p}\sqrt{1{e}_{p}^{2}}}{\sqrt{e}}\mathrm{sin}(2\Omega +2\omega ),$$
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 $V\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\sqrt{G{M}_{\text{tot}}(1+e)/D}$ so that $\sqrt{e}\simeq V\sqrt{D/G{M}_{\text{tot}}}$. Thus, the eccentricity change satisfies: $$\Delta {e}_{p}\lesssim \frac{5}{2}\frac{G{M}_{*}}{\sqrt{G{M}_{\text{HD106906}}}}\frac{{a}_{p}^{\frac{3}{2}}}{V{D}^{2}}{e}_{p}\sqrt{1e{p}^{2}}.$$(D.1)
On the other hand, the semimajor 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: $$\frac{\text{d}{i}_{p}}{\text{d}t}=\frac{3G{m}_{3}{a}_{p}^{\frac{3}{2}}\left(4{e}_{p}^{2}+1\right){R}_{x}{R}_{z}}{2{R}^{5}\sqrt{\left(1{e}_{p}^{2}\right)G{M}_{12}}}.$$
We then integrate as before to compute the change of inclination Δi_{p}. $$\begin{array}{ll}\Delta {i}_{p}=\hfill & \frac{3}{2}\frac{G{M}_{*}}{\sqrt{G{M}_{\text{HD106906}}}}\frac{{a}_{p}^{\frac{3}{2}}}{V{D}^{2}}\frac{1+4{e}_{p}^{2}}{\sqrt{1{e}_{p}^{2}}}\hfill \\ \hfill & \times (\mathrm{cos}(i)\mathrm{sin}(\Omega )\left(\mathrm{arccos}\left(\frac{1}{e}\right)+\sqrt{{e}^{2}1}\right)\hfill \\ \hfill & (\mathrm{cos}(\Omega )\mathrm{sin}(2\omega )+\mathrm{cos}(i)\mathrm{sin}(\Omega )\mathrm{cos}(2\omega ))\frac{{({e}^{2}1)}^{\frac{3}{2}}}{3{e}^{2}}).\hfill \end{array}$$
The maximum is reached for i = π∕4, Ω = π∕2 and ω = π∕2. Thus, we obtain $$\Delta {i}_{p}\lesssim \frac{G{M}_{*}}{\sqrt{G{M}_{\text{HD106906}}}}\frac{{a}_{p}^{\frac{3}{2}}}{V{D}^{2}}\frac{1+4{e}_{p}^{2}}{\sqrt{1{e}_{p}^{2}}}.$$
D.4 Longitude of the periastron change
The secularevolution of the total longitude of the periastron $\overline{{\omega}_{p}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}{\omega}_{p}+{\Omega}_{p}$ obtained at the quadrupole level writes: $$\frac{\text{d}\overline{{\omega}_{p}}}{\text{d}t}=\frac{3G{m}_{3}{a}_{p}^{\frac{3}{2}}\sqrt{\left(1{e}_{p}^{2}\right)}({R}^{2}4{R}_{x}^{2}+{R}_{y}^{2})}{2{R}^{5}\sqrt{G{M}_{12}}}.$$
We then integrate as before to compute the change of inclination $\Delta {\overline{{\omega}_{p}}}_{p}$. $$\begin{array}{ll}\Delta \overline{{\omega}_{p}}=\hfill & \frac{1}{4}\frac{G{M}_{*}}{\sqrt{G{M}_{\text{HD106906}}}}\frac{{a}_{p}^{\frac{3}{2}}}{V{D}^{2}}\sqrt{1{e}_{p}^{2}}\hfill \\ \hfill & \times \text{\hspace{0.17em}}(6{\mathrm{cos}}^{2}(i){\mathrm{cos}}^{2}(\omega )5(\mathrm{cos}(2i)3){\mathrm{cos}}^{2}(\omega )\mathrm{cos}(2\Omega )\hfill \\ \hfill & +2\mathrm{cos}(2i)(35\mathrm{cos}(2\Omega )){\mathrm{sin}}^{2}(\omega )\hfill \\ \hfill & 10\mathrm{cos}(i)\mathrm{sin}(2\omega )\mathrm{sin}(2\Omega )\text{vphantom}{)}^{2}).\hfill \end{array}$$
The maximum is reached for i = π∕2, Ω = 0 and ω = 0. Thus, we obtain $$\Delta \overline{{\omega}_{p}}\lesssim 5\frac{G{M}_{*}}{\sqrt{G{M}_{\text{HD106906}}}}\frac{{a}_{p}^{\frac{3}{2}}}{V{D}^{2}}\sqrt{1{e}_{p}^{2}}.$$
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: $${z}_{\text{max}}={a}_{p}\mathrm{sin}({i}_{p})\left(\sqrt{1{e}_{p}^{2}{\mathrm{cos}}^{2}({\omega}_{p})}+{e}_{p}\mathrm{sin}({\omega}_{p})\right).$$(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, Δ_{i}p ≫ Δe_{p}, Δi_{p} when the eccentricity tends to 1. We get: $$\begin{array}{ll}\Delta {z}_{\text{max}}\hfill & ={a}_{p}\sqrt{1{e}_{p}^{2}}\Delta {i}_{p}\hfill \\ \hfill & \lesssim \frac{G{M}_{*}}{\sqrt{G{M}_{\text{HD106906}}}}\frac{{a}_{p}^{\frac{5}{2}}}{V{D}^{2}}(1+4{e}_{p}^{2}).\hfill \end{array}$$
However, this estimate is not valid anymore when Δ_{i}p approaches π∕2, that is when sin(i_{p}) approaches 1. At this point, the estimates of Δe_{p} and $\Delta \overline{\omega}$ 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: $$\Delta {z}_{\text{max}}\lesssim {a}_{p}\mathrm{sin}\left(\tilde{{i}_{p}}\right)\left(\sqrt{1{\tilde{{e}_{p}}}^{2}{\mathrm{cos}}^{2}\left(\tilde{\omega}\right)}+\tilde{{e}_{p}}\left\mathrm{sin}\left(\tilde{\omega}\right)\right\right),$$(D.5)
where $\tilde{{i}_{p}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\mathrm{max}(\Delta {i}_{p},\frac{\pi}{2})$, $\tilde{{e}_{p}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}{e}_{p}\Delta {e}_{p}$ and $\tilde{\omega}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\mathrm{max}(\Delta {\overline{\omega}}_{p},\frac{\pi}{2})$.
Appendix E Additional materials for HD 106906 flyby 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 nonbound case, but we discarded the configurations where the eccentricity of the relative orbit is greater than 1. The resulting semimajor axis and eccentricity distributions are presented here, along with the effect of the flybys on the planet periastron, which is very similar to the nonbound 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 flyby 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 flyby closer than 0.1 pc occurred (red). 
Fig. E.3 Semimajor 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 Nbody simulations (dots) and theoretical approaches (lines), for the closer flybys, for an initial planetary apoastron of 1000 au, in the case where the two perturbers are bound. 
References
 Bailey, V., Meshkat, T., Reiter, M., et al. 2013, ApJ, 780, L4 [NASA ADS] [CrossRef] [Google Scholar]
 Baruteau, C., Bai, X., Mordasini, C., & Mollière, P. 2016, Space Sci. Rev., 205, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Beust, H. 2003, A&A, 400, 1129 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brown, A., Vallenari, A., Prusti, T., et al. 2018, A&A, 616, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brunini, A., & Fernandez, J. A. 1996, A&A, 308, 988 [NASA ADS] [Google Scholar]
 Chambers, J. E. 1999, MNRAS, 304, 793 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Chambers, J. E., Quintana, E. V., Duncan, M. J., et al. 2002, AJ, 123, 2884 [NASA ADS] [CrossRef] [Google Scholar]
 Daemgen, S., Todorov, K., Quanz, S. P., et al. 2017, A&A, 608, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 De Rosa, R. J., & Kalas, P. 2019, AJ, 157, 125 [NASA ADS] [CrossRef] [Google Scholar]
 De Zeeuw, P., Hoogerwerf, R., de Bruijne, J. H., Brown, A., & Blaauw, A. 1999, ApJ, 117, 354 [Google Scholar]
 Heggie, D. C., & Rasio, F. A. 1996, MNRAS, 282, 1064 [NASA ADS] [CrossRef] [Google Scholar]
 Jílková, L., & Zwart, S. P. 2015, MNRAS, 451, 804 [NASA ADS] [CrossRef] [Google Scholar]
 Kalas, P. G., Rajan, A., Wang, J. J., et al. 2015, ApJ, 814, 32 [NASA ADS] [CrossRef] [Google Scholar]
 Lagrange, A.M., Langlois, M., Gratton, R., et al. 2016, A&A, 586, L8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lagrange, A.M., Mathias, P., Absil, O., et al. 2019, A&A, submitted [Google Scholar]
 Levison, H. F., & Duncan, M. J. 1994, Icarus, 108, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Madsen, S., Dravins, D., & Lindegren, L. 2002, A&A, 381, 446 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nesvold, E. R., Naoz, S., & Fitzgerald, M. P. 2017, ApJ, 837, L6 [NASA ADS] [CrossRef] [Google Scholar]
 Pecaut, M. J., & Mamajek, E. E. 2016, MNRAS, 461, 794 [NASA ADS] [CrossRef] [Google Scholar]
 Pecaut, M. J., Mamajek, E. E., & Bubar, E. J. 2012, ApJ, 746, 154 [NASA ADS] [CrossRef] [Google Scholar]
 Press, W. H., Flannery, B. P., Teukolsky, S. A., et al. 1989, Numerical Recipes (Cambridge: Cambridge University Press) [Google Scholar]
 Reche, R., Beust, H., & Augereau, J.C. 2009, A&A, 493, 661 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rein, H., & Spiegel, D. S. 2014, MNRAS, 446, 1424 [Google Scholar]
 Rein, H., Hernandez, D. M., Tamayo, D., et al. 2019, MNRAS, 485, 5490 [NASA ADS] [CrossRef] [Google Scholar]
 Ribas, Á., Bouy, H., & Merín, B. 2015, A&A, 576, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rodet, L., Beust, H., Bonnefoy, M., et al. 2017, A&A, 602, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Saha, P., & Tremaine, S. 1994, ArXiv eprints [arXiv:astroph/9403057] [Google Scholar]
 Wisdom, J., & Holman, M. 1991, AJ, 102, 1528 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Nbody simulations of the 15 Myr past evolution of HD 106906 including the two flybys.
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 flybys, for ODEA and Swift HJS, as a function of the timestep 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 threebody simulations of HD 106906 flybys 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 flyby 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 flyby closer than 0.1 pc occurred (red). 

In the text 
Fig. 9 Periastron increase with respect to the distance at closest approach, from Nbody simulations (dots) and theoretical approaches (lines), for the closer flybys, 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 Nbody simulations (dots) and theoretical approaches (lines), for the closer flybys, 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 Nbody simulations (dots) and secular theoretical approach (line), for the closer flybys, 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 Nbody simulations (dots) and secular theoretical approach (line), for the closer flybys, 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 flyby. 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 flyby 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 flyby closer than 0.1 pc occurred (red). 

In the text 
Fig. E.3 Semimajor 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 Nbody simulations (dots) and theoretical approaches (lines), for the closer flybys, for an initial planetary apoastron of 1000 au, in the case where the two perturbers are bound. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.