Multiwavelength modelling of the β Leo debris disc: one, two or three planetesimal populations?*

Abstract

1 INTRODUCTION

Debris discs are distributions of dust and planetesimals with radii of 1–1000 au around main-sequence stars (see Wyatt 2008 for a recent review). The dust grains in these discs are small, and so they cannot be primordial as they would have been blown out of the system by radiation pressure on time-scales shorter than the stellar age. The dust must be continually replenished from a population of colliding planetesimals, thought to contain bodies up to ∼1 km in size (Wyatt & Dent 2002). However, as debris discs age the planetesimal population is ground down, so discs become less massive and fainter (Dominik & Decin 2003). At far-infrared and submillimetre wavelengths debris discs are optically thin, the disc-to-star contrast is favourable, and these wavelengths are sensitive to the large (up to ∼1 mm) grains that dominate the dust mass in debris discs.

The dust morphology of a debris disc can be shaped by planets in the system, so resolved images of discs help constrain models of structure and evolution of planetary systems. Resolved images can indicate that infrared excess is being produced by multiple dust populations and can also break the degeneracy between the radial location of the dust and its temperature.

The DEBRIS (Disc Emission via a Bias-free Reconnaissance in the Infrared/Sub-mm) survey (Phillips et al. 2010; Matthews et al., in preparation) is an Open Time Key Program on the Herschel Space Observatory which uses PACS (Photodetector Array Camera and Spectrometer, Poglitsch et al. 2010) and SPIRE (Spectral and Photometric Imaging REceiver, Griffin et al. 2010) to detect, resolve and characterize debris discs around a volume-limited sample of 446 A through M stars (Matthews et al. 2010). β Leo (HD 102647) was observed at 100 and 160 μm with PACS as a part of the DEBRIS survey Science Demonstration Phase (Matthews et al. 2010).

β Leo (A3V, L_* = 14.0 L_⊙) is a δ Scuti-type star at a distance of 11.1 pc. The infrared excess around this main-sequence star was first discovered by IRAS (Infrared Astronomical Satellite, Oudmaijer et al. 1992), then confirmed with Spitzer (Su et al. 2006). The excess was unresolved in mid-infrared (mid-IR) imaging (Jayawardhana et al. 2001; Akeson et al. 2009, section 2.2) although differences between the IRAS and ISO (Infrared Space Observatory) fluxes led Laureijs et al. (2002) to suggest that the disc emission may be somewhat extended in the ISO beam (52 arcsec aperture). Chen et al. (2006) obtained a Spitzer IRS (Spitzer Infrared Spectrograph) spectrum of β Leo (see Section 3.2) and found a featureless continuum spectrum consistent with dust at ∼120 K located at 19 au from the star. A very hot excess has also been partially resolved using infrared interferometry with the FLUOR (Fiber Linked Unit for Recombination) instrument at the CHARA (Center for High Angular Resolution Astronomy) array at 2 μm and BLINC (Bracewell Infrared Nulling Cryostat) at 10 μm (Akeson et al. 2009; Stock et al. 2010; see Section 2.8 for more details).

β Leo is thought to be a member of the IC2391 moving group, giving an age of 45 Myr (Nakajima, Morino & Fukagawa 2010). The age of this source determined from isochrone fitting is 50–331 Myr (Lachaume et al. 1999; Song et al. 2001). Di Folco et al. (2004) derive an age from measuring the stellar radius and suggest 100 Myr. We assume an age of 45 Myr for this paper. β Leo does not have any known companions within a few arcsecs of the star. The Washington Double Star (WDS) catalogue lists three companions with common proper motion for β Leo. These stars are located at 40 to 240 arcsec from the primary with V magnitude differences of 6.3–13 (Worley & Douglass 1997); these stars are not, however, physically associated with β Leo (Phillips et al. 2010). β Leo was also included in a high-precision radial velocity survey of early-type dwarfs with HARPS (Lagrange et al. 2009) for planetary or brown dwarf mass candidates and was found to have no companions with mass > 4.2M_jup with periods <10 d with 99.7 per cent probability.

This paper uses multiwavelength modelling of the β Leo debris disc including PACS and SPIRE observations to present a self-consistent explanation of the system. In Section 2 we present the observations, including 100- and 160-μm PACS imaging (previously published in Matthews et al. 2010) and archive Gemini MICHELLE 12- and 18-μm imaging. In Section 3 we confront the observations with models to determine the disc parameters. In Section 4 we discuss the implications of the inferred structure for the status of planet formation in this system.

2 OBSERVATIONS

2.1 Fitting the stellar photosphere

Being the fifth nearest A-type star, β Leo is bright and saturated in modern surveys such as 2MASS. In addition, the infrared excess has a hot component, which contributes to bands typically used for modelling stellar photospheres (Akeson et al. 2009; Stock et al. 2010). We use the equivalent 2MASS K_s magnitude of 1.93 from Stock et al. (2010), and include mean UBV and Hipparcos H_p photometry (ESA 1997; Mermilliod 2006). The best-fitting Castelli & Kurucz (2003) model, found by a χ² minimization method, has T_eff = 8660 K, L_★ = 14 L_⊙ and R_★ = 1.66 R_⊙. In the Multiband Imaging Photometer for Spitzer Space Telescope (MIPS) 24-μm band we predict 1171 ± 15 mJy, 1 per cent lower, but consistent with the value in Stock et al. (2010). The predicted stellar flux densities at the PACS effective wavelengths of 100 and 160 μm are 64 ± 1 and 26 ± 0.5 mJy, respectively. These are 9 per cent higher than the values presented in Matthews et al. (2010) due to improvements in the fitting procedure and data used.

2.2 Herschel PACS observations

Observations of the β Leo disc at 100 and 160 μm were taken with the ESA Herschel Space Observatory (Pilbratt et al. 2010) using the PACS (Poglitsch et al. 2010) instrument in photometry mode as a part of the Science Demonstration Phase observations for the DEBRIS survey. DEBRIS is a flux-limited survey which observes each target to a uniform depth of 1.5 mJy beam⁻¹ at 100 μm. The images were first presented in Matthews et al. (2010). β Leo was observed in a small scan-map mode (see PACS Observers’ Manual1 for details). Scan map observations had eight repeats in a single scan direction at a scan rate of 20 arcsec s⁻¹. Four 3-arcmin scan legs were performed per map with a 2-arcsec separation between legs. The total observing time was 1220 s.

These data were reduced using the Herschel Interactive Processing Environment (hipe, Ott 2010). Maps were obtained via the default PACS-naive map-making method photProject in hipe. The data were pre-filtered to remove low-frequency (1/f) noise using a boxcar filter with a width of 98 arcsec. All the bright sources in the map were masked prior to filtering to avoid ringing-type artefacts.

Photometry on the maps of β Leo (Fig. 1) using a 13-arcsec radius circular aperture centred on the emission peaks yielded fluxes of 480±30 mJy at 100 μm and 215 ± 32 mJy at 160 μm, respectively. These fluxes greatly exceed the rms noise levels in the maps of 1.4 mJy beam⁻¹ at 100 μm and 3.1 mJy beam⁻¹ at 160 μm. These values have been colour corrected using the values from the PACS Observers Manual assuming a temperature of 120 K corresponding to the temperature of the disc estimated from a blackbody fit to the spectral energy distribution (SED). The quoted errors on these fluxes do not include calibration uncertainties which are estimated to be 10 and 20 per cent at 100 and 160 μm (Poglitsch et al. 2010). These were combined in quadrature with statistical uncertainties from the rms levels in the maps. The stellar flux for fitting the photosphere is 64 and 26 mJy at 100 and 160 μm, respectively. This gives excess fluxes of 416 ± 30 mJy at 100 μm and 189 ± 32 mJy at 160 μm. These values differ from those in Matthews et al. (2010), which quoted 500 ± 50 mJy at 100 μm and 230 ± 48 mJy at 160 μm in a 20-arcsec radius aperture. This difference appears because the data presented in this paper were re-reduced using a different filter scale and a more recent version of hipe (built by the 4.2.0 developers) and a smaller aperture was used.

Figure 1

Images of the 100 (left) and 160 μm (right) emission for β Leo taken with PACS. The images have had a PSF scaled to the stellar flux (64 and 26 mJy at 100 and 160 μm, respectively) subtracted, and hence are maps of the excess emission. The pixel scale in 1 arcsec pixel⁻¹ at 100 μm and 2 arcsec pixel⁻¹ at 160 μm.

Open in new tabDownload slide

The maps of β Leo (Fig. 1) appear extended compared to the beam (indicated by the circle in the top-left corner of the images). At 100 μm the nominal beam size is 6.6 × 6.9 arcsec². At 160 μm it is 10.7 × 12.1 arcsec². We have obtained two bright point source images (Vesta and α Boo) taken in the same observing mode as the β Leo data to serve as point spread function (PSF) references. These have beam sizes of 6.6 × 6.9 arcsec² at 100 μm and 10.5 × 11.5 arcsec² at 160 μm for Vesta and 6.5 × 6.8 arcsec² at 100 μm and 10.2 × 11.7 arcsec² at 160 μm for α Boo. Fitting a Gaussian to the β Leo image gives a full width at half-maximum (FWHM) at 100 μm of 9.2 ± 0.1×10.4 ± 0.1 arcsec². At 160 μm the FWHM of a Gaussian fitted to the image is 13.6 ± 0.2 × 12.0 ± 0.2 arcsec². Fig. 1 shows the β Leo images after subtraction of the α Boo PSF scaled to the appropriate stellar flux. The PSF has an asymmetric three-lobed structure and, to ensure an accurate subtraction, was rotated to the same spacecraft angle at which the β Leo data were taken.

β Leo appears extended in all directions with respect to the PSF, with the major axis of the disc at a position angle of 125°± 15°E of N, which is the mean of the position angles of the Gaussians fitted at 100 [position angle (PA): 118°] and 160 μm (PA: 132°.5). The ratios of the semiminor to semimajor axis of the Gaussians (0.84 at 100 μm and 0.88 at 160 μm) give a mean inclination of 57°± 7° from edge on. The stellar subtraction was also done using the Vesta PSF to examine the effect of PSF variation on the residuals, and resulted in no significant change in the width of a Gaussian fitted to the disc, and no significant change in the flux in a 13-arcsec radius aperture.

2.3 MICHELLE mid-infrared observations

Mid-IR observations of β Leo taken with MICHELLE on Gemini North for β Leo were retrieved from the Gemini Science Archive [programme names: GN-2007A-C-10 (PI: Beichman) and GN-2006A-Q-10 (PI: Moerchen)]. These data are detailed in Table 1 and were taken with filters Qa (⁠ = 18.1 μm, Δλ = 1.51 μm) and N′ (⁠ = 11.3 μm, Δλ = 1.07 μm). MICHELLE has a pixel scale of 0.1005 arcsec pixel⁻¹. Imaging was taken using an ABBA chop-nod sequence with a PA of 30° for all observations. The data taken under GN-2006A-Q-10 were divided into three observation groups, with observations of a standard star from Cohen et al. (1999), HD 98118 (M0III, F_{18 m} = 4.2 Jy), taken before and after the observations of β Leo to serve as a standard star for flux calibration and to monitor the PSF. The data taken under GN-2007A-C-10 were previously published by Akeson et al. (2009). Observations of HD 109511 (K0, F_{18 m} = 1.4 Jy) preceded and followed the observations of β Leo as a standard star and PSF reference. Combining all observations, the total on-source time at 18 μm is 1751 s and at 12 μm it is 1175 s.

The data were reduced using custom routines described in Smith, Wyatt & Dent (2008). The data reduction involved determining a gain map for each observation using the mean values of each frame to construct a map of pixel responsivity. The on-source pixels were masked during this process, making this equivalent to a sky flat-field frame. A DC offset was then determined by calculating the mean pixel values in every row and every column, again masking pixels where there was source emission present. This was then subtracted from the final image to ensure a flat background. Pixels that showed high or low gain in comparison with the median response throughout the observation were masked off. To avoid errors in co-adding the data which could arise from misalignment of the images, we fitted a Gaussian with a subpixel centroid to accurately determine the centre of the image and so the position of the star. The rebinning was done using bilinear interpolation across the array.

The observations were divided into four groups (as shown in Table 1) such that flux calibration for each group was done with standards observed at a similar airmass to β Leo. The calibration levels were compared with the airmass for the standard star observations and no correlation was found, so no extinction correction was applied to the calibration factors.

For the images of β Leo in each group, a calibration factor was determined using a co-add of the two standard star observations in that group. The average calibration error for the three 18-μm groups was 8 per cent, but as seen in Table 1 the large calibration uncertainties lead to a wide range of fluxes among the three groups. The 12-μm group had a calibration error of 11 per cent. These centred, flux-calibrated images from each group were then used to produce the final images of β Leo and the standard star. The first standard (Std1) and image (Im1) in Group 1 (see Table 1) at 18 μm and Std12.1 and Im12.1 at 12 μm showed elongation in the telescope chop-nod direction, so were not used in the final co-added images, but were used when calculating the flux calibration factors.

Photometry was performed on the final co-added images of β Leo using a 1.0-arcsec radius circular aperture centred on the star. The stellar flux is expected to be 5297 and 2020 mJy at 11.3 and 18.1 μm, respectively, from fitting a model spectrum as described in Section 4. Our photometry yields fluxes of 5822 ± 476 mJy at 12 μm and 2360 ± 312 mJy at 18.6 μm including both calibration and photometric errors, giving an excess of 525 ± 476 mJy at 12 μm and 340 ± 312 mJy at 18 μm. The statistical noise was determined using an annulus with an inner radius of 7 arcsec and an outer radius of 8 arcsec centred on the star, resulting in an error averaged over all the groups of 0.31 mJy arcsec⁻² at 12 μm and 0.41 mJy arcsec⁻² at 18 μm. The IRS spectrum presented in Stock et al. (2010) gives an excess above the photosphere of 61 ± 103 mJy at 12 μm and 256 ± 53 mJy at 18 μm. The photospheric fit used here is consistent with that used in Stock et al. (2010). The MICHELLE photometry presented here therefore agrees with the IRS results, but the calibration errors are too large to detect the excess. To assess whether these mid-IR images have resolved the disc, line cuts were taken along the PA of the extension (125°E of N) seen in the PACS images (125°) for both β Leo and the standard star. These are shown in Fig. 2. β Leo shows no extension when compared with the PSF, and the excesses at both 12 and 18 μm are unresolved. This is consistent with the previous results of Akeson et al. (2009) and Moerchen et al. (2011). None of the groups shows significant extension. We also tried subtracting the standard star image scaled to the peak flux of the β Leo images, but no significant structure remains.

Figure 2

The profile of the line cuts through the total co-added mid-IR images at 125°, the position angle of the extension seen in the PACS images. The dotted lines are the final co-added image of β Leo in each group and the solid lines are the final co-added standard star in each group. The 12-μm line cut is the top image and the 18-μm one is the bottom image. β Leo shows no significant extension with respect to the PSF, indicating the excess at this wavelength that is unresolved. In fact, the observations appear narrower than the PSF, but this is not significant and is due to variation in observing conditions.

Open in new tabDownload slide

As the location of the disc is known from the resolved PACS images with peak emission at ∼5 arcsec radius, we tried convolving the co-added images with a series of Gaussians with FWHMs from 1 to 20 arcsec to find any large-scale, low surface brightness features in the outer regions of the MICHELLE images. The MICHELLE field of view is 32 × 24 arcsec², but the chop throw of 15-arcsec limits the usable region. We find no coherent features at the expected location of the debris disc in the convolved image. The surface brightness in an annulus from 0.5 to 1 arcsec in the MICHELLE 18-μm imaging of β Leo is 0.41 ± 0.34 mJy pixel⁻¹. However, the brightness in the same annulus of the PSF scaled to the flux of β Leo is 0.31 mJy pixel⁻¹.

2.4 SCUBA 2

The SCUBA 2 instrument (Sub-mm Common User Bolometer Array 2) (Holland et al. 2006) operating on the JCMT ((James Clerk Maxwell Telescope) was used to obtain fully sampled maps of β Leo, and the surrounding region at 450 and 850 μm. These data were obtained in the nights of 2010 February 16 and 17 as a part of the shared-risk observing phase of instrument commissioning. The final maps contain data from six observations obtained in each night, giving a median integration time of 3 h pixel⁻¹ in the central 1 arcmin², and a total observing time of 6 h. The ‘daisy’ pattern observing mode was used at a scanning rate of 120 arcsec s⁻¹ (Kackley et al. 2010), giving a map with a usable area of ∼9 arcmin².

The data were processed and calibrated using the smurf package in the Namaka Starlink release (Dempsey et al. 2010; Jenness et al. 2011). The data were high-pass filtered to mitigate the 1/f noise present in the data, with filter parameters set so as to retain features smaller than ∼120 arcsec. Maps were obtained at 450 and 850 μm and the data were reduced using the latest version of the smurf pipeline. β Leo was not detected in either band. Photometry on the 450-μm map gives a 3σ upper limit of 50 mJy beam⁻¹. At 850 μm the 3σ upper limit is 6 mJy beam⁻¹. The expected stellar fluxes in these wavebands are 3.1 and 0.8 mJy, respectively. Holmes et al. (2003) found an upper limit of 20 mJy at 870 μm, so these upper limits significantly improve the constraints on the SED in the sub-mm and were used in the SED fitting described in Section 3.

2.5 Spitzer MIPS observations

In addition to the photometric points at 24 and 70 μm on the SED shown in Fig. 1, we also include the partially resolved MIPS image of the disc at 24 μm, presented in Stock et al. (2010). Stock et al. (2010) observed that the 24-μm photometric point obtained using an aperture of 14.94 arcsec gives a total integrated flux of 1623 ± 33 mJy which is 2.5 per cent higher than the flux obtained with a 6.23-arcsec aperture, leading to the suggestion that the disc may be marginally resolved at this wavelength. The photometric value from the larger aperture is used for SED fitting. The MIPS 24-μm image is shown in the first panel of Fig. 3. Compared to the PSF shown in the right-hand panel of Fig. 3, the first dark Airy ring (with a radius between 5.5 and 7.5 arcsec) appears more filled in. Subtraction of the PSF scaled to the expected photospheric flux of β Leo gives residuals with an FWHM of 6.88 × 6.61 arcsec² at a PA of 118° which is larger than the FWHM of the PSF (5.69 × 5.53 arcsec²). This is again consistent with a disc that is marginally extended with an inclination of 60°± 10° from edge on and a PA of 118°± 3°. Comparing β Leo to the PSF observed for a star with hot dust at <2 au (i.e. similar to the interferometric detection) of the same spectral type (ζ Lep – 5.61×5.55 arcsec², Su et al. 2008) and to a theoretical PSF (generated using TinyTim2) also suggests that the disc is marginally resolved. If the 24-μm observations are deconvolved from the PSF (i.e. the FWHM of the PSF is subtracted in quadrature from the FWHM of the observations), it gives 3.9 arcsec for the FWHM of the residuals. The PA and inclination of these residuals are consistent with the disc resolved using PACS. The 70-μm Spitzer MIPS observations are consistent with a point source. We include the 24-μm image when considering the observational constraints for a model of the β Leo debris disc.

Figure 3

Spitzer MIPS observations of β Leo at 24 μm. The disc is detected and partially resolved compared to the PSF. The disc image is shown in the left-hand panel, the PSF is the centre panel and the lower panel shows the surface brightness profiles of the observations (solid line), and the PSF (dashed line) shows the filling in of the first dark Airy ring between 5.5 and 7.5 arcsec.

Open in new tabDownload slide

2.6 Optical observations

β Leo was imaged using the Hubble Space Telescope Advanced Camera for Surveys High Resolution Channel (HRC) coronagraph on 2004 March 25 (GO-9475, PI: Kalas). We occulted β Leo with the 1.8-arcsec diameter occulting spot approximately centred on the 1024×1024 pixel camera. Seven 140-s exposures were acquired using the F606W filter (⁠ nm, δλ = 234 nm). The A1V star HD 95418 was observed before β Leo to serve as a reference for the stellar PSF. Cosmic rays were filtered by taking the median combination of the seven β Leo exposures, as well as the six 180-s exposures on HD 95418. It should be noted that HD 95418 β UMa has an IR excess indicative of a debris disc but this is not detected in scattered light. The HD 95418 observation was registered and scaled by a factor of 1.22 to subtract the β Leo PSF. The resulting, PSF-subtracted image was then corrected for geometric distortion giving a pixel scale of 0.025 arcsec pixel⁻¹. Fig. 4 shows the resulting optical image of β Leo. We find a residual halo of light in an annulus 2.5–4.0 arcsec from the star. This halo could plausibly originate from either dust-scattered light or an imperfect PSF subtraction. The latter effect results from thermal ‘breathing’ of the telescope between the observations of the two stars, as well as a small colour mismatch between β Leo and HD 95418. The limits on surface brightness from these observations were used as constraints for the modelling described in Section 3.

Figure 4

Optical HST observations of β Leo after PSF subtraction. North is up, east is left. The green circle has diameter 8 arcsec or 88 au. The residual halo of light discussed in the text is contained interior to the green circle. Saturation columns are evident to the upper right and lower left of the star. The bright circular structure above the saturation columns is a different occulting spot in the HRC focal plane.

Open in new tabDownload slide

2.7 Herschel SPIRE observations

Photometric observations of β Leo were obtained with the SPIRE instrument on Herschel (Griffin et al. 2010), providing maps in wavebands centred at 250, 350 and 500 μm. These data were observed on the 2011 November 23 (Herschel operational day 558) using the ‘small map’ observing mode. Five repeat maps were performed, effectively giving confusion noise limited maps in all bands. The data were reduced using hipe version 6.0, build number 1985, and the standard pipeline script was used to perform the reduction. The images are shown in Fig. 5. Photometry was performed on the images. At 250 μm a 14-arcsec radius circular aperture (beam size 18.9 × 17.6 arcsec²)3 gives 51 ± 12 mJy. At 350 μm a 30-arcsec aperture (beam size 25.6 × 24.2 arcsec²) gives 18 ± 7 mJy which gives a 3σ upper limit of 39 mJy. At 500 μm a 40-arcsec aperture (beam size 38.0 × 34.6 arcsec²) gives 2.8 ± 4.1 mJy, which corresponds to a 3σ upper limit of 15.1 mJy.

Figure 5

SPIRE observations of β Leo at 250 μm. The disc is detected but not resolved at 250 μm. At 350 there is a 2σ detection, giving a 3σ upper limit on the disc and star emission of 21 mJy and the non-detection at 500 μm also gives a 3σ upper limit of 12.3 mJy, which is consistent with the fluxes expected of the SED modelling.

Open in new tabDownload slide

To assess whether the 250-μm image had resolved the disc, we fit the image with a 2D Gaussian, which has an FWHM of 22.5 ± 1.2×18.7 ± 1.1 arcsec² with a PA of 118°, consistent with the PA of the resolved disc (125°± 15°) from the PACS images. Comparison of the disc size with that of the SPIRE PSF (Neptune) found no significant extension at 250 μm.

2.8 Interferometry

Akeson et al. (2009) present interferometric observations of β Leo at 2 μm using the FLUOR instrument at the CHARA interferometer. These data show a short-baseline visibility deficit, indicating that the source is somewhat extended. Akeson et al. (2009) suggest a model in which the inferred few per cent near-infrared (NIR) excess is due to a population of dust grains within the field of view (<4.6 au), comprised of a population of blackbody grains that extends inward to the sublimation radius (0.12 au for β Leo for a sublimation temperature of 1600 K), which maximizes the 2 μm emission. This places a lower limit on the fractional dust luminosity (defined as ⁠, i.e. a measure of the IR excess for each component) for this hot population of 8.0 ± 1.1 × 10⁻⁵, which is larger than the value of 2.7 × 10⁻⁵ for the IRS emission (Chen et al. 2006). The strongest spatial constraint from the interferometry is that there is dust with a fractional luminosity of 8 × 10⁻⁵ within 4.6 au.

Stock et al. (2010) present new 10-μm nulling interferometry data from the Keck Interferometer Nuller (KIN), which suggest that there is no significant resolved emission at 10 μm in the KIN beam (FWHM of 0.50 × 0.44 arcsec², Colavita et al. 2009). Limited by its small field of view (0.6 arcsec), these observations are unable to detect extended structures beyond a radius of 3 au at the distance of β Leo (11 pc). Moreover, due to its complex transmission function, KIN observations are most sensitive to extended structures with radii from 0.1 to 1 au. Stock et al. (2010) also present an N-band (8–13 μm) nulling interferometry detection with BLINC of a null of 1.74 ± 0.3 per cent, which leads to a model-based estimate of an excess flux of 250 ± 50 mJy. BLINC is sensitive to the region between 1 and 9 au (between 0.12 and 0.8 arcsec). Stock et al. (2010) quote limits for various uniform disc models: 110 mJy for a disc between 0 and 1 au, 240 mJy for a disc from 0 to 2 au and 370 mJy for a 1–2 au ring. Stock et al (2010)’s preferred model is a ring from 2 to 3 au with a flat surface density profile, and a flux of 250 ± 50 mJy at 10 μm. This flux is used as a constraint on the SED modelling.

3 MODELLING

To determine constraints on the radial distribution of the emission seen in the resolved 100- and 160-μm images, we first considered a model of the disc structure. For example, in Section 3.1 we consider the simplest possible model. This model is composed of a single axisymmetric disc component, defined by four free parameters: inner radius (r_in), outer radius (r_out), inclination (i) and surface density distribution Σ, which is assumed to have the form Σ∝r^γ. The disc opening angle (which sets the disc height) is assumed to be 5° but this parameter is unconstrained by the modelling process. The flux from annuli in the disc at different radii was determined assuming a grain composition and a size distribution that were constrained using the emission spectrum.

In order to constrain the dust location, the width of the disc, the inclination and the surface density profile, a grid of models was run. The data used in this modelling are listed in Table 2, which indicate which observations were used to constrain the SED and the images. The model images were convolved with a PSF (image of α Boo taken in the same observing mode as the observations of β Leo and reduced in the same way) and compared with the observations using the images as well as the line cuts both in the direction of extension (PA 125°) and perpendicular to the extension (PA 35°). The use of both the line cuts allowed us to constrain the inclination simultaneously with the radial morphology. Both the wavelengths were fitted simultaneously, and a joint best fit was determined by minimizing the combined reduced χ². The reduced χ² of the fit [⁠, where ν is the number of free parameters and 1 represents a good fit] to each of these six pieces of observational data was calculated and the results were then combined linearly with equal weight in order to obtain the final chi-squared (⁠⁠).

The emission spectrum of the photosphere of β Leo was discussed in Section 2.1. To calculate the emission from the dust grains in the model they were assumed to have a size distribution with n(D) ∝D^−3.5, where D is the grain diameter, which is the standard solution for a theoretical collisional cascade (see Dohnanyi 1969), which is truncated at a minimum and maximum grain size. The minimum grain size was treated as a free parameter with a range between 0.01 and 100 times D_bl, where D_bl is the largest grain size blown out of the system by radiation pressure and depends on grain composition and stellar properties. The maximum size was fixed at 1 cm since larger grains have a negligible contribution to flux from a model with this size distribution. The grains were assumed to have a silicate core (amorphous olivine) and an accreted mantle of organic refractories produced by UV photoprocessing of ice (as used in Li & Greenberg 1997; Augereau et al. 1999). A range of compositions was also tried, with amorphous silicate fractions varying from 0 to 90 per cent by volume and with porosities (i.e. vacuum fraction of grain by volume) from 0 to 95 per cent. Dielectric constants were calculated from tabulated laboratory values (Li & Greenberg 1997; Augereau et al. 1999) using Maxwell–Garnett effective medium theory. The optical properties of the grains were calculated using Mie theory, Rayleigh–Gans theory and Geometric Optics in the appropriate size regimes (Bohren & Huffman 1983). This composition was the best fit found by minimizing χ² across a grid representing the possible compositions. However, we do not give constraints on the composition because the model used contains uncertainties in the calculation of the optical properties and we do not want to overempahize this aspect of the modelling.

These models were fitted to the IRS spectrum, BLINC point, MIPS points, PACS points and SCUBA-2 upper limits. Table 2 lists the data used in the modelling and indicates if it was used to constrain the SED, images or both. As the IRS spectrum has no obvious features to fit to, to calculate we chose six windows in the spectrum, each ∼1–2 μm wide, with a constant signal-to-noise ratio in the window. This gives 11 fluxes (six IRS, one BLINC, two MIPS, two PACS) and two upper limits with which to calculate the ⁠.

The IRS spectrum is ∼2 per cent lower than the MIPS observation at 24 μm as the MIPS observation is marginally resolved, as discussed in Section 2.5. When calculating the MIPS 24-μm photometry a larger aperture was used to ensure all the flux was included, whereas the IRS spectrum assumes a point source. A synthetic model image at 24 μm was compared with the MIPS image to check that the model reproduces the marginal extension observed. The total flux was compared to both the IRS and MIPS points and the ∼2 per cent discrepancy is a small fraction of the total flux and does not affect any conclusion about the goodness of fit of the models. There is a similar discrepancy between the BLINC flux at 10 μm and the IRS spectrum, which is again small enough not to affect the conclusions of the modelling within the uncertainties.

We considered models of increasing complexity. The initial model (described in Section 3.1) is a two-component model with a hot inner disc and a cold outer disc that has been resolved with PACS. The second model has three components, a hot inner disc, a warm component and a cold outer disc, and is described in Section 3.2. The third model consists of a single eccentric planetesimal population and is described in Section 3.3.

3.1 Two-component model

The ranges of model parameters tested for the two-component model were R_in: 5–80 au (1 au intervals); R_out: 20–150 au (5 au intervals); inclination: 0°–90° (5°) intervals, where 0° is edge on); and surface density index γ: 0 to −3.0 (0.5 intervals). The values for γ were chosen to cover possibilities such as the surface density distribution expected of grains being blown out of the system by radiation pressure (γ=−1.0), and that of the minimum mass solar nebula (MMSN) (γ=−1.5). The best-fitting two-component model to the PACS 100- and 160-μm images, when considered iteratively with the SED fitting described in Section 3, was found to be a ring between 15 ± 10 and 70 ± 5 au, with an inclination of 55°± 5° from edge on, and a surface density profile index γ=−1.5 ± 0.5, a minimum grain diameter of (3 μm), a fixed maximum grain size of 1 cm, and a composition with a silicate fraction of 20 per cent, a porosity of 20 per cent with the rest of the grain composed of organic refractories, with no ice present. The hot dust is assumed to lie between 2 and 3 au with a flat surface density profile with a silicate fraction of 60 per cent, a porosity of 20 per cent and a minimum grain diameter of 0.6 μm.

The fractional luminosity of the cold component is 3 × 10⁻⁵. As there are no obvious features in the IRS spectrum, the main constraint on the composition and size distribution comes from allowing the appropriate range of temperatures, given the constraints that are there on the grains in the region 15–70 au; blackbody grains would have to be at a radius of 19 au to achieve the observed temperature of 120 K. Although our model provides a consistent fit to both image and SED, it is not expected that the composition is uniquely constrained by this process.

To evaluate the fit of the model to both the SED and the images, we use ⁠, where and is the combined reduced χ². For a model that was a perfect fit to the observations, the reduced χ² would be 1. This best-fitting model had = 4.5 with = 1.46, = 1.58 and = 1.46. For a perfect fit, would be 3. The best-fitting SED is shown in Fig. 6. The best-fitting model image is shown on the left of Fig. 7 and can be compared directly with the observed disc (Fig. 1) which has the same colour scale.

Figure 6

Spectral energy distribution (SED) of β Leo. The photosphere of β Leo is fitted with a Kurucz model profile (L_star = 14.0 L_⊙, T_star = 8660 K) fitted to the 2MASS fluxes and shown with a dark solid line. The IRS spectrum (solid line from 5 to 30 μm) of Chen et al. (2006), MIPS fluxes (crosses), BLINC 10.1 μm (upright cross) and Gemini MICHELLE fluxes in 1-arcsec radius aperture from imaging described inSection 2.2 (diamonds) and PACS fluxes in 20-arcsec aperture from imaging described in Section 1 (squares). The upper limits are 3σ from SPIRE (350, 500 μm), SCUBA-2 (450 and 850 μm), SCUBA (higher 450 and 850 μm values) and Bolocam (1.1 mm), The excess is fitted with a two-component realistic grain fit described in Section 3.2 (dark solid line – hot component is shown with a dashed line and cold component is shown with a dotted line).

Open in new tabDownload slide

Figure 7

Model PACS images for the best-fitting two-component model – a cold disc from 15 to 70 au with a surface density profile Σ∝r^−1.5 and an unresolved hot component at 2 au. This is shown to the same colour scale and pixel size (1 arcsec at 100 μm and 2 arcsec at 160 μm) as the observations shown in Fig. 1. The 100-μm image is on the left, the 160-μm one is on the right. The residuals (observation – model) are below, with 100 μm on the left and 160 μm on the right, respectively.

Open in new tabDownload slide

There was a difference between the models preferred by the SED and by the images. The SED best fit requires a closer inner edge (r_in = 10 au) to produce enough flux in the 10–70 μm region whereas the images prefer a slightly larger inner radius (r_in = 25 au) with the same size distribution and grain composition. The best-fitting value of r_in is therefore a compromise between these two values. To constrain the error on the inner radius, a Bayesian inference method is adopted (e.g. Lay, Carlstrom & Hills 1997; Akeson et al. 2002; Pinte et al. 2008), in which each model is assigned a probability that the data are drawn from the model parameters. In cases where the Bayesian prior has a uniform probability distribution, as is the case here, this probability is ⁠, where χ² is the unreduced chi-squared. The normalization constant, P₀, is chosen so that the sum of the probabilities over all the models in the grid is unity. Once this is done for all the models in the grid, the probability distribution for a given parameter can be derived by marginalizing the eight-dimensional probability hypercube against the other seven dimensions. An example of this process for the disc inner radius is shown in Fig. 8, which shows that the distribution for possible locations of the inner edge peaks at 15 au with a range of 10–40 au.

Figure 8

The top figure shows Bayesian marginalization of the best-fitting model grid showing the probability distribution for the inner radius fit. The bottom figure shows Bayesian marginalization of the best-fitting model grid showing the probability distribution between the outer radius and the power-law index of the surface density distribution, γ. The highest probability is represented by the lightest colour, with a maximum probability of 0.53.

Open in new tabDownload slide

It is also possible to marginalize two parameters against the other six, creating a 3D plot with the probability distribution which allows us to examine the dependence of parameters upon each other. Fig. 8 also shows such a plot of the marginalization for the outer radius and γ, the exponent of the surface density power law. The outer radius and surface distribution are not well constrained by fitting the SED alone, partly due to the lack of constraints at wavelengths longer than 160 μm. The images suggest that the outer radius required is >60 au, but these parameters are degenerate as the steeper the surface density profile, the larger the outer radius needed to provide a good fit. The best-fitting choice of r_out = 70 au and surface density index γ=−1.5 represent the best compromise between these two factors. The grain size constraints on the composition and minimum grain size come solely from the SED.

Compared to the two-component disc model presented in Stock et al. (2010), which has a ring of 1 μm diameter carbonaceous grains in a ring at 2–3 au and a disc from 5 to 55 au with a flat surface density profile (γ = 0) composed of silicate grains with radii from 5 to 1000 μm. The Stock et al. (2010) model is driven by the resolved structure by BLINC and the excess level of 250 ± 50 mJy. Constraining from the excess range of 5.8–10.5 μm region (IRAC 5.8 and 8 μm, BLINC 10.5 μm), a temperature of 600 K is imposed for the inner hot component, resulting in an excess peaking at 5 μm. In comparison, the hot component in our model has a lower temperature due to a different size distribution and grain composition which means that it matches the Stock et al. (2010) model well at >5 μm but was not fitted in the IRAC bands due to differences in the photospheric model. There is a difference of 5 per cent in the photospheric fit used in the IRAC bands, but the photospheric subtraction at 10 and 24 μm produces values that are within the errors of those used by Stock et al. (2010), despite the difference in the stellar fit. When constraining the cold and warm components, we primarily consider the excess as >5 μm and refer to the Stock et al. (2010) model for a more thorough treatment of the hot excess.

The total mass in the collisional cascade, with the assumed size distribution of n(D) ∝D^−3.5, scales with ⁠. Thus scaling to the maximum grain size of 1 cm gives a dust mass of 2.1 × 10⁻⁴ M_⊕, with 3.2 × 10⁻⁸ M_⊕ of this in the hot component and the rest in the cold component. Assuming the collisional cascade extends up to bodies of D_max = 1 km, then the total mass would be 1.2 M_⊕. PSF variation can have an important effect on the observed residuals, so we repeated the modelling process using the Vesta PSFs and found no significant differences in the best-fitting model within the quoted errors, indicating that PSF variation has a negligible impact on the results.

3.1.1 Comparison with mcfost results

We also compared the results of our model to the best fit found by the mcfost code (Pinte et al. 2006), a Monte Carlo radiation transfer code in which a star radiates isotropically in space and illuminates an azimuthally symmetric parametrized disc. This was used to fit both Herschel 100- and 160-μm images and the entire SED of β Leo. A grid of approximately a million models was run with the range of parameters described above, and the joint χ² of the images and SED was minimized. The best-fitting model derived from fitting the SED and PACS images with mcfost gives qualitatively similar best-fitting parameters as those derived from our IDL modelling suite (see Section 3.2). As the two approaches give similar results, this validates the models within the constraints of the assumptions made about spatial and size distributions and compositions. The SED was weighted more heavily in calculating the χ² in the mcfost models, and the tension between the best fit for the SED and images results in a compromise with the model preferring either most of the mass in a relatively narrow ring around 30–40 au (as indicated by the SED), or a more extended disc starting further in (with smaller grains) and extending further out with a shallower surface density profile, which provides a better fit to the PACS images. As our model grid weighted the goodness of fit to the image surface brightness profiles more heavily than the SED fit, we consider the more extended disc as our best-fitting model. The 24-μm image was not included in the mcfost grid, and the extension seen in this image provides evidence against the narrow ring interpretation.

3.2 Three-component model

In Section 3.1 we describe a two-component model of the β Leo debris disc system. However, there is some compromise between the best-fitting models indicated by the SED and PACS images which leads to an uncertainty in the location of the inner edge of the disc between 15 and 30 au, as this region is within the PACS beam size. This ambiguity is due to limitations of the initial assumption of a disc with an inner radius, an outer radius and a continuous surface density profile. By considering the possibility of a more complex three-component model, we can place better constraints on the inner region of the disc. By combining the SED and surface brightness limits from the mid-IR and scattered light images, we can place limits on possible three-component models for the system and constrain any warm dust population between the hot and cold components. This three-component model has seven parameters: the radius and cross-sectional area of hot dust (r_hot, ⁠), radius and cross-sectional area of warm dust (r_warm, ⁠), radius and cross-sectional area of cold dust (r_cold, ⁠) and the composition of the warm component. We fix the composition of the hot and cold components to that inferred in Section 3.2, since we know that this results in an appropriate range of dust temperatures at a given distance to fit the PACS images. The parameters of hot dust are set by the limits from the 10-μm interferometry – here we restrict the parameter space to r_hot of 2, 5 or 8 au. The parameters of the cold dust are set by the PACS imaging as the best-fitting cold disc described in Section 3.2.

The main contribution from the cold component to the emission at 10–70 μm comes from near the inner edge of the distribution. We tried inner radii of the cold component of 15, 20, 25 and 30 au, all of which are consistent with the PACS imaging. The warm component can then be constrained through SED fitting (see Table 2), primarily considering the fit to the IRS and MIPS photometry, and then checked for consistency with the surface brightness limits from the 11.3 and 18.1 μm Gemini imaging (Section 2.3) and with the 24-μm MIPS image (Section 2.5). The best-fitting SED for the three-component model is shown in Fig. 9.

Figure 9

The SED fit for the three-component model of the β Leo disc. The SED shows the individual contributions (dotted lines) from a hot component at 2 au, a warm component at 9 au and a cold component at 30–70 au. The stellar photosphere is shown as a straight line and the cumulative disc distribution is a solid black line.

Open in new tabDownload slide

The constraints on the warm component can be examined by considering the percentage of the total 18-μm flux that is produced by each component in the disc. If we assume that the hot component is at 2 au, then for this to fit the BLINC observations at 10 μm, the warm component needs to produce 19 per cent of the total 18-μm flux, (48 mJy). The cold component with an inner edge at 15 au produces 81 per cent of the 18-μm flux, which means that there is no need for a warm component. This is shown in the first panel of Fig. 13. If the inner edge is moved to 20 au then only 70 per cent of the 18-μm emission is produced by the cold component, meaning that warm emission would have to account for 10 per cent (25 mJy) of the 18-μm emission. If the inner edge is at 30 au then the cold disc only accounts for 51 per cent (129 mJy) of the total 18-μm emission, leaving 29 per cent (74 mJy) that must come from the warm component.

Figure 13

An illustration of the possible structures of the β Leo disc based on the three different possible models of the disc. The top panel shows the two-component model, which has belt at 2 au, a gap from 2 to 15 au and an extended disc from 15 to 70 au. The middle panel shows the three-component model which has belts at 2, 9 and 30–70 au. The bottom panel shows the eccentric ring model as described in the text.

Open in new tabDownload slide

If we repeat this analysis with a hot component at 5 au (which produces 29 per cent of the 18-μm emission), we find that the inner edge of the cold component cannot be at 15 au because this produces too much of 18-μm emission, so the inner edge must be >20 au. For an inner edge at 25 au, we need a flux of 23 mJy at 18 μm (9 per cent of the total 18-μm flux) from the warm component. Similarly for an inner edge at 30 au we need a flux of 51 mJy (20 per cent of the total 18-μm emission). These possible configurations are shown in Fig. 13.

Although the above arguments imply that there are many possible hot dust radii and cold component inner edges, there are further constraints on the spectrum and many require unphysical assumptions about grain properties. There are two best fits to the three-component model. The first has dust at 2, 9 au (temperature of 160 K) and 30–70 au (γ=−2.0 for the cold component), whilst the second has dust at 5, 12 and 30–70 au (γ=−2.0). The same composition is used as that of the best-fitting two-component model (see Section 3.1). The interferometric constraints give a slight preference to the first scenario with the dust at 9 au. This gives us a = 4.45 with = 1.11, = 1.31 and = 2.03, showing that the three-component model is a better fit than the two-component model to the observations. As the three-component model has more free parameters than the two-component model (19 and 16 free parameters, respectively), it is expected to have a lower χ² than the two-component model, but as the reduced χ² takes into account the number of free parameters, ν, they can be directly compared. The significance of this decrease in χ² is discussed in Section 3.4. The SED with this fit is shown in Fig. 11.

Figure 11

Surface brightness profiles for the β Leo disc at 12 μm (top left), 18 μm (top centre), 24 μm (top right), 100 μm (bottom left) and 160 μm (bottom right) showing the fit of all the three models to the observations. These are line cuts taken parallel to the major axis of the disc (125°), summed over a width of 4 pixels. The observed profile is shown with a solid line, the model profiles with coloured dashed lines (red for the two-component model, blue for the three-component model and green for the eccentric ring model). The bottom panel of each plot shows the residuals, i.e. the observed line cut minus the model line cut to show the goodness of fit of the model.

Open in new tabDownload slide

3.3 Eccentric ring model

Two other stars that harbour hot dust at ∼1 au resolved interferometrically are HD 69830 (Beichman et al. 2006; Smith, Wyatt & Haniff 2009b) and η Corvi (Smith et al. 2009b). Both these systems can be fitted with a single continuous planetesimal population, in the form of a very eccentric (e > 0.9) ring (Wyatt et al. 2010) with a pericentre at the location of the hot dust and an apocentre corresponding to the cold belt for η Corvi, and an as yet unseen cold population for HD 69830. We used this model (Wyatt et al. 2010) as a third option to try to explain the emission seen around β Leo. This model assumes that the star formed with an eccentric planetesimal population and that this population has been evolving due to steady-state collisional erosion for the 45 Myr age of the system. The model has two dust components: cold dust coincident with the planetesimals and hot dust created in collisions at pericentre that is being removed by radiation pressure. The model parameters are the pericentre, apocentre, maximum planetesimal diameter and the size of the grains being removed by radiation pressure. The best-fitting model was found by minimizing the reduced χ² fit to the SED and the 12, 18, 100 and 160 μm images across a grid of model parameters. The range of parameters tested to find the best fit were pericentre from 0.5 to 5 au (interval 0.5 au), apocentre from 40 to 140 au (interval 10 au), and maximum planetesimal diameter between 500 and 3000 km (with an interval of 500 km). The population of grains being removed by radiation pressure were assumed to be a single size, and values from 20 to 100 per cent of the blow-out size were tested.

The SED for the best-fitting model is shown in Fig. 10. This model has a single planetesimal population with pericentres at 2 au and apocentres at 65 au. The eccentricity was therefore 0.96 and the maximum planetesimal diameter was 2000 km. The grains blown out from pericentre have a size of 2 μm. The resulting mass-loss rate given by this model is 0.005 M_⊕ Myr⁻¹, giving a mass of 0.45 M_⊕ for the parent planetesimal population. It is therefore possible to explain the β Leo disc using one continuous planetesimal population, but the possible origins of such an eccentric planetesimal population are an issue (Wyatt et al. 2010). The combined reduced χ² for this model is = 4.47, showing that the eccentric ring model has a slightly lower than the two-component model, but is a slightly worse fit than the three-component model. The significance of this change in reduced χ² is assessed in Section 3.4. The fit to the surface brightness profiles of the two-component, three-component and eccentric ring models at 12, 18, 100 and 160 μm are shown in Fig. 11. The 24-μm line cut shows that the eccentric ring model appears marginally more extended when compared with the observations. The FWHM of the eccentric ring model at this wavelength is 7.0 × 6.7 arcsec² compared to an FWHM of 6.88 × 6.61 arcsec² for the β Leo observations. This means that the eccentric ring model is slightly less favoured, but cannot be ruled out.

Figure 10

The SED of β Leo modelled using a single planetesimal population with pericentre at 2 au and apocentre at 65 au. The dashed, dotted and solid lines correspond to the contribution from the collisional cascade, blow-out grains and total emission spectrum, respectively. The diagonal solid line is the stellar spectrum, the Spitzer IRS spectrum from 10 to 40 μm is plotted with a solid line, the asterisks are the MIPS points, the triangles photometry is from Gemini Michelle and the squares PACS photometry. These are all shown with photospheric subtraction and have been colour corrected where necessary.

Open in new tabDownload slide

3.4 Ambiguity of the models

This addition of an extra parameter into each of the three-component and eccentric-ring models improves the from 4.50 to 4.45 and 4.47, respectively, where a perfect fit to the resolved images at 100 and 160 μm and the SED would have = 3.0. However, in general adding an extra parameter will improve the fit of any model, so is adding this extra parameter justified? Although the reduced χ² compensates somewhat for this as it is divided by the number of free parameters ν, when considering models with large numbers of free parameters (19 in the case of the three-component case) we need to avoid the problem of ‘overfitting’ or choosing a more complex model than is warranted by the data. We assess this using the Bayesian Information Criterion (BIC):

1

where N is the number of data points, k is the number of free parameters and is the minimum combined reduced χ² for the model. The BIC considers the fit of the model but penalizes the model for extra parameters (see Wahhaj et al. 2005; Liddle 2007). It penalizes the model more strongly for extra parameters than the reduced χ², and so it provides a more stringent test of the validity of applying more complex models. The BIC value for the two-component model is 51.1, for the three-component model it is 49.7, and for the eccentric ring model the BIC is 51.9. A lower value of the BIC is preferred; a difference of 2 between BICs of different models indicates positive evidence against the higher BIC value, and a difference of 6 indicates strong evidence against the higher BIC value. Therefore, our results indicate that there is no real preference between the three models since they all reproduce the data well. Furthermore, it is not unreasonable to assume that a realistic disc system may well be more complicated than the simple parametrizations used in these models.

These three models all provide a very good fit to the observations, but further data could allow us to resolve this ambiguity. The best way to distinguish between the possible three-component models and the two-component model would be to resolve the inner edge or the warm component. This may be possible with very deep mid-IR imaging in the N band. Resolving the disc at this wavelength would give different surface brightness profiles for the different models as shown in Fig. 12. The Spitzer image at 24 μm already suggests that the extent of the eccentric ring model may be too large. A quadratic subtraction of the FWHM of the beam from the image FWHM diameter (after PSF subtraction) indicates a diameter of ∼4-arcsec FWHM for the disc at this wavelength, significantly less than the value indicated in Fig. 13.

Figure 12

Predicted 11.3-μm surface brightness profiles for the three different models of the β Leo disc: the two-component model with a hot disc at 2–3 au and a cold disc from 15 to 70 au; a three-component model with a hot disc at 2–3 au; a warm ring at 9 au and a cold disc from 30 to 70 au and the eccentric ring model.

Open in new tabDownload slide

It may also be possible to resolve the inner edge of the cold disc in scattered light. We created simulated ACS images for all the three models to compare to Section 2.6. Mie theory was used to calculate the albedo of grains in the model and a bandpass of 0.45–0.72 μm was assumed for the F606w filter. A Henyey–Greenstein phase function (Henyey & Greenstein 1941) was used to approximate the asymmetric scattering by small particles with an asymmetry parameter (g) of 0.3 corresponding to slightly forward scattering. Due to uncertainties in scattering properties we might expect our predicted surface brightnesses to have significant uncertainties. However, the morphology for the emission (and relative level of surface brightness for different models) should be correct.

The ACS coronograph has an inner working angle of 1.8 arcsec, so in all the models the hot component is hidden behind the occulting spot. For the two-component model the peak surface brightness in scattered light would be a radius of 1.35 arcsec (i.e. behind the occulting spot) with a surface brightness of 1.7 × 10⁻³ mJy pixel⁻¹, assuming a pixel size of 25 mas. Just outside the inner working angle of the coronograph (1.9-arcsec radius) the disc would have a surface brightness of 1.3 × 10⁻³ mJy pixel⁻¹ and fall off ∝r^−3.5. The three-component model has a peak surface brightness of 2.1 × 10⁻⁴ mJy pixel⁻¹ at a radius of 2.7 arcsec and the eccentric ring model has a peak surface brightness of 4.6 × 10⁻⁴ at a radius of 3.6 arcsec. In an annulus of 2.8–3.0 arcsec the residual halo in the HST ACS observations shown in Section 2.6 has a median brightness of 5.3 × 10⁻⁴ mJy pixel⁻¹, which is consistent with the three-component model. The halo brightness decreases steeply to the sensitivity-limited sky value at 4.0 arcsec radius. If the halo were due to grain scattering, we would expect a radial brightness profile proportional to r⁻² if the star was embedded in a uniform density sheet of dust (γ = 0). The measured radial brightness profile is proportional to r⁻⁶, implying that γ∼−3. This is comparable to the dust radial profile for the south-west mid-plane of the β Pictoris dust disc (Kalas & Jewitt 1995; Golimowski et al. 2006). However, this is significantly steeper than γ=−2.0 for the three-component model.

The nature of this residual emission is unclear – it is difficult to find a disc model that can both produce sufficient scattering in this region without producing a very bad fit to the SED. Models with an outer edge at 44 au, as seen in the ACS observations or with γ < −2.0 to mimic the steep drop-off seen, are also a very bad fits to the extended emission seen with PACS at 100 and 160 μm. However, it could be reconciled with the two-component model if we are seeing a separate small grain population in the optical observations that has a steeper surface density profile. If this emission was real, then it could be used to place limits on the albedo of the grains in the eccentric ring model, or possibly rule out this scenario. Currently, the observations favour the two-component model as no peak is seen at 2.7 or 3.6 arcsec. Further optical imaging could confirm the nature and location of this emission and will possibly rule out the eccentric ring model, which has the largest surface brightness in this region.

4 DISCUSSION

The β Leo debris disc system consists of multiple dust populations. We can find a best-fitting solution with one, two or three components. These are shown in Fig. 13. The two-component model is the simplest and has dust at 2–3 and 15–70 au. There are two best fits to the three-component model. The first has dust at 2, 9 and 30–70 au, whilst the second has dust at 5, 12 and 30–70 au. The interferometric constraints give a slight preference to the first scenario. The key constraint on these models is that the warm (middle) component cannot be located inside an 8-au radius without violating the 10-μm interferometry constraints. A third possible scenario for this system is that the emission results from a single very eccentric (0.96) planetesimal population. This model consists of a steady-state collisional cascade with an apocentre coinciding with the location of the cold emission seen in the PACS images, with an additional population of hot, small grains created at collisions near pericentre which are subsequently removed from the system due to radiation pressure.

These inferred structures give us clues to the underlying dynamics of the system, as the belts are confined and separated by gaps suggesting dynamical interactions. The two-component model has a gap between 2 and 15 au which could contain planets, as could the gaps between 2 and 9 au and between 10 and 30 au in the three-component model.

4.1 Origin of the hot dust

The origin of the hot, compact dust populations is still unclear. There are five possible scenarios: (i) the small dust grains produced by collisions in the cold belt could drift towards the inner region due to PR drag; (ii) the hot dust could be evidence of a planetesimal population in a steady-state collisional cascade at a few au from the star; (iii) the dust could be transient, produced in a massive collision between planetesimals; (iv) the system could be undergoing a dynamical instability resulting in the planetesimals being thrown in from the outer belt and producing hot dust as in the Late Heavy Bombardment (Gomes et al. 2005); or (v) the grains may be produced through evaporation of comets originating from the cold planetesimal belt imaged in the sub-mm, like the Zodiacal Cloud (Nesvorný et al. 2010).

The first scenario concerning PR drag as an origin for the hot dust is improbable due to the long time-scale for PR drag compared to the shorter collisional time-scale for the outer disc (Wyatt et al. 2007). This mechanism is unlikely to produce the amount of dust observed in β Leo’s inner regions. Around late-type stars radiation pressure is inefficient at removing small grains, so these can be transported into the inner regions of the disc via Poynting–Robertson drag and stellar wind drag, even in discs with optical depths considered too high for discs to be transport dominated (τ∼ 10⁻⁴, Reidemeister et al. 2011). This is because close to the star larger grains are preferentially collisionally eliminated, causing a break in the size distribution at a critical size, below which grains are transport dominated. Since β Leo is an A-type star, however, the larger radiation pressure blow-out size will mean that all the grains remaining in the disc are larger than the critical size and hence are collisionally dominated.

The feasibility of a steady-state belt as an explanation for the hot dust can be examined by considering the maximum steady-state fractional luminosity of the belt. The fractional luminosity of the hot component from our model is 8.1 × 10⁻⁵. Using equation (18) from Wyatt (2008) for the maximum fractional luminosity of a belt at a given age,

and assuming a radius (r) for the hot component of 2.5 au, a width (dr) of 1 au, a maximum planetesimal diameter (D_c) of 60 km, planetesimal strength (⁠⁠) of 150 J kg⁻¹ (the fiducial value from Wyatt 2008), an eccentricity (e) of 0.05, stellar mass (M_*) of 2.1 M_⊙, a luminosity of 14.0 L_⊙ and an age of 45 Myr, gives a maximum fractional luminosity of 9.9 × 10⁻⁷. For the hot dust population to be considered transient, Wyatt et al. (2007) concluded that the fractional luminosity must be 1000F_max, so β Leo’s hot dust is consistent with a steady-state model, but only marginally so, and to be a steady-state phenomenon it would have to have properties that are significantly different from those for other A stars (e.g. unusually strong or large planetesimals). This implies that the third scenario is unlikely as again, a very large belt mass would be required due to the low probability of a massive collision having occurred and the short lifetime of the dust produced.

The fourth scenario also interprets the hot dust as a transient phenomenon by suggesting that the system is undergoing major dynamical perturbations in a scenario analogous to the late heavy bombardment (LHB) in the Solar system. The LHB occurred 700 Myr after the planets formed and led to a spike in the production of dust as comets were thrown in towards the star by a dynamical instability caused by the migration of the giant planets (Gomes et al. 2005). It is not known how common such events are, although Booth et al. (2009) estimate that they occur around <12 per cent of sun-like stars. Since an LHB-like event could occur at any point in a disc’s lifetime and as only a small percentage of stars undergo such an event, it is improbable that we are observing such a disturbance in the β Leo disc.

This leaves the final scenario, that the grains may be produced through disintegration of comets originating from the cold planetesimal belt imaged in the sub-mm. Since this requires neither a very high initial mass for the disc nor a low-probability event such as a collision or an LHB-like instability, it seems more likely although the dynamics of scattering and physics of dust production remain poorly understood. If the hot dust does originate from a population of scattered planetesimals then the amount of hot dust could be indicative of the configuration of the planetary system responsible for the scattering, e.g. the effect of Jupiter on scattering of planetesimals in the Solar system (Horner & Jones 2008, 2009; Horner, Jones & Chambers 2010). This implies that there could be multiple planets between the dust belts responsible for scattering planetesimals into the inner regions of the system, and so replenishing the hot dust population.

4.2 Placing β Leo into context

There are currently 16 debris discs around A-type stars that have at least one component resolved in thermal emission. These are listed in Table 3. Of these 16 discs 11 show evidence of multiple dust components, where at least one component has been imaged. The last five objects in Table 3 are A stars with evidence of multiple dust components including hot dust, but whose discs have not been resolved. Of these A stars, most are thought to be younger than β Leo, with ages in the range of tens of Myr.

The discs summarized in Table 3 show a startling diversity, with radii from a few au, i.e. ζ Lep (3–8 au, Moerchen et al. 2007a) to hundreds of au, i.e. Vega (Su et al. 2005). Some systems are confined rings with clear inner and outer edges, i.e. HR 4796 (Schneider et al. 1999) and some have very extended discs like β Pic (Kalas et al. 2000).

Fig. 14 shows a plot of the fractional luminosity versus radius of both the hot and cold disc components for all the A stars listed in Table 3 as having two disc components. Vega, Fomalhaut and β Leo are the only A stars with a hot disc component that has been resolved via interferometry (Absil et al. 2006, 2009; Akeson et al. 2009). The hot dust in other discs is either inferred from SED fitting (as for HR 4796 and HR 8799) or resolved in the mid-IR (ζ Lep, β UMa, HD 139006). This leads to potential systematic differences when comparing radii and fractional luminosities derived through different methods. These warm components could also be related to planetesimals near the ice-line (Morales et al. 2011). Although there are no clear trends in the small number of objects with known two-component discs, there are some interesting points to note.

Figure 14

Plot showing the radius and fractional luminosity of the hot and cold components for the A star debris discs listed in Table 3 as having two components.

Open in new tabDownload slide

We see two main ‘types’ of distributions between the hot and cold populations. First, discs such as Vega, Fomalhaut, β Leo and η Corvi, which have a large (>20 au) separation between the location of the hot and cold dust, and whose hot dust has a high fractional luminosity. This leads to a steep line in Fig. 14. The other type of disc, typified by Eta Tel, HD 3003 and HR 8799, show an almost flat or, for HR 8799, declining line in Fig. 14, indicating that the hot dust has a similar or lower fractional luminosity than the cold component. In the case of HR 8799 the hot and cold components are separated by a known planetary system (Marois et al. 2008) which may inhibit the transfer of material from the outer to the inner regions of the system.

Fomalhaut, Vega and β Pic show the steepest increase in fractional luminosity between cold and hot components. In all three of these cases the hot component is derived from modelling, fitting the interferometric visibility deficit seen around these objects (Absil et al. 2006, 2009, 2010). β Leo shows a gradient between the hot and cold components, similar to that of η Corvi, an F star with multiple dust populations including a cold dust disc resolved with PACS and SCUBA at 150 au (Matthews et al. 2010) and ζ Lep, which has a much more compact configuration, with its ‘cold’ disc at only 10 au.

In the Solar system the asteroid belt (radius ∼2 au, f = 0.8 × 10⁻⁷, Backman & Paresce 1993) and the Kuiper Belt (radius ∼ 30 au, f∼ 10⁻⁷, Gomes et al. 2005) would also give a flat line, similar to that of HR 8799. The radial distribution of dust in these systems could be indicative of the underlying planetary system. HR 8799 and the Solar system both have multiple planet systems between their dust belts which could influence the levels of planetesimal scattering. However, Fomalhaut and β Pic have at least one planet (Kalas et al. 2008; Lagrange et al. 2010) but still have high hot dust levels and a steep line on Fig. 14.

4.3 Implications for possible planetary systems

This diversity in disc structure and brightness seen in these discs with hot and cold components could be due to the underlying dynamics of the system, which can be affected by the presence and distribution of planets in the system. For example, in the Solar system the number of comets thrown into the inner system from the Kuiper Belt is dependent on planetary architecture (Horner et al. 2010). For the β Leo disc, the radial distributions of dust that we have derived have implications for the architecture of possible planetary systems. For the two-component model there is a gap in the disc between 2 and 15 au, in which there could be planets that are truncating the cold disc and confining the hot inner disc. For the three-component model the gaps from 3 to 9 au and from 9 to 30 au could occur again due to planets truncating the discs. The size of these gaps could be linked to clearing by a planet due to the scattering of planetesimals out from the planet’s chaotic zone. The area cleared increases for multiple planets and for planets on eccentric orbits (Bonsor, Mustill & Wyatt 2011). In the case of multiple planet systems these could influence the amount of material scattered into the inner Solar system and the possible stable locations for dust belts, so the morphology of the dust belts in the β Leo system could indicate multiple planets between the cold and hot dust belts.

5 CONCLUSIONS

We have presented detailed modelling of the debris disc around the 45-Myr-old A star β Leo. We considered multiwavelength data to construct a complete picture of this source. Resolved images taken at 100 and 160 μm using PACS on Herschel as a part of the DEBRIS survey were used to place observational constraints on the radial location of the cold dust in this system. Resolving the cold dust is key to breaking the degeneracies inherent in SED modelling. We also use detection limits from unresolved images at 12 and 18 μm from MICHELLE at Gemini, at 0.6 μm with ACS on HST, the unresolved Herschel SPIRE image at 250 μm and detailed SED modelling including all data from the literature to gain a complete picture of the disc. Modelling indicates that for a two-component model of the system consisting of a hot and cold disc, the cold disc imaged with Herschel PACS lies between 15 and 70 au, with a surface density profile Σ∝r^−1.5 at an inclination of 55° from edge on, with the hot dust at 2 au. SED fitting to observations from 5 μm to 1 mm, including the IRS spectrum (Chen et al. 2006), suggest that the grain size distribution is consistent with that for a theoretical collisional cascade with a fixed maximum grain size of 1 cm, a minimum grain size of 3 μm (⁠⁠) for a composition with a silicate fraction of 20 per cent, a porosity of 20 per cent with the rest of the grain composed of organic refractories. For the hot component the minimum grain size is 0.6 μm (⁠⁠) and the silicate fraction was 60 per cent. A three-component model indicates that another possible configuration consists of an inner edge to the cold disc at 30 au, a warm ring at 9 au and hot dust at 2 au. It is also possible to fit the observations of β Leo using a single, very eccentric (e = 0.92) planetesimal population, after Wyatt et al. (2010). This degeneracy implies that even a wealth of multiwavelength data including the resolved images may not be enough to uniquely constrain the location of the dust when there are multiple populations and the edges of the belts have not been resolved.

The age of β Leo is close to that of Fomalhaut and Vega, which also have hot dust within a few au of the star. However, its cold disc is smaller and closer to the hot emission than either of these systems, and there may be an intermediate warm component. The most analogous disc among those resolved around A stars is Eta Tel, a 12-Myr star with a resolved disc from 21 to 26 au and hot, unresolved dust at <4 au (Smith et al. 2009a), but in terms of the separation and ratio of the hot and cold components as shown in Fig. 14, β Leo most resembles HD 71155.

We have also examined the population of A stars known to have both hot and cold dust discs, and we see two main ‘types’ of distributions between the hot and cold populations. First, discs such as Vega and β Leo which have a large (>20 au) separation between the location of the hot and cold dust, and whose hot dust has a high fractional luminosity. The other type of disc, typified by Eta Tel and HR 8799, shows an almost flat or, for HR 8799, declining line in Fig. 14, indicating that the hot dust has a similar or lower fractional luminosity than the cold component. In the case of HR 8799 the hot and cold components are separated by a known planetary system (Marois et al. 2008) which may inhibit the transfer of material from the outer to inner regions of the system. However, the small number of stars in this sample and the possibility of unknown systematic issues from comparing dust detected through different methods limit the conclusions that can currently be made.

LJC would like to acknowledge the support of an STFC studentship. The authors would like to thank Andy Gibb for his help with the SCUBA-2 data reduction. The authors would like to thank James DiFrancesco for his close reading of this paper as an internal referee and his very helpful comments that much improved the paper. Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA.

REFERENCES