© ESO, 2012

1. Introduction

A debris disk around a main sequence star is a collection of small bodies left over from the planet formation process. In our solar system, the asteroid belt and Edgeworth-Kuiper belt are the two best known reservoirs of objects that remain from the planet formation process and range in size from hundreds of kilometers in diameter to meter-scale bodies (e.g. Jewitt et al. 2000; Sheppard & Trujillo 2006). Such reservoirs are highly sculpted by the evolution of the planetary system in which they form (e.g. Petit et al. 2001; Morbidelli et al. 2005; Lykawka et al. 2009), and contain objects whose accretion was stymied by the formation and migration of giant planets in the system, or simply occurred too slowly for them to grow larger. Since debris disks contain a vast number of objects on very similar orbits, they experience a continual collisional grinding which produces and continually replenishes a population of dust. This dust allows us to directly detect debris disks around other stars in two ways. The dust is heated by radiation from the central star, and therefore emits thermal radiation with a temperature characteristic of its distance from its host star (e.g. Aumann et al. 1984; Greaves et al. 2005). In addition, the smallest grains of dust can efficiently scatter the light of the host star (e.g. Smith & Terrile 1984; Kalas et al. 2005). The physical and observational properties of debris disks were defined by Lagrange et al. (2000), and their studies were recently reviewed by Wyatt (2008) and Krivov (2010) and will eventually place our solar system in context (Greaves & Wyatt 2010).

Almost all debris disks detected by the satellites IRAS, ISO and Spitzer (Bryden et al. 2006; Su et al. 2006; Trilling et al. 2008), Hubble Space Telescope (HST; Golimowski et al. 2011), and ground-based telescopes (Wyatt 2008) surround A-type and F, G, K-type stars despite several deep surveys of large samples of M stars conducted from mid-IR to submillimeter wavelengths (Plavchan et al. 2005; Lestrade et al. 2006, 2009; Gautier et al. 2007; Avenhaus et al. 2012). Currently, among the nearby M-stars, the only spatially resolved debris disk is around the very young M1 star AU Mic (Kalas et al. 2004; Liu et al. 2004; Krist et al. 2005; Wilner et al. 2012) which has been modeled by Augereau & Beust (2006) and Strubbe & Chiang (2006). In addition, there are a few candidate disks with excesses above photospheric level (e.g. Smith et al. 2006; Lestrade et al. 2006; Plavchan et al. 2009). Finally, in the cluster NGC 2547 (~40 Myr old and ~433 pc), deep Spitzer MIPS observations have revealed 11 M-stars with 24 μm excesses above photospheric level and no excess at 70 μm; these observations have been interpreted as warm dust in debris disks (Forbrich et al. 2008).

The fact that debris disks are more seldomly observed among M-stars than around higher-mass stars seems surprising at first, since all spectral types have similar detection rates of protoplanetary disks in the earlier stage of their evolution, according to observations of low density clusters like Taurus-Auriga and ρ Oph (e.g., Andrews & Williams 2005). However, in high density clusters like Orion, external photoevaporation by intense FUV radiation field can severely limit the production of planetesimals around low mass M-stars on a timescale shorter than ~10 Myr (Adams et al. 2004). Another hazard for M-stars during the first ~100 Myr is close stellar flybys, when co-eval stars are still in the expanding cluster of their birth and strongly interacting with each other. During these early close stellar flybys, planetesimals are stripped from disks, and this is more severe for disks around low mass stars in high stellar density clusters like Orion according to simulations (Lestrade et al. 2011).

Recently, Wyatt et al. (2012) have found evidence of the prevalence of debris disks in low-mass planetary systems (also Moro-Martín et al., in prep.) and suggest that this correlation could arise because such planetary systems are dynamically stable over Gyr timescales. Recent observations show that low-mass planets are more abundant among M-stars than around the other stars (Bonfils et al. 2011; Howard et al. 2012). Hence, if the correlation between debris disks and low-mass planets for G-stars applies to M-stars, then debris disks should be relatively common around them, in contrast to a paucity of detections.

However, debris disks around M-stars are harder to detect than around more massive stars at the same distance simply because they are less luminous, meaning that the dust within experiences significantly less heating. Therefore, to detect the same disk around a later type star requires deeper observations. M-star debris disks may also be less detectable because additional grain removal processes are operating. For example, a physical pecularity of M-stars is that they are structurally different from solar-type stars. Their interiors have deep convective zones −  fully convective for M3 spectral type and later −  that produce strong coronal magnetic fields responsible for their optical/radio flares and X-ray emission (Hawley et al. 2000). This activity generates also stellar winds of energetic particles (Wargelin & Drake 2001) which might dominate the circumstellar grain removal processes for a large fraction of the star lifetime (Plavchan et al. 2005).

This paper describes observations carried out as part of the key program DEBRIS (Disc Emission via a Bias-free Reconnaissance in the Infrared/Sub-mm) on the Herschel Space Observatory (Pilbratt et al. 2010). DEBRIS is an unbiased flux-limited survey to search for dust emission at λ = 100 and 160 μm toward the nearest ~89 stars of each spectral type A, F, G, K, M as evidence of debris disks (see Matthews et al. 2010; and Phillips et al. 2010, for the sample description). For selected targets, complementary Herschel observations at 70, 250, 350, 500 μm were also conducted. The first results of this program have already shown that these observations can detect disks down to much fainter levels than previously achieved, and moreover can spatially resolve debris disks at far-IR wavelengths (Matthews et al. 2010; Churcher et al. 2011; Kennedy et al. 2012a,b; Wyatt et al. 2012; Booth et al. 2012; Broekhoven-Fiene et al. 2012).

As part of this survey, we have spatially resolved a disk around the M3 spectral type star GJ 581 at λ = 70, 100, and 160 μm. Hence, this is the second resolved debris disk around an M-star, but, in contrast to the star AU Mic which is young (12 Myr, Zuckerman & Song 2004), GJ 581 is old (2–8 Gyr, see Sect. 3). Also, GJ 581 is surrounded by at least four low mass planets with minimum masses of 1.9, 15.6, 5.4, and 7.1 M_⊕, orbital radii of 0.03, 0.04, 0.07, and 0.22 AU, and eccentricities between 0.0 and 0.32, detected by radial velocity measurements (Bonfils et al. 2005; Udry et al. 2007; Mayor et al. 2009; Forveille et al. 2011). All these planets are within the tidal lock region of this M3 spectral type star (<0.25 AU) and hence are expected to be synchronously rotating and potentially undergoing atmospheric instabilities (Wordsworth et al. 2011; Kite et al. 2011). Planets GJ 581c and d are near and in the conventionally defined habitable zone (Selsis et al. 2007), respectively. The presence of one or two additional planets in the system is debated (Vogt et al. 2010; Forveille et al. 2011; Vogt et al. 2012).

In this paper, we describe the Herschel observations of GJ 581 as well as archival MIPS and IRS data from Spitzer, and NICMOS data from HST in Sect. 2. The stellar parameters of GJ 581 used are in Sect. 3. Reconnaissance of a cold debris disk around G J581 in the three PACS images at 70, 100, and 160 μm and in the presence of background sources contaminating the field is described in Sect. 4. Modeling of these images to determine the spatial distribution of the emitting dust is described in Sect. 5. The spectral energy distribution (SED) including the IRS spectrum of GJ 581 and modeling of a hypothetical second component of warm dust are described in Sect. 6. An upper limit on the brightness of scattered light using the NICMOS image is discussed in Sect. 7. Physical conditions in the disk and its relationship with the planetary system around GJ 581 are discussed in Sect. 8.

2. Observations

2.1. Herschel

GJ 581 was initially observed with Photodetector and Array Camera and Spectrometer (PACS, Poglitsch et al. 2010) on 11 August 2010 using the standard DEBRIS observing strategy, and a resolved disk was tentatively detected at 100 and 160 μm. We then acquired deeper PACS images at 100 and 160 μm on 29 July 2011, a PACS image at 70 μm (and 160 μm) on 1 August 2011, as well as SPIRE images at 250, 350 and 500 μm on 30 January 2011. These observations are summarized in Table 1.

2.1.1. PACS

The PACS observations used the mini-scan map mode with eight legs of a 3′ length, with a 4′′ separation between legs in a single scan direction at a rate of 20 arcsec s^-1, and two scan directions (70° and 110°). These data were reduced using the Herschel interactive processing environment HIPE (Ott 2010) version 7 and implement version FM6 of the flux calibration. The data were pre-filtered to remove low-frequency (1/f) noise using a box-car filter with a width of 66 arcsec at 70 and 100 μm and 102 arcsec at 160 μm. This data filtering results in the source flux density being underestimated by ~20 ± 5% as discussed in detail by Kennedy et al. (2012a). Maps were made from these filtered timelines using the photProject task in HIPE.

The pixel scales in the images presented in Fig. 1 were set to 1 arcsec at 70 and 100 μm, and 2 arcsec at 160 μm, i.e., smaller than the natural pixel scales. This enhanced sampling is possible because of the high level of redundancy provided by the scan map mode used but it comes at the cost of correlated noise between neighboring pixels. We have also made images with the natural pixel scales of 3.2 arcsec at 70 and 100 μm, and 6.4 arcsec at 160 μm to evaluate the impact on the parameter estimation in our modeling. The noise rms for the images with the natural pixel are 0.47 mJy/5.6′′ beam at 70 μm, 0.48 mJy/6.7′′ beam at 100 μm, and 0.77 mJy/11.4′′ beam at 160 μm.

Fig. 1 PACS images of GJ 581 cropped to ± 50 arcsec from the star, at 70, 100 and 160 μm from top to bottom, respectively. In the left-hand column, the raw images show that the main emission is centrally located about the star position (image center) and that there are several point-sources in the field, barely detected at 70 μm, significantly at 100 μm, and more prominently at 160 μm, as expected for submillimeter background galaxies. In Sect. 4, we show that the main emission is extended and centered on the star position, as expected for a debris disk, and mingles with a background source ~11 arcsec toward the northwest. The panels in the middle column are the photosphere-subtracted images. The panels in the right-hand column show the best subtraction (lowest residuals) of a two-point source model, which assumes that there is no debris disk around GJ 581 but an extra background source located exactly behind the star in addition to the N W source. This model is rejected because of the systematic residuals left, indicative of an extended structure, especially at 70 and 100 μm. At each wavelength, the contours levels of the three images are the same and correspond to 1, 2, 3, 9, 15σ0 (σ0 = 0.0135 mJy/1′′ pixel) at 70 μm, 1, 2, 3, 6, 9, 12, 15σ0 (σ0 = 0.0094 mJy/1′′ pixel) at 100 μm, and 1, 2, 3, 5, 7, 9, 11σ0 (σ0 = 0.0251 mJy/2′′ pixel) at 160 μm. The coordinates of the image center provided in the labels correspond to the star position at epoch of observation (July 29th–August 1st 2011). The hatched circles are the beam FWHMs: 5.6, 6.8, and 11.4 arcsec at 70, 100, and 160 μm, respectively.

2.1.2. SPIRE

Follow-up observations were taken on 30 January 2011 with SPIRE (Spectral & Photometric Imaging REceiver, Griffin et al. 2010) using the small-map mode, resulting in simultaneous 250, 350 and 500 μm images. The data were reduced using HIPE (version 7.0 build 1931), adopting the natural pixel scale of 6, 10, 14 arcsec at 250, 350 and 500 μm respectively. The noise rms are 6.1 mJy/18.2′′ beam, 7.9 mJy/24.9′′ beam, and 8.3 mJy/36.3′′ beam at 250, 350 and 500 μm, respectively, and the image at 250 μm is shown in Sect. 4.4.

2.2. Ancillary data

2.2.1. Spitzer

MIPS 70 μm observations of GJ 581 (AOR 22317568) were taken on 21 August 2007 (no 24 μm MIPS were taken) and a small measured excess, with the significance in Kóspál et al. (2009) and χ₇₀ = 2.2 in Bryden et al. (2009) with the same data, forms a tentative discovery. We re-reduced the archival data using an updated pipeline and the flux calibration summarized in Gordon et al. (2007) providing the new flux density 20.0 ± 5.3 mJy by PSF fitting to the image at the effective wavelength of 71.42 μm (color correction for T_dust = 40 K applied: 0.89¹). The uncertainty includes both statistical and calibration uncertainties and the corresponding excess ratio χ₇₀ is 2.7 with our estimate of the photospheric flux density of 5.6 mJy at 71.42 μm (see Sect. 6). Photospheric flux densities predicted for late type stars (K and M) by the Kurucz or Next Gen models have been shown to be overestimated in the mid-IR by as much as 5–10% (Gautier et al. 2007; Lawler et al. 2009). Hence, this excess can be treated as a lower limit.

The IRS observations of GJ 581 (AOR22290432) were taken on 31 August 2007, and the details of the data reduction are in Beichman et al. (2006) and Dodson-Robinson et al. (2011).

2.2.2. HST/NICMOS

GJ 581 was directly imaged with HST/NICMOS on 6 May 1998 (GO-7894; PI Todd Henry). The NICMOS data and the overall observing program are described in Krist et al. (1998) and Golimowski et al. (2004). We reanalyzed the F110W data for GJ 581 consisting of 128 s of cumulative integration on the NIC2 camera (0.076′′/pixel, 256 × 256 pixels). Target stars were not placed behind the occulting spot, near-contemporaneous observations of PSF reference stars were not made, and multiple telescope roll angles were not employed. Therefore the observations were not optimized for high-contrast imaging of low surface brightness circumstellar nebulosity. Nevertheless, we subtract the GJ 581 point-spread-function (PSF) using observations of LHS 1876 (GJ 250B) made on 24 March 1998 as part of the same scientific program, and in so doing, set constraints on the scattered light disk brightness as discussed in Sect. 8. PSF subtraction techniques, including a discussion of scattered light artifacts and other spurious features, are described in greater detailed by Krist et al. (1998).

3. Stellar parameters of GJ 581

GJ 581 (HIP 74995) lies relatively nearby (6.338 ± 0.071 pc; Phillips et al. 2010) and is classified as a star of spectral type M3.0 (Reid et al. 1995). Recent CHARA interferometric measurements of its physical radius (0.299 ± 0.010 R_⊙) imply an effective surface temperature of T_eff = 3498 ± 56 K, a bolometric luminosity of 0.01205  ±  0.00024 L_⊙ and a stellar mass of 0.28 M_⊙ (von Braun et al. 2011).

A variety of different techniques have been discussed in the literature as a means to determine the age of GJ 581, including kinematics, magnetic activity (X-ray observations), chromospheric activity, stellar color, metallicity and rotation. Leggett (1992) finds that the galactic velocities of GJ 581 are intermediate between those typical of the young and old galactic disk M-stars. Bonfils et al. (2005) conclude that the low limit on its X-ray emission, the low v sin i and the weak CaII H and K emission, taken altogether, suggest that GJ 581 is at least 2 Gyr old. Selsis et al. (2007) established an L_x/L_bol versus age relation for M- K- G-spectral type stars to estimate that the age of GJ 581 could be around 7 Gyr. Recently, Engle & Guinan (2011) have established an age-rotation period relation for M-stars and determined an age of 5.7 ± 0.8 Gyr for GJ 581. Clearly, GJ 581 is an old star well above 1 Gyr.

High contrast imaging for GJ 581 has revealed no companion of ~7 Jupiter masses or higher between 3–30 AU (Tanner et al. 2010). In addition, the limits provided by the HARPS radial velocity measurements exclude planets that are more massive than Jupiter with semimajor axes inside 6 AU (Fig. 13 in Bonfils et al. 2011).

All the parameters of GJ 581 are summarized in Table 2.

4. Herschel images of GJ 581

In Fig. 1, we present our deep PACS images of GJ 581 cropped to the region ± 50′′ from the star. At each wavelength, the main emission is close to the star position (image center) and the surrounding field is contaminated by several other sources, detected barely at 70 μm, significantly at 100 μm, and prominently at 160 μm, as is expected for submillimeter galaxies in the background. The central emission is suggestive of a spatially resolved debris disk not fully separated from the background source to the N W. In the next subsections, we analyze quantitatively these PACS images to verify this view.

Fig. 2 Radial profiles of the mean brightness of the photosphere-subtracted images. Each point of these curves was computed as the mean of the emission in an elliptical annulus centered on the peak of the main emission close to the image center. They correspond to an axisymmetric disk model inclined to 40° relative to the plane of the sky (see text). They are one pixel wide, i.e. 1 arcsec in the 70 and 100 μm images and 2 arcsec in the 160 μm image. The Gaussian with the FWHM of the beam is the profile expected for a hypothetical point source in the background aligned by chance with the star and shown for comparison.

4.1. Radial profiles of the emission in the PACS images

First, in Fig. 2, we present the radial profiles of the emission at the three wavelengths by computing the mean brightness in ellipical annuli centered on the peak of the main emission. At the three wavelengths, these radial profiles show that the emission is extended about this peak when they are compared to the Gaussian profiles of a hypothetical point source in the background.

We computed these profiles after subtracting from each image a PSF scaled to the photosphere flux density (S_photosphere = 5.8,2.8,1.1 mJy at 70, 100 and 160 μm, respectively, see Sect. 6). The emission peak position was found to be less than 2 arcsec from the star in each image, consistent with the pointing accuracy of the Herschel telescope². We have tested elliptical annuli at PA = 120° and with inclinations 0° (circular), 25°, 40° and 75°, anticipating that the disk may not be in the plane of the sky. We found that all these radial profiles were more extended than the Gaussians at the three wavelengths, and slightly more for the inclination of 40°. The imprint of the N W source at the radial distance of 11 arcsec can be seen in these profiles at 100 and 160 μm. We have tested the method by computing the radial profile of the emission of the southeast background source in the same way. We satisfactorily found that its profile matches within 1σ the Gaussian expected for a point source.

The extended emission revealedby these profiles can also be seen directly in the photosphere-subtracted images of Fig. 1 (middle column), most prominently at 70 μm because the photospheric flux density is highest at this wavelength.

4.2. Gaussian source fits

As a second approach to verify that the central emission is more extended than the PSF, we fit an elliptical 2D Gaussian to each photosphere-subtracted image with masking applied to the position of the N W source at 100 and 160 μm. At 70 μm, we found FWHM of the minor and major axes of 9.5 ± 0.7′′ and 10.7 ± 1.1′′; a position angle of 120 ± 10°; and a flux density of 18.8 ± 1.4 mJy after adding back in the photospheric contribution. At 100 μm, we found, respectively, values of 9.9 ± 1.0′′, 13.3 ± 0.5′′, 120 ± 10°, and 21.1 ± 1.5 mJy for the minor and major axes FWHM, the position angle, and the 100 μm flux density, after masking the N W source (all pixels in a 12′′ × 12′′ square centered at − 9′′ and + 6′′ from the star). Given that the FWHM of the PACS PSFs are 5.6′′ and 6.8′′ at these wavelengths, we conclude that the emission is significantly extended, as already shown by the radial profiles, and is elongated at PA ~ 120°.

The ICRS coordinates of the 70 μm Gaussian peak in this fit are 15^h10^m 25.905 ± 0.05^s and − 7°43′22.79 ± 0.7′′ and differ only by 0.3′′ from the adjusted position of the 100 μm Gaussian peak. These coordinates differ by + 0.66′′ and − 1.48′′ from the right ascension and declination of the star GJ 581 predicted with the Hipparcos astrometric parameters (Table 2). These differences are consistent with the 1σ pointing accuracy of 2′′ for the Herschel telescope. We conclude that the main emission in the PACS images at 70 and 100 μm is centered on the star position within pointing uncertainty.

The 160 μm image is the product of the coaddition of two images taken independently with the PACS 70/160 and PACS 100/160 instruments only a few days apart in 2011 (Table 1). The registration of these two images were facilitated by the negligible displacement due to proper motion over the short lapse of time and by the fortuitously small difference of 0.3′′ between the pointing positions of the two instruments as found above. Our first 160 μm image in 2010 was not coadded because no feature in the image could be used to check the registration since its signal-to-noise ratio was times lower. Although this first image was crucial in our decision to observe deeper, its use or not is inconsequential for our analysis. For the fit at 160 μm, we fixed the position of the Gaussian to the coordinates determined at 70 μm. This was necessary because the large mask (16′′ × 16′′ square) used for the N W source affected the independent determination of this position. The best fit 2D Gaussian parameters were minor and major axes of 12.8 ± 1.5′′ and 21.5 ± 2′′; a position angle of 125 ± 10°; and a flux density of 22.1 ± 5.0 mJy, after adding back in the photospheric contribution. Given the PACS PSF at 160 μm of 11.4′′, this indicates that the emission is extended at this wavelength as well.

The disk inclinations resulting from the ratios of the minor and major axes determined above, and corrected quadratically for the convolved PSFs, are: 33 ± 17° (0° is face-on), 54  ±  6°, and 71  ±  7° at 70, 100 and 160 μm respectively. Although scattered, these values are statistically consistent, since they are within 1.5σ from their weighted mean (59°), and indicate an inclined disk which has implications for the masses of the planets of the system as discussed in Sect. 8.3.

We note that the three major axes above, corrected quadratically for the convolved PSF, are very closely proportional to the wavelength, suggesting that the disk is radially broad since emission at longer wavelength probes colder dust, more distant from the central star.

The flux densities from our fits above have been scaled up to account for the flux removed by the data filtering during the reconstruction of the images; the correction factors are 16, 19 and 21% with an uncertainty of 5% estimated for point sources in the DEBRIS survey by Kennedy et al. (2012a). The uncertainty of each PACS flux density determined above is based on the quadratic sum of the statistical uncertainty in our Gaussian fit, the absolute flux calibration accuracy of 3% at 70 and 100 μm, and 5% at 160 μm, provided by the Herschel project³, and the 5% uncertainty of the correction factor for the data filtering. The uncertainty of the flux density at 160 μm is formally 2 mJy with this calculation but our fit depends to some degree on the position of the mask applied for the N W background source. So we have increased this uncertainty to 5 mJy at this wavelength based on several fits with different masks.

4.3. Considering the superposition of two point-sources

As a final test, we consider the possibility that the extended emission in the central part of the image could be caused by the superposition of two backgound sources instead of a disk. To this end, we subtracted two PSFs from each PACS image (in addition to subtracting the phostospheric emission), adjusting their flux densities in order to remove as much emission as possible. The first PSF was located at the position of the N W source, i.e. at − 9′′ and + 6′′ from the star, and the second was tested at six locations; the star position itself, as well as a half FWHM to the north, south, east, and west of the star, and half way between the star and the N W source. The lowest residuals were found after removing 1.4, 2.9 and 6.6 mJy for the first PSF at 70, 100 and 160 μm, respectively, and 4.9, 6.9 and 9.4 mJy for the second PSF at the star position at 70, 100 and 160 μm, respectively. Note that these latter flux densities are free of the photosphere contributions estimated in Sect. 6.1. Despite this removal process, there is still significant structure left in the residual images shown as the two-point source subtracted images in Fig. 1 (right-hand column). This structure can be best explained as resulting from the extended emission of the disk incompletely removed by this process. Hence, we conclude that this test rejects the possibility that the superposition of two background point sources can be responsible for the central emission.

We elaborate further by discussing the probability to find such contaminant sources in the field and their spectra. First, the probability to have one background source stronger than 6.6 mJy at 160 μm within 11′′ from the star is 18%, and to have two is only 1.8% by using the Poisson probability distribution with the mean source surface density N(S > 6.6 mJy) ~ 4000 sources/deg² at 160 μm provided by the Herschel PEP survey Berta et al. (2011, see also Sibthorpe et al. 2012). Second, the spectra of these test sources removed from the images may or may not be physical. The flux densities removed at the N W source position (1.4, 2.9 and 6.6 mJy at 70, 100 and 160 μm respectively) are consistent with the spectrum of a galaxy at z > 1.5 according to Fig. 4 of Blain et al. (2002) valid for a typical high-z galaxy enshrouded in dust (~38 K, L  ~  5  ×  10¹² L_⊙) and radiating in the far-IR and submm. However, the flux densities removed at the star position (4.9, 6.9 and 9.4 mJy at 70, 100 and 160 μm respectively) above the photospheric levels make the ratio S_{100 μm}/S_{70 μm} lower than expected for a galactic spectrum according to that work.

4.4. SPIRE images

The SPIRE image at 250 μm in Fig. 3 shows an elongated structure at PA ~ 120° which has two peaks at the 2σ level that match the positions of the star and the N W source. This structure is also extended to the S E. Although this structure is of low statistical significance, it is consistent with the emission detected at the PACS wavelengths. The two other SPIRE images, at 350 and 500 μm, are dominated by noise and no reliable structure can be recognized. A 2D Gaussian could not satisfactorily model the emission at 250 μm, and so we carried out photometry with an aperture of 36 arcsec and measured a flux density of 24 ± 6 mJy. This flux density is considered an upper limit for the disk since its emission is blended with that of the N W source.

Fig. 3 SPIRE image of GJ581 at 250 μm. Pixel size is 6′′ and the contour levels are 1 and 2σ with σ = 6.1 mJy/18.2′′ beam. The star symbol is the position of GJ 581 at the date of observation.

4.5. Spitzer MIPS image at 70 μm

A Spitzer MIPS image at 70 μm was made on 21 August 2007, four years before our PACS image (1 August 2011). We fit a Gaussian with the FWHM of 19 arcsec (70 μm MIPS beam) to the emission of the MIPS image and found the coordinates of the Gaussian peak to be at α = 15^h10^m 26.08 ± 0.17^s and δ = −7°43′23.2 ± 2.5′′ in the ICRS system. The relatively high noise of the MIPS image did not permit a solution for the FWHM. The differences in the coordinates compared with our 70 μm PACS position given in Sect. 4.2 are Δα = −2.56′′ and Δδ = + 0.45′′ (PACS minus MIPS) with an uncertainty of 3.3′′ when combining the astrometric uncertainty of the MIPS Gaussian fit (2.5′′ in both coordinates) and the pointing accuracies of 2′′ for Herschel and of 1′′ for Spitzer⁴. The predicted displacement of GJ 581 between the two epochs of observations is Δα = −4.98′′ and Δδ = −0.41′′ computed with the proper motion in Table 2. Hence, the coordinate differences of the 70 μm emission measured between the two epochs are compatible with this prediction but the star has not moved sufficiently between these epochs for us to confirm that the 70 micron emission is comoving with the star at a statistically significant level.

Fig. 4 Maps of residuals for the best fit models of the disk to the PACS images at 70, 100, and 160 μm from left to right. The post-fit residuals are between ± 3.5σ0, coded as black, blue, orange, yellow and white over this range (σ0 is the same noise rms as used for contours in Fig. 1 but the color scale is not the same as in Fig. 1). In zooming the electronic version, the contours apparent in these maps are −3, −2, −1, 1, 2, 3 σ0 (dashed contours are negative levels). The two background sources are masked in the 100 and 160 μm images as the dotted squares show.

5. Modeling the PACS images of GJ 581

5.1. Parametrized model

We fit a parametrized model of the disk to the PACS images. The model is axisymmetric and truncated by the inner radius r_in and the outer radius r_out which are free parameters in our fit. Its dust emission is optically thin, and the flux density from each element (k,l) of the grid covering the disk is (1)where B_ν(T(r_k,l)) is the Planck function that depends on the grain temperature T(r_k,l) at the radial distance r_k,l from the star, Δa is the area of the element in the grid, d is the distance to the star, and Σ_p is the coefficient of the power-law (e.g. Wyatt et al. 1999)⁵. To fit individually each PACS image in Sect. 5.2, we set the factor ϵ to unity in Eq. (1). To fit simultaneously the three PACS images at λ = 70 μm, 100 μm, and 160 μm in Sect. 5.3, we implement a gray body effect (e.g. Dent et al. 2000) by setting ϵ to unity if λ < λ₀, and ϵ = 1.0 × (λ₀/λ)^β if λ > λ₀, where λ₀ and β( > 0) are free parameters in our fit.

In our model, no assumption is made for the size distribution of the grains, their mineralogical composition and porosity. The thermal structure of the disk is taken as (2)where f_T is a scaling factor applied to the black body temperature T_BB(r) = 278(L/L_⊙)^0.25(r/1 AU)^-0.5 (K), and is a free parameter in our fit. Here we illustrate how to interpret this parameter using a simplified model of the absorption and reemission of the starlight by the grains. For dust with the absorption and emission efficiencies Q_abs and Q, a straightforward derivation shows that f_T = (Q_abs/Q)^1/4 for a grain at thermal equilibrium, ignoring that these efficiencies should be averaged over the spectrum of incoming and outgoing radiation, and integrated over the dust size distribution. If we assume that grains larger than 1 μm absorb starlight with the efficiency Q_abs = 1 and reemit it at longer wavelengths at the lower efficiency Q, then the simple interpretation is that f_T is directly related to the emission efficiency through the relationship f_T = Q ^− 1/4.

The term Σ_pr^α in Eq. (1) is the emitting cross-sectional area of the grains per unit area of the disk surface. These grains are spatially distributed according to a radial profile taken as the power-law Σ_pr^α, and their total cross-sectional area is A. Since these grains reemit with the efficiency Q, their total emission is proportional to . Hence, if α ≠ − 2, the coefficient of the power-law is (3)and the total cross-sectional area A and the power-law index α are free parameters in our fit⁶. Grains smaller than 1 μm are not important because they emit so inefficiently that their flux density is negligible (Bonsor & Wyatt 2010).

The total cross-sectional area A can be converted into mass assuming a size distribution and a mass density ρ for the material. Adopting the standard size distribution n(D) ∝ D^-3.5 for spherical particles of diameter D between D_min and D_max, the corresponding mass is (4)This model is complemented with a point source photosphere centered on the image by two free parameters (coordinates x_c and y_c) and having the flux densities estimated from the Next Gen stellar atmosphere model in Sect. 6 but lowered by ~ 20% because of the data filtering used to reconstruct the images as already mentioned in Sect. 2.1.1. This model is projected onto the sky with the inclination i and the node orientation Ω, and finally convolved by the telescope PSF provided by the PACS images of the reference star α Boo at 70 μm, 100 μm and 160 μm. Hence, our model has 9 free parameters (r_in,r_out,α,f_T,A,i,Ω,x_c,y_c).

The best fit is found by minimizing (5)computed with the residuals between the image (O) and the model (C) over all the pixels of the image, and assuming the same measurement uncertainty σ₀ for all the pixels. To compute the model, we used a grid on the sky which has a resolution of 0.5′′ at 100 μm and 70 μm, and 1′′ at 160 μm, i.e. twice as fine as the pixel size of the images of Fig. 1, and has dimensions 128 × 128 at 70 and 100 μm, and 64 × 64 at 160 μm. These dimensions can accommodate the largest disk model tested (r_out = 150 AU) extended by twice the beam FWHM. This sky grid is the same for all the models tested so that the number of degrees of freedom ν is the same for all of them. To use the probability distribution to discriminate between them, we carefully estimate the a priori uncertainty σ₀ by computing the noise rms over the image limited to the sky grid dimensions and after excluding the central part (r < 25′′) where the disk emission is located. For GJ 581, the noise rms, σ₀, is 0.0135 mJy/1′′ pixel at 70 μm, 0.0094 mJy/1′′ pixel at 100 μm, and 0.0251 mJy/2′′ pixel at 160 μm, corresponding to 0.49 mJy/5.6′′ beam, 0.50 mJy/6.7′′ beam, and 0.81 mJy/11.4′′ beam, respectively, by using the beam area π × FWHM²/4ln2.

Finally, we used the SPIRE 3σ upper limits as a constraint to reject any model with a flux density of the dust emission larger than 24 mJy at 250 μm.

5.2. Fits of individual PACS images

First, we searched for the best fit model for each individual image. The ranges of the model parameters tested were: A of 1 to 20 AU²; r_in of 3.1 to 80.0 AU; r_out of r_i to 150.0 AU; α of − 3.0 to 0.0; f_T of 1.0 to 6.0; i of 0.0 to 90.0° (0.0° is a disk seen face-on); and Ω of 0.0 to 180.0° (0.0° is north and increasing Ω is east). The range for the index α was chosen to cover possibilities such as grains being blown out of the system by radiation pressure (α = −1), and having a distribution akin to that of the minimum mass solar nebula (MMSN) (α = −1.5), as discussed in the modeling of the disk around the A-star β Leo by Churcher et al. (2011). The two background sources, N W and S E of GJ 581, were masked as shown in the residual maps of Fig. 4. Roughly 50 million models were tested in our search for the best fit.

We found the reduced 1.04, 0.99, and 0.97 for the best fits to the three images at 70, 100, and 160 μm, respectively. The numbers of degrees of freedom are ν = 408 at 70 and 100 μm, and ν = 90 at 160 μm. These reduced indicate noise-like post-fit residuals according to χ²-statistics. The residual maps in Fig. 4 do not show any systematic residuals as expected in these conditions. We stress that the uncertainty σ₀ used for Eq. (5) is the value determined a priori and is not purposefully tweaked a posteriori to make the reduced close to unity.

Fig. 5 Map of the reduced showing the correlation between the temperature factor fT and the total cross-sectional area A of the dust. The minimum is the white region in this map.

The best fit values of the parameters are in Table 3. There are significant correlations between parameters, especially between A, α, r_outer and f_T, as we found by inspecting the two-dimensional projections of the hypersurface, e.g. Fig. 5 for the pair A and f_T. To estimate the parameter uncertainties in these conditions, we have determined the lower and upper limits around the best fit value of each parameter that correspond to the fits in which the reduced are increased to 1.12 and 1.25 with all the other parameters freely adjusted. These thresholds correspond to a probability of 5% in -statistics that the reduced  of pure noise exceeds 1.12 and 1.25 for the number of degrees of freedom ν = 408 and 90, respectively. This is a standard criterium in fitting procedures. We have also inspected the corresponding residual maps and noticed nascent systematics for these degraded fits as expected. For the outer radius r_out, only the lower limits and the best fit values of the fits are given in Table 3 because the upper limits are not well constrained since any distant dust becomes very cold, even accounting for the incident interstellar radiation field as a source of heating (Lestrade et al. 2009, Fig. A1). The resulting range for r_in and r_out does not permit a conclusive estimate of the radial breadth of the GJ 581 disk.

5.3. Combined fit of the three PACS images

To consolidate these results and to break correlations between parameters, we combined the three PACS images at λ = 70, 100, and 160 μm in a single fit by setting the factor ϵ of Eq. (1) to unity if λ < λ₀ and to ϵ = 1.0 × (λ₀/λ)^β if λ > λ₀, in order to implement a gray body effect. Our search covered successively the combinations of β = 0.0, 0.5, 1.0, 1.5, 2.0 and λ₀ = 70 μm, 85 μm, 130 μm. The ranges of the other model parameters tested were the same as for the individual images in Sect. 5.2. The two background sources were masked. Note that the 160 μm image with 4 times fewer pixels has a lower weight than the two other images in this combined fit.

The best fit has the reduced of 1.03 (ν = 922) for β = 0, indicating formally no gray body effect for λ₀ < 160 μm and so a disk dominated by large grains. However, there is a high correlation between α and the pair (β, λ₀) so that these parameters cannot be properly constrained in reality. In fact, in the discussion below, we argue that small grains should be abundant in the disk by comparing timescales of collision and removal processes. The best fit values of the other parameters are in Table 3 and their uncertainties were determined as described in Sect. 5.2. The total cross-sectional area of the dust A = 2.3 AU² can be converted to the dust mass in M_⊕ for a collisional cascade, using Eq. (4) with ρ = 1.2 g/cm³ for icy grains and D_min = 1 μm. The maximum diameter D_max is unconstrained although objects larger than 10 cm contribute negligibly to the emission at the wavelengths considered in this paper. Nonetheless, the size distribution probably extends beyond 10 cm as discussed in Sect. 8.2. The inclination could be anywhere within a relatively broad range (30° < i < 70°) that matches the purely geometrical determination based of the ratio of the major and minor axes of the Gaussians fit to the central emission in Sect. 4.2. The inner radius is 25 ± 12 AU potentially providing an indication of the scale of the planetary system around GJ 581. In a similar way to the fits of the individual images, we cannot distinguish between a relatively narrow ring and a disk extending beyond 100 AU with this combined fit. The best fit value of f_T is , making the dust temperature between 50 and 30 K over the extent of the disk, despite the low luminosity of GJ 581. This factor is partially correlated with A as shown in Fig. 5 but it is clearly inconsistent with unity as we have established by forcing f_T = 1.0 in the model and found a reduced as high as 1.92 for the best fit with this constraint. An emission model including a grain size distribution instead of a fixed f_T as in our current model would be more realistic, but this can be best implemented only if the SED is finely sampled spectroscopically from mid-IR to submm (e.g. Lebreton et al. 2012).

5.4. Model with a Gaussian profile for the grain surface density

Finally, instead of a power law for the radial distribution of the grain surface density, we tested a Gaussian profile Σ_gexp( − 0.5 × ((r − r_g)/w_g)²) peaking at radius r_g and having FWHM of . The sky grid for the model, the calculation of , and the ranges of model parameters tested for A, f_T, i, and Ω were as in Sect. 5.2. Values of r_g  ranged from 2 × w_g to 150 AU, and w_g  ranged from 3.1 AU to r_g/2 (the Gaussian profile is truncated to 2w_g toward the star). We modeled the three PACS images individually and in combination. The resulting best fits have reduced of 1.03, 0.94, and 1.05 at 70, 100, and 160 μm, and of 1.04 for the combined images with β = 0 (no gray body effect). The residual maps are featureless. These best fits are statistically indistinguishable from those with the power law presented in Sects. 5.2 and 5.3. The resulting parameters are: r_g = 52 ± 15 AU, FWHM = 38 ± 15 AU, A = 2.5 ± 1.2 AU², f_T = 3.0 ± 0.5, i < 60°, and Ω = 120 ± 20°.

In this model, the inner part of the system is populated with dust making the inner radius determined with the power law surface density in Sect. 5.3 less definitive.

6. The SED and IRS spectrum of GJ 581

Fig. 6 SED of the cold disk model. The left-hand figure shows the best fit to the three PACS images only (, Table 3). The right-hand figure shows the best fit to both the three PACS images and the IRS spectrum . The modeled cold dust emission is the blue curve, the Next Gen stellar atmosphere spectrum is the gray curve, and their sum is the green curve. The IRS spectrum is in red. The upper inset zooms on the IRS wavelengths and displays spectra as Sν × (λ/24)2 on a linear scale for clarity. The lower inset zooms on the three PACS bands. In the left-hand figure, the best fit model satisfactorily fits the PACS data as shown by the lower inset but misses the IRS data as shown by the upper inset. In the right-hand figure, the best fit model partially misses the PACS data but satisfactorily fits the IRS data.

We present the SEDs of the star GJ 581 and modeled dust emission that we used to determine the fractional dust luminosity L_dust/L_∗ ~ 10^-4 of the disk. The archival IRS spectrum shows a marginally significant excess above the photospheric level that provides additional constraints on the dust emission. However, we show that a single cold disk model and a two component disk model cannot be distinguished to explain this 2-σ excess.

6.1. The SED and the fractional dust luminosity of the disk

The photometry data collected for the SEDs are summarized in Table 4. The flux densities have been color corrected when required⁷. The SED of GJ 581 is based on the NextGen stellar atmospheric model (Hauschildt et al. 1999), with the value log (g) = 5.0 and the effective temperature 3500 K, fit to the Johnson UBV and Cousins RI photometry, the JHK_s photometry from 2MASS, and the recent photometry from AKARI and WISE. Note that the flux densities of the photosphere used for our modeling in Sect. 5 were predicted from this fit (5.8, 2.8 and 1.1 mJy at 70, 100 and 160 μm, respectively). In Fig. 6 (left-hand panel), we show this SED for the star and the SED for the dust emission from our modeling. The fractional dust luminosity was determined by integrating the SED of the dust emission and is L_dust/L_∗ = 8.9 × 10^-5. This value is consistent with the fractional dust luminosity Q_absA/4πr² = 9.9  ×  10^-5 determined from the cross-sectional area of the grains A = 2.3 AU² from our fit in Table 3, using the mean disk radius r = (25 + 60)/2 = 43 AU, and assuming the absorption efficiency Q_abs = 1 for the grains larger than 1 μm. The agreement between these two independent determinations of the fractional dust luminosity provides a self-consistency check of our modeling. This fractional dust luminosity is higher than that of the Kuiper belt by several orders of magnitude.

6.2. IRS spectrum

6.2.1. Synthetic photometry

The Spitzer IRS spectrum is superimposed on the star’s SED in Fig. 6. As is standard with IRS spectra, the short wavelength module SL (7.6 − 14.2 μm) has to be adjusted to the predicted photosphere, and IRS flux densities were scaled up by the factor 1.066 for GJ 581. In Fig. 6 and insets, a small excess is apparent above the photospheric level at the longest wavelengths of the spectrum (module LL1: 20.4 − 34.9 μm).

We have carried out synthetic photometry with a rectangular bandpass between 30 and 34 μm and between 15 and 17 μm which gives the widest wavelength range while still inside of the Long-Low IRS module as prescribed in Carpenter et al. (2008, 2009). We computed the synthetic flux densities S_{31.6 μm} = 32.3 ± 1.9 mJy (IRS) and 28.4 mJy (Next Gen) yielding the 2σ excess 3.9 ± 1.9 mJy, and S_{15.96 μm} = 110.7 ± 0.85 mJy (IRS) and 109.2 mJy (Next Gen) yielding the lower significance excess 1.5 ± 0.85 mJy. We computed these synthetic flux densities as the weighted mean of the data points in these bands, and using the same weights for the corresponding Next Gen synthetic flux densities. The IRS flux density uncertainty includes an absolute calibration error of 6%. Photospheric flux densities predicted for late type stars (K and M) by the Kurucz or Next Gen models have been shown to be overestimated in the mid-IR by as much as 3–5% (Gautier et al. 2007; Lawler et al. 2009). Hence, the significance of the marginal excess at 31.6 μm is likely higher in reality. If real, this excess for the mature M-star GJ 581 is notable because, even among A-type and solar-type stars, 24  μm excesses are less frequent than 70 μm excesses and decrease with age (Rieke et al. 2005; Trilling et al. 2008; Löhne et al. 2008). In the next two sections, we investigate the implications for the system around GJ 581 if this excess is real.

Fig. 7 SED of the best fits of the cold disk and warm belt models. The green curve corresponds to the warm and cold dust emissions added to the Next Gen stellar atmosphere spectrum (gray). IRS spectrum is in red. We also show separately the best fit of the warm belt emission (rw = 0.2 AU and mw = 2.8 × 10-6 M$) (yellow) and the best fit of the cold disk emission (parameters of the combined fit are in Table 3) (blue). Insets and photometric data points are the same as in the legend of Fig. 6.

6.2.2. Modeling the IRS and PACS data with the cold disk model

First, we fit the single cold disk model of Sect. 5 simultaneously to the three PACS images and the IRS spectrum, minimizing  where  and  are the reduced  for the PACS and IRS data, respectively. With this definition, both data sets have the same weight in the fit. The best fit model thus obtained is characterized by , resulting from  and , and its SED is shown in Fig. 6 (right-hand panel). The main parameter changes are f_T = 5.5 and A = 0.8 AU², instead of f_T = 3.5 and A = 2.3 AU² in Table 3. This value of  is higher than  of the best fit in this Table and is high for the number of degrees of freedom of 1186 in χ²-statistics (probability = 1% of pure noise). It is instructive to compare the SEDs of these two fits in Fig. 6; the simultaneous fit to the PACS and IRS data in the right-hand panel does appear to be skewed to some degree. The assumption in our current model that the temperature does not depend on grain size and wavelength is a limitation. A size distribution would broaden the SED, and it may improve the ability to fit a single disk model to the flux densities of the IRS spectrum and the PACS bands simultaneously.

6.2.3. Modeling the IRS and PACS data with a two component model

We explore another possibility, a two component model in which a belt of warm dust is added to our cold disk model of Sect. 5. The model of this belt is simply based on blackbody grains (i.e., we set f_T = 1 for the warm component) located at radius r_w and having a total cross-sectional area A_w. This two parameter model is fit to the IRS spectrum alone by minimizing . We found r_w = 0.2 AU (T_dust = 191 K) and A_w = 7 × 10^-5 AU², giving a corresponding dust mass of m_w = 2.8 × 10^-6 M_$, assuming the standard grain size distribution ( ∝ D^-3.5) between 1 μm and 1 mm-sized particles and ρ = 3 g/cm³ (noting the dependence of this estimate on the unknown maximum size given in Eq. (4)). However, acceptable fits could also be found for r_w between 0.05 AU (T_dust = 382 K) and 0.4 AU (T_dust = 135 K), encompassing the orbital radii of the planets GJ 581c and GJ 581d. The IRS data alone cannot constrain f_T but if this parameter were larger than unity, the dust would be  times further out than the r_w quoted above for the corresponding T_dust. The SED of the two component model is shown in Fig. 7 where we had to decrease the cold dust cross-sectional area A by 6% from its value of Table 3 to account for the warm dust contributions to the flux densities at 70, 100 and 160 μm. The fractional dust luminosity is L_dust/L _∗  = 5.7 × 10^-5 for this warm belt shown in yellow in Fig. 7. Such a belt is comparable to the warm disk around the K0 star HD 69830 (Lisse et al. 2007). However, the proximity of the warm dust to the known planets suggests that it could be dynamically unstable (e.g. Moro-Martín et al. 2007). Nevertheless, definitive proofs are still missing to establish the reality of this warm belt in the GJ 581 system.

7. Brightness limit on scattered light around GJ 581

Our HST/NICMOS F110W image of GJ 581 after PSF subtraction (Fig. 8) is sensitive to a region from 4′′ radius (30 AU) to approximately 10′′ radius (62 AU) along PA = 120 degrees. We estimate the 3σ sensitivity to nebulosity in this region as Σ_F110W = 18.7 mag arcsec^-2. We used the radiative transfer code MCFOST (Pinte et al. 2006) to produce synthetic F110W debris disk images for pure astronomical silicate and pure water ice models that match the SED based on the geometry derived in Table 3, using the standard F110W NICMOS throughput. The maximum surface brightnesses in the 4–10′′ range at the forward scattering peak for pure water ice grains and for pure astronomical silicates are ~19.4 mag arcsec^-2 and ~21.2 mag arcsec^-2, respectively. Hence, our dust models are consistent with the nondetection of scattered light around GJ 581, but show that the disk may be detectable in deeper observations.

Fig. 8 HST/NICMOS F110W image of GJ 581 (north is up, east is left). The star was not placed behind the occulting spot available on the NIC2 camera. We subtracted the PSF using observations of GJ 250B made earlier in the same scientific program, as described in Krist et al. (1998). The circular black digital mask has 4′′ radius (30 AU) and blocks the central region where PSF subtraction artifacts are significant. Along the position angle of the disk (PA = 120 degrees) the field of view is limited to approximately 10′′ radius (62 AU). We estimate the 3σ sensitivity to nebulosity in the 4′′–10′′ radius region as ΣF110W = 18.7 mag arcsec-2. Lack of detectable scattered light at this level is consistent with the dust model derived from the far-IR PACS images.

8. Discussion

We have spatially resolved a disk around the mature M-star GJ 581 hosting four planets. This cold disk is reminiscent of the Kuiper belt in the solar system but it surrounds a low mass star (0.3 M_⊙) and has a much higher fractional dust luminosity L_dust/L_∗ of ~10^-4. It shows that debris disks can survive around M-stars beyond the first tens of Myr after the protoplanetary disk disperses, and they can be detectable although they have been elusive in searches so far.

8.1. Dust temperature in the cold disk

The factor f_T is significantly larger than unity in our analysis and indicates that the dust temperature ranges from ~50 to ~30 K from the inner to the outer radius of our modeled disk. This is about three times the black body equilibrium temperature for the dust around this low luminosity M3-type star. Values of f_T larger than unity have also been found for the debris disks around the G-type star 61 Vir (Wyatt et al. 2012) and several A-type stars (Booth et al. 2012). This is akin to disks resolved in scattered light which tend to be more extended than their sizes estimated from blackbody SED (Rodriguez & Zuckerman 2012). Values f_T > 1 are interpreted as evidence for dust grains of small sizes and/or optical properties different from blackbody spheres (Backman & Paresce 1993; Lisse et al. 2007; Bonsor & Wyatt 2010). When the SED of a debris disk is finely sampled spectroscopically from mid-IR to submm, a parameter search for composition, structure, size distribution of the grains can be conducted usefully (e.g. Lebreton et al. 2012). Such a parameter search would be degenerate for GJ 581 because of its limited photometry, and so these effects have been reduced to the single parameter f_T in our model.

Fig. 9 Ratios of radiation pressure (solid line) and stellar wind pressure (dashed line) to stellar gravity as a function of particle diameter for icy grains around GJ 581. Particles with β > 0.5 are put on hyperbolic orbits as soon as they are created and so removed from the system on orbital timescales, thus setting the blow-out limit.

8.2. Collision, Poynting-Robertson and stellar wind timescales for the GJ 581 system

In addition to gravitational forces, dust dynamics is controlled by radial forces (radiation and stellar wind pressures) and by tangential forces (Poynting-Robertson and stellar wind drags). For large dust grains, these perturbing forces act on much longer timescales than collisions, and such grains simply orbit the star until they are broken into smaller fragments in collisions with other grains. This results in a collisional cascade with a size distribution with a characteristic slope n(D) ∝ D ^− 7/2 (assuming dust grain strength is independent of size). For small dust grains, perturbing forces truncate (or at least significantly deplete) the size distribution at scales where one of the perturbing force timescales is shorter than the collisional lifetime (e.g. Wyatt et al. 2011). To ascertain the process dominating the dust removal requires a comparison of the relevant timescales.

Figure 9 shows the ratio β of the radiation pressure to stellar gravity experienced by icy dust grains of different sizes. This peaks for sizes comparable to the wavelength where the stellar spectrum peaks and is also proportional to L_∗/M_∗, where L_∗ and M_∗ are the luminosity and mass of the star (Gustafson 1994). The low luminosity of GJ 581 means that β < 0.5 for all icy grains regardless of size, and the same is true for other compositions. Since dust with β < 0.5 that is created in collisions is always placed on a bound orbit, this means that radiation pressure is not a mechanism that can be invoked to expel the dust from the system (Wyatt et al. 1999).

Figure 9 also shows the ratio β_sw of the stellar wind pressure to stellar gravity (Gustafson 1994). This depends on the stellar mass loss rate and stellar wind speed that are poorly understood and hard to measure for M-stars. Here we estimate the mass loss rate from the nondetection of X-rays from GJ 581 by ROSAT implying log L_x < 26.44 erg/s (Schmitt et al. 1995). The correlation between X-ray surface flux and mass loss rate of GKM-type stars (Wood et al. 2005) then yields the upper limit of 2 Ṁ_⊙, where the solar mass loss rate Ṁ_⊙ = 2 × 10^-14 M_⊙ yr^-1. We also consider the stellar wind speed to be ~400 km s^-1 as appropriate for GKM-type stars (Wood 2004). Finally we assumed 100% efficiency of momentum coupling between dust and the stellar wind (and used Eq. (12) of Plavchan et al. 2005). With these assumptions we found that β_sw can only be > 0.5 for dust smaller than a few nm, meaning that stellar wind pressure could only truncate the collisional cascade below the nm-scale; furthermore, stellar wind pressure would be ineffective if small grains couple inefficiently to the stellar wind (e.g. Minato et al. 2006).

Fig. 10 Dust removal timescales as a function of particle size, due to collisions (solid line), Poynting-Robertson drag (dashed line), stellar wind drag (dash-dot line), and stellar wind pressure (dotted line).

A comparison of timescales first requires an estimation of the collisional lifetime of dust grains of different sizes. Here we use Eq. (4) and the parameters for the disk found from the modeling of the combined images presented in Table 3 to derive the total mass M_tot = 2.2  ×   in M_⊕, assuming the standard size distribution between D_min = 1 μm and the diameter D_c of the largest objects and the density ρ = 1.2 g/cm³ for icy grains. The collisional lifetime is estimated using Eq. (16) of Wyatt (2008) with the additional assumptions that dust orbital eccentricities are 0.05 and that their strength is 10³ J/kg (independent of size so as to be consistent with the assumptions about the size distribution). The resulting collisional lifetime is  Myr, where D is in μm. Since we expect the cascade to extend up to sizes for which their collisional lifetime is equal to the age of the star (~5000 Myr), then as long as our assumptions apply up to large sizes we can get a rough estimate of the total mass of the disk as 0.16 M_⊕ in objects up to D_c = 0.5 km in diameter.

Figure 10 shows the timescale for dust to migrate from the inner edge of the disk at 25 AU to the star due to P-R drag. The dependence of this timescale on particle size results from a scaling ∝ 1/β which means that this has a minimum value of 60 Myr. Since this timescale is longer than the collisional lifetime at all sizes, P-R drag is not a significant loss process from the disk. Figure 10 also shows the corresponding timescale for migration due to stellar wind drag. This also includes a scaling ∝ 1/β_sw, and the efficient momentum coupling assumed here means that this timescale decreases indefinitely to smaller sizes ∝ D. As a result, stellar wind drag timescales become shorter than collisional timescales at a size of around 3 nm.

Thus, if all of the assumptions hold, we would expect the collisional cascade to extend down to 3 nm, while smaller dust is removed by stellar wind drag. However, it should be noted that there are significant uncertainties, both on the magnitude of stellar wind drag and its efficiency of coupling to small grains, and on the geometry of the dust disk which impacts the collisional lifetimes. As such this plot should be considered as representative of the kind of arguments that need to be considered when assessing the fate of material in the debris disk of GJ 581. Further study of this issue is left to later papers, but here we note that the existence (or not) of grains smaller than 1 μm is not important for the observable properties of the disk discussed in this paper, since such grains are inefficient emitters in the far-IR.

Another scenario that we have not considered in detail is that the dust is all in large mono-sized grains, in a configuration meaning that the dust collides at low enough velocities that particles bounce off each other rather than destroy each other (Heng & Tremaine 2010). Two constraints on such models are that the SED should look like a black body (since all the dust is large), and the fractional luminosity should not be large enough that collisions must necessarily occur at high velocity. Here the fractional luminosity only constrains the collision velocity at the inner edge to be > 0.3 m/s which is not sufficient to require a collisional cascade. However, although there is no evidence from our limited photometry of GJ 581 that the spectrum departs from black body shape, the resolved location of the dust shows that it is significantly hotter than black body, consistent with the presence of small grains and so incompatible with this model.

8.3. Planets and disk relationship for the GJ 581 system

First, we note that our determination of the inclination of the disk relative to the plane of the sky is 30° < i < 70° (face on disk is i = 0°). This is mostly constrained by the 160 μm image and is fairly insensitive to the masks used for the background sources. If the disk mid-plane and the orbits of the planets are coplanar, this range of inclination makes the masses of the planets of GJ 581 no more than ~ 1.6 times their measured minimum masses by radial velocity and, interestingly, ensures the long-term stability of the orbits in this system as shown in dynamical studies by Beust et al. (2008) and Mayor et al. (2009).

In our DEBRIS sample of 89 M-stars, there are only three M-stars with known planets overall (GJ 876, GJ 832 and GJ 581). GJ 581 hosts low mass planets and now has a detected disk, while GJ 876 and GJ 832 host Jupiter mass planets and have no detected disk brighter than the fractional dust luminosity 10^-5 in our survey as we shall present in a future study (Matthews et al., in prep.). Hence, using these three stars as a sample, the outcome is one disk for one low mass planet system (1/1) and no disk for two high mass planet systems (0/2). Although this is small number statistics, we note that it is suggestive that the correlation between low-mass planets and debris disks recently found for G-stars by Wyatt et al. (2012) also applies to M-stars. It is also intriguing that the only debris disk confidently detected in our current analysis surrounds the one star in the sample that hosts low-mass planets. We note that simulations by Raymond et al. (2011, 2012) suggest that a correlation might exist between low mass planets and debris as a result of planet formation processes. Note that the star AU Mic does not fall in the DEBRIS sample of the nearest M-stars, and so is not included in the statistics above; this young star (12 Myr) has a bright disk, but no known planets, although radial velocity measurements toward AU Mic are insensitive to planets with masses lower than a few Jupiters even for short orbital periods because of its high chromospheric activity (see GJ 803 in Fig. 19 of Bonfils et al. 2011).

Current programs show that a large fraction of M-stars are orbited by low-mass planets. The radial velocity survey of 102 M-stars conducted by Bonfils et al. (2011) yields the high occurence of % for low mass planets (2–10 M_⊕) around M-stars, unlike the low occurence of giant planets of ~ 2%, for orbital periods under 100 days. Transit observations in the Kepler field show that small candidate planets (2–4 R_⊕) with P < 50 days are found around 25 ± 10% of the M-stars (T_eff = 3600 − 4100 K), seven times more frequently than around the hottest stars (6600–7100 K). There is no such a dependence for larger planets (4–32 R_⊕) with P < 50 days, found uniformly around 2 ± 1% of the stars across all spectral types in the Kepler field (Howard et al. 2012). Hence, if there is a correlation between the debris disks and low-mass planets of M-stars at a level similar to that found for G stars (4/6 nearby G-stars with detectable low-mass planets have detectable disks, Wyatt et al. (2012), the high fraction of M-stars with low-mass planets would explain the detection of the disk around GJ 581 but would imply also that more disks were expected to be detected in our DEBRIS M-star sample. We defer this discussion to a paper that describes those observations in more detail, since not all observations are sensitive to disks at the same level. However, GJ 581 is at the median distance of our DEBRIS M-star sample that was observed to uniform depth, so it could simply be the brightest because of some intrinsic properties. The explanation could also be related to its multiple planetary system.

Secular perturbation theory has been applied to planetesimals in debris disks perturbed by planets in Wyatt et al. (1999) and Mustill & Wyatt (2009). The outermost planet GJ 581d (5.4 M_⊕, a_pl = 0.22 AU, and e_pl = 0.25, the highest eccentricity in the system) cannot stir the disk at a = 25 AU because the timescale for orbital crossing of planetesimals is much longer than the stellar age (Eq. (15) in Mustill & Wyatt 2009). However, a hypothetical outer planet, for example a Neptune mass planet (17_⊕) at 5 AU with a moderate orbital eccentricity of 0.2, can stir the disk at a = 25 AU in much less than the age of the system, and trigger destructive collisions of 0.5 km-sized bodies (Eq. (27) in Mustill & Wyatt 2009) to feed a collisional cascade. The most recent detection limit on m sin i from radial velocity measurements of GJ 581 over 3.3 years indicates that such a planet would not have been detected (Bonfils et al. 2011, Fig. 13), and there is a large region of parameter space of m sin i vs a over which a planet could both stir the disk and have eluded detection in radial velocity measurements.

Alternatively, mild collisions between planetesimals in a weakly excited disk could eventually form a Pluto-sized body, which in turn stirs the disk so that it produces dust (Kenyon & Bromley 2008). This self-stirring scenario is plausible since the timescale required for the formation of a Pluto-sized body at 50 AU is comparable to the age of GJ 581, even when the surface density of solids is ten times smaller than the minumim-mass solar nebula around an M3-type star (Eq. (41) Kenyon & Bromley 2008, with x_m = 0.1).

9. Conclusion

We have spatially resolved a debris disk around the M-star GJ 581 with Herschel PACS images at 70, 100 and 160 μm and modeled these observations. This is the second spatially resolved debris disk found around an M-star after AU Mic, but, in contrast, GJ 581 is much older and is X-ray quiet. Our best fit model is a disk, extending radially from 25 ± 12 AU to more than 60 AU. Such a cold disk is reminiscent of the Kuiper belt but it surrounds a low mass star (0.3 M_⊙) and its fractional dust luminosity L_dust/L _∗  of ~ 10^-4 is much higher. Also, in our best fit model, the dust temperature is found to be significantly higher than the blackbody equilibrium temperature indicating that small grains are abundant. Finally, the inclination limits of the disk make the masses of the planets small enough to ensure the long-term stability of the system according to dynamical simulations by Beust et al. (2008) and Mayor et al. (2009).

This disk complements our view of this remarkable system known to host at least four low mass, close-in planets. These planets cannot perturb sufficiently the modeled cold disk to trigger destructive collisions between planetesimals over the age of the star, but a hypothetical outer planet, for example a Neptune mass planet with an orbital radius of 5 AU and a moderate eccentricity, could replenish the system with dust. Alternatively, the self-stirring mechanism could operate for this old star causing sufficient dynamical excitation to produce the observed dust.

It is intriguing that, in our current analysis of the DEBRIS sample of 89 M-stars, the only debris disk confidently detected around a mature M-star also happens to be around the only star known to have low mass planets. This could mean that the correlation between low-mass planets and debris disks recently found for G-stars by Wyatt et al. (2012) also applies to M-stars. Then, the high fraction (~25%) of M-stars known to host low mass planets in the radial velocity and Kepler observations should make debris disks relatively common around them. If these disks have not been detected yet, it may be because searches have simply not been deep enough, or because the disk around GJ 581 is the brightest owing to some intrinsic properties; for example hosting a multiple planetary system.

Future studies and complementary observations of GJ 581 at higher angular resolution will enhance further our knowledge of this remarkable system around an M-star.

Acknowledgments

The Herschel spacecraft was designed, built, tested, and launched under a contract to ESA managed by the Herschel/Planck Project team by an industrial consortium under the overall responsibility of the prime contractor Thales Alenia Space (Cannes), and including Astrium (Friedrichshafen) responsible for the payload module and for system testing at spacecraft level, Thales Alenia Space (Turin) responsible for the service module, and Astrium (Toulouse) responsible for the telescope, with in excess of a hundred subcontractors. We thank Ben Zuckerman for comments on a draft of this article. J.F.L. gratefully acknowledges the financial support of Centre National d’Etudes Spatiales (CNES). J.H. gratefully acknowledges the financial support of the Australian government through ARC Grant DP0774000. M.B. is funded through a Space Science Enhancement Program grant from the Canadian Space Agency.

References

All Tables

All Figures

Fig. 1 PACS images of GJ 581 cropped to ± 50 arcsec from the star, at 70, 100 and 160 μm from top to bottom, respectively. In the left-hand column, the raw images show that the main emission is centrally located about the star position (image center) and that there are several point-sources in the field, barely detected at 70 μm, significantly at 100 μm, and more prominently at 160 μm, as expected for submillimeter background galaxies. In Sect. 4, we show that the main emission is extended and centered on the star position, as expected for a debris disk, and mingles with a background source ~11 arcsec toward the northwest. The panels in the middle column are the photosphere-subtracted images. The panels in the right-hand column show the best subtraction (lowest residuals) of a two-point source model, which assumes that there is no debris disk around GJ 581 but an extra background source located exactly behind the star in addition to the N W source. This model is rejected because of the systematic residuals left, indicative of an extended structure, especially at 70 and 100 μm. At each wavelength, the contours levels of the three images are the same and correspond to 1, 2, 3, 9, 15σ0 (σ0 = 0.0135 mJy/1′′ pixel) at 70 μm, 1, 2, 3, 6, 9, 12, 15σ0 (σ0 = 0.0094 mJy/1′′ pixel) at 100 μm, and 1, 2, 3, 5, 7, 9, 11σ0 (σ0 = 0.0251 mJy/2′′ pixel) at 160 μm. The coordinates of the image center provided in the labels correspond to the star position at epoch of observation (July 29th–August 1st 2011). The hatched circles are the beam FWHMs: 5.6, 6.8, and 11.4 arcsec at 70, 100, and 160 μm, respectively.

In the text

Fig. 2 Radial profiles of the mean brightness of the photosphere-subtracted images. Each point of these curves was computed as the mean of the emission in an elliptical annulus centered on the peak of the main emission close to the image center. They correspond to an axisymmetric disk model inclined to 40° relative to the plane of the sky (see text). They are one pixel wide, i.e. 1 arcsec in the 70 and 100 μm images and 2 arcsec in the 160 μm image. The Gaussian with the FWHM of the beam is the profile expected for a hypothetical point source in the background aligned by chance with the star and shown for comparison.

In the text

	Fig. 3 SPIRE image of GJ581 at 250 μm. Pixel size is 6′′ and the contour levels are 1 and 2σ with σ = 6.1 mJy/18.2′′ beam. The star symbol is the position of GJ 581 at the date of observation.
In the text

Fig. 4 Maps of residuals for the best fit models of the disk to the PACS images at 70, 100, and 160 μm from left to right. The post-fit residuals are between ± 3.5σ0, coded as black, blue, orange, yellow and white over this range (σ0 is the same noise rms as used for contours in Fig. 1 but the color scale is not the same as in Fig. 1). In zooming the electronic version, the contours apparent in these maps are −3, −2, −1, 1, 2, 3 σ0 (dashed contours are negative levels). The two background sources are masked in the 100 and 160 μm images as the dotted squares show.

In the text

	Fig. 5 Map of the reduced showing the correlation between the temperature factor fT and the total cross-sectional area A of the dust. The minimum is the white region in this map.
In the text

Fig. 6 SED of the cold disk model. The left-hand figure shows the best fit to the three PACS images only (, Table 3). The right-hand figure shows the best fit to both the three PACS images and the IRS spectrum . The modeled cold dust emission is the blue curve, the Next Gen stellar atmosphere spectrum is the gray curve, and their sum is the green curve. The IRS spectrum is in red. The upper inset zooms on the IRS wavelengths and displays spectra as Sν × (λ/24)2 on a linear scale for clarity. The lower inset zooms on the three PACS bands. In the left-hand figure, the best fit model satisfactorily fits the PACS data as shown by the lower inset but misses the IRS data as shown by the upper inset. In the right-hand figure, the best fit model partially misses the PACS data but satisfactorily fits the IRS data.

In the text

Fig. 7 SED of the best fits of the cold disk and warm belt models. The green curve corresponds to the warm and cold dust emissions added to the Next Gen stellar atmosphere spectrum (gray). IRS spectrum is in red. We also show separately the best fit of the warm belt emission (rw = 0.2 AU and mw = 2.8 × 10-6 M$) (yellow) and the best fit of the cold disk emission (parameters of the combined fit are in Table 3) (blue). Insets and photometric data points are the same as in the legend of Fig. 6.

In the text

Fig. 8 HST/NICMOS F110W image of GJ 581 (north is up, east is left). The star was not placed behind the occulting spot available on the NIC2 camera. We subtracted the PSF using observations of GJ 250B made earlier in the same scientific program, as described in Krist et al. (1998). The circular black digital mask has 4′′ radius (30 AU) and blocks the central region where PSF subtraction artifacts are significant. Along the position angle of the disk (PA = 120 degrees) the field of view is limited to approximately 10′′ radius (62 AU). We estimate the 3σ sensitivity to nebulosity in the 4′′–10′′ radius region as ΣF110W = 18.7 mag arcsec-2. Lack of detectable scattered light at this level is consistent with the dust model derived from the far-IR PACS images.

In the text

	Fig. 9 Ratios of radiation pressure (solid line) and stellar wind pressure (dashed line) to stellar gravity as a function of particle diameter for icy grains around GJ 581. Particles with β > 0.5 are put on hyperbolic orbits as soon as they are created and so removed from the system on orbital timescales, thus setting the blow-out limit.
In the text

	Fig. 10 Dust removal timescales as a function of particle size, due to collisions (solid line), Poynting-Robertson drag (dashed line), stellar wind drag (dash-dot line), and stellar wind pressure (dotted line).
In the text