Direct Imaging of the HD 35841 Debris Disk: A Polarized Dust Ring from Gemini Planet Imager and an Outer Halo from HST/STIS

We observed HD 35841 with GPI in H-band (λ_cen = 1.647 μm) using its spectroscopic ("H-Spec") and polarimetric ("H-Pol") modes as part of the Gemini Planet Imager Exoplanet Survey (GPIES; PI B. Macintosh). Details of the data sets are listed in Table 1. In both cases, the pixel scale was 14.166 ± 0.007 mas lenslet⁻¹ (De Rosa et al. 2015), a 123 mas radius focal plane mask (FPM) occulted the star, and angular differential imaging (ADI; Marois et al. 2006) was employed (the default for GPI). The average atmospheric properties for the H-Spec data set were DIMM seeing = 10 ± 02, MASS seeing = 05 ± 01, coherence time τ = 5.4 ± 1.2 ms, and airmass ranging from 1.01 to 1.06. For H-Pol, the airmass ranged from 1.08 to 1.19, but atmospheric measurements were not available from the observatory.

The spectroscopic observations divide the filter bandpass into micro-spectra that are measured by the detector and then converted into 37 spectral channels per image by the GPI Data Reduction Pipeline (DRP; Perrin et al. 2014, 2016). We used this pipeline's standard methods to assemble the raw data into 44 spectral data cubes (similar to steps taken in De Rosa et al. 2016). The star location in each channel was determined from measurements of the four fiducial "satellite" spots (Sivaramakrishnan & Oppenheimer 2006; Wang et al. 2014; Pueyo et al. 2015), as was the photometric calibration, assuming a satellite-spot-to-star flux ratio of 2.035 × 10⁻⁴ (Maire et al. 2014) and a stellar H magnitude of 7.842 ± 0.034 (Two Micron All Sky Survey [2MASS]; Cutri et al. 2003). We did not high-pass filter or smooth any of the GPI images used for our measurements and analysis. In this paper we consider only the broadband-collapsed results from the spectral data; any consideration of the disk's spectrum in reflected light is left for a future work.

We applied multiple techniques to subtract the stellar point-spread function (PSF). First, we used pyKLIP, a Python implementation (Wang et al. 2015a) of the Karhunen–Loève Image Projection (KLIP) algorithm (Soummer et al. 2012; Pueyo 2016). Subtraction was performed on a channel-by-channel basis using only angular diversity of reference images (no spectral diversity was used), and the aggressiveness of the PSF subtraction was adjusted by varying the KLIP parameters. We show aggressive and conservative reductions in Figure 1. The aggressive reduction divided each image radially into eight equal-width annuli between r = 10 and 85 pixels (px) with no azimuthal division of the annuli, used a minimum rotation threshold of 1 px to select reference images, and projected onto the first 44 KL modes. The conservative reduction was identical, except it employed only the first KL mode. The PSF-subtracted images were derotated so north is up and averaged into the final image.

Zoom In Zoom Out Reset image size

When using the aggressive pyKLIP parameters with the entire 44-image data set, we found the PSF subtraction to preferentially self-subtract the ring along its most southeast edge (Milli et al. 2012). No such effect was found for the conservative reduction. The effect possibly arose because more images were taken after transit than before transit, leading to an unequal distribution of the disk's position angle (PA) among reference images; we could avoid the bias with the aggressive parameters by using a subset of 30 images that was balanced in number of images pre- and post-transit. However, this resulted in lower S/N for the rest of the ring, so we chose to present the 44-image version in Figure 1 to better display the ring's other features and illustrate this phenomenon.

Separately, we used a modified version of the "locally optimized combination of images" (LOCI) algorithm (Lafrenière et al. 2007) on images that were median-collapsed across spectral channels. This collapse allowed faster forward modeling of the disk self-subtraction (Esposito et al. 2014) during the modeling discussed in Section 4 but did reduce S/N compared to non-collapsed reductions. The reduction presented herein used only one subtraction annulus at r = 9–120 px divided azimuthally into three subsections, with LOCI parameter values of N_δ = 0.5, W = 4 px, dr = 120 px, g = 0.625, and N_a = 160, following the conventional definitions in Lafrenière et al. (2007). To prevent speckle noise at the edge of the FPM from detrimentally biasing subtraction over the entire annulus, we set the inner radius of the region used to optimize the LOCI coefficients to 12 px instead of 9 px. We found the preferential self-subtraction of the southeast edge noted above to also occur here with a 44-image data set, so we used the more PA-balanced subset of 30 images instead. Finally, the PSF-subtracted frames were rotated to place north up and collapsed into the final median image shown in Figure 1.

The polarimetric data were created with GPI's Wollaston prism, which splits the light from the spectrograph's lenslets into two orthogonal polarization states. To combine the raw data into a set of polarization data cubes, we used the GPI DRP with the methods outlined in Esposito et al. (2016) and described in more detail in Perrin et al. (2015) and Millar-Blanchaer et al. (2015). Specific to this data set, the instrumental polarization was removed by first estimating the apparent stellar polarization in each data cube as the mean normalized difference of pixels 2–11 px from the star's location (i.e., both inside and just outside of the FPM). For each pixel in the cube, we then scaled that value by the pixel's total intensity and subtracted the product. The data cubes were also corrected for geometric distortion and smoothed with a Gaussian kernal (σ = 2 px). Combining the data cubes results in a three-dimensional Stokes cube containing the Stokes parameters {I, Q, U, V}, rotated to place north up. Finally, the Stokes cube was photometrically calibrated using the satellite spot fluxes by assuming the same satellite-spot-to-star flux ratio as we did for H-Spec and following the methods described in Hung et al. (2015). We recovered only a very low S/N total intensity detection from this H-Pol data set with 38 of field rotation, so we use only the H-Spec data for total intensity analysis.

We observed HD 35841 on 2014 November 6 with the STIS instrument on HST in its coronagraphic mode (program GO-13381, PI M. Perrin). The system was observed at two telescope roll orientations of −789 and −949 over two orbits. These orientations were chosen to align the disk's major axis, estimated by Soummer et al. (2014), perpendicular to the occulting wedge for one of the rolls, but it was ultimately offset by 24° due to scheduling constraints. At each orientation, we acquired a series of six 120 s exposures with the star centered on the 06 wide WEDGEA0.6 wedge position (hereafter A0.6), and then three longer 485 s exposures on the 1'' wide WEDGEA1.0 position (hereafter A1.0). This resulted in a total exposure time of 1440 s for separations of 03–05 but 4,350 s for separations >05.

The STIS coronagraphic mode is unfiltered ("50CCD") and sensitive from ∼0.20 to 1.03 μm, with a pivot wavelength of λ_p = 0.5754 μm (Riley 2017). The pixel scale is 50.77 mas pixel⁻¹ (Schneider et al. 2016).

To subtract the stellar PSF from the science images, we also observed HD 37002 as a reference star at a single orientation during the single orbit between visits to HD 35841. This minimized potential PSF differences caused by telescope thermal breathing. HD 37002 is an F5V star chosen for its close spectral match, similar luminosity, and on-sky proximity to HD 35841. We acquired six 110 s exposures on A0.6 and three 505 s exposures on A1.0.

We processed the A0.6 and A1.0 data sets separately using classical reference star differential imaging with the following steps. After flat-fielding via the STIS calstis pipeline and correction of the bad pixels, we registered and scaled the images of the target and reference star by minimizing the quadratic difference between each of them and the first reference image. The star center is estimated from cross-correlation of the secondary mirror struts diffraction spikes, with the absolute star center determined from a Radon transform of the first reference image (Pueyo et al. 2015). For each science frame, we subtracted either the closest reference frame or the median of all reference frames, choosing the version that minimized PSF residuals. Finally, we consolidated the PSF-subtracted frames for both wedges into one pool, rotated them to set north up, masked the areas covered by the wedges and diffraction spikes, and average-combined all of the frames using a pixel-weighted combination of their respective exposure times.

After combination, we found that some stellar background remained that was approximately azimuthally symmetric. To remove this background, we fit a sixth-order polynomial to the median radial profile measured within PA = 30°–140° (avoiding the disk), and subtracted that polynomial function from the image at all radii. We use the resulting image for all analysis, apart from one instance in Section 5.2.

We converted the final image to surface brightness units of mJy arcsec⁻² via the average "PHOTFLAM" header value of 4.116 × 10⁻¹⁹ erg cm⁻² s⁻¹ Å⁻¹ and the pixel scale. For comparison with the GPI images, we linearly interpolated the final STIS image to match the GPI plate scale, and this is the version shown in Figure 1.

We also reduced a Keck/NIRC2 ADI H-band data set from 2014 February 08 but failed to detect the disk with statistical significance. The data set comprised 97 frames of 20.0 s integrations with the 400-mas diameter coronagraph in place, totaling 139 of parallactic rotation. Unfortunately, data quality was suboptimal due to high humidity (∼70%) and variable seeing of 06–08. The compact angular extent of this disk, combined with observing conditions and a total integration time of only 32 minutes, meant the signal was not recovered from the residual speckle noise, despite attempts with classical ADI (using a median-collapsed reference PSF), LOCI, and KLIP subtractions (see Appendix A for details).

We present our total intensity GPI (spectroscopic mode) and STIS images in Figure 1. The GPI images represent the aggressive KLIP, conservative KLIP, and LOCI reductions. They show a highly inclined ring of dust with sharp inner and outer edges. We assume the ring to be approximately circular, as both ansae extend to the same projected separation (∼065 with brightness ∼3σ above the local background noise) and there is no obvious stellocentric offset along the minor axis. We consider the ansae to be the portions of the ring near its intersection with the projected major axis (i.e., the inflection point of the ellipse).

The strong brightness asymmetry in the aggressive KLIP image is a reduction artifact, as the higher KL modes preferentially self-subtract the southeast edge due to the imbalance of reference image PA's previously discussed. Therefore, we use the conservative KLIP and LOCI images as the bases for our measurements and analyses. Nevertheless, we present the aggressive KLIP image because it provides the best view of the NW back edge, which is swamped by speckled residuals in conservative PSF subtractions. Additionally, bright spots along the major axis on both sides of the star at the inner working angle are likely speckle residuals, rather than point sources or ansae of an inner ring.

Photometrically, the west edge of the ring (PA 166°–346°, measured east of north) appears consistently brighter than the east edge. From here on, we consider this west edge to be the "front edge" between the star and the observer, and the east edge to be the "back edge" behind the star. We base this on assumptions that the dust grains are primarily forward scattering, their scattering properties are constant around the ring, and the ring is approximately azimuthally symmetric. While the back edge is intrinsically fainter due to a forward-scattering phase function, it is also artificially dimmed by self-subtraction by the brighter front edge. We correct for this bias in our measurements and modeling but not in the images shown in Figure 1.

The outer edges of the ansae extend symmetrically to projected separations of r_proj ≈ 67 au (065). We detect the ring down to our inner working angle of r_proj ≈ 12 au (012) along the front edge, but the back edge is only marginally detected above the speckle noise level at r_proj ≈ 27 au (026). Interior to ∼19 au (018), the residual speckle noise is substantial and reduces the significance of our detection. The ring appears generally smooth, without clumps or vertical protrusions.

The STIS image is heavily impacted by the combined orientations of the occulting wedges in the constituent frames. Consequently, we detect just the ring ansae (and only partially in the NW). However, we also detect an asymmetric, low surface brightness component that extends at least 07 outward from the ring's outer edge and is preferentially seen west of the star. This smooth halo or "dust apron" becomes fainter with separation and reaches the background limit at r_proj ≈ 140 au (14). It is the likely source of the wing-tilt asymmetry in the Soummer et al. (2014) NICMOS 1.1 μm data, with forward scattering grains leading to preferential brightening of the disk's west side and creating an apparent deflection of the ansae toward that direction. This halo is reminiscent of similar features seen, for example, in the Fomalhaut, AU Mic, HD 32297, and HD 15745 disks (Kalas et al. 2005; Strubbe & Chiang 2006; Chiang et al. 2009; Schneider et al. 2014), but it is not sharply deflected from the main ring like the HD 61005 disk (Schneider et al. 2014).

We present the polarized intensity GPI data in Figure 2 as the radial components and of the Stokes and parameters, respectively (Schmid et al. 2006; Millar-Blanchaer et al. 2015). The signal shares roughly the same extent and shape of the total intensity ring, with greater brightness along the front edge than the back edge. The inner hole is not prominent in , which suggests it may be enhanced by ADI self-subtraction in the total intensity images; this is supported by the modeling we show later (see Figure 6). The southeast (SE) side of the disk appears brighter than the northwest (NW) side, particularly in the region where we assume the back edge to be. The image contains no clear disk signal but shows a quadrupole pattern that may result from instrumental polarization unsubtracted during reduction. We use this image to estimate noise in the data because we do not expect single scattering by circumstellar grains to generate a significant signal (Canovas et al. 2015). Dividing the image by a noise map, built from the standard deviation of pixel values in 1 px wide annuli centered on the star, we created the signal-to-noise ratio (S/N) image shown in the right-hand panel of Figure 2.

Zoom In Zoom Out Reset image size

The ring's edges appear softer in than in total intensity, particularly on the SE side. This is likely because only the total intensity ring is biased by ADI self-subtraction, which partially resembles a high-pass filter. Soft edges might indicate a vertically extended ring. In this case, light scattered by the front edge may be conflated with light scattered by the back edge for scattering angles >90°. This would affect both polarized and total intensity measurements. We discuss this possibility further in Section 4.4.

We quantitatively assessed the disk's scattering phase function by measuring its surface brightness as a function of scattering angle θ (Figure 3). These angles assume a circular ring centered on the star with an inclination of 849 (determined from modeling described in Section 4). To measure the scattering phase function, we placed apertures (radius = 2 px) on the conservative KLIP ring at a range of scattering angles from 22°, located closest to the star on the front edge, to 154°, closest to the star on the back edge (see Figure 3 inset). The ansae are at θ ≈ 90°. The NW and SE sides of the ring were measured independently.

Zoom In Zoom Out Reset image size

To estimate the self-subtraction of disk brightness by KLIP PSF subtraction, we forward modeled the effect with the "DiskFM" feature included in pyKLIP. This projects the relevant principal components onto a model of the disk's underlying brightness distribution and considers the effects of disk signal leaking into the principal components, accounting for both over- and self-subtraction of the disk. For the underlying brightness, we used the "median model" result from our MCMC fit, described in Section 4.3. We then computed the ratio of that underlying brightness to the corresponding forward-modeled brightness at every pixel to get a correction factor for each aperture, and finally multiplied each aperture measurement by said factor to get the corrected brightness values plotted in Figure 3. All of our total intensity brightness and flux measurements include these corrections.

The error bars represent 1-σ uncertainties. For each data point, they are the quadrature sum of (1) the standard deviation of mean surface brightnesses within apertures placed at the same separation but outside the disk, and (2) an assumed 15% error on the self-subtraction correction factor based on variances of measurements for different models and reductions. For measurements that are consistent with zero brightness according to our uncertainties, we plot 3-σ upper limits only.

We find the ring's brightness to be symmetric with scattering angle between its SE and NW halves. The one exception is our innermost measurement at θ ≈ 22°, for which our errors may be underestimated due to non-Gaussian noise from residual speckles close to the star. The brightness along the ring's front edge decreases by a factor of ∼20 from θ ≈ 22° to the ansae. The ring brightness along the back edge (θ > 90°) is roughly consistent with being constant in θ, although it is largely unconstrained at θ > 125°. This general shape is consistent with several other debris disks with measured phase functions (Hughes et al. 2018).

We performed similar brightness measurements on the image, also plotted in Figure 3. The uncertainties are calculated from the image as the standard deviation of mean surface brightnesses within apertures placed at the same separations as the data measurements. This assumes the noise properties are similar between and . No self-subtraction corrections are needed for our polarized intensities.

There is less variation of the polarized intensity with θ than for total intensity, as the front edge is only about 1.5–2.0 times brighter than the ansae. The back edge brightness again has relatively large uncertainties; however, it may be slightly fainter than the ansae. The SE side of the ring is preferentially brighter than the NE side, particularly on the front edge, but the asymmetry is marginal given our photometric precision. Nevertheless, these phase functions provide useful points of comparison for our models, particularly regarding grain properties.

Having brightness measurements for both total and polarized intensity, we computed their ratio to get a polarization fraction for the ring, plotted in Figure 4. The 1-σ uncertainties were propagated from the uncertainties on both sets of brightness measurements, and we exclude measurements for which we have only upper limits on the total intensity or Q_r brightness.

Zoom In Zoom Out Reset image size

The polarization fraction is generally higher to the SE than the NW but not to a significant degree, given our uncertainties. The fraction peaks near the ansae at ∼25%–30% but may be as low as a few percent at the smallest scattering angles. The location of the peak near θ = 90° (our measurement is at θ = 87°) is consistent with most predictions for scattering by micron-sized grains. Large uncertainties on brightness measurements along the back edge make for poorly constrained polarization fractions at large scattering angles; however, we see a tentative trend of the fraction decreasing as the angle increases past 90°. We do not report the fraction for θ > 130° because it is unconstrained between 0% and 100% within the 3σ uncertainties on our total intensity and brightness measurements.

The HD 35841 polarization fractions are similar to those measured for other debris disks. For example, AU Mic (Graham et al. 2007) and HD 111520 (Draper et al. 2016) both peak at 40% at the largest separations, which may be the ansae of those edge-on rings.

To compute a STIS−GPI H color, we first measured the flux within one 3 × 3 px aperture centered on each ansa in the GPI conservative KLIP image and the interpolated STIS image. We make this measurement at the ansae because they are the only places that we detect the disk with both GPI and STIS. The same aperture centers were used for both images, with locations relative to the star of (r, PA) = (058,1658) and (058, −142) for the SE and NW ansae, respectively. Both correspond to a scattering angle of 87°. The NW aperture lies just outside of the STIS image's masked region. We then subtracted a stellar STIS−GPI H color of 1.10 mag from the difference of the fluxes. The stellar color is based on the 2MASS H magnitude and an 8.88 ± 0.01 mag in the STIS 50CCD bandpass (converted from the Tycho 2 V-band value of 8.90 ± 0.01 mag; Høg et al. 2000).

We measure the STIS−GPI H color to be mag and mag for the ring's SE and NW ansae, respectively. This makes the ring slightly blue on both sides, along the lines of the optical versus near-IR colors of debris disks like AU Mic (Lomax et al. 2018 and references therein) and HD 15115 (which is blue in V − H according to Kalas et al. 2007 and Debes et al. 2008). With only two measurements, we limit our speculation as to the physical interpretation of the disk color and look forward to future visible-light observations that resolve more of the ring.

Our observations yielded no significant point sources. Based on our data, we determined limits on our sensitivity to substellar companions around HD 35841. For this, we only consider H-Spec because it achieved deeper contrast than H-Pol and lower mass limits than STIS.

We based our contrast measurements on reductions optimized for point-source detection and separate from those already presented. In this case, we used pyKLIP on the full 44-image data set, taking advantage of both angular and spectral diversity (i.e., ADI + SDI). Images were divided into nine equal-width annuli between r = 10 and 115 px that were split azimuthally into four subsections. References were restricted by a minimum rotation threshold of 1 px. We selected the first 30 KL modes among the 300 most correlated references for each target image. (More references are available now because we do not spectrally collapse the data.) The PSF was subtracted assuming two different spectral templates for a hypothetical planet: one with a flat spectrum and one with a methane-absorption spectrum (e.g., similar to that of 51 Eri b; Macintosh et al. 2015). To correct for point-source attenuation by the KLIP algorithm, we injected fake Gaussian point sources into the input images and then recovered their fluxes after PSF subtraction. The fake planet spectrum was matched to the reduction type, as either flat or methane-absorbing.

Our 5σ equivalent point-source contrast limits (Mawet et al. 2014; Wang et al. 2015b), corrected for PSF subtraction throughput, are shown in Figure 5. We translated these contrast values into planet mass limits using AMES-Cond atmosphere models (Allard et al. 2001; Baraffe et al. 2003) to convert planet luminosity to mass assuming an age of 40 Myr and "hot start" formation. At moderate separations of 08–13 (82–134 au), we can rule out planets more massive than 5 M_J and can more generally exclude 12 M_J companions or larger at 02–14 (21–144 au). These limits only apply to regions outside of the ring, however, as a planet embedded in or interior to the ring could be obscured by the dust. Therefore we are most sensitive to companions that are not coplanar with the disk.

Zoom In Zoom Out Reset image size

Our underlying disk model distributes grains in an azimuthally symmetric ring that is centered on the star. It consists of one component, which Soummer et al. (2014) found sufficient to fit the disk's SED. We chose MCFOST's "debris disk" option to define the disk's volume density profile, which follows the form

where r is the radial coordinate in the equatorial plane, z is the height above the disk midplane, and R_c is a critical radius that divides the ring into inner and outer regions with density power law indices of α_in and α_out, respectively (Augereau et al. 1999). The disk scale height varies as H(r) = H₀(r/R₀) such that the scale height is H₀ at radius R₀, while the slope of the vertical density profile is constrained by the exponent γ. We chose to set R₀ = 60 au because that radius appears by eye to pass through the middle of the ansae (i.e., it is the midpoint between the ring's inner and outer edges).

The ring's inner and outer edges are set by R_in and R_out, respectively, and tapered by a Gaussian with σ = 2 au so the volume density smoothly declines to zero (separate from Equation (1)). We found the outer radius to be weakly constrained in preliminary small-scale MCMC tests because the ring's edge gradually blends into the GPI data's background noise at ∼80 au. However, we know from the STIS image that the disk extends out to at least 110 au. Therefore, we set a conservative outer ring radius of twice this distance, R_out = 220 au, and hold it constant throughout the fitting process. The viewing geometry of the ring is set by the disk inclination i and position angle PA.

We populate the single disk component with grain particles following the power law size distribution dN(a) ∝ a^−q da, where the grain size a ranges from a minimum size a_min to a maximum size of 1 mm. We consider grains composed of three different materials³⁶ : astrosilicates (Si), amorphous carbon (aC), and water ice (H₂O). The relative abundances of these compositions are defined as fractions of the total disk mass (, , ) and are allowed to vary so long as their mass fractions sum to 1. However, all grains have pure compositions (e.g., no individual grain is 50% Si, 50% aC), are distributed in the same manner throughout the disk volume regardless of composition, and share the same size distribution and porosity within a given model. MCFOST approximates porous particles by "mixing" a grain's material composition with void (refractive index of n = 1+0i); the mixture is performed using the so-called Bruggeman mixing rule of effective medium theory to get the effective refractive index of the grains. The total dust mass M_d regulates the model's scattered-light surface brightness and thermal flux.

We do not include radiation pressure, Poynting–Robinson (P–R) drag, or gas drag effects in our model. Being relatively bright, we expect this disk's dust density to be high enough that P–R drag is negligible (Wyatt 2005). Only a non-detection of gas has been reported for this disk (Moór et al. 2011a), so we assume a standard debris disk scenario in which most of the protoplanetary gas has been removed and gas drag is also negligible. We exclude radiation pressure for simplicity, but it may have important effects on the outer disk, which we discuss later.

We derived the disk's primary morphological and grain parameters by fitting scattered-light models to our GPI total intensity and images via MCMC. Only the LOCI image of the total intensity was used in the MCMC because it allowed for faster forward modeling of the ADI self-subtraction than the KLIP image. The STIS image was not included in the fit due to its limited coverage of the ring but was used for comparison afterwards.

Uncertainty maps, calculated as the standard deviation among pixels in the data at the same projected radius from the star, were constructed at the original GPI resolution for both total intensity and . We then binned the images and uncertainty maps into 2 × 2 px bins to mitigate correlation between pixels within the same resolution element (∼3.5 px across) in our data. Ideally to remove correlations, the bin size should be equal to the size of one resolution element, but we found that binning the data to this degree removed many of the disk's defining morphological features. Alternatively, the correlations between different elements can be measured and explicitly taken into account as part of the fitting process (Czekala et al. 2015; Wolff et al. 2017).

For the H-band models, we scattered only photons with a wavelength of 1.647 μm, approximate to the central wavelength of the GPI H bandpass. We found use of a single wavelength to be a computationally efficient proxy for integrating models over the entire bandpass. To compare models to data in each iteration of the MCMC, we first constructed the models at the same resolution as the original GPI data. We then convolved them with a normalized 2D Gaussian function with a 3.8 px full-width half-maximum to approximate the GPI PSF. For total intensity only, we then forward modeled the LOCI self-subtraction using the "raw" total intensity model output by MCFOST as the underlying brightness distribution (Esposito et al. 2014, 2016). It is this forward-modeled version that we compare with the LOCI image, providing a fair comparison to the self-subtracted data. No forward modeling was necessary for the models.

Our parallel-tempered MCMC used three temperatures with 150 walkers each. We initialized walkers randomly from a uniform distribution, and then simulated each walker for 11,000 iterations (4.95 × 10⁶ samples). Parameter values were constrained by a flat prior with the limits listed in Table 2. Ultimately, the walker chains evolved significantly for ∼10,000 iterations before stabilizing (i.e., "converging"); therefore our posterior distributions are drawn from the final 1000 iterations (1.5 × 10⁵ samples) of the zeroth temperature walkers only.

Table 2.  MCMC Model Parameters

Param.	Limits	Max Lk	16%	50%	84%
i (deg)	[80, 88]	85.1	84.7	84.9	85.1
PA (deg)	[163, 167]	165.9	165.6	165.8	165.9
Rin (au)	(10, 65]	59.9	57.7	59.8	60.9
H0 (au)	(0.3, 10]	2.4	2.4	2.7	4.1
γ	(0.1, 3]	0.52	0.51	0.56	0.68
Md (M⊕)	(0.00053, 3.3)	0.11	0.11	0.15	0.19
amin (μm)	[0.1, 40]	1.5	0.16	1.5	1.6
q	(2, 6)	2.9	2.7	2.9	3.0
	(0, 1)	0.63	0.35	0.48	0.61
	(0, 1)	0.31	0.21	0.27	0.33
αout	(−6, 0)	−2.9	−3.2	−3.0	−2.8
Parameters considered upper or lower limitsa
αin	(−6, 7)	4.0	−1.6	3.8	6.2
Rc (au)	(40, 110]	51.0	45.8	53.2	57.0
porosity (%)	(0, 95)	1.4	0.50	1.5	3.3
	(0, 1)	0.058	0.079	0.24	0.42

Note.

^aThe bottom four parameters all have posterior PDFs that extend to either an upper or lower boundary imposed by our MCMC prior. Therefore, these should be formally considered lower (α_in) or upper limits (R_c, porosity, ).

Download table as:  ASCIITypeset image

The results of the MCMC are listed in Table 2 as the parameter values associated with the maximum likelihood model (i.e., "best fit") and also the 16th, 50th (i.e., median), and 84th percentiles of the marginalized posterior probability distribution functions (PDF). A corner plot for those distributions is provided in Appendix B.

We find the ring to be inclined and au wide, with dust between radii of au and 220 au. Vertically, the scale height is au at R₀, with a density profile exponent less than unity (γ = 0.51–0.68). Both parameters are weakly bimodal, favoring vertically thin disks (H₀ ≈ 2.4 au, γ ≈ 0.52) but also showing thicker disks (H₀ ≈ 4.1 au, γ ≈ 0.68) to agree nearly as well with the data. The marginalized PDF for R_c abuts our lower prior boundary of 40 au, so we only place an upper limit of 57.0 au (its 84th percentile value) on it. However, R_c < R_in suggests a sharp inner edge to the ring. This also makes α_in degenerate in most cases, so we can only place a lower limit for it at −1.6 (its 16th percentile value). The outer volume density power law α_out is better defined at . Therefore, the dust density decreases continuously from a peak near R_in to the outer edge. The PA is .

The dust properties are less constrained than the ring's morphological properties. We find the most likely minimum grain size a_min to be 1.5 μm, but there is a weaker secondary peak in the marginalized posterior at ∼0.16 μm. The blowout size (a_blow) for this star³⁷ is ∼1.6–2.1 μm; grains smaller than a_blow experience a radiation pressure force greater than the star's gravitational force and are thus ejected from the system. This additional constraint leads us to accept the larger a_min of ∼1.5 μm as the most likely value. The power law index of the grain size distribution is , and the dust mass in grains with sizes between a_min and 1 mm is approximately 0.1–0.2 M_⊕. The median mass fractions among the three compositions are 24% , 48% , and 27% , although extends to the lower prior boundary of 0%. We also note that the maximum likelihood model (which is thinner vertically) favors carbon more strongly, with mass fractions of 6% , 63% , and 31% . Grains with low porosity are strongly preferred overall, with a 99.7% confidence upper limit of <12%. We discuss some of the implications of these dust properties in Sections 5.4 and 5.5.

We present "median model" images alongside the data and data-model residuals in Figure 6. This median model is constructed using the median value of each parameter's PDF from the MCMC. As Table 2 shows, the median values are generally close to those of the maximum likelihood model, making the two models nearly indistinguishable to the eye and essentially equal in terms of (to two decimal places for GPI data and only differing by 0.01 for the STIS comparison).

Zoom In Zoom Out Reset image size

The model is a good match to the data, with little residual disk signal. A quadrupole pattern is apparent in the residual map, which is a sign of instrumental polarization that was not completely removed during data reduction. The model's signal is at least 100 times below the noise floor of our data and consistent with no disk signal. The forward-modeled total intensity agrees well with the data along the ring's front edge, but the model is comparatively faint at the ansae and the back edge. Our model grains, therefore, are more forward scattering than the real grains. The underlying total intensity model (far left panel) contains a much more vertically extended scattered-light distribution than the forward-modeled version, demonstrating how the ADI data reduction sharpens the ring's edges and generally suppresses its surface brightness.

The same underlying model recomputed for 0.575 μm scattered light (bottom of Figure 6) agrees well with the STIS image out to projected separations of ∼110 au but fails to account for halo brightness at larger separations. This is evidenced by positive residuals NW and SW of the star (the positive residuals to the NE are localized noise and not disk signal). The discrepancy may be even more pronounced near the minimum of θ where our data are incomplete. We know our model contains dust at these large separations but more scattering is required to match the observed brightness.

Thus, we propose that this outer halo brightness is produced by an additional population of grains slightly larger than a_blow that are produced in the ring and then pushed onto eccentric orbits with large apoapses by radiation pressure (Strubbe & Chiang 2006). We do not include radiation pressure directly in our models, so we do not expect them to contain scattered light from such a population. For a consistency check, however, we can approximate this eccentric dust with a separate, manually tuned model. Doing so, we find that a rudimentary model of a broad annulus at 220–450 au containing 1.7 × 10⁻³ M_⊕ of grains with a = 1.5–3.0 μm, the same inclination and composition as the ring, and H₀ = r/10 appears qualitatively similar to the outer halo in the STIS image. Its H-band brightness is also below the GPI data's sensitivity limits (for the data reductions in this work), consistent with it going undetected with GPI.

We do not explicitly fit to the measured phase functions in Figures 3 and 4; however, we do so implicitly when fitting the scattered-light images. Therefore, we can extract phase functions from our median model images, using the previously described aperture photometry method, and compare them with those measured from our data. Both are shown in Figure 7. We find that the model's total intensity and brightnesses are generally consistent with observations at all measured scattering angles, with the best agreement coming at intermediate angles of 30°–120°.

Zoom In Zoom Out Reset image size

In addition to the aperture-measured profiles, we plot the analytic scattering phase function B(θ) for total intensity calculated by MCFOST for this model. The output is in arbitrary units, so we scaled it uniformly such that it equals the observed NW brightness point at PA = 49°. Comparing this analytic profile with the aperture profile for the same model, we find B(θ) to be consistently lower at θ ≥ 60°—that is, the model's brightness is greater near the ansae and along the ring's back edge than expected, by more than 200% at some angles. We believe this excess brightness results from light scattered by the front and back edges being conflated due to viewing the inclined and vertically thick ring in projection.

As a test of this hypothesis, we calculated a "projected analytic phase function" B'(θ) that takes into account scattered light from one edge being projected onto the other edge. We first estimated the vertical distance Δz_j from the midplane of the ring at scattering angle θ_j to the midplane of the ring at the supplementary angle 180 – θ_j (i.e., the "opposite" edge). For all θ_j, we then computed the fractional density of dust contributed from the supplementary scattering angle as , where Σ = 1 at the supplementary midplane. The projected analytic phase function ends up as B'(θ_j) = B(180 − θ_j)·Σ_j, which we plot in Figure 7. It is a closer match to the measured model phase function than the original analytic phase function is, though it still underestimates the scattering somewhat (by up to 50% at 105° < θ < 125°). A localized peak occurs at θ = 90°, where the projection effect is at maximum. A second peak at θ ≈ 136° is produced by water ice preferentially scattering light incident at that angle, similar to the halo and Sun dog effects seen for sunlight in Earth's atmosphere.

We conclude that the scattering phase function measured directly from disk photometry is significantly impacted by projection effects and should not be taken at face value as the pure phase function, particularly for highly inclined and vertically extended disks. It is especially important to take this into account when comparing analytic phase functions with empirical measurements.

We modeled the disk's SED primarily as a check that our model's parameters were not ruled out by disk photometry at wavelengths outside of the near-IR. The data summarized in Table 3 comprise a spectrum from the Spitzer Space Telescope Infrared Spectrograph (IRS; Houck et al. 2004) and broadband photometric points from previous publications. This broadband photometry is composed of measurements from multiple optical and near-infrared instruments, NASA's Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010), the Multiband Imaging Photometer for Spitzer (MIPS; Rieke et al. 2004), and the Submillimeter Common-User Bolometer Array 2 (SCUBA-2; Holland et al. 2013).

In Figure 8 we plot the data alongside the SED produced from the median model from our MCMC. The SED was computed with MCFOST, assuming a stellar spectrum model from Kurucz (1993) with effective temperature of 6460 K, stellar luminosity of 2.4 L_⊙, log surface gravity of 4.0, and solar metallicity, based on SED fitting results from Moór et al. (2011b) and Chen et al. (2014), updated to the Gaia-derived distance. This is the same stellar model used throughout our modeling efforts.

Zoom In Zoom Out Reset image size

The median model SED is statistically consistent with the MIPS 160 μm point and upper limits at longer wavelengths but clearly deficient in flux at shorter wavelengths. This discrepancy may be due to another disk component not yet included in our model. This led us to model a separate inner disk component that, when added to the median model, would produce an SED similar to that observed. For simplicity, we manually tuned this inner component and consider it only a suggestion of one possible architecture for this disk.

Our best result assumes the inner component to be a narrow circular ring at r = 19–20 au containing 2.1 × 10⁻³ M_⊕ of dust. This inner component has the same inclination, fractional composition, and porosity as our median model but a larger minimum grain size of 10.0 μm and steeper size distribution with q = 4.5. Being ∼40 au closer to the star and roughly 50 times less massive than the median model, this component does not significantly impact our scattered-light models and would be indistinguishable from noise in our observations.

For comparison with the measured SED, we randomly selected 200 models from the MCMC chains to serve as the outer components and added the inner component to each to produce a distribution of two-component SEDs (Figure 8). We find that this distribution provides a good enough match to the data that our median ring model remains plausible, given a suitable inner component. The two-component SED falls within 2–3σ of all measurements apart from exceeding the 160 μm MIPS flux by 3.5σ and the 850 μm upper limit by ∼35%. The latter discrepancies stem from overproduction of flux by the median model at those wavelengths, which would be reduced if R_out is smaller than our loosely assumed 220 au and there is less cold dust in the ring as a result.

Our observations from ground- and space-based instruments, combined with MCMC modeling, have shown the HD 35841 system to include at least two, and perhaps three, debris disk components. Moving outward from the star, these are a hypothetical inner dust component, a primary dust ring, and a smooth halo extending outward from the ring. This configuration shares similarities with many other stellar systems, including our own.

The ∼57–80 au spatial scale of the primary dust ring makes it nearly two times larger in radius than the Kuiper Belt (located at ∼30–50 au; Levison et al. 2008). With HD 35841 being ∼20% more massive than the Sun and possibly hosting a narrow inner component akin to an asteroid belt, this system resembles a scaled-up version of the Solar System. We find this ring's scale height to be 4%–7% of its stellocentric radius, which is in line with measurements of 3%–10% for other disks like HR 4796A, Fomalhaut, AU Mic, and β Pic (Augereau et al. 1999; Kalas et al. 2005; Krist et al. 2005; Millar-Blanchaer et al. 2015). Despite the HD 35841 ring not being exceptionally thick, the measured scattering phase function is still significantly impacted by projection effects because the ring is highly inclined. Therefore, we reiterate that this is an important aspect to consider for phase function measurements of other disks. This issue also highlights the value of polarimetry, which avoids PSF subtraction biases and enables additional constraints on disk geometry and phase functions.

A comparison of our models with the observed SED implies the existence of a second dust component interior to the ring imaged by GPI and STIS. Our brief exploration of solutions shows this inner component could be a narrow ring at ∼20 au (02), which is just interior to our high SNR region in the GPI data but still outside of GPI's H-band inner working angle. It contains less mass than the main ring, but it will receive more incident flux, buoying the scattered-light brightness. Future direct imaging with a smaller inner working angle may resolve this dust. Interferometric observations may also help constrain it.

We note that our disk components differ from those of previous SED fits of HD 35841. Chen et al. (2014) fit blackbody rings at separations of 45 and 172 au (for a stellar distance of 96 pc) to the IRS and MIPS data. Moór et al. (2011b) used just a single infinitesimally narrow ring of modified blackbodies at ∼23 au (also for d = 96 pc). However, neither of those studies benefited from spatially resolved imaging, which requires dust at 57–80 au. Regarding their placement of dust interior to our primary ring, it is common to measure a resolved disk radius that is greater than the radius predicted by blackbody approximations, as Morales et al. (2016) demonstrated for Herschel-resolved disks. On the other hand, the outer component from Chen et al. (2014) is located just beyond the outer part of the detected halo and could represent that material.

One final morphological feature of the disk not represented in our MCMC models is the outermost part of the smooth halo extending from the ring in the STIS image. Similar features are seen in other disk images—for example, of HD 32297, HD 61005, and HD 129590 (Schneider et al. 2014; Matthews et al. 2017). These halos may be populated by ring grains that are slightly larger than a_blow and are excited onto highly eccentric orbits by radiation pressure (Strubbe & Chiang 2006). This type of eccentric grain population would cover a large surface area and scatter substantial light, but contain little mass and be relatively far from the star; thus it would contribute little to the overall SED. This is qualitatively similar to what we find for HD 35841. One could better test the link between the ring and halo populations with more holistic models that directly incorporate radiation pressure into the model physics. However, theory predicts that the bound grain population should also leave an observational signature in the halo's surface brightness radial profile, which we test below.

We measured surface brightness radial profiles for the smooth halo to see if their slopes agree with that expected for the bound, eccentric population of grains we proposed in the previous section.

To measure the radial profile, we placed apertures (radii = 2 px) along the ring's presumed major axis in the interpolated STIS image. The innermost apertures were centered at r = 51 px (74.3 au), just exterior to the ansae, and were located every 5 px out to r = 237 px (346 au) on both sides of the star. Surface brightnesses and their uncertainties were measured in the same way as the phase functions in Section 3.2. Measurements consistent with zero at less than the 1σ level are discarded, which includes all points at r > 210 au. We then fit power law functions of the form independently to the radial profiles for each side of the halo using a least-squares minimizer algorithm.

The halo radial profiles and best-fit power law functions are plotted in Figure 9. We found power law indices of α_h = −2.80 ± 0.36 in the SE and α_h = −4.18 ± 0.37 in the NW. Continuing to assume that the disk is azimuthally symmetric, we also fit a single power law to all measurements from both sides of the halo, and found an index of α_h = −3.55 ± 0.35. According to Strubbe and Chiang (2006), the surface brightnesses of collisionally dominated debris disks will vary with radius as SB ∝ r^−3.5 exterior to the dust-producing "birth ring" of planetesimals. Indices for individual sides of the halo agree with this predicted value to within 2σ and the joint index matches it nearly exactly. This is strong evidence that the halo's grains are collisionally produced in the ring and remain gravitationally bound to the star on wide and eccentric orbits. In that case, the brightness in our model from dust on circular orbits with a ≳ 74 au may simply be a proxy for this eccentric population.

Zoom In Zoom Out Reset image size

A number of our conclusions about the disk's grain properties are subsequently related to the greater disk environment, which links to planet formation and evolution. However, extracting results for the grain properties required making assumptions about the underlying scattering physics in our disk model, which in our case is based on Mie theory. Therefore, before examining our results more closely, we consider how this assumption may have affected them.

Our primary concern is that Mie theory may not accurately reproduce the scattering phase function for grains in debris disks. This would not be surprising given that disk grains, born via collision, are almost certainly not homogeneous perfect spheres like the idealized theory assumes. In particular, Min et al. (2016) found that, for equivalent grain radii and porosities, the Mie phase function generally decreases with scattering angle θ > 90°, while a phase function for irregularly shaped aggregate grains is flat or increasing. Therefore, a model exploration that is based on a Mie model may bias parameters controlling grain size, porosity, and composition away from their true values in order to match backward-scattering at large angles seen in comparison data. This would skew the resulting posterior distributions. Aggregate grains, on the other hand, would more naturally produce a backward-scattering component and may prefer different parameter values that are closer to reality. Milli et al. (2017) recently pointed out these effects for the debris ring around HR 4796A, for which a Mie model was also incompatible with the scattering phase function. It is important to keep in mind these shortcomings of Mie theory (and other simplified scattering treatments) for the following discussions of grain properties and for future disk modeling efforts.

Regarding grain structure, our models show a clear preference for a low porosity (<12% at 99.7% confidence). The disk's polarized intensity is particularly constraining in this regard, as higher porosity tends to increase polarization fraction for a given grain size and composition. Higher porosity also makes grains more forward scattering, so our view of the ring's back edge in total intensity provides additional constraints. This low porosity is in contrast to models of the AU Mic disk from Fitzgerald et al. (2007) and Graham et al. (2007) that require highly porous (80% vacuum) comet-like grains to reproduce its scattering and polarization signatures. This discrepancy may arise from differences in grain size; the Fitzgerald et al. (2007) best-fit model with 80% porosity only contained grains with 0.05 μm < a < 3.0 μm or 3 mm < a < 6 mm. Nevertheless, one interpretation of our result favoring compact grains is that little cometary activity occurs in the HD 35841 system. This may be borne out by the non-detection of CO (J = 3 − 2) from Moór et al. (2011a). Deeper searches for CO emission and other gas signatures would further investigation on this topic.

In terms of minimum grain size, a_min ≈ a_blow in most of our models. Our uniform Bayesian prior did not contain any information about a_blow, so agreement of this independent result with the fundamental physics of the system lends credence to the other parameters that we have derived from the light-scattering and polarization signatures. It also supports the case for non-porous grains, as higher porosity leads to greater surface area and a larger blowout size. Meanwhile, modeling of some debris disks, like HD 114082 (another F-type star) by Wahhaj et al. (2016), indicate an a_min several microns larger than a_blow. Perhaps this is a true dichotomy resulting from differences in grain porosity and/or shape, or maybe the difference is purely a result of inconsistent modeling methodology. A unified modeling effort applying the same methodology (or better still, a range of methodologies) to multiple disks would provide valuable insight but would also be a substantial undertaking.

Another key aspect of our model grain population is the size distribution slope, which we find to be 2.6 ≲ q ≲ 3.2 in the 99.7% confidence interval, with a median value of q = 2.9. These values lie just below average values of 3.36 (MacGregor et al. 2016) and 3.15–3.26 (Marshall et al. 2017), recently estimated for several disks based on their mm-wavelength emission. The thermal brightnesses of those disks are dominated by mm-sized grains, whereas our scattered-light brightness is dominated by micron-sized grains. The fact that our micron-appropriate q values are similar to mm-appropriate q values implies that collisional cascades proceed in a self-similar (read: single power law) way across this size range—that the physics determining particle strengths and particle velocity dispersions do not change qualitatively from millimeter to micron sizes (but see Strubbe & Chiang 2006 for why the size distribution deviates strongly from a single power law at sizes that are just above the radiation blowout limit). Future observations of the mm emission from HD 35841 would be useful for verifying that its particle size distribution is indeed characterized by a single power law from millimeters to microns.

The least constrained of our grain properties are their compositions. Though our models show at least a 2:1 preference for amorphous carbon over astrosilicates and roughly one-third of the mass in water ice, the distributions are fairly broad (see Appendix B). Given the resulting uncertainties and influence of Mie approximations, we caution against assigning too much significance to these results.

That said, we can speculate on their implications. For example, the water ice produces a backscattering peak centered around θ ≈ 135° that locally enhances the ring's back edge. Future measurements of the phase function with small uncertainties would help confirm this as a real feature. With the low grain porosity implying little cometary activity, the presence of substantial water ice in the disk would need to be explained another way. As for the other materials, the scattering properties of silicate and carbonaceous grains within a single near-IR filter band are very similar apart from albedo. Examples of both cases have been presented in studies of other debris disks. Though difficult, differentiating between a silicate-rich and a carbon-rich disk would be meaningful for the compositions of planets in the system, which presumably formed in and collected material from the same resource pool. A system abundant with carbon and water, two key materials for life on Earth, would be especially interesting from an astrobiology perspective.