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

TL;DR

This study investigates whether a recent close stellar encounter influenced the unusual debris disk and wide-orbit planet of HD 106906 by analyzing Gaia DR2 data and identifying potential perturbers within the Sco-Cen association.

Contribution

It presents the first detailed kinematic analysis linking a near-coplanar stellar flyby to the peculiar features of HD 106906, suggesting a possible dynamical influence on its planetary system.

Findings

Identified two candidate stars that had close approaches within 1 pc of HD 106906.

Both candidate stars likely formed a wide binary and had trajectories coplanar with the disk.

Close encounters may have stabilized the orbit of HD 106906 b.

Abstract

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_{\rm CA}=0.65_{-0.40}^{+0.93}$ pc ($t_{\rm CA}=-3.49_{-1.76}^{+0.90}$ Myr) and HIP 59721 at $D_{\rm CA}=0.71_{-0.11}^{+0.18}$ pc ($t_{\rm 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 ($\Delta\theta = 5.4\pm1.7$ deg and $4.2_{-1.1}^{+0.9}$ deg). 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.

Tables

Table 1. Table 1: Kinematics for HD 106906 and two candidate perturbers

Property	Unit	HD 106906	HIP 59716	HIP 59721
$\alpha$	deg	$184.47133075613\pm \left(1.48\times {10}^{-8}\right)$	$183.71107188480\pm \left(1.31\times {10}^{-8}\right)$	$183.71768277794\pm \left(1.27\times {10}^{-8}\right)$
$\delta$	deg	$-55.97558054009\pm \left(7.41\times {10}^{-9}\right)$	$-55.78991603504\pm \left(6.41\times {10}^{-9}\right)$	$-55.78439084834\pm \left(5.86\times {10}^{-9}\right)$
$\pi$	mas	$9.6774\pm 0.0429$	$8.7068\pm 0.0393$	$8.6815\pm 0.0349$
${\mu }_{\alpha }\cos \delta$	mas yr-1	$-39.014\pm 0.057$	$-35.722\pm 0.050$	$-35.184\pm 0.048$
${\mu }_{\delta }$	mas yr-1	$-12.872\pm 0.046$	$-11.041\pm 0.042$	$-11.325\pm 0.039$
RV	km s-1	$12.18\pm 0.15$aaThis work	$15.5\pm 1.1$bbChen et al. (2011)	$17.6\pm 1.7$ccSong et al. (2012)
Galactic position and motion
$X$	pc	$48.54\pm 0.22$	$53.14\pm 0.24$	$53.30\pm 0.22$
$Y$	pc	$-90.45\pm 0.40$	$-100.93\pm 0.46$	$-101.22\pm 0.41$
$Z$	pc	$11.86\pm 0.05$	$13.43\pm 0.06$	$13.48\pm 0.05$
$U$	km s-1	$-9.77\pm 0.10$	$-8.63\pm 0.51$	$-7.42\pm 0.79$
$V$	km s-1	$-20.12\pm 0.14$	$-23.09\pm 0.97$	$-24.83\pm 1.49$
$W$	km s-1	$-7.29\pm 0.05$	$-6.85\pm 0.14$	$-6.74\pm 0.20$
Galactocentric position and motion
$X_G$	pc	$-8251.38\pm 0.22$	$-8246.77\pm 0.24$	$-8246.61\pm 0.22$
$Y_G$	pc	$-90.45\pm 0.40$	$-100.93\pm 0.46$	$-101.22\pm 0.41$
$Z_G$	pc	$38.70\pm 0.05$	$40.26\pm 0.06$	$40.31\pm 0.05$
$U_G$	km s-1	$1.30\pm 0.10$	$2.45\pm 0.51$	$3.66\pm 0.79$
$V_G$	km s-1	$212.12\pm 0.14$	$209.15\pm 0.97$	$207.41\pm 1.49$
$W_G$	km s-1	$-0.01\pm 0.05$	$0.43\pm 0.14$	$0.54\pm 0.20$
HD 106906 rest frame ($-X^{\prime }$ toward Earth, $Y^{\prime }$ toward East, $Z^{\prime }$ toward North)
$X^{\prime }$	pc	$\cdots$	$11.52\pm 0.69$	$11.85\pm 0.65$
$Y^{\prime }$	pc	$\cdots$	$-0.8568\pm 0.0039$	$-0.8519\pm 0.0034$
$Z^{\prime }$	pc	$\cdots$	$0.3675\pm 0.0017$	$0.3797\pm 0.0015$
$U^{\prime }$	km s-1	$\cdots$	$3.19\pm 1.11$	$5.30\pm 1.71$
$V^{\prime }$	km s-1	$\cdots$	$-0.39\pm 0.13$	$-0.16\pm 0.13$
$W^{\prime }$	km s-1	$\cdots$	$0.13\pm 0.06$	$-0.03\pm 0.06$
HD 106906 debris disk frame (disk lies in $Y_{\mathrm{disk}}{Z}_{\mathrm{disk}}$ plane)
$X_{\mathrm{disk}}$	pc	$\cdots$	$1.136\pm 0.061$	$1.179\pm 0.057$
$Y_{\mathrm{disk}}$	pc	$\cdots$	$-0.9227\pm 0.0042$	$-0.9212\pm 0.0038$
$Z_{\mathrm{disk}}$	pc	$\cdots$	$-11.46\pm 0.69$	$-11.79\pm 0.64$
$U_{\mathrm{disk}}$	km s-1	$\cdots$	$0.30\pm 0.13$	$0.39\pm 0.17$
$V_{\mathrm{disk}}$	km s-1	$\cdots$	$-0.41\pm 0.12$	$-0.15\pm 0.11$
$W_{\mathrm{disk}}$	km s-1	$\cdots$	$-3.18\pm 1.10$	$-5.29\pm 1.70$

Full text

astrometry — planets and satellites: dynamical evolution and stability — stars: kinematics and dynamics — techniques: radial velocities — stars: individual (HD 106906, HIP 59716, HIP 59721)

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

1 Introduction

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 (catalog ) (Wagner et al., 2018), HD 141569 (catalog ) (Ardila et al., 2005; Reche et al., 2008), HD 15115 (catalog ) (Kalas et al., 2007), and $\beta$ Pictoris (catalog bet Pic) (Kalas & Jewitt, 1995; Ballering et al., 2016).

In the case of $\beta$ Pic, Kalas et al. (2000) and Larwood & Kalas (2001) explore numerical models of a non-coplanar stellar flyby ($\Delta i=30^{\circ}$) that qualitatively reproduce the disk asymmetries measured beyond the $\sim$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 $\beta$ Pic system in the last 1 Myr that could explain the asymmetry. The proximity and large apparent motion of $\beta$ 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 (catalog ), an F5V (Houk & Cowley, 1975) member of the 15 Myr-old (Pecaut & Mamajek, 2016) Lower Centaurus–Crux 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\pm 2$ MJup) 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 $\sim$50 au radius ring inclined to the line of sight by $\sim$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 $\beta$ Pic (Kalas et al., 2015). The planet is oriented $\sim$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 (catalog ) 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.

2 Systemic Velocity of HD 106906

HD 106906 consists of a $1.37\,M_{\odot}$ primary and $1.34\,M_{\odot}$ secondary with a projected separation of between 0.36–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\pm 1.7$ km s*-1*, Gontcharov, 2006; $15.1\pm 0.3$ km s*-1*, Chen et al., 2011), and two are from a combination of multiple measurements ($8.4$ km s*-1*, $N=22$, Evans, 1967; $11.1$ km s*-1*, $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 $\chi^{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 emcee (Foreman-Mackey et al., 2013) to fit the seven orbital elements: period $P$, eccentricity $e$, argument of periapsis $\omega$, epoch of periapsis $\tau$, primary semi-amplitude $K_{1}$, mass ratio $q$, and the systemic velocity $\gamma$. 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\pm 0.001$ days) and eccentricity ($e=0.669\pm 0.002$) were well constrained, and are consistent with an independent analysis by Lagrange et al. (2019). We found a systemic velocity of $\gamma=12.18\pm 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.

3 Close stellar encounters with HD 106906

3.1 Identifying candidate perturbers

The revised systemic radial velocity was combined with astrometry from Gaia DR2 (Gaia Collaboration et al., 2018) to compute the position and kinematics of HD 106906 given in Table 1. 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 Wright & Mamajek (2018), augmented with members reported in Song et al. (2012). Astrometric measurements were obtained from Gaia DR2, while radial velocities came from a number of sources, either from various literature sources (Levato et al., 1996; Torres et al., 2006; Gontcharov, 2006; Kharchenko et al., 2007; Chen et al., 2011; Song et al., 2012; Dahm et al., 2012; Kordopatis et al., 2013; Desidera et al., 2015), 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_{\rm CA}$) were computed for each of the $10^{5}$ trials. This process was repeated for each star within the sample.

3.2 HIP 59716 and HIP 59721

\startlongtable

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 (catalog ) (F5V, Houk & Cowley, 1975) at $t_{\rm CA}=-3.49_{-1.76}^{+0.90}$ Myr, $D_{\rm CA}=0.65_{-0.40}^{+0.93}$ pc; and HIP 59721 (catalog ) (G9V, Torres et al., 2006) at $t_{\rm CA}=-2.18_{-1.04}^{+0.54}$ Myr, $D_{\rm 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_{\odot}$ 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 (catalog ) to the Sun 70  kya at $0.25^{+0.11}_{-0.07}$  pc (Mamajek et al., 2015), 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\arcsec$ 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 ($\Delta\pi=0.025\pm 0.053$ mas), and the proper motion difference ($0.61\pm 0.07$ mas yr*-1*) is less than the circular motion of a face-on 2,750 au orbit ($\sim 1.7$ mas yr*-1*). 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 ($\pi=7.56\pm 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*-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$\arcsec$ 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\arcsec$, had a parallax within 0.078 mas, and had total proper motion within 0.68 mas yr*-1*, implying a chance alignment probability of $\sim$10*-7*.

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^{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_{\rm CA}=-3.5_{-1.8}^{+0.9}$ Myr, $D_{\rm CA}=0.63_{-0.39}^{+0.92}$ pc for HIP 59716 and $t_{\rm CA}=-2.2_{-1.0}^{+0.6}$ Myr, $D_{\rm 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 $\sim$2 % of the simulations. The closest approach distances were not significantly different for the bound simulations, increasing slightly for HIP 59716 to $D_{\rm CA}=0.42_{-0.17}^{+0.30}$ pc, and decreasing slightly for HIP 59721 to $D_{\rm 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/\rho$ conversion factor peaking at unity for a uniform eccentricity distribution (Dupuy & Liu, 2011). Only in the most extreme cases where $a/\rho\sim 3$ and $e\sim 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.

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^{\prime}$ pointing toward the current location of the Sun and $Y^{\prime}$ and $Z^{\prime}$ pointing toward East and North from the perspective of the Sun. In this coordinate system, the plane of the disk is rotated $15^{\circ}$ about the $X^{\prime}$ axis and $85^{\circ}$ about the $Y^{\prime}$ axis. The second is rotated about the first such that the disk lies in one of the planes of the coordinate system ($Y_{\rm disk}Z_{\rm disk}$). A schematic diagram is shown in Figure 4.

The trajectories of both stars were almost coplanar with the current plane of the debris disk: $\Delta\theta=5\fdg 4\pm 1\fdg 7$ ($5\fdg 2_{-1.8}^{+1.5}$) for HIP 59716 and $\Delta\theta=4\fdg 2_{-1.1}^{+0.9}$ ($4\fdg 8\pm 1\fdg 0$) for HIP 59721 (Fig. 1, 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\pm 1.1$ km s*-1* and $5.3\pm 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.

4 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_{\rm CA}$ and $t_{\rm CA}$ using the linear approximation described previously. We assume that the radial velocities are perfectly measured. The correlations between $D_{\rm CA}$ and $t_{\rm CA}$ and the radial velocities of the three stars are shown in Figure 5. $D_{\rm CA}$ and $t_{\rm 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_{\rm CA}$ and $t_{\rm CA}$ follow lines of constant relative velocity. Increasing the relative velocity between HD 106906 and HIP 59716 causes a significant decrease in $D_{\rm CA}$ and an increase in $t_{\rm CA}$ (a more recent flyby). Interestingly, if HIP 59716’s radial velocity was similar to that of HIP 59721, the median value of $D_{\rm 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.

5 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 (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\sim 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. 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 ($<$104 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–21 Myr. Assuming a younger age of $\sim$10 Myr, the formation of the planet and migration into mean-motion resonance with the inner binary has to occur in $\sim$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 much more rapidly (e.g., Boss, 2011).

For the older age estimate of $\sim$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\times 10^{4}$ 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 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_{\odot}$ particle at $(0,0)$ pc, accounting for the combined mass of the HD 106906 binary, around which a 11 $M_{\rm Jup}$ particle was placed with an orbital eccentricity of $\log_{10}(1-e)$ drawn from $\mathcal{U}(-4,-1)$ (where $\mathcal{U}$ denotes the uniform distribution), a semi-major axis with a corresponding periapsis distance ($r_{\rm 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_{\odot}$, initialized at $y=-3$ pc, $\log_{10}(x)$ drawn from $\mathcal{U}(-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_{\rm peri}$ and 3,330 (18,615) led to a decrease in $r_{\rm peri}$, with 1,541 (1,525) resulting in the planet being ejected from the system. The correlation between the change in periapsis distance ($\Delta r_{\rm 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 $\Delta r_{\rm peri}$ and $D_{\rm CA}$, with closer encounters leading to a larger $\Delta r_{\rm peri}$ for a given initial eccentricity. The presence of a debris ring at $\sim$50 au around the HD 106906 binary suggests that the $r_{\rm peri}$ for the planet is now $>50$ au, otherwise the disk would be disrupted with each periastron passage for relative inclinations $\lesssim 10$ deg (Jílková & Zwart, 2015). Periapsis distances interior to the radius of the debris ring are possible for higher relative inclinations ($\gtrsim 40$ deg) without disrupting the disk. These $N$-body simulations suggest that $r_{\rm peri}$ for a highly eccentric planet can be plausibly raised from $\sim$1.5 to $\sim$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.

6 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 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*-1* significantly increases the probability of a closest approach distance within 0.1 pc (Fig. 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.