HD 106906: A Case Study for External Perturbations of a Debris Disk

Circumstellar debris disks are produced by the rocky and icy material leftover from the formation of the star and any planets in the system. To date, over 1700 debris disks have been detected via the excess infrared emission in their star's spectral energy distribution (Cotten & Song 2016), and over 40 have been resolved with optical or infrared imaging (http://circumstellardisks.org). The architecture of the underlying planetary system can leave a distinct imprint on the morphology of a debris disk (Mouillet et al. 1997; Wyatt et al. 1999; Matthews et al. 2014; Nesvold & Kuchner 2015b; Lee & Chiang 2016; Nesvold et al. 2016).

Modeling debris disk morphology is often focused on the effects of a planetary-mass perturber orbiting interior to or within the disk (Mouillet et al. 1997; Chiang et al. 2009; Nesvold & Kuchner 2015a; Pearce & Wyatt 2015). Nonetheless, exterior companions have been detected and inferred for several systems (Rodriguez & Zuckerman 2012; Bailey et al. 2014; Mawet et al. 2015), and dynamical modeling suggests that an exterior perturber can also excite eccentricities of the particles in a debris disk, via Laplace–Lagrange secular perturbations in the near-coplanar case (Thébault et al. 2012) or Kozai–Lidov perturbations for mutually inclined companions and disks (Nesvold et al. 2016), inducing asymmetries in the disk and triggering a collisional cascade. Collisions between the disk particles produce smaller dust grains whose thermal emission or scattered light can then be spatially resolved with infrared or optical imaging (Wyatt 2008).

HD 106906 is an ideal example of a system with an exterior perturber, as it harbors a confirmed exterior companion with a model-atmosphere-derived mass of $11\,{M}_{\mathrm{Jup}}$ at a projected separation of 650 au outside a resolved disk (Bailey et al. 2014). Scattered-light imaging of the disk with the Gemini Planet Imager (GPI), the Hubble Space Telescope's Advanced Camera for Surveys (HST/ACS), and SPHERE has revealed that the disk is a ring viewed nearly edge-on (inclination $\sim 85^\circ$), with a inner region cleared of small dust grains (Kalas et al. 2015; Lagrange et al. 2016). These observations noted four major features of the disk morphology and system geometry:

• 1.  
  The position angle (PA) of the disk is oriented $\sim 21^\circ$ counterclockwise from the PA of the companion, which constrains the orbit of the companion relative to the disk.
• 2.  
  The inner disk has little to no vertical extension. While Kalas et al. (2015) tentatively suggested the presence of a "warp" in the disk's vertical structure on the west side of the disk, this warp was not confirmed by Lagrange et al. (2016). This lack of vertical extension indicates that the inclinations of the disk particles have not been excited.
• 3.  
  The east side of the disk is brighter than the west side in GPI and SPHERE near-infrared images.
• 4.  
  Kalas et al. (2015) observed a faint extension on the west side of the disk out to nearly 500 au, but only diffuse emission on the east side.

These latter two features indicate that the disk may be an eccentric ring, which will exhibit a brightness asymmetry toward the pericenter side ("pericenter glow"; Wyatt et al. 1999) and a faint, extended tail toward the apocenter side (Lee & Chiang 2016).

We modeled the HD 106906 system to demonstrate that the observed exterior companion can shape the disk into a flat, eccentric ring, and that all four of these morphological features can be reproduced without invoking the presence of a second companion. We also used the observed features and asymmetries of the HD 106906 disk to place constraints on the orbit of the observed companion. In Section 2, we describe the collisional and dynamical simulations we performed on the parent bodies and dust grains in the HD 106906 disk. In Section 3, we present the simulated brightness images produced by our simulations for comparison with observations. In 4, we discuss the implication of these results and show how they can be used to constrain the orbit of HD 106906b. In Section 5, we summarize our conclusions and suggest opportunities for future work.

Given that collisions between particles in a disk with sufficiently high optical depth (${L}_{\mathrm{IR}}/{L}_{* }\approx 1.4\times {10}^{-3}$ for HD 106906 Chen et al. 2011) will both produce the small grains seen in observations and may affect the dynamics of the disk, we simulated the HD 106906 system using the Superparticle-Method Algorithm for Collisions in Kuiper belts and debris disks (SMACK; Nesvold et al. 2013). We then recorded the dust-producing encounters between parent bodies, simulated the orbits of the generated dust grains under the influence of radiative forces following the method of Lee & Chiang (2016), then simulated the surface brightness of the dust using a Henyey–Greenstein scattering phase function (Henyey & Greenstein 1941).

2.1. SMACK Model

SMACK is based on the N-body integrator REBOUND (Rein & Liu 2012), but approximates each particle in the integrator as a collection of bodies with a range of sizes between 1 mm and 10 cm in diameter, traveling on the same orbit. This group of bodies is called a "superparticle" and is characterized by a size distribution, position, and velocity. The superparticles act as test particles in the integration, and orbit the star under the influence of perturbations by any planets in the simulation. Each superparticle is approximated as a sphere with some finite radius. When REBOUND detects that two superparticles are overlapping in space, SMACK statistically calculates the number of bodies within each superparticle that will collide and fragment, removes these bodies from their size distributions, and redistributes the fragments. SMACK also corrects the trajectories of the parent superparticles to conserve angular momentum and energy, compensating for the kinetic energy lost to fragmentation.

The parameters for the SMACK simulation of the HD 106906 system described in this work are listed in Table 1. The masses of the star and companion were 2.5 ${M}_{\odot }$ and 11 ${M}_{\mathrm{Jup}}$, respectively. The initial semimajor axis range of 65–85 au for the 10,000 superparticles in the simulated disk was chosen in anticipation that the disk would spread during the 15 Myr course of the simulation. The radial extent of the HD 106906 ring as observed in scattered-light imaging is $\sim 50\mbox{--}100\,\mathrm{au}$ (Kalas et al. 2015; Lagrange et al. 2016). The orbital parameters of the companion were chosen such that the gravitational perturbations from the companion would excite the eccentricities of the disk particles without causing a vertical extension of the disk on the timescale of the system's age, and such that the simulated companion's orbit could reproduce the position of the observed companion.

Table 1.  Initial Conditions of the Disk and Companion for the Simulation

Parameter	Initial Disk Values	HD 106906b
Semimajor Axis (au)	65–85	700
Eccentricity	0–0.01	0.7
Inclination (°)	0–0.29	8.5
Longitude of Nodes (°)	0–360	90
Argument of Pericenter (°)	0–360	−90
Optical Depth	$1.4\times {10}^{-3}$	⋯
Density (g cm−3)	1.0	⋯

Download table as:  ASCIITypeset image

2.2. SMACK versus Collisionless N-body

To measure the effects of collisions on the dynamics of the disk particles, we also performed a collisionless N-body simulation of the disk using the Wisdom–Holman integrator of REBOUND with collision detection and resolution turned off. The collisionless N-body simulation used the same companion and disk parameters as the SMACK simulation (Table 1), with 10,000 test particles to represent the disk. Figure 1 shows the time evolution of the simulated disk's average eccentricity, inclination, longitude of nodes, and argument of pericenter. Although there are small variations, most notably in the average eccentricity, between the SMACK simulation and the collisionless N-body simulation, the maximum difference for each parameter is $\leqslant 10 \%$, indicating that fragmenting collisions have a minimal effect on the dynamics of this system.

Zoom In Zoom Out Reset image size

2.3. Dust Model

To generate the simulated images of the dust grains in the HD 106906 system, we adapted the method of Lee & Chiang (2016), which extended the dust orbit calculations of Wyatt et al. (1999) to include estimates of the surface brightness. Our SMACK simulation output the locations of dust-producing collisions during the 15 Myr simulation, as well as the orbits of the parent bodies producing the dust. We selected the first 10⁴ dust production events occurring after time t = 5 Myr. For each dust production event, we generated 10 dust orbits, each with a β value randomly chosen from a power-law distribution with an index 3/2 (where $\beta \approx {F}_{\mathrm{rad}}/{F}_{\mathrm{grav}}$ represents the ratio of the radiative and gravitational forces acting on a dust grain). The maximum possible value for the β value was set by the parent body's orbit:

$$\begin{eqnarray}&&{\beta }_{\max }=\displaystyle \frac{1-{e}_{{\rm{p}}}}{2(1+{e}_{{\rm{p}}}\cos {f}_{{\rm{p}}})},\end{eqnarray} \tag{ 1 }$$

where ${e}_{{\rm{p}}}$ and ${f}_{{\rm{p}}}$ are the parent body's eccentricity and true anomaly, respectively. The semimajor axis a, eccentricity e, and argument of pericenter ω of each dust orbit are given by the parent body's orbit (specifically ${a}_{{\rm{p}}}$, ${e}_{{\rm{p}}}$, and ${f}_{{\rm{p}}}$) and the β value assigned to the orbit:

$$\begin{eqnarray}&&a=\displaystyle \frac{{a}_{{\rm{p}}}(1-{e}_{{\rm{p}}}^{2})(1-\beta )}{1-{e}_{{\rm{p}}}^{2}-2\beta (1-{e}_{{\rm{p}}}\cos {f}_{{\rm{p}}})},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&e=\displaystyle \frac{\sqrt{{e}_{{\rm{p}}}^{2}+2\beta {e}_{{\rm{p}}}\cos {f}_{{\rm{p}}}+{\beta }^{2}}}{1-\beta },\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\omega ={\omega }_{{\rm{p}}}+\arctan \left(\displaystyle \frac{\beta \sin {f}_{{\rm{p}}}}{{e}_{{\rm{p}}}+\beta \cos {f}_{{\rm{p}}}}\right).\end{eqnarray} \tag{ 4 }$$

The inclination i and longitude of nodes Ω of the dust orbit were set to be equal to the corresponding values for the parent body, ${i}_{{\rm{p}}}$ and ${{\rm{\Omega }}}_{{\rm{p}}}$, respectively. For each dust orbit, we generated 10 dust grains, and assigned each a mean anomaly selected randomly from a uniform distribution between 0° and 360°. Thus, each dust production event produced 100 final dust grain locations.

After constructing the dust population from the SMACK results, we simulated the surface brightness of the dust using $\phi (g,\theta )/{\beta }^{2}{r}^{2}$, where $\phi (g,\theta )$ is the Henyey–Greenstein scattering phase function with asymmetry parameter g, θ is the angle between the dust grain and the observer's line-of-sight (with the vertex at the star), and r is the distance between the dust grain and the star. Following Lee & Chiang (2016), we used g = 0.5.

Figure 2 shows the simulated brightness of the dust produced by SMACK after 5 Myr, scaled for comparison with Figure 1 of Lagrange et al. (2016) and Figure 3 of Kalas et al. (2015). The SMACK-produced dust population exhibits a brightness asymmetry at pericenter (Figure 2(a)) as well as a faint extension on the apocenter side (Figure 2(b)), although it does not reproduce the diffuse emission on the eastern side of the disk suggested by Kalas et al. (2015).

Zoom In Zoom Out Reset image size

The scattered-light brightness enhancement at pericenter is a signature of an eccentric ring (Wyatt et al. 1999; Pan et al. 2016), indicating that the orbits of the disk particles in the simulation, initially assigned eccentricities $\lt 0.01$ and random longitudes of node and arguments of pericenter, have become more eccentric and apsidally aligned due to their secular resonance with the distant companion (Li et al. 2014; see Section 4). The extended "tail" seen on the apocenter side of the disk is also a signature of an eccentric ring, in which high-eccentricity dust grains are produced near the ring's pericenter on apsidally aligned orbits, and then observed as they travel to and from their distant apocenters (Lee & Chiang 2016).

The relatively low inclination of the companion relative to the disk results in a flat, narrow ring, with no appreciable vertical extension after 5 Myr, and the observed relative inclination between the companion and the disk ($\sim 21^\circ$) is reproduced by our choice of 5° relative inclination and the 85° viewing inclination. The observed exterior companion is therefore able to reproduce four of the observed morphological features of the disk with no requirement for a second companion.

The physical mechanisms described above allow us to place constraints on the orbit of the companion. The relative PA of the companion and the disk's line of nodes on the sky is related to the inclination of the companion's orbit relative to the plane of the disk, but this relationship is complex and also depends on the longitude of nodes and argument of pericenter of the companion relative to the disk. In addition, these angles change with time, as the gravitational perturbations from the companion cause the orbits of the disk particles to precess coherently together (see Section 4). Instead, we can use the disk's morphology to place constraints on the companion's orbit. For example, the pericenter of the companion cannot be too close to the outer edge of the disk or the companion's chaotic zone will truncate the disk. As a rough estimate, we can calculate the relationship between the radius of the inner edge of the companion's chaotic zone, ${r}_{{\rm{z}}}$, and the companion's pericenter distance, ${r}_{\mathrm{pl}}$, using the analytically derived classical chaotic zone relationship for circular orbits (Wisdom 1980), $({r}_{\mathrm{pl}}-{r}_{{\rm{z}}})/{r}_{\mathrm{pl}}=1.3\,{\mu }^{2/7}$, where μ is the companion-to-star mass ratio. If we set the inner edge of the companion's chaotic zone to be the outer edge of the disk, ${r}_{{\rm{z}}}=100\,\mathrm{au}$, the minimum pericenter location for the companion is ${r}_{\mathrm{pl}}\approx 138\,\mathrm{au}$. This constraint contains a degeneracy between the companion's eccentricity and semimajor axis, ${a}_{\mathrm{pl}}(1-{e}_{\mathrm{pl}})\gtrsim 138\,\mathrm{au}$.

There also exists an upper limit on the companion's semimajor axis, as the secular timescale must be less than the age of the system for the companion's secular perturbations to produce the observed asymmetries in the disk. If we constrain the secular timescale to be at least 10 Myr, the companion's semimajor axis and eccentricity are constrained by ${a}_{\mathrm{pl}}\sqrt{1-{e}_{\mathrm{pl}}^{2}}\lesssim 661\,\mathrm{au}$ (see Section 4).

We can use the disk's vertical extent to place an upper limit on the mutual inclination between the plane of the disk and the companion's orbit. We ran three collisionless N-body simulations of the disk perturbed by an 11 ${M}_{\mathrm{Jup}}$ companion with the same orbital parameters as in the SMACK simulation, but varying the mutual inclination between the companion and disk to 85, 20°, or 30°. The N-body particles represented the parent bodies in the disk. We simulated the production of dust grains by generating and recording one dust orbit matching the location and velocity of each parent body every 10 years during each simulation. We then generated 100 dust grains from each dust orbit and produced simulated brightness images using the procedure described in Section 2.3. Figure 3 shows the simulated brightness of each disk, viewed 5° from edge-on, at 10 Myr, as well as the location of the parent bodies in each disk. Perturbations from a higher companion inclination produce a larger vertical extent in the disk. Resolved images of the system show a flat disk, indicating that the companion's inclination relative to the disk must be $i\lesssim 20^\circ$. Our SMACK simulation demonstrates that the observed exterior companion can excite the necessary eccentricities in the ring within the age of the system without creating a significant vertical extension if the companion has a moderately large eccentricity ($e\approx 0.7$) but a small inclination ($i\approx$ 85). Given these orbital parameters, the companion's semimajor axis would need to be ∼700 au to match its observed projected position.

Zoom In Zoom Out Reset image size

At time t = 0 year in our SMACK simulation of the HD 106906 system, the disk particles have eccentricity $\leqslant 0.01$ and longitudes of node and arguments of pericenter distributed randomly between 0° and 360°, forming a circular belt. By time t = 5 Myr, gravitational perturbations from the companion have increased the average eccentricity of the particles to ∼0.18 (Figure 4). Increasing the average eccentricity of the particles alone would only produce a broader circular disk, but the secular resonance formed with the companion also causes the particles' longitudes of node and arguments of pericenter to converge. In other words, the orbits of the disk particles begin to apsidally align, producing a coherent eccentric ring. Figure 4 illustrates this with plots of the time evolution over 50 Myr of the longitudes of node, arguments of pericenter, eccentricities, and inclinations of 10 randomly chosen disk particles in a collisionless N-body simulation with the same system parameters as the SMACK simulation described in Table 1. We used 10,000 particles to represent the disk, simulated with the Wisdom–Holman integrator in REBOUND with collisions turned off. The orbits of these 10 randomly chosen particles become roughly apsidally aligned within ∼4 Myr. This behavior is consistent with the hierarchical nearly coplanar secular evolution, investigated in Li et al. (2014), which showed that the resonance angle for low-inclination companions is the sum of the longitude of nodes (Ω) and the argument of pericenter (ω). The test particles, although they were initially assigned random values of Ω and ω, are captured into resonance with the companion, which both pumps up their eccentricity and aligns their orbits, forming an eccentric disk.

Zoom In Zoom Out Reset image size

The secular precession timescale of particles in the disk, due to an external planetary-mass companion with mass ${M}_{\mathrm{pl}}$ and eccentricity ${e}_{\mathrm{pl}}$, is defined as (Naoz 2016)

$$\begin{eqnarray}&&{t}_{\sec }\sim \displaystyle \frac{16}{30\pi }\displaystyle \frac{{M}_{* }+{M}_{\mathrm{pl}}}{{M}_{\mathrm{pl}}}\displaystyle \frac{{P}_{\mathrm{pl}}^{2}}{{P}_{{\rm{d}}}}{(1-{e}_{\mathrm{pl}}^{2})}^{3/2},\end{eqnarray} \tag{ 5 }$$

where ${P}_{\mathrm{pl}}$ is the companion's period around the star (mass M_*) and ${P}_{{\rm{d}}}$ is the period of particles in the disk. The eccentric companion induces gravitational perturbations that result in the disk particles orbiting in a resonance, where the disk particles' longitude of pericenter $\varpi =\omega +{\rm{\Omega }}$ is the resonant angle (Li et al. 2014).

The precession timescale can be used to place an upper limit on the companion's semimajor axis. In order to perturb the entire disk (down to its inner edge at ∼50 au) within the age of the system ${t}_{\mathrm{age}}\approx 10\,\mathrm{Myr}$, the inner edge of the disk must have experienced at least one half-cycle of secular perturbation, so the secular timescale at 50 au must be $\tfrac{1}{2}{t}_{\sec }\lesssim {t}_{\mathrm{age}}$. Using a stellar mass of 2.5 ${M}_{\odot }$ and a companion mass of 11 ${M}_{\mathrm{Jup}}$, and Equation (5), this yields

$$\begin{eqnarray}&&{a}_{\mathrm{pl}}\sqrt{1-{e}_{\mathrm{pl}}^{2}}\lesssim 661\,\mathrm{au}.\end{eqnarray} \tag{ 6 }$$

We have shown that the observed exterior companion in the HD 106906 system can shape the system's debris disk into a flat, eccentric, dust-producing ring and reproduce its observed morphological features and asymmetries. Our SMACK simulations also allow us to place constraints on the orbit of the companion using the morphology of the disk.

While we have demonstrated that the perturbations from the observed, exterior companion can excite eccentricities in the HD 106906 ring, alternative mechanisms for eccentricity excitation also exist. For example, a second companion on an eccentric orbit interior to the debris ring could force an eccentricity on the ring. Future simulations investigating the plausibility of this scenario may be able to constrain the orbit of the outer companion based on stability requirements.

Constraining the orbit of HD 106906b could have implications for its formation scenario. Prior to the publication of resolved images of the disk, it was suggested (using N-body simulations) that the companion formed interior to the disk and was scattered onto a highly eccentric orbit (Jílková & Portegies Zwart 2015). This study concluded that the disk can survive perturbations by a companion with an apocenter distance of 650 au and a pericenter distance interior to the disk if the companion's inclination is $\gtrsim 10^\circ$. However, this resulted in a significantly vertically perturbed disk by 10 Myr, regardless of the companion's inclination. Our simulations indicate that a companion with an orbit completely exterior to the disk can reproduce the observed asymmetries without vertically extending the disk, supporting the scenario in which the companion formed in situ.

A more thorough exploration of the parameter space may be able to place further constraints on the companion's orbit using the observed geometry of the system. The methodology we presented in this work can also be generalized to other debris disk observations to explore whether their observed asymmetries could be explained by the presence of an undetected distant exterior companion, and to investigate how the morphologies of these disks could constrain the orbits of their exterior perturbers. Other debris disk systems with exterior massive perturbers are likely not uncommon; surveys indicate that $\sim 25 \%$ of debris disks exist in binary or triple star systems (Rodriguez & Zuckerman 2012). Planetary-mass exterior companions like HD 106906b may be responsible for the asymmetries in observed debris disks such as HD 61005 (Hines et al. 2007; Esposito et al. 2016) and HD 15115 (Rodigas et al. 2012; Schneider et al. 2014; MacGregor et al. 2015), for example.

Numerical simulations were performed on the Memex High Performance Computing Cluster at the Carnegie Institution for Science. E.N. was supported by the Carnegie DTM Postdoctoral Fellowship. S.N. acknowledges partial support from a Sloan Foundation Fellowship. The authors wish to thank the anonymous referee for a prompt and helpful review.