A Near-coplanar Stellar Flyby of the Planet Host Star HD 106906

Close encounters by external stellar perturbers have the potential to significantly modify the architecture of planetary systems (Laughlin & Adams 1998; Kenyon & Bromley 2004). Stars passing very close to our solar system have been invoked to explain the formation of the Oort cloud of comets (Duncan et al. 1987), comet showers (Hills 1981), the disruption of the Kuiper Belt (Kobayashi et al. 2005), and the distant detached orbits of dwarf planets such as 90377 Sedna (Brown et al. 2004), as well as the hypothetical Planet Nine (Bromley & Kenyon 2016). Such close stellar encounters have also been invoked to explain the orbital properties of extrasolar planets (Zakamska & Tremaine 2004; Malmberg & Davies 2009; Parker & Quanz 2011; Pfalzner et al. 2018), as well as the observed structure of disks around HD 100453 (Wagner et al. 2018), HD 141569 (Ardila et al. 2005; Reche et al. 2008), HD 15115 (Kalas et al. 2007), and β Pictoris (Kalas & Jewitt 1995; Ballering et al. 2016).

In the case of β Pic, Kalas et al. (2000) and Larwood & Kalas (2001) 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 (Kalas & Jewitt 1995). Using astrometry from the Hipparcos catalog, Kalas et al. (2001) 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 (Gaia Collaboration et al. 2018) 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 (Houk & Cowley 1975) member of the 15 Myr old (Pecaut & Mamajek 2016) Lower Centaurus–Crux (LCC) subgroup of the Scorpius–Centaurus OB association (Chen et al. 2011). 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 (Bailey et al. 2014) exterior to a debris disk that was inferred from the spectral energy distribution of the star (Chen et al. 2005). 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° (Kalas et al. 2015; Lagrange et al. 2016). 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 (Kalas et al. 2015). 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 (Kalas et al. 2015). HD 106906 is also a spectroscopic binary with a mass ratio near unity (Lagrange et al. 2019).

The hypothetical dynamical history of this system has been further studied by Jílková & Zwart (2015), Nesvold et al. (2017), and Rodet et al. (2017). Rodet et al. (2017) 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.

HD 106906 consists of a 1.37 M_⊙ primary and 1.34 M_⊙ secondary with a projected separation of between 0.36 and 0.58 au (Rodet et al. 2017). 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⁻¹, Gontcharov 2006; 15.1 ± 0.3 km s⁻¹, Chen et al. 2011), and two are from a combination of multiple measurements (8.4 km s⁻¹, N = 22, Evans 1967; 11.1 km s⁻¹, N = 2, Nordström et al. 2004). 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; Mayor et al. 2003) 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 (Husser et al. 2013) to both stars which can then be used as a template to fit the Keplerian motion of both stars in the remaining epochs (Figure 1). 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 χ² 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.

Zoom In Zoom Out Reset image size

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 emcee (Foreman-Mackey et al. 2013) to fit the seven orbital elements: period P, eccentricity e, argument of periapsis ω, epoch of periapsis τ, primary semi-amplitude K₁, 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 Lagrange et al. (2019). We found a systemic velocity of γ = 12.18 ± 0.15 km s⁻¹ 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.

3.1. Identifying Candidate Perturbers

3.2. HIP 59716 and HIP 59721

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, Houk & Cowley 1975) at ${t}_{\mathrm{CA}}=-{3.49}_{-1.76}^{+0.90}$ Myr, ${D}_{\mathrm{CA}}={0.65}_{-0.40}^{+0.93}$ pc and HIP 59721 (G9V, Torres et al. 2006) at ${t}_{\mathrm{CA}}=-{2.18}_{-1.04}^{+0.54}$ Myr, ${D}_{\mathrm{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 3.3). 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 2. These distances are comparable to the closest approach of WISE J072003.20-084651.2 to the Sun 70 kya at ${0.25}_{-0.07}^{+0.11}$ pc (Mamajek et al. 2015), although the solar system is significantly more compact than the HD 106906 system.

Zoom In Zoom Out Reset image size

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 (Mason et al. 2001) 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⁻¹) is less than the circular motion of a face-on 2750 au orbit (∼1.7 mas yr⁻¹). The similar kinematics of this group of stars had been previously noted. Using Hipparcos astrometry, Shaya & Olling (2010) 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; van Leeuwen 2007).

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⁻¹), 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⁻¹. 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⁸ stars and compared their positions, parallax, and proper motions in the sky plane to that of HIP 59716. Only 20 of the 10⁸ stars were within 24'', had a parallax within 0.078 mas, and had total proper motion within 0.68 mas yr⁻¹, implying a chance alignment probability of ∼10⁻⁷.

3.3. 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 REBOUND package (Rein & Liu 2012; Rein & Spiegel 2015) to predict the position of the three stars within a Galactic potential (Bovy 2015) over the previous 10 Myr. We initialized 10⁴ 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 yr. The closest approach times and distances were unchanged from the linear approximation ${t}_{\mathrm{CA}}=-{3.5}_{-1.8}^{+0.9}$ Myr, ${D}_{\mathrm{CA}}={0.63}_{-0.39}^{+0.92}$ pc for HIP 59716 and ${t}_{\mathrm{CA}}=-{2.2}_{-1.0}^{+0.6}$ Myr, ${D}_{\mathrm{CA}}={0.62}_{-0.17}^{+0.21}$ pc for HIP 59721 (marginalized distributions shown in Figure 2).

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}_{\mathrm{CA}}={0.42}_{-0.17}^{+0.30}$ pc, and decreasing slightly for HIP 59721 to ${D}_{\mathrm{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 2750 au suggests a comparable semimajor axis with the a/ρ conversion factor peaking at unity for a uniform eccentricity distribution (Dupuy & Liu 2011). 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.

Zoom In Zoom Out Reset image size

3.4. Flyby Geometry

The position and motion for the three stars in the Galactic and Galactocentric coordinate systems are given in Table 1. 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). The distribution of closest approach locations in the second coordinate system is shown in Figure 3, and a schematic diagram of the flyby is shown in Figure 4.

Zoom In Zoom Out Reset image size

The trajectories of both stars were almost coplanar with the current plane of the debris disk: Δθ = 54 ± 17 ($5\buildrel{\circ}\over{.} {2}_{-1.8}^{+1.5}$) for HIP 59716 and ${\rm{\Delta }}\theta =4\buildrel{\circ}\over{.} {2}_{-1.1}^{+0.9}$ (48 ± 10) for HIP 59721 (Figure 1, second column; the 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⁻¹ and 5.3 ± 1.7 km s⁻¹ 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.

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⁸ 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 distances, 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⁴ 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 5. 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.

Zoom In Zoom Out Reset image size

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 (Kalas et al. 2015; Lagrange et al. 2016) and the planet at a large projected separation (Bailey et al. 2014). 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 (Kalas et al. 2015). Several theories have been suggested to explain the observed asymmetry including perturbation by the planet on its current orbit (Jílková & Zwart 2015; Nesvold et al. 2017) and a scattering event where the planet interacted with either an unseen low-mass companion (Kalas et al. 2015) or the inner binary (Rodet et al. 2017) 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; Jiang & Tremaine 2010; Mamajek et al. 2013) interior to which the gravitational interaction between the stars is stronger than the Galactic potential. In a recent study, Rodet et al. (2017) 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 semimajor axis of the planet to rapidly oscillate before being ejected from the inner system. Rodet et al. (2017) 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 (Batygin & Brown 2016), 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 (Bromley & Kenyon 2016).

In the scenario described by Rodet et al. (2017), 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⁴ 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 (Pecaut & Mamajek 2016), the plausible ages for HD 106906 range from 9 to 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., Lambrechts & Johansen 2014), although formation via gravitational instability would occur more rapidly (e.g., Boss 2011).

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., Ansdell et al. 2017). 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 Rodet et al. (2017).

5.1. Flyby Simulations

We ran 5 × 10⁴ REBOUND 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 semimajor 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 ${\mathrm{log}}_{10}(1-e)$ drawn from ${ \mathcal U }(-4,-1)$ (where ${ \mathcal U }$ denotes the uniform distribution), a semimajor axis with a corresponding periapsis distance (r_peri) of 1.5 au, and a mean anomaly M drawn from ${ \mathcal U }(0,2\pi )$ 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, ${\mathrm{log}}_{10}(x)$ drawn from ${ \mathcal U }(-2,0.5)$ pc, with a velocity of ${dy}/{dt}=3$ km s⁻¹. In half of the simulations the perturber was instead initialized at y = 3 pc, with a velocity of dy/dt = −3 km s⁻¹, 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 (6385) caused an increase in ${r}_{\mathrm{peri}}$ and 3330 (18,615) led to a decrease in r_peri with 1541 (1525) 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 6. 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° (Jílková & Zwart 2015). Periapsis distances interior to the radius of the debris ring are possible for higher relative inclinations (≳40°) 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 coplanar 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.

Zoom In Zoom Out Reset image size

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 Rodet et al. (2017). 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⁻¹ significantly increases the probability of a closest approach distance within 0.1 pc (Figure 5). 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.

Software: Astropy (The Astropy Collaboration et al. 2013), Matplotlib (Hunter 2007), galpy (Bovy 2015), REBOUND (Rein & Liu 2012).