Multiband Polarimetric Imaging of HR 4796A with the Gemini Planet Imager

The Gemini Planet Imager (GPI) is a high-contrast imaging instrument built for the Gemini Observatory; it has been operating at the Gemini South Telescope (Macintosh et al. 2006). GPI has an integral field polarimetry mode as well as a integral field spectrograph mode. In the polarimeter mode, light is sampled in the pupil plane by a lenslet array followed by a polarizing beam splitter. The raw image consists of an array of spatial resolution elements or spaxels, with pairs of spots of orthogonal linear polarizations. The light is modulated by a half-wave plate (HWP) between each exposure, such that I, Q, and U images can be constructed from a sequence of exposures.

We observed HR4796A with GPI on 2014 March 22. A summary of the observations are listed in Table 1. Thirty-five 60 s exposures were taken in J-band (${\lambda }_{c}=1.24\mu {\rm{m}}$) polarimetry mode with 65° of field rotation, followed by thirty-eight 60 s exposures in the K₁-band (${\lambda }_{c}=2.05\mu {\rm{m}}$) mode with 438 of field rotation under good seeing conditions. The half-wave plate was rotated between position angles 0°, 22°, 45°, 68° throughout each sequence. We additionally used H-band (λ_c = 1.65μm) polarimetry mode data whose acquisition and reduction is described in Perrin et al. (2014).

The raw data were reduced using the GPI Data Reduction Pipeline (Perrin et al. 2014). The raw images were dark-subtracted, flexure-corrected, destriped, and corrected for bad pixels. The orthogonal polarization spots were then extracted from each raw image to form polarization cubes, each with two frames of orthogonal polarization. The cubes were then divided by a polarized flat field. Bad pixels were identified and replaced via interpolation. The star's position was measured using the position of reference satellite spots diffracted from starlight behind the coronagraph (Wang et al. 2014).

Data taken in GPI's polarimetry mode are particularly suited for both polarized and angular differential imaging, both of which suppress the starlight and improve the contrast by orders of magnitude. For polarized differential imaging (PDI), we subtracted the two frames of orthogonal polarization for each datacube, removing the majority of the starlight, which has a randomly oriented polarization. Stokes cubes were then constructed from the resultant frames using a singular value decomposition method (Perrin et al. 2015) to recover Stokes parameters from the data and HWP-modulated time-variable measurement matrices. The mean stellar polarization was corrected for by first measuring the average polarized intensity (Stokes parameter I) inside the focal plane mask. The total intensity image scaled by the measured intensity was then subtracted from the linearly polarized intensity image. The final image that was fit to in subsequent modeling described in Section 3 was a radial polarization image, a combination of the Q and U images that gives the polarization in the radial direction. For an optically thin single scattering disk, all of the signal is expected to lie in this radial polarization.

Another advantage of using these polarization data is that the sum of the linear polarization states can be combined and processed using an angular differential imaging (ADI; Marois et al. 2006) algorithm, to produce a PSF subtracted total intensity image. For each data cube, we combined the two linear polarization states to form a series of total intensity images to correspond to each polarization image. We then used a Python implementation of Karhunen–Loeve Image Projection (KLIP; Soummer et al. 2012), PyKLIP (Wang et al. 2015), to perform this angular differential imaging. We optimized the size and number of subtraction regions, as well as the number of basis vectors subtracted, to minimize PSF self-subtraction of the disk by making measurements of the signal-to-noise ratio (S/N) at various points along the disk as a function of KLIP parameters. Our measurements indicated the optimal parameters were one basis vector and one subtraction.

Our forward model as described in Section 2 requires a convolution of our model with a PSF. The PSF for GPI is challenging to model, due to its complex structure and variability (Wang et al. 2014). As such, rather than use a Gaussian or Airy function, we used a PSF that extracted from the four satellite spots dispersed in each image of HR4796A. Since polarimetric frames are broadband—and therefore have overlapping satellite spots in a single frame—we extracted the PSF structure from satellite spots of observations taken in GPI's integral field spectrograph (IFS) mode. We elected to use the HR4796A satellite spots even though the field is noisier than that of observations of other stars because the of the dependence of the PSF shape on the stellar spectrum. The stellar spectra would affect the relative weights of the extracted satellite spots at different wavelengths. We used thirty-six 60 s H-band frames taken on 2016 March 18, fifty-eight 30 s J-band frames taken on 2016 March 23, and thirty-six 60 s K₁-band frames taken on 2016 March 13.

These data were also reduced using the GPI Data Reduction Pipeline (Perrin et al. 2014). The raw images were dark-subtracted, flexure-corrected, destriped, and corrected for bad pixels. The spectra for each spaxel were then extracted to form 3D data cubes. The data cubes were then further corrected for bad pixels and distortion.

In order to estimate the PSF, we first summed the spectral mode images along the wavelength axis and the polarimetric frames on the polarized axis, both giving estimates of the total intensity across the bandpass. We took the median image of these flattened spectral mode images and the median of the polarimetric mode images. We high-pass filtered both median images with a box size of 6 spaxels, in order to remove large-scale structure from the main image of the star behind the occulting PSF. This box size was chosen to optimize the uniformity of the background structure surrounding the star. Each spectral channel was linearly combined with a weight to find the least-square difference with the polarimetric satellite spot. The need for different weighting parameters stems from the difference in throughput between spectral and polarimetric mode. To get an image of the PSF, we registered the spectral satellite spots of each wavelength, multiplied each one by the weights we had fitted for, and summed them along the wavelength axis. The PSF is highly asymmetric, with lobes at four locations around the core. Thus, for our model convolution, we azimuthally averaged this PSF because each image of the disk is derotated for sky rotation. making the final PSF a combination of rotated PSFs from individual exposures.

In order to extract the geometric parameters and brightness function of the disk, we fit a geometric model to the data. By fitting a model generated purely from an arbitrary description of phase and geometric parameters, we remain agnostic to any assumptions about the physical forces on the dust grains, the orbital grain distribution, or the properties of the dust grains. In this procedure, we additionally use KLIP forward modeling (Pueyo 2016) to account for self-subtraction of the disk brightness in total intensity.

We selected our preferred prescription for the disk by minimizing the number of parameters needed to achieve comparable χ² values. We found that modeling the disk as an ellipse (as opposed to an offset circle) added extra parameters that did not improve the χ² sufficiently to warrant the more complicated disk. We therefore modeled the ring as a series of nested circles offset from the star. We fit for Ω, the position angle of the major axis, and i, the inclination of the nested circles.

We constructed the model by mapping each pixel location $(x,y)$ in the sky plane to a radius $r(x,y)=\sqrt{{x}_{\mathrm{disk}}{\left(x,y\right)}^{2}+{y}_{\mathrm{disk}}{\left(x,y\right)}^{2}}$ in the disk plane and a $\theta (x,y)={\tan }^{-1}(x/y)$ in the sky plane. The intensity of each pixel is then

$$\begin{eqnarray}&&I(x,y)=\displaystyle \frac{{B}_{r}(r){B}_{\theta }(\theta )}{{r}_{* }^{2}},\end{eqnarray} \tag{ 1 }$$

where B_r is a radial profile, B_θ is the azimuthal intensity profile, and r_* is a unitless factor that scales with the distance between the location (x, y) and the star. Here, B_θ is a periodic spline interpolation with varying numbers of knots, with the intensity at every knot as a free parameter. Transforming B_θ to B_ϕ, where ϕ is the scattering angle, gives the SPF of each model disk. As all of the resulting phase curves are normalized in our analysis in Section 5, the units on I(x, y) are arbitrary.

The knots are evenly spaced in the sky plane along the disk, as shown in the blue points in Figure 1. By spacing the points with a separation larger than the scale of the PSF, we minimized the effects of spatially correlated noise residuals. To find the optimal number of spline points, we used the dust modeling package MCFOST (Pinte et al. 2006a) to generate a model with similar geometry to HR4796A and a known two-component Henyey–Greenstein phase function with the dust modeling package MCFOST. We then injected this model into a separate polarization data set with no disk detection. We used our modeling procedure to recover this artificial disk. Transforming the recovered B_θ to B_ϕ gave a curve that could be directly compared to a theoretical B_ϕ scattering phase function. Using a minimum number of 13 knots, we were able to recover the scattering phase function to the 1% level at all observable scattering angles. In this test injection, the fit was invariant to the location of the spline points as long as there were a sufficient number of points to fully describe the shape. However, this test assumes that the disk has a smooth phase function at all points—which, as will be shown in Section 5, may not be the case. It also assumes that there is no large-scale noise structure that the model would fit to, which fortunately is the case for most of the images we have modeled. Though the locations of the spline points are not densely sampled in ϕ versus B_ϕ space, as long as the intrinsic SPF is smooth, the spline will recover its shape at all accessible scattering angles. This is in contrast to extractions of an SPF that use aperture photometry to sample to the brightness B_θ, which can only be described by brightnesses at the discrete locations of the apertures and suffer from self-subtraction bias.

Zoom In Zoom Out Reset image size

The radial profile B_r is a broken power law:

$$\begin{eqnarray}{B}_{r}(r)=\left\{\begin{array}{ll}{\left(\displaystyle \frac{r}{{r}_{0}}\right)}^{-{\gamma }_{\mathrm{in}}} & \mathrm{if}\ {r}_{\mathrm{in}}\lt r\lt {r}_{0}\\ {\left(\displaystyle \frac{r}{{r}_{0}}\right)}^{-{\gamma }_{\mathrm{out}}} & \mathrm{if}\ {r}_{0}\lt r\lt {r}_{\mathrm{out}}\\ 0 & \mathrm{otherwise}\end{array}\right.,\end{eqnarray} \tag{ 2 }$$

where r₀ represents the central radius in milliarcseconds, γ_in the inner radial profile, and γ_out the outer radial profile, as free parameters. The radial profile was found to be very sharp, to the point that the γ factors were degenerate with the inner and outer radii. To reduce the number of parameters and avoid unbounded parameters, r_in and r_out were fixed at 70 and 100 au, respectively. These radii were selected by deprojecting the disk and finding the radii where the S/N of the disk falls below 10%.

We then used our model disks to fit to the J- and K₁-band total and polarized intensity images, as well as the H-band polarized intensity image. The H-band total intensity did not have enough field rotation to reliably be forward modeled. We created model images with the above parameters, which we then convolved with our derived convolution PSF. In order to simulate the effects on the final data product due to the KLIP PSF subtraction, we developed the DiskFM module for PyKLIP (Wang et al. 2014) specifically for forward modeling extended objects based on the mathematical framework presented in Pueyo (2016). Due to this extra step of modeling, we fixed the geometrical parameters of the total intensity disks to those found from their polarized intensity counterparts and only fit the scattering phase spline function. This is a natural choice because the total and polarized intensity images are generated from the same raw data images.

We fit each disk's geometric parameters independently from band to band, to account for various physical and nonphysical effects. The position angle of the line of nodes (Ω) and inclination (i) could differ between bands, as there is some uncertainty in the rotation of the instrument relative to north. The radial profile parameters were fitted separately to reflect possible changes between the distributions of differently sized dust grains due to differing effects of radiation pressure and gravitational forces.

We independently fit the by using a linear least-squares algorithm to optimize the χ² using the uncertainty maps. We then used the resultant parameters to seed a fit using an ensemble Markov chain Monte Carlo (MCMC) using the emcee package (Foreman-Mackey et al. 2013). The final geometric parameters are shown in Table 2, with their error bars derived from the distributions of the final walkers. After fitting for the geometrical parameters, we fixed all of the geometrical parameters for each model disk and fit only the spline parameters.

The data and best-fit models for the polarized intensity data are shown in the left and middle columns of Figure 2. The third column shows the difference between the final and data images divided by our noise map. The residuals for the H- and K₁- band model are consistent with the estimated data uncertainties. Some residual structure may be seen in the J-band image northeast ansa, which we will later discuss in Section 4.2.

Zoom In Zoom Out Reset image size

Fits to the total intensity data are shown in Figure 3. The northwest portion of the model in the J-band image overfitted to speckle noise, most evident in the image of the model without forward modeling, shows an unphysical dip in intensity. This portion of the disk shows a dip in intensity where it would be expected that the intensity would be smooth. Because of this, we have decided to omit the J-band polarized fraction and total intensity curves from the phase curve analysis in Section 5.

Zoom In Zoom Out Reset image size

The final distributions of the MCMC walkers for the K₁-band fit are shown in Figure 4, and the best-fit parameters with 3σ variance for all bands are listed in Table 2. It is evident from both of the walker distributions in Figure 4 that the variance of the final parameters are unrealistically small, most likely due to some model mismatch. In the final values for the PA and ω (the direction of the offset), listed in Table 2, we have included the variance of the image from true north of −1° ± .001° found by Konopacky et al. (2014). Calculations of the radius in milliarcseconds shown in the last line of the table have included the error in assumed plate scale of 14.14 ± 10⁻⁵ milliarcseconds.

Zoom In Zoom Out Reset image size

Table 2 also shows comparisons to parameters found by Milli et al. (2017) in H band and Schneider et al. (2014) in the F25ND3 filter. Milli et al. (2017) found their geometric parameters by fitting radial profiles to cuts of the image and fitting ellipses through the maximal radial values of every profile. To compare with our circular model, Table 2 shows the average distance of every point along the ellipse to the star with their best-fit parameters of a semimajor axis $a=1066\,\mathrm{mas}$ and ellipticity e = 0.07. As the J-band image has strong residual structure that is likely driving the fit parameters, it is more useful to compare parameters between K₁- and H-band parameters. The geometry parameters of the position angle and inclination in these bands are consistent not only with each other but also with Milli et al. (2017). The radii found in Milli et al. (2017) are consistent with our derived K₁- and J-band models, but not with the H-band model. This may be due to biasing of fit by noise close to the focal plane mask in the H band. Overall, the most consistent and reliable geometric measurements come from the K₁band.

The inner radial profile exponents γ_in are large, indicating an unresolved inner edge. The error bars are unrealistically small due to nonuniformity in the radial profile, with one ansa forcing a steeper radial profile and the other forcing a broader radial profile. This effect is most visible in the J-band residual image in the northern portion of the disk in the upper left panel of Figure 2. There are residuals between the data and the model exactly at the midplane, and a visual inspection of the J-band data and model show that this is likely due to the data's radial profile being sharper in that region. A more direct representation of the radial profile can be seen in Figure 5, in which we plotted the intensity radially along the cuts in the directions shown in Figure 5(d). The K₁-band radial profiles in Figure 5(a) show little systematic deviation between the model and data. In the H-band radial profiles, the radial profiles near the ansae are well-fit, though there is evidently noise at small scattering angles near the focal plane mask. This is a likely explanation for the small radius in the H-band fit. Figure 5(c) shows that the southwest radial fits are good in the J band, but the peaks of the model cuts are systematically lower than those of the data in the NE region. The radial profile fit is forced by the inner and outer sides of the profiles, lowering the peak. As the model's radial profile is uniform about the disk, this indicates that the radial profile in the northeast half of the disk is sharper than in the southeast half. Such an effect would be seen most evidently in the J band, as it has the smallest PSF and highest resolution. Qualitatively, the narrowing of the disk in the northwest side is consistent with the Olofsson et al. (2019) measurements of polarized SPHERE/ZIMPOL data. We refer the reader to Olofsson et al. (2019) for an in-depth discussion of the physical mechanisms that may be causing this effect.

Zoom In Zoom Out Reset image size

The polarized intensity curves are shown in Figure 6, normalized at a scattering angle of 90°. The scattering phase functions are strongly forward-scattering, with both the polarized and total intensity phase curves peaking at the smallest scattering angles. The NE and SW curves in the H and K₁ bands are symmetric, while the NE ansa of the J-band image has a bump at a scattering angle of 55° due to the residual structure seen in the images at this scattering angle. While the phase curves have similar behavior from 70° − 120°, the heights of the peaks vary with wavelength. The K₁-band phase curve (λ_c = 2.05μm) evidently has a sharper forward scattering peak than the H band's (λ_c = 1.65μm), which is sharper than the J band's (λ_c = 1.12μm). Polarized intensity phase curves taken by ZIMPOL at λ_c = 0.74μm show that this effect extends to smaller wavelengths, with the phase curve similarly decreasing from 80° − 120° but plateauing from 13° to 80° (Milli et al. 2019). The source of this chromaticity is unknown, as it is plausible for the effect to be caused by a different spatial distribution of multiple grain populations or chromatic effects of a single dust population. Since we are only analyzing the polarized intensity and not the polarized fraction of the J- and H-band data sets, consistency cannot be checked for chromatic effects with modeling.

Zoom In Zoom Out Reset image size

The K₁-band total intensity and polarized fraction curves are shown in the middle and bottom panels, respectively in Figure 7. The total intensity curves were normalized at 1 at 90°, while the polarized fraction is unitless. Though it is challenging to measure the polarized and total intensities in physical units, the unnormalized curves can be divided to calculate the polarized fraction. Consistent with phase curves in a similar band in Milli et al. (2017), the total intensity phase curve exhibits a forward-scattering peak and a flat distribution that rises at scattering angles larger than 70°. The polarized fraction curve peaks at ∼40° at 50%, consistent with the lower limit found by Hinkley et al. (2009) of 44% as well as the peak polarization of 50% at a scattering angle of 50° found by Perrin et al. (2014).

Zoom In Zoom Out Reset image size

We used the MCFOST package (Pinte et al. 2006b) to generate theoretical Mie and DHS phase functions (Min et al. 2005) to fit to our measured phase functions. We modeled to our highest-fidelity curves: the southwest K₁-band total intensity and the polarized fraction. We used MCFOST to compute total intensity and polarized fraction phase curves using a given set of parameters at the central wavelength of the K₁-band filter. Because the change in grain properties over the K₁-band filter is small for most materials, integrating over the whole band did not significantly affect the morphology of the curve. For the total intensity curves, we compared the data to a scaled model where we found the scaling factor by taking the ratio of the model and data curves at every scattering angle and taking the median of those ratios.

Using the scaled total intensity curve and the polarized fraction, we then computed reduced ${\chi }_{\nu }^{2}$ values for each curve. As the profiles were generated from a previous fitting procedure, we expected the errors to be correlated, but that was ignored in this χ² calculation. In the total intensity fit, we excluded regions at scattering angles smaller than 20° and larger than 140°, as the data are unreliable close to the focal plane mask (shown in Figure 7). We truncated the polarized fraction curve past 120°, as the signal of the polarized intensity is undetected. The locations of these cut-off scattering angles with respect to the disk are shown in Figure 1.

We ran a grid search over the minimum grain sizes a_min, the exponent of the power law that describes the grain size distribution a_exp and the grain composition. We assumed a grain size distribution of:

$$\begin{eqnarray}\displaystyle \frac{{dN}}{{da}}\sim \left\{\begin{array}{ll}{a}^{-{a}_{\exp }} & \mathrm{if}\ {a}_{\min }\lt a\lt {a}_{\max }.\\ 0, & \mathrm{otherwise}.\end{array}\right.\end{eqnarray} \tag{ 3 }$$

We parameterized the grain composition in terms of the real and imaginary component indices of refraction of the dust grains. By doing so, we remain agnostic to the chemical composition of the grains. We also eliminate the need for the porosity parameter, whose effects are captured by the real and imaginary indices of the dust grain population, assuming a uniform effective medium.

The ranges of our fitting parameters are shown in Table 3. The phase curves were integrated over a range from a_min to a_max. The maximum grain size, a_max, was fixed at 1 mm due to the sharp power law that dictates there are few large grains for any of the proposed grain size distribution. The limits of the real and imaginary indices of refraction were gained from the ranges of the indices for physical grain compositions at the K₁-band wavelength. Measured real and imaginary indices for a variety of different materials at K₁band are shown in Figure 11. Whereas the usually assumed exponent a_exp for the grain size power law distribution is usually assumed to be 3.5, following Mathis et al. (1977), we fit over 10 different power laws.

We evaluated both DHS and Mie models for the grid defined in Table 3. We examined the resultant curves for each model, using the metrics of the lowest χ_ν for the total intensity curve, lowest χ_ν for the polarized intensity curve, and lowest χ_ν for both curves simultaneously. The ideal model needs the three distinctive properties of the HR4796A model: a strong forward-scattering peak in the total intensity curve, a gradual increase in the total intensity curve at the backscattering side, and a peak in polarized intensity at 40°.

We found that neither model could simultaneously reproduce all of the features of both the SPF and polarized phase function. Figure 8 shows the best-fit models for the simultaneous χ_ν. Though not a close model in total intensity, the DHS model is able to reproduce the features of the forward-scattering peak as well as the shape of the curve in the backscattered direction. On the other hand, the DHS model is unable to reproduce the peak in polarized fraction at 40°.

Zoom In Zoom Out Reset image size

While the Mie model generates a polarized fraction curve with a peak closer to that of the data, the Mie model fails to recover the magnitude of the peak in total intensity as well as the increase in intensity at backscattering angles, exhibiting instead a flat distribution at scattering angles greater than 40°.

We computed the goodness-of-fit metrics for the total intensity and polarized fraction phase curves independently of each other by calculating the ${\chi }_{\nu }^{2}$ of each while ignoring the other. The models with the lowest χ_ν of the total intensity phase curves are shown in Figure 9, with the parameters listed in Table 4. In this case, the best-fit Mie model is able to reproduce the backscattering increase, but cannot produce a forward-scattering peak sharp enough to match the model. The DHS model has a good fit to the overall curve, with a ${\chi }_{\nu }^{2}$ under 1. We compare the DHS ${\chi }_{\nu }^{2}$ of the total intensity–only fit versus the ${\chi }_{\nu }^{2}$ of the combined polarized fraction and total intensity fits (shown in rows 5 and 4, respectively). This comparison reveals that the total intensity fit has an improvement on the total intensity χ², but a drastically worse polarized fraction.

Zoom In Zoom Out Reset image size

Table 4.  Best Fit Parameters for Different ${\chi }_{\nu }^{2}$ for Different Grain Models

Grain Model	Metric ${\chi }_{\nu }^{2}$	Summed ${\chi }_{\nu }^{2}$	Tot Intensity ${\chi }_{\nu }^{2}$	Pol Frac ${\chi }_{\nu }^{2}$	${a}_{\min }$	aexp	n	k
Mie	Sum	12.4	13.3	11.8	1.9	2.9	3.7	3.72
Mie	Tot Intensity	155.8	3.3	363.1	0.3	5.7	1.1	10.0
Mie	Pol Frac	165.2	293.6	3.0	3.7	3.7	2.0	2.7 × 10−2
DHS	Sum	11.5	5.2	20.4	13.9	3.3	3.4	3.7
DHS	Tot Intensity	41.2	0.8	96.2	26.8	4.1	1.1	1.4
DHS	Pol Frac	815.3	1458.8	3.0	1.9	2.9	3.4	3.0 × 10−2

Note. The bolded row indicates the ${\chi }_{\nu }^{2}$ that each set of parameters was optimized for. The third column shows the sum of the total intensity and polarized fraction ${\chi }_{\nu }^{2}$ values. The fourth and fifth columns are the ${\chi }_{\nu }^{2}$ values for the total and polarized fraction, respectively. The second column describes the metric over which the best-fit parameters were derived. The second and fifth rows list the best-fit parameters in the sixth to ninth columns for the best summed ${\chi }_{\nu }^{2}$. The third and sixth columns are the parameters for the best total intensity ${\chi }_{\nu }^{2}$, and the fourth and sixth list the best-fit parameters for the best polarized intensity.

Download table as:  ASCIITypeset image

The curves produced by fitting only to the polarized fraction are shown in Figure 10. Both curves have overall structures similar to the data phase curves, with a peak at 40°, but they exhibit an unexpectedly jagged curve. Images produced by MCFOST with these phase curves do not visibly show any of this roughness, given the pixel sampling and PSF convolution. Our model, constructed assuming a smooth phase curve, would therefore be unable to detect any extra structure on the curve without overfitting.

Zoom In Zoom Out Reset image size

In order to further evaluate the generated DHS model in comparison to more physical models, we compared the phase space of likely indices of refraction we derived from our fits to indices of various other materials in Figure 11. Following Bruggeman mixing rules, mixtures of two or more materials result in indices intermediate to the indices of the materials being mixed. A mixture of any dust compositions would lie somewhere along the semilinear track traced out by the materials already shown. Porosity, essentially a mixture of void with a dust grain composition, would additionally move any point along the same track.

Zoom In Zoom Out Reset image size

The polarized fraction DHS best fits, boxed in red, occupy a part of parameter space that is not only far from any pure dust grain composition, but would also be far from any mixture with any porosity. The parameter space of decent total intensity fits using DHS is fairly broad and overlaps with the track of physical compositions. However, the lack of overlap between the polarized fraction fits and the total intensity fits precludes any confident conclusions about the grain composition derived from the DHS fits.

Both the Mie and DHS models are meant to be substitutes for more realistic—but more computationally expensive—models of aggregate dust grains. These aggregate dust grain models get exponentially more expensive with grain size. Our models that produce the smallest ${\chi }_{\nu }^{2}$ values all exhibit large grain sizes of 2–26 μm, for which aggregate models have not been extensively generated. This analysis questions the validity of Mie and DHS models in producing meaningful results in this limit. The phase curves for HR4796A are unlike other phase curves in the defining features of the sharp total intensity forward-scattering peak at 25°, the modest backscattering peak, and the polarized intensity peak at 25°. Neither model was able to fully produce all three features simultaneously.

There are a number of other ways that the dust population model can be improved in future work regarding the parameterization of the size distribution. Most obviously, the sharp cutoffs of our dust grain size distribution at the minimum and maximum grain sizes are unphysical. Given the steep outer power-law a_exp and our large maximum grain size, it is unlikely that increasing the maximum grain size cutoff would appreciably affect the resultant model. On the other hand, creating a more gradual distribution of grains rather than one that sharply cuts off at the minimum grain size would likely affect the model phase curves. It is also likely that the jagged polarized fraction phase curves would be smoothed by the inclusion of smaller grains, but not without affecting the goodness of fit to the DHS total intensity model. Another major possibility is that not only is there a mixture of grains with different compositions, but there are also multiple dust grain populations with different size distributions.