© ESO, 2014

1. Introduction

The debris disc around β Pictoris is one of the most studied planetary systems. Smith & Terrile (1984) discovered the debris disc when the Infrared Astronomical Satellite (IRAS) detected an infrared excess associated with the star (Aumann 1985), and since then it has been studied from the visible and ultraviolet to the millimeter wavelengths. In optical/near-infrared scattered light, the debris disc appears as a bright dust disc seen almost edge-on (e.g. Kalas & Jewitt 1995). It can be decomposed into two components (Augereau et al. 2001): an inner warped disc detected up to around 150 AU and an outer main disc extending up to 1000 AU (Smith & Terrile 1987). Currie et al. (2013) reported the detection of the disc with the NAOS-CONICA (NaCo) instrument at the Very Large Telescope (VLT) in the L′-band, but did not investigate and discuss its properties. Observations in this wavelength range still dominated by scattered light from the star, but close to the transition with thermal emission, could bring precious additional constraints on the dust properties. In addition, adaptive optics (AO) corrections from ground-based telescopes perform better at longer wavelengths, with Strehl ratios typically between 60% and 80%, while the angular resolution (~0.1″) is still sufficient to probe the disc inner regions, typically within a few tens of astronomical units. These regions are critical because the disc interacts with potential planets therein. For instance, Mouillet et al. (1997) attributed the observed inner warp of the disc to the gravitational perturbation of a massive body located on an inclined orbit. This perturber, discovered later by Lagrange et al. (2009), was indeed a giant planet orbiting within 10 AU from the star (Chauvin et al. 2012), and the inclination of its orbit with respect to the main outer disc was confirmed (Lagrange et al. 2012a). The position angle (PA) of the main disc was found to be 29.0° ± 0.2°, whereas the inner warp is offset by .

Imaging debris discs in the L′-band centred around 3.8μm is a difficult task. The thermal emission from the telescope/instrument and the sky emission, called hereafter the background emission, is much larger than at shorter wavelengths, while the star to disc contrast remains the same¹. Quasi-static speckles from imperfect optics prevail inside 1′′ whereas background emission is the dominant noise source beyond that separation. Therefore, high contrast and high angular resolution techniques are required, as well as specific data reduction tools. Innovative methods such as angular differential imaging (ADI, Marois et al. 2006), locally-optimized combination of images (LOCI, Lafrenière et al. 2007), and principal components analysis (PCA, Soummer et al. 2012; Amara & Quanz 2012) were developed to detect and characterize faint companions. They are well suited to image point sources. For discs, they are also used to reveal high spatial frequency features (Buenzli et al. 2010; Thalmann et al. 2011; Boccaletti et al. 2012; Lagrange et al. 2012a or Lagrange et al. 2012b), but any extended structure is strongly affected by these image processing techniques (Milli et al. 2012; Esposito et al. 2014). In this paper, we present new images of the disc obtained with VLT/NaCo in the L′-band (Sect. 2). We detail the morphology from 3.8″ down to 0.4′′ (Sect. 3), and perform forward modelling to correct for the biases induced by ADI (Sect. 4). Finally, we discuss our results in Sect. 5, before concluding in Sect. 6.

2. Observations and data reduction

2.1. Observations

The star was observed on 31 January 2013 during the science verification observing run of NaCo’s Annular Groove Phase Mask (AGPM) coronagraphic mode (Absil et al. 2013). The AGPM coronagraph operates in the L′-band at 3.8μm with an undersized Lyot stop blocking the light from the secondary mirror and telescope spiders (Mawet et al. 2013). Because the signal-to-noise ratio (S/N) was not high enough on this individual data set to retrieve the disc morphology with sufficient accuracy², we have used six additional data sets that had been recorded for orbit monitoring. These additional data were obtained in pupil-tracking mode, the star being saturated. A summary of the seven data sets with the average observing conditions is shown in Table 1.

2.2. Data reduction of individual data sets

For the AGPM observations, the sky subtraction was different from that described in Absil et al. (2013). The main source of noise that needs to be overcome to reveal faint extended emission beyond 1′′ comes from the high and variable background at this wavelength. The standard way to evaluate the background in high contrast coronagraphic imaging is to open the AO loop and perform an offset of the telescope to measure an nearby empty sky region. This technique does not provide a sufficient accuracy in our case because the telescope thermal emission significantly changed, depending on the AO loop status (open or closed). Therefore, the cosmetic treatment of the raw data only consisted of dark subtraction, flat-field correction, and bad-pixel interpolation. The background is removed in a later stage, during the ADI subtraction process, along with the stellar halo. The data were binned every 100 s, resulting in a data cube of 90 images, and the images were re-centred with the same procedure as described in Absil et al. (2013), leading to an uncertainty in the AGPM centring of 0.1 mas and an uncertainty in the star centring with respect to the AGPM centre of 8.5 mas.

In the case of the six non-coronagraphic data sets, the sky was evaluated using a different telescope offset position and subtracted from the raw images. The same cosmetic treatment was then applied. Because of the small telescope offsets, the total field of view is limited to a radius of 3.8″. The star centre position is obtained from a fit to the low-flux level wings of the saturated image and has an uncertainty of 7 mas, obtained with the same method described in the appendix of Lagrange et al. (2012a).

Then, the star subtraction algorithm consisted of a PCA performed in two to five concentric annuli between 0′′ and 3.8′′. This method uses a truncated basis of eigenvectors created by a Karhunen-Loève transform of the initial set of images, to perform the subtraction of the point-spread function (PSF). The number of principal components retained varied for each data set between 1 and 4. It was adjusted visually to maximize the disc SNR on each individual reduction. We did not do any frame selection to maximize the exposure time, thus the signal from the disc.

2.3. Combination of the data and subtraction of β Pictoris b

On the individual data sets, the residual thermal noise is non-Gaussian at L′ and shows large patterns of the size of several resolution elements. For instance, for resolution elements between 0.5′′ and 1′′, the Shapiro-Wilk test (Shapiro & Wilk 1965) yields a p-value below 2% for all individual data sets, except for the AGPM image (p-value of 90%), and this is highly dependent on the number of principal components retained for each data set. In order to increase the disc SNR, we combined several uncorrelated data taken at different epochs. This represents the first approach of that kind in high contrast imaging. Combining seven data sets reduces the standard deviation by a factor 2.7 for the resolution elements between 1′′ and 2′′, which is also what is expected from the central limit theorem assuming uncorrelated noise ().

For each data set, an AO-corrected unsaturated PSF recorded before and/or after the science observations was used as a photometric calibrator. The six final images were individually scaled to the total flux of the central star β Pictoris, cropped to the same field of view of radius 3.8′′, and averaged. The total flux of the star was calculated by performing aperture photometry on the AO-corrected unsaturated PSF. The flux increases with the aperture radius until convergence is reached for a radius of 1.35′′, therefore, we used this radius to make sure we encircled the total star flux and are not subject to the PSF variability due to the AO correction. Note that the disc is well below the noise level in the unsaturated PSF and does not bias the photometry of the star.

As visible in Fig. 1, the presence of β Pictoris b prevents from detecting the disc inside 0.6″ on the south-west (SW) extension. The planet orbital motion between 2009 and 2013 induces a bright dot smeared along 5.5 px, or 150 mas. To disentangle any disc structure or additional companions from the signature of β Pictoris b, we subtracted the planet signal from the seven data sets. We used the negative fake planet technique described in Lagrange et al. (2010) and Bonnefoy et al. (2011) to retrieve the photometry and astrometry of the planet at each epoch. We used the same optimization area as described in Absil et al. (2013). The photometry and astrometry derived for the position of the planet in each individual data set are listed in Table 2. The astrometry error budget includes the statistical error in the companion position calculated with negative fake planets as in Absil et al. (2013), the calibration errors from the platescale, true north and rotator offset, and the error in the star centre described in Sect. 2.2. The photometry error budget includes the statistical error plus a contribution accounting for photometric variations of 0.05 mag for a photometric night (ESO data quality control standards).

Fig. 1 PCA-reduced image obtained after combination of the seven individual images, colour scale-adjusted to show the smearing of the planet (green arrow), which prevents disc signal extraction at this location. The green circle has a radius of 0.4′′ and is centred on the star (green square). The black arrow has a position angle of 29° and a total length of 6″.

Similarly, disc non-axisymmetric features or additional point sources in orbital motion around the star would also smear an arc in the image. Assuming these features are in keplerian motion around the star, we explored the amplitude of this effect and show the results in Fig. 2 for four different semi-major axis. A feature orbiting at 30 AU from the star can smear an arc of 240 mas, and this value goes down to 170 mas for a 60 AU-orbit and 135 mas for a 100 AU-orbit. Therefore, we cannot exclude that our image contains disc features smeared over a few hundred mas that would otherwise have appeared more compact.

Fig. 2 Projected distance smeared by a potential disc feature or planet in keplerian motion around the star during the time span between the first and last observations (3.1 yrs).

3. Characterization of the disc

3.1. Correction for biases due to ADI

Fig. 3 Two different disc reductions: cADI-disc (top) and PCA-disc (bottom). The black arrow shows the apparent midplane, has a position angle of 29° (equal to that of the main disc), and a total length of 6′′. The colour scale is identical for both images but not linear (square root). The green square marks the position of the star and the green circle has a radius of 0.4″.

The signal of extended objects may be heavily influenced by ADI (Milli et al. 2012). Therefore, a careful calibration was performed. We used three approaches, including:

• 1.

  iterations as in Lagrangeet al. (2012a) in order to obtainimages that have the least biases from ADI;

• 2.

  forward modelling to constrain the dust density distribution as in Thalmann et al. (2014), presented in Sect. 4.2; and

• 3.

  injection of model discs in the data to calibrate the morphological parameters measured in the images, as in Boccaletti et al. (2012). This is presented in Sect. 5.

The iterative approach was already implemented in Lagrange et al. (2012a) and is called cADI-disc, where cADI refers to classical ADI and consists of a simple subtraction of the median of the cube of images to subtract the star. The resulting image is shown in Fig. 3 (top). Here, we also adapted this technique for the PCA, namely we ran a first reduction to get an image of the disc. After applying a mask to this image to set the regions unaffected by the disc to 0, we subtracted this image from the initial data cube. We therefore obtained a library that ideally does not contain any astrophysical signal and we applied a PCA on our initial data cube using this reference library. The reduced image is shown in Fig. 3 (bottom). The disc in the cADI-disc image is slightly thicker than in the PCA-disc image because disc self-subtraction is higher with the latter algorithm. We will therefore use the cADI-disc image to extract the surface brightness profiles of the disc. Speckle noise is weaker in the PCA-disc image, so we will use this image to extract the morphological parameters of the disc other than the vertical thickness. It should be noted that this iterative technique minimizes disc self-subtraction but it does not totally prevent it, in particular, in faint vertically extended parts of the disc, far from the midplane. The pixels of the image that contain the disc signal at all observed parallactic angles (because the disc is vertically too extended) still suffer from disc self-subtraction. This affects essentially the pixels inside 0.8″ and is very clearly visible on the PCA-reduced image of a disc model (Fig. 10, middle, described later in Sect. 4.3). Space-based observations that do not use field rotation for PSF subtraction indeed reveal a faint emission at all azimuths in the sky plane at a radius of 0.5″ (e.g. Fig. 2 of Golimowski et al. 2006).

3.2. Disc morphology

The disc is detected unambiguously from underlying speckle and thermal noise up to the maximum field of view available, i.e. 3.8″, and down to 0.4″. This is the closest separation ever reached on the disc. The SNR map of the disc is shown in Fig. 4. In the subsequent sections, we use the terminology spine and apparent midplane of the disc to refer to the curve joining the brightest pixels of the disc in each vertical profile, and to the line going through the star with a PA of 29°, corresponding to the best estimate of the PA of the main outer disc at Ks, as seen in Sect. 1.

Fig. 4 S/N map of the PCA-disc image. The colour scale is linear and spans from 0 to 10. The green circle indicates the start of the 5σ detection at a radius of 0.4″.

3.2.1. Offset and bow of the disc

The disc is clearly vertically offset from the star within 2″: the star lies below the spine. This is illustrated in Fig. 5, in which each vertical profile has been normalized to its peak value to highlight the position of the spine. At 1″, the spine is about 2 AU above the midplane of the disc. This is also clearly visible on the isophotes of the disc in Fig. 6.

Moreover, these two figures show that the disc is bowed towards the north-west (NW), i.e. above the apparent midplane, between −3″ and 3″. This bow of the spine in the innermost regions was already marginally detected by Boccaletti et al. (2009) in the H and Ks band, and was reported in the optical by Golimowski et al. (2006) who mentioned it for the SW extension. This bow can be reproduced by anisotropic scattering of the dust, combined with a slight inclination of the disc with respect to an edge-on view, and it will be discussed later. Two additional observations support this interpretation. First, the isophotes are more widely spaced above the apparent midplane (NW side) than below (SE side). This was previously observed by Golimowski et al. (2006). Second, the disc appears thicker above the spine than below it (see Figs. 3 and 6).

Fig. 5 PCA-disc image where each vertical profile has been normalized by the maximum spine brightness to enhance the vertical position of the spine. The two green circles have a radius of 0.4″ and 2″ respectively. The PA of the main disc (29°) as measured by Lagrange et al. (2012a) is indicated by a solid red line, whereas the best fit PA of 30.8° (NE) and 211.0° (SW) as measured from the L′ images are shown with dashed lines. The waviness of the spine is discussed in Sect. 5.2.

Fig. 6 Isophotes corresponding to the cADI-disc image of Fig. 3. The black arrow passes through the star (green square), has a position angle of 29°, and a total length of 6″.

3.2.2. Position angle of the disc

We extracted the PA of the disc along with the position and maximum brightness of the spine from our cADI-disc and PCA-disc images with the same technique as described in Lagrange et al. (2012a) for the Ks image. In the Ks band, Lagrange et al. (2012a) found a PA of 29.3° ± 0.2° for the main disc by fitting a Lorentzian profile beyond 4.8′′. To measure the PA of the warped inner component, Lagrange et al. (2012a) decomposed the vertical profiles into two Lorentzian profiles and found in that case a PA of 29.0° ± 0.2° and for the main and warped discs (see Table 2 in Lagrange et al. 2012a).

Unlike in the Ks image, the field of view available in our L′ image is too small to access a region dominated only by the main disc, which would allow us to determine the PA of this outer component and derotate the images. If we try to find the best PA of the disc by nulling the spine vertical position on average between 3″ and 3.8″, we get a PA value for the north-east (NE) and SW extensions of 30.8° ± 0.6° and 211.0° ± 0.7°. The same technique applied to our Ks images yields a similar PA within error bars (30.2° ± 0.5° and 210.5° ± 0.6°), confirming that the L′ image is compatible with the Ks image and confirming that our view of the disc is a superimposition of the two contributions of the outer main disc and the inner warped disc. Figure 5 shows the lines of PA 29°/ 209° and 30.8°/ 211.0° on top of the spine location. This shows that within 3.5″ the optimal PA of 30.8°/ 211.0° (red dashed line) is indeed the most compatible with the L′ data. At the largest separations beyond 3.5″, the disc tends to align again on the PA of the main disc (red solid line). We show in Sect. 5.4 that this is perfectly compatible with a two-component disc if the inner component is brighter than the outer component within 3.5″.

The dynamical range is not high enough, however, to allow us to fit a two-component Lorentzian profile and separate the contributions of the main disc from that of the warped disc. Therefore we decided to use the PA of 29° measured from the Ks image as an a priori to derotate our L′ image.

3.2.3. Position of the spine

The spine position as measured from the Lorentzian fit is displayed in Fig. 7. The bow of the disc is clearly confirmed by this plot and its extension can be determined: the spine is a concave function of the separation between −3.″ and 3.5″. On the NE side, the spine of the disc goes below 0 beyond 2.5″ and the slope changes sign beyond 3″. This fact is also well reproduced by the model of the disc presented in Sect. 5.4. The spine of the disc versus separation is not a smooth curve but present ripples. We show in Sect. 5.2 that these small-amplitude ripples are due to the background and speckle noise and a smooth curve is very likely, given the error bars, however, we note a larger ripple at 0.8″ on the NE extension, which is very clear in Fig. 5. This deviation is at 2σ above the noise level and its significance will be discussed in Sect. 5.2.

Fig. 7 Departure of the spine with respect to the apparent midplane (defined in Sect. 3.2 as the plane at 29° PA) as measured on the PCA-disc image. NE is to the left. The uncertainty, shown here in red, is the 3σ error in the centroid position of the Lorentzian profile fitted to each vertical profile, after setting a radially-dependent noise level to each pixel. The waviness of the curve is due to the background and remaining speckle noise, as shown in Sect. 5.2. We note, however, a larger ripple at 0.8″ on the NE side.

3.2.4. Brightness of the spine

The surface brightness of the spine as a function of the separation is plotted in Fig. 8. The noise is evaluated as the standard deviation measured azimuthally in the regions where the disc is not detected. The two extensions of the disc show no overall brightness asymmetry within error bars inside 1.5″. There is one localized asymmetry slightly inside 2″: the SW extension is locally 0.5 mag brighter than its NE counterpart. A localized brightness asymmetry was also detected in mid-infrared images of the disc by Telesco et al. (2005) at wavelengths beyond 8μm, but at a larger separation, namely 2.7″. In the optical, the SW extension is also slightly brighter than the NE extension between 50 and 100 AU or 2.5″ and 5″ (Golimowski et al. 2006).

The surface brightness profile appears smooth between 0.5″ and 3.8″, compatible with a single power-law dependance with the separation. The inflection seen in the optical at 2″ by Golimowski et al. (2006) is not detected in our data. A linear regression to the NE and SW extensions between 0.5″ and 3.7″ yields slopes of −2.77 ± 0.18 and −2.57 ± 0.16. The linear fit is shown as a dotted line in Fig. 8. The steeper slope of the NE extension explains why it appears slightly fainter beyond 1.5″. The error bar shown in Fig. 8 only includes the measurement error. Although the image was corrected for self-subtraction by applying the iteration technique described in Sect. 3.1, significant self-subtraction could still occur inside 0.8″, therefore, the brightness inside 0.8″ should be considered as a lower bound. To overcome this difficulty, we performed forward modelling with an innovative approach.

Fig. 8 Surface brightness of the spine, measured on the cADI-disc image from the fit of a Lorentzian profile. The red and purple shaded areas show the 3σ error in the amplitude of the Lorentzian function fitted to each vertical profile.

4. Forward modelling

4.1. Modelling philosophy

We used scattered light disc models to interpret the observed features and to disentangle ADI artifacts from real features. The disc models were generated with the GRaTeR code (Augereau et al. 1999). For each of the seven individual cubes of frames, the disc model is rotated to the appropriate parallactic angles of the initial frames and then subtracted. The resulting cubes are re-reduced using the same PCA algorithm as described previously. The seven reduced images are then combined to obtain one single disc-subtracted image. These steps are repeated iteratively by varying the free parameters of the disc model until a merit function is minimized. The minimization algorithm is a downhill simplex method or amoeba. For each minimization, three different sets of initial conditions were explored, all representing physically acceptable conditions to reduce the risk of finding local minima. We found that all sets agreed within less than 5%. The merit function is a reduced chi squared computed in the part of the image where the disc is detected³. Thalmann et al. (2014) used this forward modelling approach as well to study the disc of LkCa 15.

4.2. Assumptions

We neglected thermal emission over scattering by the dust particles. We discuss this assumption further in Sect. 5 where we show that, even if the thermal emission might not be negligible depending on the exact grain properties (especially their composition), the thermal emission cannot fully dominate the scattered light above 0.5″ and the main conclusion concerning the dust distribution will be unchanged. To limit the parameter space, we assumed an azimuthally symmetric distribution of optically thin dust with constant effective scattering cross-section. We used a Henyey-Greenstein phase function, whose single parameter g is not computed from theoretical dust optical properties but is adjusted to the data. We did not attempt to model the warp in a first approach but used a single population model similar to the single dust population presented in Ahmic et al. (2009).

We followed the work of Augereau et al. (2001), called A01 hereafter, who derived the surface density of the parent belt (PB) of planetesimals in the disc from dynamical modelling. Their model locates the PB between 50 and 120 AU and can explain both the surface brightness of the disc in scattered light and its infrared emission. It can also reproduce the recent submillimetric observations at 870 μm obtained with ALMA with a resolution of ~0.6″ (Dent et al. 2014) very well. We therefore assumed a dust distribution matching that of the PB within this range. Beyond a radius r_out = 115 AU, we extrapolated the dust surface density with power laws to account for the grains that are blown away by radiation pressure (but still bound to the star) and that populate the regions outside the outer edge of the parent belt. Pantin et al. (1997) argues that the dust density should decrease for distances larger than 117 AU following a power law with an index within the range [–1.5; –2.3] according to visible observations. This is in agreement with theoretical expectations (e.g. Thébault & Wu 2008). We therefore tested three different laws of indices p_out = −1.5, − 1.7 and −2.3. Little is known about the presence of dust within 50 AU. Resolved mid-infrared images by Pantin et al. (1997) and modelling of the spectral energy density by Li & Greenberg (1998) suggest that most of the dust lies beyond 30 AU but do not exclude an inner dust population around 10 AU. This is supported by the work of Okamoto et al. (2004) who revealed what may be several planetesimal belts of amorphous silicate grains. In this context, we consider three different scenarios for the dust distribution within 50 AU. In the first scenario, little dust lies inside 50 AU and the surface density follows that of the PB. In the second and third scenario, the dust surface density follows a power law of index p_in = 0.5 and 1.5 respectively.

All in all, we ended up with 3 × 3 = 9 models for the dust surface density profile. They are summarized in Fig. 9. We had to further assume a vertical distribution Z(r,z) for the dust to convert our surface density into volume density and compute the scattered light images. To do so, we followed Ahmic et al. (2009) who assumed that the dust vertical volumetric density follows an exponential profile with an exponent γ = 0.5, a scale height ξ₀ at the radius r_out, and a flaring exponent β = 1.5 so that .

Fig. 9 Description of the dust surface density used in the nine different models, based on the surface density derived by Augereau et al. (2001) for the parent bodies (PB, black line). Each of the nine models is made of three parts delimited by the two vertical dotted lines at 50 AU and 115 AU: an inner region (purple curves) following either a power law or the PB distribution, a common intermediate region following the PB distribution, and an outer region following a power law (green curves).

We observed one degeneracy in the remaining free parameters of the disc. The inclination of the disc i and the anisotropic scattering parameter g are degenerate. More inclined discs require a smaller g parameter to explain the offset and bow of the spine. For this reason, we considered explicitly distinct inclination values i of 89°, 88°, 87°, 86°, and 85°. We did not consider an edge-on disc (i = 90°) because the offset of the spine can only be reproduced assuming anisotropic scattering combined with a slight inclination of the disc. Likewise, we did not explore an inclination below i = 85° because it produces a disc significantly thicker than observed. The remaining free parameters for the model are therefore the PA, ξ₀, | g |, and a flux scaling factor.

4.3. Results

We show in Table 3 the best parameters corresponding to the most relevant cases in the (i,p_in,p_out) exploration space. They are relatively robust to the different scenarios investigated. The values of the best PA all lie between 30.6° and 30.9°, for instance, and for a given inclination, we observe little dispersion in the parameters | g | and ξ₀. Overall, the best χ² is obtained with an inclination i of 86°, and the three best χ² are obtained with an index p_in of 1.5, which is an indication that the dust surface density is not as steep as that of the PB within 50 AU. The best model is shown in Fig. 10 (top image), while the middle image shows the disc model after PCA-reduction without noise to illustrate the biases of ADI. We also show the residuals after subtraction of the model from the data (bottom image). Disc residuals are still present in this model-subtracted image inside 1″, mainly above the apparent midplane, and this will be discussed in the next section.

For all disc models with an inclination i = 86°, we observe that an inner dust distribution following that of the PB systematically yields the worst χ², and a value of p_in = 1.5 always obtains the best figures of merit. The different values of p_out yield very similar results. For a higher inclination of 88°, models are more anisotropic and with larger scale heights. This corresponds to the degeneracy between i and g mentioned earlier and to the fact that a vertically thicker disc is needed when the inclination is closer to edge-on. For such an inclination, we can reasonably rule out the less steep inner dust distribution p_in = 0.5 because it systematically leads to a higher χ².

Fig. 10 From top to bottom: initial disc model, disc after PCA reduction without noise, residuals remaining after the disc model was subtracted from the data cubes and the cubes were re-reduced. The colour scale is linear and two times larger for the initial model than for the other two images. The green circles have a radius of 0.4″ and the black arrows indicate the apparent midplane (PA = 29°) with a total length of 6″.

We conclude from this analysis that the PA of the disc as measured within the 3.8″ field of view lies between 30.9° and 30.6°. This is in agreement with the value of the PA measured from our single-Lorentzian fit: the PA of the disc as measured inside 3.8″ is greater than that measured at larger distances where the warp component is not detected.

The disc scatters light anisotropically with a value of the Henyey-Greenstein parameter g dependent on the inclination of the disc. For a disc inclined by 89°, a relatively strong anisotropic scattering with | g | = 0.6 is necessary to reproduce the bow of the disc. There is a moderate anisotropy (| g | = 0.4 to 0.5) for a disc inclined by 88°, and only a slight anisotropy (| g | = 0.3 to 0.4) for an inclination of 86°. If this anisotropy is because of the enhanced forward-scattering efficiency of the dust, as is often claimed (Golimowski et al. 2006), then the NW side is the closest to the Earth. These g values are compatible with those derived in the near-infrared for other debris discs, for instance | g | = 0.3 ± 0.03 for HD 181327 (Schneider et al. 2006), and Rodigas et al. (2012) could also explain the bow of HD 15115 detected at L′ with an anisotropic scattering parametrized by | g | = 0.5.

5. Discussion

5.1. Possible scenarios to explain the missing flux inside 1′′

To test further the agreement of our models with the data, we injected a fake disc (our best-fit model for i = 86°, p_in = 1.5 and p_out = −2.3) at a PA 90° away from the disc in the data cubes. The result is shown in Fig. 11. Visually, the model looks very similar to the observed disc beyond 1″ but appears too faint inside that region, as already seen in the model-subtracted image of Fig. 10. Quantitatively, the reduced fake disc is not detected inside 0.5″ and lacks about 0.3 mag/arcsec² at 0.8″.

Fig. 11 Image obtained after injection of the best disc model at 90° from the real disc and reducing the data again. The field of view is 3.8″ and the real disc extends from NE to SW.

We discuss three possibilities to explain the excess emission within 1″ compared to our best scattered light model:

• 1.

  There is more dust inside 20 AU than predictedby our model. This warm dust would therefore enhance thescattered light from the star and could also contribute to thermalemission at L′. Wahhaj et al. (2003) estimate from thermal imaging at 17.9μm that there is indeed an inner ring at 14 AU. Pantin et al. (1997) showed that mid-infrared observations are compatible with a dust depletion around 20 AU, but the surface density could increase again inside 10 AU (see the solid red curve in Fig. 14 described later in Sect. 5.3). Augereau et al. (2001) also proposed a population of hot grains within 20 AU to explain the 12 μm flux. Last, interferometric measurements showed that the star exhibits a near-infrared excess at H and Ks (Di Folco et al. 2004; Defrère et al. 2012) that cannot be entirely explained by scattered light from the edge-on disc within the interferometric field of view. All these observations are supported by the fact that none of our models can explain the strong brightness of the disc within 1″. To investigate this possibility, we used GRaTeR to compute the amount of scattering and thermal emission occurring along the line of sight for our best-fit dust distribution model, assuming two different types of grain populations: porous astronomical silicates (Draine 2003) or more realistic aggregates⁴ made of a silicate core coated by an organic refractory mantle, with water ice partly filling the holes created by the porosity, if the temperature is less than the water sublimation temperature (Augereau et al. 1999). In both cases, the grain size distribution follows the traditional a^-3.5 power-law decrease resulting from a collisional cascade with a minimum size a_min. The result is shown in Fig. 12. It confirms that porous silicates cannot explain the missing flux because scattering predominates by a factor over 100 beyond 0.2″. For the population of more complex aggregates, we cannot exclude a contribution of thermal emission dominant up to 0.5″ and still significant between 0.5″ and 1″, especially if the grain population is micronic or sub-micronic, which is expected from the population in the inner 20 AU (Augereau et al. 2001). At 0.8″, a ratio between thermal emission and scattered light of 0.3 is enough to explain the 0.3 mag missing flux, which is obtained for a grain population of minimum size below 0.8 μm (Fig. 12).

  Fig. 12 Flux ratio between the dust thermal emission and the scattered light integrated along the line of sight assuming our best-fit dust distribution and either porous astronomical silicates (left) or the more complex dust aggregates described in the text (right).

• 2.

  The Henyey-Greenstein scattering phase function we used does not reproduce the scattering properties of the dust grains well. The theoretical phase function calculated in the frame of the Mie theory, assuming hard silicate spheres of a few microns, predicts a much more peaked forward-scattering behaviour (see for instance Mulders et al. 2013). If we assume the flux within 1″ is mostly due to forward scattering events occurring around 100 AU, then the scattering phase angles are below 11°, thus well within the range where the Mie theory predicts a strong peak, and this could explain our missing flux. This may also explain why the excess emission in Fig. 10 (bottom) lies above the disc midplane. The scattering phase function is very sensitive to the shape and structure of the grains, however, and it remains to be proven that such theoretical phase functions are realistic for porous fluffy grains that can significantly deviate from a spherical shape (Augereau et al. 1999; Schneider et al. 2006).

• 3.

  Our forward modelling approach cannot constrain the dust distribution at short separation. The optimization region where the χ² is minimized starts at 0.4″ or 7.7 AU, so any change in the disc brightness or morphology within 0.4″ does not affect the χ². Pixels at 0.4″ also have a smaller weight in the χ² because the noise is higher at that separation, which could introduce a bias towards more external regions. This explanation alone cannot explain the missing flux observed up to 1″. An additional explanation might be the fact that the real and fake disc interfere during the star subtraction process, because the field rotation is close to 90° in some data sets. This would reinforce self-subtraction in the inner regions but both the fake and real disc would be similarly self-subtracted, whereas we do not notice a flux change in the real disc between Figs. 11 and 3.

Those three scenarios do not exclude each other. Based on the observational evidence mentioned above, especially the mid-infrared imaging, we favor the first scenario, but a combination of the three is very likely.

5.2. The reality of the ripples

We measured the departure of the spine with respect to the apparent midplane on the fake disc with the same technique as that applied to the real disc, and show the result in Fig. 13 (black solid line). On the same graph, we plot the measurement as performed on the real disc (dashed line), and performed on the initial model (red line).

First, the measurement performed on the fake disc injected at 90° also shows ripples with the same amplitude as those measured on the real disc. Therefore, the ripples detected longwards 1″ on each extension in Fig. 5 are likely to be caused by the remaining noise. Inside 1″, the fake disc spine measurements vary significantly before and after PCA reduction. In particular, we also detect a ripple on the left side at about 0.8″ with more than half the amplitude of that detected on the real disc. This is because of a very bright speckle located at a PA of about 110° (see Fig. 11, green arrow, also visible in Fig. 3), very close to the fake disc PA of 120.8°. This speckle clearly tends to shift the spine of the fake disc upwards, as visible in Fig. 11. We should emphasize that at this separation the reduced fake disc is 0.3 mag fainter than the real disc, so the measurement is more sensitive to the underlying speckle noise. Therefore, the current SNR does not allow us to rule out the fact that the deviation of the spine detected at 0.8″ might be due to the residual noise. This deviation is, however, twice as large as that induced by remaining speckles and is detected on three of the seven individual images (those with the best SNR), which may not be a coincidence.

Fig. 13 Departure of the spine with respect to the apparent midplane as measured on the initial disc model (top right image of Fig. 10), on the fake disc injected at 90° (Fig. 11) and on the real disc. NE is to the left for the real disc and SE is to the left for the fake disc.

Fig. 14 Comparison between the P97 and A01 models. The vertical dotted line indicates the separation of 115 AU beyond which the dust surface density inverted from the 12 μm observations was interpolated with power laws.

Second, we note some significant variations between the departure of the spine as measured on the initial disc model and on the real data. While the value is, on average, close to 0 and 2 AU for the NE and SW side, the slope is too steep for the SW side and not steep enough for the NE side. We will describe this behaviour in more detail in Sect. 5.4 with a two-component disc model.

5.3. Comparison with mid-infrared observations

In mid-infrared around 10 μm, thermal emission by dust is the dominant process over scattering in the inner 100 AU and observations are not limited by the contrast to the star while they still provide a resolution of 5 AU (Pantin et al. 1997; Okamoto et al. 2004; Telesco et al. 2005). It is therefore relevant to compare the dust distribution derived from our images to that derived from mid-infrared observations. Using observations at 12 μm from the TIMMI camera, Pantin et al. (1997) proposed a model of the dust distribution in the inner 100 AU, by inverting the integral equation relating the spatial flux to the dust surface density. Their choice of grain composition and size distribution led to the dust surface density displayed in Fig. 14 (plain red line), later referred to as P97. This profile corresponds to the average between NE and SW surface density presented in Fig. 5 of Pantin et al. (1997), after correcting for the revised distance of β Pictoris of 19.3 pc. Again, building a scattered light image using this distribution requires extrapolating the surface density beyond 100 AU, and we followed the same approach as described previously with different power laws of indices p_out = −1.5, − 1.7, and −2.3. The comparison between the P97 and A01 models is shown in Fig. 14. The dust surface density from P97 models predicts a narrower peak density and varies significantly from our models within the inner 20 AU, predicting more dust within 10 AU than our (A01, p_in = 0.5) and (A01, p_in = 1.5) models. We chose the same dust vertical distribution as our models, e.g. an exponential profile of index γ = 0.5, a flaring exponent β = 1.5, and a scale height ξ₀ at r_out = 115 AU, and run our minimization algorithm to subtract the best parametric model from our data. The results are shown in Table 4. The P97 model does not explain the data better than our models with a χ² degraded by 3.6% and 14% on average for an inclination of 88° and 86°, and the disc-subtracted image still displays some disc residuals within 1″. We can conclude that the A01 models better reproduce our observations beyond 0.5″ or 10 AU. Moreover, our estimate of the parameters PA, ξ₀ and | g | is robust to the underlying density distribution because we find similar values as in Table 3. Little can be said inside 0.4″ on the reality of the rise in dust surface density at about 8 AU in P97 (see Fig. 14) for two reasons: we are probably not dominated by scattering at this separation, especially if such a warm dust population is present, and the exact scattering properties of the dust at small phase angles are still uncertain.

5.4. Disc model with an inner and outer component

Fig. 15 Schematic of the disc model including a main and warped component. For clarity, the disc is shown edge-on, on a vertical cut. The red line indicates the region where most of the dust would lie in the A01, pin = 1.5 model.

Among the various disc parameters compatible with the data, it clearly appears that the PA differs significantly from the value measured in scattered light at larger separations beyond 5″. In an attempt to better understand this discrepancy and take advantage of our knowledge of the disc from previous modelling, we used a more sophisticated model to take the warped component of the disc into account. We therefore introduced a two-component model, as sketched in Fig. 15, with a main outer disc and a warped inner disc. We have too limited a field of view and lack a signal from the main outer disc to constrain the surface density of the outer component and lead the same forward modelling approach as for the single-component model presented in Sect. 4.2. This will be the purpose of a future analysis using multiple-wavelength data and scattered-light images with a larger field of view. Instead, we performed a qualitative analysis to understand if this more complex model can better explain the position of the spine. We assumed that the warped inner disc has a surface density identical to our A01, p_in = 1.5 model until the separation r_warp, and extrapolated with a power law of index p_out = −1.7 beyond r_warp. The main outer disc surface density is 0 inside r_warp and follows our A01, p_in = 1.5, p_out = −2.3 model beyond.

We further assumed that each of the two components can be described with the same parameters β = 1.5, γ = 0.5, inclination i = 86°, anisotropy coefficient | g | = 0.36 and scale height ξ₀ = 5.1 AU at 115 AU, taken from the best model of Table 3. The specific parameters of each components are:

• for the outer main disc, the inner radiusr_warp and the PA of 29°;

• for the inner warped disc, the outer radius r_warp, and the PA offset by φ = 4° with respect to the outer component.

In this model, the only free parameters are now the radius r_warp and the flux scaling factor. To constrain r_warp, we used the flux ratio between the main and the warped disc as measured from the two-component fit of vertical profiles by Golimowski et al. (2006) on the images from the ACS instrument aboard the Hubble Space Telescope. A ratio of 2.5 was measured at a separation of 80 AU. If we assume that this ratio does not depend on the wavelength, this translates into a value of r_warp of about 60 AU as shown in Fig. 16. Interestingly, Fig. 16 shows us that the main source of emission seen within 2.8″ comes from light scattered by the warped disc rather than by the main disc. The brightness ratio between the two components of the disc is indeed smaller than 1 inside 2.8″.

Fig. 16 Brightness ratio between the main outer disc and the warped inner disc for different values of rwarp. The dashed line indicates the ratio of 2.5 measured by Golimowski et al. (2006) at 80 AU or 4.15″. The vertical dotted lines indicate the different values of rwarp.

Using this model, we then performed our single-Lorentzian fit again after derotating the image by the PA of the main disc, i.e. 29° (Fig. 17). We also display the fit done on each individual components of the model. The black curve shows the same trend as our measured spine, namely an overall offset above the apparent midplane inside 2′′ with a larger offset on the NE than on the SW extension. We clearly see that the spine position of the double-component model lies in between that of the single-component models and tends to align with that of the main outer component beyond 3″. Within 3′′, the departure of the spine with respect to the apparent midplane is steeper on the NE side than on the SW side, as seen in our data in Fig. 7.

This model also produces a subtle brightness asymmetry of 0.1 mag/arcsec² between the NE and SW extension inside 2′′. Indeed, the spine position of the main and warped components (green and red curves in Fig. 17) being closer to each other on the SW side, the brightness of the summed contribution is slightly larger on the SW side than on the NE side. This difference is compatible with our measurements and can explain the slight asymmetry noted between 1.5′′ and 2′′ in Fig. 8.

Fig. 17 Departure of the spine with respect to the apparent midplane as measured on a disc model presenting an outer disc aligned with a PA of 29°, plus an inner warped component extending up to 60 AU. The black curve reproduces the behaviour of the measured departure of the spine well, as shown in Fig. 7.

All in all, this two-component model can explain the offset and bow of the disc. The presence of the warped component creates an asymmetry with a departure of the spine greater on the NE side than on the SW side, as observed in the data.

6. Conclusion

• 1.

  We have obtained L′-band images of β Pictoris that reveal the disc in its innermost regions inside 1″ dominated by light scattered from the warped inner component. The surface brightness distribution shows no overall asymmetry between the two extensions of the disc within error bars. The spine is offset and bowed due to anisotropic scattering combined with a small inclination of the disc estimated around 86°. A large ripple is detected in the spine position on the NE side at 0.8″ with a confidence level of 2σ, although we cannot rule out that it comes from the residual speckle noise. If real, this deviation of the spine might be the sign of a gravitational perturbation. We do not detect any obvious point-source within the disc after subtracting the signal from β Pictoris b, but two aspects should be highlighted. First, the brightness of the disc hinders the detection of faint point-sources. Second, as a result of our strategy to combine multi-epoch data, any planet on an orbit smaller than that of β Pictoris b would smear an arc on the combined image of the seven epochs and its flux would be diluted over more than one resolution element. A paper dedicated to the detection limits on potential planetary companions located in the disc midplane is in preparation. The paper will make use of these combined data sets to go deeper than in Absil et al. (2013) and will be complemented by radial velocity constraints.

• 2.

  The observations are compatible with the presence of dust within 10 AU from the star, although degeneracies in the modelling prevents us from accurately constraining the inner dust distribution. Our single-component disc model explains the overall morphology of the disc well. We favour a best model inclined by 86° with an anisotropic scattering coefficient | g | of 0.36. However, this model cannot explain the strong brightness of the innermost parts of the disc inside 1″. Different scenarios can be proposed to explain why our models fail to predict this strong brightness, the most likely being either that thermal emission starts to be significant inside 0.5″ or that enhanced forward scattering occurs for small phase angles. Additional observations in scattered light and at shorter wavelengths insensitive to thermal emission, but with a sufficient resolution and contrast to probe the inner 0.5″ are necessary to disentangle those two scenarios.

• 3.

  An investigation of a two-component disc model, mimicking a warped disc extending up to 60 AU and offset by 4° with respect to a main outer disc extending from 60 to 115 AU shows that within 2.8″ the scattered light mainly comes from the warped inner disc. The measured position of the spine is well reproduced by such a model, especially the bow and offset. A multi-wavelength modelling with a larger field of view is now needed to confirm our understanding of the dust distribution and bring additional constraints on the parameters derived in this study.

Acknowledgments

We would like to thank the referee Dr. A. Brandecker for his constructive comments that helped to significantly improve the quality of the article. J.M. acknowledges financial support from the ESO studentship program and the Labex OSUG 2020 and O.A. acknowledges financial support from a FNRS Research Associate during part of this work. Last, the research leading to these results also received funding from the European Research Council under the European Union’s Seventh Framework Programme (ERC Grant Agreement No. 337569) and from the French Community of Belgium through an ARC grant for Concerted Research Actions.

References

All Tables

All Figures

	Fig. 1 PCA-reduced image obtained after combination of the seven individual images, colour scale-adjusted to show the smearing of the planet (green arrow), which prevents disc signal extraction at this location. The green circle has a radius of 0.4′′ and is centred on the star (green square). The black arrow has a position angle of 29° and a total length of 6″.
In the text

	Fig. 2 Projected distance smeared by a potential disc feature or planet in keplerian motion around the star during the time span between the first and last observations (3.1 yrs).
In the text

	Fig. 3 Two different disc reductions: cADI-disc (top) and PCA-disc (bottom). The black arrow shows the apparent midplane, has a position angle of 29° (equal to that of the main disc), and a total length of 6′′. The colour scale is identical for both images but not linear (square root). The green square marks the position of the star and the green circle has a radius of 0.4″.
In the text

	Fig. 4 S/N map of the PCA-disc image. The colour scale is linear and spans from 0 to 10. The green circle indicates the start of the 5σ detection at a radius of 0.4″.
In the text

Fig. 5 PCA-disc image where each vertical profile has been normalized by the maximum spine brightness to enhance the vertical position of the spine. The two green circles have a radius of 0.4″ and 2″ respectively. The PA of the main disc (29°) as measured by Lagrange et al. (2012a) is indicated by a solid red line, whereas the best fit PA of 30.8° (NE) and 211.0° (SW) as measured from the L′ images are shown with dashed lines. The waviness of the spine is discussed in Sect. 5.2.

In the text

	Fig. 6 Isophotes corresponding to the cADI-disc image of Fig. 3. The black arrow passes through the star (green square), has a position angle of 29°, and a total length of 6″.
In the text

Fig. 7 Departure of the spine with respect to the apparent midplane (defined in Sect. 3.2 as the plane at 29° PA) as measured on the PCA-disc image. NE is to the left. The uncertainty, shown here in red, is the 3σ error in the centroid position of the Lorentzian profile fitted to each vertical profile, after setting a radially-dependent noise level to each pixel. The waviness of the curve is due to the background and remaining speckle noise, as shown in Sect. 5.2. We note, however, a larger ripple at 0.8″ on the NE side.

In the text

	Fig. 8 Surface brightness of the spine, measured on the cADI-disc image from the fit of a Lorentzian profile. The red and purple shaded areas show the 3σ error in the amplitude of the Lorentzian function fitted to each vertical profile.
In the text

Fig. 9 Description of the dust surface density used in the nine different models, based on the surface density derived by Augereau et al. (2001) for the parent bodies (PB, black line). Each of the nine models is made of three parts delimited by the two vertical dotted lines at 50 AU and 115 AU: an inner region (purple curves) following either a power law or the PB distribution, a common intermediate region following the PB distribution, and an outer region following a power law (green curves).

In the text

Fig. 10 From top to bottom: initial disc model, disc after PCA reduction without noise, residuals remaining after the disc model was subtracted from the data cubes and the cubes were re-reduced. The colour scale is linear and two times larger for the initial model than for the other two images. The green circles have a radius of 0.4″ and the black arrows indicate the apparent midplane (PA = 29°) with a total length of 6″.

In the text

	Fig. 11 Image obtained after injection of the best disc model at 90° from the real disc and reducing the data again. The field of view is 3.8″ and the real disc extends from NE to SW.
In the text

	Fig. 12 Flux ratio between the dust thermal emission and the scattered light integrated along the line of sight assuming our best-fit dust distribution and either porous astronomical silicates (left) or the more complex dust aggregates described in the text (right).
In the text

	Fig. 13 Departure of the spine with respect to the apparent midplane as measured on the initial disc model (top right image of Fig. 10), on the fake disc injected at 90° (Fig. 11) and on the real disc. NE is to the left for the real disc and SE is to the left for the fake disc.
In the text

	Fig. 14 Comparison between the P97 and A01 models. The vertical dotted line indicates the separation of 115 AU beyond which the dust surface density inverted from the 12 μm observations was interpolated with power laws.
In the text

	Fig. 15 Schematic of the disc model including a main and warped component. For clarity, the disc is shown edge-on, on a vertical cut. The red line indicates the region where most of the dust would lie in the A01, pin = 1.5 model.
In the text

	Fig. 16 Brightness ratio between the main outer disc and the warped inner disc for different values of rwarp. The dashed line indicates the ratio of 2.5 measured by Golimowski et al. (2006) at 80 AU or 4.15″. The vertical dotted lines indicate the different values of rwarp.
In the text

	Fig. 17 Departure of the spine with respect to the apparent midplane as measured on a disc model presenting an outer disc aligned with a PA of 29°, plus an inner warped component extending up to 60 AU. The black curve reproduces the behaviour of the measured departure of the spine well, as shown in Fig. 7.
In the text