A Near-coplanar Stellar Flyby of the Planet Host Star HD 106906
Robert J. De Rosa Paul Kalas




Robert J. De Rosa
rderosa@stanford.edu

0000-0002-4918-0247]Robert J. De Rosa
Kavli Institute for Particle Astrophysics and Cosmology, Stanford University, Stanford, CA 94305, USA
Department of Astronomy, University of California, Berkeley, CA 94720, USA

0000-0002-6221-5360]Paul Kalas
Department of Astronomy, University of California, Berkeley, CA 94720, USA
SETI Institute, 189 Bernardo Ave., Mountain View, 94043, USA


We present an investigation into the kinematics of HD 106906 using the newly released Gaia DR2 catalog to search for close encounters with other members of the Scorpius–Centaurus (Sco–Cen) association. HD 106906 is an eccentric spectroscopic binary that hosts both a large asymmetric debris disk extending out to at least 500 au and a directly imaged planetary-mass companion at a projected separation of 738 au. The cause of the asymmetry in the debris disk and the unusually wide separation of the planet is not currently known. Using a combination of Gaia DR2 astrometry and ground-based radial velocities, we explore the hypothesis that a close encounter with another cluster member within the last 15 Myr is responsible for the present configuration of the system. Out of 461 stars analyzed, we identified two candidate perturbers that had a median closest approach (CA) distance within 1 pc of HD 106906: HIP 59716 at D_ CA=0.65_-0.40^+0.93 pc (t_ CA=-3.49_-1.76^+0.90 Myr) and HIP 59721 at D_ CA=0.71_-0.11^+0.18 pc (t_ CA=-2.18_-1.04^+0.54 Myr), with the two stars likely forming a wide physical binary. The trajectories of both stars relative to HD 106906 are almost coplanar with the inner disk (Δθ = 54±1.7 and 42_-1.1^+0.9). These two stars are the best candidates of the currently known members of Sco–Cen for having a dynamically important close encounter with HD 106906, which may have stabilized the orbit of HD 106906 b in the scenario where the planet formed in the inner system and attained high eccentricity by interaction with the central binary.






§ INTRODUCTION


Close encounters by external stellar perturbers have the potential to significantly modify the architecture of planetary systems <cit.>. Stars passing very close to our solar system have been invoked to explain the formation of the Oort cloud of comets <cit.>, comet showers <cit.>, the disruption of the Kuiper Belt <cit.>, and the distant detached orbits of dwarf planets such as 90377 Sedna <cit.>, as well as the hypothetical Planet Nine <cit.>.  Such close stellar encounters have also been invoked to explain the orbital properties of extrasolar planets <cit.>, as well as the observed structure of disks around HD 100453 <cit.>, HD 141569 <cit.>,  HD 15115 <cit.>, and [bet Pic]β Pictoris <cit.>. 

In the case of β Pic, <cit.> and <cit.> explore numerical models of a non-coplanar stellar flyby (Δ i=30^∘) that qualitatively reproduce the disk asymmetries measured beyond the ∼200 au projected radius; the northeast disk extension is radially extended and vertically flat whereas the southwest disk extension is radially truncated and vertically extended <cit.>.  Using astrometry from the Hipparcos catalog, <cit.> identified several candidate perturbers that had a close encounter with the β Pic system in the last 1 Myr that could explain the asymmetry. The proximity and large apparent motion of β Pic were conducive for this analysis given the precision of the measurements available at the time. The order of magnitude improvement in the astrometric precision of the recently published Gaia DR2 catalog <cit.> makes it feasible to extend this analysis to more distant stars that also show evidence of dynamical perturbation.

In this paper we present an investigation into the dynamical history of HD 106906, an F5V <cit.> member of the 15 Myr-old <cit.> Lower Centaurus–Crux subgroup of the Scorpius–Centaurus OB association <cit.>. HD 106906 rose to prominence with the discovery of a planetary-mass companion (11±2 M_Jup) at a wide projected separation of 738 au <cit.>, exterior to a debris disk that was inferred from the spectral energy distribution of the star <cit.>. The inner disk was later resolved in near-infrared scattered light as a ∼50 au radius ring inclined to the line of sight by ∼85^∘ <cit.>.  However, visible light Hubble Space Telescope (HST) data showed a striking asymmetry in the morphology of the outer disk that resembles the disk asymmetry of β Pic <cit.>. The planet is oriented ∼21^∘ from the position angle of the disk midplane, suggesting that the planet's orbit is not coplanar with the disk <cit.>. HD 106906 is also a spectroscopic binary with a mass ratio near unity <cit.>. 

The hypothetical dynamical history of this system has been further studied by <cit.>, <cit.>, and <cit.>. <cit.> quantified a scenario where HD 106906 b originally formed in a disk near the binary star, migrated inward, encountered an unstable mean-motion resonance with the binary, was ejected into a high-eccentricity orbit, and then had its periastron raised into a stable region by an external stellar perturber. Here we search for potential stellar perturbers consistent with this scenario using the exquisite precision provided by the Gaia DR2 astrometry and ground-based radial velocities measurements.



§ SYSTEMIC VELOCITY OF HD 106906

HD 106906 consists of a 1.37 M_⊙ primary and 1.34 M_⊙ secondary with a projected separation of between 0.36–0.58 au <cit.>. The full orbital characterization of the system—most critically the systemic velocity—has not yet been published. There are four radial velocities for this system within the literature, two are instantaneous (10.2±1.7 km s^-1, ; 15.1±0.3 km s^-1, ), and two are from a combination of multiple measurements (8.4 km s^-1, N=22, ; 11.1 km s^-1, N=2, ). Given the scatter in these measurements, possibly caused by the orbit of the binary, we opted instead to use archival spectroscopic observations of HD 106906 to measure the systemic velocity.

Calibrated optical spectra of HD 106906 from the High Accuracy Radial velocity Planet Searcher (HARPS; ) instrument were obtained from the European Southern Observatory (ESO) Archive (program IDs 192.C-0224, 098.C-0739, 099.C-0205) for 76 epochs from 2014 January  to 2017 June. The radial velocities computed automatically by the HARPS Data Reduction Software were strongly biased by the blending of the rotationally broadened lines of both stars in all but four of the epochs, and were therefore not used in the subsequent analysis. Instead, we use the epoch in which the lines of the two stars are the most separated (2017 May 30) to fit a model atmosphere <cit.> to both stars which can then be used as a template to fit the Keplerian motion of both stars in the remaining epochs (Figure <ref>). The fit was restricted to seven lines/groups of lines centered on 5040, 5370, 5448, 5529, 6141, 6441, and 6680 Å. The best-fit model was found via χ^2 minimization by varying the temperature, surface gravity, radial velocity and rotational velocity of both stars, the flux ratio, and the metallicity. For each trial the model atmosphere was generated via linear interpolation, rotationally broadened, Doppler shifted, and smoothed to the resolution of HARPS. The HARPS and model atmosphere spectra were then continuum normalized with a linear fit to the continuum.


The resulting template for each star was used to fit the Doppler shift induced by the Keplerian motion of each star over the three-year baseline. We used the affine-invariant Markov chain Monte Carlo package  <cit.> to fit the seven orbital elements: period P, eccentricity e, argument of periapsis ω, epoch of periapsis τ, primary semi-amplitude K_1, mass ratio q, and the systemic velocity γ. At each step in the chain, the radial velocity of both stars at each epoch was predicted from the orbital elements. The templates constructed previously were Doppler shifted to the corresponding velocity, summed, and compared to the observed spectrum at that epoch. Mass functions that were significantly discrepant with the apparent brightness of the system were excluded. The period (P=49.233±0.001 days) and eccentricity (e=0.669±0.002) were well constrained, and are consistent with an independent analysis by <cit.>. We found a systemic velocity of γ=12.18±0.15 km s^-1 after marginalizing over the remaining parameters, the error being a combination of statistical error from the fit and an assumed systematic uncertainty to account for biases in the fitting process estimated by repeating the fitting process using subsets of the lines described previously.



§ CLOSE STELLAR ENCOUNTERS WITH HD 106906




 §.§ Identifying candidate perturbers

The revised systemic radial velocity was combined with astrometry from Gaia DR2 <cit.> to compute the position and kinematics of HD 106906 given in Table <ref>. Uncertainties were propagated in a Monte Carlo fashion, drawn from a multivariate normal distribution using the correlation coefficients given in the Gaia catalog for the astrometry, and from a normal distribution for the radial velocity. One potential source of error in the Gaia astrometry is the motion of the photocenter relative to the center of mass. All sources within the Gaia catalog were fit with a five-parameter single star model. The magnitude of this error should be small given the comparable mass and flux ratio of the two components.

A catalog of Sco–Cen members was constructed following <cit.>, augmented with members reported in <cit.>. Astrometric measurements were obtained from Gaia DR2, while radial velocities came from a number of sources, either from various literature sources <cit.>, or from the Gaia catalog. The final sample consisted of 461 members that had both a Gaia DR2 parallax and either a ground-based or Gaia radial velocity measurements. 

Candidate perturbers were identified through a linear approximation of the space motion of each of the sample stars relative to HD 106906. In this approximation, the velocity of each star is constant and gravitational forces are not considered. Uncertainties in the true position and motion of each star were propagated in a Monte Carlo fashion. A trial position, proper motion, and parallax were drawn from a multivariate normal distribution using the correlation coefficients in the Gaia catalog, and a radial velocity was drawn from a normal distribution. Closest approach times and distances (D_ CA) were computed for each of the 10^5 trials. This process was repeated for each star within the sample.



 §.§ HIP 59716 and HIP 59721

rcccc

Kinematics for HD 106906 and two candidate perturbers
5
1
0pt
Property     Unit     HD 106906     HIP 59716     HIP 59721
α     deg     184.47133075613± (1.48×10^-8)     183.71107188480±(1.31× 10^-8)     183.71768277794±(1.27× 10^-8) 

δ     deg     -55.97558054009± (7.41×10^-9)     -55.78991603504±(6.41× 10^-9)     -55.78439084834±(5.86× 10^-9) 

π     mas     9.6774±0.0429     8.7068±0.0393     8.6815±0.0349 

μ_αcosδ     mas yr^-1     -39.014±0.057     -35.722±0.050     -35.184±0.048 

μ_δ     mas yr^-1     -12.872±0.046     -11.041±0.042     -11.325±0.039

RV     km s^-1     12.18±0.15a     15.5±1.1b     17.6±1.7c

5cGalactic position and motion

X     pc     48.54±0.22     53.14±0.24     53.30±0.22 

Y     pc     -90.45±0.40     -100.93±0.46     -101.22±0.41 

Z     pc     11.86±0.05     13.43±0.06     13.48±0.05 

U     km s^-1     -9.77±0.10     -8.63±0.51     -7.42±0.79 

V     km s^-1     -20.12±0.14     -23.09±0.97     -24.83±1.49 

W     km s^-1     -7.29±0.05     -6.85±0.14     -6.74±0.20 

5cGalactocentric position and motion

X_G     pc     -8251.38±0.22     -8246.77±0.24     -8246.61±0.22 

Y_G     pc     -90.45±0.40     -100.93±0.46     -101.22±0.41 

Z_G     pc     38.70±0.05     40.26±0.06     40.31±0.05 

U_G     km s^-1     1.30±0.10     2.45±0.51     3.66±0.79 

V_G     km s^-1     212.12±0.14     209.15±0.97     207.41±1.49 

W_G     km s^-1     -0.01±0.05     0.43±0.14     0.54±0.20 

5cHD 106906 rest frame (-X^' toward Earth, Y^' toward East, Z^' toward North)

X^'     pc         11.52±0.69     11.85±0.65 

Y^'     pc         -0.8568±0.0039     -0.8519±0.0034 

Z^'     pc         0.3675±0.0017     0.3797±0.0015 

U^'     km s^-1         3.19±1.11     5.30±1.71 

V^'     km s^-1         -0.39±0.13     -0.16±0.13 

W^'     km s^-1         0.13±0.06     -0.03±0.06 

5cHD 106906 debris disk frame (disk lies in Y_ diskZ_ disk plane)

X_ disk     pc                   1.136±0.061     1.179±0.057 

Y_ disk     pc                   -0.9227±0.0042     -0.9212±0.0038 

Z_ disk     pc                   -11.46±0.69     -11.79±0.64 

U_ disk     km s^-1         0.30±0.13     0.39±0.17 

V_ disk     km s^-1         -0.41±0.12     -0.15±0.11 

W_ disk     km s^-1         -3.18±1.10     -5.29±1.70 

aThis work
b<cit.>
c<cit.>






Two candidate perturbers, likely forming a wide physical binary, were identified as having a closest approach distance within 1 pc of HD 106906 within the last 15 Myr: HIP 59716 (F5V, ) at t_ CA=-3.49_-1.76^+0.90 Myr, D_ CA=0.65_-0.40^+0.93 pc; and HIP 59721 (G9V, ) at t_ CA=-2.18_-1.04^+0.54 Myr, D_ CA=0.71_-0.11^+0.18 pc. The closest approach times and distances are different for the two stars despite them likely forming a wide binary due to the large uncertainty on the relative radial positions and velocities, and the fact that this calculation does not incorporate gravitational forces (see Section <ref>). We estimate masses of 1.37 and 1.22 M_⊙ for the two stars based on their spectral type. The closest approach times and distances for the two stars is shown in Figure <ref>. These distances are comparable to the closest approach of WISE J072003.20-084651.2 to the Sun 70  kya at 0.25^+0.11_-0.07  pc <cit.>, although the solar system is significantly more compact than the HD 106906 system. 

HIP 59716 and HIP 59721 are presently separated from each other by 24 on the sky and are listed in the Washington Double Star Catalog <cit.> as being a physical binary based on the similarity of their parallax. The Gaia DR2 astrometry is consistent with this assessment; the parallax of the two stars are indistinguishable (Δπ=0.025±0.053 mas), and the proper motion difference (0.61±0.07 mas yr^-1) is less than the circular motion of a face-on 2,750 au orbit (∼1.7 mas yr^-1). The similar kinematics of this group of stars had been previously noted. Using Hipparcos astrometry, <cit.> estimated a probability of 92 % that HD 106906 and HIP 59716 were bound. No assessment was made regarding the nature of HIP 59721 due to the poor quality of the parallax measurement in the Hipparcos catalog (π=7.56±5.84 mas, ).

While the Gaia astrometry provides exquisite constraints on the position and motion of the two stars in the sky plane (projected separation to within 0.006 au, relative motion to within 0.008 au yr^-1), the position and motion in the radial direction is less well constrained; the relative distance is known to within 140,000 au, and the relative radial velocity to within 0.3 au yr^-1. This precision is insufficient to confirm that the two stars are bound. Instead, we make a probabilistic argument based on the chance of finding two stars in LCC with a projected separation within 24 that share a common parallax and proper motion. Using a simplistic assumption that the Galactic positions and motions of LCC members within our sample follow normal distributions, we simulated 10^8 stars and compared their positions, parallax, and proper motions in the sky plane to that of HIP 59716. Only 20 of the 10^8 stars were within 24, had a parallax within 0.078 mas, and had total proper motion within 0.68 mas yr^-1, implying a chance alignment probability of ∼10^-7.



 §.§ N-body simulations


The proximity of HD 106906 and HIP 59716 at closest approach, and between HIP 59716 and HIP 59721 at the present epoch, motivated us to investigate the effect of  the gravitational interaction between the stars and the influence of the Galactic potential. We used the N-body  package <cit.> to predict the position of the three stars within a Galactic potential <cit.> over the previous 10 Myr. We initialized 10^4 simulations to sample the uncertainties on the astrometry and radial velocity for each star. Each simulation was integrated backwards in time, with the absolute positions and velocities recorded every 100 years. The closest approach times and distances were unchanged from the linear approximation t_ CA=-3.5_-1.8^+0.9 Myr, D_ CA=0.63_-0.39^+0.92 pc for HIP 59716 and t_ CA=-2.2_-1.0^+0.6 Myr, D_ CA=0.62_-0.17^+0.21 pc for HIP 59721 (marginalized distributions shown in Figure <ref>).

The consistency between the simulations with and without gravity was expected. The difference in the Galactic potential experienced by the three stars over a small fraction of their Galactic orbital period is negligible, and the large uncertainty on the relative radial velocity and distance between HIP 59716 and HIP 59721 meant that the pair were only bound in ∼2 % of the simulations. The closest approach distances were not significantly different for the bound simulations, increasing slightly for HIP 59716 to D_ CA=0.42_-0.17^+0.30 pc, and decreasing slightly for HIP 59721 to D_ CA=0.42_-0.18^+0.30 pc. We note that these values are derived from a small number of simulations (212 of 10,000), and may not constitute a representative sample of bound orbit trajectories. The projected separation between the two stars of 2,750 au suggests a comparable semi-major axis with the a/ρ conversion factor peaking at unity for a uniform eccentricity distribution <cit.>. Only in the most extreme cases where a/ρ∼ 3 and e∼ 1 does the apocenter distance, where eccentric binaries spend more of their time, become comparable to the closest approach distance. For the purposes of this study, we assume that the dynamics of the binary system—if it is indeed bound—does not have a large effect on the flyby of HD 106906.



 §.§ Flyby geometry

The position and motion for the three stars in the Galactic and Galactocentric coordinate systems are given in Table <ref>. We defined two new coordinate systems with HD 106906 at the origin (in terms of both position and velocity). The first has -X^' pointing toward the current location of the Sun and Y^' and Z^' pointing toward East and North from the perspective of the Sun. In this coordinate system, the plane of the disk is rotated 15^∘ about the X^' axis and 85^∘ about the Y^' axis. The second is rotated  about the first such that the disk lies in one of the planes of the coordinate system (Y_ diskZ_ disk). A schematic diagram is shown in Figure <ref>. 

The trajectories of both stars were almost coplanar with the current plane of the debris disk: Δθ = 54±17 (52_-1.8^+1.5) for HIP 59716 and Δθ = 42_-1.1^+0.9 (48±10) for HIP 59721 (Fig. <ref>, second column; values in parentheses were calculated from the N-body simulation). The relative inclination between the flyby and the orbit of the planet cannot be constrained due to the lack of multi-epoch astrometry to measure the planet’s orbit. We calculated a relative velocity of 3.3±1.1 km s^-1 and 5.3±1.7 km s^-1 at the time of closest approach in the N-body simulation for HIP 59716 and HIP 59721, respectively. The velocity vector is almost entirely in the disk plane, consistent with a coplanar encounter.



§ EFFECT OF MEASUREMENT UNCERTAINTIES



The relative three-dimensional motion between HD 106906 and the two candidate perturbers is almost entirely in the radial direction, and as such our uncertainty on the closest approach distance is almost entirely due to the uncertainty in the relative radial positions and motions of the three stars. While the uncertainty in the relative radial positions of the three stars is a factor of 10^8 larger than the relative tangential positions, the dominant source of uncertainty on the predicted closest approach distance is the relative radial velocity as this error accumulates as the positions of each star are traced back in time. 

To explore the effects of varying the relative radial velocities on the closest approach differences, we repeated the traceback analysis for HIP 59716 and HIP 59721 for a range of plausible radial velocities. For each combination of radial velocities, we drew 10^4 random variates from the multivariate normal distribution that describe the Gaia astrometry and their correlation for each star and determined D_ CA and t_ CA using the linear approximation described previously. We assume that the radial velocities are perfectly measured. The correlations between D_ CA and t_ CA and the radial velocities of the three stars are shown in Figure <ref>. D_ CA and t_ CA are only sensitive to the relative radial velocity between HD 106906 and the two other stars rather than their absolute velocities; lines of constant D_ CA and t_ CA follow lines of constant relative velocity. Increasing the relative velocity between HD 106906 and HIP 59716 causes a significant decrease in D_ CA and an increase in t_ CA (a more recent flyby). Interestingly, if HIP 59716’s radial velocity was similar to that of HIP 59721, the median value of D_ CA would be near the global minimum found in this grid search. HIP 59716 and HIP 59721 both have relatively large uncertainties on their radial velocities, and periodic measurements will be required in order to rule out the presence of binary companion to either of the two stars.



§ MAGNITUDE OF DYNAMICAL INTERACTION

There are two peculiar features of the HD 106906 system that may be explained by the dynamical perturbation caused by a passing star: the asymmetric debris disk <cit.> and the planet at a large projected separation <cit.>. The debris disk has two components: an inner ring at a radius of 50 au that exhibits a strong brightness asymmetry, and a sharp feature resolved with HST extending at least 500 au westward along the disk plane with no counterpart seen on the eastern side <cit.>. Several theories have been suggested to explain the observed asymmetry including perturbation by the planet on its current orbit <cit.>, and a scattering event where the planet interacted with either an unseen low-mass companion <cit.> or the inner binary <cit.> and then migrated through the disk. The presence of the planet at a projected separation of 738 au is also a mystery. Canonical planet formation theories suggest that massive planets cannot form at such wide separations. Instead, the planet must have either formed closer in and later migrated or was scattered to its current separation, or formed more like a binary star in the initial fragmentation of the molecular cloud as HD 106906 itself was being formed.

The flyby of HIP 59716 and HIP 59721 may have been important in the context of both the observed asymmetry in the debris disk and the current wide separation of the planet. The two stars passed well within the tidal radius of the HD 106906 binary (1.9 pc; ) interior to which the gravitational interaction between the stars is stronger than the Galactic potential. In a recent study, <cit.> postulated that the present orbit of the planet could be explained by a combination of dynamical interaction with an eccentric inner binary (we find e∼0.67 based on our orbit fit discussed previously) and a close encounter with a passing star.

In this scenario, the planet formed at a location within the circumbinary disk that was stable to the dynamical effects of the binary. After formation, the planet migrated inward until it was temporarily trapped in a 1:6 mean-motion resonance with the inner binary. Dynamical interaction with the inner binary cause the semi-major axis of the planet to rapidly oscillate before being ejected from the inner system. <cit.> suggests that the very eccentric (or hyperbolic) orbit of the ejected planet could then be stabilized by an interaction with a passing star that decreases the eccentricity of the orbit and raises the periapsis distance out of the chaotic zone surrounding the binary. Without a stabilizing encounter the planet would enter the chaotic zone surrounding the inner binary each time it goes through periapsis, rapidly leading to the complete ejection of the planet from the system.

The configuration of the flyby required to increase the periapsis distance out of the chaotic zone depends on the mass and closest approach distance of the passing star, the inclination of the encounter, and the angle between the velocity vector of the planet and passing star. Such a scenario has also been suggested for the hypothesized ninth planet within our own solar system <cit.>, with the planet being ejected from the inner solar system due to interactions with the other gas giants, and then stabilized by the gravitational influence of a passing star <cit.>.

In the scenario described by <cit.>, the planet is rapidly ejected from the inner system after reaching the 1:6 mean-motion resonance with the inner eccentric binary. The window for a stabilizing timescale is short (<10^4 yr); either the planet is ejected immediately on a hyperbolic trajectory, or it first achieves a highly eccentric orbit that evolves into a hyperbolic trajectory after a few passes through the chaotic zone surrounding the binary. Given a mean age of 15 Myr for LCC members, and an intrinsic scatter of 6 Myr <cit.>, the plausible ages for HD 106906 range from 9–21 Myr. Assuming a younger age of ∼10 Myr, the formation of the planet and migration into mean-motion resonance with the inner binary has to occur in ∼5 Myr for the planet to be ejected as the perturbing stars were passing by. These timescales are similar to those proposed for planet formation via pebble accretion (e.g., ), although formation via gravitational instability would occur much more rapidly (e.g., ). 

For the older age estimate of ∼20 Myr, the planet would need to spend over ten million years outside of the chaotic zone surrounding the binary before entering the mean-motion resonance approximately four million years ago. These timings are inconsistent with the requirement of a massive circumstellar disk to enable planetary migration; disks around the majority of young stars dissipate within a few million years (e.g., ). In this case, if a stabilizing encounter did occur, it is more likely to have been with another cluster member in the more distant past, rather than with HIP 59716 and HIP 59721. Future work could search for a dynamical process where the ejection of HD 106906 b occurs at a later epoch than calculated by <cit.>.  



 §.§ Flyby simulations


We ran 5 × 10^4  N-body simulations to explore the dynamical impact of a stellar flyby on a planet with a highly eccentric (e>0.9) orbit. For these simulations we assumed that HIP 59716 and HIP 59721 form a physical binary, and treated them as a single particle of their combined mass given current uncertainties about their orbit. This is a reasonable approximation when the flyby distance is significantly larger than the binary separation, where the gravitational influence of the two stars would not be resolved by the orbiting planet, but may not be valid in cases where the minimum separation between the planet and perturbing stars is comparable to their orbital semi-major axis. We note that these simulations do not model the ejection mechanism itself, rather the planet is initialized on a highly eccentric orbit to simulate the orbital configuration after an ejection has occurred. A more detailed modeling effort will be needed to combine an ejection event with a stellar flyby into a single simulation.

Each simulation was initialized with a 2.71 M_⊙ particle at (0, 0) pc, accounting for the combined mass of the HD 106906 binary, around which a 11 M_ Jup particle was placed with an orbital eccentricity of log_10(1-e) drawn from 𝒰(-4, -1) (where 𝒰 denotes the uniform distribution), a semi-major axis with a corresponding periapsis distance (r_ peri) of 1.5 au, and a mean anomaly M drawn from 𝒰(0, 2π) rad so that the flyby occurs at a different phase of the orbit within each simulation. The two candidate perturbers were modeled by a single particle with a mass of 2.59 M_⊙, initialized at y = -3 pc, log_10(x) drawn from 𝒰(-2, 0.5) pc, with a velocity of dy/dt = 3 km s^-1. In half of the simulations the perturber was instead initialized at y = 3 pc, with a velocity of dy/dt = -3 km s^-1, to account for the unknown direction of the orbit of the planet. A coplanar encounter was assumed to limit the size of phase space to explore. Each simulation was advanced either until the perturbers had reached y = 3 pc (or -3 pc), or for 2 Myr, whichever occurred first.

Of the 25,000 prograde (retrograde) flybys; 21,670 (6,385) caused an increase in r_ peri and 3,330 (18,615) led to a decrease in r_ peri, with 1,541 (1,525) resulting in the planet being ejected from the system. The correlation between the change in periapsis distance (Δ r_ peri) and the closest approach distance between the two stars for the simulations in which the eccentricity decreased is shown in Figure <ref>. There is a clear correlation between Δ r_ peri and D_ CA, with closer encounters leading to a larger Δ r_ peri for a given initial eccentricity. The presence of a debris ring at ∼50 au around the HD 106906 binary suggests that the r_ peri for the planet is now >50 au, otherwise the disk would be disrupted with each periastron passage for relative inclinations ≲10 deg <cit.>. Periapsis distances interior to the radius of the debris ring are possible for higher relative inclinations (≳40 deg) without disrupting the disk. These N-body simulations suggest that r_ peri for a highly eccentric planet can be plausibly raised from ∼1.5 to ∼100 au with a co-planar stellar flyby within 0.2 pc. Future observational and theoretical work is needed to more precisely constrain the geometry of the encounter and its dynamical consequences.



§ CONCLUSIONS

HIP 59716 and HIP 59721 are the best candidates of the currently known members of Sco–Cen for a dynamically important close encounter with HD 106906 within the last 15 Myr. The flyby of these two stars fulfill many of the criteria for the stabilization scenario described in <cit.>. Their trajectories are almost coplanar with the debris disk in its current orientation, their velocities relative to HD 106906 at closest approach are low (the change in velocity of the orbiting planet being inversely proportional to the relative velocity of the passing star at closest approach), and the distribution of closest approach distances for HIP 59716 is consistent with a dynamically significant encounter within 0.5 pc.

The biggest source of uncertainty in the closest approach distances are the relative radial velocities. Future spectroscopic observations of the two candidate perturbers will be essential to precisely determine their radial velocities. For example, increasing the relative velocity between HD 106906 and HIP 59716 by 1 km s^-1 significantly increases the probability of a closest approach distance within 0.1 pc (Fig. <ref>). The astrometry for each star will also be improved with  upcoming Gaia data releases. Not only will the precision of the measurements improve, but the photocenter motion of short period binaries will also be accounted for. With these data in hand, a more precise determination of the kinematics of the three stars can be made.

We wish to thank Eric Nielsen, Ian Czekala, Ruth Murray-Clay, and Anne-Marie Lagrange  for useful discussions relating to this work. We also wish to thank the referee for the comments that helped improve the quality of this work. The authors were supported in part by NSF AST-1518332, NASA NNX15AC89G and NNX15AD95G. This work benefited from NASA's Nexus for Exoplanet System Science (NExSS) research coordination network sponsored by NASA’s Science Mission Directorate. This work has made use of data from the European Space Agency (ESA) mission Gaia (<https://www.cosmos.esa.int/gaia>), processed by the Gaia Data Processing and Analysis Consortium (DPAC; <https://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 and the VizieR catalog access tool, both operated at the CDS, Strasbourg, France.

Astropy <cit.>,  
          Matplotlib <cit.>,
          galpy <cit.>,
          REBOUND <cit.>



