© ESO, 2012

1. Introduction

The evolution of planets in the solar system is intimately connected to the existence of a reservoir of planetesimals in its outer regions. The depletion of the Kuiper belt and the main asteroid belt that occurred during the late heavy bombardment not only delivered large amounts of water, volatiles and carbonaceous material onto the inner planets, but it also reduced the later rate of catastrophic impacts onto the Earth, opening the path to the emergence of life (Morbidelli et al. 2000; Greaves & Wyatt 2010; Hartogh et al. 2011). Although detecting Kuiper belt-like objects around nearby planetary systems will remain an unachievable goal for the foreseeable future, their collisional erosion produces circumstellar disks of dust responsible for a characteristic excess emission detectable at infrared and sub-millimeter wavelengths. Recent studies using the Herschel Space Observatory reveal that as much as  ~30% of nearby Main-Sequence F, G, K stars are surrounded by cold debris rings analogous to the Kuiper belt (latest Herschel/DUNES open time key program results, Eiroa et al., in prep.). However, the cold temperature and the large grain sizes make it impossible to detect solid-state features in their spectrum, thus preventing unambiguous identification of the dust composition. We are thus left with model-dependent methods to constrain the dust composition from color and polarimetric measurements, and from the spectral energy distribution (hereafter SED).

Another uncertainty in planet formation theories lies in the lack of detailed observational constraints for the mutual dust and gas dissipation timescales in circumstellar disks, or more precisely, on how the gas-to-dust mass ratio evolves with time and location in disks at the very early stages of planet formation. The study of stellar clusters of different ages shows that, statistically, a population of young stars loses its massive dust disks in only a few million years on route to the main sequence (Carpenter et al. 2005), although, for individual objects, the transition from an optically thick to an optically thin disk is expected to occur more quickly (a million years or less, Cieza et al. 2007; Currie & Sicilia-Aguilar 2011). However not much is known about the characteristic evolutionary timescales for the gas in these disks. Observing young stars at different evolutionary stages can help elucidate the extent to which gas and dust simultaneously dissipate in disks. This requires observations of both the dust and gas components and, at the same time, modeling of the continuum and line emission.

As part of the GASPS (Gas in Protoplanetary Systems) Herschel Open Time Key Programme, we observed HD 181327, a young (~12 Myr) F5.5V star located at 51.8 pc (Holmberg et al. 2009, Table 1) belonging to the β Pictoris moving group (Zuckerman & Song 2004; Mamajek et al. 2004). HD 181327 was identified as a debris disk hosting main sequence star through its SED with a fractional infrared luminosity L_IR/L_⋆ = 0.2% (Mannings & Barlow 1998). Schneider et al. (2006), from NICMOS coronagraphic observations in scattered light, discovered a ring-like disk of circumstellar dust located at radius  ~89 AU from its star. The ring is inclined by 31.7 ± 1.6° from face-on and shows an apparent brightness assymetry with respect to the minor axis that is well explained by a strong directionally preferential scattering. No photocenter nor pericenter offset is seen in the ring relative to the position of the central star and the disk appears axisymmetric. Recent reprocessing of these data (Schneider et al., in prep.) confirms the annulus shape of the disk, and suggests a narrow radial distribution of the dust in the system (FWHM = 24.5 AU). A low surface brightness diffuse halo is seen in the NICMOS image and in complementary HST/ACS images at a distance of  ~4′′ ≤ r ≤ 9′′ from the central star, which may correspond to a population of very weakly bound grains originating in the main ring. Mid-infrared images at 18.3 μm by Chen et al. (2008) suggest azimuthal asymmetries in the density profile of the ring with respect to the apparent minor axis. These 18.3 μm observations also show that little to no emission (≲ 15%) comes from inside the ring resolved in scattered light. Therefore, the density profile of the HD 181327 dust disk should be well constrained by the scattered light and thermal emission images, allowing a detailed study of the dust properties. Interestingly, from IRS and MIPS-SED data, Chen et al. (2008) also inferred the presence of crystalline ice in the debris disk.

In this paper we perform SED and line fitting of new Herschel/PACS observations of HD 181327 (presented in Sect. 2), supplemented by data from the literature. We also present new ATCA imaging of the disk (Sect. 2). In Sect. 3, we derive surface density profiles from the newly reduced HST/NICMOS data, which well constrain the ring geometry. We describe our dust model and present the results of the SED fitting in Sects. 3 and 4. We assess the gas content in this young debris disk in Sect. 5. Finally we discuss the results in Sects. 6–8 and present our conclusions in Sect. 9.

2. Observations and data reduction

2.1. Herschel/PACS continuum

We obtained new far-IR photometry and spectroscopy of HD 181327 using the PACS instrument onboard the Herschel Space Observatory (Pilbratt et al. 2010; Poglitsch et al. 2010). The observations were carried out in two modes: 1) 159 s of point-source chop-nod mode imaging at 70 and 160 μm (obsid: 1342183658) and 2) two 276 s scan map imaging observations at 100 and 160 μm (obsids: 1342209057-8). The scan maps were executed with the medium scan speed (20″/s), at scan angles of 70° and 110°. The two cross scan maps at different scan angles were combined to increase the signal-to-noise by cutting down on streaking in the image background along the scan direction.

We reduced the photometry with the Herschel Interactive Processing Environment (HIPE, Ott 2010) version 4.2. The uncertainties in the absolute flux calibration for this version are 10% at 70  and 100 μm and 20% at 160 μm (see the PACS Scan Map Astronomical Observation Templates release note¹). For the 70′′ and 160 μm chop-node images, we performed aperture photometry using apertures of 7″ and 11″, respectively. For the 100 and 160 μm scan map observations, we used apertures of radii 12″ and 24″, respectively. An annulus for sky background estimation was placed from 10″−20″ beyond the aperture in each case, and an aperture correction was applied to the measured fluxes. Both estimations at 160 μm are fully compatible within the error bars; the 160 μm scan-map flux is used in the rest of the study. Since the source is marginally resolved at 70 and 100 μm (see Sect. 2.1), we calculated our own aperture corrections following the procedure described in Poglitsch et al. (2010). The photometry results are presented in Table 2.

Fig. 1 PACS 70 μm imaging of (left) HD 181327 and (right) a PSF reference star, HD 148387. Both images were acquired in chop-nod mode. The field of view of each image is 25″ × 25″. Brightness contours at 5, 10, 25, 50 and 75 times the background noise in the HD 181327 image are overlaid on both images.

HD 181327 is more extended than the Herschel point-spread function (PSF) in both the 70 μm and 100 μm images. The 70 μm image of HD 181327 appears in Fig. 1, along with the 70 μm image of a PSF reference star (HD 148387) observed using the same chop-nod mode. The upper panel of Fig. 2 shows the azimuthally-averaged radial brightness profiles for HD 181327 and HD 148387. In the case of the 100 μm scan map observation of HD 181327, we used a different PSF reference observation, i.e. a scan map of α Boo. The 100 μm α Boo radial profile appears in Fig. 2, bottom panel. At both wavelengths, the azimuthally averaged radial profile of HD 181327 is clearly extended beyond that of the PSF reference. The FWHM of the HD 181327 70 μm profile is $\mathrm{6}\stackrel{\mathrm{″}}{.}\mathrm{44}$, compared to a FWHM of $\mathrm{5}\stackrel{\mathrm{″}}{.}\mathrm{28}$ for HD 148387. At 100 μm, the FWHM for HD 181327 is $\mathrm{7}\stackrel{\mathrm{″}}{.}\mathrm{58}$, compared to $\mathrm{6}\stackrel{\mathrm{″}}{.}\mathrm{94}$ for α Boo.

Fig. 2 Radial profiles of HD 181327 at 70 μm (top), compared to HD 148387 as a PSF reference, and at 100 μm (bottom), compared to α Boo as a PSF reference. Both plots show the azimuthally averaged, normalized radial brightness profile vs. radius from the star. Each error bar is the standard deviation of the brightness values at that radius. The FWHM values for the profiles are given on the plot.

2.2. ATCA 3.2 mm observations

The Australia Telescope Compact Array (ATCA) was used to observe HD 181327 at 3.2 mm in May and September 2010. The observations were carried out in continuum double sideband mode, with 2046 × 1 MHz channels per sideband with an effective total bandwidth of just over 2 GHz. The ATCA has five antenna with 3-mm receivers, and our observations were carried out with the compact H75 array configuration with baselines of 31–98 m. The quasars PKS 1815-553 and PKS 1253-055 were used for the gain and bandpass calibrations respectively, and Uranus was used for the absolute flux calibration. The ATCA calibration uncertainty is estimated to be about 20% in the 3-mm band.

Fig. 3 ATCA 93.969 GHz (3.192 mm) map of HD 181327 with contours at −2 (dashed red line), 2 and 3σ (solid blue lines), where σ is the image rms given by 1.1 × 10-4 Jy/beam. The dashed blue area represents the beam. The Gaussian fit is overplot with contours at 2 and 3σ (solid green lines). The red and green crosses denote respectively the star position and the Gaussian fit center.

The data were calibrated and reduced with the MIRIAD package (Sault et al. 1995). The frequency pair for the observations were centered at 93 and 95 GHz, and during the reduction process the data was averaged over 64 channels to improve the signal-to-noise, shifting the frequency pairs to 92.968 and 94.969 GHz. The final CLEANed map combined the two sidebands, resulting in an image at 93.969 GHz (3192.5 μm) presented in Fig. 3. Our total on-source integration time was 770 minutes. We determined an integrated flux of 7.2 × 10^-4 Jy, with an rms of σ = 1.1 × 10^-4 Jy/beam. Including the 20% calibration uncertainty, this yields an uncertainty of 2.5 × 10^-4 Jy. The disk is marginally resolved. Assuming a Gaussian profile, we find a PA of  ~102°, consistent with that seen on the NICMOS image (107 ± 2°, Table 1), and along which the FWHM is  ~7.8″ (beam: 3.5″ × 2.8″, PA = 85.4°). The center of the Gaussian fit is also shifted by 0.6′′ toward the West with respect to the star position². Structures consistent with a belt at  ~100 AU are detected along the semi-minor axis with a peak flux of 4.01 × 10^-4 Jy/beam on the NW side. We note a possible brightness asymmetry with respect to the SE side.

2.3. Herschel/PACS spectroscopy

We obtained far-infrared spectroscopy using the PACS integral field unit which provides a 5  ×  5 array of $\mathrm{9}\stackrel{\mathrm{″}}{.}\mathrm{4}\mathrm{\times }\mathrm{9}\stackrel{\mathrm{″}}{.}\mathrm{4}$ spaxels. The data were taken in chop-nod mode using the line scan and range scan modes. The spectra cover eight small wavelength ranges (~1.8  μm wide) targeting particular emission lines, including [O i] at 63 μm and [CII] at 158 μm. The spectral resolution of the data varies between R ≈ 1100 and 3400. Additional details on the observations appear in Mathews et al. (2010).

The data were calibrated with the Herschel Interactive Processing Environment (HIPE) version 4.2. The wavelength-dependent flux and aperture corrections provided by the PACS development team were applied, giving an uncertainty in the absolute flux calibration of 30% (see the PACS Spectroscopy Performance and Calibration Document²). Significant emission was detected only in the central spaxel. No emission lines are seen, only continuum flux. The central spaxel spectrum centered on the [O i] 63 μm line appears in Fig. 4.

The continua were binned at half the instrumental FWHM and fitted with 1st degree polynomials, using the fluxes and errors produced by our HIPE reduction script. The HIPE uncertainties appeared to be slightly too small, however. New flux random errors were calculated from the rms of the data minus the continuum fits, over wavelength ranges spanning 6 to 14 times the instrumental FWHM, centered on the expected line rest wavelengths. The continuum fluxes at the rest wavelengths of the expected emission lines appear in Table 3. To determine the upper limit on the integrated line flux (F_int) for each available emission line, we integrated the data over a small wavelength range, either  ±1 or  ±1.5 resolution elements, centered on the rest wavelength of the expected emission line. The flux errors were propagated to determine the 1σ upper limit on each integrated line flux. The 3σ integrated line flux upper limits appear in Table 3 as well.

Fig. 4 PACS spectrum of HD181327, centered on the wavelength of the [O i] 63 μm emission line. The linear fit to the continuum flux is overplotted with the red line. The [O i] 63 μm emission line was not detected.

3. Model setup

The dust and gas models developed in Sects. 4 and 5 rely on some assumptions on the spatial distribution of the dust and the gas, as well as on the dust properties. These are detailed below in Sects. 3.1 and 3.4. Additional data collected from the literature and stellar properties are summarized in Sect. 3.2.

3.1. New surface density profile for the HD 181327 ring

Our spatially resolved imaging constrains the location of the dust. We advantageously impose this constraint in order to examine the dust properties in details and to get an upper limit on the gas mass from the Herschel non-detections of gaseous emission lines. The disk around HD 181327 was spatially resolved in scattered light with HST/NICMOS 1.1 μm imaging by Schneider et al. (2006), which alleviates many uncertainties on the disk structure. The NICMOS images revealed a ring-like disk with a peak density position r₀ = 1.71′′ (~88.6 AU, according to the new distance to the star, 51.8 pc, by Holmberg et al. 2009), with an azimuthally medianed photometric FWHM of  ~0.538′′ (~28 AU).

Fig. 5 New HST/NICMOS 1.1 μm surface brightness profiles (top, nominal and  ±1σ profiles), and surface density profiles resulting from direct inversion of the brightness profiles (bottom), assuming two possible values for the asymmetry parameter. The color-coding is the same for both panels. The striped areas identify regions dominated by PSF-substraction residuals where the density has been extrapolated.

Recently, Schneider et al. (in prep.) improved the data reduction of the original HST/NICMOS observations of HD 181327 thanks to PSF subtraction template images better matched to the stellar PSF structure in the HD 181327 images. This results in an overall improvement of the signal-to-noise, especially inside the ring, and leads to a slight sharpening of the belt. The new profiles, displayed in the top panel of Fig. 5, show a peak at a distance of 88.6 AU and a surface brightness fainter by  ~21% (0.76 mJy arcsec^-2 at the peak position) compared to the original Schneider et al. (2006) study.

We derive the surface density profile Σ(r) by inverting the new scattered light surface brightness profile. We employ the methodology described in Augereau & Beust (2006), assuming a disk inclination of i = 31.7° from face-on. The inversion procedure yields the product σ_scaΣ(r), where σ_sca is the mean scattering cross section of the grains at 1.1 μm. This is equivalent to the product of the vertical optical depth τ_⊥(r), multiplied by the mean albedo ω at that wavelength (σ_scaΣ(r) = ωτ_⊥(r), see Eqs. (2), (3) in Schneider et al. 2006). The inversion solution depends on both the anisotropic scattering phase function, and on the disk vertical profile. Nevertheless, Schneider et al. (2006) showed that the density profile is largely independent of the shape for the vertical dust distribution in the range of acceptable vertical scale heights H₀. We here assume a Gaussian vertical profile with a scale height H₀ = 6 AU at the observed peak radius of the ring r₀ = 88.6 AU (H₀/r₀ ~ 0.068, consistent with the theoretical lower limit of 0.04 from Thébault 2009). The observed phase function of the disk at 1.1 μm is well-represented by a Henyey-Greenstein scattering phase function with asymmetry parameter |g_HG|  = 0.30 ± 0.03. This value was discussed by Chen et al. (2008) in light of asymmetries identified in T-ReCS thermal emission images in the Q_a-band. Chen et al. (2008) suggest that the asymmetries may have led to an overestimate of |g_HG|, which may be as small as 0.16 ± 0.04. We thus perform the inversion for both |g_HG|  = 0.16 and |g_HG| = 0.30.

The results are displayed in the bottom panel of Fig. 5, where we also show the impact of the surface brightness uncertainties on the inferred surface density profile. We find an actual surface density peak at 89.5 AU from the star. A higher |g_HG| leads to a globally higher surface density profile. The surface brightness measured interior to  ~1.5′′ (~78 AU) is not reliable because of PSF subtraction residuals. We fit the profiles from 79 to 83 AU obtained for both |g_HG| values by a very steep power-law function of the form  ~r^+10±1. We extrapolate this profile interior to 78 AU as shown in Fig. 5 (striped area). This extrapolation is supported by the mid-IR results of Chen et al. (2008), who show that the inner part of the disk is largely cleared toward the center, containing little to no dust.

Outside of the belt peak position, the improved data reduction results in a profile that falls off as  ~r^−4.7±0.3 from  ~92 to 162 AU, a steeper fall-off than the Schneider et al. (2006) original profiles (~r^-3). The dust is thus more tightly confined to the parent-belt than previously estimated, and the profile is also steeper than theoretical expectations. The exact value of the outer slope can only slightly impact the modelling results, as most of the emission arises from the inner part of the belt where the temperature is larger, except possibly for the longer-wavelengths.

Eventually, the differences between the several radial profiles in Fig. 5 essentially translate into a scaling factor that impacts the total disk mass estimate. We estimate that the uncertainty on the asymmetry parameter adds an additional  ~4% to the uncertainty on our final mass values. The density profile obtained for |g_HG| = 0.30 is used in the rest of this study. We stress that, in order to focus on the other disk properties, we make the assumption that this profile is valid for any dust grain constrained by the SED, i.e. for grain sizes ranging from the micron to the millimetre scale. This assumption assumes no radial segregation in dust sizes and it remains reasonable as long as the flux at any wavelength arises predominantly from regions close to the density peak position (see Sect. 6.1).

3.2. Spectral energy distribution

The Herschel and ATCA data presented in this paper are complemented by IR and sub-mm data collected from the literature. These are listed in Table 4. The shortest wavelength flux we use is the 12 μm IRAS photometry for which the emission shows no excess from the stellar photosphere level. This observation is documented in Schneider et al. (2006) together with IR-excess from the IRAS Faint Source Catalogue and obtained with the Spitzer/MIPS instrument. These are complemented by ground-based mid-IR measurements performed with the T-ReCS imager on Gemini-S by Chen et al. (2008), and 13 points extracted from the 5–35 μm Spitzer/IRS spectrum that we reduced using the c2d Legacy team pipeline (Lahuis et al. 2006).

Additional MIPS observations were obtained in SED-mode by Chen et al. (2006) between 54.5 and 95.5 μm, well-constraining the shape of the disk SED around its maximum emission wavelength. Far-IR excess was reported by Moór et al. (2006) using the ISO database. The shape of the SED in the far-IR domain (100–200 μm) is well constrained by our new Herschel measurements.

Because the uncertainties on the dust optical constants at millimetre wavelengths are high, we chose not to use the ATCA flux to constrain the dust model (see Sect. 6.1 for a discussion on the compatibility of the model to these observations). Therefore, the sub-mm regime is only represented by a measurement with the LABOCA bolometer array at the 12-m telescope APEX (Nilsson et al. 2009, λ ~ 870 μm), which thus represents a key piece of information. These archival data, together with our new observations, constitute the 49 measurements that are displayed in Fig. 6.

Fig. 6 Observed spectral energy distribution of HD 181327 (see Tables 2–4).

3.3. Star properties and dik geometry

We derive excesses from a comparison of the SED with a synthetic stellar spectrum extracted from the “NextGen” model atmosphere grid (Hauschildt et al. 1999) and scaled to match the visible magnitude V = 7.0 mag (Bessell 2000), assuming log (g) = 4.5, T_eff = 6600 K, which leads to a stellar luminosity of 3.33 L_⊙ and a mass of 1.36 M_⊙. This stellar model shows good agreement with the Spitzer/IRS spectrum at wavelengths shorter than 10–15 μm, longward of which the excess starts dominating over the stellar photospheric emission. The star properties and disk geometry adopted in this study are summarized in Table 1.

3.4. Dust properties: size distribution and composition

For both the dust and gas modeling, we make the assumption that the grain population has a single differential size distribution valid anywhere in the ring. This is a central hypothesis that necessarily limits the uniqueness of the solution. We discuss its impact in Sect. 6. The differential size distribution follows a classical power law dn(a) ∝ a^κda, from a minimum size a_min to a maximum size a_max. The minimum grain size a_min and the power law exponent κ are free parameters that are constrained by the SED fitting (Sect. 4). The maximum grain size a_max is a fixed parameter, sufficiently large so that it does not affect the light emission properties at the wavelengths considered (we chose a_max ~ 8   mm). For standard κ values (close to −3.5 or smaller), the mass of solid material in a debris disk can be “arbitrarily” large depending on the assumed higher cut-off value of the size distribution (that can extend up to planetesimal sizes). In this article, we define the dust mass M_dust as the total mass of material contained in grains smaller than 1 mm. This definition has been followed within the GASPS consortium because it allows direct comparison between different studies independently from the model assumptions.

We adopt a complex grain model that bears some similarities to those succesfully employed by several authors for fitting SEDs of young debris disks (Li & Greenberg 1998; Augereau et al. 1999; Li & Lunine 2003a,b), and for fitting scattered light disk colors (Köhler et al. 2008). We consider icy porous aggregates, with silicate cores and organic refractory mantles (or simplified versions of this model; see Sect. 4), to mimic fluffy particles made of essentially unalterred interstellar material (the so-called “cold-coagulation” dust model in Li & Lunine 2003a,b). Assuming cosmic abundances, the volume of the carbonaceous mantle V_C is estimated to be  ~2.12 times that of the silicate core V_Si in the interstellar medium, a value similar to the one derived from the mass spectra of comet Halley in the solar system (V_C/V_Si ~ 2, Greenberg 1998).

To calculate the optical properties, we use a number of representative materials with well established optical indexes. These are either derived from laboratory measurements or from observational and analytical arguments. The materials are: (1) the astronomical silicates (astroSi) from Draine (2003), a theoretical material consisting of amorphous silicates polluted by carbonaceous materials and other metals (ρ = 3.5  g   cm^-3); (2) the ACAR sample of amorphous carbon from Zubko et al. (1996) as organic refractories (ρ = 1.95  g   cm^-3); (3) the H₂O-dominated ice from Li & Greenberg (1998) (amorphous ice, ρ = 1.2  g   cm^-3); and (4) the crystalline water ice from Warren (1984) (ρ = 0.92  g   cm^-3).

The optical indexes of an aggregate are obtained using a standard effective medium theory (EMT), the Bruggeman mixing rule (Bohren & Huffman 1983), and following the methodology of Li & Greenberg (1998) and Augereau et al. (1999). The optical properties are computed with Mie theory valid for hard spheres. The EMT-Mie combination has been shown to be a good approximation for the absorption and scattering efficiencies, provided that grain shapes are close to ellipsoids such that the porosity can effectively be seen as the volume fraction of vacuum (Shen et al. 2009, and references therein). In fact, Shen et al. (2009) show that only the polarization and the phase functions are significantly sensitive to grain shapes. We define the composition of a grain through the volume fractions of its different components: $\begin{array}{ccc}{\mathit{v}}_{\mathrm{Si}}\mathrm{=}\frac{{\mathit{V}}_{\mathrm{Si}}}{{\mathit{V}}_{\mathrm{mat}}}\mathit{,}\hspace{0.17em}{\mathit{v}}_{\mathrm{C}}\mathrm{=}\frac{{\mathit{V}}_{\mathrm{C}}}{{\mathit{V}}_{\mathrm{mat}}}\mathit{,}\hspace{0.17em}{\mathit{v}}_{\mathrm{ice}}\mathrm{=}\frac{{\mathit{V}}_{\mathrm{ice}}}{{\mathit{V}}_{\mathrm{mat}}}\mathit{,}\hspace{0.17em}\mathrm{𝒫}\mathrm{=}\frac{{\mathit{V}}_{\mathrm{vac}}}{{\mathit{V}}_{\mathrm{mat}}\mathrm{+}{\mathit{V}}_{\mathrm{vac}}}& & \end{array}$(1)where V_Si, V_C, V_ice and V_vac are the total volume of astronomical silicates, carbonaceous material, ice and vacuum in the grain, respectively, and V_mat = V_Si + V_C + V_ice  is the total volume of solid material. denotes the porosity. Should one of the materials within the grains reach its sublimation temperature, it is replaced by vacuum and the optical properties are adjusted accordingly in the model. In most cases however, the disk is too cold for any material to sublimate. Figure 7 shows the mass opacity (absorption cross section per unit mass) calculated using Eq. (4) of Draine (2006) for several prototypical dust compositions, with the size distributions given in Table 6.

Fig. 7 Mass opacity of spherical grains computed with the Mie theory and Bruggeman mixing rule for several prototypical compositions detailed in Table 6 and their associated best-fit size distributions.

4. Modeling the dust continuum emission

Fig. 8 Best fit models (Table 6) for several 1- or 2-material mixtures (top) and 3- to 4-material mixtures (bottom). Red crosses: photometric data (Fig. 6). Solid grey line: synthetic stellar spectrum.

We use the Grater code (Augereau et al. 1999) to compute a grid of models with a goal of reproducing the HD 181327 spectral energy distribution. Grater is a SED, image and interferometric data fitter, specially designed to efficiently model optically thin disks with parametric grain size and radial distributions, or distributions from dynamical simulations. The code accounts for both the scattered light and continuum emission of dust grains in thermal equilibrium with a star.

In the present case, the spatial dust distribution is constrained by the resolved HST images (Sect. 3.1, Fig. 5). The surface density profile derived from inverting the observed surface brightness was imposed in the model. Our study of the disk continuum emission focuses on the properties of the dust particles, namely their size distribution and composition. To that end, we computed large grids of models reflecting various possible dust properties, starting from the simplest composition (astronomical silicates) and progressively increasing the model complexity until obtaining a good match to the SED. For each model, the dust mass M_dust and the quality of the fitting are obtained from linear least-square fit between the observed and the synthetic SEDs. The χ² listed in this section are obtained using the 48 data points presented in Tables 2–4, and with the numbers of degrees of freedom listed in Table 6. We subsequently use a statistical (Bayesian) inference method to derive the probability distribution of each parameter (see Appendix A).

Fig. 9 Bayesian probabilities for the models in Table 6: minimum grain size, size distribution power law index, and total dust mass.

In total we consider  ~300 compositions, with 77 × 50 grain size distributions for each, leading to  ~1   150   000 models. All the parameters of the grain models are summarized in Table 5.

4.1. First steps toward a complex dust model

We start by fitting the HD 181327 SED with a single- or two-material grain model, which show the need for a more complex dust model. We consider several mixtures of the materials presented in Sect. 3.4 to assess the broad properties of the grains. The model then includes 3 or 4 free parameters: (1) the minimum grain size a_min; (2) the slope of the size distribution κ; (3) the total mass of the disk M_dust; and in the case of 2-component grains: (4) the volume fraction of either carbon (v_C) or ice (v_ice), or the porosity (). The detailed results of the investigation are summarized in Table 6, the best-fitting SEDs are displayed in Fig. 8 (top panel) and probability curves are shown in Fig. 9.

4.1.1. Silicates or organic refractories?

The most basic grain model to consider is the case of compact spheres made of pure astronomical silicates (astroSi), a configuration that is often used by default in many modeling approaches. We let a_min, κ and M_dust vary and we compute a grid of χ² on which we apply a Bayesian analysis. We find that pure astronomical silicate grains would have a_min = 0.99 ± 0.13   μm and κ = −3.43 ± 0.04, leading to a total dust mass M_dust = 0.35 ± 0.03   M_⊕ for grains smaller than 1 mm (Table 6). These values are consistent with the ones typically found for debris disks, but they actually lead to very poor fits, as illustrated in Fig. 8 (top panel) where the model with the smallest χ² (${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}\mathrm{=}\mathrm{10.3}$) is displayed. Although the model is in fairly good agreement with the data over the full spectral range and agrees well with the sub-mm measurement, it fails to reproduce the shape of the emission peak at  ~50−100 μm  and the mid-IR fluxes, shifting the dominant part of the disk to colder temperatures. It overestimates by 30–40% the excess at far-IR wavelengths (λ ~100−160 μm), in particular the Herschel/PACS measurements.

We take a further step and consider mixtures of amorphous silicates and carbonaceous materials. We ran a Bayesian analysis and find that the actual best composition corresponds to pure carbon. Qualitatively, the effect on the SED is not significant (see Fig. 8, top panel), even if it leads to a better minimum chi-square (${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}\mathrm{=}\mathrm{7.75}$). This model gets closer to the measured fluxes at the peak emission (λ ~ 50−100 μm), but does not provide an overall good fit. The addition of carbonaceous material slightly shifts the best fit parameters of the grain size distribution (Table 6), resulting in slightly larger a_min values and slightly steeper size distributions

4.1.2. Icy grains?

Chen et al. (2008) suggested that HD 181327 may contain a significant fraction of water in the form of crystalline ice, producing a broad peak in the combined Spitzer/IRS (5–35 μm) and MIPS-SED (60–77 μm) spectrum. To assess the presence of water ice in the debris disk of HD 181327, we model the disk with grains made of a mixture of astronomical silicates and water ice, either in the form of amorphous ice or crystalline ice.

The results of the statistical analysis are documented in Table 6. The quality of the fit is only slightly improved compared to the pure silicate model: the best models have ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}\mathrm{=}\mathrm{8.90}$ for amorphous ice, ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}\mathrm{=}\mathrm{10.4}$ for crystalline ice. In fact, the addition of water ice, crystalline or amorphous, does not modify the general shape of the synthetic SED (see Fig. 8, top panel). The smallest χ² calculated with amorphous ice is slightly smaller than the one of crystalline ice, but in both cases the values are too high to accept these models.

4.1.3. Porous grains?

As a final step in our preliminary study of the dust composition, we consider a grid of models consisting of silicate grains with different porosities, from compact silicates to extremely fluffy silicate grains (). The marginal probability distribution for the porosity reveals a bimodal probability distribution: the most probable models are the ones with either no porosity (discussed in Sect. 4.1.1), or very high porosity (). Both solutions are presented in Table 6. A visual inspection of Fig. 8 resolves the ambiguity. The compact grain model better reproduces the minimum and maximum wavelength ends of the SED, and in particular fits the sub-mm point, but it fails at reproducing the intermediate wavelength range around the emission peak. On the other hand, the highly porous grain model is orders of magnitude below the observed value at 870 μm, but it is within the error bars in the other regions of the SED, especially at Herschel wavelengths. In sum, none of these models lead to overall good fits to the whole SED.

4.2. Construction of an advanced grain model

From what precedes, we conclude that none of the 2-component mixtures considered are sufficient to provide a satisfactory fit to the SED. This is essentially a consequence of the disk surface density being known a priori in this model; it cannot be adjusted to compensate for the emission properties of the dust grains. We note, however, that the addition of one component to silicates, should it be carbonaceous material, amorphous ice, or porosity, provides some slight improvements. This suggests that their combined effect could produce a more realistic model.

We evaluate the improvement provided by a 3-component mixture on the fitting of the SED. We showed in Sect. 4.1.1 that it is not clear how to disentangle the respective contribution of silicates and organic materials by looking at the SED. Instead of including an additional free parameter, we consider grains (approximately) in fixed Si/C ratio described in Sect. 3.4 (v_C/v_Si = 2), and we restrict the discussion to the ice fraction and the porosity of the grains. We thus assume that grains are made of silicate and carbonaceous material, which are either icy (v_ice parameter) or porous ( parameter). The results are reported in Table 6. For porous ice-free grains, the best model is found for a porosity % and leads to some improvement to the fit with ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}\mathrm{=}\mathrm{6.05}$. If we now consider only compact icy grains, we can achieve a very reasonable ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}\mathrm{=}\mathrm{4.48}$ with comparison to previous models, assuming an amorphous ice fraction v_ice = 0.7. Overall, both the addition of either water ice or vacuum to silicate + carbonaceous grains do contribute to improving the fit to the SED of HD 181327. This can can be seen in Fig. 8 (bottom panel). Nevertheless, the fits are still unsatisfactory (${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}$ much larger than 1), and justify an even more advanced dust model.

We next consider a 4-component grain model inspired by the “cold coagulation dust model” by Li & Lunine (2003a,b), which has been successfully employed for modeling debris disks (Li & Greenberg 1998; Augereau et al. 1999; Li & Lunine 2003a,b; Köhler et al. 2008). We keep the Si/C volume ratio fixed (v_Si:v_C = 1:2) and let both the ice fraction and the porosity vary, in addition to the other usual free parameters in this study, namely the minimum grain size a_min, the slope of the size distribution κ, and the total mass of the disk M_dust. This leads to a five-free-parameter disk model to fit the 48 data points in the SED. We compute three grids of models corresponding to the three surface density profiles shown in Fig. 5 for |g_HG| = 0.30 (i.e. the nominal and  ± 1σ profiles). The Bayesian analysis is then applied to the 3 grids separately and the probabilities for each point of the parameter space are calculated as the sum of the probabilities for the 3 profiles. The detailed results are documented in Table 6.

Fig. 10 Bayesian probabilities for the grain composition in the case of 4-material dust mixtures.

Figure 10 presents the results of the Bayesian analysis on the composition of the dust grains, with a 2D-map showing the probability as a function of and v_ice for the grid models, after projection of the probabilities over the other parameters. It clearly illustrates the need for a composition within a limited range of possible values. In particular, the fraction of H₂O-dominated ices is well defined: the probability peak is obtained for an ice fraction v_ice = 0.67 ± 0.07. The porosity is not so strongly constrained but we find a peak of probability for . The best fit model in terms of χ² is found for grains made of 23% carbonaceous materials, 12% silicates, 65% ice, and with a porosity . As can be seen in the bottom panel of Fig. 8, for such grains, the synthetic SED very reliably describes the observed SED from the mid-IR to the sub-mm, with ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}\mathrm{=}\mathrm{1.64}$, a value significantly better than those obtained for simpler grain models. The grain size distribution derived from the statistical analysis has a_min = 0.81 ± 0.31 μm and κ = −3.41 ± 0.09, a power law index very consistent with that predicted for theoretical collisional cascades (κ = −3.5). The mass in grains smaller than 1 mm, M_dust = 0.051 ± 0.016   M_⊕, may look surprisingly small at first glance compared to the other models, but it results from the significant fractions of ice and vacuum in the grains, implying larger mass opacities (Fig. 7). The temperature of the grains at the surface density peak position (89.5 AU) ranges from 86 K for the smallest grains to 40 K for  ~8 mm-sized grains.

To summarize, our approach, based on a fixed surface density profile constrained by resolved imagery of the disk, clearly highlights the need to consider the presence of both ice and porosity to find an overall good fit to the SED. Icy grains are naturally expected at such a large distance from a solar-type star and the fluffiness of the grains is likely a by-product of their formation process. We applied a similar cold-coagulation dust model as Li & Lunine (2003a) did fitting the very young HD 141569A debris disk, except that we used no a priori assumption on the ice fraction, but with arguing the grain composition must reflect the primordial abundances of condensable elements, they chose v_ice ≈ 30% for cold regions and they find a best-fitting porosity . Hence we propose a dust model that is especially ice-rich with a degree of porosity that is not unexpected. Further comparisons with the literature and implications on the nature of the dust are discussed in more details in Sect. 6.

5. Limits on the gas content

At 12 Myr-old, HD 181327 offers a chance to assess the amount of gas in a young debris disk, at a time when, according to planet formation theories, the accretion of gas by giant planets has already occurred while the formation of rocky planets might be ongoing.

We use the upper limits on the oxygen and ionized carbon fine-structure lines ([O   I]  at 63 and 145 μm, and  [C   II]  at 158 μm, see Table 3) to constrain the amount of gas present in the disk. We solve the photochemistry and energy balance in the disk using the photochemical code ProDiMo³. The abundance of 71 species is computed in steady-state, under the influence of adsorption, desorption (thermal, photo and cosmic-ray induced), and molecular hydrogen formation on grain surfaces.

The gas density profile is supposed to match the one of the dust, as presented in Sect. 3.1. The total mass of solid (defined as the mass in grains from 0.81 μm to 8 mm) is constrained by the best dust model of Sect. 4.2 (M_solid = 4.9    ×  10^-7   M_⊙ = 0.164   M_⊕) and we let the total atomic+molecular gas mass M_gas vary from a gas-poor disk (M_gas/M_solid = 10^-2) to a gas-rich disk (M_gas/M_solid = 10³). The parameters of the model are summarized in Table 7. One of the main parameters controlling the gas energy balance is the amount of PAHs (f_PAH), which remains unconstrained in the disk of HD 181327 because their abundance is too small to produce detectable PAH emission features in the IR. We adopted two “extreme” values for f_PAH: a low PAH abundance (f_PAH = 10^-5) and a high PAH abundance (f_PAH = 0.1)⁴.

The non-LTE line fluxes for the  [O   I] ,  [C   II] , CO J = 3  →  2, and CO J = 2 →  1 transitions are given in Fig. 11 and Appendix B. Among all the lines, only the CO J = 3 →  2 and J = 2 →  1 line fluxes behave monotonically, increasing in flux with increasing disk mass. The reasons for the apparent erratic behaviour of the fine-structure lines come from the non-linear nature of several physico-chemical phenomena. For example, the  [C   II]  flux depends on the abundance of C⁺ and on the excitation of the first energy level. The C⁺ abundance first increases with total gas mass until CO self-shielding against photodissociation is efficient enough such that a significant fraction of C⁺ is converted into neutral C and CO. The presence of CO also changes the heating and cooling balance, since the CO absorption of IR stellar photons becomes the main source of heating. The change from a disk with no CO self-shielding to a CO self-shielded disk manifests itself as a jump in CO fluxes between a disk with gas-to-solid ratio of 1 to a disk with gas-to-solid ratio of 10. The abundance of PAHs has a weak influence on the gas temperature because the small UV flux from a F5/F6 star without excess accretion luminosity cannot ionize PAHs. On the other extreme, when the gas mass is very small (at gas-to-solid mass ratio below 1), the PAH and molecular absorption heating are nonexistent. Photoelectric heating by large grains dominates and explains why the line fluxes are similar for the high and low PAH abundance models.

Line fluxes from models with f_PAH = 10^-5 are consistent with the PACS upper limits for all gas-to-solid mass ratios (upper part of Table B.1). Only in the case f_PAH = 0.1, the disk gas mass is constrained by [OI] at 63 μm to be below  ~17 M_⊕ (M_gas/M_solid ≲ 10², lower part of Table B.1).

Eventually, even if the integration of the  [O   I]  63 and 145 μm, and  [C   II]  lines we performed was not sufficient to derive strong constraints on the gas-to-solid mass ratio in the disk, the combination of the GRaTeR and ProDiMo codes proved efficient to make predictions on the fluxes in the CO J = 3 →  2 and J = 2 →  1 lines. These lines will be detectable for instance by ALMA, finally allowing to assess the gas content of the debris disk.

Fig. 11 Expected gas line fluxes as a function of gas-to-solid mass ratio for either a low-PAH abundance (thin lines) or a high-PAH abundance (thick lines). The horizontal lines denote the PACS 3σ upper-limits and do not include the 30% flux calibration uncertainty. The gas-to-solid mass ratios assume a solid mass of 0.164 M⊕.

6. Additional constraints on the disk

The surface density profile used in this study comes from a new reduction of the scattered light images of Schneider et al. (2006) obtained with HST/NICMOS (see Sect. 3.1). This ensures by definition that the shape of any synthetic scattered light surface brightness profile in this study perfectly agrees with the observed profile at 1.1   μm. However, the consistency is limited to the shape, and the agreement in terms of absolute flux is not guaranteed. This is because our fitting approach in Sect. 4 consists of adjusting the disk mass to get a best fit to the SED, which is equivalent to scaling up or down the surface density profile. A reason for adopting this strategy comes from the fact that the scattered light images do not directly provide the surface density Σ(r), but its product by the mean scattering cross section σ_sca at the observing wavelength (1.1 μm). The fitting process to the SED, on the other hand, yields an independent, absolute scaling of the surface density profile, thus allowing one to obtain an initial estimate of σ_sca. We find σ_sca = 3.7   μm² in the NICMOS 1.1   μm band.

A second estimate of σ_sca is obtained by calculating the theoretical scattering cross section for compact spherical grains given the grain properties (composition, size distribution) that were inferred in Sect. 4: ${\mathit{\sigma }}_{\mathrm{sca}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{=}{{\mathrm{\int }}^{}}_{\mathit{a}\mathrm{=}{\mathit{a}}_{\min }}^{{\mathit{a}}_{\max }}\mathit{\pi }{\mathit{a}}^{\mathrm{2}}{\mathit{Q}}_{\mathrm{sca}}\mathrm{\left(}\mathit{a,\lambda }\mathrm{\right)}\mathrm{d}\mathit{n}\mathrm{\left(}\mathit{a}\mathrm{\right)}$ where Q_sca is the scattering efficiency at wavelength λ and for a grain radius a. For the best model (in the Bayesian sense, Table 6) we find σ_sca = 16.8   μm² at 1.1   μm, a factor of about 4.5 higher than the first estimate. This implies a peak flux density of 3.14 mJy arcsec^-1 inconsistent with the 0.76 mJy arcsec^-1 measured in the NICMOS profile. In other words, this means that our model predicts an albedo (ω = 59% at 1.1   μm) that is a factor of about 4.5 larger than the observed value. This suggests that smaller grains than those required to fit the SED dominate the scattering behavior of the disk. It is worth noting that such a conclusion is not specific to HD 181327, but is also observed in other debris disks. Recent examples include HD 207129 for which Krist et al. (2010) predict an albedo (based on dust properties also obtained from a fit to the SED) that is an order of magnitude larger than the value suggested by the scattered light image. A similar inconsistency is noticed for HD 92945 by Golimowski et al. (2011) who point out that different observing diagnostics suggest different minimum grain size.

It is important to keep in mind that although our model suggests the grains are very porous and hence likely to have an irregular, fractal structure, we used hard spheres (Mie theory) to calculate both the thermal emission and scattered light optical properties. Unlike the thermal behavior, light scattering will likely strongly depend on the exact shape of the dust aggregates, and on the size of their elemental constituents, two features that are not accounted for by the Mie theory (Voshchinnikov et al. 2007). Another piece of evidence that Mie theory might be a limiting factor for a self-consistent model of HD 181327 comes from the observed scattering asymmetry parameter g_HG that is a measure of the direction of the light scattered by a particle. For our best model (Sect. 4.2) the Mie theory yields a mean scattering asymmetry parameter |g_HG|  = 0.95 at λ = 1.1   μm, much larger than the |g_HG|  = 0.30 measured on the NICMOS image. Only a population of much smaller particles could produce such a low scattering asymmetry parameter if we stick to the Mie formalism, as already noted by Schneider et al. (2006). An analogous problem arises from the modeling of the debris disk of HD 92945 (Golimowski et al. 2011). Our micron-sized grains are highly anisotropic, essentially due to their spherical shapes. Incorporating vacuum homogeneously inside the grains does not suffice to decrease |g_HG| and a more accurate description of the shape of fluffy aggregates is eventually needed. In fact, the scattering efficiencies and asymmetry parameters of fluffy aggregates are known to be correctly described by an EMT-Mie model (porous spheres) only in the Rayleigh regime, i.e. when the size of the inclusion (vacuum here) is small compared to the wavelength (see e.g. Voshchinnikov et al. 2007). This is not the case for the relatively large grain sizes and porosity when the wavelength becomes small.

6.1. A narrower parent belt?

Unlike scattering, we expect the continuum thermal emission to be largely independent of the exact shape of the grains. Our best model of HD 181327’s disk produces a sub-mm flux of 39 mJy at 870 μm, consistent within 2σ with the LABOCA measurement (52 ± 6 mJy). The ATCA measurement at 3.2 mm was not model-fitted, because of large uncertainties on the optical constants at millimetre wavelengths. As a sanity check, we calculate the flux in the ATCA band using a power-law extrapolation of the SED. In the far-IR to mm regime, the model emission declines as λ^{−2.51 ± 0.165}. This results in a flux equal to $\mathrm{1.}{\mathrm{6}}_{-0.3}^{\mathrm{+}\mathrm{0.4}}$ mJy, which is within  ~3   σ compatible with the measured flux of 0.72 ± 0.25 mJy. Thus, the model remains valid for the longest wavelength observations (coldest emitting particles), showing that the ring resolved in the near-IR with the HST is compatible with the spatial distribution of the largest grains, at least at first order.

A limit of our approach is the assumption that the grain size distribution is the same everywhere in the disk, while it is known that only the smallest grains can populate the outer disk due to radiation pressure. The maximum size and power law slope of the dust size distribution should vary with the distance from the star (e.g. Fig. 4 in Augereau et al. 2001), whereas the model keeps both large grains in the outer disk and a constant slope for the size distribution. Such a radial variation of the size distribution is suggested by a comparison of a synthetic image at 70 μm with the PACS 70 μm map (Sect. 2.1). The synthetic image appears broader than the observed one with a FWHM of 7.68″, versus 6.44″ for the PACS image. Our interpretation is that the PACS band reveals larger grains that experience less the force due to the radiation pressure and thus are more bound to the parent belt. Assuming all grains are equally distributed results in putting too many large grains in the outer disk: a more complex model should account for the true dynamics of the particles.

6.2. Emission from inside the belt?

The region of the disk interior to 75 AU cannot be reliably assessed from the NICMOS data because of the coronagraph and artefacts from the PSF-substraction interior to 1.5″ (78 AU). Therefore we cannot exclude the existence of an additional dust component interior to the main belt, that would be invisible in the surface density profile we use. If micron-sized grains populate the inner disk, they would radiate mainly in the mid-IR, which might affect the results of our study. From the Gemini/T-ReCS images of Chen et al. (2008), we estimate that less than  ~15% of the light emitted by the disk at 18.3 μm originates from a region interior to 1″ (Fig. 2 in Chen et al. 2008).

We perform a simple test of the robustness of our results. Suppose that a belt of black body dust particles is orbiting the star at a distance r < 1″. We calculate a black body spectrum than peaks at 18.3 μm (T = 150 K) and we adjust it such that its maximum flux represents 15% of the 53 mJy attributed to the disk at that wavelength. Then we subtract the black body flux from the observed SED of Table 4 to simulate an observation of the system excluding the hypothetical inner component. We reproduce the approach of Sect. 4.2 with this new SED. Interestingly the results of the Bayesian analysis are not significantly impacted. To account for the deficit of flux at mid-IR wavelengths, only the minimum grain size needs to be slightly adjusted to eliminate the smallest grains: a_min = 1.0   μm and κ = −3.43, while we have a_min = 0.89   μm and κ = −3.42 for our best fit (Table 6). The composition remains unchanged for the best model (${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}\mathrm{=}\mathrm{1.52}$) and no change is induced that would improve the consistency of the model with the scattered light flux (Sect. 6). Consequently we cannot reject the possible existence of an inner component to the disk, but none of our conclusions would suffer from its existence.

Fig. 12 βpr = |FP |/| FG| as a function of grain radius. The horizontal dashed line denotes the limit between bound (β < 1/2) and unbound (β > 1/2) grains.

7. Grain dynamics and ice origin

7.1. Blowout size

Grains in the HD 181327 debris disk for which radiation pressure overcomes gravity are put on unbound orbits and eventually ejected from the system on a short time scale (see for instance Krivov 2010, and ref. therein). The balance between radiation pressure (F_P) and stellar gravity (F_G) forces is evaluated with the distance independent, but size-dependent, β_pr factor, given by $\begin{array}{ccc}{\mathit{\beta }}_{\mathrm{pr}}\mathrm{=}\frac{\mathrm{\left|}{\mathit{F}}_{\mathrm{P}}\mathrm{\right|}}{\mathrm{\left|}{\mathit{F}}_{\mathrm{G}}\mathrm{\right|}}\mathrm{=}\frac{\mathit{\sigma }}{\mathit{m}}\frac{\mathrm{⟨}{\mathit{Q}}_{\mathrm{pr}}\mathrm{⟩}{\mathit{L}}_{\mathit{\star }}}{\mathrm{4}\mathit{\pi cG}{\mathit{M}}_{\mathit{\star }}}& & \end{array}$(2)where  ⟨Q_pr⟩  is the pressure efficiency averaged over the stellar spectrum, m the grain mass, σ the grain geometrical cross section, c the speed of light, G the gravitational constant, and L_⋆ and M_⋆ documented in Table 1. Assuming that the grains are spherical, Eq. (2) yields a critical size below which the grains originating from the parent belt are ejected from the system: a_blow = a(β_pr = 1/2) (assuming circular orbits for the parent bodies).

In Fig. 12, we represent the β_pr ratios as a function of grain radius for different compositions. Our best 4-component grain model has a blowout size a_blow = 4.83 μm, that is 5 times larger than the inferred minimum grain size in the system. Models implying no porosity nor carbon, on the other hand, have blowout sizes of the order of the inferred value of a_min, more consistent with the theoretical argument that grains below the blowout size should be expelled from the system. The response of silicate grains to the stellar irradiation and gravity fields is only slightly affected by the presence of ice, leading to similar shape and moderated values for β_pr(a). Organic refractory grains on the other hand are highly sensitive to the radiation pressure and they can reach high β_pr values. Figure 12 shows that porosity has the effect of flattening the β_pr(a) function, pushing the critical value β_pr = 1/2 to larger grains. Hence the addition of porosity to a mixture of carbon and silicates produces a strong effect: β_pr not only reaches a high value at its maximum due to carbon, but the β_pr(a) function decreases slowly such that the critical value β_pr = 1/2 is reached only for large grain sizes.

Decreasing the amount of carbon in the grains (it was imposed somehow arbitrarily, see Sect. 4.2) would mechanically decrease their blowout size. Pure silicate grains with the same porosity and ice content have a_blow = 0.89 μm, but they do not provide good fits to the SED. We computed additional grid of models with alternative v_C/v_Si values. The results confirm the conclusion of Sect. 4.1.1: all Si+C+H₂O+porosity models lead to similarly good fits to the SED, as long as a sufficient fraction of carbon is considered. For instance, if silicates represent only 10% of the C+Si mixture, then v_ice remains at 0.65 and the porosity is reduced to 55%, with M_dust = 0.043   M_⊕, a_min = 0.89   μm and ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}\mathrm{=}\mathrm{1.66}$. All these values remain within the uncertainties of the original model. Models involving carbon-poor grains do not provide good fits to the SED, but as long as v_C/v_Si ≳ 0.5, the results remain unchanged. In summary, the discrepancy between a_min and a_blow does not lie on the hypothesis we made on the composition of the grains.

Our estimation of a_blow is based on the assumption that the grains are spherical and homogeneous, but the actual response of inhomogeneous aggregates to radiation pressure cannot be that simple. Saija et al. (2003) for instance show that the fluffiness of the grains tends to decrease their radiation pressure cross-section, as if aggregates behaved like their elemental constituents. This would make grains below the theoretical blowout size less sensitive to radiation pressure, and allow them to remain bound to the star. Thus the apparent inconsistency may be interpreted as another signpost of the complex structure of the grains. We emphasize that fixing a_min as equal to a_blow calculated with the (EMT-)Mie theory would not be a relevant approach when dealing with porous grains.

Another possibility is that grains are not in dynamical equilibrium (i.e. an imbalance between production and ejection), following for instance a major collisional event, such as an asteroid break-up. An excess of dust grains from all sizes would have been produced, including grains below the blowout size. These grains could have survived destructive collisions and radiation pressure blowing for long enough to be observed. Of course, observing such an event is unlikely, because dust typical lifetime is only of the order of 1000 years, which is small with comparison to the age of the star (see Sect. 7.2). But this “rare event hypothesis” would be consistent the scenario proposed by Chen et al. (2008) to justify the brightness asymmetry they claim to observe in thermal light.

7.2. Time scales and origin of the ice

Collisions occurring between grains in a debris disk are a key process as they generally result in the production of smaller dust particles, either by fragmentation of the target, or excavation. Any phenomenon longer than the characteristic collisional time scale will be considered negligible in the evolution of the dust population. Based on the principle that any particle will cross the disk midplane twice during one orbit, Thébault & Augereau (2007) proposed a simplified relation for the collisional time scale, valid for the smallest grains: ${\mathit{t}}_{\mathrm{col}}\mathrm{=}{\mathit{t}}_{\mathrm{orbit}}\mathit{/}{\mathit{\tau }}_{\mathrm{\perp }}^{\mathrm{geo}}$, with t_orbit the orbital period at the surface density peak, and ${\mathit{\tau }}_{\mathrm{\perp }}^{\mathrm{geo}}$ the geometrical vertical optical depth. Here the orbital period at 89.5 AU is  ~850 years, leading to t_col ~ 7600 years for the best model. This is more than three orders of magnitude smaller than the age of the system and two orders of magnitude smaller than the characteristic time-scale for Poynting-Robertson drag: t_pr = 0.1 Myr for β_pr = 1/2 grains (Wyatt 2008). The evolution of the HD 181327 debris disk is clearly dominated by collisions requiring a mechanism to generate dust on a 10³−10⁴ years characteristic time scale.

Attempts to model time-dependent chemistry with ProDiMo reveals that the grain water ice mantle would be photo-evaporated in a few hundred years (for  ~1 μm grains), to a few thousand years (for  ~10 μm grains), a similar (or shorter) time scale to the dust mutual collision, requiring as well an ice reprocessing mechanism. This provides a consistent view that the icy dust grains are maintained through collisions in a reservoir of icy planetesimals. We highlight that, while detecting the mid-infrared spectral features in the emission of a debris disk is very challenging – unambiguous detection have only been made in brighter circumstellar disks (see, e.g. Malfait et al. 1999) – the sole analysis of the continuum SED alone can provide clues on the amount of ice in dust disks, thus on the nature of the parent bodies, often referred to as the leftovers of planet formation.

In the solar system Kuiper belt, the study of icy planetesimals, Kuiper belt objects (KBOs) in particular, is a prolific topic (for a recent review see e.g. Jewitt 2010). The presence of ice, together with the porous nature of the objects, is inferred from estimates of their densities. Unambiguous detection of water ice has been made possible by the spectroscopic near-IR signature of crystalline ice, for instance by Delsanti et al. (2010) who completed their study of the plutino Orcus by radiative transfer considerations concluding that it contains a mixture dominated by amorphous and crystalline water ice. Confronted with space weathering, crystalline ice should be amorphidized in a short timescale. This implies that ice needs to be stored (and possibly produced) in the interior of the object before being brought to the surface, possibly through cryovolcanic events. The nature of water ice in KBOs remains unclear, but the vast majority of these objects are too small to experience such geological processes (Dumas et al. 2011,and references therein). Thus KBOs covered by crystalline ice can only constitute particular examples, and the water ice covering the collisionally produced dust in the solar system is more likely amorphous. The young age of HD 181327 allows the possibility that primordial crystalline ice has survived in its debris disk. Our study tilts in favor of the amorphous nature of that ice, suggesting that either ice was predominantly incorporated in an amorphous form in HD 181327’s Kuiper belt objects, or that the amorphization timescale is shorter than  ~10 Myr in this distant belt.

8. Discussion

The properties of the HD 181327 debris disk can be compared to those of the disks surrounding two of its co-moving “sister” stars, β Pictoris (A3V, 1.7 M_⊙) and AU Microscopii (M1V, 0.5 M_⊙), which formed in the same molecular cloud from the same material and are coeval. Of course comparisons are limited by the specific models used to study each object. The relatively large porosity found here should be put in the context of the grain properties found for AU Mic (Graham et al. 2007; Fitzgerald et al. 2007). Assuming a Si+C+H₂O mixture, only highly porous grains (), reproduces the scattered light profiles (exhibiting strong forward scattering), the high polarization and the thermal SED. In that model, the Mie theory remains valid to represent the scattering properties of porous grains, because the grains filling the outer disk are very small and stay in the Rayleigh regime. A total dust mass of 0.01   M_⊕ is found, that essentially reflects the contribution from the inner mm-sized grains, concentrated in a  ~40 AU belt. A similar porous grain model was used by Li & Greenberg (1998) to model β Pictoris  yielding a dust mass  ~0.4 M_⊕, contained in the  ~100 AU-wide outer disk. Overall, “cold coagulation-like” dust models prove effective at reproducing observations of these young debris disks. A trend appears for increasing dust mass with spectral type as could be intuitively expected, but despite this difference of scale, the three stars are orbited by a similar material.

Given the young age of these objects, we are investigating the conditions holding within disks just after the planetary embryos may have formed, the gas dispersed, and possible accretion onto outer giant planets happened. To explain the belt-like structure of the HD 181327 disk, the most straightforward explanation is to invoke the dynamical influence of a planet inside the belt. The Wisdom (1980) analytical criteria gives the width of the chaotic zone a planet mean-motion resonances would create. Assuming circular orbits for the planetesimals, as suggested by the sharpness of the belt, this yields a range of possible values for the planet semi-major axis depending on its mass. Assuming the inner edge of the disk is located at 78 AU from the star, the yet undiscovered planet would lie from 73 AU for a Neptune-mass planet, down to 59 AU for a 10 Jupiter-mass planet.

The material at that stage appears as very hydrated providing a considerable reservoir of water in the outer system. The existence of such a reservoir is hypothesized in scenarios of the solar system formation to explain the unexpectedly large water content on Earth for instance (Morbidelli et al. 2000). Indeed, terrestrial planets at this age (12 Myr) are likely too hot to hold on to water and later impacts of icy comet-like bodies from the outer regions are a possible mechanism for delivering water and volatiles to planets in the habitable zone. Thus, not only HD 181327’s debris disk reveals that the system has already grown planetary embryos – including maybe the bricks of rocky planets in the inner disk – but it likely holds the conditions that were prevailing before the solar system late heavy bombardment that depleted the Kuiper belt and determined terrestrial planets history.

9. Summary and conclusions

We have presented new Herschel/PACS and ATCA imaging of the HD 181327 debris disk, a 12 Myr old star belonging to the β Pictoris moving group. The 70 and 100 μm images as well as the 3.2 mm observation reveal a marginally resolved extended emission. Fluxes in the three PACS photometric channels were derived and complemented by continuum fluxes obtained from PACS spectroscopic observations. Together, these data provide precise information on the continuum thermal emission of the debris disk, but detailed models such as the one developed in this study are required to interpret the precise nature of the emitting material. Degeneracies between the grain properties and disk’s spatial structure are broken in our model by making use of a reprocessed HST/NICMOS coronagraphic image revealing that the dust is bound to a narrow belt peaking at 89.5 AU from the central star. We focus our study on the dust properties (composition, size distribution) using a radiative transfer code, and run a large grid of debris disk models for subsequent statistical (Bayesian) analysis. We conclude that the thermal light spectrum of the disk is well reproduced assuming a mixture of silicates, carbonaceous material and amorphous ice, as well as vacuum to mimic a significant porosity. The observed SED can be very reliably reproduced by a population of grains larger than  ~1 μm, which follows a power-law distribution matching the ideal collisional equilibrium case and represents  ~0.05 M_⊕ (1.5 × 10^-7   M_⊙) in grains with sizes up to 1 mm. Extending the size distribution to the parent bodies of the observed dust grains leads to a mass of  ~50 M_⊕ for planetesimals up to 1 km. The composition we infer for the dust grains must be a good representation of a considerable reservoir of material hidden in unseen planetesimals.

In this study, we attempted to use Herschel/PACS spectroscopy to detect for the first time far-infrared emission lines in the debris disk. The non-detections we obtain were used to investigate the gas content of the debris disk with our photochemical disk model ProDiMo based on the dust model found with GRaTeR. The main uncertainty lies in the amount of PAHs in the disk. The current non-detections of [OI] and [CII] lines alone do not provide unambiguous upper limits to the gas content. Only in the case of a high PAH abundance, we can set an upper limit to the gas mass of  ~17   M_⊕ (one Neptune mass). Coupling with other tracers, in particular with the CO lines accessible to ALMA, offers much better prospects to reach lower limits on the low gas content of debris disks. However, steady-state chemistry is likely not ideally suited to debris disk modeling, and in future studies more appropriate prescriptions, time-dependent chemistry in particular, will be explored.

We note that both the low albedo estimated from the scattered light images and the dust dynamics are not well modelled by spherical grains filled with vacuum, and that the grains are likely complex fluffy aggregates. Though their thermal behavior is correctly described by the EMT-Mie theory, the scattering of light at short wavelengths and their sensitivity to radiation pressure require a more accurate description of the exact shape of the aggregates. Our statistical analysis highlights that ignoring the porosity of the grains or setting the minimum grain size as equal to the blowout size would result in incorrect interpretation of the observations, in particular because the grain opacity is strongly dependent on the exact grain composition.

We stress some necessary limitations of the models. First, we use a grain size distribution that is independent of the distance to the star, whereas it is known that only the smallest grains can populate the outer parts of the disk under the effect of radiation pressure. In this paper we discussed clues that the observed scattered light profile may not be a perfect representation of that seen in thermal light, by Herschel or ATCA. The long wavelength emission emanating from the largest grains might originate from an even narrower belt. Further studies will need to account for the radial differentiation of grain sizes with the objective to fit altogether resolved images from the near-IR to the millimetre domain. Secondly, the Dohnanyi-like power-law size distribution is known to be an insufficient description; if the real distribution implies wavy patterns as proposed by Thébault & Augereau (2007), then adjusting a power-law likely puts too many/too few small grains in the disk. Dynamical models are required to understand better the general outcome of the collisional evolution of the disk.

Acknowledgments

We wish to thank Paul Smith and Karl Stapelfeldt for their helpful comments regarding the Spitzer/MIPS-SED data, and Johan Olofsson for reducing the Spitzer/IRS spectrum. We are grateful to the anonymous referee who suggested improvement to the manuscript. We also thank the Programme National de Planétologie (PNP) and the CNES for supporting part of this research. J.L. and J.C.A. thank the French National Research Agency (ANR) for financial support through contract ANR-2010 BLAN-0505-01 (EXOZODI). This study is based, in part, on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute (STScI), which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with program # 10177. Support for program # 10177 was provided by NASA through a grant from the STScI.

References

Appendix A: Statistical analysis of the results

The strength of the GRaTeR code, in comparison with an optimization algorithm that would provide a unique best-fit solution and possibly miss local minima, is that all the models are stored enabling fast and easy post-processing of the results.

To find the best solution in the grid of models, we use a classical least-square minimization approach with n_free =  3 to 4 free parameters depending on the grain composition and N = 48, which results in 43 to 45 degrees of freedom (d.o.f.) (Table 6). The free parameters are κ, a_min and M_dust, and depending on the situation v_C, v_ice or . M_dust is handled separately because it is automatically adjusted to reproduce the surface density profile in the code.

Let a be the vector of parameters. Here . For each set of parameters, we can define the quantity: $\begin{array}{ccc}{\mathit{\chi }}^{\mathrm{2}}\mathrm{\left(}{a}\mathrm{\right)}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}\frac{\mathrm{(}{\mathit{D}}_{\mathit{i}}\mathrm{-}{\mathit{M}}_{\mathit{i}}{\mathrm{)}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathit{i}}^{\mathrm{2}}}& & \end{array}$(A.1)where D_i and M_i are the observed and modeled values, respectively, and σ_i is the measured uncertainty. For each χ² value, we can define the likelihood of the data, i.e. the probability of the data given the parameters of a model. This assumes the errors are distributed normally and are uncorrelated. $\begin{array}{ccc}\mathit{P}\mathrm{\left(}{D}\mathrm{\right|}{a}\mathrm{)}\mathrm{=}\frac{\mathrm{1}}{\mathit{C}}\hspace{0.17em}\exp \left(\frac{\mathrm{-}{\mathit{\chi }}^{\mathrm{2}}\mathrm{\left(}{a}\mathrm{\right)}}{\mathrm{2}}\right)& & \end{array}$(A.2)where C is a normalization constant and D =  {D(i)} _i = 1...N. In order to evaluate the effect of each individual parameter a_i on the quality of the fit, we implement a Bayesian inference method. We define the posterior probability of a as the probability distribution of the parameters a given the data D. It has to be noted that this formulation implicitly assumes that the model is suited to describe the observation. Bayes’ theorem says that this quantity can be written:

$\begin{array}{ccc}\mathit{P}\mathrm{\left(}{a}\mathrm{\right|}{D}\mathrm{)}\mathrm{=}\frac{\mathit{P}\mathrm{\left(}{a}\mathrm{\right)}\mathit{P}\mathrm{\left(}{D}\mathrm{\right|}{a}\mathrm{)}}{\mathit{P}\mathrm{\left(}{D}\mathrm{\right)}}& & \end{array}$(A.3)where P(D) is the prior probability of the data under all possible hypothesis (a normalization constant in practice), and P(a) is the prior probability distribution of the parameters a. Here we use uniform priors for all the parameters, meaning that we do not favor any realization of the parameters prior to obtaining the data. That fundamental theorem provides a simple way to independently assess the parameters of a model. Indeed, the marginal probability of each parameter a_i can be computed regardless of the values of the other parameters aj ≠ i: $\mathit{P}\mathrm{\left(}{\mathit{a}}_{\mathrm{i}}\mathrm{\right|}{D}\mathrm{)}\mathrm{=}\underset{}{{\mathrm{\int }}^{{\mathit{a}}_{\mathrm{0}}}}\mathit{...}\underset{}{{\mathrm{\int }}^{{\mathit{a}}_{\mathit{i}\mathrm{-}\mathrm{1}}}}\underset{}{{\mathrm{\int }}^{{\mathit{a}}_{\mathit{i}\mathrm{+}\mathrm{1}}}}\mathit{...}\underset{}{{\mathrm{\int }}^{{\mathit{a}}_{{\mathit{i}}_{\max }}}}\mathit{P}\mathrm{\left(}{a}\mathrm{j}\ne \mathrm{i}\mathrm{\right|}{D}\mathrm{)}\mathrm{d}{\mathit{a}}_{\mathrm{0}}\mathit{...}\mathrm{d}{\mathit{a}}_{{\mathrm{i}}_{\max }}\mathit{.}$(A.4)Then we obtain a probability law as a function of the realization of each individual parameter. Finally, we estimate the parameters by calculating the mathematical expectation and variance. These are the values presented in Table 6 where the uncertainties are 3σ, $\begin{array}{ccc}\mathit{E}\mathrm{\left(}{\mathit{a}}_{\mathit{i}}\mathrm{\right)}& \mathrm{=}& \underset{{\mathit{a}}_{{\mathit{i}}_{\min }}}{\overset{{\mathit{a}}_{{\mathit{i}}_{\max }}}{\mathrm{\int }}}{\mathit{a}}_{\mathit{i}}\mathit{P}\mathrm{\left(}{\mathit{a}}_{\mathit{i}}\mathrm{\right|}{D}\mathrm{)}\\ {\mathit{\sigma }}^{\mathrm{2}}\mathrm{\left(}{\mathit{a}}_{\mathit{i}}\mathrm{\right)}& \mathrm{=}& \mathit{E}\mathrm{\left(}{\mathit{a}}_{\mathit{i}}^{\mathrm{2}}\mathrm{\right)}\mathrm{-}\mathit{E}\mathrm{(}{\mathit{a}}_{\mathit{i}}{\mathrm{)}}^{\mathrm{2}}\mathit{.}\end{array}$

Appendix B: Gas line fluxes (see Table B.1.)

All Tables

All Figures

	Fig. 1 PACS 70 μm imaging of (left) HD 181327 and (right) a PSF reference star, HD 148387. Both images were acquired in chop-nod mode. The field of view of each image is 25″ × 25″. Brightness contours at 5, 10, 25, 50 and 75 times the background noise in the HD 181327 image are overlaid on both images.
In the text

	Fig. 2 Radial profiles of HD 181327 at 70 μm (top), compared to HD 148387 as a PSF reference, and at 100 μm (bottom), compared to α Boo as a PSF reference. Both plots show the azimuthally averaged, normalized radial brightness profile vs. radius from the star. Each error bar is the standard deviation of the brightness values at that radius. The FWHM values for the profiles are given on the plot.
In the text

	Fig. 3 ATCA 93.969 GHz (3.192 mm) map of HD 181327 with contours at −2 (dashed red line), 2 and 3σ (solid blue lines), where σ is the image rms given by 1.1 × 10-4 Jy/beam. The dashed blue area represents the beam. The Gaussian fit is overplot with contours at 2 and 3σ (solid green lines). The red and green crosses denote respectively the star position and the Gaussian fit center.
In the text

	Fig. 4 PACS spectrum of HD181327, centered on the wavelength of the [O i] 63 μm emission line. The linear fit to the continuum flux is overplotted with the red line. The [O i] 63 μm emission line was not detected.
In the text

Fig. 5 New HST/NICMOS 1.1 μm surface brightness profiles (top, nominal and  ±1σ profiles), and surface density profiles resulting from direct inversion of the brightness profiles (bottom), assuming two possible values for the asymmetry parameter. The color-coding is the same for both panels. The striped areas identify regions dominated by PSF-substraction residuals where the density has been extrapolated.

In the text

	Fig. 6 Observed spectral energy distribution of HD 181327 (see Tables 2–4).
In the text

	Fig. 7 Mass opacity of spherical grains computed with the Mie theory and Bruggeman mixing rule for several prototypical compositions detailed in Table 6 and their associated best-fit size distributions.
In the text

	Fig. 8 Best fit models (Table 6) for several 1- or 2-material mixtures (top) and 3- to 4-material mixtures (bottom). Red crosses: photometric data (Fig. 6). Solid grey line: synthetic stellar spectrum.
In the text

	Fig. 9 Bayesian probabilities for the models in Table 6: minimum grain size, size distribution power law index, and total dust mass.
In the text

	Fig. 10 Bayesian probabilities for the grain composition in the case of 4-material dust mixtures.
In the text

	Fig. 11 Expected gas line fluxes as a function of gas-to-solid mass ratio for either a low-PAH abundance (thin lines) or a high-PAH abundance (thick lines). The horizontal lines denote the PACS 3σ upper-limits and do not include the 30% flux calibration uncertainty. The gas-to-solid mass ratios assume a solid mass of 0.164 M⊕.
In the text

	Fig. 12 βpr = |FP |/| FG| as a function of grain radius. The horizontal dashed line denotes the limit between bound (β < 1/2) and unbound (β > 1/2) grains.
In the text