An Exo–Kuiper Belt with an Extended Halo around HD 191089 in Scattered Light

2.1.1. HST/NICMOS (2006)

On 2006 May 27, HD 191089 was observed using HST/NICMOS (Proposal ID: 10527; PI: D. Hines) with the NIC2-CORON aperture and F110W filter (λ_c = 1.12 μm, inner working angle (IWA): 03, pixel scale: 75.65 mas pixel⁻¹; Viana et al. 2009) with two telescope orientations each observing the target for 2303.67 s (see Table 2 for the observation log). These HST/NICMOS data, totaling 16 frames, were previously presented in Soummer et al. (2014). We perform an updated reduction including application of two different point-spread function (PSF) subtraction algorithms to remove the starlight and speckle noise and reveal the debris disk around HD 191089.

We obtain the scattered-light image of the system using the multireference differential imaging (MRDI) technique. Specifically, we retrieve 849 F110W exposures of 70 diskless reference stars in the ALICE archive of the NICMOS observations (PI: R. Soummer; Choquet et al. 2014; Hagan et al. 2018).⁴⁴ For each target exposure, we first decompose its 10% closest ALICE images in correlation (i.e., L²-norm sense and a total of 85 images) using both the Karhunen–Loève image projection (KLIP; Soummer et al. 2012) and nonnegative matrix factorization (NMF;⁴⁵ Ren et al. 2018b) data reduction methods, then model the target with these components. The disk then resides in the residual image when the empirical PSF model is subtracted from the target exposure. The final NICMOS disk image is then the elementwise mean of the derotated individual residual exposures. We list the reduction parameters in Table 3. ⁴⁶

To calibrate the NICMOS disk image and obtain a measurement of its surface brightness, we multiply the reduced data by the calibrated F110W PHOTFNU parameter⁴⁷ F_ν = 1.21 × 10⁻⁶ Jy s count⁻¹, then divide the data by the NICMOS pixel area on sky to obtain the surface brightness data in units of Jy arcsec⁻². See Figure 1 for the results.

Zoom In Zoom Out Reset image size

2.1.2. HST/STIS (2014)

On 2014 July 19 and August 13 (two visits each), HD 191089 was observed using HST/STIS (Program ID: 13381; PI: M. Perrin) using the 50CORON aperture (λ_c = 0.58 μm, pixel scale: 50.72 mas pixel⁻¹; Riley et al. 2018), totaling 76 frames. Each visit was performed with a different telescope orientation, with position angles of −8423, −6123, −223, and 2559 for the y-axis in the images (N to E). These telescope orientations are selected to obtain 360° azimuthal coverage of the disk down to the occulting mask (similar to Schneider et al. 2014). In each visit, we first obtained 16 short 32 s exposures on the WEDGEA0.6 position to probe the inner region down to a half-width of 03. Then we obtained three longer 483.3 s exposures on the WEDGEA1.0 position (half-width: 05) to deeply probe the exterior region of the disk.

To perform PSF subtraction, we also observed the reference star HD 196081 (selected for color, brightness, and on-sky proximity matches to HD 191089⁴⁸ ) using the same aperture positions, with one visit interleaved between the two science (HD 191089) visits at each epoch. With the numerical mask created by Debes et al. (2017) to mask out the STIS occulters, we perform multiple PSF subtractions using different approaches. We first apply the classical reference differential imaging (cRDI) technique to subtract the starlight in each exposure by minimizing the residual variation in the region of the coronagraphically unapodized diffraction spikes (excluding where the disk resides). We also obtain the STIS reference exposures from the STIS PSF archive created in Ren et al. (2017) for MRDI reduction: for each target exposure, we perform NMF reduction using the 10% most-correlated references in the STIS archive (i.e., 85 images for WEDGEA1.0 and 45 images for WEDGEA0.6). The final STIS disk image is then the elementwise mean of the individual derotated PSF-subtracted exposures.

To calibrate the STIS image in physical surface brightness units, we convert the PHOTFLAM photometric parameter in the raw FITS file for the HD 191089 observations (F_λ = 4.15 ×10⁻¹⁹ erg cm⁻² Å⁻¹ count⁻¹) to F_ν = 4.56 × 10⁻⁷ Jy s count⁻¹ using the conversion equation in Appendix B.2.1 of Viana et al. (2009),

$$\begin{eqnarray*}&&{F}_{\nu }=\displaystyle \frac{{\lambda }_{{\rm{c}}}^{2}{F}_{\lambda }}{3\times {10}^{-13}},\end{eqnarray*}$$

where λ_c = 0.58 μm is the pivot wavelength of STIS. We then multiply our combined STIS image in count s⁻¹ pixel⁻¹ by F_ν and divide it by the STIS pixel area on sky to obtain the surface brightness data in units of Jy arcsec⁻².

2.1.3. Gemini/GPI (2015)

On 2015 September 1, HD 191089 was observed using Gemini/GPI in H-band (λ_c = 1.65 μm, pixel scale: 14.166 ± 0.007 mas pixel⁻¹; De Rosa et al. 2015) polarimetric mode (H-Pol; Perrin et al. 2015) during the Gemini Planet Imager Exoplanet Survey (PI: B. Macintosh; Macintosh et al. 2014). We took 28 exposures, each with 88.74 s integration time, with a total field rotation of 1022. The airmass ranged from 1.002 to 1.01, the differential motion image monitoring seeing measurement was 116 ± 017, and the multi-aperture scintillation sensor seeing was 11 ± 03.

To obtain the Stokes cube ({I, Q, U, V}) for the HD 191089 debris disk, we followed the recipes described in Perrin et al. (2014) and Millar-Blanchaer et al. (2015) and reduced the raw exposures using the GPI data reduction pipeline (DRP; Perrin et al. 2014, 2016) and the automated data processing architecture (Data Cruncher; Wang et al. 2018). The Q and U components in the traditional Stokes cube were then transformed to the local Stokes cube ({${{ \mathcal Q }}_{\phi }$, ${{ \mathcal U }}_{\phi }$}), with ${{ \mathcal Q }}_{\phi }$ representing the polarized light perpendicular or parallel to the radial direction and ${{ \mathcal U }}_{\phi }$ at ±45° from it (Monnier et al. 2019). On the local Stokes maps of HD 191089, we notice two similar low spatial frequency octopole structures with a rotation of ∼45°. However, for a ${{ \mathcal U }}_{\phi }$ map, we do not expect any signal from an optically thin disk with single scattering events on the dust (see Canovas et al. 2015 for a discussion of the validity and exceptions). Given that we observe similar structures in the other GPI polarimetry observations, we expect such a structure to be one of the instrumental artifacts (T. E. Esposito et al. 2019, in preparation). To reduce the systematic errors induced by it, we fit an octopole model using the ${{ \mathcal U }}_{\phi }$ map and remove it from the ${{ \mathcal U }}_{\phi }$ map, then rotate the model 45° and remove it from the ${{ \mathcal Q }}_{\phi }$ map.

We flux-calibrate these data following the procedure described in Hung et al. (2015, 2016). For HD 191089, we first adopt an H-band flux of F_⋆ = 3.749 ± 0.119 Jy (Skrutskie et al. 2006), then we adopt the satellite-to-star ratio R = 2.035 × 10⁻⁴ from the GPI DRP (Wang et al. 2014; Perrin et al. 2014) and average satellite spot total flux $S=(1.12\pm 0.17)\times {10}^{3}\,\mathrm{count}\,{{\rm{s}}}^{-1}$ in the FITS file header. Combining these, we obtain a conversion factor of

$$\begin{eqnarray*}&&{F}_{\nu }=\displaystyle \frac{{{RF}}_{\star }}{S}=(6.87\pm 1.06)\times {10}^{-7}\,\mathrm{Jy}\,\ {\rm{s}}\,\ {\mathrm{count}}^{-1}.\end{eqnarray*}$$

We apply that conversion factor and normalize the local Stokes maps by the exposure time and pixel area to obtain the disk surface brightness data in units of Jy arcsec⁻². The flux-calibrated data are then geometrically corrected and smoothed by convolving with a Gaussian kernel (σ = 14.166 mas, i.e., the scale of 1 GPI pixel; Millar-Blanchaer et al. 2016) to remove the high spatial frequency noise that impacts regions smaller than the Nyquist-sampled PSF of GPI.

Based on different PSF subtraction methods for each data set, we estimate the uncertainties as follows.

NMF and cRDI (STIS, NICMOS). We estimate the noise in these reductions from the ensemble of science frames to probe the temporal variations from frame to frame (16 NICMOS frames, 64 and 12 STIS frames at the two different positions on the mask; see Table 2). We proceed as follows. After subtraction of the PSF, we compute the pixelwise standard deviation across the science frames to obtain the typical noise map per frame, which is used to account for the noise added by PSF subtraction. We then replicate this noise map for N_frame times and derotate each with the same angle as each science frame. We obtain the final noise map by computing the square root of the quadratic sum of these derotated noise maps. See Figure 2 for the S/N maps.

Zoom In Zoom Out Reset image size

KLIP (NICMOS). For the NICMOS-KLIP reduction, we use the ALICE library of reference stars that are processed identically to HD 191089 to estimate the residual speckle noise. For each HD 191089 image (16 total), we first select 25% of the most-correlated images in the reference library (i.e., 212 reference images). In this way, we obtain 413 reference images that are correlated with at least one HD 191089 image (i.e., ∼50% of the entire library). We split these 413 images into 25 groups, each containing 16 images (with the leftovers randomly discarded). In each group, the 16 images are treated as the mock target images to simulate 16 HD 191089 nondetection images. For one mock target image, we first identify the real target that was observed, then we remove the images that are taken on the same real target from the reference library to avoid self-subtraction. We then use the updated reference library to perform KLIP subtraction of the mock image following the identical procedure as for HD 191089 (i.e., 19 eigenmodes from the 84 most-correlated references). The 16 reduced mock images are then derotated using the same orientation angles as the ones in the HD 191089 observations. For each group, we take the elementwise mean of the 16 rotated reduced mock target images as the mock result for one realization of HD 191089 nondetection. We then obtain a total of 25 realizations of nondetections using the 25 full groups from the 413 mock images. We take the elementwise standard deviation from the 25 mock nondetections as the noise map for our NICMOS-KLIP reduction.

Polarimetry (GPI). For the GPI ${{ \mathcal Q }}_{\phi }$ map, we used the convolved ${{ \mathcal U }}_{\phi }$ image as a proxy for the uncertainty (Millar-Blanchaer et al. 2015). To obtain the noise map, for each angular separation to the star, we calculate the standard deviation of an annulus with 3 pixel width in the convolved ${{ \mathcal U }}_{\phi }$ image as the noise. This noise map provides a reasonable estimation of PSF subtraction residuals, photon and detector noise, and residual instrumental polarization (Millar-Blanchaer et al. 2015). See Figure 3 for the convolved ${{ \mathcal Q }}_{\phi }$ and ${{ \mathcal U }}_{\phi }$ images used for signal-to-noise ratio (S/N) calculation.

Zoom In Zoom Out Reset image size

The HD 191089 debris disk is detected in all of our HST and GPI observations. Here we compare the quality of the different reductions and discuss their relative merits for measuring different quantities of interest.

STIS. The disk is detected in the STIS data with high morphological and photometric fidelity in both the cRDI and NMF reductions (Figure 4). The system shows a bright parent belt, surrounded by a faint and diffuse halo detected up to ∼6'' from the star. These features are detected with a dynamic range of 3 orders of magnitude, with a peak ring surface brightness of ∼1 mJy arcsec⁻² and the diffuse halo detected down to a few μJy arcsec⁻². Theoretically, unlike KLIP reduction, both the cRDI and NMF reduction methods are free from oversubtraction caused by overfitting disk features using the references. This was discussed in Ren et al. (2018b) and is confirmed in Ren et al. (2018a), where the NMF method is able to successfully retrieve the spiral arms for the MWC 758 system.

Zoom In Zoom Out Reset image size

For HD 191089, the bright disk and excellent cRDI quality enable a quantitative comparison between cRDI and NMF. In Figure 4, we present the reduction results from different methods for both the STIS and NICMOS data for a comparison between the methods. For the STIS data, the ring surface brightness is consistent within ±10% between cRDI and NMF. The NMF achieves higher S/Ns along the major and minor axes of the disk by a factor of ∼3, i.e., the regions containing coronagraphically unapodized diffraction spikes, which are marked by red ellipses in Figure 5.

Zoom In Zoom Out Reset image size

NICMOS. The halo component discovered in STIS is confirmed in the NICMOS-NMF image. The NICMOS-KLIP image in Figure 4 is able to recover the halo that was not observed in the original discovery image (Soummer et al. 2014). The new reduction is obtained with a larger field of view and fewer KLIP components. In retrospect, it is not surprising that the original reduction by Soummer et al. (2014) did not detect the halo: the KLIP method is based on principal component analysis, which requires mean subtraction for each individual image. When modeling the target with KLIP, the reduced image also has zero mean, which offsets the faint halo as negative background; thus, the halo cannot be recovered with positive signals. In addition, to maximize the removal of the starlight, a large number of KLIP components was used in Soummer et al. (2014). In this case, KLIP also removed the halo because its extended diffuse structure resembles PSF wing.

The NMF independently confirms the existence of the halo in the NICMOS data while recovering the ring regardless of increasing component number. The NMF component basis is nonnegative—thus, not orthogonal—and it does not perform direct projection of vectors as KLIP, which falls into the overfitting regime, but instead searches for the best nonnegative combination of nonnegative NMF components. The halo does not resemble the NMF components that are used to model the PSF wings, and thus it remains in the residual after PSF subtraction.

In the NICMOS results in Figure 4, the NIMCOS-NMF image has about two times the ring surface brightness of the NICMOS-KLIP image, supporting the expected behavior of the two methods: KLIP's overfitting because of direct vector projection cannot be avoided, even with a small number of eigenmodes, and the mean-subtraction offset reduces the overall flux of the system. In contrast, NMF is expected to preserve the surface brightness of the NICMOS disk as for the STIS data.

GPI in total intensity. We attempted to detect the disk in total intensity from this same data set with MRDI: for each polarization-direction pair of H-Pol exposures, we derive a single total intensity image. For starlight subtraction, we select the 100 most-correlated GPI H-band PSFs from the library of 15,847 exposures, then perform KLIP and NMF reductions. We present the S/N maps of the null detections in Figure 5: we do not detect the disk or any point source using either KLIP or NMF. We hypothesize that a larger or better reference library is needed to obtain the best match for HD 191089.

To estimate an upper limit on the surface brightness for the disk, we calculate the noise map for the GPI-KLIP image using the standard deviation across the individual reduced images. We estimate a 1σ uncertainty of ∼3 mJy arcsec⁻² at the minor axis and ∼1 mJy arcsec⁻² along the major axis. Assuming gray scattering for the dust in the ring seen with the NICMOS F110W filter and GPI H band, the NICMOS-NMF result will produce surface brightnesses of ∼2.2 and ∼1.0 mJy arcsec⁻², respectively. Therefore, even the disk is not removed by the reduction methods; it is below a detection threshold of 1σ (which corresponds to a total S/N of 5 for extended structure spanning ∼25 pixels; Debes et al. 2019). In this way, its detectability is beyond the limit of the current methods with the current GPI H-band PSF library.

For the ring, the GPI ${{ \mathcal Q }}_{\phi }$ map provides the highest spatial sampling, with STIS and NICMOS offering total intensity observations at different wavelengths. Using the GPI ${{ \mathcal Q }}_{\phi }$ map, we obtain the geometric structure for the ring in Section 3.3.1, including inclination, semimajor axis, position angle, ring center position, and eccentricity of the deprojected ring. With the geometric information and the parallactic distance to HD 191089, we are able to calculate the average intensity of the disk as a function of physical distance to the host star (i.e., the radial profile).

For the halo, the STIS total intensity image is able to probe the largest spatial extent with a larger field of view and high sensitivity. Assuming that the halo is coplanar with the ring, we can measure the radial profile for the halo, which will help to identify whether the halo is a geometric extension of the ring by comparing their outer power-law indices. We use the cRDI reduction for this purpose, since it covers a larger field of view than the NMF reduction obtained from the fixed-width STIS archive in Ren et al. (2017).

For both the ring and the halo, if we assume that they are coplanar, then we can measure the disk surface brightness as a function of scattering angle—the SPF. The SPF is related to properties of the dust and therefore provides insights into the composition, size distribution, and minimum dust size for the system, even though this information is degenerate. To extract the information for the dust, we adopt the spatial distribution information from GPI ${{ \mathcal Q }}_{\phi }$ measurements, then use radiative transfer modeling tools to model the ring. Given that the ring is well resolved with all three instruments, the SPFs at different wavelengths are expected to help constrain the dust size and composition.

The spatial distribution of the dust in the debris disk can be parameterized using a three-dimensional function in cylindrical coordinates: the radial distribution in the midplane and the vertical distribution along the axis perpendicular to the midplane. Along the radial direction, the dust in the system follows a combination of two power laws,

$$\begin{eqnarray}&&\rho (r)\propto {\left\{{\left(\displaystyle \frac{r}{{r}_{{\rm{c}}}}\right)}^{-2{\alpha }_{\mathrm{in}}}+{\left(\displaystyle \frac{r}{{r}_{{\rm{c}}}}\right)}^{-2{\alpha }_{\mathrm{out}}}\right\}}^{-\tfrac{1}{2}},\end{eqnarray} \tag{ 1 }$$

where α_in > 0 and α_out < 0 approximate the midplane dust density power-law indices interior and exterior to r = r_c (Augereau et al. 1999). Along the normal direction of the disk midplane, the dust follows a Gaussian dispersion form, i.e.,

$$\begin{eqnarray}&&Z(r,z)\propto \exp \left[-{\left(\displaystyle \frac{z}{\zeta (r)}\right)}^{2}\right],\end{eqnarray} \tag{ 2 }$$

with

$$\begin{eqnarray}&&\zeta (r)={{hr}}^{\beta },\end{eqnarray} \tag{ 3 }$$

where

$$\begin{eqnarray}&&\beta =1\end{eqnarray} \tag{ 4 }$$

for a nonflared debris disk. For the HD 191089 system, h = 0.04 is adopted from the vertical structure study (Thébault 2009).

The power-law index for the surface density radial profile is then (Augereau et al. 1999)

$$\begin{eqnarray}&&{\rm{\Gamma }}=\alpha +\beta ,\end{eqnarray} \tag{ 5 }$$

where α is α_in or α_out in Equation (1), with Γ being Γ_in or Γ_out as the corresponding indices.

Since illumination decreases as a function of distance from the star, our images must be corrected for illumination effects before the surface density can be measured. The relationship between the surface brightness power-law index, γ, and the surface density power-law index, Γ, is

$$\begin{eqnarray}&&\gamma ={\rm{\Gamma }}-2,\end{eqnarray} \tag{ 6 }$$

where γ is γ_in or γ_out, corresponding with Γ being Γ_in or Γ_out.

The radial distribution of the system in Equation (1) is supplemented with two extra parameters, r_in and r_out, which are the cutoff radii. They are introduced to describe the clearing of materials interior and exterior to the disk, i.e., when r < r_in or r > r_out, ρ(r) = 0.

3.3.1. Ellipse Parameters

3.3.2. Radial Distribution

We measure the surface brightness power-law indices by averaging the flux density at the same stellocentric radius in the highest-S/N regions. Although this approach implicitly assumes isotropic scattering, which is a likely incorrect assumption, we tested it by performing the fit in narrow wedges (e.g., along the minor or major axis) and obtained consistent results, indicating that the scattering anisotropy does not change our results significantly. We plot the surface brightness radial profiles for the ring and halo in Figure 6 and summarize the results in Table 4. For both the ring and the halo, we define the inner and outer cutoff radii, r_in and r_out, as the radii for which the corresponding average surface brightness is consistent with zero at the 1σ level.

Zoom In Zoom Out Reset image size

Table 4.  Disk Morphology Parameters

Parameter	Ring	Halo	Meaning
Instrument	GPI	STIS
θinca	${{59}^{^\circ }}_{-{2}^{^\circ }}^{+{4}^{^\circ }}$	59°b	Disk midplane inclination from face-on
θPAa	${{70}^{^\circ }}_{-{3}^{^\circ }}^{+{4}^{^\circ }}$	70°b	Position angle for major axis (N to E)
γinc	4.9 ± 0.2	⋯d	Surface brightness power-law indices
γoutc	−6.1 ± 0.2	−2.68 ± 0.04
rcc	43.6 ± 0.3 au	⋯d	rc in Equation (1)
rinc	26±4 au	⋯d	Inner and outer clearing radii
routc	78±14 au	640±130 au
rcenter	45.6 ± 0.2 au	⋯d	Peak position and FWHM of a Gaussian ring
FWHM	24.9 ± 0.4 au	⋯d

Notes.

^aFitted with Debris Ring Analyzer (Stark et al. 2014). ^bValues adopted from the GPI results. ^cFitted for Equations (1) and (6). ^dUnconstrained from the STIS image.

Download table as:  ASCIITypeset image

We measure the following parameters for the ring from the GPI ${{ \mathcal Q }}_{\phi }$ observation: γ_in = 4.9 ± 0.2, r_c = 43.6 ± 0.3 au, and γ_out = −6.1 ± 0.2.⁵⁰ We also fit a Gaussian ring to the deprojected ${{ \mathcal Q }}_{\phi }$ surface density map, and the ring is centered at r_center = 45.6 ± 0.2 au with a 24.9 ± 0.4 au FWHM. We do not report the results from STIS or NICMOS, since they have low spatial sampling and are noisy near the corresponding IWAs.

We measure the geometric parameters for the halo from the STIS data because they cover the largest field of view. We measure a surface brightness power-law index of γ_out = −2.68 ± 0.04. The power-law indices measured at different position angles are also consistent with the integrated measurement within 1σ. Therefore, we report the integrated radial profiles to reduce systematic uncertainty.

In the above calculation, we have assumed a flat disk (i.e., h ≪ 1); however, the disk is not perfectly flat, and the line of sight passes through different radii at different heights. Under this scenario, for a nonflared disk, we calculate that the radial separation will be modified by multiplicative factors of

$$\begin{eqnarray}&&{\left(1\pm 0.5h\tan {\theta }_{\mathrm{inc}}\right)}^{-1}.\end{eqnarray} \tag{ 7 }$$

For the HD 191089 system, these line-of-sight intersections increase the uncertainties by ∼3% for h = 0.04. Therefore, the approximation of a flat disk will not bias the results for a disk with small scale height. For a continuous vertical distribution of the dust that is following a Gaussian decay, this effect is then an upper limit, since most of the scatterers are close to the midplane. We thus ignore this effect, given its model-dependent minor impact on the uncertainties.

We find that the radial power-law indices, γ_out, for the ring and halo differ by >10σ. Despite the physical connection between the ring and the halo, although the two power-law indices are measured at different wavelengths, when the radial distribution of the dust is wavelength-independent, the indices are indicating that the halo is not a geometrical extension of the ring's outer part.

Hughes et al. (2018) summarized the SPFs for different systems including zodiacal dust and debris disks and found a tentative universal SPF trend for the dust in debris disks. To further investigate the similarities and differences of SPFs in different systems, we first derive the uncertainty for the scattering angles and the SPFs in Appendix A, then measure the empirical SPFs for the ring and halo for our HD 191089 observations.

For the ring and halo, we present in Figure 7 the SPFs averaged for both sides (i.e., the NE and SW sides) using STIS-cRDI and NICMOS-NMF results to minimize oversubtraction. We observe different trends of the SPFs between the ring and the halo: in the STIS data, the halo is more forward- and backward-scattering than the ring; in the halo, the backward scattering is less strong than the forward trend. In addition, for the phase functions of the ring, the GPI ${{ \mathcal Q }}_{\phi }$ polarized-light image and STIS and NICMOS total intensity images have similar trends.

Zoom In Zoom Out Reset image size

We measure the polarization fraction for the ring using the GPI ${{ \mathcal Q }}_{\phi }$ and NICMOS-NMF images. Given the fact that we cannot recover the ring in total intensity with GPI H-band observations at 1.65 μm (Section 2.3), we instead use the 1.12 μm NICMOS-NMF observation that is both close to the GPI wavelengths and has less of an oversubtraction effect. The polarization fractions derived from GPI ${{ \mathcal Q }}_{\phi }$ and NICMOS-NMF images are around 20%–40%, with no clear trend (Figure 8). The polarization fraction is possibly peaking at ∼110°; however, it cannot be as firmly established as in previous measurements (e.g., Perrin et al. 2014; Frattin et al. 2019; Milli et al. 2019). To better constrain the polarized fraction values, we need H-band total intensity observations to rule out the wavelength-dependent effect.

Zoom In Zoom Out Reset image size

To compare SPFs in the STIS image, we first normalize the SPFs by dividing their average surface brightness at 90° ± 10° scattering angle. We then divide the normalized halo SPF by that of the ring to illustrate the difference. In Figure 9, the halo is likely both more forward- and backward-scattering in the probed scattering angles.

Zoom In Zoom Out Reset image size

To compare the HD 191089 SPFs with the ones in the literature, we present the normalized SPFs in linear scale for selected samples including both solar system objects (Saturn's D68 and G rings; Hedman & Stark 2015) and circumstellar disk systems (HD 181327: Stark et al. 2014; HR 4796 A: Milli et al. 2017) in Figure 10. Comparing with the previous studies, the HD 191089 ring SPF lies between the Saturn rings and the other samples, while its halo SPF lies above the Saturn rings.

Zoom In Zoom Out Reset image size

Using the HST disk images, we calculate the color of the dust as follows. We first bin the 1.12 μm NICMOS-NMF image and the 0.58 μm STIS-NMF image to two images with a pixel size of ∼150 mas to reduce correlated noise. We then divide the binned images by the corresponding pysynphot (STScI Development Team 2013) NICMOS or STIS counts of a T = 6450 K blackbody to obtain the magnitude per pixel. We calculate the difference of the two magnitude maps and present the radial profile along the ansae of the system in Figure 11.

Zoom In Zoom Out Reset image size

For the ring, the dust scatters ∼25% more light (Δmag ≈0.25) within the NICMOS F110W passband than in the STIS 50CCD, showing a red scattering property. For the halo, Δmag ≈ −1. Dust in the halo is expected to be generated in the ring through collisional cascade, then the smaller dust that is more sensitive to radiation pressure migrates outward to form the halo (e.g., Strubbe & Chiang 2006; Thébault & Wu 2008). Assuming that scattered-light images primarily probe the cross sections of the dust whose sizes are comparable to the observing wavelength, this red-to-blue trend from the ring to the halo in Figure 11 is consistent with this scenario.

We model the HD 191089 ring with the MCFOST (Pinte et al. 2006, 2009) radiative transfer modeling code and describe the dust using a distribution of hollow spheres (DHS; Min et al. 2003, 2005, 2007). As a derivative of the Mie theory, where the dust grains are assumed to be spherical (Mie 1908), DHS adopts vacuum centers for these spherical dust grains. To approximate small irregularly shaped dust, the only additional parameter in DHS from Mie—maximum vacuum fraction (f_max)—parameterizes the central vacuum fraction in the dust that is uniformly ranging from zero to f_max. The DHS has been used to successfully reproduce the SPF of linearly polarized scattered light with incident unpolarized light for quartz particles in the laboratory (Min et al. 2005), characterize the spectral features of the interstellar medium (e.g., Min et al. 2007; Poteet et al. 2015), and better fit the scattering properties of the HR 4796A circumstellar disk system (Milli et al. 2015).

In our study, we adopt the geometrical parameters derived in Section 3.2 for the ring and perform radiative transfer modeling with MCFOST using the DHS theory to probe the dust properties (e.g., dust mass, minimum dust size, composition, porosity, and maximum void fraction). Due to the complex nature of the constituent dust, DHS is a computationally intensive technique, and we thus adopt parallel computation in the Python environment and use the DebrisDiskFM package (Ren & Perrin 2018).⁵¹ The framework is developed to efficiently explore debris disk properties through radiative transfer modeling.

DebrisDiskFM is based on two software codes. We use MCFOST (version 3.0.33) to generate disk model images using given input parameters. We use emcee (version 3.0rc1; Foreman-Mackey et al. 2013), which makes use of a Markov chain Monte Carlo (MCMC) strategy with affine invariant ensemble samplers (Goodman & Weare 2010), to obtain the posterior distributions for these parameters. With the two software codes, we use DebrisDiskFM to distribute posterior calculations among multiple computation nodes in a computer cluster,⁵² with each node calculating its own MCFOST models in parallel. Specifically, for each combination of input parameters, we store the model in a unique folder, with the folder named by the hashed string for the array of input parameters. We also append the folder name with a hashed random number, preventing multiple nodes from simultaneously accessing the same folder and causing errors.

We model the STIS-cRDI and GPI ${{ \mathcal Q }}_{\phi }$ rings, since they cover the largest wavelength range and have higher data quality. To study the dust properties, we assume the dust is made of different types of grains, each with a pure composition. We adopt the three compositions in Esposito et al. (2018): the amorphous silicate dust (i.e., "astronomical silicates," denoted by "Si"; Draine & Lee 1984), amorphous carbonaceous dust (denoted by "C"; Rouleau & Martin 1991), and H₂O-dominated ice described in Li & Greenberg (1998) to model the β Pic disk (denoted by "ice"). The size of the dust, a, follows a power-law distribution with index q, i.e.,

$$\begin{eqnarray}&&{dN}(a)\propto {a}^{-q}\,{da}.\end{eqnarray} \tag{ 8 }$$

The distribution is truncated at a minimum size of a_min, and we set a maximum limit of the dust size, a_max = 1000 μm. For the HD 191089 system, we set q = 3.5 for the expected dust size distribution for debris disks undergoing collisional cascade (e.g., theory and simulation: Dohnanyi 1969; Pan & Schlichting 2012; observation: MacGregor et al. 2016; Esposito et al. 2018).

In this paper, we probe the following seven parameters of interest through disk modeling:

• 1.  
  disk mass, M_disk, which generally controls the overall brightness of the disk at different wavelengths;
• 2.  
  porosity;
• 3.  
  mass fraction for "astronomical silicates," f(Si);
• 4.  
  mass fraction for amorphous carbonaceous dust, f(C);
• 5.  
  mass fraction for ice, f(ice);
• 6.  
  minimum dust size, a_min; and
• 7.  
  maximum void fraction for DHS, f_max.

Given the fact that the composition parameters are interconnected, i.e., f(Si) + f(C) + f(ice) = 1, we only explicitly sample f(Si) and f(C). We also set the lower limit for a_min to be 0.5 μm based on the calculation of the blowout sizes for different dust compositions by Arnold et al. (2019). In the implementation of MCMC modeling of the system using the DebrisDiskFM framework, for a given set of parameters, we first generate two parameter files for MCFOST to represent the spatial sampling and field of view for the three instruments. We then perform radiative transfer modeling with MCFOST, using the DHS theory, to calculate the images for the three instruments at their central wavelengths in Table 2.

To simulate instrument responses, we convolve the disk models with TinyTim PSFs (Krist et al. 2011)⁵³ for STIS and a two-dimensional Gaussian profile for GPI (FWHM = 53.8 mas, corresponding to 3.8 times the pixel scale for GPI to match the GPI PSF; Esposito et al. 2018).

We compare the STIS model directly with the STIS-cRDI image. We compare the GPI model with the GPI ${{ \mathcal Q }}_{\phi }$ image by first converting the Stokes Q and U models to a ${{ \mathcal Q }}_{\phi }$ model, then comparing the PSF-convolved ${{ \mathcal Q }}_{\phi }$ model with the observation. With the models and observations, we maximize the following log-likelihood function:

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}{ \mathcal L }\left(\theta | {X}_{\mathrm{obs}}\right) & = & -\displaystyle \frac{1}{2}\displaystyle \sum _{i=1}^{N}{\left(\displaystyle \frac{{X}_{\mathrm{obs},i}-{X}_{\mathrm{model},i}}{{\sigma }_{\mathrm{obs},i}}\right)}^{2}\\ & & -\,\displaystyle \sum _{i=1}^{N}\mathrm{log}{\sigma }_{\mathrm{obs},i}-\displaystyle \frac{N}{2}\mathrm{log}(2\pi ).\end{array}\end{eqnarray} \tag{ 9 }$$

In the above equation, X is a flattened image with N pixels. Subscripts "obs" and "model" denote observation and model, respectively, and σ_obs,i is the uncertainty for X_obs,i at the ith pixel. We only focus on the disk region (i.e., between r_in and r_out determined by GPI ${{ \mathcal Q }}_{\phi }$ image) to minimize the influence from the halo. In this paper, we assume that the pixel noise follows a Gaussian distribution, and that all pixels are independent of each other (see, e.g., Wolff et al. 2017, for a proper treatment of correlated pixels).

4.2.1. GPI ${{ \mathcal Q }}_{\phi }$ Image

One-step GPI. While the ring is resolved with all three instruments, we first only fit the GPI ${{ \mathcal Q }}_{\phi }$ image to investigate the dust properties. Using the flat priors presented in Table 5, we assign 60 chains to explore the parameter space and run the MCMC modeling procedure for 9000 steps. We discard the first 2000 steps that are identified as the burn-in stage and calculate the posterior distributions using the last 7000 steps (a total of 4.2 × 10⁵ models).

Table 5.  Independent Dust Properties Retrieved from Ring Image Modeling

Parametera	Priorb	Posteriorc
Image		GPI ${{ \mathcal Q }}_{\phi }$	STIS
${\mathrm{log}}_{10}{M}_{\mathrm{disk}}$ (M⊙)	(−12, −4)	$-{7.462}_{-0.010}^{+0.015}$	$-{6.79}_{-0.04}^{+0.03}$
Porosity	(0, 1)	${0.600}_{-0.008}^{+0.012}$	${0.01}_{-0.01}^{+0.01}$
f(Si)	(0, 1)	${0.005}_{-0.005}^{+0.016}$	${0.50}_{-0.05}^{+0.04}$
f(C)	(0, 1)	${0.004}_{-0.004}^{+0.011}$	${0.50}_{-0.04}^{+0.05}$
f(ice)d		${0.989}_{-0.010}^{+0.009}$	0.003e
${\mathrm{log}}_{10}{a}_{\min }$ (μm)	(−0.3, 2)	$-{0.115}_{-0.007}^{+0.013}$	${0.24}_{-0.03}^{+0.01}$
fmax	(0, 1)	${0.213}_{-0.014}^{+0.026}$	0.001e
${\chi }_{\nu }^{2}$ (GPI best fit)		1.05	240
${\chi }_{\nu }^{2}$ (STIS best fit)		12	10
q	3.5f
${\mathrm{log}}_{10}{a}_{\max }$ (μm)	3f
Iteration number	9000f
Burn-in	2000f

Notes.

^aThe morphological parameters in Table 4 are adopted. ^bThe parameters are limited to three decimal digits, with uniform sampling in the prior range. For M_disk and a_min, they are log-uniformly sampled. ^c16th, 50th, and 84th percentiles. ^dThe mass fraction for ice is $f(\mathrm{ice})=1-f(\mathrm{Si})-f({\rm{C}})$; thus, only Si and C are explicitly sampled. ^e95th percentile. ^fKept fixed.

Download table as:  ASCIITypeset image

Table 5 reports the credible intervals from the posterior distributions in Figure 13 (generated by corner; Foreman-Mackey 2016). We notice the small uncertainties in the results; these uncertainties are underestimated, since the correlated noise is ignored in our likelihood function (e.g., Wolff et al. 2017). In addition, we also argue that these small uncertainties are also the results from the limited dust models and that various scenarios may lead to small uncertainties (e.g., spatial distribution asymmetries, non-power-law surface density distribution); thus, the inclusion of correlated noise may still not lead to physical interpretations of the retrieved parameters.

4.2.2. STIS Image

4.2.3. Implications

We present the results from independent GPI and STIS modelings and discuss the implications as follows.

SED. We generate the SEDs corresponding to the best-fit parameters for the STIS and GPI ${{ \mathcal Q }}_{\phi }$ images with MCFOST and compare them with the observed one in Figure 12. We notice that the SED from the STIS best fit overpredicts the emission from the system, while the SED from the GPI ${{ \mathcal Q }}_{\phi }$ best fit underpredicts the emission. Even though different sizes of dust dominate the SED and the scattered-light images, the inability to reproduce the observed SED adds more evidence that our study using DHS models cannot provide a consistent model satisfying all observables. For instance, the true dust albedo differs from the one predicted in DHS models.

Zoom In Zoom Out Reset image size

Mie versus DHS. Although we do not obtain a set of parameters that is able to explain both data sets, we can still focus on the f_max parameter (i.e., maximum void fraction for the DHS theory). We argue that it is the most profound parameter that is added to the DHS theory from the Mie theory; the f_max parameter approximates dust grains from spheres to aggregate-like structures (Min et al. 2003, 2005, 2007), and we believe it is one of the keys to understanding the scattered-light properties of the dust.

For the STIS image in total intensity, we obtain f_max ≈ 0. In DHS theory, this value corresponds to the Mie theory scenario. For the GPI ${{ \mathcal Q }}_{\phi }$ image in polarized light, we obtain f_max ≈ 0.21. In this scenario, the DHS theory is preferred to the Mie theory, and the f_max parameter is smaller than what is in the interstellar dust (e.g., 0.7; Min et al. 2007), suggesting that dust properties are different under various environments. Given that the two data sets differ in two fundamental ways (different wavelengths and total intensity/polarization observation), we cannot determine which is the root cause of this discrepancy. However, this discrepancy in the f_max parameter is still informative, since it is an approximation of the dust structure—the discrepancy indicates that neither Mie nor DHS is able to well approximate the structure of the dust seen in scattered-light images.

Distinct parameters. In the radiative transfer modeling of the scattered-light images, we retrieve distinct sets of parameters in Table 5 and Figure 13. Although the best fits are not statistically consistent with each other, we categorize the difference in the dust properties in the ring into two groups.

• 1.  
  Physically possible: minimum dust size. The retrieved minimum dust size is ∼0.8 μm for GPI ${{ \mathcal Q }}_{\phi }$ and ∼1.7 μm for STIS. Both values are within a factor of ∼2, and the inclusion of correlated noise in the likelihood function may resolve the discrepancy. The values are reasonably consistent with blowout size calculations, although we caution that these are themselves highly dependent on disk composition, porosity, and aggregate structure—the blowout size for dust can vary by an order of magnitude for different composition and porosity (see the Arnold et al. 2019 calculation for HD 181327, a star similar to HD 191089).
• 2.  
  Physically impossible: porosity and composition fraction. Specifically, these parameters take values in limited ranges (i.e., from 0% to 100%), and the discrepancy is currently at the ∼50% level. Even though the discrepancies for these parameters can be alleviated with larger uncertainties, their physical meaning is then uninformative. A possible solution is to increase the uncertainties to ∼50% (assuming the inclusion of correlated noise is able to achieve it); however, that would render these parameters meaningless, since the large uncertainties would not reject any values in the physically plausible values (i.e., from 0% to 100%).

We note that the above parameter values are based only on MCMC fitting results, and neither may be correct given the limitation of DHS or Mie models. Specifically, these models are optimized for spectral fitting (e.g., Min et al. 2007; Poteet et al. 2015), but neither model is able to properly constrain the composition from scattered-light images (e.g., Milli et al. 2019). See Section 5.2 for more discussion on the dust properties retrieved from radiative transfer modeling.

4.2.4. Other Attempts

In addition to the above modeling efforts, we have performed separate modeling attempts with loosened prior constraints.

For GPI. Set the lower limit in the prior for a_min to be 0.01 μm, and keep q unconstrained. We observed a steeper q ≈ 4.15 with a_min ≈ 0.02 μm to describe the GPI ${{ \mathcal Q }}_{\phi }$ image, but it still does not recover the STIS flux density. In addition, although this smaller dust is not blown out by radiation pressure, its collisional cascade suppliers (slightly larger dust) are blown out (e.g., Burns et al. 1979; Silsbee & Draine 2016; Arnold et al. 2019); thus, this scenario is not stable.

For GPI and STIS. Simultaneously model the STIS and GPI images with the conditions above (i.e., as in the previous bullet point). However, the best fits indicate a bimodal distribution, with one better recovering the STIS image and the other better recovering the GPI image. In the former model, the GPI ${{ \mathcal Q }}_{\phi }$ model displays negative polarization at the smallest scattering angles, failing to properly recover the observation; in the latter, the STIS model does not recover the ansae in the observed data.

Based on our modeling efforts, we conclude that the STIS and GPI data sets cannot be consistently reproduced with a single model. As a result, we present our separate models for the GPI ${{ \mathcal Q }}_{\phi }$ image and the STIS image in Figures 13 and 14. See Section 5.2 for more discussion on applying the DHS theory to disk modeling.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Although our radiative transfer modeling efforts cannot explain the scattering properties of the ring, we are confident in the results of our geometrical analysis, since they are based only on the surface brightness distribution of the system.

5.1.1. Ring Measurables

Ring clearing radii (r_in and r_out). Churcher et al. (2011) observed the ring at 18.3 μm, reporting a dust belt from 26 to 84 au⁵⁴ , which is consistent with our fitting results of r_in = 26 ± 4 au and r_out = 78 ± 14 au at the 1σ level. The position angle of the major axis, as well as the inclination, is better constrained with our high spatial resolution data in scattered light with GPI.

Brightness asymmetry. For the ring, we are not able to find brightness asymmetry beyond ∼10% or 1σ with the GPI and STIS data. Although a tentative ∼20% asymmetry was observed at the 1.8σ level in the 18.3 μm observation by Churcher et al. (2011), if the dust follows the same spatial distribution at these wavelengths, the 18.3 μm emission asymmetry is likely a result of statistical or instrumental fluctuation.

Planet perturber. In the GPI H-band total intensity HD 191089 observations, we did not detect any point source (Figure 5). We report 5σ point-source contrast limits of ∼1 × 10⁻⁵–∼3 × 10⁻⁶ between 03 and 08 with the forward-modeling planet detection method to correct for selfsubtraction and oversubtraction in Ruffio et al. (2017). Using these contrast limits, if a planet is shepherding the ring, for a system with an age of 22 Myr and using the evolution tracks in Spiegel & Burrows (2012), its mass is expected be smaller than ∼5 M_Jup.

Using the Morrison & Malhotra (2015) analysis for the outer edge of a planet's chaotic zone and assuming this outer edge is the inner edge of the ring of a debris disk, the upper limit on planet mass can be translated to the lower limit on the semimajor axis of the planet's orbit:

$$\begin{eqnarray}&&{a}_{{\rm{p}}}=\displaystyle \frac{{r}_{{\rm{i}}{\rm{n}}}}{1+1.7{\left(5{M}_{{\rm{J}}{\rm{u}}{\rm{p}}}/{M}_{\star }\right)}^{0.31}}.\end{eqnarray} \tag{ 10 }$$

When we substitute the measured r_in into the above equation, we obtain a lower limit of a_p = 20 ± 3 au for the semimajor axis of the planet.

5.1.2. Ring as an Exo–Kuiper Belt

The spatial extent of the ring (r_center = 45.6 ± 0.2 au, FWHM = 24.9 ± 0.4 au for a Gaussian ring) resembles that of the solar system's Kuiper Belt (30–50 au; Stern & Colwell 1997; Bannister et al. 2018). To investigate the scenario where the ring is an extrasolar version of the Kuiper Belt, which is perturbed by a corresponding planet (i.e., Neptune), we first deproject the S/N map of the GPI observation to a face-on view, then search for possible mean-motion resonance orbits for the inner and outer radii of the disk. Based on Kepler's third law, we search for the period ratios corresponding to different combinations of stellocentric radii in Figure 15, where the period ratios are mapped to simple fraction values that correspond to the strongest resonance orbits of a hypothetical potential well.

Zoom In Zoom Out Reset image size

We set the boundary ranges to be consistent with pixelwise S/N ≈ 1. Among different combinations of the strongest resonances, we obtain four pairs of 4:3 and 5:2 resonance orbits to resemble the extent of the solar system Neptunian resonance orbits in Chiang et al. (2003). For the other resonance pairs, their boundary ranges are too narrow to cover the GPI disk, these pairs cannot be resolved because of the limitation of instrumental spatial resolution, or these solutions are consistent with the ranges of the 4:3 and 5:2 pairs; therefore, they are not presented or analyzed in this paper.

With the four resonance radii pairs, we are able to compute the location of the potential well (i.e., 1:1 resonance) at r_(1:1) = 060 ± 002, corresponding to a stellocentric distance of r_(1:1) = 29.9 ± 1.2 au. The mean resonance radii and the hypothesized 1:1 orbit of the potential well are shown in Figure 15. To confirm these resonances, deeper high-resolution and high-S/N observations are needed to firmly establish the edges of the ring. If the 1:1 gravitational potential is caused by a planet, it is likely of small mass and requires the future LUVOIR or HabEx missions for observation.

5.1.3. Halo: Radial Distribution

Overall distribution. In the STIS data, the halo extends to r_out = 640 ± 130 au, with a surface density power-law index of Γ_out = −0.68 ± 0.04. This power-law distribution index for the surface density profile is shallower than −1.5, i.e., the classical expectation for the halo of debris disks (e.g., Strubbe & Chiang 2006; Thébault & Wu 2008). It is also shallower than the power-law index of −1, which is expected for the steady-state radial motion of the dust (Jewitt & Meech 1987). Using the dust size–free approximation in Jewitt & Meech (1987), a gravity-dominated slowdown of the radial motion of the dust corresponds to a power-law index of −0.5, and a constant acceleration results in an index of −1.5.

A power-law index of $-{0.68}_{-0.04}^{+0.04}$ between the two values (−0.5, −1.5) is thus caused by the joint effect from multiple force sources. In debris disks, the −1.5 power-law index has already taken into account both the slowdown from gravity and the acceleration from radiation pressure; the $-{0.68}_{-0.04}^{+0.04}$ power-law index is thus calling for additional slowdown sources, and the slowdown by the interstellar medium is a plausible candidate. In fact, a similar surface density profile has been observed in the outskirts of the HR 4796 A halo, which has an index of −0.7 and a large-scale structure that is strongly suggestive of interaction with the interstellar medium (Schneider et al. 2018).

Following the Jewitt & Meech (1987) analytic derivation relating surface density to the radial speed of the dust, we assume the size distribution of the dust is independent of its stellocentric distance. Under this assumption, if the dust has an outward radial speed of v(r) ∝ r^x, where r is the radial separation, then the surface density will be Γ(r) ∝ r^−x−1. In this way, the dust in the halo of HD 191089 has a radial speed of v(r) ∝ r^{−0.32 ± 0.04}.

Local distribution. The surface density power law of the HD 191089 halo deviates from the classical model at the >10σ level, calling for detailed investigation of the local variation of the halo. As an attempt to investigate the variation of the surface density distribution at different stellocentric distances, we compute the power-law indices of the surface density radial profile for the halo at different radii in Figure 16. Assuming the dust size distribution is independent of its stellocentric distance, we also present the Jewitt & Meech (1987) analytical derivation between radial speed and surface density power-law index under difference scenarios. At different stellocentric distances, we observe the following.

• 1.  
  Interior to ∼200 au, the radial speed of dust decreases as stellocentric distance increases. This indicates the decrease of the net inward force that slows down the outward motion of the dust.
• 2.  
  Between ∼200 and ∼300 au, the radial speed reaches a constant, then increases as stellocentric distance increases. At ∼300 au, the surface density power-law index reaches that for the classical model for the debris disk halo (e.g., Strubbe & Chiang 2006; Thébault & Wu 2008).
• 3.  
  Exterior to ∼300 au, albeit with large uncertainty, the radial speed marginally increases, then reaches a constant as stellocentric distance increases.

Zoom In Zoom Out Reset image size

Given the complex two-dimensional residual structure in the single SPF-corrected distribution of the halo (Section 5.1.4), we do not further discuss the trend of the local distribution power-law indices. In addition, our analysis is based on the assumption that dust size does not vary as a function of stellocentric separation; however, the assumption is invalid for collision-dominated debris disks (e.g., Strubbe & Chiang 2006; Thébault & Wu 2008). Therefore, a full dynamical modeling of the halo is needed to better explain the local variations of the halo.

5.1.4. Halo: Surface Density Variation

In the classical model for the debris disk halo, the dominant dust size for optical depth decreases as stellocentric distance increases (e.g., Strubbe & Chiang 2006; Thébault & Wu 2008). Using the STIS observations of the HD 181327 halo, Stark et al. (2014) found a consistency between the observed SPF change (under Mie theory) and the classical model. However, we do not find a clear SPF variation trend for the HD 191089 halo. To investigate the HD 191089 halo, we adopt the averaged halo SPF from measurement and explore the two-dimensional surface density variation for the HD 191089 halo.

To investigate the deviation of the scattering properties from the same SPF for the halo in STIS at different stellocentric radii, we first scale the whole halo by the surface brightness radial distribution as measured for Figure 6, thus eliminating both the distance-dependent illumination and radial density distribution factors. We then divide the image by the empirically averaged SPF for the halo in Figure 10 based on the scattering angles for each pixel. The scaled STIS image is then deprojected to a face-on view, rotated to align the major axis with the x-axis, and subtracted by the median to show the first-order deviation from an identical SPF in all of the halo. Based on the quality of the NICMOS-NMF data, only the γ_out = −2.68 correction is applied.

The two-dimensional deviation from one SPF for both the STIS-cRDI and NICMOS-NMF are shown in Figure 17. We observe overdensity regions in the NE and SW sides of the STIS data at the ∼25% level, with the STIS NE region likely matching the NICMOS NE overdensity region. The underdensity region to the NW region in the STIS image is likely influenced by the truncation of signal by STIS's Wedge B (Figure 2).

Zoom In Zoom Out Reset image size

Possible explanations for the deviation from a uniform SPF in the halo are as follows: (1) when the scattering properties of the dust are the same in the halo, the deviations are corresponding to local surface density variations; (2) when there is no density variation, the dust's scattering properties are different; or (3) both effects are jointly affecting the SPFs in the halo.

Stellar encounter. Based on the complicated structure of the halo, we investigate the scenario of whether the halo was created by stellar encounter events (e.g., De Rosa & Kalas 2019). In the current epoch, the star at a separation of 114 to the southwest of HD 191089 in the STIS field of view (partially seen at the bottom right corner in Figure 1(a)) is a background star. It is identified by Gaia DR2 with Source ID 6847146784384527872 at d = 1.06 ± 0.06 kpc (Gaia Collaboration et al. 2018); thus, it is not responsible for creating the halo.

To trace the positions of the nearby stars in the past, we retrieve 44 stars that are within 5 pc from HD 191089 using the Gaia DR2 archive and use the proper motions to linearly propagate the locations of the 44 stars in the past (∼0.5 Myr ago). We notice three stars that have the nearest projected approach from ∼1.2 to ∼1.5 pc, which happened between ∼0.2 and ∼0.35 Myr from now.⁵⁵ For star encounter events, the closest approaches are typically smaller than ∼200 au (e.g., Pfalzner 2003; Pfalzner et al. 2018), and distinct features such as spiral arms dissipate beyond ∼1000 yr (Pfalzner 2003). Therefore, if a star encounter event created the halo for the HD 191089 system, it should happen early in a cluster environment, where close encounters are more frequent (e.g., the solar system; Pfalzner et al. 2018). Under this mechanism, the halo would have been dissipated.

In this paper, based on previous efforts in debris disk modeling with DHS (e.g., Min et al. 2010; Milli et al. 2017), we adopted the DHS theory to model the dust in the HD 191089 ring for the GPI ${{ \mathcal Q }}_{\phi }$ image in polarized light and the STIS image in total intensity. We cannot yet interpret the disk images with one model across different instruments. Although we have considered a limited combination of compositions, they span the range from refractory (carbon) to pure ice and even void (through porosity). The models consider a broad range of refractive index that encompasses most standard astronomical compositions. Therefore, it is unlikely that our failure to find a good fit is solely due to not trying another composition.

In our modeling results, the STIS image favors Mie theory, and the GPI ${{ \mathcal Q }}_{\phi }$ image favors DHS theory with a maximum void fraction of ∼21%. The discrepancy in this parameter, which is the only additional parameter in DHS from Mie, indicates that the shape of the dust cannot be well approximated by either theory. In addition, although previous modeling attempts (e.g., Rodigas et al. 2015; Choquet et al. 2017) have encountered that simple models cannot reproduce a total intensity observation at multiple wavelengths, the analysis performed in this paper adds another dimension—polarization observation—to the complexity of disk modeling.

For some of the dust parameters (e.g., porosity, mass fraction of compositions), the discrepancies may be mathematically resolved using larger uncertainties by taking into account the correlated noise (e.g., Wolff et al. 2017). However, this resolution does not change the best-fit values, and in this way, it would make these parameters less informative, since they can only take values in a limited range (i.e., from 0% to 100%). The lack of meaning for these retrieved parameters further supports the fact that current models (i.e., Mie, DHS) cannot properly depict the debris disk images obtained at different wavelengths and with different observational techniques.

The advances in dust descriptions may solve the discrepancies in the radiative transfer modeling of debris disks, including using laboratory measurements such as the Amsterdam-Granada Light Scattering Database (Muñoz et al. 2012) or adopting advanced models for dust shapes and optical properties (e.g., discrete-dipole approximation: Purcell & Pennypacker 1973; Draine & Flatau 1994; Min et al. 2006; Rayleigh–Gans–Debye (RGD) theory: Sorensen 2001; the T-matrix method: Mishchenko et al. 1996; Gaussian random spheres: Muinonen et al. 1996; aggregation of small particles: Kempf et al. 1999; Tazaki et al. 2016; Tazaki & Tanaka 2018; Arnold et al. 2019). Given the more realistic descriptions of dust properties, the latter may help to resolve the discrepancies encountered with current dust models. For example, Tazaki et al. (2016) calculated the SPFs for aggregates using the RGD theory and found that backward scattering was underestimated in previous studies with simple models; Arnold et al. (2019) calculated blowout size for aggregates and found that the size can vary by as much as an order of magnitude for different particle models. Although these treatments may resolve the discrepancies, MCMC retrievals of dust properties for these advanced descriptions of dust are currently limited by computational power.

For the HD 191089 system, the contribution from the halo in the STIS data is also calling for a more complex structural model. In addition, if the halo is not coplanar with the ring, then the halo will bias the SPF of the ring and add another dimension of complexity to the problem.

To quantify the coordinate values, the x and y coordinates can be directly measured in the image of the system with (x, y) ± (δx, δy), where δ denotes the uncertainty of the parameters in this paper. In this section, given the other measured quantities of the system (i.e., inclination and position angle of the semimajor axis), we obtain the z coordinates for this system.

To determine the z coordinates, we first set up the mathematical representation of the system. Let O = (0, 0, 0) be the origin of the 3D Cartesian coordinate system, with O placed at the geometric center of the debris disk, and unit vector $\hat{x}=(1,0,0)$ pointing W, $\hat{y}=(0,1,0)$, N, $\hat{z}\,=(0,0,1)$, pointing toward the observer. Let ${\hat{n}}_{\mathrm{disk}}=(a,b,c)$ denote the unit normal vector of the midplane of the debris disk system. Then, all of the points on the disk midplane, which also contains the origin O, satisfy

$$\begin{eqnarray}&&{ax}+{by}+{cz}=0,\end{eqnarray} \tag{ 11 }$$

with a² + b² + c² = 1.⁵⁶

The inclination of the system, which is denoted by ${\theta }_{\mathrm{inc}}\pm \delta {\theta }_{\mathrm{inc}}\in [0^\circ ,90^\circ ]$ and defined as the dihedral angle between the disk midplane and the xOy-plane, satisfies

$$\begin{eqnarray*}&&\cos {\theta }_{\mathrm{inc}}=\displaystyle \frac{{\hat{n}}_{\mathrm{disk}}\cdot {\hat{n}}_{{xOy}}}{\sqrt{| | {\hat{n}}_{\mathrm{disk}}| {| }^{2}| | {\hat{n}}_{{xOy}}| {| }^{2}}},\end{eqnarray*}$$

where · is the dot product between two vectors, and ${\hat{n}}_{{xOy}}=\hat{z}$ is the unit normal vector for the xOy-plane. The above equation becomes

$$\begin{eqnarray}&&\cos {\theta }_{\mathrm{inc}}=\displaystyle \frac{(a,b,c)\cdot (0,0,1)}{\sqrt{({a}^{2}+{b}^{2}+{c}^{2})({0}^{2}+{0}^{2}+{1}^{2})}}=c.\end{eqnarray} \tag{ 12 }$$

The position angle of the system, which is denoted by ${\theta }_{\mathrm{PA}}\pm \delta {\theta }_{\mathrm{PA}}\in [0^\circ ,180^\circ ]$, is defined as the angle from N to the intersecting line between the system and the xOy-plane. For the intersecting line, it satisfies ax + by = 0, since the points on it are represented as (x, y, 0). Let the mathematical slope angle of the line be θ_slope, which is defined as the counterclockwise angle from $\hat{x}$ to the line. Then, we have the relationship between the mathematical slope angle and astronomical position angle,

$$\begin{eqnarray*}&&{\theta }_{\mathrm{PA}}={\theta }_{\mathrm{slope}}-90^\circ ,\end{eqnarray*}$$

where

$$\begin{eqnarray*}&&\tan {\theta }_{\mathrm{slope}}=\displaystyle \frac{{dy}}{{dx}}=-\displaystyle \frac{a}{b};\end{eqnarray*}$$

therefore, we have

$$\begin{eqnarray}&&\tan {\theta }_{\mathrm{PA}}=\tan ({\theta }_{\mathrm{slope}}-90^\circ )=-\cot {\theta }_{\mathrm{slope}}=\displaystyle \frac{b}{a}.\end{eqnarray} \tag{ 13 }$$

For the points on the disk midplane, we can substitute Equations (12) and (13) into Equation (11); then we have the z-coordinate of the points as

$$\begin{eqnarray}\begin{array}{rcl}z & = & -\displaystyle \frac{1}{c}({ax}+{by})\\ & = & -\sec {\theta }_{\mathrm{inc}}a\left(x+\displaystyle \frac{b}{a}y\right)\\ & = & \sec {\theta }_{\mathrm{inc}}\sin {\theta }_{\mathrm{inc}}\cos {\theta }_{\mathrm{PA}}(x+y\tan {\theta }_{\mathrm{PA}})\\ & = & \tan {\theta }_{\mathrm{inc}}(x\cos {\theta }_{\mathrm{PA}}+y\sin {\theta }_{\mathrm{PA}}).\end{array}\end{eqnarray} \tag{ 14 }$$

Assuming the measured parameters are independent, the squared uncertainty for z is therefore

$$\begin{eqnarray}\begin{array}{rcl}{\delta }^{2}z & = & \ {\delta }^{2}\left[\tan {\theta }_{\mathrm{inc}}(x\cos {\theta }_{\mathrm{PA}}+y\sin {\theta }_{\mathrm{PA}})\right]\\ & = & {\left(x\cos {\theta }_{\mathrm{PA}}+y\sin {\theta }_{\mathrm{PA}}\right)}^{2}{\delta }^{2}(\tan {\theta }_{\mathrm{inc}})\\ & & +\,{\tan }^{2}{\theta }_{\mathrm{inc}}{\delta }^{2}(x\cos {\theta }_{\mathrm{PA}}+y\sin {\theta }_{\mathrm{PA}})\\ & = & {\left(x\cos {\theta }_{\mathrm{PA}}+y\sin {\theta }_{\mathrm{PA}}\right)}^{2}{\sec }^{4}{\theta }_{\mathrm{inc}}{\delta }^{2}({\theta }_{\mathrm{inc}})\\ & & +\,{\tan }^{2}{\theta }_{\mathrm{inc}}\{[{x}^{2}{\sin }^{2}{\theta }_{\mathrm{PA}}{\delta }^{2}({\theta }_{\mathrm{PA}})+{\cos }^{2}{\theta }_{\mathrm{PA}}{\delta }^{2}(x)]\\ & & +\,[{y}^{2}{\cos }^{2}{\theta }_{\mathrm{PA}}{\delta }^{2}({\theta }_{\mathrm{PA}})+{\sin }^{2}{\theta }_{\mathrm{PA}}{\delta }^{2}(y)]\}\\ & = & {\tan }^{2}{\theta }_{\mathrm{inc}}{\cos }^{2}{\theta }_{\mathrm{PA}}{\delta }^{2}(x)+{\tan }^{2}{\theta }_{\mathrm{inc}}{\sin }^{2}{\theta }_{\mathrm{PA}}{\delta }^{2}(y)\\ & & +\,{\sec }^{4}{\theta }_{\mathrm{inc}}{\left(x\cos {\theta }_{\mathrm{PA}}+y\sin {\theta }_{\mathrm{PA}}\right)}^{2}{\delta }^{2}({\theta }_{\mathrm{inc}})\\ & & +\,{\tan }^{2}{\theta }_{\mathrm{inc}}({x}^{2}{\sin }^{2}{\theta }_{\mathrm{PA}}+{y}^{2}{\cos }^{2}{\theta }_{\mathrm{PA}}){\delta }^{2}({\theta }_{\mathrm{PA}}).\end{array}\end{eqnarray} \tag{ 15 }$$

Now, by combining the (x, y) coordinates and the inclination and position angle of the system, we can obtain z ± δz from Equation (14) and the square root of Equation (15). In this paper, for the HD 191089 disk, the 1σ uncertainties for x and y are estimated to be 0.33 pixel,⁵⁷ and the uncertainties for θ_inc and θ_PA are obtained from the GPI ${{ \mathcal Q }}_{\phi }$ image using the Debris Ring Analyzer package by Stark et al. (2014).

From the (x, y, z) coordinates of the dust in the debris disk system, and given the position of the star at (x₀, y₀, z₀) ± (δx₀, δy₀, δz₀), we can then measure the scattering angle for the photons. A photon, which is emitted from the star and then interacts with the material at (x, y, z), has an original direction of ${\boldsymbol{r}}=(x,y,z)-({x}_{0},{y}_{0},{z}_{0})$. The photon, when collected by the observer, has a final direction of $\hat{z}$. The scattering angle of this photon, which is defined as the angle between ${\boldsymbol{r}}$ and $\hat{z}$, is thus

$$\begin{eqnarray}\begin{array}{rcl}{\theta }_{\mathrm{scatter}} & = & \arccos \left(\displaystyle \frac{{\boldsymbol{r}}\cdot \hat{z}}{\sqrt{| | {\boldsymbol{r}}| {| }^{2}| | \hat{z}| {| }^{2}}}\right)\\ & = & \arccos \left[\displaystyle \frac{z-{z}_{0}}{\sqrt{{\left(x-{x}_{0}\right)}^{2}+{\left(y-{y}_{0}\right)}^{2}+{\left(z-{z}_{0}\right)}^{2}}}\right].\end{array}\end{eqnarray} \tag{ 16 }$$

Substituting Equation (14) into the above equation, we have the scattering angle at the (x, y) position in the detector frame (i.e., on the xOy-plane),

$$\begin{eqnarray}&&{\theta }_{\mathrm{scatter}}=\arccos \left\{\displaystyle \frac{\tan {\theta }_{\mathrm{inc}}(x\cos {\theta }_{\mathrm{PA}}+y\sin {\theta }_{\mathrm{PA}})-{z}_{0}}{\sqrt{{\left(x-{x}_{0}\right)}^{2}+{\left(y-{y}_{0}\right)}^{2}+{\left[\tan {\theta }_{\mathrm{inc}}(x\cos {\theta }_{\mathrm{PA}}+y\sin {\theta }_{\mathrm{PA}})-{z}_{0}\right]}^{2}}}\right\}.\end{eqnarray} \tag{ 17 }$$

Assuming the measured parameters are independent, the corresponding squared uncertainty for θ_scatter is then

$$\begin{eqnarray*}\begin{array}{rcl}{\delta }^{2}{\theta }_{\mathrm{scatter}} & = & \ \displaystyle \frac{1}{1-{\left[\tfrac{z-{z}_{0}}{\sqrt{{\left(x-{x}_{0}\right)}^{2}+{\left(y-{y}_{0}\right)}^{2}+{\left(z-{z}_{0}\right)}^{2}}}\right]}^{2}}\\ & & \times \,{\delta }^{2}\left[\displaystyle \frac{z-{z}_{0}}{\sqrt{{\left(x-{x}_{0}\right)}^{2}+{\left(y-{y}_{0}\right)}^{2}+{\left(z-{z}_{0}\right)}^{2}}}\right],\end{array}\end{eqnarray*}$$

denoting $v^{\prime} \equiv v-{v}_{0}$ for $v\in \{x,y,z\}$ and thus ${\delta }^{2}(v^{\prime} )={\delta }^{2}(v)+{\delta }^{2}({v}_{0})$; then, the above equation becomes

$$\begin{eqnarray}&&\begin{array}{l}{\delta }^{2}{\theta }_{\mathrm{scatter}}=\displaystyle \frac{{x}^{{\prime} 2}+{y}^{{\prime} 2}+{z}^{{\prime} 2}}{{x}^{{\prime} 2}+{y}^{{\prime} 2}}{\delta }^{2}\left(\displaystyle \frac{z^{\prime} }{\sqrt{{x}^{{\prime} 2}+{y}^{{\prime} 2}+{z}^{{\prime} 2}}}\right)\\ =\ \displaystyle \frac{{x}^{{\prime} 2}+{y}^{{\prime} 2}+{z}^{{\prime} 2}}{{x}^{{\prime} 2}+{y}^{{\prime} 2}}\displaystyle \frac{{z}^{{\prime} 2}}{{x}^{{\prime} 2}+{y}^{{\prime} 2}+{z}^{{\prime} 2}}\\ \quad \times \ \left[\displaystyle \frac{{\delta }^{2}(z^{\prime} )}{{z}^{{\prime} 2}}+\displaystyle \frac{{\delta }^{2}(\sqrt{{x}^{{\prime} 2}+{y}^{{\prime} 2}+{z}^{{\prime} 2}})}{{x}^{{\prime} 2}+{y}^{{\prime} 2}+{z}^{{\prime} 2}}\right]\\ =\ \displaystyle \frac{{z}^{{\prime} 2}}{{x}^{{\prime} 2}\,+\,{y}^{{\prime} 2}}\left\{\displaystyle \frac{{\delta }^{2}(z^{\prime} )}{{z}^{{\prime} 2}}+\displaystyle \frac{1}{{x}^{{\prime} 2}\,+\,{y}^{{\prime} 2}\,+\,{z}^{{\prime} 2}}\right.\\ \quad \times \,\left.\left[\displaystyle \frac{1}{4({x}^{{\prime} 2}+{y}^{{\prime} 2}+{z}^{{\prime} 2})}{\delta }^{2}({x}^{{\prime} 2}+{y}^{{\prime} 2}+{z}^{{\prime} 2})\right]\right\}\\ =\ \displaystyle \frac{{z}^{{\prime} 2}}{{x}^{{\prime} 2}\,+\,{y}^{{\prime} 2}}\left[\displaystyle \frac{{\delta }^{2}(z^{\prime} )}{{z}^{{\prime} 2}}+\displaystyle \frac{{x}^{{\prime} 2}{\delta }^{2}(x^{\prime} )+{y}^{{\prime} 2}{\delta }^{2}(y^{\prime} )+{z}^{{\prime} 2}{\delta }^{2}(z^{\prime} )}{{\left({x}^{{\prime} 2}+{y}^{{\prime} 2}+{z}^{{\prime} 2}\right)}^{2}}\right].\end{array}\end{eqnarray} \tag{ 18 }$$

Substituting the value and uncertainty of z from Equations (14) and (15) into Equations (16) and (18), we can obtain the value and uncertainty for the scattering angle, θ_scatter ± δθ_scatter.

In this paper, the input uncertainties are obtained from the GPI ${{ \mathcal Q }}_{\phi }$ image measured with the Debris Ring Analyzer package (Stark et al. 2014). The SPF is then the scattering-angle dependence of the flux density of the system at specific radial separations. The original measurements are then averaged to reduce measurement errors, and the final SPF is obtained by correcting the limb-brightening effect (by dividing the observed SPF by that of an isotropic disk model; Milli et al.2017).

In probability theory, for n independent and identically distributed real-valued random variables (${X}_{1},{X}_{2},\,\cdots ,\,{X}_{n}\in {\mathbb{R}}$), the empirical cumulative distribution function (ECDF) is defined as⁵⁸

$$\begin{eqnarray}&&{\hat{F}}_{n}(x)=\displaystyle \frac{1}{n}\displaystyle \sum _{i=1}^{n}{{\bf{1}}}_{[{X}_{i},\infty )}(x),\end{eqnarray} \tag{ 19 }$$

where ${{\bf{1}}}_{[{X}_{i},\infty )}(x)$ is an indicator function that is equal to 1 only when ${X}_{i}\leqslant x\lt \infty$ (otherwise ${{\bf{1}}}_{[{X}_{i},\infty )}(x)=0$). For the ECDF, the Glivenko–Cantelli theorem (e.g., Chung 2001) describes its asymptotic relationship with the cumulative distribution function (CDF) of a random variable $X\in {\mathbb{R}}$, which is denoted by F_X(x), that ${\hat{F}}_{n}(x)$ converges to F_X(x) uniformly to F_X(x) as $n\to \infty$ almost surely, i.e.,

$$\begin{eqnarray}&&| | {\hat{F}}_{n}-{F}_{X}| {| }_{\infty }\equiv \mathop{\sup }\limits_{x\in {\mathbb{R}}}| {\hat{F}}_{n}(x)-{F}_{X}(x)| \mathop{\to }\limits^{{\rm{a}}.{\rm{s}}.}0.\end{eqnarray} \tag{ 20 }$$

In statistics, for a random variable $X\in {\mathbb{R}}$ with CDF F_X(x), the PIT states that Y = F_X(x) is uniformly distributed between 0 and 1 (Rosenblatt 1952), in the sense that

$$\begin{eqnarray}\begin{array}{rcl}{F}_{Y}(y) & = & \Pr (Y\leqslant y)\\ & = & \Pr [{F}_{X}(x)\leqslant y]\\ & = & \Pr [X\leqslant {F}_{X}^{-1}(y)]\\ & = & {F}_{X}[{F}_{X}^{-1}(y)]\\ & = & y,\end{array}\end{eqnarray} \tag{ 21 }$$

which is the CDF of a random variable that is uniformly distributed between zero and 1, and F_Y(y) is the CDF of the random variable Y. Based on this property, the PIT is used to sample distributions, especially the ones that do not have parametric expressions.

To use the marginal posterior distribution from the previous MCMC run as the prior for the next run, we transfer the information between the two MCMC runs by combining the Glivenko–Cantelli theorem and the PIT.

First, convert discrete points to a continuous distribution. Based on the Glivenko–Cantelli theorem, for a specific random variable X with a large number of samples, we can treat its marginal ECDF from the previous MCMC run as its CDF, then use the ECDF as the prior for the next MCMC run.

Second, sample from an ECDF. We first sample a standard uniform random quantile variable $Y\in [0,1]$, then we find the corresponding empirical quantile in the given ECDF, i.e., ${\hat{F}}_{n}^{-1}(Y)$. Using the PIT in Equation (21), we have

$$\begin{eqnarray}&&{\hat{F}}_{n}^{-1}(Y)\sim {\hat{F}}_{n}.\end{eqnarray} \tag{ 22 }$$

Third, combining the Glivenko–Cantelli theorem in Equation (20) with the quantile distribution in Equation (22), we have

$$\begin{eqnarray}&&{\hat{F}}_{n}^{-1}(Y)\sim {F}_{X},\end{eqnarray} \tag{ 23 }$$

i.e., for a standard uniform random variable Y, its corresponding quantile for the ECDF of a random variable X follows the distribution of X.

In this section, we establish the connection between the GPI ${{ \mathcal Q }}_{\phi }$ image and the STIS total intensity image through radiative transfer modeling. We first obtain the posterior distribution of the disk parameters by radiative transfer modeling the GPI ${{ \mathcal Q }}_{\phi }$ image with MCMC (Section 4.2.1). We then calculate the marginal distributions from the MCMC posterior values⁵⁹ and use the PIT to treat them as the priors when modeling the STIS image.

The posteriors with the PIT approach for the STIS image are presented in Figure 18. In this section, our purpose is to demonstrate the statistical deviation of the posteriors from the GPI best-fit values. We thus constrain the PIT sampling ranges to be between the 2.5th and 97.5th percentiles, which is analogous to the conventional definition of two-tailed statistical significance: p-value >0.05.

Zoom In Zoom Out Reset image size

In this paper, we have ignored the correlated spatial noise in the images in our MCMC modeling; however, correlated noise is expected to increase the uncertainty of the extracted parameters (e.g., Czekala et al. 2015; Wolff et al. 2017). To better quantify the statistical deviation of the two sets of disk parameters that are extracted from the two disk images (Table 5), rigorous treatment of the correlated noise is necessary (e.g., Wolff et al. 2017).