THE AU MIC DEBRIS DISK: FAR-INFRARED AND SUBMILLIMETER RESOLVED IMAGING

Large area photometric mapping observations were performed using the PACS and SPIRE cameras at 70, 160, 250, 350, and 500 μm. The large scan-map mode was used and nearly perpendicular cross-linked scans were performed with both instruments. The total observing time was 3.8 hr. The images at each observed wavelength are shown in Figure 2.

Zoom In Zoom Out Reset image size

The PACS observations used scan legs of 74 with a scan-leg spacing of 38''. Medium scan rate maps containing 11 scan legs were repeated 11 times in each scanning direction to achieve the required map depth. Data were obtained at both 70 and 160 μm simultaneously.

The nominal cross-linked large map observing parameters were used for the SPIRE observations, with a smaller map of 4 arcmin by 4 arcmin being executed. In total, 10 repeat maps were used to reach the required depth.

The data were processed using version 13.0 of the Herschel interactive pipeline environment (HIPE, Ott 2010). The standard pipeline processing steps were used for both the PACS and SPIRE data. Versions 69 and 13.1 of the PACS and SPIRE calibration products were applied respectively. The final maps have pixel scales of 1, 2, 6, 10, and 14'' at 70, 160, 250, 350, and 500 μm respectively.

As part of the data processing, the PACS timelines were first masked at the source location, as well as in other areas of bright emission, and then high-pass filtered to reduce the impact of 1/f noise. The filtered data were then converted to maps using the "photProject" task. Maps were likewise made from the SPIRE data using the "naiveMapper" task. No filtering was performed on the SPIRE data.

Observations of α Tau and Neptune, adjusted to the correct spacecraft position angle, were used as the model PACS and SPIRE point-spread functions (PSFs) respectively, and used for model image convolution.

Observational data presented in this paper were also taken using the SCUBA-2 camera (Holland et al. 2013) on the JCMT. The data were obtained both as part of the SONS JCMT Legacy survey (Phillips et al. 2010) and a PI program. The observations were taken with the constant speed DAISY pattern (Holland et al. 2013), which maximizes exposure time and provides uniform coverage in the central 3' diameter region of a field. The total integration time was just over 5 hr, split into 10 separate ∼30 minute observations. Observing conditions were generally excellent with precipitable water vapour levels less than 1 mm, corresponding to zenith sky opacities of around 1.0 and 0.2 at 450 and 850 μm respectively (equivalent to JCMT weather "grade 1"; ${\tau }_{225\;\mathrm{GHz}}$ of $\;\lt 0.05$). The data were calibrated in flux density against the primary calibrator Uranus and also secondary calibrators CRL 618 and CRL 2688 from the JCMT calibrator list (Dempsey et al. 2013), with estimated calibration uncertainties amounting to 10% at 450 μm and 5% at 850 μm.

The SCUBA-2 data were reduced using the Dynamic Iterative Map-Maker within the STARLINK SMURF package (Chapin et al. 2013) called from the ORAC-DR automated pipeline (Cavanagh et al. 2008). The map maker used a configuration file optimized for known position, compact sources. It adopts the technique of "zero-masking" in which the map is constrained to a mean value of zero (in this case outside a radius of 60'' from the center of the field), for all but the final interation of the map maker (Chapin et al. 2013). The technique not only helps convergence in the iterative part of the map-making process but suppresses the large-scale ripples that can produce ringing artefacts. The data are also high-pass filtered at 1 Hz, corresponding to a spatial cut-off of ∼150'' for a typical DAISY scanning speed of 155'' s⁻¹. The filtering removes residual low-frequency (large spatial scale) noise and, along with the "zero-masking" technique, produces flat, uniform final images largely devoid of gradients and artefacts (Chapin et al. 2013).

To account for the attenuation of the signal as a result of the time series filtering, the pipeline re-makes each map with a fake 10 Jy Gaussian added to the raw data, but offset from the nominal map center by 30'' to avoid contamination with any detected source. The amplitude of the Gaussian in the output map gives the signal attenuation, and this correction is applied along with the flux conversion factor derived from the calibrator observations. The final images were made by coadding the 10 maps using inverse-variance weighting, re-gridded with 1 arcsec pixels at both wavelengths. The final images at both wavelengths have been smoothed with a 7'' FWHM Gaussian to improve the signal-to-noise ratio. The FWHMs of the primary beam are 79 and 130 at 450 μm and 850 μm, respectively.

We use a simple model to interpret the Herschel, SCUBA-2, and ALMA observations. Preliminary tests find that the Herschel 70 μm image is significantly more extended than the parent belt imaged with ALMA at 1.3 mm, so the basic requirement of the model is that it reconciles the different extent of these images and produces a blackbody-like spectrum. As is already known, the solution is that the greater radial extent at 70 μm originates in the halo of small grains also seen in scattered light, and that these are not seen with ALMA because these grains emit inefficiently at wavelengths significantly larger than their physical size. Thus, our model comprises two components; the first is the parent bodies for which we use the ALMA modeling results of MacGregor et al. (2013), and the second is a halo whose basic properties are to be informed by previous theory work and determined from the modeling.

With only six photometric points in the spectrum, and a only a few beams of resolution in the new Herschel and JCMT images, our modeling approach is physically motivated but simple. It has been used previously to model many Herschel-resolved debris disks (e.g., Kennedy et al. 2012a, 2012b; Matthews et al. 2014a). For each disk component, a single azimuthally symmetric 3D dust distribution is used, which is simply a small-scale height disk with a power-law surface density dependence between the inner and outer radii. At each radial location, the disk emission is assumed to arise from a modified blackbody, with a power-law radial temperature dependence. This approach is therefore largely empirical, and the derived blackbody parameters can be subsequently compared with more detailed dust models to draw conclusions about the dust properties, in particular, their typical size. The limited resolution also precludes derivation of comparative radial profiles, as was done for HR 8799's halo (Matthews et al. 2014a).

Practically, we generate high-resolution images at each observed wavelength, as viewed from a specific direction to produce the desired disk geometry. We can include an arbitrary number of disk components, as well as point sources to model unresolved sources such as warm inner disk regions or background sources. Each image is then convolved with the appropriate PSF, which may be an observation of a bright calibration source (i.e., PACS, SPIRE), or simply a Gaussian with specified width(s) and position angle (i.e., JCMT, ALMA), and then resampled to the resolution of the observation. We model the ALMA image and not the visibilities, which is acceptable given the good uv coverage and dynamic range of the ALMA data. The models are fit to the data using a combination of by-eye variation of parameters and least-squares minimization. Since we do not search all parameter space, our best-fit models are not necessarily unique, but must be considered a reasonable interpretation.

Our model is guided by previous modeling work. The parent belt is modeled using the structure derived by MacGregor et al. (2013) based on ALMA data, here using a surface density distribution between 8.8 and 40 AU with a radial power law dependence of ${\rm{\Sigma }}\propto {r}^{2.8}$. The inner edge is not tightly constrained. MacGregor et al. find that the inner edge of the parent belt could be as far out as 21 AU, suggesting a much narrower disk width, though we note that Schneider et al. (2014) find evidence of a stello-symmetric warp in the disk with an outer edge of 15 AU, which supports the idea that some component of the disk has an inner edge closer than 21 AU. Note that we are calling this component the parent belt for convenience, and that the extent inward of 40 AU may be due to stellar wind drag from a true source region that is much narrower, but likely still concentrated at ∼40 AU. The emission properties of grains in this belt are assumed to be blackbody-like, and follow ${T}_{\mathrm{PB}}={T}_{0,\mathrm{PB}}/\sqrt{r}$, meaning that the surface brightness profile matches the power-law found by MacGregor et al. (2013). Because radiation and wind forces are relatively weak for AU Mic, small dust can reside in the planetesimal belt and ${T}_{0,\mathrm{PB}}$ can be higher than expected for blackbodies, so we leave it as a free parameter. We also add the unresolved point source component at the stellar position seen in the ALMA data (MacGregor et al. 2013), though this source is too faint to be detected at any other wavelength (see Section 3.2).

To reproduce the Herschel 70 μm image, we include the halo component, with a surface density profile fixed to ${\rm{\Sigma }}\propto {r}^{-1.5}$ (e.g., Strubbe & Chiang 2006), which extends from 40 AU to 140 AU. (The outer edge is poorly constrained by the data and fixed at 140 AU.) The surface density of this component is forced to join smoothly to the planetesimal belt, and the grain properties are allowed to vary via their temperatures and emission spectra. Their spectra are modified blackbodies (i.e., a Planck function multiplied by ${({\lambda }_{0}/\lambda )}^{{\beta }_{\lambda }}$ for $\lambda \gt {\lambda }_{0}$), but we reduce the number of parameters by enforcing a temperature law suited for small grains, ${T}_{\mathrm{halo}}={T}_{0,\mathrm{halo}}{r}^{-1/3}$. We vary ${\lambda }_{0}$, ${\beta }_{\lambda }$, and ${T}_{0,\mathrm{halo}}$ as free parameters. As noted above, we expect ${\lambda }_{0}\lesssim 100$ μm and ${\beta }_{\lambda }\gt 0$, so that the halo is detected at 70 μm, but not at 1.3 mm with ALMA. These values are consistent with measured values of ${\lambda }_{0}$ and β for other disks (i.e., Booth et al. 2013; Matthews et al. 2014a).

Finally, we include the point source visible to the west of AU Mic in Figure 2. This source is included in part to ensure it does not bias the model at 160 μm, but is also used to derive a flux density to subtract from the 160 μm aperture flux derived above. We also include the central point source seen in the 1.3 mm ALMA image. With this simple model, we have four main parameters; one for the parent belt temperature and three for the halo component temperature and spectrum. There are in addition eight R.A./decl. offsets for the Herschel, JCMT, and ALMA images, but these are relatively well constrained. The aim of the model fitting is therefore essentially to find the relative weight of the two components at each wavelength, thereby empirically deriving a coarse spectrum for each in a way that is independent of assumptions about grain properties. We found that the relatively poor resolution (excepting the ALMA data) means that the parameters are poorly constrained due to degeneracies, in particular for the halo component. However, all solutions we found are sufficiently similar that the ultimate interpretation is the same.

Figure 3 shows the original images, the parent belt and halo model components, and the model-subtracted residuals (for just the planetesimal belt and the full model) at wavelengths of 70, 160, 450, and 1300 μm. First, to make it clear that the halo is detected, the fourth column shows the data with only the parent belt subtracted. In this column, the brightness of the parent belt is 1.8 and 1.4 times brighter than in the final model at 70 and 160 μm. The 70 μm residuals clearly show that the disk is more extended than the parent belt, and that a component with greater extent (i.e., the halo) is required. The fifth (rightmost) column shows the residuals with the addition of this component, and that our parent belt + halo model reproduces the data well. Thus, the halo contributes significant emission at 70 μm, but little at longer wavelengths.

At all wavelengths, few residuals above the 3σ level are seen. The negative residuals at 160 μm appear to lie along the scan directions (at roughly $\pm 45^\circ$ relative to N) suggesting that they are artefacts of the imaging. The positive residual to the SE lies beyond the disk extent and is not seen on the opposite side of the star, so should not be considered as a deficiency of the model. The ALMA residual map has a bright point source offset from the stellar position. This is a background source identified by MacGregor et al. (2013) as well. It is located $0\buildrel{\prime\prime}\over{.} 3$ west and $1\buildrel{\prime\prime}\over{.} 8$ south of the star and is detected with a flux density of 0.15 mJy (a 5σ detection).

The main derived model parameters are ${T}_{0,\mathrm{PB}}=245$ K, ${T}_{0,\mathrm{halo}}=161$ K, ${\lambda }_{0}=12$ μm, and ${\beta }_{\lambda }=1$. Because the structure of the parent belt is fixed, ${T}_{0,\mathrm{PB}}$ is moderately well constrained by the contribution required at 70 μm. The halo properties, which are our primary interest, are less well constrained because is it only strongly detected at 70 μm. However, the conclusion that ${\lambda }_{0}$ is shorter than ∼70 μm, and that ${\beta }_{\lambda }\gtrsim 1$ is robust because the halo is required by the images to have little contribution at 160 μm and beyond. Given the conclusions of previous studies that the halo is populated by small dust on eccentric orbits these parameters are consistent with our expectations. The specific value of ${\lambda }_{0}=12$ μm suggests that ∼1 μm sized grains dominate the emission, but a poor constraint on this parameter means that the grains could be larger or smaller, though by no more than an order of magnitude. We return to the likely size of blowout grains when considering small-grain dynamics in Section 4.2.

Figure 7 shows the star-subtracted SED and the disk model, and here the relative contributions of the two components as a function of wavelength can again be seen. We note that the ALMA fluxes (1300 μm) lie somewhat above the fit, but are consistent with the model at the $2-\sigma$ level. We discuss the spectral index of the disk in Section 4.3 below. As noted above, at 70 μm the emission from the planetesimal belt and halo components is similar, but at longer wavelengths the planetesimal belt dominates. The planetesimal belt is dominated by emission at 40 AU, where the grains have a temperature of 39 K. The temperature of the halo component at that separation from the star is slightly higher (50 K), but poorly constrained due to the lack of resolution and degeneracy with ${\lambda }_{0}$ and ${\beta }_{\lambda }$. The fact that these two components yield a composite spectrum that is a simple blackbody is a reminder of the power of resolved imaging, since much information can be hidden within a single component SED.

MacGregor et al. (2013) found evidence of an excess above the photosphere at the position of the star in their high resolution ALMA imaging, and we include this component in the ALMA model image. The residuals from fitting our dust model to the Herschel and JCMT images in Figure 6 show no sign of any other unresolved excess above the stellar photosphere, so the best we can do is place an upper limit on the emission from such a centralized component of the point-source sensitivity at a given wavelength (i.e., $3-\sigma$, MacGregor et al. 2013). It is also possible, given the considerably lower resolution of our observations, that such emission, if present, has been incorporated into the belt emission and removed, which adds uncertainty to the flux density of the belt.

Zoom In Zoom Out Reset image size

From Figure 4 of Cranmer et al. (2013), the coronal model of AU Mic predicts a flux of a few μJy around 100 μm with a relatively flat distribution shortward of 1300 μm where the ALMA detection was made. Our point source sensitivity is on the order of 1 mJy at 70 μm and ∼10 mJy at 160 μm, meaning that, even in the absence of the disk, we would not be able to detect a coronal thermal heating contribution to the total flux in the Herschel data.

A direct estimate of the mass of dust in AU Mic's debris disk can be made from its submillimeter flux densities. Debris disks are optically thin at these wavelengths, and so the mass of the disk is directly proportional to the emission, following the relation:

$$\begin{eqnarray}&&{M}_{d}=\displaystyle \frac{{F}_{\nu }\ {d}^{2}}{{\kappa }_{\nu }\ {B}_{\nu }({T}_{d})}\end{eqnarray} \tag{ 1 }$$

In the Rayleigh–Jeans limit, Equation (1) reduces to

$$\begin{eqnarray}&&{M}_{d}[{M}_{\oplus }]=5.8\times {10}^{-10}\displaystyle \frac{{F}_{\nu }[\mathrm{mJy}]\ {(d[\mathrm{pc}])}^{2}\ {(\lambda [\mu {\rm{m}}])}^{2}}{{\kappa }_{\nu }[{\mathrm{cm}}^{2}\;{{\rm{g}}}^{-1}]\ {T}_{d}[{\rm{K}}]}.\end{eqnarray} \tag{ 2 }$$

While the AU Mic SED supports a single temperature fit of 53 K, our image modeling requires a separate planetesimal belt of 39 K and a warmer halo, which contributes negligible flux density to the total emission at 850 μm (see Figure 7). This temperature is comparable to the 40 K adopted by Liu et al. (2004) in the absence of an SED temperature fit. The color corrected flux density measured with SCUBA-2 at 850 μm (12.5 mJy) is within the $1-\sigma$ limit of the value measured by Liu et al. (2004) with SCUBA (14.4 ± 1.8 mJy), although the SCUBA measurement relied on the single bolometer "photometry" mode, rather than a mapping technique, which could have missed some emission.

Zoom In Zoom Out Reset image size

Taking 12.5 mJy as the 850 μm flux density of the disk yields a slightly lower 1 mm grain dust mass of 0.01 ${M}_{\oplus }$ relative to the Liu et al. (2004) measurement, for the same dust opacity of 1.7 cm² g⁻¹. The dominant sources of uncertainty in the mass are the flux density and the temperature of the planetesimal component of the disk, amounting to a 20% uncertainty in the disk mass.

The SED fit to the disk around AU Mic requires only a single temperature component. The emission is well fit (${\chi }^{2}$ is minimized) from mid-IR through submillimeter wavelengths by a pure blackbody with a temperature of 53 ± 2 K. Figure 7 shows, however, that the halo and planetesimal belt components could still be segregated in temperature, with the parent bodies of the belt colder (∼39 K) than the halo (see also Fitzgerald et al. 2007).

Figure 8 shows the distribution of temperature with grain size for three grain compositions at the location of the planetesimal belt. The temperature distributions were calculated as in Augereau et al. (1999) and Wyatt & Dent (2002). We used three different compositions: (1) a mix of one part amorphous astronomical silicates to two parts organics, (2) pure astronomical silicates, and (3) pure ice (Li & Greenberg 1997). The models shown in Figure 8 are not porous. Porous grains do not have the peak in temperature near 1 μm, instead decreasing steadily to 20–25 K from maximum temperatures similar to the non-porous case. It is clear that pure water ice grains are ruled out at 39 K for all grain sizes, although we note that other compositional mixes could be modeled to satisfy the temperature distribution we have derived.

Zoom In Zoom Out Reset image size

Orbiting dust is subject to forces from both radiation and stellar winds. Around early-type stars, radiation dominates, while stellar winds are thought to be more important for late-type stars such as AU Mic. These forces are commonly split into radial and tangential components. The radial component is called "radiation pressure" or the "radiation force," and effectively reduces the stellar mass seen by the particle, with an analogous effect for the stellar wind. The tangential component is called Poynting-Robertson drag, and causes dust grains to spiral into the star, again with an analogous (and much stronger) effect for stellar wind. See Burns et al. (1979) for a detailed review.

The small grain halo seen previously in scattered light and confirmed here in the far-IR is similar in morphology to those seen around other debris disks (e.g., β Pic: Pantin et al. 1997; Augereau et al. 2001; HR 8799: Su et al. 2009; Matthews et al. 2014a). In other systems, the halo is attributed to the effect of the radial component of the radiation force from the star, which increases the eccentricity of smallest bound grains. The key measure of this effect is the parameter β, the ratio of the radiation force to the gravitational attraction of the star (distinct from the millimeter-wave spectral slope ${\beta }_{\lambda }$ used above). Grains with $\beta \gt 1$ are unbound, and grains liberated from a parent body on a circular orbit are unbound if $\beta \gt 0.5$. Figure 9 shows β as a function of grain diameter, calculated for amorphous silicates + organics, and ice as described above. The solid curves show β for the effect of radiation force on solid grains (left panel) and porous grains (right panel). It is clear that a range of β is possible due to the poorly constrained grain compositions. Also, the assumption of Strubbe & Chiang (2006) that the radiation force blowout grain size exists (i.e., that $\beta \gt 0.5$ for any size) is not well founded, particularly if the grains are porous. Essentially, the low luminosity of an M-type star does not produce a strong enough radiation force to significantly affect the small grains. Even taking into account the flaring of the star is not enough to remove grains via the radiation force (Augereau & Beust 2006).

Zoom In Zoom Out Reset image size

Unlike early-type stars, an M star can have a relatively strong stellar wind (Plavchan et al. 2005). The force from this particle wind causes an effect similar to that of the radiation force and is able to reproduce the scattered light halo around AU Mic (Augereau & Beust 2006). The dotted lines in Figure 9 demonstrate the effective ${\beta }_{\mathrm{SW}}$ parameter for the stellar wind force (using the prescription of Strubbe & Chiang 2006). The left panel shows that for compact grains the wind blowout size is very small, <0.1 μm for the favored AU Mic mass-loss rate of $10\times {\dot{M}}_{\odot }$ (Cully et al. 1994, based on EUV flaring). Such small grains are extremely poor scatterers of light, even at optical wavelengths, so they would appear very faint. Augereau & Beust (2006) get around this problem by including flares, which increase the time-averaged mass-loss rate to $300\times {\dot{M}}_{\odot }$. Such a high mass-loss rate would, however, make the system transport-dominated rather than collision-dominated, i.e., stellar wind drag would also become important, filling the region interior to 30 AU with small dust grains in a manner contrary to the scattered light observations (Strubbe & Chiang 2006). Lim & White (1996) also point out that radio flares would be obscured if the stellar winds of M stars were many orders of magnitude more massive than the Sun's.

The right panel of Figure 9 shows that porosity can also play an important role. By increasing the grain porosity, the effect of wind force increases, in turn increasing the blowout grain size. For p = 0.9, the blowout size is an order of magnitude larger for the same mass-loss rate, and the scattering efficiency of the grains greater, thus providing the likely solution. For a stellar wind rate of ∼100 ${\dot{M}}_{\odot }$, the size of halo grains is roughly 0.1–1 μm. These grains have emission properties consistent with those derived from the image modeling above. This conclusion of high porosity compares well with Graham et al. (2007), Fitzgerald et al. (2007), and Shen et al. (2009) who all find that porosity is necessary to explain the scattering and polarization properties of the grains, although Shen et al. (2009) find that a composition of random aggregates require a porosity of just ∼0.6 compared to the 0.9–0.94 required by Graham et al. (2007) for the Mie theory applied to spherical grains or aggregates.

The Strubbe & Chiang (2006) model of the disk as a collision-dominated, narrow birth ring from which the smallest grains are blown out does a good job of explaining the halo and lack of small grains interior to 30 AU. The recent ALMA results of MacGregor et al. (2013), however, show that there are indeed larger grains interior to the presumed birth ring. Given the 1300 μm observing wavelength, those observations should be dominated by grains with $a\sim \lambda /2\pi$, or 200 μm. Drag forces, which can act on larger grains over long timescales, are therefore a likely cause of this interior emission. Augereau & Beust (2006) concluded that Poynting–Robertson drag is negligible for AU Mic and that stellar wind drag may evacuate the inner regions (see also Strubbe & Chiang 2006). However, more recent work suggests that stellar wind drag can fill the interior regions to some degree, without seriously violating the scattered light constraints (Schüppler et al. 2015). Alternatively, the radially increasing surface density of this disk is reminiscent of a "self-stirred" disk that collisionally grinds down from the inside out (Kenyon & Bromley 2002; Kennedy & Wyatt 2010). Schüppler et al. (2015) also find that the self-stirred model is plausible for AU Mic's disk.

The data compiled in this work provide a long lever arm to measure the spectral index, α where ${F}_{\nu }\propto {\lambda }^{-\alpha }$ is the emission from submillimeter through millimeter wavelengths. As noted above, the 1300 μm fluxes from the SMA and ALMA both lie above the nominal "best fit" line of the models produced by fits to the SED and the images. Comparisons of the 350/1300 spectral index and the 450/1300 spectral index yields α values ranging from ∼1.5 to close to 2.0, which was the inferred spectral index of 350 versus 1300 (SMA) by Wilner et al. (2012). It is clear from the SED that these wavelengths are in the Rayleigh–Jeans part of the spectrum, which was not definite based on the data available to Wilner et al. (2012) and Gáspár et al. (2012).

A spectral index of ∼2 is lower than most other debris disks. Based on the compiled data of Gáspár et al. (2012), the three A stars in that study (β Pictoris, Vega, and Fomalhaut) have the highest value of α, which could suggest that the size distribution is different around A stars, being more consistent with a collisional cascade. Gáspár et al. (2012) show that for a collisional quasi steady-state at a single temperature, the Rayleigh–Jeans tail of the SED should yield a spectral index of 2.65, for a self-similar collisional differential size distribution index of 3.65. The shallower slopes measured by Wilner et al. (2012) and a comparison of the 450 μm flux density from SCUBA-2 of Table 2 (which actually seems to lie somewhat below the SED fit) and the elevated flux densities measured by the SMA and ALMA at 1300 μm suggest a size distribution index of 3.0, i.e., that the size distribution is not the product of a collisional cascade. The shallow slope may however arise from the way that the parent belt and halo components add to yield the total spectrum (Figure 7, Fitzgerald et al. 2007). Such a shallow distribution constraint has also been noted for two IR bright disks (HD 377 and HD 104680) not detected with the VLA at 9 mm (Greaves et al. 2012). In AU Mic, it appears the grains may be very large ($\gt \gt 1$ mm). Longer wavelength observations of the disk beyond 1300 μm would be very beneficial to further constrain the size distribution.

The spectral slope is also dependent on composition. Schüppler et al. (2015) present more detailed modeling of the AU Mic disk, including investigation of the impact of compositional dependence on the derived spectral slope. They find a best fit for a combination of silicate, carbon, and vacuum in equal measure, with little change in the spectral slope fit with small additions of ice or variations of these materials. Their test also provides clear evidence for porous grains. As in our analysis, their SED fits also underestimate the 1300 μm flux densities from SMA and ALMA, though, as in our case, this result should be interpreted with caution since the exploration of the parameter space is limited.