An Exo-Kuiper Belt and An Extended Halo around HD 191089 in Scattered Light

TL;DR

This study presents detailed scattered light imaging of the HD 191089 debris disk, revealing a Kuiper Belt-like ring and an extended, bluer halo, with analysis of dust scattering properties and modeling challenges.

Contribution

It provides high-resolution images of the debris disk's structure and investigates dust scattering, highlighting the need for more complex models for accurate interpretation.

Findings

The debris disk has a Kuiper Belt-like ring at ~46 au and a halo extending to ~640 au.

The halo is bluer and exhibits more forward- and backward-scattering than the ring.

Radiative transfer modeling shows inconsistencies, suggesting more realistic dust models are required.

Abstract

We have obtained Hubble Space Telescope STIS and NICMOS, and Gemini/GPI scattered light images of the HD 191089 debris disk. We identify two spatial components: a ring resembling Kuiper Belt in radial extent (FWHM: ${\sim}$25 au, centered at ${\sim}$46 au), and a halo extending to ${\sim}$640 au. We find that the halo is significantly bluer than the ring, consistent with the scenario that the ring serves as the "birth ring" for the smaller dust in the halo. We measure the scattering phase functions in the 30{\deg}-150{\deg} scattering angle range and find the halo dust is both more forward- and backward-scattering than the ring dust. We measure a surface density power law index of -0.68${\pm}$0.04 for the halo, which indicates the slow-down of the radial outward motion of the dust. Using radiative transfer modeling, we attempt to simultaneously reproduce the (visible) total and (near-infrared) polarized intensity images of the birth ring. Our modeling leads to mutually inconsistent results, indicating that more complex models, such as the inclusion of more realistic aggregate particles, are needed.

Tables

Table 1. Table 1: System properties

Properties	HD 191089	Reference
Distance (pc)	$50.14\pm 0.11$	1
RA (J2000)	20 09 05.215	1
Dec (J2000)	-26 13 26.520	1
Spectral Type	F5V	2, 3
$M_{\star }$ ($M_{\odot }$)	$1.4\pm 0.1$	4
$T_{\mathrm{eff}}$ (K)	$6450$	1
$V$ (mag)	$7.18$	5
$J$ (mag)	$6.321$	6
$H$ (mag)	$6.091$	6
Association	$\beta$ Pic Moving Group	7
Age (Myr)	$22\pm 6$	8
$L_{\mathrm{dust}}/L_{\star }$	$\left(14.2\pm 0.5\right)\times {10}^{-4}$	9
Proper Motion (RA)	$40.17\pm 0.07$ mas yr-1	1
Proper Motion (Dec)	$-67.38\pm 0.05$ mas yr-1	1
Radial Velocity	$-5.4\pm 0.4$ km s-1	1

Table 2. Table 2: Observation Log

Instrument	Filter	${\lambda }_c$aaFor STIS, ${\lambda }_c$ is the pivot wavelength.	Pixel Scale	IWAbbIWA for STIS is the half-width of the wedge-shaped occulter.	$T_{\exp }$	$N_{\mathrm{frame}}$	$\Delta {\theta }_{\mathrm{PA}}$	UT Date
		($\mu$m)	(mas pixel-1)	(″)	(s)		(°)
NICMOS	F110W	1.12	75.65	0.3	4607.34	16	30.0	2006 May 27
STIS	50CCD	0.58	50.72	0.3	2048.00	64	109.8	2014 Jul 19, 2014 Aug 13
				0.5	5799.60	12
GPI	$H$-Pol	1.65	14.166	0.123	2484.72	28	102.2	2015 Sep 01

Table 3. Table 3: Reduction Parameters

	NICMOS	STIS	GPI
Classical RDI	NaaY: performed; N: not performed or unavailable.	YaaY: performed; N: not performed or unavailable.	N
Reference (Classical)	N	HD 196081	N
PSF References #	85	45 (A0.6)	100
		84 (A1.0)
KLIP Truncation #	19	N	20
NMF TruncationbbThe NMF reductions become stable with more than 10 components. #	10	10 (A0.6)	20
		10 (A1.0)
Polarimetry	N	N	Y

Table 4. Table 4: Disk Morphology Parameters

Parameter	Ring	Halo	Meaning
Instrument	GPI	STIS
${\theta }_{\mathrm{inc}}$aaFitted with Debris Ring Analyzer (Stark et al., 2014).	$59_{}^{\circ }{}_{-2^{\circ }}{}^{+4^{\circ }}$	${59}^{\circ }$ ccValues adopted from the GPI results.	Disk mid-plane inclination from face-on.
${\theta }_{\mathrm{PA}}$aaFitted with Debris Ring Analyzer (Stark et al., 2014).	$70_{}^{\circ }{}_{-3^{\circ }}{}^{+4^{\circ }}$	${70}^{\circ }$ ccValues adopted from the GPI results.	Position angle for major axis (N to E).
${\gamma }_{\mathrm{in}}$bbFitted for Equations (1) & (6).	$4.9\pm 0.2$	$\cdots$ddUnconstrained from the STIS image.	Surface brightness power law indices.
${\gamma }_{\mathrm{out}}$bbFitted for Equations (1) & (6).	$-6.1\pm 0.2$	$-2.68\pm 0.04$
$r_c$bbFitted for Equations (1) & (6).	$43.6\pm 0.3$ au	$\cdots$ddUnconstrained from the STIS image.	$r_c$ in Equation (1).
$r_{\mathrm{in}}$bbFitted for Equations (1) & (6).	$26\pm 4$ au	$\cdots$ddUnconstrained from the STIS image.	Inner and outer clearing radii.
$r_{\mathrm{out}}$bbFitted for Equations (1) & (6).	$78\pm 14$ au	$640\pm 130$ au
$r_{\mathrm{center}}$	$45.6\pm 0.2$ au	$\cdots$ddUnconstrained from the STIS image.	Peak position and FWHM of a Gaussian ring.
FWHM	$24.9\pm 0.4$ au	$\cdots$ddUnconstrained from the STIS image.

Table 5. Table 5: Independent Dust Properties Retrieved from Ring Image Modeling

ParameteraaThe morphological parameters in Table 4 are adopted.	PriorbbThe parameters are limited to 3 decimal digits, with uniform sampling in the prior range. For $M_{\mathrm{disk}}$ and $a_{\min }$, they are log-uniformly sampled.	Posteriorcc16th, 50th, and 84th percentiles.
Image		GPI ${𝒬}_{\varphi }$	STIS
${\log }_{10}{M}_{\mathrm{disk}}$ ($M_{\odot }$)	$(-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\left(\mathrm{Si}\right)$	(0, 1)	${0.005}_{-0.005}^{+0.016}$	${0.50}_{-0.05}^{+0.04}$
$f\left(C\right)$	(0, 1)	${0.004}_{-0.004}^{+0.011}$	${0.50}_{-0.04}^{+0.05}$
$f\left(\mathrm{ice}\right)$ddThe mass fraction for ice is $f\left(\mathrm{ice}\right)=1-f\left(\mathrm{Si}\right)-f\left(C\right)$, thus only Si and C are explicitly sampled.		${0.989}_{-0.010}^{+0.009}$	$0.003$ee$95$th percentile.
${\log }_{10}{a}_{\min }$ ($\mu$m)	$(-0.3,2)$	$-{0.115}_{-0.007}^{+0.013}$	${0.24}_{-0.03}^{+0.01}$
$f_{\max }$	(0, 1)	${0.213}_{-0.014}^{+0.026}$	$0.001$ee$95$th percentile.
${\chi }_{\nu }^2$ (GPI best fit)		1.05	240
${\chi }_{\nu }^2$ (STIS best fit)		12	10
$q$	3.5ffKept fixed.
${\log }_{10}{a}_{\max }$ ($\mu$m)	3ffKept fixed.
Iteration Number	$9000$ffKept fixed.
Burn-in	$2000$ffKept fixed.

Equations

$$F_{\nu }=\frac{{\lambda }_{\mathrm{c}}^2{F}_{\lambda }}{3\times 10^{-13}},$$

$$F_{\nu }=\frac{R{F}_{\star }}{S}=(6.87\pm 1.06)\times 10^{-7}\mathrm{Jy}\mathrm{s}\mathrm{coun}{\mathrm{t}}^{-1}.$$

$$\rho (r)\propto {\left\{{\left(\frac{r}{r_{\mathrm{c}}}\right)}^{-2{\alpha }_{\mathrm{in}}}+{\left(\frac{r}{r_{\mathrm{c}}}\right)}^{-2{\alpha }_{\mathrm{out}}}\right\}}^{-\frac{1}{2}},$$

$$Z(r,z)\propto \exp \left[-{\left(\frac{z}{\zeta (r)}\right)}^2\right],$$

$$\zeta (r)=h{r}^{\beta },$$

$$\beta =1$$

$$\Gamma =\alpha +\beta ,$$

$$\gamma =\Gamma -2,$$

$$(1\pm 0.5h\tan {\theta }_{\mathrm{inc}}{)}^{-1}.$$

$$\mathrm{d}N(a)\propto a^{-q}\mathrm{d}a.$$

$$\log \mathcal{L}\left(\theta \mid X_{\mathrm{obs}}\right)=$$

$$-\sum _{i=1}^N\log {\sigma }_{\mathrm{obs},i}-\frac{N}{2}\log (2\pi ).$$

$$a_{\mathrm{p}}=\frac{r_{\mathrm{in}}}{1+1.7(5M_{\mathrm{Jupiter}}/M_{\star }{)}^{0.31}}.$$

$$ax+by+cz=0,$$

$$\cos {\theta }_{\mathrm{inc}}=\frac{{\hat{n}}_{\mathrm{disk}}\cdot {\hat{n}}_{xOy}}{\sqrt{||{\hat{n}}_{\mathrm{disk}}|{|}^2||{\hat{n}}_{xOy}|{|}^2}},$$

$$\cos {\theta }_{\mathrm{inc}}=\frac{(a,b,c)\cdot (0,0,1)}{\sqrt{(a^2+b^2+c^2)(0^2+0^2+1^2)}}=c.$$

$${\theta }_{\mathrm{PA}}={\theta }_{\mathrm{slope}}-90^{\circ },$$

$$\tan {\theta }_{\mathrm{slope}}=\frac{\mathrm{d}y}{\mathrm{d}x}=-\frac{a}{b},$$

$$\tan {\theta }_{\mathrm{PA}}=\tan ({\theta }_{\mathrm{slope}}-90^{\circ })=-\cot {\theta }_{\mathrm{slope}}=\frac{b}{a}.$$

$$z$$

$$=-\sec {\theta }_{\mathrm{inc}}a\left(x+\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}}).$$

$${\delta }^{2}z=$$

$$=$$

$$+{\tan }^2{\theta }_{\mathrm{inc}}{\delta }^2(x\cos {\theta }_{\mathrm{PA}}+y\sin {\theta }_{\mathrm{PA}})$$

$$=$$

$$+{\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)]\}$$

$$=$$

$$+{\sec }^4{\theta }_{\mathrm{inc}}(x\cos {\theta }_{\mathrm{PA}}+y\sin {\theta }_{\mathrm{PA}}{)}^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}}).$$

$${\theta }_{\mathrm{scatter}}$$

$$=\arccos \left[\frac{z-z_0}{\sqrt{(x-x_0{)}^2+(y-y_0{)}^2+(z-z_0{)}^2}}\right].$$

$${\theta }_{\mathrm{scatter}}=\arccos \left\{\frac{\tan {\theta }_{\mathrm{inc}}(x\cos {\theta }_{\mathrm{PA}}+y\sin {\theta }_{\mathrm{PA}})-z_0}{\sqrt{(x-x_0{)}^2+(y-y_0{)}^2+{\left[\tan {\theta }_{\mathrm{inc}}(x\cos {\theta }_{\mathrm{PA}}+y\sin {\theta }_{\mathrm{PA}})-z_0\right]}^2}}\right\}.$$

$${\delta }^2{\theta }_{\mathrm{scatter}}=$$

$${\delta }^2\left[\frac{z-z_0}{\sqrt{(x-x_0{)}^2+(y-y_0{)}^2+(z-z_0{)}^2}}\right],$$

Full text

stars: imaging — stars: individual: HD 191089 — techniques: image processing — radiative transfer — protoplanetary disks

††facilities: HST (NICMOS, STIS), Gemini:South (Gemini Planet Imager)††software: corner.py (Foreman-Mackey, 2016), DebrisDiskFM (Ren & Perrin, 2018), Debris Ring Analyzer (Stark et al., 2014), emcee (Version 3.0rc1: Foreman-Mackey et al., 2013), MCFOST (Pinte et al., 2006, 2009), nmf_imaging (Ren, 2018), pysynphot (STScI Development Team, 2013).

\newdateformat

monthyearday\THEYEAR \monthname[\THEMONTH] \twodigit\THEDAY

\acceptjournallThe Astrophysical Journal

1 Introduction

Debris disks, the extrasolar analogs of the asteroid belt and Kuiper Belt, have been detected around ${\sim}20\%$ of the nearest stars (A-type stars: Thureau et al., 2014; FGK stars: Eiroa et al., 2013; Montesinos et al., 2016; Sibthorpe et al., 2018). They are expected to be the results from the grinding down of larger dust (Wyatt, 2008), however the diversity in observables such as morphology and surface brightness suggests that they are shaped by a variety of mechanisms (e.g., Artymowicz & Clampin, 1997; Stark et al., 2014; Lee & Chiang, 2016). Imaging studies of debris disks in scattered light use not only space-based instruments (e.g., STIS: Schneider et al., 2009, 2014, 2016, 2018; Konishi et al., 2016; NICMOS: Soummer et al., 2014; Choquet et al., 2016, 2017, 2018) that offer the best telescope stability, but also extreme adaptive optics–equipped ground-based instruments (e.g., GPI: Perrin et al., 2015; Kalas et al., 2015; Hung et al., 2015; Draper et al., 2016; Millar-Blanchaer et al., 2015, 2016; Esposito et al., 2018; SPHERE: Boccaletti et al., 2015; Lagrange et al., 2016; Wahhaj et al., 2016; Feldt et al., 2017; Milli et al., 2017; Matthews et al., 2017; Engler et al., 2017; Sissa et al., 2018; Olofsson et al., 2018; Milli et al., 2019) that provide the best angular resolution and probe closer-in regions of the disks.

Multi-wavelength studies can provide complementary insights in understanding circumstellar disks, since different wavelengths probe distinct regions and parameter space for a disk (Ertel et al., 2012; Sicilia-Aguilar et al., 2016). In this paper, we focus on using scattered light observations to understand one of these systems. In previous studies, the combination of both space- and ground-based instruments has been implemented to study both protoplanetary (e.g., PDS 66: Wolff et al., 2016) and debris disks (e.g., 49 Ceti: Choquet et al., 2017, HD 35841: Esposito et al., 2018), and those observations are interpreted using radiative transfer codes (e.g., Augereau & Beust, 2006; Milli et al., 2015; Wolff et al., 2017; Esposito et al., 2018). We perform such an analysis for the debris disk surrounding HD 191089 to study its specific properties via measurement and radiative transfer modeling in this paper.

We list the properties of the system in Table 1: HD 191089 is an F5V star with $T_{\text{eff}}=6450$ K located at $50.14\pm 0.11$ pc (Gaia DR2: Gaia Collaboration et al., 2018). Moór et al. (2006) identified it as a member of the Beta Pictoris moving group with space velocities compatible with the kinematics of the group, and Shkolnik et al. (2017) estimated the age of this group to be $22\pm 6$ Myr based on the consensus of the group members.

Before a resolved scattered light image was reported, HD 191089 was first identified by Mannings & Barlow (1998) as a debris disk candidate based on IRAS infrared excess. Chen et al. (2014) suggested a two-temperature model to explain the Spectral Energy Distribution (SED) of the system, while Kennedy & Wyatt (2014) argued for one temperature. The latest SED analysis including CSO, Herschel, and JCMT photometry up to 850 $\mu$m seems to confirm the latter hypothesis with a best fit obtained using a single-component disk: assuming the dust behaves as a blackbody, SED analysis suggests a dust mass ${\sim}0.037M_{\earth}$, a temperature of $89$ K, and a radius ${\sim}17$ au (Holland et al., 2017). However, as is commonly seen for many disks, the radius derived assuming blackbody dust is several times smaller than the radius observed in resolved images (e.g., SEDs: Mittal et al., 2015; images: Hughes et al., 2018, and references therein). This indicates that the dust grains are not simple black bodies, and the smallest dust are not efficient emitters at long wavelengths.

The HD 191089 disk was first resolved using Gemini/T-ReCS at $18.3\ \micron$ by Churcher et al. (2011), with the region interior to $28$ au reported to have little emission. The disk was then detected in scattered light in a re-analysis of the archival 2006 HST/NICMOS observations by Soummer et al. (2014), with the apparent disk extent and orientation consistent with Churcher et al. (2011).

To further characterize the debris disk, we observed the target with HST/STIS and Gemini/GPI. By carrying out a multi-wavelength study, we aim to understand: (1) the spatial distribution of the dust; (2) the scattered light color of the dust; (3) the scattering phase functions of the dust for the different spatial components of the system (if available), and how they suit the trends of the current observed debris disks; (4) the dust properties, including size distribution, structure, and compositional information; and (5) whether a universal description of the dust is able to explain the observations across different wavelengths and observational techniques.

The structure of this paper is as follows. In Section 2, we describe our HD 191089 observations and the data reduction procedure. In Section 3, we describe measurables derived from the observations. In Section 4, we describe our radiative transfer modeling efforts in studying the disk. In Sections 5 and 6, we discuss our findings and provide concluding remarks.

2 Observations & Data Reduction

2.1 Data Observation and Reduction

2.1.1 HST/NICMOS (2006)

HD 191089 was observed using HST/NICMOS (Proposal ID: 10527, PI: D. Hines) with the NIC2-CORON aperture and F110W filter ($\lambda_{\text{c}}=1.12\ \micron$, inner working angle [IWA]: $0\farcs 3$, pixel scale: $75.65$ mas pixel*-1*, Viana et al., 2009) on 2006 May 27, with two telescope orientations each observing the target for $2303.67$ s. These HST/NICMOS data, totaling 16 frames, have previously been 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 star light and speckle noise and reveal the debris disk around HD 191089.

We obtain the scattered light image of the system using the Multi Reference 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)111https://archive.stsci.edu/prepds/alice/. For each target exposure, we first decompose its 10% closest ALICE images in correlation (i.e., $L^{2}$-norm sense , and a total of 85 images) using both the Karhunen-Loève Image Projection (KLIP, Soummer et al., 2012) and the Non-negative Matrix Factorization (NMF222 Python-based nmf_imaging package: Ren (2018)., Ren et al., 2018a) 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 element-wise mean of the derotated individual residual exposures333The ${\sim}0.9\%$ $x$- and $y$-direction scale difference in Schneider et al. (2003) is ignored..

To calibrate the NICMOS disk image and obtain a measurement of its surface brightness, we multiply the reduced data by the calibrated F110W PHOTFNU parameter444http://www.stsci.edu/hst/nicmos/performance/photometry/postncs_keywords.html $F_{\nu}=1.21\times 10^{-6}$ Jy s count*-1*, then divide the data by the NICMOS pixel area on-sky to obtain the surface brightness data in units of Jy arcsec*-2*.

2.1.2 HST/STIS (2014)

HD 191089 was observed using HST/STIS (Program ID: 13381, PI: M. Perrin) using the 50CORON aperture ($\lambda_{\text{c}}=0.58\ \micron$, pixel scale: $50.72$ mas pixel*-1*, Riley et al., 2018) on 2014 July 19 and 2014 August 13 (2 visits each), totaling 76 frames. Each visit was performed with a different telescope orientation, with position angles of $-84\fdg 23,-61\fdg 23,-2\fdg 23$, and $25\fdg 59$ 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 $0\farcs 3$. Then we obtained 3 longer 483.3 s exposures on the WEDGEA1.0 position (half-width: $0\farcs 5$) 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 191089555http://www.stsci.edu/hst/phase2-public/13381.pro) 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 star light 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 element-wise 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_{\lambda}=4.15\times 10^{-19}$ erg cm*-2* Å*-1* count*-1*) to $F_{\nu}=4.56\times 10^{-7}$ Jy s count*-1* using the conversion equation in Appendix B.2.1 of Viana et al. (2009):

[TABLE]

where $\lambda_{\rm c}=0.58\ \mu$m is the pivot wavelength of STIS. We then multiply our combined STIS image in count s*-1* pixel*-1* by $F_{\nu}$, and divide it by the STIS pixel area on-sky to obtain the surface brightness data in units of Jy arcsec*-2*.

2.1.3 Gemini/GPI (2015)

HD 191089 was observed using Gemini/GPI in $H$-band ($\lambda_{\text{c}}=1.65\ \micron$, pixel scale: $14.166\pm 0.007$ mas pixel*-1*, De Rosa et al., 2015) polarimetric mode (“$H$-Pol”, Perrin et al., 2015) on 2015 September 01 during the Gemini Planet Imager Exoplanet Survey (GPIES; PI: B. Macintosh, Macintosh et al., 2014). We took 28 exposures, each with $88.74$ s integration time with a total field rotation of $102\fdg 2$. The airmass ranged from $1.002$ to $1.01$, the differential motion image monitoring (DIMM) seeing measurement was $1\farcs 16\pm 0\farcs 17$, and the multi-aperture scintillation sensor (MASS) seeing was $1\farcs 1\pm 0\farcs 3$.

To obtain the Stokes cube ({$I,Q,U,V$}) for the HD 191089 debris disk, we follow 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 ${\pm}45^{\circ}$ from it (Monnier et al., 2019). On the local Stokes maps of HD 101089, we notice two similar low spatial frequency octopole structures with a rotation of ${\sim}45^{\circ}$. 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 is to be one of the instrumental artifacts (Esposito et al., 2019). 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^{\circ}$ and remove it from the $\mathcal{Q}_{\phi}$ map.

We flux calibrate these data following the procedure described in Hung et al. (2015, 2016): we first adopt for HD 191089 an $H$-band flux of $F_{\star}=3.749\pm 0.119\,\rm{Jy}$ (Skrutskie et al., 2006), then adopt the satellite-to-star ratio $R=2.035\times 10^{-4}$ from 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}\,\rm{count\,s^{-1}}$ in the FITS file header. Combining these, we obtain a conversion factor of

[TABLE]

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*-2*. The flux-calibrated data are then geometrically corrected and smoothed by convolving with a Gaussian kernel ($\sigma=14.166$ mas, i.e., the scale of one GPI pixel: Millar-Blanchaer et al., 2016) to remove the high spatial frequency noise that impacts regions smaller than the Nyquist-sampled point-source PSF of GPI.

2.2 Noise Estimation

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

NMF & 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 pixel-wise 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_{\rm 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.

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., ${\sim}50\%$ of the entire library). We split these $413$ images into 25 groups each containing 16 images (with the left-overs randomly discarded). In each group, the $16$ images are treated as the mock target images to simulate 16 HD 191089 non-detection images. For one mock target image, we first identify the real target that was observed, then 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 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 element-wise mean of the 16 rotated reduced mock target images as the mock result for one realization of HD 191089 non-detection. We then obtain a total of $25$ realizations of non-detections using the $25$ full groups from the $413$ mock images. We take the element-wise standard deviation from the $25$ mock non-detections 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 SNR calculation.

2.3 Comparison of the Reduction Methods

The HD 191089 debris disk is detected in all of our HST and GPI observations. We here 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 the NMF reductions (Figure 4). The system shows a bright parent belt, surrounded by a faint and diffuse halo detected up to ${\sim}6\arcsec$ from the star. These features are detected with a dynamic range of 3 orders of magnitude, with the peak ring surface brightness of ${\sim}1$ mJy arcsec*-2* and the diffuse halo detected down to a few $\mu$Jy arcsec*-2*. Theoretically, unlike KLIP reduction, both cRDI and NMF reduction methods are free from over-subtraction caused by over-fitting disk features using the references. This was discussed in Ren et al. (2018a) and is confirmed in Ren et al. (2018b) where the NMF method is able to successfully retrieve the spiral arms for the MWC 758 system.

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

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 star light, 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.

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 non-negative — 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 non-negative combination of non-negative 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 ${\sim}2$ times the ring surface brightness of the NICMOS-KLIP image, supporting the expected behavior of the two methods: KLIP’s over-fitting because of direct vector projection cannot be avoided even with a small number of eigen-modes, 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 dataset 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 SNR 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 upperlimit 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\sigma$ uncertainty of ${\sim}3$ mJy arcsec*-2* at the minor axis, and ${\sim}1$ mJy arcsec*-2* 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 brightness of ${\sim}2.2$ mJy arcsec*-2* and ${\sim}1.0$ mJy arcsec*-2*, respectively. Therefore, even the disk is not removed by the reduction methods, it is below a detection threshold of $1\sigma$ (which corresponds to a total SNR of 5 for extended structure spanning ${\sim}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.

3 Disk Morphology & Measurement

3.1 Strategy to Measure Disk and Dust Properties

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, as well as the 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 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 they are coplanar, then we can measure the disk surface brightness as a function of scattering angle—the scattering phase function (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.

3.2 Mathematical Description

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

[TABLE]

where $\alpha_{\rm in}>0$ and $\alpha_{\rm out}<0$ approximate the mid-plane dust density power law indices interior and exterior to $r=r_{\rm c}$ (Augereau et al., 1999). Along the normal direction of the disk mid-plane, the dust follows a Gaussian dispersion form, i.e.,

[TABLE]

with

[TABLE]

where

[TABLE]

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

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

[TABLE]

where $\alpha$ is $\alpha_{\rm in}$ or $\alpha_{\rm out}$ in Equation (1) , with $\Gamma$ being $\Gamma_{\rm in}$ or $\Gamma_{\rm 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, $\gamma$, and the surface density power law index, $\Gamma$, is

[TABLE]

where $\gamma$ is $\gamma_{\rm in}$ or $\gamma_{\rm out}$, corresponding with $\Gamma$ being $\Gamma_{\rm in}$ or $\Gamma_{\rm out}$.

The radial distribution of the system in Equation (1) is supplemented with two extra parameters $r_{\rm in}$ and $r_{\rm 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_{\rm in}$ or $r>r_{\rm out}$, $\rho(r)=0$.

3.3 Disk Morphology

3.3.1 Ellipse Parameters

We measure the geometric parameters for the disk from the GPI $\mathcal{Q}_{\phi}$ observations, because it has the smallest pixels, and is less biased by the reduction methods.

We assume that the peak radial polarized surface density matches the peak radial particle density, and use the Debris Ring Analyzer package from Stark et al. (2014) to fit the peak intensity of the ring in $10^{\circ}$ azimuthal wedges by minimizing the $\chi^{2}$ value between the observation and an ellipse model, and quantify the uncertainties by assessing the change of the $\chi^{2}$ values for the corresponding degrees of freedom on a grid of the explored parameters (Choquet et al., 2018).

We measure an inclination666The uncertainties calculated in this paper are $1\sigma$ unless otherwise specified. of $\theta_{\rm inc}={59^{\circ}}_{-2^{\circ}}^{+4^{\circ}}$ from face-on, and a semimajor axis of $45_{-1}^{+2}$ au , and a position angle of $\theta_{\rm PA}={70^{\circ}}_{-3^{\circ}}^{+4^{\circ}}$ from North to East for the major axis. There is no significant offset between the location of the star and the center of the ring ($3\sigma$ upper limit: 8 au). The $3\sigma$ upper limit for the eccentricity of the deprojected ellipse is $0.3$.

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 SNR 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 the 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_{\rm in}$ and $r_{\rm out}$, as the radii for which the corresponding average surface brightness is consistent with zero at the $1\sigma$ level.

We measure the following parameters for the ring from the GPI $\mathcal{Q}_{\phi}$ observation: $\gamma_{\rm in}=4.9\pm 0.2$, $r_{\rm c}=43.6\pm 0.3$ au, and $\gamma_{\rm out}=-6.1\pm 0.2$.777For the input data, we have taken both $x$- and $y$-uncertainty into account with the orthogonal distance regression method (Boggs et al., 1989). We also fit a Gaussian ring to the deprojected $\mathcal{Q}_{\phi}$ surface density map, and the ring is centered at $r_{\rm center}=45.6\pm 0.2$ au with a $24.9\pm 0.4$ au full width at half maximum (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 $\gamma_{\rm out}=-2.68\pm 0.04$. The power law indices measured at different position angles are also consistent with the integrated measurement within $1\sigma$. Therefore, we report the integrated radial profiles to reduce systematic uncertainty.

In the above calculation, we have assumed a flat disk (i.e., $h\ll 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 non-flared disk, we calculate that the radial separation will be modified by multiplicative factors of

[TABLE]

For the HD 191089 system, these line-of-sight intersections increase the uncertainties by ${\sim}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 mid-plane. We thus ignore this effect given its model-dependent minor impact on the uncertainties.

We find that the radial power law indices, $\gamma_{\rm out}$, for the ring and the halo differ by ${>}10\sigma$. 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.

3.4 Scattering Phase Function (SPF)

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 A, then measure the empirical SPFs for the ring and the halo for our HD 191089 observations.

For the ring and the halo, we present Figure 7 the SPFs averaged for both sides (i.e., the NE and SW sides) using STIS-cRDI and NICMOS-NMF results to minimize over-subtraction. 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.

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\ \micron$ (Section 2.3), we instead use the $1.12\ \micron$ NICMOS-NMF observation that is both close to the GPI wavelengths and has less over-subtraction 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 ${\sim}110^{\circ}$, however it cannot be as firmly established as in previous measurements (e.g., Perrin et al., 2014; Milli et al., 2019; Frattin et al., 2019). To better constrain the polarized fraction values, we need $H$-band total intensity observations to rule out the wavelength-dependent effect.

To compare SPFs in the STIS image, we first normalize the SPFs by dividing their average surface brightness at $90^{\circ}\pm 10^{\circ}$ 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 scattering and backward scattering in the probed scattering angles.

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.

3.5 Disk Color

Using the HST disk images, we calculate the color of the dust as follows: we first bin the $1.12\,\mu$m NICMOS-NMF image and the $0.58\,\mu$m STIS-NMF image to two images with pixel size of ${\sim}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.

For the ring, the dust scatters ${\sim}25\%$ more light ($\Delta$mag $\approx 0.25$) within the NICMOS F110W passband than that in STIS 50CCD, showing a red scattering property. For the halo, $\Delta$mag $\approx-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 outwards 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.

4 Disk Modeling

4.1 Radiative Transfer Modeling Tool

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_{\rm max}$)—parameterizes the central vacuum fraction in the dust that is uniformly ranging from 0 to $f_{\rm max}$. DHS has been used to successfully reproduce the scattering phase function of linearly polarized scattered light with incident unpolarized light for quartz particles in 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, maximum void fraction.). DHS is a computationally intensive technique due to the complex nature of the constituent dust, we thus adopt parallel computation in the Python environment and use the DebrisDiskFM package (Ren & Perrin, 2018)888https://github.com/seawander/DebrisDiskFM. 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 cluster999We performed the calculations on the Maryland Advanced Research Computing Center (MARCC): https://www.marcc.jhu.edu., 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 the input parameters. We also append the folder name by a hashed random number, preventing multiple nodes simultaneously accessing the same folder and causing errors.

4.2 Modeling the Ring

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”, Draine & Lee, 1984, denoted by “Si”), amorphous carbonaceous dust (Rouleau & Martin, 1991, denoted by “C”), and ${\rm H_{2}O}$-dominated ice described in Li & Greenberg (1998) to model the $\beta$ Pic disk (denoted by “ice”). The size of the dust, $a$, follows a power-law distribution with index $q$, i.e.,

[TABLE]

The distribution is truncated at a minimum size of $a_{\rm min}$, and we set a maximum limit of the dust size, $a_{\rm max}=1000$ $\mu$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:

• •

Disk mass, $M_{\rm disk}$, which generally controls the overall brightness of the disk at different wavelengths;

• •

Porosity;

• •

Mass fraction for “astronomical silicates”, $f({\rm Si})$;

• •

Mass fraction for amorphous carbonaceous dust, $f({\rm C})$;

• •

Mass fraction for ice, $f({\rm ice})$;

• •

Minimum dust size, $a_{\rm min}$;

• •

Maximum void fraction for DHS, $f_{\rm max}$.

Given the fact that the composition parameters are interconnected, i.e., $f({\rm Si})+f({\rm C})+f({\rm ice})=1$, we only explicitly sample $f({\rm Si})$ and $f({\rm C})$. We also set the lower limit for $a_{\rm min}$ to be $0.5\,\mu$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)101010http://www.stsci.edu/software/tinytim/ for STIS, and with a 2-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 compare the PSF-convolved $\mathcal{Q}_{\phi}$ model with the observation. With the models and observations, we maximize the following log likelihood function:

[TABLE]

In the above equation, $X$ is a flattened image with $N$ pixels. Subscripts obs and model denote observation and model, respectively. $\sigma_{{\rm obs},i}$ is the uncertainty for $X_{{\rm obs},i}$ at the $i$-th pixel. We only focus on the disk region (i.e., between $r_{\rm in}$ and $r_{\rm out}$ determined by GPI $\mathcal{Q}_{\phi}$ image) to minimize the influence from the halo. In this paper, we assume 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 burn in stage, and calculate the posterior distributions using the last $7000$ steps (a total of $4.2\times 10^{5}$ models).

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, 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

One-Step STIS

We fit the STIS image using the same priors as for the GPI image. We run MCMC modeling for $9000$ steps with $60$ chains, and discard the first $2000$ as burn-in steps as for the GPI $\mathcal{Q}_{\phi}$ image. For this approach, we also present the posterior results in Table 5 and Figure 13.

Two-Step MCMC

Noticing the two distinct sets of MCMC posteriors for the two images, we try to establish the connection between the GPI image and STIS image with a two-step MCMC fitting: we obtain the posterior distributions for the seven variables from the GPI $\mathcal{Q}_{\phi}$ image, then use them as the priors to fit the STIS image. In this way, we use the GPI posterior distribution as the null hypothesis, and test it on the STIS image.

We use Probability Integral Transform (PIT: B) to draw samples from the GPI posterior distributions. To find the posterior ranges from the STIS image, we only explore the GPI posteriors between their 2.5th and 97.5th percentiles as an analogy to the conventional level of significance for one parameter (i.e., $p$-value $\geq 0.05$ for double-tailed distribution). We run MCMC modeling for $3500$ steps with $60$ chains, and discard the first $1000$ as burn in stage.

For the STIS modeling with the PIT approach, we generate the posterior distributions and overplot them on the one-step GPI data in Figure 18 (Appendix B.4). From the following, we notice the statistical deviation of the STIS posteriors from the GPI posteriors: (1) the STIS posteriors using PIT are adjacent to the 2.5th- or 97.5th-percentile boundaries of the GPI posteriors, indicating the trend of drifting away from the null hypothesis; and (2) the STIS posteriors using PIT have extremely narrow credible intervals, indicating there is no statistically preferred solution within the explored intervals.

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 over-predicts the emission from the system, while the SED from the GPI $\mathcal{Q}_{\phi}$ best fit under-predicts the emission. Even though different sizes of dust dominate the SED and the scattered light images, the inability to reproduce the observed SED adds another 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.

Mie vs DHS

Although we do not obtain a set of parameters that is able to explain both datasets, we can still focus on the $f_{\rm 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_{\rm 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_{\rm max}\approx 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_{\rm max}\approx 0.21$. In this scenario, the DHS theory is preferred to the Mie theory, and the $f_{\rm max}$ parameter is smaller than what is in the interstellar dust (e.g., 0.7: Min et al., 2007), suggesting dust properties are different under various environments. Given that the two datasets differ in two fundamental way (different wavelengths, and total intensity/polarization observation), we cannot determine which is the root cause of this discrepancy. However, this discrepancy in the $f_{\rm 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:

• •

Physically Possible: minimum dust size. The retrieved minimum dust size is ${\sim}0.8\,\mu$m for GPI $\mathcal{Q}_{\phi}$and ${\sim}1.7\,\mu$m for STIS. Both values are within a factor of ${\sim}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 blow out 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).

• •

Physically Impossible: porosity and composition fraction. Specifically, these parameters takes values in limited ranges (i.e., from $0\%$ to $100\%$), and the discrepancy is currently at ${\sim}50\%$ level. Even though the discrepancies for these parameters can be alleviated with larger uncertainties, their physical meaning are then uninformative—a possible solution is to increase the uncertainties to ${\sim}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 loosen prior constraints.

For GPI

Set the lower limit in the prior for $a_{\rm min}$ to be 0.01 $\mu$m, and keep $q$ unconstrained. We observed a steeper $q\approx 4.15$ with $a_{\rm min}\approx 0.02\,\mu$m to describe the GPI $\mathcal{Q}_{\phi}$ image, but it still does not recover the STIS flux density. In addition, although these 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 & 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 datasets 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.

5 Discussion

5.1 Spatial Distribution

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_{\rm in}$ and $r_{\rm out}$): Churcher et al. (2011) observed the ring at $18.3\micron$, reporting a dust belt from 26 to 84 au111111Updated with the Gaia DR2 distance to HD 191089., which is consistent with our fitting results of $r_{\rm in}=26\pm 4$ au and $r_{\rm out}=78\pm 14$ au at $1\sigma$ 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 ${\sim}10\%$ or $1\sigma$ with the GPI and STIS data. Although a tentative ${\sim}20\%$ asymmetry was observed at the $1.8\sigma$ level in the $18.3\mu$m observation by Churcher et al. (2011), if the dust follows the same spatial distribution at these wavelengths, the $18.3\ \micron$ emission asymmetry is likely resulted from 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\sigma$ point-source contrast limits of ${\sim}1\times 10^{-5}$–${\sim}3\times 10^{-6}$ between $0\farcs 3$ and $0\farcs 8$ with the forward modeling planet detection method to correct for self- and over-subtraction 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 ${\sim}5~{}M_{\rm Jupiter}$.

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 semi-major axis of the planet’s orbit:

[TABLE]

When we substitute the measured $r_{\rm in}$ into the above equation, we obtain a lower limit of $a_{\rm p}=20\pm 3$ au for the semi-major axis of the planet.

5.1.2 Ring as an Exo–Kuiper Belt

The spatial extent of the ring ($r_{\rm center}=45.6\pm 0.2$ au, FWHM $=24.9\pm 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 SNR 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 with the strongest resonance orbits of a hypothetical potential well.

We set the boundary ranges to be consistent with pixel-wise ${\rm SNR}\approx 1$: among different combinations of the strongest resonances, we obtain 4 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 either too narrow to cover the GPI disk, or 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 4 resonance radii pairs, we are able to compute the location of the potential well (i.e., $1:1$ resonance) at $r_{(1:1)}=0\farcs 60\pm 0\farcs 02$, correponding with a stellocentric distance of $r_{(1:1)}=29.9\pm 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 SNR 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_{\rm out}=640\pm 130$ au, with a surface density power-law index of $\Gamma_{\rm out}=-0.68\pm 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 slow-down of the radial motion of the dust corresponds a power law index of $-0.5$, and a constant acceleration results into 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 slow-down from gravity and the acceleration from radiation pressure; the $-0.68_{-0.04}^{+0.04}$ power law index is thus calling for additional slow-down sources, and the slow-down by 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 has 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)\propto r^{x}$, where $r$ is the radial separation, then the surface density will be $\Gamma(r)\propto r^{-x-1}$. In this way, the dust in the halo of HD 191089 has a radial speed of $v(r)\propto r^{-0.32\pm 0.04}$.

Local Distribution

The surface density power law of the HD 191089 halo deviates from the classical model at ${>}10\sigma$ level, calling for detailed investigation for the local variation of the halo. As an attempt to investigate the variation of the surface density distribution at different stellocentric distance, 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 that

• •

Interior to ${\sim}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.

• •

Between ${\sim}200$ au and ${\sim}300$ au, the radial speed reaches a constant then increases as stellocentric distance increases. At ${\sim}300$ au, the surface density power law index reaches that for the classical model for debris disk halo (e.g., Strubbe & Chiang, 2006; Thébault & Wu, 2008).

• •

Exterior to ${\sim}300$ au, albeit with large uncertainty, the radial speed marginally increases then reaches a constant as stellocentric distance increases.

Given the complex 2-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 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 2-dimensional surface density variation for the HD 191089 halo.

To investigate the deviation of the scattering properties from a 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 and 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 $\gamma_{\rm out}=-2.68$ correction is applied.

The 2-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 side of the STIS data at ${\sim}25\%$ level, with the STIS NE region likely matching the NICMOS NE overdensity region. The under-density region to the NW region in the STIS image is likely influenced by the truncation of signal by STIS’s Wedge B (Figure 2).

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 with local surface density variations; (2) when there is no density variation, then 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 $11\farcs 4$ to the southwest of HD 191089 in the STIS field of view (partially seen at the bottom right corner in Figure 1a) is a background star. It is identified by Gaia DR2 with Source ID 6847146784384527872 at $d=1.06\pm 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 (${\sim}0.5$ Myr ago). We notice three stars that have the nearest projected approach from ${\sim}1.2$ pc to ${\sim}1.5$ pc, which happened between ${\sim}0.2$ Myr and ${\sim}0.35$ Myr from now121212The results do not change using the $329$ stars within $10$ pc.. For star encounter events, the closest approach are typically smaller than ${\sim}200$ au (e.g., Pfalzner, 2003; Pfalzner et al., 2018), and distinct features such as spiral arms dissipate beyond ${\sim}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.

5.2 Dust Properties from Radiative Transfer Modeling

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 encompass 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 maximum void fraction of ${\sim}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 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 of correlated noise (e.g., Wolff et al., 2017). However, this resolution does not change the best-fit values, and in this way it would yield 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 wavelength and 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; $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 scattering phase functions 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 as large 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 co-planar with the ring, then the halo will bias the SPF of the ring and add another dimension of complexity to the problem.

6 Summary

In this paper, we report our detection and characterization of the HD 191089 debris disk by combining space- and ground-based instruments: HST/NICMOS, HST/STIS, and Gemini/GPI. Using these three instruments, we are able to study the disk in scattered light at three different wavelengths: $0.58\,\micron$ and $1.12\,\micron$ in total intensity, and $1.65\,\micron$ ($H$-band) in polarized intensity. In the scattered light images, we are able to identify two components in the debris disk system: a ring, and a fainter fan-like halo structure. For the STIS and GPI $\mathcal{Q}_{\phi}$ images, we implement radiative transfer modeling to retrieve the dust information. Assuming the ring and the halo are coplanar, we summarize our findings as follows.

Measurement:

• •

The HD 191089 system has two spatial components: one exo–Kuiper Belt ring from $26\pm 4$ au to $78\pm 14$ au, and an halo extending to $640\pm 130$ au. The center of ring does not have a significant offset from that of the star. The halo has an overall radial surface density power law index of $-0.68\pm 0.04$, with local variations indicative of the interaction with the interstellar medium.

• •

The ring has an inclination of ${59^{\circ}}_{-2^{\circ}}^{+4^{\circ}}$, enabling the scattering phase function measurements for the two components from ${\sim}30^{\circ}$ to ${\sim}150^{\circ}$ in the STIS data. In the range of scattering angles probed by our observations, both forward and backward scattering are stronger for the dust in the halo than that in the ring.

• •

The polarization fraction curve calculated using $1.65\,\mu$m GPI $\mathcal{Q}_{\phi}$ and $1.12\,\mu$m NICMOS-NMF images does not have a clear trend, however it is possibly peaking at ${\sim}110^{\circ}$. Our result may be influenced by wavelength, and it can be better constrained with future total intensity observations in $H$-band.

• •

From the color of the dust derived from the NICMOS and STIS observations, the dust in the ring is redder than that in the halo. This is consistent with larger dust in the ring, which is also consistent with theoretical simulations that the ring serves as the “birth ring” for the smaller dust in the halo.

• •

In comparison with an identical SPF trend (Hughes et al., 2018), the SPFs of the dust in the HD 191089 system are likely deviating from the trend.

• •

If the ring is shaped by the strongest orbital resonances, the gravitational well is likely at $29.9\pm 1.2$ au, resulting a $4:3$ resonance with the inner edge, and a $5:2$ resonance with the outer edge.

Radiative Transfer Modeling:

We use DHS theory to model the HD 191089 ring images observed with STIS and GPI.

• •

Most of the extracted parameters are statistically not consistent with each other (e.g., composition, porosity). Specifically, the GPI $\mathcal{Q}_{\phi}$ ring favors DHS theory with ${\sim}21\%$ maximum void fraction, while the STIS ring favors Mie theory (i.e., DHS with $0\%$ void fraction). This maximum void fraction parameter approximates the structure of the dust, which is the only parameter that is added between the two theories. The discrepancy for this parameter thus suggests that neither DHS nor Mie is a good approximation of the dust structure. However, both values are smaller than the best fit for interstellar dust (0.7: Min et al., 2007), suggesting that dust properties are different in different enviroments.

• •

The discrepant dust parameters retrieved in our DHS radiative transfer modeling of the ring may be mathematically resolved with larger uncertainties. However, for the parameters that takes limited range of values from $0\%$ to $100\%$—e.g., porosity and mass fractions of compositions—large uncertainties will render these parameters less informative on dust properties. Advanced description of dust models such as aggregates are expected to physically resolve the discrepancies, however such MCMC analyses are currently limited by computational power.

We appreciate the suggestions from the anonymous referee, which significantly improved this paper. B.R. thanks the useful discussions with Xinyu Lu and the comments from Kevin Schlaufman which improved the paper. E.C. acknowledges support from NASA through Hubble Fellowship grant HST-HF2-51355 awarded by STScI, operated by AURA, Inc. under contract NAS5-26555, and support from HST-AR-12652, for research carried out at the Jet Propulsion Laboratory, California Institute of Technology. T.E. was supported in part by NASA Grants NNX15AD95G/NEXSS, NNX15AC89G, and NSF AST-1518332. C.P. acknowledges funding from the Australian Research Council via FT170100040 and DP180104235. G.D. acknowledges support from NSF grants NNX15AD95G/NEXSS, AST-1413718 and AST-1616479. This research has made use of data reprocessed as part of the ALICE program, which was supported by NASA through grants HST-AR-12652 (PI: R. Soummer), HST-GO-11136 (PI: D. Golimowski), HST-GO-13855 (PI: É. Choquet), HST-GO-13331 (PI: L. Pueyo), and STScI Director’s Discretionary Research funds, and was conducted at STScI which is operated by AURA under NASA contract NAS5-26555. The input images to ALICE processing are from the recalibrated NICMOS data products produced by the Legacy Archive project, “A Legacy Archive PSF Library And Circumstellar Environments (LAPLACE) Investigation,” (HST-AR-11279, PI: G. Schneider). Based on observations obtained at the Gemini Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), National Research Council (Canada), CONICYT (Chile), Ministerio de Ciencia, Tecnología e Innovación Productiva (Argentina), Ministério da Ciência, Tecnologia e Inovação (Brazil), and Korea Astronomy and Space Science Institute (Republic of Korea). This research project (or part of this research project) was conducted using computational resources (and/or scientific computing services) at the Maryland Advanced Research Computing Center (MARCC).

Appendix A The Scattering Phase Function

A.1 3D Cartesian Coordinates of the Disk System

To quantify the coordinate values, the $x$ and $y$ coordinates can be directly measured in the image of the system with $(x,y)\pm(\delta x,\delta y)$, where $\delta$ 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 semi-major 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 towards West, $\hat{y}=(0,1,0)$ towards North, $\hat{z}=(0,0,1)$ pointing towards the observer. Let $\hat{n}_{\rm disk}=(a,b,c)$ denote the unit normal vector of the mid-plane of the debris disk system. Then, all the points on the disk mid-plane, which also contains the origin $O$, satisfies

[TABLE]

with $a^{2}+b^{2}+c^{2}=1$.131313Note: the $a$ symbol in this section is not the dust size used in the main text.

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

[TABLE]

where $\cdot$ 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

[TABLE]

The position angle of the system, which is denoted by $\theta_{\rm PA}\pm\delta\theta_{\rm PA}\in[0^{\circ},180^{\circ}]$, is defined as the angle from North 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 $\theta_{\rm slope}$, which is defined as the counter-clockwise angle from $\hat{x}$ to the line. Then, we have the relationship between the mathematical slope angle and astronomical position angle,

[TABLE]

where

[TABLE]

therefore we have,

[TABLE]

For the points on the disk mid-plane, we can substitute Equations (A2) and (A3) into Equation (A1), then we have the $z$-coordinate of the points as

[TABLE]

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

[TABLE]

Now combining the $(x,y)$ coordinates and the inclination and position angle of the system, we can obtain $z\pm\delta z$ from Equation (A4) and the square-root of Equation (A5). In this paper, for the HD 191089 disk, the $1\sigma$ uncertainties for $x$ and $y$ are estimated to be $0.33$ pixel141414Scaled from a conservative $3\sigma$ uncertainty of 1 pixel, and the uncertainties for $\theta_{\rm inc}$ and $\theta_{\rm PA}$ are obtained from the GPI $\mathcal{Q}_{\phi}$ image using the Debris Ring Analyzer package by Stark et al. (2014).

A.2 Scattering Angle

From the $(x,y,z)$ coordinates of the dust in the debris disk system, and given the position of the star at $(x_{0},y_{0},z_{0})\pm(\delta x_{0},\delta y_{0},\delta z_{0})$, 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 $\vec{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 and the angle between $\vec{r}$ and $\hat{z}$, is thus

[TABLE]

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

[TABLE]

Assuming the measured parameters are independent, the corresponding squared uncertainty for $\theta_{\rm scatter}$ is then

[TABLE]

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 {widetext}

[TABLE]

Substituting the value and uncertainty of $z$ from Equations (A4) and (A5) into Equations (A6) and (A8), then we can obtain the value and uncertainty for the sacttering angle, $\theta_{\rm scatter}\pm\delta\theta_{\rm 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).

Appendix B Sampling from Posterior Distributions

B.1 Glivenko-Cantelli Theorem

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 as151515Note: the $x$ and $y$ symbols in B are statistical variables, not Cartesian coordinates.

[TABLE]

where $\mathbf{1}_{[X_{i},\infty)}(x)$ is an indicator function which is equal to 1 only when $X_{i}\leq x<\infty$ (otherwise $\mathbf{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\rightarrow\infty$ almost surely, i.e.,

[TABLE]

B.2 Probability Integral Transform (PIT)

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

[TABLE]

which is the CDF of a random variable that is uniformly distributed between 0 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.

B.3 Posterior as Prior

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 find the corresponding empirical quantile in the given ECDF, i.e., $\hat{F}_{n}^{-1}(Y)$. Using the PIT in Equation (B3), we have

[TABLE]

Third, combining the Glivenko-Cantelli Theorem in Equation (B2) with the quantile distribution in Equation (B4), we have

[TABLE]

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$.

B.4 GPI Posteriors as STIS Priors

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 values161616Note: the correlation of the parameters are ignored., 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.

In this paper, we have ignored the correlated spatial noises 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).