A Deep Polarimetric Study of the Asymmetrical Debris Disk HD 106906

We obtained deep polarimetric data of HD 106906 in the H band (peak λ = 1.647 μm) from GPIES. These data have a greater total integration time (2564.79 s compared to 709.94 s) and field rotation (203 compared to 71) than the polarimetric data presented in Kalas et al. (2015). We also obtained polarimetric J- (peak λ = 1.235 μm) and K1-band (peak λ = 2.048 μm) data through one of the Gemini Observatory's Large and Long Programs (PI: Christine Chen). The H-band data were taken on 2015 June 30 using the GPI's polarimetric (H-Pol) mode. The J- and K1-band data were later taken over a 4 day period from 2016 March 24 to 28, again in polarimetric (J-Pol and K1-Pol) mode. Due to a malfunction in the GPI's apodizer wheel at that time, all three observations used the coronagraph apodizing mask that was optimized for the H band ("APOD H G6205"). Further details about the polarimetric observations can be found in Table 1.

Table 1. Details of Observations with GPI for Each Band

Mode	Band	Date	texp (s)	N	ΔPA (deg)	MASS Seeing
Pol	J	2016 Mar 26	59.65	54	35.2	077
Pol	H	2015 Jun 30	59.65	43	20.3	027
Pol	K1	2016 Mar 28	88.74	40	36.5	133

Note. Here t_exp = the exposure time for each frame, N = the number of frames, and ΔPA = the total parallactic angle rotation.

Download table as:  ASCIITypeset image

Each of these sets of data were reduced through the GPI Data Reduction Pipeline (Perrin et al. 2014). In GPI's polarization mode, a pair of spots is created by each lenslet, representing both orthogonal polarization states Q and U. The first step in the reduction process is to create 3D cubes from the raw data, where the three dimensions contain an image for each orthogonal polarization state. The raw data are dark-subtracted, bad pixel–corrected, and cleaned of correlated detector noise, then combined into polarization data cubes. The images are then corrected for geometric distortion. The position of the star is found by measuring a set of fiducial satellite spots (Wang et al. 2014).

Next, the data cubes are combined to create the final composite data cube. In this final step, the instrumental polarization is measured from the apparent stellar polarization in each data cube by measuring the average fractional polarization inside the occulting mask, which is then subtracted from each cube (Millar-Blanchaer et al. 2015). Additionally, each cube is smoothed with a Gaussian kernel (FWHM of 1 pixel), and the satellite spot–to–star flux is measured in order to determine a photometric calibration factor. After this measurement, the 3D cubes are combined to create a single radial Stokes image, consisting of Stokes [I, Q_ϕ , U_ϕ , V]. Finally, using the stellar magnitude and photometric calibration factor for each filter+apodizer combination, the data were calibrated from ADU coadd⁻¹ to Jy arcsec⁻² (Hung et al. 2016).

The final images in each filter can be seen in Figure 1, where the left column represents the data in Q_ϕ , and the right column represents the data in U_ϕ . The white circles represent the size of the focal plane mask (FPM), with a radius of 009, 012, and 015, respectively.

Zoom In Zoom Out Reset image size

For this section, we calculated the peak surface brightness profile in each wavelength in order to quantify the surface brightness asymmetry. To do this, we find the peak brightness along the spine of the disk. After rotating the disk to be horizontal, we find the maximum ("peak") surface brightness along vertical slices. The peak surface brightness is then binned by 3 pixels along the horizontal axis and averaged. These results can be found in Figure 4, showing the peak surface brightness in mJy arcsec⁻² as a function of projected separation from the star in arcseconds. The error bars are calculated using the noise map in each filter.

Zoom In Zoom Out Reset image size

In all three bands, we observe the same general shape of the surface brightness profiles. Within ∼050, the surface brightness peaks closest to the star and decreases steeply toward larger stellar separations. Approximately between 050 and 070, the surface brightness briefly plateaus before again decreasing more gradually out to ∼11.

Looking at the relative peak surface brightness between each band, we confirm that there is a brightness asymmetry in all three wavelengths across the disk. Table 2 shows the percentage by which the SE extension is brighter than the NW extension, segmented by radial separation from the star. While the asymmetry differs slightly between each band, all three exhibit an SE extension that is roughly 10%–30% brighter than the NW extension. We find that the H band exhibits the greatest asymmetry with 30.9% ± 4.4% brighter flux in the SE extension between the inner working angle (IWA) and 040, about 10% greater than that measured in Kalas et al. (2015), and shows another peak of 27.0% ± 3.9% between 060 and 080. The asymmetries in the J and K1 bands are more modest, with a smaller asymmetry along the disk compared to the H band. While the J and K1 bands have larger uncertainties, because these three bands are very close in wavelength, we do not expect that the true surface brightness asymmetry varies too significantly between them. The asymmetry seen in the H band is more likely a better representation of the true brightness asymmetry observed at these near-IR wavelengths due to its much higher S/N. Regardless, this analysis shows that a modest surface brightness asymmetry exists along the disk in all three bands, and that the previously resolved asymmetry does not result from an artifact.

Table 2. Brightness Asymmetry between the NW and SE Extensions for Each Band as a Function of Stellar Separation

Band	IWA–040	040–060	060–080	080–10
J	26.9% ± 9.7%	9.7% ± 12.6%	22.3% ± 10.7%	18.3% ± 15.0%
H	30.9% ± 4.4%	18.9% ± 4.1%	27.0% ± 3.9%	21.0% ± 3.9%
K1	23.4% ± 11.2%	14.5% ± 11.2%	18.0% ± 12.3%	23.9% ± 17.2%

Note. Units are in a percentage of how much brighter the SE extension is compared to the NW extension.

Download table as:  ASCIITypeset image

Comparing the flux of the disk to the flux of the star in different bands tells us about the disk scattered-light color. Figure 5 shows the disk color as a function of the difference in the disk surface brightness minus the stellar magnitude between each band. We use the 2MASS J-, H-, and K-band magnitudes for HD 106906AB (6.946 ± 0.03, 6.759 ± 0.038, and 6.683 ± 0.026 mag, respectively). To measure the color, the disk flux in each band is integrated between 030 and 085 and over ∼015 vertically on both sides of the disk. The total flux is then converted to mag arcsec⁻² and compared between two bands, where the stellar magnitudes in each band are then subtracted (i.e., Δ(M_disk arcsec⁻² − M_star) for H − K1 = (H − K1)_disk − (H − K1)_star). Here a value above zero indicates a red disk color, a value below zero indicates a blue disk color, and a value of zero indicates a neutral disk color. We find that the disk predominantly exhibits a blue color, where calculating Δ(M_disk arcsec⁻² − M_star) between each band gives J − H = −0.38 ± 0.06 mag, J − K1 = −0.52 ± 0.07 mag, and H − K1 = −0.14 ± 0.05 mag. The difference between J − K1 and H − K1 also suggests that the disk color becomes increasingly more neutral toward longer wavelengths, consistent with disk color trends seen with other debris disks, such as HD 15115 and HR 4796A (Rodigas et al. 2012, 2015).

Zoom In Zoom Out Reset image size

The disk color gives some information about the dust grains in the disk, as color is associated with certain dust grain properties, such as the minimum dust grain size and porosity. For example, at a wavelength of 1.6 μm (∼H band) versus 2.2 μm (∼K band) and a porosity of zero, numerical simulations have shown that a strong blue disk color could be produced by a submicron minimum dust grain size, whereas a neutral disk color could be produced by a larger minimum dust grain size of several microns. A porosity greater than zero can also increase the estimated minimum dust grain size as well (Boccaletti et al. 2003; Kirchschlager & Wolf 2013). The fact that our H − K1 color measurement is only weakly blue and close to neutral suggests that the disk may have a larger minimum dust grain size. This being said, the minimum dust grain size is difficult to untangle and quantify from the disk color alone, as dust composition and the scattering phase function (SPF) also have an effect. What is more important is that we do not see a difference in color between the two extensions, implying that the dust grain properties do not significantly vary on either side of the disk.

To quantify the structure of the disk, we first measured the vertical height via the FWHM of the disk as a function of radial separation from the star. We chose to do this with the H band, as it has the highest S/N. We first rotate our image so that the disk's major axis is horizontal. Next, we measure the surface brightness along vertical cross sections of the disk at various radial separations from the star until the whole disk has been covered. We then fit a Gaussian to each vertical brightness profile and use this to extract the FWHM at each cross section. The results can be seen in Figure 6, where error bars are taken from fitting the Gaussian to the data.

Zoom In Zoom Out Reset image size

On both sides, we measure the FWHM out to r ≈ 11. We find the weighted mean FWHM (weights = 1/σ²) to be 016, represented by the gray dashed line in Figure 6. Comparing this value to the instrumental FWHM of the PSF (∼0054) confirms that the disk is fully resolved. This value of 016 is also ∼18% larger than the value measured in Kalas et al. (2015) of ∼013, although still consistent within the scatter of our FWHM data points. Subtracting the instrumental FWHM from our weighted mean in quadrature gives us an intrinsic FWHM of ∼015, or 15.6 au. While the FWHM fluctuates around the weighted mean, there appear to be no significant trends in the FWHM with stellocentric distance. This is consistent with a radially narrow ring, as projection effects would cause an increase in the vertical FWHM with distance if the disk were broad. A narrow ring also aligns with the conclusions made by Kalas et al. (2015) and Lagrange et al. (2016); however, we note that the vertical FWHM is only an indirect measure of the radial width, and we cannot make a firm conclusion based on it alone.

While the vertical height tells us about the disk structure, measuring the vertical offset, which traces the disk spine, may be even more revealing. The disk spine refers to where the peak surface brightness along the disk is located relative to the star in the vertical direction. To measure this, we refer back to our Gaussian fitting (Section 3.2), where the spine location is related to the mean value of the Gaussian. We can then compare the vertical offset between the spine location and the star as a function of radial separation. These results can be seen as the blue data points in Figure 7, where the error bars are taken from the Gaussian fitting. At first glance, the vertical offset looks to be asymmetric, where past 100 au (∼10), the vertical offset dips below −2 au in the SE extension, whereas this is not the case for the NW extension.

Zoom In Zoom Out Reset image size

To investigate this asymmetry further, we generate a simple inclined, circular, narrow ring model, following the same procedure used by Duchêne et al. (2020). Our free parameters include the ring radius (R_d ), inclination (i), PA measured from east of north, and offsets of the ring center from the star along the major and minor axes (δ_x and δ_y ). Here we are assuming that the star is located perfectly at the center of our data; however, it is important to note that some small systematic uncertainties may be present. A best-fitting model is then found by using a Markov Chain Monte Carlo (MCMC) fit via the Python package emcee (Foreman-Mackey et al. 2013); for more details about this method, see Section 4.2. Our best-fitting model can be seen in Figure 7, represented by the orange line. Through this modeling, we find ${R}_{d}={101.84}_{-1.61}^{+1.75}$ au, $i={85.57}_{-0\buildrel{\prime\prime}\over{.} 15}^{+0\buildrel{\prime\prime}\over{.} 14}$, and $\mathrm{PA}={104.31}_{-0\buildrel{\prime\prime}\over{.} 30}^{+0\buildrel{\prime\prime}\over{.} 16}$, where our inclination and PA are consistent with those derived in Kalas et al. (2015) and Lagrange et al. (2016). More interestingly, we find ${\delta }_{x}={16.48}_{-1.76}^{+2.04}$ and ${\delta }_{y}=-{2.99}_{-0.29}^{+0.28}$ au, which implies a significant eccentricity, placing the star closer to the SE extension. While possible systematic uncertainties may lower the significance of the small offset along the minor axis, this is not the case for the large offset along the major axis. Therefore, we can use δ_x to place a lower limit on an eccentricity of ≳0.16. Because modeling the spine suggests an inherently asymmetrical disk geometry, a circular ring model is no longer appropriate, where modeling with a proper ellipse will be needed to further probe the eccentricity and disk orientation.

To evaluate whether or not the x-offset is reasonable, we estimate the expected surface brightness asymmetry from a 16.48 au offset by comparing 1/r² (where r is the distance from the star to disk) and the SPF between the SE and NW at the same radial separations from the star. We first calculate r on both sides between 045 and 050, roughly halfway across the disk, where the SPF would not be affected by limb brightening near the ansae or the noise present close to the star. Using a disk radius of 101.84 au and a stellar offset of 16.48 au toward the SE extension, we obtain an SE/NW 1/r² ratio of ∼1.4. This means that based solely on the relative distances between the star and the disk, we would expect the SE extension to be about 1.4 times brighter than the NW extension at those separations, about 12% greater than what we measure with our H-band data (see Table 2). However, the difference between the SPF at those separations can also affect the expected brightness asymmetry. We therefore calculate the scattering angles between 045 and 050 and assume that the SPF in the HD 106906 disk conforms with the generic SPF observed in the solar system and other debris disks (Hughes et al. 2018) to estimate the SPF ratio between the SE and NW. Doing so, using our calculated scattering angles of ∼25°–27° for the NW and ∼29°–33° for the SE, we find that the NW extension has an SPF ∼1.25 times greater than the SE extension. Taking this into account, we expect that the greater SPF in the NW extension partially cancels out the 1/r² effect and therefore decreases the expected brightness asymmetry. While these are only rough estimates, as we do not know the exact SPF of this disk, they show that we would not expect a large brightness asymmetry despite the large offset derived from our spine fitting, and that the more moderate brightness asymmetry observed with our data is reasonable.

The 3D radiative-transfer code MCFOST is designed to model circumstellar disks around young objects and can create fully polarized, scattered-light, synthetic observations (Pinte et al. 2006). For our models, the "debris disk" option is chosen within MCFOST, where the volume profile is defined by the following equation (Augereau et al. 1999):

$$\begin{eqnarray}&&\rho (r,z)\ \propto \ \displaystyle \frac{\exp \left(-{\left(\tfrac{\left|z\right|}{h(r)}\right)}^{\gamma }\right)}{{\left[{\left(\tfrac{r}{{Rc}}\right)}^{-2{\alpha }_{\mathrm{in}}}+{\left(\tfrac{r}{{Rc}}\right)}^{-2{\alpha }_{\mathrm{out}}}\right]}^{1/2}},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&h(r)={H}_{0}\left(\displaystyle \frac{r}{{R}_{0}}\right),\end{eqnarray} \tag{ 2 }$$

In Equation (1), r is the radial distance from the star in the plane of the disk, and R_c is the critical radius inside and outside of which the surface density approaches power laws with indices α_in and α_out, respectively. In addition, z is the height above the disk midplane, and h(r) describes the scale height profile, with a vertical exponent of γ. Equation (2) defines h(r), where H₀ is the scale height at a reference radius, R₀.

To model the dust grain properties, Mie theory is used, where the dust grains are assumed to be spherical and compact. This idealized model, which does not capture the complexities expected in real dust grains, has implications that will be discussed later in this paper. However, for now, this assumption allows for easy and fast scattering computations. In addition, the dust grain size follows a power-law distribution of ${dN}(a)\propto {a}^{-q}{da}$, where a ranges from the minimum dust grain size ${a}_{\min }$ up to a size of 1 mm, and q is the power law for the dust grain size distribution. We also explore dust composition, looking exclusively at the volume fractions of astrosilicates, amorphous carbon, and water ice. The optical indices for each material are obtained from the literature (Si, Draine & Lee 1984; aC, Rouleau & Martin 1991; H₂O, Li & Greenberg 1998). The dust compositions (Si, aC, and H₂O) are set using effective medium theory, in that they make up one dust grain species (e.g., a single grain could consist of 50% Si, 25% aC, and 25% H₂O). First, the porosity is fixed for all dust grains that assume a single power-law size distribution. Then, the composition of the dust grains is allowed to vary as the relative fractions of each component are kept as free parameters.

The MCFOST code creates each model at a user-specified wavelength. For our purposes, we choose 1.647 μm (where the H filter throughput peaks). Here MCFOST produces a 4D fits image, which includes the final model in Stokes I, Q, U, and V intensities. The Q and U maps for each model are converted to Stokes Q_ϕ and convolved with an instrumental PSF. Both our H-band data and models are masked at radial separations greater than 12 and smaller than 012 to ensure that noisy pixels are left out.

To determine the best-fitting values for each parameter explored, the models must be compared to the data. To do this, we again utilize the Python package emcee (Foreman-Mackey et al. 2013). It uses MCMC, which allows us to find confidence intervals and best estimates for each parameter defining the density structure and dust grain properties of the disk.

In this process, we model the H-band data alone, as they have the highest S/N. To begin with, properties that are related to the density structure of the disk are first allowed to vary over a range of priors. The parameters that would most strongly affect the disk brightness are kept constant; this includes dust grain properties such as the minimum grain size and composition. This is done to limit parameter space and therefore speed up computation time. A list of density structural parameters, the initial value for each, and their prior range can be viewed in Table 3.

Table 3. The Density Structural Parameters and Dust Grain Properties Chosen to Be Constrained, Along with the Initial Value and Prior Limits Chosen for MCMC Analysis

Param.	Initial Val.	Limits
H0 [au]	5.0	[1, 10]
γ	2.0	[0.1, 3]
Rc [au]	100	[50, 200]
αin	1.0	[0.5, 15]
αout	−1.0	[−5, −0.5]
x* [au]	0.0	[−40, 40]
z* [au]	0.0	[−10, 10]
Mdust [M⊕]	0.033	[0.00033, 3.33]
Porosity	0.5	[0, 1]
VSi	0.5	[0, 1]
VaC	0.5	[0, 1]
${V}_{{{\rm{H}}}_{2}{\rm{O}}}$	0.5	[0, 1]
${a}_{\min }$ [μm]	1.5	[0.3, 4]
q	3.5	[2, 6]

Download table as:  ASCIITypeset image

Our free parameters for the density structure include H₀, γ, R_c, α_in, and α_out. We also include stellar offset parameters, x_* and z_*, where x_* is the offset along the disk major axis and z_* is the offset in the vertical direction relative to the disk. The parameter, x_*, is related to δ_x from our disk spine fitting in Section 3.3, where the main difference is that the star itself is being offset instead of the inclined ring, meaning that an offset toward the SE extension would be negative rather than positive. The parameter, z_*, differs from δ_y , as it represents an offset above or below the disk midplane, where δ_y is the projected offset along the minor axis. While our parameterized models do not account for geometrical asymmetries (the synthetic disk remains circular), so as to not make our models overly complex, a stellar offset can still be used to create a surface brightness asymmetry. Important structural parameters, such as inclination and PA, were not considered in this analysis, as they are already well defined from previous work, as well as our spine fitting, to be ∼85° and ∼104°, respectively (Kalas et al. 2015; Lagrange et al. 2016). Therefore, we fix these parameters in order to reduce the parameter space and computational time. Additionally, the inner radius is set to 50 au, based on observations from Kalas et al. (2015), and the outer radius to 200 au (roughly twice the outer radius observed in our data), where R₀ is fixed at 100 au.

After we model the density structure, we then hold these values constant at the median (50th) percentile values and instead vary the dust grain parameters. A list of the dust grain parameters, the initial value for each, and their prior range can be viewed in Table 3. Here M_dust is the total dust mass in Earth masses, porosity defines how fluffy or compact the dust grains are (where zero porosity corresponds to a perfectly compact dust grain), V_Si is the volume fraction of astronomical silicate, V_aC is the volume fraction of amorphous carbon, and ${V}_{{{\rm{H}}}_{2}{\rm{O}}}$ is the volume fraction of water ice. Finally, a_min is the minimum dust grain size, and q is the power law for the distribution of dust grain sizes.

For each run, 200 walkers are deployed to explore the parameter space, starting at the initial values as stated in Table 3. These walkers move around freely in order to find the maximum log-likelihood values for each parameter, as well as their posterior distributions. The likelihoods between our models and H-band data are calculated pixel by pixel using the χ² function. This occurs for up to several thousand iterations until after the parameters have clearly converged or the log likelihood for the walkers is no longer increasing.

For the MCMC analysis of the density structural parameters, emcee was allowed to run for a few thousand iterations, until well after convergence, to make sure that we had the best results. We then reran the MCMC analysis, where this time, the density structural parameters were kept fixed using the median percentile values (50%) obtained in the previous run. This analysis was allowed to run again until the log likelihood was no longer noticeably changing. The results for both the density structure and dust grain parameters can be found in Table 4, which lists the 16th, 50th, and 84th percentiles. The final model, created from the median percentile values from each parameter's confidence intervals, can be seen in Figure 8.

Zoom In Zoom Out Reset image size

Table 4. MCMC Results for Density Structural and Dust Grain Parameters

Param.	16%	50%	84%
H0 [au]	4.07	4.34	4.38
γ	0.708	0.715	0.723
Rc [au]	68.92	72.21	75.31
αin	0.79	1.03	1.32
αout	−2.34	−2.26	−2.16
x* [au]	−3.63	−3.58	−3.35
z* [au]	0.25	0.27	0.29
Mdust [M⊕]	0.012	0.02	0.11
Porosity	0.02	0.1	0.29
VSi	0.19	0.39	0.53
VaC	0.06	0.38	0.70
${V}_{{{\rm{H}}}_{2}{\rm{O}}}$	0.04	0.20	0.46
${a}_{\min }$ [μm]	0.91	3.17	3.72
q	2.99	3.19	3.30

Note. This includes the values for the 16th, 50th, and 84th percentiles.

Download table as:  ASCIITypeset image

From our MCMC modeling, we obtain an H₀ of ${4.34}_{-0.03}^{+0.04}$ au with a vertical profile exponent of $\gamma ={0.715}_{-0.01}^{+0.01}$. This would give the scale height to disk radius a ratio of h/r ≈ 4%, which is consistent with the h/r found for other debris disks such as β Pic, Fomalhaut, and AU Mic (Kalas et al. 2005; Krist et al. 2005; Millar-Blanchaer et al. 2015). However, the shallow γ appears to be in contradiction to the Gaussian profile used in Section 3.2, driven by the shape of our data, where a shallow γ would suggest a more Lorentzian profile. Given our limited S/N and model limitations, we note that it may be difficult to firmly establish the nature of the vertical profile. Next, we obtain a critical radius, R_c , of ${72.21}_{-3.30}^{+3.10}$ au, which lies roughly halfway between the set inner radius (50 au) and our disk radius of 101.84 au derived in Section 3.3. We also find surface density power laws of ${\alpha }_{\mathrm{in}}={1.03}_{-0.24}^{+0.29}$ and ${\alpha }_{\mathrm{out}}=-{2.26}_{-0.08}^{+0.11}$. These shallow power laws give the disk a dust density distribution ranging from a fixed inner radius of 50 au to a half-maximum at ∼111 au and an R_peak (radius of peak density) of ∼64 au, consistent with that derived in Lagrange et al. (2016). However, this is in contradiction to the assumption of a radially narrow ring based on our vertical FWHM profile, as these values give ΔR/R = 0.95, consistent with a radially broad disk. Although projection effects may still be present within our limited S/N data, another possibility is that the inner radius lies beyond 50 au. Given that the near-edge configuration of the disk makes it difficult to constrain R_in, we can only put an upper limit on ΔR/R of 0.95. Future modeling and analysis of polarized and spectral data may be able to help shed more light on the location of the inner radius.

For the stellar offset parameter, z_*, we find values that are very close to zero (0.25–0.29 au), meaning that the star lies within the midplane of the disk, as we would expect. Surprisingly, the MCMC favors a stellar offset of only ∼3.58 au toward the SE extension of the disk, which is much smaller than the stellar offset of ∼16.5 au predicted from our spine fitting. An offset this small is unable to match the surface brightness asymmetry observed in the H band. This can be seen in Figure 8, where there is a large residual remaining in the SE extension. One reason why the MCMC might prefer models with a smaller stellocentric offset could be the lack of a radial extent asymmetry in our data, as a larger offset would create a truncated SE extension. This would be an issue caused by the fact that, despite having a stellocentric offset, our models are geometrically symmetric circles.

While our confidence intervals for the density structure parameters were well constrained, constraining the dust grain parameters proved to be much more difficult. We obtain a total dust mass of ${0.02}_{-0.008}^{+0.09}$ M_⊕. This is on the same order as the masses derived from millimeter (∼0.05 M_⊕; Kral et al. 2020) and mid-IR observations (0.067 M_⊕; Chen et al. 2011); however, our confidence interval is large and spans two probability peaks at ∼0.016 and ∼0.11 M_⊕. We obtain a porosity of ${0.10}_{-0.08}^{+0.19}$, showing a significant preference for grains with a low porosity or compact grains. The composition is nearly unconstrained, with all three compositions (Si, aC, and H₂O) having very large confidence intervals and multiple peaks in their posterior distributions. The MCMC results slightly favor larger fractions of Si and aC, with smaller fractions of H₂O; however, we cannot make strong conclusions about the composition given their broad posteriors.

For the minimum dust grain size, we obtain ${a}_{\min }\,={3.17}_{-2.26}^{+0.55}$ μm, where we have a large probability peak at 0.91 μm and a second, smaller peak between 3 and 4 μm. While a minimum dust grain size of 0.91 μm is more consistent with the calculated blowout size for the system (∼0.85 μm), ¹⁸ a larger size may be more consistent with the disk color measured in Section 3.1. To test this, we compute Δ(M_disk arcsec⁻² − M_star) over the same integrated distances as in Section 3.1 for our median model at all three wavelengths. This is done for an a_min of both 0.91 and 3.17 μm. We find that an a_min of 3.17 μm is more consistent with the observed disk color (H − K1 = −0.15 mag), compared to our model with an a_min of 0.91 μm, which exhibits a much stronger blue color (H − K1 = −0.66 mag). Therefore, we find that the data favor a larger minimum dust grain size, although a_min could be further informed in future work using the degree of linear polarization.

Finally, we obtain a power law for the dust grain size distribution of $q={3.19}_{-.20}^{+0.11}$. This is slightly smaller than the predicted size distribution slope of q ≈ 3.5 for a collisional cascade (Dohnanyi 1969; Marshall et al. 2017); however, other debris disks have been found to deviate from this value. Additionally, our obtained value for q is comparable to the size distribution slope found for the millimeter observations of 3.36 ± 0.13 (Kral et al. 2020) within 1σ.

Comparing our H-band data with the best-fit model, we see that while the NW extension seems to be modeled well with no lingering residuals, as stated before, the SE extension features a significant residual along the disk. We find this to be more evidence for a geometrically asymmetric and likely eccentric disk, as our symmetric/circular models are unable to fully model both the disk structure and surface brightness asymmetry simultaneously. This may have affected our final modeling results, since we would need a model that has an asymmetric dust density profile to more accurately model our asymmetric disk. However, we leave this more complex modeling for future work.

There are several limitations with our MCFOST models that could impact our results. The first, as discussed above in Section 4.3, is that our models are not inherently eccentric. While we are able to offset the star, we assume the disk to be circular with a symmetric dust density profile, which would not be the case for a truly eccentric disk. This assumption has likely caused our final model to underestimate the surface brightness asymmetry, as seen in the residual in Figure 8. Following the same steps as in Section 3.3, we compute 1/r² and the SPF between 045 and 050, this time using a stellar offset of 3.58 au. Doing so, we find that the effects from 1/r² and the SPF at these radial separations almost completely cancel each other out, meaning that no brightness asymmetry would be seen with an offset this small. While a larger offset is needed to reach the observed brightness asymmetry, increasing the stellar offset in our symmetric models creates a more ill-fitting model, as this creates a radial extent asymmetry not observed in our data. This leads us to believe that in order to better represent the disk, we need a model that is eccentric with an asymmetric dust density profile.

Another limitation is that our models assume a smooth dust mass distribution defined with only one power law, which may not necessarily be realistic, as certain sizes of dust grains may be missing or the distribution may more closely resemble a function composed of more than one power law (Strubbe & Chiang 2006). High-resolution Atacama Large Millimeter/submillimeter Array data would be useful in this case to constrain the dust mass distribution with millimeter observations, which can then be used within MCFOST to create a more realistic model.

The final limitation comes from the use of Mie theory in our models. As mentioned before, Mie theory uses compact spherical grains, which makes computations much faster and easier but is unrealistic. In reality, dust grains are not perfectly spherical but rather irregularly shaped aggregates. This may strongly affect the accuracy of our results for the dust grain properties, as they rely heavily on this assumption and how light is scattered. For example, Min et al. (2015) found that the SPF differed between Mie dust grains and realistic aggregates for scattering angles greater than 90° (back side of the disk). Where for Mie grains, the SPF generally decreased, they found that for aggregates, the SPF either increased or stayed flat. Another discrepancy due to Mie theory is the effect this has on the blowout size. Arnold et al. (2019) calculated a_blow for aggregates with a porosity of 76.4% and found that, in general, it was two to three times higher than for Mie compact grains. In addition, multiple studies of other debris disks have also found that Mie theory is unable to match the observations (Milli et al. 2017; Arriaga et al. 2020; Duchêne et al. 2020). While we leave more complex modeling with realistic aggregates for future work, we believe our results are still useful in comparing to previous work done on HD 106906 and other debris disks.

In this section, we compare our results to previous simulation work, including both works focused specifically on HD 106906 and more general disk–planet interaction studies.

The empirical analysis of our data supports a significant eccentricity; this could be due to interaction with HD 106906b, either as it was ejected from the disk or due to secular orbital effects. There are two previous studies that consider these scenarios using particle simulations (Nesvold et al. 2017; Rodet et al. 2017). While Rodet et al. (2017) examined the morphological outcomes of the disk from both the planet being ejected and secular effects, Nesvold et al. (2017) solely examined the secular effects.

Rodet et al. (2017) found that the ejection of the planet from the disk was able to reproduce the asymmetries on larger scales, as seen with HST, where particles are blown out to high eccentricities on the right side of the disk. However, the ejection was not able to reproduce the SE/NW asymmetry on smaller scales, as seen with GPI in this work and Kalas et al. (2015). On the other hand, the effect of the planet on an eccentric, inclined orbit was able to reproduce the smaller-scale asymmetries with a larger density in the SE than the NW. This also gives rise to a needle-like halo at larger distances, but not as strongly as the ejection scenario, suggesting that the two effects together could explain the observations on large and small scales.

In contrast, Nesvold et al. (2017) conducted an in-depth study of the effect on the disk from the planet, initially located at 700 au (in situ) with an adopted eccentricity of 0.7, as it orbits at different inclinations. Through this alone, they were able to recreate both the small- and large-scale asymmetries, where a needle-like structure was seen up to 500 au, while closer in, the SE extension is brighter than the NW extension. The simulation also predicts that the planet should induce a large eccentricity on the disk, up to ∼0.3, which agrees well with the large eccentricity (e ≳ 0.16) we derive from our spine fitting. One caveat is that this model predicts a greater surface brightness asymmetry than we observe in our GPI data, particularly past 50 au (∼50% brighter SE extension compared to ∼18%–27% derived in Section 3.1); however, this is taking into account the total scattered light compared to polarized light.

In a more general case, Lee & Chiang (2016) looked at a number of different debris disk morphologies induced by a single eccentric planet companion. One of these morphologies is the "Needle," which is described as an eccentric and vertically thin disk viewed nearly edge-on. In their analysis, they found that a single planet with an eccentricity of 0.7 creates a radial length asymmetry, where one side extends much farther than the other (the needle), as well as a brightness asymmetry. In this case, the shorter side of the disk is brighter closest to the star, but eventually, at larger distances, the longer side becomes brighter. This is consistent with observations of the HD 106906 debris disk, as the "Needle" is observed with HST at large radial separations in the NW extension, where a brightness asymmetry is observed with GPI at smaller radial separations. While in these simulations, the planet has a mass of 10 M_⊕ and is orbiting interior to the disk, the results produced are similar to the simulations done with a larger planet on an eccentric orbit outside of the disk for HD 106906. While we cannot fully discriminate against a planet perturber interior to the disk, HD 106906b is still a likely cause for the observed disk morphology.

These simulations show that a single eccentric planet on both an internal and an external orbit can cause a debris disk to become eccentric, as well as produce a brightness asymmetry consistent with what we observe with our GPI data. Additionally, Lee & Chiang (2016) and Nesvold et al. (2017) showed that this is sufficient to create the needle-like halo that is observed with HST. The recent constraints of HD 106906b's orbit also indicate that it most likely has a misaligned and eccentric orbit, with a periastron distance that comes reasonably close to the star (Nguyen et al. 2021). This is similar to the conditions used for the planet's orbit in Nesvold et al. (2017), which are able to reproduce both the small- and large-scale asymmetries, providing further evidence for the planet being the disk's perturber. Based on simulation work and our observations, dynamical interaction with HD 106906b remains a strong candidate for the source of the disk's asymmetries.

Another possibility for the observed surface brightness asymmetry is an interaction with the ISM (Debes et al. 2009). If a debris disk passes through a dense cloud of interstellar gas, small grains can be stripped from the disk into the halo due to ram pressure. This has been used as one attempt to explain the morphological features in other debris disks, such as the "wing" feature for HD 32297 and HD 61005 (Schneider et al. 2014), as well as the "needle" feature for HD 15115 (Rodigas et al. 2012).

An ISM interaction comes up as a scenario for HD 106906 because of the disk's asymmetrical needle-like halo, where passing through a dense ISM cloud would cause small dust grains to be blown into the NW extension of the halo, as observed in the optical with HST (Kalas et al. 2015). For this to occur, we would expect the proper motion of the binary to be in the opposite direction of the needle (SE direction). Here HD 106906AB has a proper motion of −39.014 mas yr⁻¹ in R.A. and −12.872 mas yr⁻¹ in decl. (Gaia Collaboration et al. 2018), meaning that the proper motion is pointing in the SW direction. This is inconsistent with what we would expect for an ISM interaction. In addition, we do not see a difference in the color along the disk, where we might expect an ISM interaction to cause the NW extension of the disk to be more blue due to the direction of small grains blown and distributed toward the NW.

Taking these inconsistencies into consideration, we determine that an ISM interaction is unlikely to be the source of the disk's asymmetry. Based on our findings of an asymmetric structure and consistent disk color for both extensions, the source of asymmetry should be an event that changes the dust density distribution without greatly changing the dust grain size distribution. In this case, we find that dynamical interaction with the planet companion is more consistent with observations.