Issue 
A&A
Volume 602, June 2017



Article Number  A12  
Number of page(s)  15  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201630269  
Published online  24 May 2017 
Origin of the wideorbit circumbinary giant planet HD 106906
A dynamical scenario and its impact on the disk
^{1} Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France
email: laetitia.rodet@univgrenoblealpes.fr
^{2} Instituto de Astronomia, Geofísica e Ciências Atmosféricas, Universidade de São Paulo, Rua do Matão, 1226, Cidade Universitária, 05508900, São Paulo − SP, Brazil
^{3} Laboratoire d’Astrophysique de Bordeaux, Univ. Bordeaux, CNRS, B18N, Allée Geoffroy SaintHilaire, 33615 Pessac, France
Received: 16 December 2016
Accepted: 4 March 2017
Context. A giant planet has been recently resolved at a projected distance of 730 au from the tight pair of young (~13 Myr) intermediatemass stars HD 106906AB in the Lower Centaurus Crux (LCC) group. The stars are surrounded by a debris disk which displays a ringlike morphology and strong asymmetries at multiple scales.
Aims. We aim to study the likelihood of a scenario where the planet formed closer to the stars in the disk, underwent inward diskinduced migration, and got scattered away by the binary star before being stabilized by a close encounter (flyby).
Methods. We performed semianalytical calculations and numerical simulations (Swift_HJS package) to model the interactions between the planet and the two stars. We accounted for the migration as a simple force. We studied the LCC kinematics to set constraints on the local density of stars, and therefore on the flyby likelihood. We performed Nbody simulations to determine the effects of the planet trajectories (ejection and secular effects) onto the disk morphology.
Results. The combination of the migration and meanmotion resonances with the binary star (often 1:6) can eject the planet. Nonetheless, we estimate that the flyby hypothesis decreases the scenario probability to less than 10^{7} for a derived local density of stars of 0.11 stars/pc^{3}. We show that the concomitant effect of the planet and stars trajectories induce spiralfeatures in the disk which may correspond to the observed asymmetries. Moreover, the present disk shape suggests that the planet is on an eccentric orbit.
Conclusions. The scenario we explored is a natural hypothesis if the planet formed within a disk. Conversely, its low probability of occurrence and the fact that HD 106906 b shares some characteristics with other systems in ScoCen (e.g., HIP 78530, in terms of mass ratio and separation) may indicate an alternative formation pathway for those objects.
Key words: methods: numerical / celestial mechanics / planetary systems / planets and satellites: dynamical evolution and stability / planetdisk interactions
© ESO, 2017
1. Introduction
More than 3500 exoplanets have been found in the last three decades^{1}, but few among them have been detected to be hundreds of astronomical units (au) from their star. As the development of direct imaging reveals more of those wide planetarymass companions, classical theories of planet formation fail at explaining their origin. In the two scenarios, core accretion (Pollack et al. 1996) and gravitational instability (Boss 1997), the planets form within the primordial gas disk. However, the limited extent of the disk (see e.g. Fig. 5 in LiemanSifry et al. 2016) does not enable the formation of a giant planet far away from its star. Thus, when the star around which orbits the very wide and massive HD 106906AB b turned out to be a binary star (Lagrange et al. 2016b), it has been suggested that dynamical interactions could account for the current position of the planet (Lagrange et al. 2016b; Wu et al. 2016).
The planet HD 106906 (or also HIP 59960) is located at a distance of 103 ± 4 pc (Van Leeuwen 2007) and belongs to the Lower Centaurus 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 17 Myr, with an agespread of about 10 Myr (Pecaut et al. 2012). In recent years, high contrast imaging has revealed the circumstellar environment of HD 106906AB: an 11 ± 2 M_{J} planet located at 732 ± 30 au in projected separation (Bailey et al. 2013) and an asymmetric debris disk nearly viewed edgeon, imaged by SPHERE (Lagrange et al. 2016a), GPI and HST (Kalas et al. 2015) and MagAO (Wu et al. 2016). More recently, the binary nature of HD 106906 was inferred thanks to observations with the instruments HARPS and PIONIER (Lagrange et al. 2016b). It turns out to be a 13 ± 2 Myr old SB2 binary consisting of two F5 Vtype stars with very similar masses. Table 1 summarizes the key characteristics of the system components. No further information is known about the orbit of the planet, which must have an orbital period of at least 3000 years. The binary orbit is also not much constrained yet, but given its short orbital period (<100 days), it will presumably be better known in the near future.
Key characteristics of the HD 106906 system.
The edgeon debris disk has an unusual shape: its luminous intensity has a very asymmetric profile. The longest peak, pointing westward, extends up to 550 au, while the east edge reaches 120 au only (see Figs. 1 and 3 of Kalas et al.2015). Conversely, below 120 au, the disk is more luminous on its east side than on its west side. This reversed asymmetry might suggest the presence of a spiral density wave extending over the whole disk, and viewed edgeon from the Earth. Finally, a large cavity splits the disk into two debris belts. Chen et al. (2014) modeled the stars’ excess emission and suggested 13.7 and 46 au for the radii of the belts. The latter likely corresponds to the one imaged by Lagrange et al. (2016a) and Kalas et al. (2015) at ~50 au.
Despite the richness of the observations, the geometry and kinematics of the whole system are strongly underconstrained. If the actual planetbinary distance is less than 1000 au, then the orbit inclination with respect to the plane of the disk must be significant (20 degrees). However, a coplanar configuration cannot be excluded, but the separation should then be around 3000 au. In any case, the large separation between the planet and the central binary, as well as the possible misalignment between the planet orbit and the debris disk, challenges classical mechanisms of planet formation.
According to current theories, planet formation takes place in the primordial gaseous disk. However, as we mentioned above, forming a giant planet via core accretion or gravitational instability at 700 au or more from any central star appears very unlikely, 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. The disk asymmetries (in particular the suspected spiral structure) indicate strong ongoing dynamical interaction with the dust. This may suggest that the planet did not form where it resides today, but may have formed inside and be scattered afterwards. The recently discovered binary nature of HD 106906AB is indeed a source of potentially strong dynamical perturbations that could trigger planet scattering.
The purpose of this paper is to investigate both analytically and numerically the scenario that could have lead to the presentday characteristics of the HD 106906 system starting from a planet formation within the circumbinary disk. As viscosityinduced migration tends to make the planet move inwards, in Sect. 2 we will study the likelihood of an ejection via interactions with the binary, and we will then discuss in Sect. 3 how the planet could have stabilized on such a wide orbit. Finally, in Sect. 4 we briefly analyze the effect of this scattering scenario on the disk and the processes that could have shaped it as it currently appears. Numerical simulations in our analysis have been performed using the Swift_HJS symplectic integration package (Beust 2003), a variant of the Swift package developed by Levison & Duncan (1994), but dedicated to multiple stellar systems.
2. Ejection
2.1. Basic scenario
We investigated how HD 106906 b, supposed initially orbiting the binary on a nearly coplanar orbit, could have been ejected from the disk via dynamical interactions. When it is located far away enough from its host stars, a circumbinary planet may have a very stable orbit. On the other hand, if it migrates too close to the binary, it undergoes a close encounter with the stars and can be ejected.
The binary is surrounded by a chaotic zone where no stable circumbinary orbit is possible. Dvorak (1986) uses a semianalytical approach to compute the upper critical orbit (lower radius of the stable zone) and lower critical orbit (upper radius of the chaotic zone) around two stars of same masses for different binary eccentricities, and found that this gap size typically ranges between two and three times the semimajor axis of the binary orbit. Numerical results for this mass ratio are missing, so that we computed the limits of the gap with our Swift_ HJS package and compared them to the semianalytical approach in Fig. 1. In each simulation, the evolutions of 10 000 test particles have been studied during 10^{5} orbital periods of the binary. The particles have been randomly chosen with semimajor axes between 1.5 and four times the binary semimajor axis a_{B}, eccentricities between 0 and 0.1, and inclinations with respect to the binary orbital plan between 0 and 3°. The time step has been chosen to be 1/20 of the binary orbital period.
Fig. 1
Chaotic zone (in dark gray) as a function of the binary eccentricity, for binary components of same masses. The lighter part designates a critical zone, where some test particles can survive. The red lines represent the lower and upper critical orbit parabolic fits found by Dvorak (1986) in its study of circumbinary planet stability. The 1:6 commensurability is the strongest outside the chaotic zone (see Sect. 2.2) for e_{B} ≥ 0.4. 
Artymowicz & Lubow (1994) showed that this chaotic zone also affects the gas of the disk, with gap sizes similar to the values given by our algorithm (Beust 2003). Consequently, as the migration necessarily stops at the inner edge of the disk, the planet should never reach the chaotic zone this way. It will remain confined close to the lower critical orbit, where it may never be ejected. Meanmotion resonances (hereafter MMR) may help overcoming this difficulty. During its inward migration, the planet is likely to cross MMRs with the binary. It may then be captured by the resonance and furthermore undergo an eccentricity increase that could drive its periastron well inside the chaotic zone.
2.2. Meanmotion resonances
Nested orbits are in a configuration of MMR when their orbital periods are commensurable. For fixed masses and neglecting the precession, this is fully controlled by the semimajor axis ratio a_{B}/a (subscript B refers to the binary): the orbits are said to be in a p + q: p resonance when (1)
where p, q are integers, and T and m designate respectively Keplerian periods and masses. Resonances are described using the characteristic angle (2)
where λ designates the mean longitude and ω the periastron longitude. σ represents the longitude of a conjunction between the binary and the planet, where all three bodies are aligned, measured from the line of apsides of the planet. If σ stops circulating and begins to oscillate around an equilibrium position (libration), it means that the conjunctions repeatedly occur roughly at the same places on the planet orbit: the system is locked in the resonant configuration. If the resonant conjunction occurs in the location where the interacting bodies are sufficiently far away from each other (like in the NeptunePluto case), then the resonance acts as a stabilizing mechanism that prevents close encounters. MMRs are nevertheless known to enhance eccentricities. If the eccentricities are too highly excited, then the conjunctions may no longer occur at safe locations, often causing instability. For a review on MMRs, see Morbidelli (2002).
The way a MMR can enhance the eccentricity of the planet can be studied in a semianalytical way. Details about this procedure are given in Beust (2016), Beust & Morbidelli (2000) and Yoshikawa (1989). Basically, if we restrict the study to orbits with negligible σlibration, the interaction Hamiltonian can be averaged over the motion of the binary for constant σ. This gives a one degree of freedom autonomous Hamiltonian. Phasespace diagrams with level curves of this Hamiltonian can then be drawn in (ν ≡ ω−ω_{B},e) space to explore the overall dynamics. To adapt the method to this unusual case where the inner bodies have similar masses, we calculated the resonant Hamiltonian of a planet orbiting the center of mass of a binary with binary mass parameter μ ≡ m_{2}/m_{B} (where m_{2} is the mass of the second star): (3)
where G is the gravitational constant, r_{B} ≡ R_{2}−R_{1} and r ≡ R−(μR_{1} + (1−μ)R_{2}) if R, R_{1} and R_{2} are respectively the position vectors of the planet, the first and the second component of the binary. We could then perform the integration over the orbital motions and derive the phase space diagram for the interesting commensurabilities. The result is displayed in Fig. 2 in the 1:6 MMR case, for a binary eccentricity of e_{B} = 0.4. Most of the level curves of the Hamiltonian exhibit important change in the planet eccentricity; therefore, starting at low eccentricity, the resonant interaction can drive the planet to higher eccentricity regime and cause it to cross the chaotic zone (indicated in red on the figure) at periastron, leading to ejection.
Fig. 2
Isocontour in the (ν = ω−ω_{B}, e) phase space of the average interaction Hamiltonian of a test particle trapped in 1:6 meanmotion resonance with a binary eccentricity of e_{B} = 0.4, assuming a binary mass parameter of μ = 1/2. Each curve represents a trajectory in the (ν, e) space. Above the red line, the planet has part of its orbit in the chaotic zone. 
Our choice of focusing on the 1:6 meanmotion resonance should not be surprising. Indeed, according to Fig. 1, it is the lowest order resonance that lies outside the chaotic zone for e_{B} ≥ 0.4: it occurs at a/a_{B} ≃ 3.3. Any lower order (thus potentially stronger) resonance such as 1:2, 1:3, etc. falls inside the chaotic zone, and could not be reached by the planet according to our scenario. Moreover, the topology of the diagram depends on the binary eccentricity: the higher it is, the higher is the change of eccentricities depicted by the level curves. And those curves are flat for a circular binary orbit.
However, the semianalytical study is not sufficient here to study the dynamical route that leads to ejection. Indeed, libration of the resonant angle σ and chaos on short timescale, not taken into account in the computation of the phasespace diagram, are not negligible for a binary with mass parameter close to 1/2. We thus performed numerical simulations of 10^{5} binary orbital periods of dynamical evolution for different binary eccentricities and different initial angular conditions, to study the stability of different ratio of MMR. All runs were performed starting with a semimajor axis close to the resonant value, with a time step set to 1/20 of the binary orbital period.
Only a few resonances located outside the chaotic zone are finally able to trigger ejection: the 1:6 and the 1:7 one. The simulations allowed to check not only the ability of the resonances to generate instability, but also the time needed to eject the planet, as well as the typical width of the starting resonant zone that leads to ejection, which is typically 0.01 au. Table 2 summarizes the results obtained with various e_{B} values and μ = 1/2 with the 1:6 resonance.
Effect of the 1:6 meanmotion resonance and ejection duration for different eccentricities of the binary, starting with a planet eccentricity of e = 0.05.
The simulations confirm that resonance stability depends on binary eccentricity e_{B}, and that the resonance gets weaker when the order of the resonance  q  increases. An important result is that whenever ejection occurs, it happens within a very short timescale, always much shorter than the typical time needed (≲1 Myr) to form the planet from the gaseous disk. Our first conclusion is thus that the planet cannot have formed within the resonance. This validates the idea outlined above that it first formed at larger distances in a more stable position, and furthermore migrated inwards and was possibly trapped in a meanmotion resonance before being ejected. In the following, we investigate this scenario.
2.3. Migration
In recent decades, planet migration has become an unavoidable ingredient to explain the configuration of some planetary systems. Due to tidal interactions with the primordial gas disk, giant planets (mass > 10 M_{⊕}) undergo first a type I, and furthermore a type II migration once they have created a gap in the disk (Baruteau et al. 2014). It consists of a drift that can be directed toward the star, whose characteristic timescale depends on the position, and characteristics of the planet and on the viscous properties of the disk.
We have assumed that the planet has approximately reached its final mass when it arrives at the location of unstable MMRs, that is between 1 and 2 au from the stars. The characteristic time of migration varies in inverse proportion to the quantity α_{ν}h^{2}Σ, where α_{ν} is the viscosity parameter, h the aspect ratio and Σ the surface density (Lin & Papaloizou 1986). However, not only the values of those quantities are unknown in HD 106906 primordial disk, but also this simple dependency does not seem to match nor the known planetary population (Mordasini et al. 2009) neither the results inferred by hydrodynamical simulations (Dürmann & Kley 2015). Taking this fact into account, estimating the mass of the primordial disk to be around 0.6% of the stellar mass (Andrews et al. 2013) and varying the viscosity parameter and the aspect ratio around the observed values (e.g., Pinte et al. 2015), we obtained a large range of migration timescales. To obtain the largest overview without trying every single velocity, we choose to run our tests with four different migration velocities at 2 au: 10^{3}, 10^{4}, 10^{5} and 10^{6} au/yr.
Simulating the whole process of diskinduced migration in the circumbinary environment is beyond the scope of the present paper. Using a hydrodynamic code, Nelson (2003) computed the migration of a planet in a circumbinary disk and show that it was likely to get locked into a meanmotion resonance. As their stars had very different masses, their results can not be applied here, so we choose to add to the SWIFT_HJS code an additional extraforce that mimics the migration mechanisms they observed. This force is designed in such a way that its secular effect averaged over the orbital motion of the planet just induces the desired steadystate semimajor axis drift da/ dt = v_{mig}, v_{mig} being a fixed arbitrary migration velocity, and has no effect on the eccentricity nor on the longitude of periastron. Further details about the choice of the force are provided in Appendix A. We derive: (4)
where (e_{r},e_{θ}) are the 2D cylindrical radial and orthoradial unit vectors in the local referential frame attached to the planet’s motion. Thus, F_{mig} depends on the planet position via the radius r, the vector e_{θ} and the planet mean angular motion n = 2π/T. The parameter v_{mig} is set at the beginning of the simulation, according to the timescale we want for the migration. We note that with the above convention, inward migration corresponds to v_{mig}< 0. Of course, the migration is implicitly assumed to hold as long as the planet moves inside the disk.
Whether migration would inhibit or enhance the effects of MMR is not a straightforward issue. Resonance trapping induced by type II migration was found to exist for some commensurabilities between two protoplanets orbiting a star (Snellgrove et al.2001; Nelson & Papaloizou2002). But MMRs with a binary are more difficult to predict, especially those located near the chaotic zone, like those we are focusing on here. Moreover, the expected lifetime of the gas disk is roughly three million years around massive stars (Haisch Jr et al.2001; Ribas et al.2015), so that the formation, migration and hypothetical ejection must all occur by this time.
First unstable resonance and corresponding ejection time for different eccentricities of the binary and different migration velocities, starting with e = 0.05 and a = 2 au.
We thus performed numerical simulations, where the planet was initially put outside (~2 au) the zone of interest. Whether the planet formed just outside the critical zone or whether it migrated toward there is irrelevant, only the values of the orbital elements and migration velocity at the entrance of the zone of interest matter to conclude on the possibility of ejection. Migration was added using the additional force depicted in Eq. (4), with the diverse migration velocity prescriptions described above. The simulations were pursued until the planet gets captured in a meanmotion resonance and furthermore ejected, or until it reaches the inner edge of the disk, that is, the chaotic zone, in the case of no resonant capture. Again, the time step has been taken to be ≲T_{B}/ 20. The main result is that migration, regardless of its velocity or of the binary eccentricity, always leads to a resonant trapping followed by an ejection after a reasonable amount of time spent in the resonance.
Fig. 3
Evolution of the planet semimajor axis with respect to time for a binary eccentricity of e_{B} = 0.2 and a 10^{5} au/yr migration velocity. The semimajor axis of the binary has been set to a_{B} = 0.4 au. The plot illustrates the migration, then ejection, of the planet after it has been trapped into a 1:6 resonance. The effect of the 1:7 resonance, weaker, is also visible on the plot. As the planet has a perturbed Keplerian motion around the binary, the exact locations of MMRs are not straightforward to derive (see Appendix B). 
In Fig. 3, an example of the effect of both migration and resonances is visible via the evolution of the semimajor axis of the planet. The figure illustrates the full dynamical evolution corresponding to v_{mig} = 10^{5} au/yr and e_{B} = 0.2. In addition to high frequency oscillations that illustrate the chaotic nature of the dynamics, we see a gradual semimajor decrease at a speed corresponding to the initial prescription, followed by a capture in the 1:6 MMR resonance that finally leads to ejection. Interestingly, we note a temporary trapping in the 1:7 resonance than occurs before the final capture in the 1:6. The 1:7 resonance appears not to be strong enough to fully inhibit the migration, while the 1:6 does.
Table 3 summarizes the ejection times obtained in the various cases tested. Comparing Tables 3 and 2, we note that migration, despite causing important smallscale variability of the semimajor axis, enhances resonant instabilities. However, this efficiency is probably overestimated because of the simplicity of our migration model. Deeper analysis of the diskplanet interaction close to the resonance would be needed. Moreover, close to its inner edge, the disk is strongly shaped by the binary and some eccentric ringlike features may affect the protoplanet migration (Mutter et al. 2016).
We may now summarize the analysis that has been conducted in this section by reviewing the time evolution of this tentative ejection process. The formation of a giant planet takes a variable amount of time depending on the process and the location: from several periods if formed via gravitational instabilities to a million periods if formed via core accretion (Chabrier et al. 2014). Consequently, in order for HD 106906 b to acquire its mass, it must have formed in a relatively stable location over the timescale involved, at least at a distance of 2 au. However, as giant planets are believed to form beyond the snow line, whose location is estimated to ~10 au around ~3 M_{⊙} star (Kennedy & Kenyon 2008), the stability of the planet formation position is a priori ensured. After a substantial growth of the planet, migration occurs, whose strength depends on the primordial disk characteristics, and pushes the planet into a less stable zone. For the planet to be ejected, it has to enter the zone of destabilizing resonances (1:6, 1:7), which lies around 1.5 au (Fig. 3). All in all, if a_{formation}/v_{mig} is inferior to the disk lifeexpectancy, the scattering of the planet is a natural outcome in a system with binary mass ratio close to unity.
3. Stabilization
3.1. The idea of a close encounter
In the previous subsection, we demonstrate that a giant planet which formed reasonably close to the binary is likely to undergo an ejection. However, ejection does not imply stabilization on a distant orbit around the binary, as it is most likely the case for HD 106906 b. Eventually, the planet follows an hyperbolic trajectory and does not need more than 10 000 years to completely fly away from its host star. Indeed, suppose that the planet gets ejected on a stillbound orbit via a close encounter with the binary: the orbit may have a very distant apoastron, but its periastron will necessarily lies in the region where it originates, that is the immediate vicinity of the binary. Therefore, after one orbital period, the planet is back at periastron and undergoes a new violent encounter with the binary that is likely to cause ejection. Such episodes have been actually recorded in our simulations.
Thus, in order to stop the ejection process and stabilize the planet orbit, an additional dynamical process is needed to lower its eccentricity and increase its periastron. In the absence of other wide companion of similar mass orbiting the binary (Lagrange et al. 2016a), a close encounter with a passing star is a natural candidate. Recalling that the ScoCen association must have had a more important stellar density several million years ago, this event might have occurred with non negligible probability.
The impact of dynamical interactions on planetary systems in open clusters has been studied intensively since the discovery of the first exoplanets. An effective cross section has been computed by Laughlin & Adams (1998), that characterizes the minimal encounter distance needed to raise the eccentricity of a Jovian planet at 5 au from 0 to over 0.5. They found ⟨ σ ⟩ = (230au)^{2}, which gives a stellar encounter rate of about 0.01 disruptive encounter in our system lifetime. More precisely, Parker & Quanz (2012) conducted Nbody simulation to observe the planet orbital elements after a flyby, and found a probability between 20 and 25% that a 30 au planet undergoes at least a 10% eccentricity change in a ten million year period. In our case, the situation seems even easier, because we want to modify the orbit of an unstable planet already far from its star, thus with a trajectory that can easily be swayed. However, the encounter needs to happen at the right time of the planet life, during the ~1000 years that would last the ejection. Moreover, the encounter should be weak enough not to definitely eject the planet, but strong enough to circularize the orbit to a reasonable eccentricity. We note that weak encounters are more likely to occur than strong ones.
3.2. Probability of a stabilizing flyby
Fig. 4
Example of a coplanar configuration where a passing star (in red) stabilizes a wide unstable planet orbit. Before the flyby, the planet orbit still has a very low periastron, and after it gets much wider, thanks to the interaction with the passing star. We recall that according to Kepler’s laws, the planet spends most of its time near apoastron, so that any flyby is likely to occur when the planet is at or near this point. 
Of course, not all flyby geometries will generate the desired effect. The flyby is entirely defined by the mass of the passing star M_{∗}, the closest approach (or periastron) to the binary p_{∗}, the velocity of the passing star at closest approach v_{∗}, the inclination i_{∗} of the passing star orbit with respect to the planet orbit, its longitude of ascending node Ω_{∗} measured from the line of apsides of the planet orbit, and the argument of periastron ω_{∗} with respect to the line of nodes. A scheme of the effect of a stellar flyby is sketched in Fig. 4 in the coplanar case, in a configuration voluntarily favorable to a restabilization: when the planet is at the apoastron of a wide unstable orbit. In fact, the apoastron is also the most likely position of the planet, as it spend there most of its time.
Figure 5 shows the results of a parametric study limited to coplanar flybys (we studied the inclined cases as well) for a given angle ω_{∗} (45 degrees), in (p_{∗},v_{∗}) 2D parameter space, for three different M_{∗} values (0.1, 1 and 5 M_{⊙}) and assuming the planet was at the apoastron of a very wide unstable orbit before the encounter (such as in Fig. 4). The planet is 1000 au away from the binary when the flyby occurs, this is why a cutof can be seen around au. In each case, the gray area represents the zone in parameter space that is reachable (plausible v_{∗}) and actually causes a significant periastron increase of the planet. In this peculiar configuration, taking into account the distribution of p_{∗} and v_{∗} (see below), a stabilization is very likely.
Fig. 5
Area in the disruptive star phase space which succeed to raise the periastron from 1 au to over 2 au in the case of a coplanar encounter of periastron argument ω_{∗}−ω = 45 degrees. v_{∗} designates the maximum relative velocity of the disruptive star, p_{∗} designates its smaller distance from the binary. 
As a more general approach, the probability of a convenient encounter can be estimated by an integration over the relevant flyby parameters. Taking a homogeneous distribution of stars in the cluster with characteristic distance d, and a Gaussian distribution of relative velocities with dispersion σ, the number of flybys that would raise the planet periastron above q_{stable} is (5)
where q_{f} is the final periastron reached by the planet after the flyby perturbation, τ_{ejection} is the characteristic time of ejection, τ_{cluster} ≡ d/σ is the characteristic time in the cluster (timescale needed to have a convenient flyby), ℋe is the Heaviside function and n_{v} the velocity distribution of unbound stars (6)This gives the probability of having a stabilizing flyby, for a given mass M_{∗} of the passing star. Apart from the role of M_{∗} (see Fig. 5), this probability is strongly though indirectly dependent on the orbital parameters of the planet before the flyby, that is on the state of advancement of the planet ejection. It is higher when the planet lies initially on a wide, unstable but still bound orbit (as in Fig. 4). On the other hand, it is nearly zero as long as the planet is still close to the binary (i.e., before ejection) and if it is already on a hyperbolic trajectory. In order to compute analytically the value of ℋe(q_{f}−q_{stable}) for every set of parameters (p_{∗},v_{∗},i_{∗},Ω_{∗},ω_{∗}) given any initial planet position and velocity, we assume a linear trajectory for the perturber. The direction of the velocity change caused by this approximated encounter can then be analytically derived, as well as the new planet orbit. In the computation, we assumed a velocity dispersion of σ = 0.2 au/yr (1 km s^{1}) (Madsen et al. 2002), and the order of magnitude of the characteristic time of ejection τ_{ejection} has been set to 10^{3} yr.
The most critical dependence of our formula (5) is on the local distance between stars d. The present and past density of the LCC is not known. Therefore, we attempted to determine it through a kinematic study in Appendix D. From 141 stars for which complete data could be retrieved, we could trace back the density of LCC through time. The results show that the early density was roughly 1.7 times the present density, evaluated around 0.05 star/pc^{3} in the close neighborhood of HD 106906. Moreover, the contribution of field stars (not related to LCC) has been estimated to be similar to the contribution of LCC. From this piece of information, we derived that the present local density is lower than ≈0.11 star/pc^{3}. This density, consistent with the density of the solar neighborhood (Reid et al. 2002), corresponds to d ~ 2 pc. If our scenario happened in such an environment, the probability of a close encounter (p_{∗}< 5000 au) just following the planet ejection is below 1 × 10^{7}.
Nevertheless, our estimate of the LCC density is based on a small number of luminous (and mostly earlytype) stars for which the kinematics can be inferred. In our case, the flyby of any object more massive than the planet can stabilize the orbit and impact our probabilities. Therefore, we considered the extreme case where neighboring bodies in the cluster are separated by d = 0.1 pc, a density similar to the one taken in Laughlin & Adams (1998) and Parker & Quanz (2012). We report the probabilities for that high density and for the case of a 1 M_{⊙} perturber in Table 4, for different initial conditions. We note that the number of encounters for any d> 5000 au roughly scale with d^{3}, so that lowerdensity results can be easily retrieved from the table.
Number of close encounters with a 1 M_{⊙} star raising the planet periastron above a given value (2, 50 or 150 au) depending on the trajectory of the planet before the flyby.
3.3. Conclusions
Table 4 shows that the probability of a stabilizing flyby remains low. As outlined above, the most favorable case corresponds to initially wide elliptical orbits before encounter. However, most of the time, the planet is directly ejected on an hyperbolic orbit instead of a wide elliptical orbit. And even if this occurs, the subsequent periastron passages in the vicinity of the binary quickly lead to a definitive ejection.
The probabilities have been computed for a 1 M_{⊙} perturber only, less than half of our system 2.7 M_{⊙} star. Though the perturber to host star mass ratio do matter to evaluate the flyby impact (e.g., Jílková et al. 2016), the 1 M_{⊙} results give an upper bound that accounts for the encounter with lighter stars, and a rough estimate for encounters with heavier stars (see Fig. 5), which are less abundant.
We therefore conclude that while our scenario uses generic ingredients (migration, MMR, flyby), it is not very likely to happen because of the low probability of a flybyassisted stabilization. An indirect proof could be provided if we could see traces of planet ejection on the disk. Moreover, constraints on the presentday orbit of HD 106906 b would certainly help refining this scenario: a very high planet eccentricity could raise its likelihood, but the secular effect of such a planet passage in the disk every thousands of years would have big consequences on the disk morphology.
Fig. 6
Nbody simulation showing the consequence of the ejection of a 11 M_{J} planet through a disk. From left to right, the snapshots have been taken 0, 1000 and 10 000 years after the ejection. The planet starts on an hyperbolic orbit similar to what we observe in the simulations we performed: a = 10 au and e = 1.1 (corresponding to periastron q = 1 au). Above is a spatial representation of the top view of the disk (the planet trajectory is depicted in black), below is the density along the y axis, integrated over the x and z axis. 
4. Debris disk
In this section, we investigate the consequences of our scattering scenario on the disk particles repartition, to check whether it matches the observations (shortdistance asymmetry, longdistance asymmetry and extended inner cavity).
4.1. Ejection through the disk
An essential part of the scenario we outline in this paper is the violent scattering of the planet by the binary. Most of the time, the planet switches directly from a close orbit around the binary to a fast hyperbolic trajectory toward the edge of the system. As of yet we did not mention the effect of such an ejection on the debris disk surrounding HD 106906AB. The passage of a ~10 M_{Jup} planet across the disk should presumably induce drastic perturbations on it. In order to investigate this issue, we ran a Nbody simulation with 10 000 test particles, neglecting the interactions between them to access the first order of perturbation. The particles have been randomly chosen with semimajor axes between 5 and 100 au, eccentricities between 0 and 0.05, and inclinations with respect to the binary orbital plan between 0 and 2°. As the main effect of the ejection is due to close encounters between the planet and the disk particles, we use the package Swift_RMVS (Levison & Duncan 1994) that is designed to handle such trajectories. However, this package is not devised to work in multiple stellar system, so that the binary will be here approximated by a single star. The binary effect on the dust being negligible above 5 au for the duration of the perturbation (approximately ten times the planet ejection time, that is 10 000 years), this approximation has almost no consequences on the final dust distribution. Time steps have been set to at most 1/20 of the particles orbital periods, but Swift_RMVS automatically adjusts them to manage close encounters.
The result is displayed in Fig. 6. After the initial spirallike propagation of the eccentricity disturbance created by the planet, the disk homogenizes on an oblong asymmetric shape that could possible match the needle we observe up to ~500 au. In the case where the planet is first scattered on a wide eccentric orbit before being ejected, the process gives eccentricity to some test particles, but the effect is negligible compared to the effect of the ejection that comes next. However, in any case, the asymmetry might not last forever. Orbital precession induced by the inner binary (not taken into account in our simulation) should finally randomize the longitudes of periastron of the particles on a much longer timescale and restore the initial axisymmetric disk shape. For a particle orbiting the binary at 100 au, the precession period (see Appendix C) due to the binary is ~4 × 10^{7} yr. Of course it is shorter closer to the star, but this remains comparable or larger than the age of the system except in the innermost parts of the disk. Hence still observing the asymmetry today at 500 au should not be surprising even if was created a long time ago. However, our mechanism cannot explain the reversed asymmetry at shorter distance. This inversion presumably corresponds to a spiral density wave extending across the disk that needs a steadystate perturbation to be sustained over a long enough timescale.
4.2. The effect of a stellar encounter
In Sect. 3, we discussed the possibility that a stellar flyby could have stabilized the planet midejection. The effect of such encounters on a disk has been studied intensively (for example in Breslau et al. 2014; Jílková et al. 2016). This effect is of course very dependent on the mass ratio of the stars and on the encounter periastron and eccentricity.
In fact, most encounters that would stabilize an unstable planet are compatible with the current shape of the disk. We can, for example, consider the case of a 1 M_{⊙} star perturber. The majority of the suitable encounters have periastrons superior to 1000 au (see Fig. 5). According to the computations of Jílková et al. (2016), this and the high mass of our star implies that all the disk particles will remain bound. Indeed, the transfer radius, that is the minimum radius where capture is possible, is well superior to the observed limit of the debris disk. For our disk to be depleted, the transfer radius should be inferior to ≈100 au, which corresponds to an encounter periastron around 250 au. Thus, though the problem is strongly underconstrained, our scenario is likely to be compatible with the existence of the disk.
4.3. Secular carving
The secular action of the planet orbiting the binary on its present day large stabilized orbit is an obvious longterm source of perturbation on the disk. We note that we make here a clear distinction between the initial, short term perturbation triggered by the planet on the disk during its ejection process, which effect has been described in the previous subsection, and the longterm secular action of the planet as it moves on its distant bound orbit. It is known that eccentric companions (planets or substellar) orbiting at large distance a star surrounded by a disk create spiral density waves within the disk (Augereau & Papaloizou 2004). To a lesser extent, binaries do the same on circumbinary disks (Mutter et al. 2016). The following study nevertheless shows that the asymmetry currently observed in the HD 106906 disk cannot be due to the sole action of the binary, but rather requires an outer source of perturbation like the planet, that enhances the density waves induced by the binary.
We investigate here the secular action of the planet on the disk, combined with that of the binary, using simulations with our Swift_HJS package. Of course with only a projected position, our knowledge of the current orbit of the planet is sparse. Some orbital configurations may nevertheless be ruled out as they would lead to a destruction of the disk. Jilkova & Zwart (2015) studied intensively the impact of each orbital configuration on the disk via the percentage of particles that remain bound n_{bound} and the fraction of bound particles that suit the observation constraints f_{d/b}. Although nor the disk neither the binary was resolved at that time, their conclusion still can be used, at least on a qualitative level. They showed that a planet periastron larger than 50 au or an inclination larger than five to ten degrees is enough to keep a relatively good agreement with the observations (f_{d/b}> 0.5) without completely depleting the disk. However, to better match the observations (f_{d/b}> 2/3), the periastron must lie outside the outer radius of the disk. The maximal inclination is constrained by the observation, that is about twenty degrees. No further constraints can be provided by the simulation of Jilkova and Zwart to rule out any inclination between zero and 20 degrees if the planet orbit does not go across the disk. They point out that Kozailike mechanisms can lead to some wobbles in inclination, but small enough for the disk to remain in a nearly coplanar state.
Assuming that the planet fulfills these requirements, we compute the asymmetries induced on the disk and compare the result to the observation. The disk was initially made of 10 000 test particles with same initial conditions than in the previous subsection. The result of a typical run is displayed in Fig. 7. Basically, if the periastron of the planet is close enough to the outer edge of the disk, it generates an important asymmetry in the disk within a timescale of between five and ten million years (~10^{8} binary periods, ~10^{3} planetary orbit). In Fig. 8, the density profile has been computed along the x axis. The resulting plot displays an asymmetry similar to the observations: the east side (in blue) is brighter than the west side on short scale, but its density drops well above the west side density. The shape of the perturbation resembles a circular arc, but it actually consists of two overlapping spiral arcs, one driven by the planet and the other one created by the binary.
Fig. 7
Top view of the evolution of the debris disk after ten million years of perturbation by a planet on a coplanar orbit whose coplanar orbit has a periastron of 200 au and an apoastron of 1000 au. The color scale represents the relative density. Strong asymmetries can be seen. 
Fig. 8
Density along the x axis, integrated over the y and z axis, obtained from Fig. 7. The gray zone marks approximately the cavity that we observe today. The asymmetry seems to reverse when we get farther to the stars. 
This issue can be studied analytically. The approach is analogous to the study without binary, as conducted in Wyatt (2005). Consider a test particle orbiting the binary. Suppose that the planet has a Keplerian orbit around the center of mass of the binary, so that the system is hierarchical. The instantaneous Hamiltonian controlling its motion can be written as (7)where H_{Kep} = −Gm_{B}/ 2a is the pure Keplerian Hamiltonian, and where the two remaining terms constitute the disturbing function, one part arising from the binary, and the other part from the planet. For a binary of mass parameter μ, these independent perturbations are written
where, in a frame whose origin is at the center of mass of the binary, r is the position vector of the particle, r_{B} is the radius vector between the two individual stars, m_{p} is the mass of the planet, and r_{p} is its position vector. More generally, B subscribed quantities will refer to parameters of the binary, p subscribed quantities to the planet, while unsubscribed parameters will correspond to the orbiting particle.
Fig. 9
2D representation of a debris disk after 0, 10^{6} and 10^{7} years under the thirdorder approximation of the influence of the binary and a planet whose periastron is 200 au and apoastron is 1000 au. The color scale represents the local number of particles. At first, only the edge of the disk is affected, but after 10^{7} years, two spirals components appear within the disk. 
Both terms of the disturbing function are then expanded in ascending powers of the semimajor axis ratios a_{B}/a and a/a_{p}, truncated to some finite order (three here) and averaged independently over all orbital motions, assuming implicitly that the particle is not locked in any meanmotion resonance with the planet or with the binary. Higher orders terms of the disturbing function will be neglected on initial examination, but their influence will be studied in a forthcoming paper. The secular evolution of the particle’s orbital elements is then derived via Lagrange equations. Details on this procedure are given in Appendix C.
Starting from a disk made of particles on circular orbits, we use this theory to compute their instantaneous polar coordinates (r(t),θ(t)) in the disk and compute theoretical synthetic images. The result is shown in Fig. 9, which must be compared with Fig. 7. We note the presence of a circular arc very similar to the one obtained in the numerical simulation. This peculiar shape is due to the combination of two spiral waves winded in opposite senses, induced by the planet and the binary via differential precession and eccentricity excitation on the disk particles.
The test particles precession velocities and periastrons are represented in Fig. 10 as a function of their semimajor axis. In the inner part of the disk, the precession is dominated by the binary, so that the speed of the orbital precession decreases with increasing semimajor axis. The results is a trailing spiral wave that can be seen in Fig. 9. Conversely, in the outer part of the disk, the precession is mostly due to the planet, so that its is now an increasing function of the semimajor axis. This creates a leading spiral density wave. The superposition of both spirals in the intermediate region generates the observed circular arc. The exact location of this arc corresponds to the periastron of the particles whose semimajor axis minimizes the precession velocity, that is, around 55 au. Moreover, we see in Fig. 10 that all particles from a ~ 60 au to ~100 au have the same periastron. The combination of the two effects enhances the density of the arc, as can be observed in Figs. 7 and 9.
Fig. 10
Precession velocities (above) and periastrons (below) of Fig. 9 test particles with respect to their semimajor axis (that does not change with time) after ten million years evolution. Above, the red curve displays the precession induced by the binary, the blue curve by the planet and the black curve depicts both contributions. 
The two spirals are, however, not fully independent. Orbital precession of the particles actually has no visible effect on their global distribution in the disk as long as their orbits are circular. In Appendix C, we show that due to the small size of its orbit and its mass ratio close to one, the binary has very little influence on the eccentricity of the particles compared to the planet, even in the inner part of the disk. In fact, while the outer spiral is fully due to the planet, in the inner spiral, the eccentricity oscillations are also driven by the planet, while the precession is controlled by the binary. Moreover, the contrast of the density wave highly depends on the planet orbital shape (the amplitude of the eccentricity oscillations is roughly proportional to e_{p} within our approximation). For example, if the planet apoastron is 1000 au, its periastron should be less than 500 au (e> 0.3) in order to create a significant asymmetry as the one we observe, in a reasonable timescale.
5. Discussion
5.1. Disk cavity
In the previous sections, we did not study the origin of the large cavity observed within the disk by Lagrange et al. (2016a) and Kalas et al. (2015). It is possible that one or more unseen planets could have carved and sustained this cavity. In that case, one of those unseen planet may be responsible for the ejection of the known one, instead of one of the binary star. Those planet(s), if on eccentric orbits, could also influence the shape of the disk (Lee & Chiang 2016). Therefore, we ran Nbody simulations with the Swift_RMVS package (the same setting as in the Debris Disk section) to quantify at first order the minimum mass of a single planet, checking if it can carve the observed cavity between the two belts of debris surrounding the pair of stars. This assumes that if one planet alone is responsible for gap, its mass will be higher than in the case where multiple planets are considered. The end result must reproduce the inner edge of the cavity at 13 Myr located between 10 and 15 au (inferred from the IR excess modeling, Chen et al. 2014). We assume that the outer edge of the cavity around corresponds to the separation of the ring (65 au) measured on the SPHERE images (Lagrange et al. 2016a). The simulations give a minimum mass of 30 M_{J} for a single noneccentric planet located at 30 au, which is well above the detection limits in Fig. 6 of Lagrange et al. (2016a). However, the disk is viewed edgeon, so that the coronagraph used during the observation (radius of 93 mas or 9.5 au) hides part of the orbital plan. In Fig. 11 we computed a 2D detectionlimits map from the data published in Lagrange et al. (2016a). The map confirms that a small zone around the coronagraph has detection sensitivity above 30 M_{J}. An additional giant planet on a 30 au circular orbit will spend 20% of its time (ten years) in this blind zone and therefore could have been missed. For the case of an eccentric orbit, the mass of the perturber could be only 1 M_{J}. This is too low compared to the known planet mass to produce an ejection, but high enough to have a noticeable effect on the disk morphology.
Fig. 11
Minimum mass of planets (in Jupiter masses) that can be detected into the H2 data published in Lagrange et al. (2016a) around HD 106906. The contrasts has been translated into masses using Baraffe et al. (2003) model adapted to the SPHERE filters. 
5.2. Alternative scenarios
Our scenario makes use of standard ingredients (resonances, migration, scattering, flyby) envisioned or observed in young planetary systems (Baruteau et al. 2014) and account for all known components of the system. Nevertheless, the low probability of occurrence we estimate in Sect. 3 because of the need for a nearby star flyby at the right time makes the scenario implausible. If we suppose that the planet was in fact on a stable orbit before the flyby, then this flyby event could have happened at any time, and not necessarily during the early age of the system. However, the probability for a flyby to have a significant effect on the planet without ejecting it decreases dramatically when the planet gets closer to its host star. Taking the data from Parker & Quanz (2012), we can expect a probability of around 0.1 for a disruption superior to 10% on eccentricity without ejection for a 30 au Jovian planet in the system lifetime. Among the disruptive encounters, it is then hard to tell how many would put the planet on a suitable orbit (apoastron greater than 700 au, that is e> 0.75). Plus, such a change of orbit would lead to a very small planet periastron, which will strongly deplete the disk (see Sect. 4.2).
Alternatively, the planet could have been stolen from another system. Indeed, captured planets tend to have eccentric orbit (Malmberg et al. 2011). However, for the final orbit to be so wide, the initial orbit must also have been wide (Jílková et al. 2016). All in all, such a scenario would only turn over the problem, as we would have to account for the wide initial orbit on the first place.
Conversely, the disk could replace the flyby in our scenario. Indeed, to follow the idea of Kikuchi et al. (2014), the planet could have been accelerated by the gas at its apoastron after a first scattering, and its orbit could have been rendered stable this way. It is interesting to note that some of the circumstellar disks of ~2.5 M_{⊙} stars recently resolved with ALMA at high angular resolution shows gas extending up beyond the separation of HD 106906 b (e.g., Walsh et al. 2016, and references therein). The total mass of HD 106906 A and B is around 2.7 M_{⊙} and it is therefore possible that the binary bore such an extended primordial disk that would have circularized the orbit of the ejected HD 106906 b.
Before the discovery of HD 106906AB binary status that indicates strong gravitational interactions, Bailey et al. (2013) suggested that it may have formed in situ. On the one hand, the existence of extended protoplanetary disks implies that our planet may have formed in HD 106906AB primordial disk. On the other hand, HD 106906 b is not the only planetarymass companion detected at very large projected separation, and such bodies have usually no known scatterers in their environment (see Bryan et al. 2016, even though their study was conducted over a small number of systems less wide than HD 106906 and with lighter stars). Among the systems harboring a planetarymass companion of similar separation and mass ratio, we can name HIP 77900^{2} (Aller et al. 2013), HIP 78530 (Lafrenière et al. 2011), or the triple system Ross 458 (Goldman et al. 2010). In Fig. 12, we represented the wide young planetarymass companions discovered by direct imaging. We note that HD 106906 b has the lowest planet/star mass ratio above 100 au. The proximity of HIP 78530A b and HIP 77900A b (two brown dwarfs that are also part of ScoCen) in that diagram, could indicate that HD 106906AB b formed in situ (within the disk, or like a stellar companion).
Fig. 12
Planet mass ratios with respect to projected separation. Only planets that belong to young systems (<0.1 Gyr) are displayed, with the exception of the circumbinary planet Ross 458 c. HD 106906 b has the lowest mass ratio beyond 100 au. HIP 78530 b and HIP 77900 b are its two closest neighbors in the diagram. Being identified as brown dwarfs, they may have formed in situ by cloud collapse. Data retrieved from the exoplanets.eu database. 
6. Conclusion
We have shown that HD 106906AB b could have formed within the primordial disk and be scattered away on a wide orbit during the first ten million years of the system life. This scenario involves the combination of diskinduced migration and meanmotion resonances with the binary. However, if the scattering is likely to occur, the stabilization of the planet on its current wide orbit is delicate, and requires more than gravitational interactions with the binary. A flyby scenario has been suggested, but the stabilization only occurs for a restricted part of the overall encounters trajectories. The low density (<0.11 stars/pc^{3}) that we estimated for the LCC makes a close encounter even more unlikely.
The disk has multiple features, that each could be explained within the frame of our scenario, but also outside of this frame. Two spiral density waves are created if the planet have for the last ten million years had an eccentric orbit with periastron around 200 au. A needle extending to 500 au could have been created by the ejection of the planet, but a smaller needle could be provoked by an eccentric and inclined outer orbit (see Fig. 7) or by an eccentric inner orbit (Lee & Chiang 2016). Nesvold et al. (2017) also studied the secular effect of an eccentric, inclined outer orbit for HD 106906 b in a recent paper, and could produce asymmetries whose brightness repartition is consistent with the observations.
The scenario we explored builds on the observed components of the system (disk, binary star) and on the hypothesis that the planet could not have formed via core accretion or gravitational instability at several hundreds of au. Nevertheless, the low probability of occurrence of our scenario demands that we reconsider those assumptions. Alternative hypothesis like the circularization of the planet orbit via the interaction with the disk gas or in situ formation could explain the present architecture of the system. But this requires that the disk extends up to the separation of the planet and contains enough gas at that separation. Recent high quality images of circumstellar disks extending beyond 700 au around massive stars and the close properties of other systems in ScoCen (HIP 78530A b and HIP 77900A b) argue for this alternative formation pathway.
Finally, we note that many of the methods depicted here are easily generalizable to other circumbinary environment. Nbody simulations with a simple migration force can be applied on any circumbinary planet to have a quick overlook of the stability of its early trajectory. Flyby may not be the most efficient process to stabilize a planet, because of the rarity of suitable close encounters. Destabilization by a flyby is much more probable. Finally, ejection, outer and inner orbits can create huge asymmetries in the disk during the first ten million years of a system. In particular, an inner orbit enhances the dynamical perturbations created by an outer orbit by speeding up the precession, while the outer orbit if eccentric can enhance the perturbations created by the inner orbit by providing eccentricity to the disk.
Contrary to HD 106906 b, HIP 77900 b has not been confirmed by the common proper motion test. Nonetheless, Aller et al. (2013) argue that lowgravity features in HIP 77900 b spectrum is compatible with the object being a member of ScoCen, and therefore a plausible companion to HIP 77900 A.
Acknowledgments
We thank the anonymous referee for reviewing our work and for insightful 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. This work has made use of data from the European Space Agency (ESA) mission Gaia (http://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, http://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research has made use of the SIMBAD database operated at the CDS (Strasbourg, France). P.A.B.G. acknowledges financial support from the São Paulo State Science Foundation (FAPESP). We thank Cecilia Lazzoni and François Ménard for fruitful discussions.
References
 Aller, K. M., Kraus, A. L., Liu, M. C., et al. 2013, ApJ, 773, 63 [NASA ADS] [CrossRef] [Google Scholar]
 Andrews, S. M., Rosenfeld, K. A., Kraus, A. L., & Wilner, D. J. 2013, ApJ, 771, 129 [NASA ADS] [CrossRef] [Google Scholar]
 Artymowicz, P., & Lubow, S. H. 1994, ApJ, 421, 651 [NASA ADS] [CrossRef] [Google Scholar]
 Augereau, J.C., & Papaloizou, J. 2004, A&A, 414, 1153 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 BailerJones, C. A. L. 2015, PASP, 127, 994 [NASA ADS] [CrossRef] [Google Scholar]
 Bailey, V., Meshkat, T., Reiter, M., et al. 2013, ApJ, 780, L4 [NASA ADS] [CrossRef] [Google Scholar]
 Baraffe, I., Chabrier, G., Barman, T. S., Allard, F., & Hauschildt, P. 2003, A&A, 402, 701 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 BarbierBrossat, M., & Figon, P. 2000, A&AS, 142, 217 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Baruteau, C., Crida, A., Paardekooper, S.J., et al. 2014, Protostars and Planets VI, 667 [Google Scholar]
 Beust, H. 2003, A&A, 400, 1129 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Beust, H. 2016, A&A, 590, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Beust, H., & Morbidelli, A. 2000, Icarus, 143, 170 [NASA ADS] [CrossRef] [Google Scholar]
 Boss, A. P. 1997, Science, 276, 1836 [NASA ADS] [CrossRef] [Google Scholar]
 Breslau, A., Steinhausen, M., Vincke, K., & Pfalzner, S. 2014, A&A, 565, A130 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bryan, M. L., Bowler, B. P., Knutson, H. A., et al. 2016, ApJ, 827, 100 [NASA ADS] [CrossRef] [Google Scholar]
 Chabrier, G., Johansen, A., Janson, M., & Rafikov, R. 2014, Protostars and Planets VI, 619 [Google Scholar]
 Chen, C. H., Mamajek, E. E., Bitner, M. A., et al. 2011, ApJ, 738, 122 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, C., Mittal, T., Kuchner, M., et al. 2014, VizieR Online Data Catalog: J/ApJS/211/25 [Google Scholar]
 de Bruijne, J. H. J. 1999, MNRAS, 310, 585 [NASA ADS] [CrossRef] [Google Scholar]
 Desidera, S., Covino, E., Messina, S., et al. 2015, A&A, 573, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 De Zeeuw, P., Hoogerwerf, R. V., de Bruijne, J. H., Brown, A., & Blaauw, A. 1999, ApJ, 117, 354 [Google Scholar]
 Duflot, M., Figon, P., & Meyssonnier, N. 1995, A&AS, 114, 269 [Google Scholar]
 Duriez, L. 1992, Mod. Meth. Cel. Mech., 13, 9 [Google Scholar]
 Dürmann, C., & Kley, W. 2015, A&A, 574, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dvorak, R. 1986, A&A, 167, 379 [NASA ADS] [Google Scholar]
 Galli, P. A. B., Teixeira, R., Ducourant, C., Bertout, C., & BenevidesSoares, P. 2012, A&A, 538, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Galli, P. A. B., Moraux, E., Bouy, H., et al. 2017, A&A, 598, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Girard, T. M., van Altena, W. F., Zacharias, N., et al. 2011, AJ, 142, 15 [NASA ADS] [CrossRef] [Google Scholar]
 Goldman, B., Marsat, S., Henning, T., Clemens, C., & Greiner, J. 2010, MNRAS, 405, 1140 [NASA ADS] [Google Scholar]
 Gontcharov, G. A. 2006, Astron. Lett., 32, 759 [NASA ADS] [CrossRef] [Google Scholar]
 Haisch Jr, K. E., Lada, E. A., & Lada, C. J. 2001, ApJ, 553, L153 [NASA ADS] [CrossRef] [Google Scholar]
 Holmberg, J., Nordström, B., & Andersen, J. 2007, A&A, 475, 519 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
 Jílková, L., & Zwart, S. P. 2015, MNRAS, 451, 804 [NASA ADS] [CrossRef] [Google Scholar]
 Jílková, L., Hamers, A. S., Hammer, M., & Zwart, S. P. 2016, MNRAS, 457, 4218 [NASA ADS] [CrossRef] [Google Scholar]
 Johnson, D. R. H., & Soderblom, D. R. 1987, AJ, 93, 864 [NASA ADS] [CrossRef] [Google Scholar]
 Kalas, P. G., Rajan, A., Wang, J. J., et al. 2015, ApJ, 814, 32 [NASA ADS] [CrossRef] [Google Scholar]
 Kennedy, G. M., & Kenyon, S. J. 2008, ApJ, 673, 502 [NASA ADS] [CrossRef] [Google Scholar]
 Kikuchi, A., Higuchi, A., & Ida, S. 2014, ApJ, 797, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Kordopatis, G., Gilmore, G., Steinmetz, M., et al. 2013, AJ, 146, 134 [Google Scholar]
 Lafrenière, D., Jayawardhana, R., Janson, M., et al. 2011, ApJ, 730, 42 [NASA ADS] [CrossRef] [Google Scholar]
 Lagrange, A.M., Langlois, M., Gratton, R., et al. 2016a, A&A, 586, L8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lagrange, A.M., Mathias, P., Absil, O., et al. 2016b, A&A, submitted [Google Scholar]
 Laskar, J., & Boué, G. 2010, A&A, 522, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Laughlin, G., & Adams, F. C. 1998, ApJ, 508, L171 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, E. J., & Chiang, E. 2016, ApJ, 827, 125 [NASA ADS] [CrossRef] [Google Scholar]
 Levison, H. F., & Duncan, M. J. 1994, Icarus, 108, 18 [NASA ADS] [CrossRef] [Google Scholar]
 LiemanSifry, J., Hughes, A. M., Carpenter, J. M., et al. 2016, ApJ, 828, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Lin, D., & Papaloizou, J. 1986, ApJ, 309, 846 [NASA ADS] [CrossRef] [Google Scholar]
 Lindegren, L., Lammers, U., Bastian, U., et al. 2016, A&A, 595, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Madsen, S., Dravins, D., & Lindegren, L. 2002, A&A, 381, 446 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Malmberg, D., Davies, M. B., & Heggie, D. C. 2011, MNRAS, 411, 859 [Google Scholar]
 Mermilliod, J.C., Mayor, M., & Udry, S. 2009, A&A, 498, 949 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Morbidelli, A. 2002, Modern Celestial Mechanics, Advances in Astronomy and Astrophysics (CRC) [Google Scholar]
 Mordasini, C., Alibert, Y., Benz, W., & Naef, D. 2009, A&A, 501, 1161 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mutter, M. M., Pierens, A., & Nelson, R. P. 2016, MNRAS, in press [Google Scholar]
 Nelson, R. P. 2003, MNRAS, 345, 233 [NASA ADS] [CrossRef] [Google Scholar]
 Nelson, R. P., & Papaloizou, J. C. 2002, MNRAS, 333, L26 [NASA ADS] [CrossRef] [Google Scholar]
 Nesvold, E. R., Naoz, S., & Fitzgerald, M. P. 2017, ApJ, 837, L6 [NASA ADS] [CrossRef] [Google Scholar]
 Parker, R. J., & Quanz, S. P. 2012, MNRAS, 419, 2448 [Google Scholar]
 Pecaut, M. J., & Mamajek, E. E. 2016, MNRAS, 461, 794 [Google Scholar]
 Pecaut, M. J., Mamajek, E. E., & Bubar, E. J. 2012, ApJ, 746, 154 [NASA ADS] [CrossRef] [Google Scholar]
 Pinte, C., Dent, W. R., Menard, F., et al. 2015, ApJ, 816, 25 [Google Scholar]
 Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62 [NASA ADS] [CrossRef] [Google Scholar]
 Preibisch, T., & Mamajek, E. 2008, in Handbook of Star Forming Regions, Vol. II, ed. B. Reipurth, 235 [Google Scholar]
 Reid, I. N., Gizis, J. E., & Hawley, S. L. 2002, AJ, 124, 2721 [NASA ADS] [CrossRef] [Google Scholar]
 Ribas, Á., Bouy, H., & Merín, B. 2015, A&A, 576, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Robin, A. C., Reylé, C., Derrière, S., & Picaud, S. 2003, A&A, 409, 523 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Roeser, S., Demleitner, M., & Schilbach, E. 2010, AJ, 139, 2440 [NASA ADS] [CrossRef] [Google Scholar]
 Snellgrove, M., Papaloizou, J., & Nelson, R. 2001, A&A, 374, 1092 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Song, I., Zuckerman, B., & Bessell, M. S. 2012, AJ, 144, 8 [NASA ADS] [CrossRef] [Google Scholar]
 Torres, C. A. O., Quast, G. R., da Silva, L., et al. 2006, A&A, 460, 695 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Van Leeuwen, F. 2007, A&A, 474, 653 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Walsh, C., Juhász, A., Meeus, G., et al. 2016, ApJ, 831, 200 [NASA ADS] [CrossRef] [Google Scholar]
 Wenger, M., Ochsenbein, F., Egret, D., et al. 2000, A&AS, 143, 9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wilson, R. E. 1953, General catalogue of stellar radial velocities (Carnegie Institute Washington D.C. Publication) [Google Scholar]
 Wu, Y.L., Close, L. M., Bailey, V. P., et al. 2016, ApJ, 823, 24 [NASA ADS] [CrossRef] [Google Scholar]
 Wyatt, M. C. 2005, A&A, 440, 937 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Yoshikawa, M. 1989, A&A, 213, 436 [NASA ADS] [Google Scholar]
 Zacharias, N., Finch, C. T., Girard, T. M., et al. 2012, VizieR Online Data Catalog: I/322 [Google Scholar]
Appendix A: Ad hoc force to account for type II migration
We searched for a coplanar migration force F_{mig} that induces a constant variation of the planet average semimajor axis, but no change in the planet average eccentricity, nor in the periastron longitude. Let F_{mig} = F_{r}e_{r} + F_{θ}e_{θ} be the description of the force in the local referential attached to the planet movement. Gauss equations are (Duriez 1992): where e_{x} and e_{y} are the vectors in the fixed frame and . We want to assume a simple form for F_{r} and F_{θ} that could then be easily averaged over time. The simplest positiondependent force would be F_{r} = A(1 + ccosu) et F_{θ} = B(1 + dcosu), where A, B, c and d are unknown functions of (a,e), constant at first order over a oneperiod integration. Our conditions are then summarized to: where is the mean motion. Taking A = 0, , any c and d = 3e/ 2, we finally obtain Eq. (4).
Appendix B: Location of meanmotion resonances
Equation (1) gives the semimajor axis of a resonant circumbinary planet when its orbit is purely Keplerian. When we take into account the perturbation caused by the binary on the planet orbit, the commensurability of periods that characterizes MMRs can not be easily associated with a semimajor axis, mainly due to orbital precession.
The movement of a circumbinary planet (binary of mass parameter μ) is controlled by the Hamiltonian
where H_{Kep} = −Gm_{B}/ (2a) is the Keplerian Hamiltonian, and where H_{bin} is given by Eq. (8). If the planet orbits at sufficiently large distance from the binary, H_{bin} is a perturbative term that triggers orbital evolution of the planet. This can be investigated analytically via a truncated expansion of H_{bin} in ascending powers of a_{B}/a, and an averaging over both orbital motions. To lowest order, this yields (B.3)Strictly speaking, this approximation is not valid at the exact location of MMRs, as the motions of both orbits are not independent anymore, but it gives a good insight of the perturbation of the planet orbit when it is near the MMRs. Moreover, numerical verifications show that this order two approximation is still relevant for a ≥ 3a_{B}, and could thus be made to study the 1:6 resonance. Lagrange equations (Duriez 1992) then give where n is the Keplerian meanmotion. Thus, if T_{0}(a) is the Keplerian period 2π/n, then the period of the mean longitude T_{λ} is (B.6)The MMR configuration is characterized by the steadiness of σ = (p + q) /qλ_{B}−p/qλ−ω. However, in our study, the planet orbit remains almost circular until ejection, so that the planet line of apsides is not a good reference. Taking the binary line of apsides (constant in time) as the new reference, the resonance characterization writes T_{λ} = p/ (p + q)T_{B}. All in all, the resonant location a_{res} satisfies (B.7)
Appendix C: Spiral density wave
As mentioned in the text, the motion of a particle moving in a circumbinary disk is controlled by the Hamiltonian H_{Kep} + H_{bin} + H_{pla}, where H_{bin} and H_{pla} are perturbative terms given by Eqs. (8) and (9). Following the approach of Wyatt (2005), these terms are then expanded in ascending powers of a_{B}/a and a/a_{p}, truncated to some finite order and averaged over the orbital motion of both orbits (see Laskar & Boué 2010). To second order and third order, the result is In HD 106906 configuration, the binary mass parameter is very close to 1/2, and the semimajor axis of the binary is very small compared to the distance between the binary and the relevant part of the disk, between 50 and 100 au. Thus, the binary part of U_{3} can be neglected. Using Lagrange equations, we derive the equations of evolution: where M is the mean anomaly, is the mean motion and a is a constant of motion in the secular regime. As we want to study the evolution of an initially almost circular particle orbit, we note that we cannot neglect the planetary part of U_{3}, because of the 1 /e factor in dω/ dt. The two first equation are coupled, Eq. (C.5) will be solved in a second phase after injection of their solution. These equations are nonetheless irregular for small eccentricity regime. Thus, we will use the complex variable z = eexp(iω) to render them regular (Wyatt 2005). Moreover, from Eq. (C.4), we can deduce that the eccentricity is maximum when ω = ω_{p}. This information, combined with the initial value of the Hamiltonian, allows us to compute the maximal eccentricity as a function of a. These maximums prove themselves to be less than 0.5 in any case, so that we can linearize the system in z for an easier solving. It yields (C.6)where We now solve the system and get the eccentricity, precession and mean anomaly as a function of time. For null initial eccentricity, it is written If we represent the motion of z on the complex plane, we get exactly the circle depicted in Fig. 2 of Wyatt (2005). These formula were used to generate Fig. 9.
Appendix D: Density of stars around HD 106906
The first step to investigate the density of stars around the HD 106906 system is to build a complete list of known members in the LCC subgroup of the ScoCen association. Our list of LCC members is based on previous surveys of this region (De Zeeuw et al. 1999; Preibisch & Mamajek 2008; Song et al. 2012; Pecaut & Mamajek 2016) and consists of 369 stars. In the following, we estimate the current density of stars around the planetary system and its evolution in time. Thus, our methodology requires prior knowledge of the distances, proper motions and radial velocities for the individual stars in our sample.
The TychoGaia Astrometric Solution (TGAS, Lindegren et al. 2016) from the Gaia data release 1 provides trigonometric parallaxes and proper motions for 203 stars in our sample. To access more proper motion data, we also searched for this information in the PPMXL (Roeser et al. 2010), UCAC4 (Zacharias et al. 2012) and SPM4 (Girard et al. 2011) catalogs. Doing so, we find proper motions for 368 stars of the sample. We use the TGAS proper motions for the 203 stars and take the weighted mean of the multiple measurements given by the other catalogs (PPMXL, UCAC4 and SPM4) for the remaining 165 stars. Then, we searched the SIMBAD/CDS databases (Wenger et al. 2000) for radial velocity information using the data mining tools available on the site. The radial velocities that we use in this work come from Wilson (1953), Duflot et al. (1995), BarbierBrossat & Figon (2000), Torres et al. (2006), Gontcharov (2006), Holmberg et al. (2007), Mermilliod et al. (2009), Chen et al. (2011), Song et al. (2012), Kordopatis et al. (2013) and Desidera et al. (2015). We found radial velocity for 184 stars of our sample.
We apply the methodology developed by BailerJones (2015) to convert parallaxes into distances (see Sect. 7 of his paper). The systematic errors of about 0.3 mas in the TGAS parallaxes reported by Lindegren et al. (2016) were added quadratically to the parallax uncertainties. The threedimensional position of the stars are calculated from the individual distances in a XYZ grid where X points to the Galactic center, Y points in the direction of Galactic rotation, and Z points to the Galactic north pole. The reference system has its origin at the Sun. Then, we use the procedure described in Johnson & Soderblom (1987) to compute the UVW components of the spatial velocity for each star that are given in the same reference system. We perform a 3σ clipping in the distribution of proper motions, parallaxes, radial velocities and spatial velocities to remove obvious outliers. This procedure reduces the dataset to a total of 312 stars, but only 141 stars in this sample exhibit published radial velocities and 102 stars have complete data (proper motions, radial velocities and parallaxes). Based on this subset of 102 stars we calculate a revised mean spatial velocity of the LCC association of (U,V,W) = (−8.5,−21.1,−6.3) ± (0.2,0.2,0.2) km s^{1} (not corrected for the solar motion).
We note that 39 stars in the sample of 141 stars with known radial velocities do not have published parallaxes in the TGAS catalog. Individual parallaxes (and distances) can be inferred for these stars from the movingcluster method under the assumption that they are comoving. This method uses proper motions, radial velocities and the convergent point position of the moving group to derive individual parallaxes for group members (Galli et al. 2012). We emphasize that the soderived kinematic parallaxes are meaningful and provide valuable information in this work to increase the number of stars with measured parallax in our sample. We adopt the space motion listed above and the formalism described in Sect. 2 of Galli et al. (2017) to estimate the convergent point position and the kinematic parallaxes for each group member. Using a velocity dispersion of σ_{v} = 1.5 km s^{1} and distance estimate of 120 pc for the LCC association (see e.g., de Bruijne 1999) we find a convergent point solution located at (α_{cp},δ_{cp}) = (104.8°,−37.2°) ± (1.0°,0.8°) with chisquared statistics and correlation coefficient of ρ = −0.98. To gain confidence in the soderived kinematic parallaxes we compare our results with the trigonometric parallaxes from the TGAS catalog for the stars in common. We find a mean difference of 0.1 mas and rms of 0.6 mas, that are significantly smaller than the typical error on the kinematic parallaxes (~0.8 mas) derived from the movingcluster method in this analysis. This confirms the good agreement between the two datasets. Thus, the final sample with complete information (proper motion, radial velocity and parallax) that we use in this work to estimate the early density of stars around HD 106906 consists of 141 stars.
Fig. D.1
Evolution of the density of stars for different radii around the HD 106906 planetary system. The colored regions indicate the upper and lower limits (at the 1σ level) for the density of stars at a given radius. 
In a subsequent analysis, we consider the present day location of the 141 stars and use the UVW spatial velocity for each star to calculate their XYZ positions backward in time. We compute the stellar positions as a function of time in steps of 0.1 Myr from t = 0 (current position) to t = −14.0 Myr. The latter value is chosen to be consistent with an upper limit for the age estimate of the HD 106906 system as derived by Pecaut et al. (2012) from different evolutionary models. Then, we count the number of stars in the vicinity of HD 106906 for different radii (r = 5,10,15,...,30 pc) and determine the density of stars around the target. Figure D.1 illustrates the results of this investigation. Our analysis indicates that the early density of stars around HD 106906 (at t = −7.7 Myr) was higher than the current value by a factor of about 1.7 for r = 5 pc. At this stage it is important to mention that our result for density of stars is restricted to known members of the LCC association with complete data in our sample for which we can calculate spatial velocities and compute their positions back in time. As soon as new data (parallaxes and radial velocities) from the upcoming surveys (e.g., Gaia) become available and other group members are identified, a more refined analysis of this scenario will be made possible.
One alternative approach to better constrain the density of stars around HD 106906 consists of investigating the contribution of field stars (not related to the LCC association) in our solution. In this context, we use the model of stellar population synthesis from Robin et al. (2003) to simulate a catalog of pseudostars and their intrinsic properties (e.g., distances, spectral types, ages, magnitudes, etc.) in the direction of the HD 106906 system. We run the model with a distance range from 0 to 300 pc and a solid angle of 20 deg^{2} centered around the target. These values are chosen to include known members of the LCC association that is clearly spread in angular extent and exhibits significant depth effects along the line of sight. We do not constrain our simulations in magnitudes and spectral types to get a more complete picture of the stellar population in this region. We use a distance step of 0.5 pc in our simulations that is the minimum value that can be used in the model. The synthetic stars are all supposed to be at the same coordinates. So, we run a number of 1000 simulations to generate random coordinates for the simulated stars and use them (together with the distances provided by the model) to calculate the stellar threedimensional positions in the XYZ grid. Figure D.2 shows the density of stars around HD 106906 for different radii obtained from our sample of LCC stars, the pseudostars from our simulations and a combined result that includes both (cluster + field). Although this analysis cannot be extrapolated backward in time (as in Fig. D.1), it yields a more refined value for the current (t = 0) density of stars. However, we emphasize that the results obtained for small radii around the target (i.e., r ≤ 5 pc) are calculated with a small number of stars (typically, less than ten stars) and they should be regarded with caution. Thus, we conclude that the presentday density of stars in the vicinity of the HD 106906 system for r> 5 pc is ≤0.07 stars/pc^{3} (within the 1σ error bars). We infer from the results presented in Fig. D.1 an upper limit of ~0.11 stars/pc^{3} for the density of stars around HD 106906, a result that will need further confirmation as soon as more data becomes available.
Fig. D.2
Current density of stars in the vicinity of the HD 106906 system inferred from LCC cluster members and field stars. The colored region indicates the upper/lower limits (at the 1σ level) for the final density of stars (cluster + field) at different radii around the target. 
All Tables
Effect of the 1:6 meanmotion resonance and ejection duration for different eccentricities of the binary, starting with a planet eccentricity of e = 0.05.
First unstable resonance and corresponding ejection time for different eccentricities of the binary and different migration velocities, starting with e = 0.05 and a = 2 au.
Number of close encounters with a 1 M_{⊙} star raising the planet periastron above a given value (2, 50 or 150 au) depending on the trajectory of the planet before the flyby.
All Figures
Fig. 1
Chaotic zone (in dark gray) as a function of the binary eccentricity, for binary components of same masses. The lighter part designates a critical zone, where some test particles can survive. The red lines represent the lower and upper critical orbit parabolic fits found by Dvorak (1986) in its study of circumbinary planet stability. The 1:6 commensurability is the strongest outside the chaotic zone (see Sect. 2.2) for e_{B} ≥ 0.4. 

In the text 
Fig. 2
Isocontour in the (ν = ω−ω_{B}, e) phase space of the average interaction Hamiltonian of a test particle trapped in 1:6 meanmotion resonance with a binary eccentricity of e_{B} = 0.4, assuming a binary mass parameter of μ = 1/2. Each curve represents a trajectory in the (ν, e) space. Above the red line, the planet has part of its orbit in the chaotic zone. 

In the text 
Fig. 3
Evolution of the planet semimajor axis with respect to time for a binary eccentricity of e_{B} = 0.2 and a 10^{5} au/yr migration velocity. The semimajor axis of the binary has been set to a_{B} = 0.4 au. The plot illustrates the migration, then ejection, of the planet after it has been trapped into a 1:6 resonance. The effect of the 1:7 resonance, weaker, is also visible on the plot. As the planet has a perturbed Keplerian motion around the binary, the exact locations of MMRs are not straightforward to derive (see Appendix B). 

In the text 
Fig. 4
Example of a coplanar configuration where a passing star (in red) stabilizes a wide unstable planet orbit. Before the flyby, the planet orbit still has a very low periastron, and after it gets much wider, thanks to the interaction with the passing star. We recall that according to Kepler’s laws, the planet spends most of its time near apoastron, so that any flyby is likely to occur when the planet is at or near this point. 

In the text 
Fig. 5
Area in the disruptive star phase space which succeed to raise the periastron from 1 au to over 2 au in the case of a coplanar encounter of periastron argument ω_{∗}−ω = 45 degrees. v_{∗} designates the maximum relative velocity of the disruptive star, p_{∗} designates its smaller distance from the binary. 

In the text 
Fig. 6
Nbody simulation showing the consequence of the ejection of a 11 M_{J} planet through a disk. From left to right, the snapshots have been taken 0, 1000 and 10 000 years after the ejection. The planet starts on an hyperbolic orbit similar to what we observe in the simulations we performed: a = 10 au and e = 1.1 (corresponding to periastron q = 1 au). Above is a spatial representation of the top view of the disk (the planet trajectory is depicted in black), below is the density along the y axis, integrated over the x and z axis. 

In the text 
Fig. 7
Top view of the evolution of the debris disk after ten million years of perturbation by a planet on a coplanar orbit whose coplanar orbit has a periastron of 200 au and an apoastron of 1000 au. The color scale represents the relative density. Strong asymmetries can be seen. 

In the text 
Fig. 8
Density along the x axis, integrated over the y and z axis, obtained from Fig. 7. The gray zone marks approximately the cavity that we observe today. The asymmetry seems to reverse when we get farther to the stars. 

In the text 
Fig. 9
2D representation of a debris disk after 0, 10^{6} and 10^{7} years under the thirdorder approximation of the influence of the binary and a planet whose periastron is 200 au and apoastron is 1000 au. The color scale represents the local number of particles. At first, only the edge of the disk is affected, but after 10^{7} years, two spirals components appear within the disk. 

In the text 
Fig. 10
Precession velocities (above) and periastrons (below) of Fig. 9 test particles with respect to their semimajor axis (that does not change with time) after ten million years evolution. Above, the red curve displays the precession induced by the binary, the blue curve by the planet and the black curve depicts both contributions. 

In the text 
Fig. 11
Minimum mass of planets (in Jupiter masses) that can be detected into the H2 data published in Lagrange et al. (2016a) around HD 106906. The contrasts has been translated into masses using Baraffe et al. (2003) model adapted to the SPHERE filters. 

In the text 
Fig. 12
Planet mass ratios with respect to projected separation. Only planets that belong to young systems (<0.1 Gyr) are displayed, with the exception of the circumbinary planet Ross 458 c. HD 106906 b has the lowest mass ratio beyond 100 au. HIP 78530 b and HIP 77900 b are its two closest neighbors in the diagram. Being identified as brown dwarfs, they may have formed in situ by cloud collapse. Data retrieved from the exoplanets.eu database. 

In the text 
Fig. D.1
Evolution of the density of stars for different radii around the HD 106906 planetary system. The colored regions indicate the upper and lower limits (at the 1σ level) for the density of stars at a given radius. 

In the text 
Fig. D.2
Current density of stars in the vicinity of the HD 106906 system inferred from LCC cluster members and field stars. The colored region indicates the upper/lower limits (at the 1σ level) for the final density of stars (cluster + field) at different radii around the target. 

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.