© ESO 2019

1. Introduction

Due to the discovery of hundreds of debris discs by the Infrared Astronomical Satellite (IRAS; e.g. Aumann 1985) and thousands of exoplanets by the current large missions like Kepler (e.g. Borucki et al. 2010) we know that planetary systems are not rare in the solar neighbourhood.

There is growing evidence from the infrared and submillimetre surveys of the last few years that debris systems with multiple planetesimal belts around young main sequence stars are common as well (Carpenter et al. 2009; Moór et al. 2011; Morales et al. 2013; Chen et al. 2014; Jang-Condell et al. 2015). Their architectures could be similar to that of the solar system, harbouring two prominent rings of debris: the main asteroid belt, and the Edgeworth-Kuiper belt.

Several hypotheses were recently developed to explain the formation of such systems. The most favourable of them are disc sculpting by planets into distinct planetesimal belts (e.g. Dipierro et al. 2016; Geiler & Krivov 2017, and references therein) and the interaction of dust grains with the gas in the disc (e.g. Lyra & Kuchner 2013). The former mechanism assumes the presence of planets formed in the protoplanetary disc. In analogy to the solar system, the planets carve ring-like structures in the disc removing the gas and debris from their immediate vicinities, the so-called chaotic zones. They create a wide gap depleted of debris material which separates the regions of the disc stillcontaining dust and planetesimals. Evidence for this hypothesis could be provided by the planetary systems HR 8799 and HD 95086. In HR 8799, four planets (Marois et al. 2010) reside in the hole between the warm and cold dust belts (Su et al. 2009). A similar disc configuration with a warm and a cold component (Su et al. 2015) and a planet located in between them is detected in the HD 95086 system (Rameau et al. 2013).

Multiple ring structures can also originate from local instabilities induced by the gas in dusty discs. Using the hydrodynamical simulation of the gas–disc interaction, Lyra & Kuchner (2013) showed that the photoelectric heating of the gas by locally concentrated dust grains can create a pressure maximum resulting in a further dust concentration at this location. However, in order for this mechanism to work, a strong coupling of dust to the gas is required and therefore a disc dust-to-gas ratio lower than unity. This condition might be fulfilled in the so-called transitional discs such as HD 141569A (Thi et al. 2014) or in the gas-rich debris discs such as HIP 73145 (Moór et al. 2017) in which multiple ring structures were observed (Perrot et al. 2016; Konishi et al. 2016; Feldt et al. 2017).

The existence of double planetesimal belts is often inferred by analysing the stellar spectral energy distribution (SED). The SEDs of such systems show an infrared (IR) excess on top of the stellar photosphere which can be well fitted by a combination of two Planck’s laws with different temperatures. This implies that at least two spatially separated populations of dust exist in the system (e.g. Morales et al. 2011; Ballering et al. 2013): a so-called hot/warm component located closer than 10 au to the host star and a cold component usually between 70 and 150 au away.

Generally, the SED modelling can provide neither the actual number of belts in the system nor an accurate estimate of their location andspatial extent due to degeneracy between orbital distance and grain properties. Multiple cold belts, for instance, do not necessarily show up as distinct features in the SED. The actual distribution of the dust can only be revealed with direct imaging allowing high spatial resolution and high contrast (e.g. Matthews et al. 2017). Usually, the innermost component is inaccessible for direct imaging because it is faint and located too close to the star. To date, only a few examples of debris discs with multiple belts have been spatially resolved in scattered light; for example HD 92945 (Golimowski et al. 2011), HIP 73145 (Feldt et al. 2017), HIP 67497 (Bonnefoy et al. 2017) and HD 141569A (Perrot et al. 2016). In all these cases, the positions of the belts measured by direct imaging are different from those inferred from the SED modelling.

For a successful detection of multiple belts, a low inclination of the system is an advantage. All the above-mentioned objects have an inclination i lower than 80°. Obviously, it is much more difficult to identify more than one belt in edge-on discs. In this work, we investigate the presence of two belts of planetesimals in the well-known edge-on disc around HD 15115.

HD 15115 (HIP 11360) is a nearby zero-age main sequence F4IV star (Harlan 1974). According to the Gaia DR1 catalogue (Gaia Collaboration 2016), HD 15115 has a parallax angle of 20.76 ± 0.41 mas corresponding to a distance of 48.16 ± 0.95 pc which we apply in this work. This value is slightly different from the former determination of 45.2 ± 1.3 pc (van Leeuwen 2007). The age of the star is very uncertain with a range from 25 Myr, as a member of the β Pictoris moving group (BPMG; Mamajek & Bell 2014), although the BPMG membership of HD 15115 has been questioned (Debes et al. 2008), to 100–500 Myr, based on different age indicators (Debes et al. 2008).

The IR excess indicating thermal emission from circumstellar dust was first detected with IRAS (Silverstone 2000). Using the Hubble Space Telescope (HST) Advanced Camera for Surveys (ACS) and the Keck near-IR camera NIRC2, Kalas et al. (2007) subsequently discovered a highly asymmetrical debris disc in the optical and in the near-IR (H band). The disc was found to be nearly edge-on extending along the east-west direction with a west side detected out to a stellocentric radius of ≈580 au (≈12′′) and the east side to at least 340 au (≈7′′).

The complexity of the inner disc structure was noted by Debes et al. (2008) who imaged HD 15115 with the HST Near-Infrared and Multi-Object Spectrometer (NICMOS) at 1.1 μm and claimed a possible “warp” on the western side of the disc at around 1′′. They also established a wavelength dependence of the brightness asymmetries. Kalas et al. (2007) measured an overall blue V –H colour for the west side, which together with the shape of the disc gave it its nickname: “the blue needle”. Debes et al. (2008) detected a relatively intense variation of disc colours as a function of distance from the central star: at 2′′ they measured a neutral colour up to 1.1 μm and a blue spectrum for λ > 1.1 μm, while at 1′′, two disc extensions show a red F110W - H colour. This suggests the presence of several populations of grains with various properties, such as grain size or composition, at different separations. The red colour of the inner disc could be due to the strong absorption features of olivine at 1 μm (Debes et al. 2008). Rodigas et al. (2012) found, on the contrary, that the disc possesses mostly a spectrally neutral K_S –L′ colour (hence, “the grey needle”).

Analysing the data obtained with PISCES on the Large Binocular Telescope (LBT) and LBTI/LMIRcam, these latter authors showed that the east-west asymmetry is still present in the K_S band but the disc appears to be almost symmetrical in the L′ band. This symmetry was confirmed by Mazoyer et al. (2014) who detected a ring-shaped structure with a radius of R ≈ 1.99′′ and a position angle (PA) of 98.8° in H and K_S band images from the Gemini Near-Infrared Coronagraphic Imager (NICI) archival data sets. They also found that the disc has an inclination of ~ 86°, shows no noticeable centre offset with respect to the star, and probably has an inner cavity. This work provided the first unbiased view of the disc structure though with some discrepancies with the results of recent observationsof HD 15115 with HST (Schneider et al. 2014). Schneider et al. (2014) confirm the bow-like structure measured by Rodigas et al. (2012) and Mazoyer et al. (2014) but obtain a different inclination of ~ 80° for an inner “half-ring” while they estimate that the extended disc is edge-on (Schneider et al. 2014, see their Table 5). The authors also note a bifurcation of the scattered light away from the major axis on the east disc side and a significantly larger vertical extension of the disc on the west side compared to the east.

These puzzles could not be solved with the 1.3 mm observations using the Submillimeter Array (SMA; MacGregor et al. 2015). The resolved millimetre emission showed a rather symmetric distribution of large grains in the disc with a tentative detection of the asymmetric extension towards the west. MacGregor et al. (2015) fitted these data with a model of a broad planetesimal belt extending from R_in = 43 to R_out = 113 au.

Such a broad debris belt could be a combination of two or even more narrower planetesimal belts remaining unresolved in the SMA data. Previous studies have already classified the HD 15115 disc as a double belt system based on the SED modelling. Moór et al. (2011) used the disc fluxes at λ = 24 μm, λ = 70 μm, and λ = 160 μm (Table C.1) as well as the Spitzer Infrared Spectrograph (IRS) 5–35 μm spectrum for their fit and found two belts with radii R_{warm dust} = 4 au and R_{cold dust} = 42 au. Chen et al. (2014) derived R_{warm dust} = 5 au and R_{cold dust} = 175 au based on the fits to the same Spitzer IRS spectrum and the IR excess at λ = 13 μm, λ = 24 μm, and λ = 31 μm.

In this paper, we present new observations of HD 15115 in the broad bands J and H carried out with SPHERE (Spectro-Polarimetric High-contrast Exoplanet REsearch; Beuzit et al. 2008) at the Very Large Telescope (VLT) in Chile. The paper is organised as follows. In Sect. 2 we describe the observations and present the data. Section 3 is dedicated to the methods of data reduction and Sect. 4 to the data analysis. The modelling of the disc geometry using both the total intensity and the polarimetric data is presented in Sect. 5. Section 6 describes the measurement of stellar and disc fluxes and compares the results from different datasets. In Sect. 7, we discuss our data analysis and possible indications for multiple belts in the HD 15115 system from earlier observations and finally summarise our results in Sect. 8.

2. Observations

The HD 15115 observations (098.C-0686 and 096.C-0640) were performed with the Infra-Red Dual-beam Imager and Spectrograph (IRDIS, Dohlen et al. 2006), the near-infrared imaging sub-instrument of SPHERE. We used the pupil-stabilised mode of IRDIS in order to allow for angular differential imaging (ADI; Marois et al. 2006)while the field of view rotates. The observations of HD 15115 with the broad-band filter BB_H (λ_c = 1626 nm, Δλ = 291 nm, hereafter H band) in both channels of the instrument were conducted on Oct 29, 2015, and with the broad-band filter BB_J (λ_c = 1258 nm, Δλ = 197 nm, hereafter J band) on Nov 19, 2015. To block the stellar light we used an apodized Lyot coronagraph (N_ALC_YJH_S) including a focal mask (ALC2) with a diameter of 185 mas (inner working angle of 95 mas, Boccaletti et al. 2018). The detector integration time (DIT) and the total exposure time together with the observing conditions are summarised in Table 1. To measure thestellar flux, several non-coronagraphic frames were recorded before and after each observation with DIT ≈ 0.84 s in H band and DIT = 2.00 s in J band using a neutral density filter ND2.0 with a ~ 10⁻³–10⁻² transmissivity to prevent the saturation of the detector.

At the beginning of the science observation a “centre frame” was taken using the deformable mirror waffle mode (Langlois et al. 2013). This frame provides the accurate measurement of the star position when using the coronagraph. The stability of the star position during the observation is ensured by the Differential Tip-Tilt Sensor (DTTS, Baudoz et al. 2010).

For the polarimetric measurements on Oct 3, 2016, we used the dual-polarization imaging mode (DPI) of IRDIS (Langlois et al. 2010), for which the field of view is stabilised. The data were taken in the J band with the same apodized Lyot coronagraph as in the imaging mode. We performed ten QU cycles consisting of consecutive measurements of Stokes parameters Q, −Q, U, −U setting DIT to 32 s per frame. The total exposure time of the polarimetric observations was 5120 s. As for the first observations, the stellar flux in the J band is obtained with an off-axis image of the star using a neutral density filter.

A raw frame with imaging or polarimetric data recorded with IRDIS has a format of 2048 × 1024 pixels. In Classical Imaging (CI), raw frames contain two images side-by-side (left and right channel) simultaneously acquired with a broad-band filter in the common path. In DPI mode, the left and right images contain the intensities of orthogonal linear polarizations, also taken simultaneously with a common broad-band filter. The common optical path of both channels ensures that the two beams have similar optical aberrations allowing an efficient suppression of speckles within the AO (Adaptive Optics) control radius.

The pixel scale of the IRDIS detector is 12.25 mas and the field of view (FOV) is 11′′ × 11′′.

3. Data reduction

3.1. Classical imaging data

To produce the calibration data, raw frames (flats, darks and sky background) are processed with the SPHERE data-reduction and handling (DRH) pipeline esorex (Pavlov et al. 2008). A custom routine is used to perform the basic calibration of the science frames: a bad pixels removal, the subtraction of the sky (instrumental) background, and the flat-field correction. The calibrated science frames with the size of 2048 × 1024 pixels are then split in two halves to create two temporal data cubes (the left and right channels of the instrument) containing 192 (H band) or 384 (J band) frames of 1024 × 1024 pixels each.

The position of the star is determined by fitting a 2D Gaussian function to the four satellite spots in the “centre frame”. The star coordinates are derived from the intersection point of two lines connecting the centres of two opposite spots. These coordinates are used for re-centring all frames in the data cube.

In the next step, reduced and centred data cubes are processed with several custom algorithms making use of ADI technique: cADI (classical ADI; Marois et al. 2008), KLIP (Karhunen-Loève Image Projection; Soummer et al. 2012), and LOCI (Locally Optimized Combination of Images; Lafrenière et al. 2007).

The final median-combined images obtained with cADI, after the stellar light had been removed, are shown in Figs. 1a and d.

Similar results were obtained with the KLIP algorithm (Figs. 1b and e). This method uses a principal component analysis to construct a reference image for the stellar light from the projection of the science image onto an orthogonal basis given by the Karhunen-Loève (KL) transform of the same science images. We performed KLIP data reductions with 3, 5, 15, and 25 KL vectors (or eigenimages) forming the basis for a reference image. The entire images were used without a frame selection to avoid large disc-flux self-subtraction effects. Visual inspection of the final images showed that the larger number of vectors yields a very strong self-subtraction of the disc flux.

For comparison, we also performed the LOCI data reduction (Figs. 1c and f). The LOCI algorithm provides a high detectability performance similar to the KLIP procedure. The reference image is constructed for each frame individually as a linear combination of other images in the data cube which have a minimal difference in the parallactic angle to the processed frame. The requested minimum difference in field rotation, in units of arc length, is set by parameter N_δ . The optimisation algorithm works in annular segments with the size specified by the parameter N_A which gives the number of “PSF cores” that should fit in the segment area. The frames must cover a large rotation angle to preserve the flux of an extended structure like an edge-on debris disc. For the LOCI reduction we adopted parameters N_δ = 1.25 × FWHM and N_A equal to 1000 and 300.

To correct for the IRDIS True North (TN) offset, all frames are additionally rotated clockwise by 1.7° (Maire et al. 2016) during the ADI procedure.

Reduced frames of the stellar flux measurement are centred by fitting a 2D Gaussian function with a FWHM of ~ 3 pixels to the stellar profile. The sky subtraction is not applied to these frames. Instead, we measure the instrumental background in the median frames of both IRDIS channels and subtract it. The median frames are then divided by the mean transmissivity of the neutral density filter ND2.0 equal to 3.3 × 10⁻³ and 1.5 × 10⁻² in the J and H bands, respectively (see SPHERE manual¹). The final frames that are used to measure the stellar count rate (see Sect. 6) are the mean frames of two flux measurements performed before and after each observation.

3.2. Polarimetric data

We used theDRH pipeline esorex to produce the calibration data. The raw polarimetric science frames are calibrated with the same custom procedure (subtraction of the sky or instrumental background, flat-field correction, bad pixel removal) as applied to the total intensity data.

The Stokes Q and U parameters were obtained with the double-difference method from calibrated data using a custom routine. We applied the same basic data-reduction steps as described in Engler et al. (2017) where the procedure is explained in detail. For the frame centring we have used the same method as for the centring of the total intensity data set (see Sect. 3.1). After that, each Q and U frame was visuallyinspected to check the correct position of the centre of the right channel frame relative to the left channel frame.

The final Stokes Q and U parameters are calculated as follows:

where Q⁺ is measured along the axis pointing to the north, U⁺ is measured along the axis pointing to the northeast, and Q⁻ and U⁻ are measured along the axes perpendicular to Q⁺ and U⁺, respectively.

The Q and U images contain the linear polarization components of the total polarized flux received from the debris disc. To generate the image of the disc showing the total polarized flux, the instrumental coordinate system is transformed to the polar coordinate system where the linear polarization components are represented by the azimuthally oriented Stokes Q_φ parameter and U_φ parameter, the linear polarization component at a 45° angle to the Q_φ parameter:

where φ is the polar angle measured east of north in the coordinate system centred on the star. Here we adopt the sign convention for Q_φ and U_φ parameters which is consistent with Q_r and U_r parameters defined in Schmid et al. (2006).

The instrumental polarization signal is corrected through a forced normalization of the total counts in frames of opposite polarization states.

Finally, the Q_φ and U_φ images are rotated clockwise by1.7° to perform the astrometric TN correction (see Sect. 3.1).

Figure 2 shows the Q_φ image with the detected polarized light. The debris disc appears as a thin almost symmetric stripe oriented in the east-west direction and extending beyond 2.5′′ (120 au). There is no detection in the U_φ image, as expected. Therefore the Q_φ image shows the total intensity of the polarized light.

4. Data analysis

4.1. Disc morphology

4.1.1. Total intensity data

All data reductions (cADI, KLIP, LOCI) reveal an almost edge-on disc with the bright side to the north of the disc major axis. The famous bow-like structure, previously reported in different studies (Rodigas et al. 2012; Mazoyer et al. 2014; Schneider et al. 2014) and corresponding to the brightest side of the highly inclined ring, is clearly visible in both J and H. Additionally, our images unveil the disc ansae symmetrically located to the east and west of the star. In some images, parts of the fainter side of the disc are also visible near the ansae. The brighter disc side is detected with a signal-to-noise ratio (S/N) higher than 4 (see S/N map in Fig. A.1) up to a distance of ~ 4′′ eastwards and 6′′ on the west side (up to the end of frame).

4.1.2. Polarized intensity data

Figure 3 shows the binned Q_φ image after binning by 10 and 2 pixels along the x- and y-axis, respectively. The smaller vertical binning is applied to preserve the spatial information provided by the high-resolution imaging. In this figure we define an x–y coordinate system (Fig. 3a) where the star is at the origin, + x is the coordinate along the major axis westwards (θ_disc), − x is eastwards (θ_disc − 180°), and the y-axis is perpendicular to this with the positive axis towards the north.

There are three surface brightness (SB) peaks on the east side within r = 1.5′′ in azimuthal direction: at A(−0.80′′, 0.29′′), B(−1.30′′, 0.05′′), and C(−0.75′′, − 0.16′′) indicated in Fig. 3b. However, the points A and C are inside the AO correction area dominated by strong stellar residuals and cannot be considered as a detection. A straight line going through point B and the star intersects the small region with increased SB on the west side shown by arrow D. This region is mirror-symmetrically located to the point B with respect to the star.

Such a distribution of the polarized SB conforms with the highly inclined double belt system in which the outer and inner belts have slightly different PA and inclination from our toy model of the system. Figure 4 shows the synthetic images of the polarized intensity in a debris disc with single and double belt in comparison to the data. To create synthetic images, the model² with radius of the outer disc R₀ = 2′′ and inclination of 86° was used. The other parameters were arbitrarily chosen, for instance, the radius of the inner belt was set to 1.5′′ . In order to have a similar level and pattern of noise in the simulated images and in the data, both models (single and double belt) were inserted into the Q_φ image. The disc flux was previously removed from this Q_φ image. This was done by replacement of the image area containing disc flux through the same area of the mean of two images obtained by rotation of the Q_φ image by 15° clockwise (first image) and by 15° counterclockwise (second image). The model images were binned and smoothed in the same way as was done to obtain the Q_φ image presented in Fig. 3. The visual comparison between model images (Figs. 4a and c) and the Q_φ image (Fig. 4b) suggests that the HD 15115 debris disc could be a double-ring system because the SB distribution observed in particular on the east side of the star can be explained by the presence of a second inner belt. We also note that the single-belt model (Fig. 4c) lacks polarized flux inside 1.5′′.

Motivated by this comparison, we fitted an ellipse to the SB peaks A, B, C, and D and found that, if such an inner belt were to exist, it would possibly have a deprojected radius of r ≈ 1.25′′ ± 0.07′′ and a PA θ_{inn ring} = 276° ± 3°, and thus a difference in PA up to 6° with respectto the outer belt.

The goal of the data analysis presented in the following sections is to find out whether there is some additional evidence for the double-belt structure suspected from the polarimetric observation. For this purpose, we explore the position of the disc spine (Sect. 4.3), compare themodels for the spatial distribution of dust in an optically thin disc with the total and polarized intensity data (Sect. 5), and investigate the variation of the disc SB along the disc axis calculated in apertures of different sizes (Sect. 6).

Fig. 1. Total intensity images of the HD 15115 debris disc obtained with different data-reduction techniques in the H band (left column) and in the J band (right column). The number of modes in KLIP reductions is 10. The star position is marked with a yellow asterisk. The FOV of each displayed image is 8.2′′ × 3.1′′.

Fig. 2. Polarized intensity images of the HD 15115 debris disc in J band. Images show the Stokes Qφ parameter smoothed by a Gaussian kernel with σ = 2 pixels (panela) and with an additional 4 × 4 binning (panel b). The white arrow points to the local surface brightness peak responsible for the impression of the inner belt (see also Fig. 3). The star position is marked with a yellow asterisk. The colour scale shows the flux in counts per binned pixel. The surface brightness contours in panel b are given for 5 and 15 counts per binned pixel. The FOV of each displayed image is 8.2′′ × 3.1′′.

Fig. 3. Qφ images binned by 10 pixels along the x-axis and by 2 pixels along the y-axis and smoothed via a Gaussian kernel with σ = 1 pixel (panel a) and with σ = 2 pixels (panel b) to highlight the structure and position of both rings relative to each other. The central part of the images, dominated by strong residuals, is masked out. The position of the star is indicated by a yellow asterisk. The x-axis (interrupted white line going through the star position) shows the outer disc major axis at PA = 279.0°. The FOV of the displayed panel a is 9.8′′× 0.74′′. The letters “A”, “B”, “C” and “D” are used to label the surface brightness peaks discussed in Sect. 4.1. Colour bar shows flux in counts per binned pixel.

Fig. 4. Comparison of the double (column a) and single (column c) belt models with the Qφ image (column b). All images are binned by 10 pixels along the x-axis and by 2 pixels along the y-axis and additionally smoothed via a Gaussian kernel with σ = 2 pixels (top row) and with σ = 1 pixel (bottom row).

4.2. Position angle and size of the disc

We estimated the disc radius, the PA, and the inclination from the H band imaging data reduced with LOCI and KLIP algorithms where the disc front and back sides are visible. We derived the R_disc = 2.01′′ ± 0.05′′, θ_disc = 278.7° ± 0.2° and the disc inclination i = 86.5° ± 0.5° by fitting an ellipse to the ring contours and the SB peaks at the east and west disc edges.

We also derived the PA and the radius of the disc by determining the centroids of bright ansae at the radial separation r ≈ 2′′ to the east and west of the star in the polarimetric data in J band (Fig. B.1). This location corresponds to the radius of the symmetrical ring detected by Mazoyer et al. (2014) in the H band data which is evident in our data as well. For an edge-on disc or ring with a characteristic radius R_disc, the maximum polarized flux is expected at the locations of the planetesimal belt or at scattering angles close to 90°. This is a natural outcome of the interplay between the angular distribution of the polarized flux, which is described by the polarimetric phase function, and the radial distribution of the grain number density, which is maximal at R_disc (see also Sect. 5.3). This provides a relatively simple method to directly measure the extent of the planetesimal belt and its position on the sky. Using this approach we obtain for the HD 15115 outer disc the radius R_disc = 1.96′′ ± 0.05′′ (94.39 ± 2.41 au) along the disc major axis and the PA θ_disc = 279.0° ± 0.2°. After taking into account the PA measured in the total intensity data, this yields a mean PA of the disc θ_disc = 278.9° ± 0.1°.

These values are in good agreement with the disc radius R_disc = 1.99′′, PA θ_disc = 278.8° ± 0.4° and the disc inclination i = 86.2° found by Mazoyer et al. (2014) in the NICI data. Our measurement of the PA agrees also very well with results from previous studies of HD 15115 at different wavelengths (Kalas et al. 2007; Schneider et al. 2014). We do not see any evidence in our data of an east/west misalignment as noted by Debes et al. (2008) who obtained PA θ_disc = 278.9° ± 0.2° for the west side and PA θ_disc = 276.6° ± 1.3° for the east side at 1.1 μm. It should be noted that the NICMOS data have much lower resolution compared to the IRDIS data and the PA measurement of Debes et al. (2008) could be affected by the strong stellar residuals around the coronagraph.

4.3. Spine of the disc

The geometric parameters of the disc can be inferred from the spine which we define as the location of the maximum intensity of a Moffat function fitted to the disc cross-sections perpendicular to the midplane. The Moffat function is given by the following equation (Trujillo et al. 2001):

where a_M is the peak flux located at a vertical distance y₀ from the disc major axis and α and β are two constant parameters determined from the fit.

We measured the spine position of the brighter side of the disc for several data reductions (cADI, KLIP, LOCI) in H and J bands (with a 4 × 4 binning) and derived the mean vertical offset as plotted in Fig. 5. The y₀ values, which are larger than y₀ = 100 mas or smaller than y₀ = −25 mas, are considered as poor fits due to a low S/N and are excluded from this analysis. The error bars show the standard deviation of the offset values obtained in all data sets without taking these outliers into account.

The spine has a well defined bow-like shape as already described in Mazoyer et al. (2014; see their Fig. 5) but here with a much smaller dispersion. The smaller error bars are due to the improved angular resolution and contrast provided by SPHERE compared to NICI. We also note that the spine dispersion is relatively small for the radial separations between 1.2′′ and 3.7′′ on the west side and between 1.0′′ and 2.4′′ on the east side. Between r ≈ 1.2′′ and r ≈ 2.2′′ the spine is roughly symmetric with respect to x₀. The smallest -offset is approximately 45 ± 5 mas (2.2 au) around |x| = 2.0′′ ± 0.1′′. For smaller radial separations the vertical offset is increasing and achieves a maximum of ~ 140 mas at x ≈ ±1.0′′ in our data. For radial separations x > 2.0′′ the spine deviates from the major axis by 45–60 mas towards the north. The same behaviour is observed on the east side, although with a lower S/N. A similar spine curve which does not reach the disc major axis at the location of the planetesimal belt was measured for another edge-on debris disc around F star HIP 79977 (Engler et al. 2017) and A star β Pictoris (Milli et al. 2014). Such a course of the spine curve is expected (see Sect. 5.4) and explained by the anisotropic scattering of the dust (Mazoyer et al. 2014).

Fig. 5. Spine vertical distance from the disc major axis measured as a vertical offset of the Moffat profile centre. The plot traces the mean location of the spine in all available images (cADI, KLIP, LOCI reductions of the total intensity data) in H and J bands. The star position is at (0, 0). The black line shows the best-fit ellipse to the data points between |x| ≥ 1.2′′ and |x| ≤ 2.0′′.

Spine diagnostics I

The turning point of the spine at ~ 2′′ clearly indicates the location of the edge of the parent belt of the planetesimals. The geometry of the disc bow traced by the accurately measured spine allows us to put constraints not only on the characteristic radius of the belt but also on the disc inclination. We derive these parameters by fitting an ellipse to the spine between |x| ≥ 1.2′′ and |x| ≤ 2.0′′. The centre of the ellipse is set to the star position. The best-fit ellipse (Fig. 5) has a semi-major axis of 2.06′′ ± 0.06′′ and a semi-minor axis of 0.15′′ ± 0.01′′ yielding a disc inclination of 85.8° ± 0.7°. This inclination estimate is in good agreement with i = 86.2° found by Mazoyer et al. (2014) with the same methodology.

5. Disc modelling

5.1. One-belt model

To put constraints on the basic geometric parameters of the disc, we generate a set of models from a parameter grid and compare the synthetic images for the scattered light with the KLIP disc images in J and H bands (Figs. 1b and e) by evaluating the reduced for the model fit (Eq. 7).

The synthetic image of scattered light is created for each model assuming a rotationally symmetric distribution of the dust. In accordance with the theory of the “birth” ring, small dust grains are generated in a collisional cascade triggered by large planetesimals residing in a narrow ring at a distance R₀ from the star. Therefore, the radial distribution of the grain number density has a peak at the location of the parent belt and can be described by the power laws r^α_in (α_in > 0) for r < R₀ and r^α_out (α_out < 0) for r > R₀ to account for the decrease of the grain number density inside and outside of the ring (e.g. Augereau et al. 2001). Adopting a Gaussian function for a vertical distribution, the grain number density distribution n(r, h) in the disc can be written in cylindrical coordinates (r, h) as: (5)

where H(r) is the disc scale height which scales with a radius like , H₀ is the scale height of the vertical profile of the disc at R₀, and H₀∕R₀ defines the aspect ratio of the disc. For simplicity, we set the disc flare index to γ = 1.

Furthermore, we assume that the disc is optically thin, that is, that the single scattering dominates over multiple scattering of photons, and adopt the Henyey-Greenstein (HG) phase function for the angular dependence of the scattered radiation (Henyey & Greenstein 1941): (6)

where parameter g is the average of the cosine of the scattering angle for the phase function. There is g = 0 for the isotropic scattering, forward scattering grains have 0 < g ≤ 1, while for − 1 ≤ g < 0 the scattering is peaked backwards. In this work,we do not consider g < 0 assuming only forward scattering for the dust grains. We do not introduce the size distribution of dust grains to calculate the asymmetry parameter g and the cross section for scattering σ_sca averaged over all particle sizes (Engler et al. 2017). In our model, these parameters are constant through the disc, which is a reasonable approximation if we take into account that disc images preferentially probe grain sizes comparable to the wavelength of observation.

The parameter space investigated in this work is based on the results of a previous modelling of the HD 15115 disc and includes the following parameter sets (a grid of 2160 models in total):

• power-law index α_in: 2, 3, 4, 5, 6, 7, 8, 9, 10;

• power-law index α_out: -2, -3, -4, -5, -6, -7, -8;

• asymmetry parameter g: 0.2, 0.3, 0.4, 0.5, 0.6;

• disc aspect ratio H₀∕R₀: 0.005, 0.0075, 0.01, 0.02, 0.03, 0.05.

The radius and inclination of the disc in all models are set to 2′′ (96 au) and 86°, respectively.

Each model is first rotated and inserted into an empty data cube to mimic the actual position of the disc in the data. Then each frame of this cube is convolved with the instrumental PSF measured from the data and the cube is processed through the KLIP algorithm. Thereby, we use the forward-modelling approach to account for the effects of the disc self-subtraction which always arise when the ADI technique is used (Soummer et al. 2012). In forward modelling, the synthetic images of models are projected on the same vector basis as determined during the processing of the original data cube. Figure 6 shows an example of a synthetic model image (Fig. 6a) convolved with the PSF (Fig. 6b) and then processed through the KLIP algorithm (Fig. 6c) using an identical vector basis which was applied to obtain the KLIP images presented in Figs. 1b and e. The self-subtraction of the disc flux caused by the KLIP post-processing is clearly visible in the form of dark stripes on both sides of the disc in the data and the models.

To estimate the goodness of the model fit to the KLIP images, we define the reduced metric as follows: (7)

where N_data is the number of binned pixels used to perform the minimization, F_{i, data} is the flux measured in pixel i with an uncertainty σ_{i, data}, F_{i, model}(p) is the modelled flux of the i pixel and A is a scaling factor. The degree of freedom of the fit is denoted by ν = N_data − N_par, where N_par represents the number of variable parameters which is equal to five in this case.

A metric is evaluated for each KLIP-processed model of the grid with respect to the data presented in Figs. 1b and e. We applied a 4 × 4-pixel binning of the data and the models. Because of the pronounced brightness asymmetry between both disc sides, the west and east extensions were fitted separately in the first step. We performed the minimization in rectangular areas that are different for the west side (1′′ < x < 4′′ and − 0.2′′ < y < 0.2′′) and for the east side ( − 3′′ < x < −1′′ and − 0.2′′ < y < 0.2′′) for reasons related to the S/N. We also considered both disc sides simultaneously using the same areas. The error in intensity is calculated as the standard deviation of the flux distribution in binned pixels in concentric annuli excluding image areas which contain disc flux.

We consider the model providing the minimum as a best-fit model. Additionally, we determine the parameter sets of all well-fitting models bounded by a threshold of (Thalmann et al. 2014) to take the systematics into account. The threshold is derived with corresponding to a 1σ deviation of the χ² distribution with ν degrees of freedom. Tables 2 and 3 provide the best-fit model parameters (Cols. 2, 4, 6) as well as the range of parameters of well-fitting models satisfying the threshold noted in the column headers (Cols. 3, 5, 7). The parameter values of the best-fitting models can differ significantly between the east and west sides. The flux distribution on the west side can be well described with the HG asymmetry parameter g = 0.3 and the relatively shallow radial power-law indices α_in = 3 and α_out = −4 (H band data) or α_out = −3 (J band data), whereas the east side is best represented with g = 0.5 and a much steeper index α_in = 8 (H band data) and α_out = −6. In the J band data, there is no large difference for the inner radial index α_in between the east and the west side. The aspect ratio of 0.03 between disc radius and width seems to give a good fit to both disc sides, in particular when modelling J band data. As a final remark on the modelling results, we note that the χ² -minimization method is based on the assumptions that the flux measurement errors have a normal Gaussian distribution and the pixel values are independent. These assumptions might not be warranted in the considered data set and some other well-fitting parameter combinations might exist outside the quoted ranges.

5.2. Two-belts model

The analysis of polarimetric data (Sect. 4.1) suggests that the observed morphology of the HD 15115 disc could be explained by the presence of an inner belt of planetesimals with a radius of R < 2′′. If such an inner component exists it is not showing as an obvious feature in the total intensity data in particular because of the high inclination of the HD15115 disc. The phase function in total intensity being different than in polarimetry, an inner component will dominate the disc image at very small phase angles precisely where the ADI post-processing strongly attenuates the disc signal. To further test the disc geometry we consider both the one- and two-belts models and compare their KLIP post-processed signatures.

For this purpose, we use the one-belt model, which gives a reasonable fit to the H band data ( for the west side and for the east side). The model is generated with the following parameters: R₀ = 2′′, i = 86°, g = 0.3, α_in = 5, α_out = −4, H₀ ∕R₀ = 0.01 and PA θ_{out belt} = 279°. Figure 6 shows the appearance of the modelled belt (Fig. 6a), the effect of the blurring instrumental PSF (Fig. 6b), and the KLIP processing (Fig. 6c) on it.

For ease of comparison, when modelling the two-belts configuration, we use the same parameters for the outer ring as described in the previous paragraph. The radius of the inner belt is set to check text R_{inn belt} = 1.3′′ in accordance with the location of the SB peaks (see Sect. 6.3). For the remaining parameters of the inner belt, we adopt the mean values of the investigated parameter space: the asymmetry parameter g = 0.4, radial power-law indices α_in = 5 and α_out = −5, and the aspect ratio of 0.01. The intensity ratio between two belts is another free parameter in the two-belts model. For simplicity, we set this parameter to 1, assuming that the grain number density in the outer belt is higher by a factor of to take into account the light dilution with the distance from the star. For the analysis presented in this section, we consider two models, each consisting of two belts: a non-coplanar model where the inner belt has a different inclination and PA from the outer belt, and a coplanar model where the inner belt has the same inclination and PA as the outer belt.

In a non-coplanar model, the outer belt (PA θ_{out belt} = 279°) has an inclination of i = 86°; the inner belt has a PA of θ_{inn belt} = 276° as measured from the polarimetric data (see Sect. 4.1) and an inclination of i = 80°. Such inclination was estimated by Schneider et al. (2014) for the inner “half-ring” component of the HD 15115 disc (see note to their Table 5). The inner ring is shown in the top panel of Fig. 7. The other panels of this figure present the non-coplanar model before convolution with the instrumental PSF (Fig. 7b) and after convolution (Fig. 7c), and the image obtained with the KLIP algorithm (Fig. 7d).

The effect of a misalignment of two belts is clearly visible in Figs. 7b and c. Because the relative PA and inclination of belts are non-zero, we expect the maximum SB of the disc to be shifted to the west side (to the right in this figure). The PSF convolution tends to blur the separation between the two belts (Fig. 7c), which however, is reinforced afterwards (but also attenuated in intensity) owing to the high pass filtering produced by the KLIP processing (Fig. 7d).

Figure 8 shows the coplanar two-belts model before convolution with the PSF (Fig. 8a) and after KLIP processing (Fig. 8b). The disc morphology in Fig. 8a resembles that of a single ring (cf. Fig. 6a), as expected. The two-belts model is slightly brighter in the central part and appears to be thicker in the KLIP-processed image (Fig. 8b) compared to the one-belt model (cf. Fig. 6c).

The comparison between KLIP images of all models (Figs. 6c, 7d and 8b) illustrates the challenge of detection of multiple components in an edge-on disc. To see the split of the belts in a KLIP-reduced image, we need to resolve the system down to very small angular separations from the star.

In order to test how well the models with two belts fit the H band data, we calculated the reduced for the models presented in Figs. 7d and 8b setting the intensity ratio between theinner and outer belt equal to 0.5, 1, 2, 3 and 4. A simultaneous fit of both disc sides yields the in the range3.38 and 3.67 for the non-coplanar model and in the range 3.68 and 3.89 for the coplanar model (ν = 1060, N_par = 10). The minimum is obtained forthe non-coplanar model with an intensity ratio of 3. This is the smallest value of all tested models (including single belt models)indicating a higher conformity with the data. The number of free parameters in this model is however larger than in the one-beltmodels, and it is possible that this is a result of data over-fitting. Therefore, we compare the best one-belt model (see Col. 6 in Table 2 for the model parameters) with the non-coplanar two-belts model usingthe Bayesian information criterion (BIC) to check whether the additional parameters are warranted. This criterion introduces a penalty term for the additional parameters which increase the likelihood function of the model . We define the BIC as follows (Schwarz 1978): (8)

where is a set of the best-fitting parameters for the model providing the maximized value of the likelihood function . Under the assumption that the model errors are normally distributed, this function can be estimated as: (9)

The model with the lowest BIC is preferred. We estimate the BIC = 3901+ const for the model with one belt and BIC = 3653+ const for the non-coplanar two-belts model. The latter BIC is smaller, leading us to conclude that the double belt model provides a better fit to the data; additional parameters are justified. The change in BIC of 248 is greater than 10, which is in support of the model with two belts (Raftery 1995).

For all other variants of the fitting procedure (each side separately or using the J band data), the considered models with two belts are among the 20 best of all the models we tested.

The visual comparison between the best-fitting one-belt model, the non-coplanar and coplanar two-belts models, and the data is shown in Fig. 9. The models are inserted into the KLIP-processed H band image (Fig. 9a) where the disc flux was removed in the same way as was done with the Q_φ image (see Sect. 4.1). Because all models are rotationally symmetric, their intensity is adjusted to match the intensity of the west disc side for an easier comparison with the asymmetric HD 15115 disc. The one-belt model (Fig. 9b) and the non-coplanar two-belts model (Fig. 9c) are comparable with the data (Fig. 9a). There is no significant visual difference between these two models, they both reveal the back side of the disc, for instance. On the contrary, the two-belts coplanar model (Fig. 9d) does not show the disc opening (back side is invisible).

We conclude that the two-belts model with the inner ring parameters R_{inn belt} = 1.3′′, i = 80° and PA θ_{inn belt} = 276° provides a better fit to the H band data than the coplanar model or single-belt model according to the method andthe BIC. However, we have to keep in mind that this conclusion is based on the χ² calculation and is therefore subject to the limitations mentioned in Sect. 5. This means that none of the considered models can be absolutely ruled out. Higher spatial resolution and higher S/N of the data, especially at small angular separations, are necessary to differentiate between these models.

Fig. 6. Example of synthetic images of one-belt model. Panel a: model before convolution with the instrumental PSF. Panel b: model convolved with the instrumental PSF. Panel c: model after processing with the KLIP algorithm. The flux is shown in arbitrary units.

Fig. 7. Synthetic images of the two-belts model with an inner belt (shown separately in panel a) which is non-coplanar with the outer belt. The outer belt has the same geometry and asymmetry parameter g as the one-belt model demonstrated in Fig. 6. Panels b and c: two-belts model before and after convolution with the instrumental PSF, respectively. Panel d: KLIP-processed model. The flux is shown in arbitrary units.

Fig. 8. Synthetic images of the two-belts model with an inner belt coplanar with the outer belt. The model parameters are the same as in the model shown in Fig. 7 except for the inclination and PA of the inner belt which are the same as for the outer belt. Panel a: model before convolution with the PSF. Panel b: model after processing with the KLIP algorithm. The flux is shown in arbitrary units.

Fig. 9. Comparison of the HD 15115 disc (panel a) with the models inserted in the image (H band, KLIP reduction) where the disc wasremoved. Panel b: model with one belt. Panel c: two-belts model with an inner belt which has a different PA and inclination with respect to the outer belt. Panel d: two-belts model with an inner ring coplanar with the outer ring. In the processed images of the one-belt model (panel b) and non-coplanar two-belts model (panel c) the back side of the disc is visible, whereas it is not recognisable in the image of the coplanar two-belts model (panel d). Colour bar shows flux in arbitrary units.

5.3. Polarized light model

Since we have polarimetric data in the J band, we would like to test whether the two-belts model conforms with the polarimetric data set.

Polarimetry is particularly useful to study the structure of an edge-on debris disc. Due to the different angular dependence of the scattered phase function and polarized phase function on scattering angle, the distribution of the disc SB and therefore the appearance of the disc are different in the images of scattered and scattered polarized light. For a nearly edge-on disc, the polarized disc flux has a maximum at scattering angles close to θ = 90° (e.g. Engler et al. 2017) while the scattered flux peaks at θ = 0° due to diffraction of light. The polarimetric measurement traces the extent of the planetesimal belt because we measure a peak polarized SB of the disc at the location of the belt ansae. Therefore, in the case of the edge-on disc with multiple planetesimal belts, we should measure several SB peaks along the disc axis. These peaks indicate the presence of the inner belts and their extent.

In orderto demonstrate this behaviour, we simulated images of polarized light from an edge-on disc consisting of either one or two belts (Figs. 10a and b). For these synthetic images we used the model described in Engler et al. (2017). In this model, we assume that the degree of linear polarization of the scattered light p_sca has an angular dependence with the scattering angle θ such as thatfor Rayleigh scattering:

where p_m is the maximum fractional polarization achieved at a scattering angle of θ = 90°.

The parameters of the one-belt model correspond to the best-fit model of the J band (see Col. 6 in Table 3). The two-belts model is the non-coplanar model shown in Fig. 7 (see Sect. 5.2 for parameters). Both models were convolved with the PSF.

Comparison of Figs. 10a and b clearly demonstrates that polarimetric measurement helps to identify a multiple belt system. The polarized intensity image of such a system reveals several SB peaks symmetrically located with respect to the star on both disc sides. In this case, the radial SB profile measured along the disc major axis reflects the location and the number of the planetesimal belts (see insets in Figs. 10a and b).

Compared to the disc images in scattered light (see Figs. 7c or 8a), the polarimetric data offer the better opportunity to discover the inner components of the edge-on debris disc. In the total intensity images, the inner components are invisible. They can show up after the data post-processing, for example with the KLIP algorithm, but at very small radial separations from the star; they therefore remain inaccessible (Figs. 7d and 8b).

The polarimetric data of HD 15115 in the J band (Fig. 10d) could be interpreted as a two-belt system. However, for a meaningful conclusion concerning the disc structure based on the modelling results, a higher S/N of the polarimetric data is needed. Therefore we do not perform a fitting of the disc models to the Q_φ image to explore the whole parameter space but compare qualitatively (visually) and quantitatively (-value) only the one-belt model and the non-coplanar two-belts model discussed in Sect. 5.2 with the polarimetric data (Fig. 10d).

To do this, we inserted both models into the Q_φ image which contains no disc flux and was already used for the illustration of the toy model (Sect. 4.1). The visual comparison of the resulting images reveals that the one-belt model (Fig. 10c) lacks the flux in the inner part of the disc in comparison to the data and to the two-belts model (Fig. 10e) and therefore the latter matches the Q_φ image better.

Using the same areas as for the total intensity data and the same binning as well, we also found a marginal preference for the two-belts model with contrary to the one-belt model with .

Fig. 10. Model images of debris disc in polarized light (top row) and comparison of the 4 × 4 binned Qφ image (panel d) with the one-belt model (panel c) and the two-belts model (panel e) in the bottom row. The insets in the top row show the radial SB profiles of the one-belt model (panel a) and the two-belts model (panel b). All models are convolved with the instrumental PSF. Colour bar shows flux in counts per binned pixel.

5.4. Spine diagnostics II

As discussed in Sect. 4.3, the spine measurement can be used as a diagnostic tool for the estimation of the disc radius and inclination. However, one should take into account the degeneracy between these parameters and other geometrical parameters describing the radial distribution of the dust such as α_in , α_out, the HG asymmetry parameter g, and aspect ratio H₀∕R₀ (Engler et al. 2017). The course of the spine curve is defined by a combination of all these parameters. Therefore, in the following we investigate how the shape of the spine curve is affected if one of the parameters α_in , α_out, g, H₀ ∕R₀ or inclination varies. For this analysis, we use the KLIP processed grid models with the fixed radius of the belt (R₀ = 2′′) as described in Sect. 5.

Figure 11 illustrates the degeneracy between α_in , α_out, g, H₀ ∕R₀, and disc inclination. In each panel, one parameter takes three different values while the others remain constant. As can be seen in Fig. 11a, the steeper the dust number density law inside the belt (increasing α_in ), the higher the spine elevation above the midplane inside the parent belt, whereas the steeper fall-off of the grain number density (decreasing α_out) shifts the spine closer to the midplane (Fig. 11b). On the contrary, the spine is displaced upwards when the asymmetry parameter g increases (Fig. 11c). The impact of the disc aspect ratio H₀∕R₀ on the position of the spine is similar to that of the asymmetry parameter – the larger the disc scale height, the larger the spine elevation (Fig. 11d). The vertical offset of the spine from the disc axis can also be affected by the disc inclination. The smaller the inclination, the larger the vertical displacement of the spine and the steeper the curve slope in the inner and outer parts of the disc (Fig. 11e).

Interestingly, we note that the curvature of the spine in the outer region (beyond r = 2′′) changes with decreasing α_out: the spine is bowed for α_out ≤ −6 while it is more like a straight line for larger values. There might be a trend for a bow in the data for r > 3′′ but the S/N of the disc is getting small at these separations. This bow could be partially reproduced by some particular model (Fig. 11b).

Figure 11 clearly demonstrates that the spine is fully determined by a combination of all parameters once the disc radius is fixed. Therefore, as an alternative to calculating a 2D χ² on scattered light images, which implies a large number of degrees of freedom, we used the spine itself and derived for angular separations in the range 1–3′′ on the west side of the disc as already done in Boccaletti et al. (2012) for the HD 32297 debris disc. We define the according to

where N_data is the number of data points between 1′′ and 3′′, y_{i data} are the vertical offsets of the spine measured from the data with uncertainties σ_{i data} (see Sect. 4.3), and y_{i model} are the vertical offsets of the model spine.

The minimal value is obtained for α_in = 4, α_out = −7, g = 0.6, H₀ ∕R₀ = 0.03 and incl = 86°. However, there are other combinations of parameters, in particular the model with α_in = 3, α_out = −7, g = 0.5, H₀ ∕R₀ = 0.01 and incl = 86° with , which provide a similar good fit to the data. As can be seen in Fig. 11f, the spines of these two models show a high degree of similarity with the data. For comparison purposes, the spine of the two-belts model (), where the inner belt is inclined at 80° (see Sect. 5.2), is also presented (blue asterisks in Fig. 11f). Investigating a larger parameter space with the two-belts model should allow to find a better match to the spine curve.

Fig. 11. Investigation of the parameter degeneracies between αin, αout , g, H0 ∕R0 , and inclination. Red dots mark the location of the spine measured in all 4 × 4 binned imaging data as a vertical offset of the Moffat profile peak (see Sect. 4.3 for detailed description). Panel a: impact of varying αin on the spine location of the model with fixed parameters g = 0.4, αout = −4, H0 ∕R0 = 0.03, incl = 86°. Panel b: impact of varying αout on the spine location of the model with fixed parameters g = 0.4, αin = 6, H0 ∕R0 = 0.03, incl = 86°. Panel c: impact of varying g on the spine location of the model with fixed parameters αin = 6, αout = −4, H0 ∕R0 = 0.03, incl = 86°. Panel d: impact of varying H0∕R0 on the spine elevation of the model with fixed parameters αin = 6, αout = −4, g = 0.4, incl = 86°. Panel e: impact of varying inclination on the spine location of the model with fixed parameters αin = 6, αout = −4, g = 0.4, H0 ∕R0 = 0.03. Panel f: comparison of the two best-fitting one-belt models and the two-belts model with the data. The parameters of the model best tracking the spine are g = 0.6, αin = 4, αout = −7, H0 ∕R0 = 0.03, incl = 86°. The second best model has g = 0.5, αin = 3, αout = −7, H0 ∕R0 = 0.01 and incl = 86°.

6. Photometry

We used the instrumental PSFs for the flux normalisation, in both total intensity and polarimetry. In classical imaging, the two IRDIS channels are summed up for photometric measurements. The stellar flux is measured in a circular aperture centred on the PSF with a radius of r = 0.6′′ (r = 50 pixels) after the subtraction of the background (see Sect. 3.1). Taking the transmission of the neutral density filter ND2.0 and DITs (see Sect. 2) into account, we obtain a stellar count rate of (7.27 ± 0.80) ⋅ 10⁷ ADU s⁻¹ in J band and (6.52 ± 0.78) × 10⁷ ADU s⁻¹ in H band. To measure the disc flux from the total intensity data, we first calculated the correction factor to account for the ADI-induced self-subtraction using the best-fit models obtained for both bands (see Col. 6 in Tables 2 and 3). The ratio of the unprocessed and the ADI-processed model images provides the flux reduction factor for each pixel (in the binned images for each binned pixel) inherent to the ADI technique. Hereafter, the disc flux is corrected for the ADI-induced flux reduction. The total disc flux is integrated in two rectangular areas containing pixels with (x, y) coordinates lying within the bounds 1′′ < |x| < 4′′ and − 0.14′′ < y < 0.26′′. We obtain 2900 ± 300 ADU/s and 4100 ± 500 ADU s⁻¹ in J and H, respectively. This yields a brightness ratio between the disc and the star of (4.0 ± 0.6) × 10⁻⁵ in the J band and (6.3 ± 1.0) × 10⁻⁵ in the H band for scattered light.

The polarized flux measured in the same areas in J band is 240 ± 70 ADU s⁻¹. This flux has to be corrected for the reduced instrument throughput in DPI (Langlois et al. 2010, SPHERE User manual). Since the flux measurement in DPI is performed with all instrument components placed in the optical path, we estimate the correction factor by comparing the total stellar flux received in both DPI channels with the stellar flux measured in both CI channels (see Sect. 2). Calibrating the polarimetric measurement in this way, we assume that the stellar light is unpolarised and that the polarimetric efficiency of the instrument did not significantly change during the observation. We obtain a flux reduction factor in DPI of 3.8 for the left channel and 3.5 for the right channel. Taking into account the correction factor of 3.5, the actual total polarized flux of the disc amounts to 840 ± 240 ADU s⁻¹. With this value, the polarized light contrast of the disc in J band is (1.2 ± 0.6) × 10⁻⁵. The ratio between the intensity of the polarized light and total intensity of the scattered light is 0.29 ± 0.09.

In the total intensity we measured that the east side flux is ~ 40% of the west side flux in both J and H band data (more precise ~38 and ~ 44%, respectively) while in polarimetry F_east = 0.8 F_west.

6.1. Disc colour

Total intensity and polarimetric SB profiles are displayed in Fig. 12 for comparison. Here, the disc flux is measured in square 0.25′′× 0.25′′ (20 × 20 pixels) apertures centred on the disc axis. This aperture size was chosen for a better comparison with the radial SB profiles of the HD 15115 disc in V and H bands presented by Kalas et al. (2007; their Fig. 3) and in 1.1 μm band published by Debes et al. (2008; their Fig. 3).

For stellocentric distances ranging between 1.5′′ and 3.0′′ we find that the west side is ~ 0.2–0.5 mag arcsec⁻² brighter than the east side in the J band while this difference is as large as ~ 0.7–1 mag⋅arcsec⁻² in the H band.

The steeper fall-off of the SB in the H band ( ∝ r^−3.3 in the range 1.2–4′′ on the west side and ∝ r^−5.9 in the range 1.2–3′′ on the east side) compared to the J band ( ∝ r^−2.6 in the range 1.2–4′′ on the west side and ∝ r^−3.5 in the range 1.2–3′′ on the east side) produces a change of colour at a particular separation to the star. We observe a red or grey colour (J–H ≥ 0) at short separations, while the colour turns to blue (J–H < 0) at large separations. This colour change occurs between 2.5 and 3.5′′ on the west side and between 2 and 3′′ on the east side. For large separations, J–H can be as large as − 0.7 mag arcsec⁻².

It is interesting to note that the polarized light SB profiles of the west and east side show a high degree of similarity: a slow growth between 1.2′′ and 2′′ and a fall-off ∝ r^−3.4 for r > 2′′. The SB profiles are truncated because of the low S/N of the polarimetric data at large separations.

6.2. Polarization fraction

The polarization fraction (PF) of scattered light in the J band is measured in square 0.25′′× 0.25′′ apertures along the disc axis (Fig. 13) and is calculated from the ratio between the SB of the polarized and the total intensity data (corrected for the ADI bias). The radial PF profiles of both disc sides show a similar behaviour, increasing from PF ≈ 0.03 at r = 1.25′′ up to ~ 15–20% between r = 2′′ and r = 3′′, thus close to the edge of the planetesimal belt. This value is consistent with the degree of polarization in the K band measured by Tamura et al. (2006) along the midplane of the β Pictoris debris disc. A similar fractional polarization (~10–25%) is obtained by Asensio-Torres et al. (2016) for another edge-on disc HD 32297 in the H band.

At radial separations of r ≳ 3′′, the S/N of our polarimetric data is too low to make a conclusion about the polarization degree of scattered light. It is also not possible to determine whether the PF further increases or decreases, or remains approximately constant with a growing radial distance as was measured for an edge-on debris disc around M star AU Mic (Graham et al. 2007).

Fig. 12. Radial SB profiles of the west and east disc sides in the J and H band. All data points show the disc magnitude measured in square 0.25′′× 0.25′′ (20 × 20 pixels) apertures, centred on the disc major axis, relative to the stellar magnitude in the corresponding filter. Circles indicate the relative magnitude of the total intensity of scattered light (Stokes I) and diamonds mark the magnitude of the polarized flux (Stokes Qφ). Errors are calculated as a standard deviation of the flux distribution in the concentric annuli.

Fig. 13. Polarization fraction of the scattered light measured in square 0.25′′× 0.25′′ apertures along the disc axis using total intensity (Stokes I) and polarized flux (Stokes Qφ) profiles in the J band plotted in Fig. 12.

6.3. Surface brightness profiles

The aperture size 0.25′′× 0.25′′ (20 × 20 pixels) used in the previous sections is too large to see SB variations on smaller scales while the angular resolution provided by IRDIS allows for a more detailed examination of the flux distribution inside r = 2′′. In Fig. 14, we used a finer sampling of 4 × 4 pixels (0.049′′× 0.049′′) together with a σ = 1 pixel and σ = 2 pixels Gaussian smoothing for classical imaging and polarimetric data, respectively. The errors were computed as the standard deviation of the flux distribution in concentric annuli centred on the star and do not include systematic uncertainty.

The SB profiles measured in all data sets show a peak at a radial separation of r ≈ 2′′ from the star on both disc sides (denoted as “F1” in Fig. 14c). This peak is clearly attributed to the outer planetesimal ring. All three observations show remarkably consistent results for the radius of this belt: for the east side we measure R_east = 2.00′′ ± 0.05′′ and for the west side R_west = 1.91′′ ± 0.05′′. The small difference in the belt radius between both sides indicates an offset of ~ 0.045′′ of the belt towards the east with respect to the star.

Further inside the ring, the SB profiles exhibit another peak at r ≈ 1.2−1.3′′ also on both disc sides in classical imaging data of J and H bands (denoted as “F2” in Fig. 14c). In polarized intensity data, a SB peak is detected only on the east side at the same radial separation. The location of this peak corresponds to the increased SB at the position indicated by the arrow “B” in Fig. 2. We observe some kind of correlation between these peaks and think that they could mark the ansae of the inner belt. There is a third SB peak on the east side of the disc (“F3” in Fig. 14c) but at a lower significance level.

The increased SB at the locations “F1”, “F2”, and “F3” is also visible in the disc images in Fig. 15. For an easier comparison of Figs. 14 and 15 and better identification of the SB peak locations, the images in Fig. 15 are presented with the same x-axis scale for 1 arcsec as in Fig. 14. The vertical dotted lines going through the positions of the SB peaks in Fig. 14 are also shown in Fig. 15.

Fig. 14. Radial SB profiles measured along the disc major axis in 0.049′′× 0.049′′ square apertures in the cADI image of the J band data (panel a), the cADI image of the H band data (panel b), and in the Qφ image (panelc, J band). All data sets show the SB peaks at r = 2′′ (“F1”) and r ≈ 1.3′′ (“F2”) symmetrically located with respect to the star, except the west disc side in the Qφ image. Images of the disc in Fig. 15 are presented with the same x-axis scale for 1 arcsec as in this figure for ease in identification of the SB peak locations highlighted with the vertical dashed lines.

Fig. 15. Images of the HD 15115 disc illustrating the position of SB peaks plotted in Fig. 14. The original data were 4 × 4 binned to reduce the noise and smoothed via a Gaussian kernel with σ = 2 pixels. The position of the star in each image is marked by a yellow dotted line and corresponds to the position of the star in Fig. 14 (yellow line). Images are presented in linear scale: J band in [−25, 85], H band in [−120, 400], J band polarimetry in [−10, 34].

6.4. Phase function

In the total intensity data, in particular in the H band dataset, the back side, although faint, is detected. The observed part of the planetesimal ring shown in the inset in Fig. 16 spans the range of the scattering angles between 35° and 140°. Measuring the azimuthal variation of the SB along the rim of the ring, we can constrain the scattering phase function. For the calculation of the scattering angles, we need to make some general assumptions on the grain properties and geometry of the disc. We assume that grains scatter stellar radiation preferentially in the forward direction; which would mean that the northern part of the disc is the near side. We also consider that the circumstellar material is confined to a narrow ring with a radius of R₀ = 1.95′′ (94 au). Using the cADI total intensity image in the H band obtained with selected frames (only frames with a seeing condition below 0.9′′ were used), we performed the aperture photometry of the west disc side in circular apertures with a radius of 6 pixels. The disc area within each aperture was attributed a unique scattering angle which was calculated for the near side of the disc according to: (10)

and for the back side of the disc (11)

where (x_c, y_c) are coordinates of the aperture centre in the coordinate system shown in Fig. 3a. The apertures were centred on the spine of the front and back side of the disc defined as the brightest points of the vertical cross-sections perpendicular to the disc major axis.

To correct for the flux loss caused by the ADI data reduction, the measured flux values were multiplied with a scaling factor derived from the modelling of ADI bias (see Sect. 6). The errors of this measurement are estimated as one standard deviation of the flux distribution measured in the apertures of the same size which were placed at the same radial distance from the star but at different PAs.

The obtained intensity profile is shown in Fig. 16. The flux values are arbitrarily normalized so that the intensity at 83° equals 0.15. As Fig. 16 shows, the scattering intensity steeply increases with a decreasing scattering angle for θ < 50° and gives rise to a relatively broad diffraction peak. The measured curve is rather flat between 60° and 90°, whereas for larger angles corresponding to the fainter back side of the disc (θ > 90°) the scattering intensity is decreasing. For comparison, we measured the intensity profiles of several models performing the same aperture photometry. Instead of plotting each model profile individually, we show the range of possible model values with a shaded red area in Fig. 16.

We usedthe χ²-minimization algorithm to find a HG function (Eq. (6)) that matches the intensity profile and fitted the data also with a linear combination of two HG functions following the same approach as in Milli et al. (2017). According to the -criterion, the better fit f(θ) is given by a linear combination of two HG functions (Engler et al. 2017): (12)

with g₁ = 0.94 ± 0.06, g₂ = 0.21 ± 0.01 and scaling factor w = 0.89 ± 0.10 (blue curve in Fig. 16). For this combination, the reduced metric is 1.22. This is smaller than of 1.43 which is obtained for the best-fit of the HG function with the single asymmetry parameter g = 0.57 ± 0.0 (purple curve). However, the visual comparison of both fits in Fig. 16 shows that they are nearly identical and thus two additional free parameters of the combined HG functions do not significantlyimprove the fit. The lower BIC (Sect. 5.2) for the HG function with the single g parameter also indicates that this model should be considered as the best-fitting one.

In Fig. 16, we also compare the measured intensity profile with three scattering phase functions derived by Min et al. (2010) for dust grains with radii equal to 5, 15, and 50 μm which are covered by small regolith particles³. The grains are treated as micro-asteroids reflecting light in the geometric optics regime. Their phase functions are obtained for a narrow grain size distribution and include the diffraction by spherical particles and the reflectance function according to the analytical model for the bidirectional reflectance (Hapke 1981). For a comparison, all three functions in Fig. 16 are normalized to the data point at θ = 83°.

At large scattering angles (θ > 90°) the measured phase function of the HD 15115 disc does not follow the phase functions obtained by Min et al. (2010). When comparing the diffraction parts of the phase functions (θ < 80°), the scattering function for grains with a size of 5 μm (solid black curve) fits best the measured intensity profile. Assuming that this grain size is the average size ⟨a⟩ of the grain size distribution, we derive the minimum grain size a_min = 3 μm from the relation a_min = 3∕5 ⟨a⟩ which is valid for the dust grains in a collisional cascade equilibrium expected at the location of the parent body ring. This minimum grain size is in agreement with the estimated blowout size for HD 15115 of ~ 1–3 μm assuming a stellar luminosity of 3.5 L_⊙, a mass of 1.3 M_⊙, and a grain average density in the range from 1 to 3 g cm⁻³ (Rodigas et al. 2012; Hahn 2010).

It is interesting to note that the comparison of micro-asteroids scattering functions with the observed phase functions of the debris discs around HR 4796 A and Fomalhaut shows the compliance with the curves for grains with radii of a = 15 μm (Milli et al. 2017) and a = 50 μm (Min et al. 2010), respectively. This could suggest that the grains in the HR 4796 A and Fomalhaut discs are larger on average than the dust grains we observe on the west side of the HD 15115 debris disc.

Fig. 16. Normalized scattering intensity (red diamonds) measured in circular apertures along the disc spine on the west side shown in the inset image. The blue and purple curves show HG fits to the data points. The black curves show scattering phase functions obtained by Min et al. (2010) for three different sizes of grains a covered withregolith particles. The shaded grey area indicates 1σ-uncertainties of the measurement. The shaded red area shows the range of values measured from the models.

7. Discussion

7.1. Unpolarized versus polarized light

As mentioned in Sect. 6, the polarized intensity (Stokes Q_ϕ ) measured for the east side of the disc is a factor of 1.25 lower compared to the west side. In the total intensity data (Stokes I) of the J band, we obtain for the east side an intensity which is a factor of 2.67 lower than on the west side. The different factors for the total and polarized intensities, when comparing the east side with the west side, can indicate a difference in the grain number density, grain properties, or both.

In the case of HD 15115, this could be the result of an asymmetric loss of small particles due to the eastward motion of the system through the interstellar medium as suggested by Debes et al. (2009). Such motion could produce a surplus of smaller grains alongthe line of sight on the west side in front of the planetesimal belt. With polarimetry, we trace predominantly 90° scattering and therefore the ansae of the belt. The total intensity strongly favours forward scattering from grains in front of the star and the parent belt. This could cause the strong Stokes I signal in the west where the particles produce a lot of “extra” scattered intensity I but only a little “extra” Stokes Q_ϕ signal.

An alternative interpretation of the disc brightness asymmetry was recently proposed by Lee & Chiang (2016). They presented a model explaining the “needle”-like morphologies observed in the edge-on discs by a planet perturbation. Their model requires an eccentric parent belt with a brighter disc side much closer to the star and much shorter in length than the fainter side. Such a model is in contradiction with our data. We observe a symmetric parent belt (see Fig. 14) in both J and H bands, in total and in polarized intensity. If the brightness asymmetry of the HD 15115 disc is caused by planets, then the perturbation mechanism should be more complex than the model suggested by Lee & Chiang (2016).

7.2. Colour

The photometric analysis (Figs. 12 and 14) suggests that inside the planetesimal belt (r ≲ 2′′) the disc colour is red and becomes bluer with increasing distance from the star (r > 2′′). A similar colour behaviour has previously been observed for the HD 15115 disc by Debes et al. (2008) and Kalas et al. (2007). Debes et al. (2008) argued that the change in the dust scattering efficiency beyond r ≈ 2′′ could be due to the absorption features of dust material at 1.6 μm and inside r ≈ 2′′ due to absorption of olivine at ~1 μm, for instance. Another possibility to explain the complex colour of the disc would be the existence of dust populations with different average grain sizes varying with distance from the star. At the radial position of the “birth” ring we would expect to find dust grains larger on average than the blowout size a_blow. Beyond the “birth” ring small dust particles might dominate the scattering cross section due to the radiation blow-out of the smallest dust grains (Stark et al. 2014; Krivov 2010, and references therein). Assuming that the collisional cascade power law for the grain size distribution dn(a) ∝ a^−3.5 da is valid at the location of the parent body ring, we would expect that this distribution becomes steeper, and therefore the disc colour bluer, with distance from the star.

7.3. Evidence for multiple belts?

Our images show an axisymmetrical disc in agreement with the observations of Kalas et al. (2007) and Rodigas et al. (2012). Mazoyer et al. (2014) have shown that in H and K_S filters, the disc has a symmetric ring shape with a radius of R₀ ≈ 2′′ and an inclination of ~ 86° well matching the geometrical parameters which we derive for the ring observed in the broad bands H and J (Fig. 1). This ring could be an outer planetesimal belt of a multiple ring system.

The disc morphology seen in our polarimetric data (Sect. 4.1, Fig. 3) suggests the presence of a second inner planetesimal ring which could have a slightly different PA and inclination compared to the outer belt. The structural clumps that we see in the Q_φ image (Fig. 2b) and ripples of the SB profile along the disc axis, in particular on the east side of the disc (Fig. 14) are similar to the ripples in the SB profile of the polarized light detected by Tamura et al. (2006) in the β Pic debris disc. The 2 μm polarimetry of the β Pic disc has revealed many details in the inner part of the disc, including the previously observed warp. In analysing the disc structure in the region between 50 and 121 au, Tamura et al. (2006) concluded that the ripples of the polarimetric SB profile may indicate a multiple ring system. We think that we see the same effect in the SB profile of the polarized light in Fig. 14c. It is obvious that the outer SB peaks (denoted as “F1” in Fig. 14c) mark the location of the outer ring (R ≈ 2′′) of the HD 15115 disc. Therefore the other SB peaks, which we measure further inside the disc, may also indicate inner rings.

The modelling of the total intensity and polarized intensity data in both J and H bands does not rule out the two-belts model. According to the χ² criterion forthe goodness-of-fit, which we have calculated for more than 2000 models, the two-belts model provides a good fit to the total intensity data, being always among the 20 best models. For the H band dataset, it is even the best model when fitting both disc sides simultaneously (see Table 2).

In Appendix C, we show that the two-belts model is also consistent with the spectroscopic and photometricdata of HD 15115, if the modified blackbody function (Backman & Paresce 1993) is used to fit the disc SED.

There are some indications from previous studies that support the idea of HD 15115 as a multiple belt system:

• As already mentioned in Sect. 1, Moór et al. (2011) and Chen et al. (2014) claimed that the disc SED canbe better fitted with a combination of two blackbody models, which implies the presence of two planetesimal belts. However, both authors found that an inner belt has a small radius of 4 or 5 au. Therefore this belt cannot be the same inner belt suggested by us.

• Debes et al. (2008) and Rodigas et al. (2012) measured an increasing PA of the eastern disc side with decreasing distance to the star for r < 1.6′′ (see their Figs. 2a and 4a, respectively) which is in contradiction to a bow-like disc shape reported by Rodigas et al. (2012) and Schneider et al. (2014). Both authors suggested a nearly edge-on disc model with a high HG scattering asymmetry parameter to explain the “bowing” of the brighter (northern) disc side. In accordance with this model, the PA of the east disc side should decrease following the course of the disc spine when approaching the star. This means that the vertical offset of the spine from the disc axis towards the north should increase for smaller separations, as is expected for a single belt with a high inclination and a high asymmetry parameter for forward scattering (see Fig. 12 in Schneider et al. 2014 for an illustration). The measured vertical offset of the spine is decreasing in observations at 1.1 μm (Debes et al. 2008, their Fig. 2a) and in K_S band for r < 1.6′′ (Rodigas et al. 2012, their Fig. 4a). These measurements could be affected by the lower spatial resolution of the PISCES and NICMOS data and by their very special structure of the stellar halo with more low spatial frequencies, which could mimic a disc structure. However, our data also show a small decrease of the spine offset between r = 1.5′′ and r ≈ 1.25′′ on the east side in the total intensity images of both J and H bands (see Fig. 5).

Furthermore, Rodigas et al. (2012) measured a clear break in the radial SB power law at r ≈ 1.4′′ on the east side of the disc: from a power law index of −1.1 for r > 1.4′′ to a power law index of −5.27 for 1.2′′ < r < 1.4′′ (see their Fig. 5a). The power law break of the radial SB profile is often connected with the location of the planetesimal belt. The rapid increase of the SB could be due to stellar residuals or due to an additional source of emission inside r ≈ 1.5′′ which could be a planetesimal ring. In Fig. 17, we reproduce Fig. 1a from Rodigas et al. (2012) presenting a K_S band image of the debris disc; it shows the increasing SB below the disc axis close to the coronagraph only on the east side (indicated with a white arrow in Fig. 17). We see a similar feature in our H band data at the same position (see Fig. 18). It might be responsible for the small dip in the spine curve (see Fig. 5) mentioned above.

• The inner belt with a slightly different PA from the outer belt can easily explain the bifurcation of the east disc side observed by Schneider et al. (2014) and Debes et al. (2008, see their Fig. 1; feature A is explained as a contamination from the PSF subtraction but it roughly coincides with the bifurcated structure). The radial blowout of the smallest dust grains from two planetesimal belts at different PAs can create the same morphology.

It is also possible that the SB peak, which we observe to the east at r ~1.3′′ from the star (indicated with a white arrow in Fig. 2b and with the arrow “B” in Fig. 3b) and attribute to the east edge of the inner belt, is actually associated with some clumpy feature of the outer planetesimal belt. Assuming that this dust clump moves on the Keplerian orbit, we could observe its orbital movement when its deprojected radial position changes by ~ 0.2′′ (clump width). Such a position shift on a circular orbit of the radius r = 2′′ for a 1.3 M_⊙ star would be observable in about 18.5 yr.

Fig. 17. KS band image of the HD 15115 debris disc as published by Rodigas et al. (2012).

Fig. 18. H band image of the HD 15115 debris disc (cADI, this paper).

8. Summary

In this paper, we present classical imaging and polarimetric observations of the HD 15115 debris disc known as the “blue needle” with SPHERE-IRDIS in the broadbands J and H. The total intensity images (Stokes I, Fig. 1) in both bands were extracted with ADI using cADI, KLIP, and LOCI processing algorithms. Images of the polarized scattered light (Stokes Q_φ , Fig. 2) in the broad band J were obtained using the double difference method.

Our results from the analysis of the disc geometry and photometry can be summarised as follows:

• The overall disc structure observed in both bands is an axisymmetric, nearly edge-on debris ring consistent with the detection by Mazoyer et al. (2014). The radius of the disc towards the east is R_east = 2.00′′ ± 0.05′′ (~96 au) and towards the west it is R_west = 1.91′′ ± 0.05′′ (~92 au). By fitting an ellipse to the spine curve we derive a disc inclination of 85.8° ± 0.7°.

• The PA of the disc major axis measured in the total intensity and polarized intensity data is θ_disc = 278.9° ± 0.1°.

• The disc photometry and radial SB profiles along the disc major axis show the east-west asymmetry. In two rectangular areas bounded by the coordinates 1′′ < |x| < 4′′ and − 0.14′′ < y < 0.26′′ in the total intensity data of J and H bands, we measure ~2.5 times less flux on the east side compared to the west side. In the same areas of the polarimetric data taken in the J band, this difference is smaller and the polarized flux on the east side amounts to ~80% of the flux measured on the west side (Sect. 6). This might indicate a different average grain size in the disc to the east compared to the disc to the west due to the loss of small dust grains towards the west.

• Based on the forward modelling results, we measure a change in the disc colour between 2′′ and 3.5′′ from red (for r < 2′′) to predominantly blue (for r > 3.5′′). This colour difference may indicate different grain size distributions at the location of the parent body belt and in the outskirts where the smallest grains are swept out of the disc.

• The image of the polarized light from the disc exhibits some SB peaks which can also be seen in the radial SB profile along the disc axis (Fig. 14c). Similar SB peaks can be recovered in the total intensity data (Figs. 14a and b). They may be indicative for the presence of a second (inner) planetesimal ring. Based on the J band polarimetry, we estimate that such a ring would have a radius of between 1.3′′ and 1.5′′ and a slightly different PA from the outer ring.

• We have performed extensive model calculations to investigate the compliance of the two-belts model with the total intensity and polarized intensity data. According to the χ² criteria, both the one-belt and the two-belts models provide similar goodness of the fit. Therefore, the existence of the inner ring cannot be rule out. Taking into account the relatively low S/N of our data, no robust conclusion can be made.

To prove the presence of the second belt in the inner regions, it would be valuable to re-observe the HD 15115 disc with ZIMPOL/SPHERE polarimetry. An observation at shorter wavelengths would provide a required high resolution of ~ 25–30 mas when observing in Very Broad Band (VBB, λ_c = 735 nm) or I band (λ_c = 790 nm). The Strehl ratio is higher in the I band but the count rate is ~40% lower compared to the VBB. Under good observing conditions (seeing ≤ 0.7, coherence time ≥ 3.5 ms) a high contrast of 10⁻⁶ can be achieved. Polarimetry as a differential technique has also a decisive advantage over ADI in that it does not suffer from the disc self-subtraction, which could be crucial if we want to study the close-in disc structure.

Acknowledgements

We would like to thank the referee for many thoughtful comments which helped to improve this paper. This work issupported by the Swiss National Science Foundation through the grant number 200020 - 162630. We acknowledge the financial support from the Programme National de Planétologie (PNP) of CNRS-INSU co-funded by the CNES, and by the Agence Nationale de la Recherche (ANR-14-CE33-0018). J.O. acknowledges financial support from the ICM (Iniciativa Científica Milenio) via the Núcleo Milenio de Formación Planetaria grant, from the Universidad de Valparaíso, and from Fondecyt (grant 1180395). J.M. acknowledges support from the NASA through the NASA Hubble Fellowship grant HST-HF2-51414.001 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555.

Appendix A: Signal-to-noise ratio

Figure A.1 demonstrates the S/N of the total intensity images in J and H bands.

Fig. A.1. S/N map for different data-reduction techniques in the H band (left column) and in the J band (right column). The number of modes in KLIP reductions is 10. The star position is marked with a white cross.

Appendix B: Position angle of the disc

Fig. B.1. Qφ image illustrating the method to determine the PA and size of the outer disc. The image was 4 × 4 binned and smoothed via a Gaussian kernel with σ = 1 pixel. The brightest pixels of the regions inside the yellow circles indicate the location of the planetesimal belt which is schematically drawn as a yellow ellipse. For the clarity, the ellipse was shifted upwards (white dotted lines show the parallel shift of the brightest pixels). The position of star is marked by a yellow asterisk. The white line represents the disc major axis. The colour bar shows the counts per binned pixel.

Appendix C: Consistency of the two-belts configuration with the disc SED

To check the consistency of the two-belts model with the available spectroscopic and photometric data of the HD 15115 disc, we fitted the disc SED with a single temperature blackbody (BB) emission B_ν (λ, T_dust) and with the modified Planck function (MBB, Backman & Paresce 1993):

where B_ν(λ, T_dust) is the Planck function defined by the wavelength λ and the dusttemperature T_dust, H is the Heaviside step function, λ₀ is the characteristical wavelength, above which the absorption efficiency of dust grains scales like with β being the opacity index. The free parameters λ₀ and β allow to account for the steeper falloff of the debris disc spectra at longer wavelengths than the Planck function and to have a higher dust temperature than the blackbody temperature. The free parameter of the MBB model are the dust temperature T_dust, the opacity index β, the characteristic wavelength λ₀ and the scaling factor. In the BB model, the belt radius and thus the dust temperature are fixed. The only free parameter of this model is the scaling factor for the Planck function. We adopted a belt radius of 2′′ (96 au) with T_dust = 39 K for the one belt model, and additionally a radius of the second belt of 1.3′′ (61 au) with T_{dust inn} = 48 K for the two-belts model. To fit the thermal emission from two planetesimal belts, we used a linear combination of two BB or two MBB functions, assuming the same optical properties and average grain size in both belts.

The disc SED is constructed using the sampled 5–35 μm Spitzer/IRS spectrum (Moór et al. 2011) and fluxes compiled from the literature (Table C.1). The stellar photosphere emission of a star is calculated with L_* = 3.5 L_⊙ and T_eff = 6780 K (Moór et al. 2011).

The derived best fits for each model are plotted in Fig. C.1. The single and double temperature BB models (blue lines) cannot reproduce the data well. The contribution of the lower temperature component in the two-belts BB model to the derived best fit (blue solid line) is negligible. Contrary to that, both MBB models, with one and two belts (red lines), respectively, fit the data very well, except for one data point at 13 μm (Chen et al. 2014) which also disagrees with the Spitzer/IRS spectrum (AOR 10885632) retrieved from the archive.

Fig. C.1. Disc SED and the best fits for the BB and MBB models with one and two belts. Red asterisks show the sampled 5–35 μm Spitzer/IRSspectrum. Violet asterisks represent the photometric data specified in Table C.1.

References

All Tables

All Figures

	Fig. 1. Total intensity images of the HD 15115 debris disc obtained with different data-reduction techniques in the H band (left column) and in the J band (right column). The number of modes in KLIP reductions is 10. The star position is marked with a yellow asterisk. The FOV of each displayed image is 8.2′′ × 3.1′′.
In the text

Fig. 2. Polarized intensity images of the HD 15115 debris disc in J band. Images show the Stokes Qφ parameter smoothed by a Gaussian kernel with σ = 2 pixels (panela) and with an additional 4 × 4 binning (panel b). The white arrow points to the local surface brightness peak responsible for the impression of the inner belt (see also Fig. 3). The star position is marked with a yellow asterisk. The colour scale shows the flux in counts per binned pixel. The surface brightness contours in panel b are given for 5 and 15 counts per binned pixel. The FOV of each displayed image is 8.2′′ × 3.1′′.

In the text

Fig. 3. Qφ images binned by 10 pixels along the x-axis and by 2 pixels along the y-axis and smoothed via a Gaussian kernel with σ = 1 pixel (panel a) and with σ = 2 pixels (panel b) to highlight the structure and position of both rings relative to each other. The central part of the images, dominated by strong residuals, is masked out. The position of the star is indicated by a yellow asterisk. The x-axis (interrupted white line going through the star position) shows the outer disc major axis at PA = 279.0°. The FOV of the displayed panel a is 9.8′′× 0.74′′. The letters “A”, “B”, “C” and “D” are used to label the surface brightness peaks discussed in Sect. 4.1. Colour bar shows flux in counts per binned pixel.

In the text

	Fig. 4. Comparison of the double (column a) and single (column c) belt models with the Qφ image (column b). All images are binned by 10 pixels along the x-axis and by 2 pixels along the y-axis and additionally smoothed via a Gaussian kernel with σ = 2 pixels (top row) and with σ = 1 pixel (bottom row).
In the text

	Fig. 5. Spine vertical distance from the disc major axis measured as a vertical offset of the Moffat profile centre. The plot traces the mean location of the spine in all available images (cADI, KLIP, LOCI reductions of the total intensity data) in H and J bands. The star position is at (0, 0). The black line shows the best-fit ellipse to the data points between |x| ≥ 1.2′′ and |x| ≤ 2.0′′.
In the text

	Fig. 6. Example of synthetic images of one-belt model. Panel a: model before convolution with the instrumental PSF. Panel b: model convolved with the instrumental PSF. Panel c: model after processing with the KLIP algorithm. The flux is shown in arbitrary units.
In the text

Fig. 7. Synthetic images of the two-belts model with an inner belt (shown separately in panel a) which is non-coplanar with the outer belt. The outer belt has the same geometry and asymmetry parameter g as the one-belt model demonstrated in Fig. 6. Panels b and c: two-belts model before and after convolution with the instrumental PSF, respectively. Panel d: KLIP-processed model. The flux is shown in arbitrary units.

In the text

	Fig. 8. Synthetic images of the two-belts model with an inner belt coplanar with the outer belt. The model parameters are the same as in the model shown in Fig. 7 except for the inclination and PA of the inner belt which are the same as for the outer belt. Panel a: model before convolution with the PSF. Panel b: model after processing with the KLIP algorithm. The flux is shown in arbitrary units.
In the text

Fig. 9. Comparison of the HD 15115 disc (panel a) with the models inserted in the image (H band, KLIP reduction) where the disc wasremoved. Panel b: model with one belt. Panel c: two-belts model with an inner belt which has a different PA and inclination with respect to the outer belt. Panel d: two-belts model with an inner ring coplanar with the outer ring. In the processed images of the one-belt model (panel b) and non-coplanar two-belts model (panel c) the back side of the disc is visible, whereas it is not recognisable in the image of the coplanar two-belts model (panel d). Colour bar shows flux in arbitrary units.

In the text

Fig. 10. Model images of debris disc in polarized light (top row) and comparison of the 4 × 4 binned Qφ image (panel d) with the one-belt model (panel c) and the two-belts model (panel e) in the bottom row. The insets in the top row show the radial SB profiles of the one-belt model (panel a) and the two-belts model (panel b). All models are convolved with the instrumental PSF. Colour bar shows flux in counts per binned pixel.

In the text

Fig. 11. Investigation of the parameter degeneracies between αin, αout , g, H0 ∕R0 , and inclination. Red dots mark the location of the spine measured in all 4 × 4 binned imaging data as a vertical offset of the Moffat profile peak (see Sect. 4.3 for detailed description). Panel a: impact of varying αin on the spine location of the model with fixed parameters g = 0.4, αout = −4, H0 ∕R0 = 0.03, incl = 86°. Panel b: impact of varying αout on the spine location of the model with fixed parameters g = 0.4, αin = 6, H0 ∕R0 = 0.03, incl = 86°. Panel c: impact of varying g on the spine location of the model with fixed parameters αin = 6, αout = −4, H0 ∕R0 = 0.03, incl = 86°. Panel d: impact of varying H0∕R0 on the spine elevation of the model with fixed parameters αin = 6, αout = −4, g = 0.4, incl = 86°. Panel e: impact of varying inclination on the spine location of the model with fixed parameters αin = 6, αout = −4, g = 0.4, H0 ∕R0 = 0.03. Panel f: comparison of the two best-fitting one-belt models and the two-belts model with the data. The parameters of the model best tracking the spine are g = 0.6, αin = 4, αout = −7, H0 ∕R0 = 0.03, incl = 86°. The second best model has g = 0.5, αin = 3, αout = −7, H0 ∕R0 = 0.01 and incl = 86°.

In the text

Fig. 12. Radial SB profiles of the west and east disc sides in the J and H band. All data points show the disc magnitude measured in square 0.25′′× 0.25′′ (20 × 20 pixels) apertures, centred on the disc major axis, relative to the stellar magnitude in the corresponding filter. Circles indicate the relative magnitude of the total intensity of scattered light (Stokes I) and diamonds mark the magnitude of the polarized flux (Stokes Qφ). Errors are calculated as a standard deviation of the flux distribution in the concentric annuli.

In the text

	Fig. 13. Polarization fraction of the scattered light measured in square 0.25′′× 0.25′′ apertures along the disc axis using total intensity (Stokes I) and polarized flux (Stokes Qφ) profiles in the J band plotted in Fig. 12.
In the text

Fig. 14. Radial SB profiles measured along the disc major axis in 0.049′′× 0.049′′ square apertures in the cADI image of the J band data (panel a), the cADI image of the H band data (panel b), and in the Qφ image (panelc, J band). All data sets show the SB peaks at r = 2′′ (“F1”) and r ≈ 1.3′′ (“F2”) symmetrically located with respect to the star, except the west disc side in the Qφ image. Images of the disc in Fig. 15 are presented with the same x-axis scale for 1 arcsec as in this figure for ease in identification of the SB peak locations highlighted with the vertical dashed lines.

In the text

Fig. 15. Images of the HD 15115 disc illustrating the position of SB peaks plotted in Fig. 14. The original data were 4 × 4 binned to reduce the noise and smoothed via a Gaussian kernel with σ = 2 pixels. The position of the star in each image is marked by a yellow dotted line and corresponds to the position of the star in Fig. 14 (yellow line). Images are presented in linear scale: J band in [−25, 85], H band in [−120, 400], J band polarimetry in [−10, 34].

In the text

Fig. 16. Normalized scattering intensity (red diamonds) measured in circular apertures along the disc spine on the west side shown in the inset image. The blue and purple curves show HG fits to the data points. The black curves show scattering phase functions obtained by Min et al. (2010) for three different sizes of grains a covered withregolith particles. The shaded grey area indicates 1σ-uncertainties of the measurement. The shaded red area shows the range of values measured from the models.

In the text

	Fig. 17. KS band image of the HD 15115 debris disc as published by Rodigas et al. (2012).
In the text

	Fig. 18. H band image of the HD 15115 debris disc (cADI, this paper).
In the text

	Fig. A.1. S/N map for different data-reduction techniques in the H band (left column) and in the J band (right column). The number of modes in KLIP reductions is 10. The star position is marked with a white cross.
In the text

Fig. B.1. Qφ image illustrating the method to determine the PA and size of the outer disc. The image was 4 × 4 binned and smoothed via a Gaussian kernel with σ = 1 pixel. The brightest pixels of the regions inside the yellow circles indicate the location of the planetesimal belt which is schematically drawn as a yellow ellipse. For the clarity, the ellipse was shifted upwards (white dotted lines show the parallel shift of the brightest pixels). The position of star is marked by a yellow asterisk. The white line represents the disc major axis. The colour bar shows the counts per binned pixel.

In the text

	Fig. C.1. Disc SED and the best fits for the BB and MBB models with one and two belts. Red asterisks show the sampled 5–35 μm Spitzer/IRSspectrum. Violet asterisks represent the photometric data specified in Table C.1.
In the text