ODEA: Orbital Dynamics in a complex Evolving Architecture -- Application to the planetary system HD 106906

Mixed-variable symplectic integrators are widely used in orbital dynamics. However, they have been developed for Solar system-type 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 F-type stars surrounded by a debris disk and a planetary-mass companion on a wide orbit. 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. 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 fly-bys that occurred on system HD 106906, recently evidenced by De Rosa&Kalas (2019), 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. 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 fly-bys 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. ODEA is a very good choice for the study of non-Solar 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 fly-by on system HD 106906.


Symplectic algorithms
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 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.

HD 106906
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 or more 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.

Structure of the code: Swift HJS
Let us consider the gravitational N-body problem, with masses (m i ) i=1,..,N , positions (r i ) i=1,..,N and impulsions 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 HamiltonianH. The latter is chosen 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 properH in orbital mechanics, the key idea is to split the Hamiltonian into two integrable parts: 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.H 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 (r k , p k ) i=1,..,N are designed with a Jacobi-like approach: r k is the relative position of the center of mass of orbit k's satellites with respect to that of its centers, and p k is the relative conjugate momentum. The first coordinates r 1 and p 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 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 Hamiltonian describing 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 B k ≡ 1/m k ∂H B /∂r k .

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 a Kep k = GM k /r 2 ij , 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 min k T k /20, where if orbit k is bound or if its smallest approach has not yet occurred, or 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.

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 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 objective estimate of the relevance of the numerical scheme. However, its computation is tedious (grows as N 4 , see Appendix). 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 B k /a Kep k for each orbit k, where a Kep k and a B k are the accelerations r 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.

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 consider the 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.

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-Systemlike 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. Fig. 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 the Appendix). 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 by Table 2, 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 to resolve 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 −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 fly-by configuration. The Kepler solver used to integrate H A handles bound or unbound orbits as well.
3. Application to system HD 106906

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 subsection 3.3 that the effect of a fly-by is inversely proportional to 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). (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 study at that point.

N-body simulations performed by De Rosa & Kalas
The initialization of the simulations is designed with a Monte-Carlo approach, following De Rosa & Kalas. The 3×6 parameters and their respective precision are the 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 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

Objects
Timescale (yr) Host binary star period 10 −1 Planet period 10 3 Duration of the fly-by 10 5 Perturbers binary period 10 6 Time of the fly-by 3.10 7 yr ago Age of the system (15 ± 6).10 7 10 2 10 1 10 0 10 1 Distance at closest approach (pc) 10 4 10 5 10 6 Distance at closest approach (au) 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). 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 1,000 yr, with outputs every 1,000 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. 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 Fig. 14 and 15 in the appendix.
The distributions of the time and velocities of the perturber at closest approach are represented on 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.

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 1,000 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 remained 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 1,000 au, and 95 % for an apoastron of 3,000 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 1,000 (a = 500.5 au, e = 0.998) or 3,000 au (a = 1500.5 au, e = 0.9993).
The necessary energy to completely eject the planet is 1 2 GM HD106906 /a p , 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 2 /yr 2 . A proportion of 2.10 −3 less corresponds to an elliptic trajectory with apoastron 1,000 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 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 timescale 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.

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 1,000 or 3,000 au apoastron, a 0.1 pc encounter is not enough to significantly raise the periastron: a closer fly-by is required. For the 1,000 au apoastron case, the distance at closest approach must be less than 0.01 pc, that is 2, 000 au. For the 3,000 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 such distances, 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 1,000 and 3,000 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 degrees, 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 minimal altitude 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.

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: 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): 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): (1 − ep)(1 + e p ) 10 1 Distance at closest approach (pc) 10 4 10 3 10 2 10 1 10 0 10 1 10 2 10 3 10 4 10 5 0 10 5 10 4 10 3 10 2 10 1 10 0 10 1 10 2 10 3 10 4 Increase of planet periastron (AU) HIP 59716 HIP 59721 Bound Ejected Secular Impulse (Circular) Impulse (Apoastron) 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 1,000 au. The grey part corresponds to a periastron change inferior to +1 au, which will not secure the planet stability.
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: 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 −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).
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 10 2 10 1 Distance at closest approach (pc) 10 5 10 4 10 3 10 2 10 1 10 0 10 1 10 2 10 3 10 4 10 5 0 10 5 10 4 10 3 10 2 10 1 10 0 10 1 10 2 10 3 10 4 10 5 10 6 Increase of planet periastron (AU) HIP 59716 HIP 59721 Bound Ejected Secular Impulse (Circular) Impulse (Apoastron) 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 3,000 au. The grey part corresponds to a periastron change inferior to +1 au, which will not secure the planet stability. . 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 1,000 au. The grey line indicates the projected elevation of the planet. significant elevation above the disk plane. For any Keplerian orbit, the maximum elevation z max above the reference plane is given by: 1 − e 2 p cos 2 (ω p ) + e p | sin(ω p )| .
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 the appendix. The resulting maximal altitude is represented on Fig.11 and 12. 10 1 Distance at closest approach (pc)

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.

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 (1,000 au), with a medium relative inclination (∼ 45 • ) and ran a simulation with the three massive bodies (HD 106906 ABb and the perturbers) and 1,000 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 years around the fly-by epoch, with a time step of 1 yr. The resulting disk is represented on 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.

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 to study 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 ≤ 3,000 au), which is one order 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 efficiently 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. Fig. 14 and 15 represents the distribution of the coordinates of the bodies in the simulations. Fig. 16 and 17 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.    10 1 Distance at closest approach (pc) 10 4 10 3 10 2 10 1 10 0 10 1 10 2 10 3 10 4 10 5 0 10 5 10 4 10 3 10 2 10 1 10 0 10 1 10 2 10 3 10 4 Increase of planet periastron (AU) Bound Secular Impulse (Circular) Impulse (Apoastron) Fig. 17. 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 1,000 au, in the case where the two perturbers are bound.