Investigating the presence of two belts in the HD 15115 system

TL;DR

This study uses near-infrared observations and modeling to analyze the debris disk around HD 15115, investigating the potential existence of a second inner belt and characterizing the disk's structure, brightness, and color variations.

Contribution

First detailed near-IR imaging and modeling of HD 15115's debris disk, exploring evidence for a second inner belt and analyzing disk properties in scattered and polarized light.

Findings

Detected an axisymmetric belt at ~2 arcseconds with specific inclination and position angle.

Observed brightness asymmetry with the west side brighter than the east in total intensity.

Identified a color gradient from red to blue with increasing distance, indicating grain size variation.

Abstract

We present new observations of the edge-on debris disk around HD 15115 (F star at 48.2 pc) obtained in the near-IR. We search for observational evidence for a second inner planetesimal ring in the system. We obtained total intensity and polarimetric data in the broad bands J and H and processed the data with differential imaging techniques achieving an angular resolution of about 40 mas. We observe an axisymmetric planetesimal belt with a radius of $\sim$2$''$, an inclination of $85.8^{\circ} \pm 0.7^{\circ}$ and position angle of 278.9$^{\circ} \pm 0.1^{\circ}$. A grid of models describing the spatial distribution of the grains in the disk is generated to constrain the geometric parameters of the disk and to explore the presence of a second belt. We perform a photometric analysis of the data and compare disk brightness in two bands in scattered and in polarized light. The analysis shows that the west side is $\sim$2.5 times brighter in total intensity than the east side in both bands, while for polarized light in the J band this ratio is only 1.25. The maximum polarization fraction is 15--20% at $r\sim$2.5$''$. We also find that the J - H color of the disk appears to be red for the radial separations $r\lesssim2''$ and is getting bluer for the larger separations. This apparent change of disk color from red to blue with an increasing radial separation could be explained by the decreasing average grain size with distance. The presence of an additional inner belt slightly inclined with respect to the main planetesimal belt is suspected from the polarized intensity image but the analysis and modeling presented here cannot establish a firm conclusion due to the faintness of the disk and its high inclination.

Tables

Table 1. Table 1: Summary of observations.

Observation	Instrument	Filter	Field	Integration Time			Observing Conditions
date	mode		rotation [∘]	DIT [s]	Total [min]		Airmass	Seeing [”]	Coh. time [ms]
2015-10-29	IRDIS CI	BB_H	23.2	16	51.2		$1.18\pm 0.02$	$0.92\pm 0.09$	$1.7\pm 0.2$
2015-11-19	IRDIS CI	BB_J	23.5	8	51.2		$1.18\pm 0.02$	$0.78\pm 0.10$	$1.6\pm 0.3$
2016-10-03	IRDIS DPI	BB_J	stabilized	32	85.3		$1.20\pm 0.03$	$0.87\pm 0.18$	$2.2\pm 0.4$

Table 2. Table 2: Best-fit model parameters of the H band data.

Optimized	East side ($\nu =425$)		West side ($\nu =635$)		Both sides ($\nu =1065$)
parameter	Best fit	Models with	Best fit	Models with	Best fit1	Models with
	${\chi }_{\nu ,\min }^2=1.40$		${\chi }_{\nu ,\min }^2=1.53$		${\chi }_{\nu ,\min }^2=3.63$
HG parameter $g$	$0.5$	0.5, 0.6	0.3	0.3	0.3	0.3-0.4
Inner radial index ${\alpha }_{\mathrm{in}}$	9	6-9	2	2-4	2	2-3
Outer radial index ${\alpha }_{\mathrm{out}}$	-6	-6	-4	-4	-4	-4
Aspect ratio $H_0/R_0$	0.05	0.03–0.05	0.01	0.005–0.03	0.01	0.005–0.03

Table 3. Table 3: Best-fit model parameters of the J band data.

Optimized	East side ($\nu =425$)		West side ($\nu =635$)		Both sides ($\nu =1065$)
parameter	Best fit	Models with	Best fit	Models with	Best fit	Models with
	${\chi }_{\nu ,\min }^2=1.67$		${\chi }_{\nu ,\min }^2=2.22$		${\chi }_{\nu ,\min }^2=4.24$
HG parameter $g$	0.4	0.4, 0.5	0.3	0.3	0.3	0.3
Inner radial index ${\alpha }_{\mathrm{in}}$	3	2-6	3	2-5	2	2-3
Outer radial index ${\alpha }_{\mathrm{out}}$	-5	-5, -6	-3	-3	-3	-3
Aspect ratio $H_0/R_0$	0.03	0.03–0.05	0.03	0.01–0.05	0.03	0.075–0.03

Table 4. Table 4: Photometric measurements of HD 15115.

Wavelength	Flux density	Error	Instrument	Ref.
($\mu$m)	(mJy)	(mJy)
13	170.2	2.9	Spitzer/IRS	1
24	58.3	2.3	Spitzer/MIPS	2
31	85.6	2.3	Spitzer/IRS	1
70	451.9	32.6	Spitzer/MIPS	2
160	217.3	27.8	Spitzer/MIPS	2
850	4.9	1.6	JCMT/SCUBA	3

Equations

$$Q=0.5\ast (Q^{+}-Q^{-})$$

$$U=0.5\ast (U^{+}-U^{-}),$$

$$Q_{\phi }=-Q\cos 2\phi -U\sin 2\phi$$

$$U_{\phi }=Q\sin 2\phi -U\cos 2\phi$$

$$f_M(y)=a_M{\left[1+{\left(\frac{y-y_0}{\alpha }\right)}^2\right]}^{-\beta },$$

$$\begin{gathered} n(r,h)\sim {\left({\left(\frac{r}{R_0}\right)}^{-2{\alpha }_{in}}+ \\ {\left(\frac{r}{R_0}\right)}^{-2{\alpha }_{out}}\right)}^{-1/2}\cdot \exp \left[-{\left(\frac{|h|}{H(r)}\right)}^2\right], \end{gathered}$$

$$f(\theta )=\frac{1-g^2}{4\pi (1+g^2-2gcos(\theta ){)}^{3/2}},$$

$${\chi }_{\nu }^2=\frac{1}{\nu }\sum _{i=1}^{N_{\mathrm{data}}}\frac{{\left[F_{i,\mathrm{data}}-A\cdot F_{i,\mathrm{model}}(\text{@vec}p)\right]}^2}{{\sigma }_{i,\mathrm{data}}^2}$$

$$BIC(\mathcal{M})=-2\ln (\mathcal{L}(\widehat{\text{@vec}p}))+N_{\mathrm{par}}\ln (N_{\mathrm{data}}),$$

$$\ln (\mathcal{L})=-\frac{1}{2}{\chi }^2+const$$

$$p_{\mathrm{sca}}(\theta )=p_m\frac{1-{\cos }^2\theta }{1+{\cos }^2\theta },$$

$${\chi }_{\mathrm{spine}}^2=\frac{1}{N_{\mathrm{data}}}\sum _{i=1}^{N_{\mathrm{data}}}\frac{{\left[y_{\mathrm{i}\mathrm{data}}-y_{\mathrm{i}\mathrm{model}}\right]}^2}{{\sigma }_{\mathrm{i}\mathrm{data}}^2},$$

$$\theta =\arcsin \frac{\sqrt{x_c^2+y_c^2}}{R_0}$$

$${\theta }_{\mathrm{back}}=180^{\circ }-\arcsin \frac{\sqrt{x_c^2+y_c^2}}{R_0},$$

$$f(\theta )=w{f}_{\mathrm{HG}1}(\theta ,g_1)+(1-w)f_{\mathrm{HG}2}(\theta ,g_2),$$

$$B_{\nu ,\mathrm{mod}}(\lambda ,T_{\mathrm{dust}})=B_{\nu }(\lambda ,T_{\mathrm{dust}})\times (H({\lambda }_0-\lambda )+H(\lambda -{\lambda }_0)({\lambda }_0/\lambda {)}^{\beta }),$$

Full text

*Aims. *We present new observations of the edge-on debris disk around HD 15115 (F star at 48.2 pc) obtained in the near-IR. We search for observational evidence for a second inner planetesimal ring in the system.

*Methods. *We obtained total intensity and polarimetric data in the broad bands J and H and processed the data with differential imaging techniques achieving an angular resolution of about 40 mas. A grid of models describing the spatial distribution of the grains in the disk is generated to constrain the geometric parameters of the disk and to explore the presence of a second belt. We perform a photometric analysis of the data and compare disk brightness in two bands in scattered and in polarized light.

Results. We observe an axisymmetric planetesimal belt with an inclination of $85.8^{\circ}\pm 0.7^{\circ}$ and position angle of $\sim$278.9${}^{\circ}\pm 0.1^{\circ}$. The photometric analysis shows that the west side is $\sim$2.5 times brighter in total intensity than the east side in both bands, while for polarized light in the J band this ratio is only 1.25. We also find that the J - H color of the disk appears to be red for the radial separations $r\lesssim 2^{\prime\prime}$ and is getting bluer for the larger separations. The maximum polarization fraction is 15–20% at $r\sim$2.5′′. The polarized intensity image shows some structural features inside the belt which can be interpreted as an additional inner belt.

*Conclusions. *The apparent change of disk color from red to blue with an increasing radial separation from the star could be explained by the decreasing average grain size with distance.

The presence of an inner belt slightly inclined with respect to the main planetesimal belt is suspected from the data but the analysis and modeling presented here cannot establish a firm conclusion due to the faintness of the disk and its high inclination.

1 Introduction

1.1 Debris disks with multiple belts

Due to the discovery of hundreds of debris disks 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 a growing evidence from the infrared and submillimeter 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 the disk 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 disk (e.g. Lyra & Kuchner, 2013). The former mechanism assumes the presence of planets formed in the protoplanetary disk. In analogy to the solar system, the planets carve ring-like structures in the disk 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 disk still containing 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 disk 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 disks. Using the hydrodynamical simulation of the gas-disk 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 disk dust-to-gas ratio lower than unity. This condition might be fulfilled in the so-called transitional disks such as HD 141569A (Thi et al., 2014) or in the gas-rich debris disks 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 modeling can provide neither the actual number of belts in the system nor an accurate estimate of their location and spatial 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 disks 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 modeling.

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^{\circ}$. Obviously, it is much more difficult to identify more than one belt in edge-on disks. In this work, we investigate the presence of two belts of planetesimals in the well-known edge-on disk around HD 15115.

1.2 Previous observations of 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\pm 0.41$ mas corresponding to a distance of $48.16\pm 0.95$ pc which we apply in this work. This value is slightly different from the former determination of $45.2\pm 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 $\beta$ 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 disk in the optical and in the near-IR (H band). The disk was found to be nearly edge-on extending along the east-west direction with a west side detected out to a stellocentric radius of $\approx 580$ au ($\approx 12^{\prime\prime}$) and the east side to at least 340 au ($\approx 7^{\prime\prime}$).

The complexity of the inner disk 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 $\mu$m and claimed a possible “warp” on the western side of the disk at around $1^{\prime\prime}$. They also established a wavelength dependence of the brightness asymmetries. Kalas et al. (2007) measured an overall blue V - H color for the west side, which together with the shape of the disk gave it its nickname: “the blue needle”. Debes et al. (2008) detected a relatively intense variation of disk colors as a function of distance from the central star: at $2^{\prime\prime}$ they measured a neutral color up to 1.1 $\mu$m and a blue spectrum for $\lambda>1.1$ $\mu$m, while at $1^{\prime\prime}$ two disk extensions show a red F110W - H color. This suggests the presence of several populations of grains with various properties, such as grain size or composition, at different separations. The red color of the inner disk could be due to the strong absorption features of olivine at 1 $\mu$m (Debes et al., 2008). Rodigas et al. (2012) found, on the contrary, that the disk possesses mostly a spectrally neutral KS - L*′* color (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 KS band but the disk 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\approx 1.99^{\prime\prime}$ and a position angle (PA) of $98.8^{\circ}$ in H and KS band images from the Gemini Near-Infrared Coronagraphic Imager (NICI) archival data sets. They also found that the disk has an inclination of $\sim$86∘, shows no noticeable center offset with respect to the star and probably has an inner cavity. This work provided the first unbiased view of the disk structure though with some discrepancies with the results of recent observations of 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 $\sim 80^{\circ}$ for an inner “half-ring” while they estimate that the extended disk 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 disk side and a significantly larger vertical extension of the disk 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 millimeter emission showed a rather symmetric distribution of large grains in the disk 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_{\rm in}=43$ to $R_{\rm 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 disk as a double belt system based on the SED modeling. Moór et al. (2011) used the disk fluxes at $\lambda=24$ $\mu$m, $\lambda=70$ $\mu$m, and $\lambda=160$ $\mu$m (Table 4) as well as the Spitzer Infrared Spectrograph (IRS) 5-35 $\mu$m spectrum for their fit and found two belts with radii $R_{\rm warm\,dust}=4$ au and $R_{\rm cold\,dust}=42$ au. Chen et al. (2014) derived $R_{\rm warm\,dust}=5$ au and $R_{\rm cold\,dust}=175$ au based on the fits to the same Spitzer IRS spectrum and the IR excess at $\lambda=13$ $\mu$m, $\lambda=24$ $\mu$m, and $\lambda=31$ $\mu$m.

In this paper, we present new observations of HD 15115 in the broadbands 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 organized 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 modeling of the disk geometry using both the total intensity and the polarimetric data is presented in Sect. 5. Section 6 describes the measurement of stellar and disk 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 summarize 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-stabilized 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 ($\lambda_{c}=1626$ nm, $\Delta\lambda=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 ($\lambda_{c}=1258$ nm, $\Delta\lambda=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 summarized in Table 1. To measure the stellar flux, several non-coronagraphic frames were recorded before and after each observation with DIT $\approx$ 0.84 s in H band and DIT = 2.00 s in J band using a neutral density filter ND2.0 with a $\sim$10${}^{-3}-10^{-2}$ transmissivity to prevent the saturation of the detector.

At the beginning of the science observation a “center 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 stabilized. 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 32s per frame. The total exposure time of the polarimetric observations was 5120s. 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$\times$1024 pixels. In Classical Imaging (CI), raw frames contain two images side-by-side (left and right channel) simultaneously acquired with a broadband 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 broadband 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^{\prime\prime}\times 11^{\prime\prime}$.

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 cali- bration 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$\times$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$\times$1024 pixels each.

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

In the next step, reduced and centered 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 1d.

Similar results were obtained with the KLIP algorithm (Figs. 1b and 1e). 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 disk-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 disk flux.

For comparison, we also performed the LOCI data reduction (Figs. 1c and 1f). 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_{\delta}$. The optimization 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 disk. For the LOCI reduction we adopted parameters $N_{\delta}=1.25\times$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^{\circ}$ (Maire et al., 2016) during the ADI procedure.

Reduced frames of the stellar flux measurement are centered by fitting a 2D Gaussian function with a FWHM of $\sim$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\cdot 10^{-3}$ and $1.5\cdot 10^{-2}$ in the J and H bands, respectively (see SPHERE manual111https://www.eso.org/sci/facilities/paranal/instruments/sphere/doc.html). 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 the DRH 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 centering we have used the same method as for the centering of the total intensity data set (see Sect. 3.1). After that, each $Q$ and $U$ frame was visually inspected to check the correct position of the center of the right channel frame relative to the left channel frame.

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

[TABLE]

[TABLE]

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 disk. To generate the image of the disk 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_{\varphi}$ parameter and $U_{\varphi}$ parameter, the linear polarization component at a 45∘ angle to the $Q_{\varphi}$ parameter:

[TABLE]

[TABLE]

where $\varphi$ is the polar angle measured east of north in the coordinate system centered on the star. Here we adopt the sign convention for $Q_{\varphi}$ and $U_{\varphi}$ parameters which is consistent with $Q_{\rm r}$ and $U_{\rm 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_{\varphi}$ and $U_{\varphi}$ images are rotated clockwise by $1.7^{\circ}$ to perform the astrometric TN correction (see Sect. 3.1).

Figure 2 shows the $Q_{\varphi}$ image with the detected polarized light. The debris disk 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_{\varphi}$ image, as expected. Therefore the $Q_{\varphi}$ image shows the total intensity of the polarized light.

4 Data analysis

4.1 Disk morphology

Total intensity data

All data reductions (cADI, KLIP, LOCI) reveal an almost edge-on disk with the bright side to the north of the disk 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 disk ansae symmetrically located to the east and west of the star. In some images, parts of the fainter side of the disk are also visible near the ansae. The brighter disk side is detected with a signal-to-noise ratio (S/N) higher than 4 (see S/N map in Appendix, Fig. 19) up to a radial distance of $\sim$4*′′* eastwards and $6^{\prime\prime}$ on the west side (up to the end of frame).

Polarized intensity data

Figure 3 shows the binned $Q_{\varphi}$ 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 ($\theta_{\rm disk}$), $-x$ is eastwards ($\theta_{\rm disk}-180^{\circ}$), and the $y$- axis is perpendicular to this with the positive axis towards the north.

There are three surface brightness peaks on the east side within $r=1.5^{\prime\prime}$ in azimuthal direction: at A$(-0.80^{\prime\prime},\,0.29^{\prime\prime})$, B$(-1.30^{\prime\prime},\,0.05^{\prime\prime})$, and C$(-0.75^{\prime\prime},\,-0.16^{\prime\prime})$ 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 surface brightness 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 surface brightness 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 disk with single and double belt in comparison to the data. To create synthetic images, the model222The detailed description of the model used can be found in Sect. 5. with radius of the disk $R_{0}=2^{\prime\prime}$ and inclination of $86^{\circ}$ was used. The other parameters were arbitrarily chosen, for instance, the radius of the inner belt was set to $1.5^{\prime\prime}$. 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_{\varphi}$ image. The disk flux was previously removed from this $Q_{\varphi}$ image. This was done by replacement of the image area containing disk flux through the same area of the mean of two images obtained by rotation of the $Q_{\varphi}$ image by $15^{\circ}$ clockwise (first image) and by $15^{\circ}$ counterclockwise (second image). The model images were binned and smoothed in the same way as was done to obtain the $Q_{\varphi}$ image presented in Fig. 3. The visual comparison between model images (Fig. 4a and 4c) and the $Q_{\varphi}$ image (Fig. 4b) suggests that the HD 15115 debris disk 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\approx 1.25^{\prime\prime}\pm 0.07^{\prime\prime}$ and a PA $\theta_{\rm inn\,ring}=276^{\circ}\pm 3^{\circ}$, and thus a difference in PA up to 6*∘* with respect to 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 disk spine (Sect. 4.3), compare the models for the spatial distribution of dust in an optically thin disk with the total and polarized intensity data (Sect. 5), and investigate the variation of the disk surface brightness along the disk axis calculated in apertures of different sizes (Sect. 6).

4.2 Position angle and size of the disk

We estimated the disk radius, the PA, and the inclination from the H band imaging data reduced with LOCI and KLIP algorithms where the disk front and back sides are visible. We derived the $R_{\rm disk}=2.01^{\prime\prime}\pm 0.05^{\prime\prime}$, $\theta_{\rm disk}=278.7^{\circ}\pm 0.2^{\circ}$ and the disk inclination $i=86.5^{\circ}\pm 0.5^{\circ}$ by fitting an ellipse to the ring contours and the surface brightness peaks at the east and west disk edges.

We also derived the PA and the radius of the disk by determining the centroids of bright ansae at the radial separation $r\approx 2^{\prime\prime}$ to the east and west of the star in the polarimetric data in J band (Fig. 20). This location corresponds to the radius of the symmetrical ring $R_{\rm disk}=1.99^{\prime\prime}$ detected by Mazoyer et al. (2014) in the H band data which is evident in our data as well. For an edge-on disk or ring with a characteristic radius $R_{\rm disk}$, the maximum polarized flux is expected at the locations of the planetesimal belt or at scattering angles close to $90^{\circ}$. 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_{\rm disk}$ (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 disk the radius $R_{\rm disk}=1.96^{\prime\prime}\pm 0.05^{\prime\prime}$ ($94.39\pm 2.41$ au) along the disk major axis and the PA $\theta_{\rm disk}=279.0^{\circ}\pm 0.2^{\circ}$. After taking into account the PA measured in the total intensity data, this yields a mean PA of the disk $\theta_{\rm disk}=278.9^{\circ}\pm 0.1^{\circ}$.

These values are in good agreement with the disk radius $R_{\rm disk}=1.99^{\prime\prime}$, PA $\theta_{\rm disk}=278.8^{\circ}\pm 0.4^{\circ}$ and the disk inclination $i=86.2^{\circ}$ 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 $\theta_{\rm disk}=278.9^{\circ}\pm 0.2^{\circ}$ for the west side and PA $\theta_{\rm disk}=276.6^{\circ}\pm 1.3^{\circ}$ for the east side at 1.1 $\mu$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 disk

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

[TABLE]

where $a_{M}$ is the peak flux located at a vertical distance $y_{0}$ from the disk major axis and $\alpha$ and $\beta$ are two constant parameters determined from the fit.

We measured the spine position of the brighter side of the disk for several data reductions (cADI, KLIP, LOCI) in H and J bands (with a 4$\times$4 binning) and derived the mean vertical offset $\left\langle y_{0}(x)\right\rangle$ as plotted in Fig. 5. The $y_{0}$ values, which are larger than $y_{0}=100$ mas or smaller than $y_{0}=-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^{\prime\prime}$ and $3.7^{\prime\prime}$ on the west side and between $1.0^{\prime\prime}$ and $2.4^{\prime\prime}$ on the east side.

Between $r\approx 1.2^{\prime\prime}$ and $r\approx 2.2^{\prime\prime}$ the spine is roughly symmetric with respect to $x_{0}$. The smallest $\left\langle y_{0}(x)\right\rangle$-offset is approximately $45\pm 5$ mas (2.2 au) around $|x|=2.0^{\prime\prime}$ $\pm 0.1^{\prime\prime}$. For smaller radial separations the vertical offset is increasing and achieves a maximum of $\sim$140 mas at $x\approx\pm 1.0^{\prime\prime}$ in our data. For radial separations $x>2.0^{\prime\prime}$ the spine deviates from the major axis by 45–60 mas towards the north. The same behavior is observed on the east side, although with a lower S/N. A similar spine curve which does not reach the disk major axis at the location of the planetesimal belt was measured for another edge-on debris disk around F star HIP 79977 (Engler et al., 2017) and A star $\beta$ 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).

Spine diagnostics I

The turning point of the spine at $\sim$2*′′* clearly indicates the location of the edge of the parent belt of the planetesimals. The geometry of the disk 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 disk inclination. We derive these parameters by fitting an ellipse to the spine between $|x|\geq 1.2^{\prime\prime}$ and $|x|\leq 2.0^{\prime\prime}$. The center of the ellipse is set to the star position. The best-fit ellipse (Fig. 5) has a semi-major axis of $2.06^{\prime\prime}\pm 0.06^{\prime\prime}$ and a semi-minor axis of $0.15^{\prime\prime}\pm 0.01^{\prime\prime}$ yielding a disk inclination of $85.8^{\circ}\pm 0.7^{\circ}$. This inclination estimate is in good agreement with $i=86.2^{\circ}$ found by Mazoyer et al. (2014) with the same methodology.

5 Disk modeling

5.1 One belt model

To put constraints on the basic geometric parameters of the disk, we generate a set of models from a parameter grid and compare the synthetic images for the scattered light with the KLIP disk images in J and H bands (Figs. 1b and 1e) by evaluating the reduced $\chi^{2}_{\nu}$ for the model fit (Eq. 3).

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_{0}$ 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^{\alpha_{\rm in}}$ $(\alpha_{\rm in}>0)$ for $r<R_{0}$ and $r^{\alpha_{\rm out}}$ $(\alpha_{\rm out}<0)$ for $r>R_{0}$ 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 disk can be written in cylindrical coordinates $(r,h)$ as:

[TABLE]

where $H(r)$ is the disk scale height which scales with a radius like $H(r)=H_{0}\,{\left(\frac{r}{R_{0}}\right)}^{\gamma}$, $H_{0}$ is the scale height of the vertical profile of the disk at $R_{0}$, and $H_{0}/R_{0}$ defines the aspect ratio of the disk. For simplicity, we set the disk flare index to $\gamma=1$.

Further, we assume that the disk 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):

[TABLE]

where $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\leq 1$, while for $-1\leq 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 $\sigma_{\rm sca}$ averaged over all particle sizes (Engler et al., 2017). In our model, these parameters are constant through the disk, which is a reasonable approximation if we take into account that disk 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 modeling of the HD 15115 disk and includes the following parameter sets (a grid of 2160 models in total):

• •

power-law index $\alpha_{\rm in}$: 2, 3, 4, 5, 6, 7, 8, 9, 10

• •

power-law index $\alpha_{\rm out}$: -2, -3, -4, -5, -6, -7, -8

• •

asymmetry parameter $g$: 0.2, 0.3, 0.4, 0.5, 0.6

• •

disk aspect ratio $H_{0}/R_{0}$: 0.005, 0.0075, 0.01, 0.02, 0.03, 0.05.

The radius and inclination of the disk 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 disk 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-modeling approach to account for the effects of the disk self-subtraction which always arise when the ADI technique is used (Soummer et al., 2012). In forward modeling, 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 1e. The self-subtraction of the disk flux caused by the KLIP post-processing is clearly visible in the form of dark stripes on both sides of the disk in the data and the models.

To estimate the goodness of the model fit to the KLIP images, we define the reduced $\chi_{\nu}^{2}$ metric as follows:

[TABLE]

where $N_{\rm data}$ is the number of binned pixels used to perform the $\chi^{2}_{\nu}$ minimization, $F_{i,\,{\rm data}}$ is the flux measured in pixel $i$ with an uncertainty $\sigma_{i,\,{\rm data}}$, $F_{i,\,{\rm model}}(\@vec{p})$ is the modeled flux of the $i$ pixel and $A$ is a scaling factor. The degree of freedom of the fit is denoted by $\nu=N_{\rm data}-N_{\rm par}$, where $N_{\rm par}$ represents the number of variable parameters $\@vec{p}=(p_{1},p_{2},...,p_{N_{\rm par}})$ which is equal to five in this case.

A $\chi^{2}_{\nu}$ metric is evaluated for each KLIP-processed model of the grid with respect to the data presented in Figs. 1b and 1e. We applied a $4\times 4$-pixel binning of the data and the models. Because of the pronounced brightness asymmetry between both disk 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^{\prime\prime}<x<4^{\prime\prime}$ and $-0.2^{\prime\prime}<y<0.2^{\prime\prime}$) and for the east side ($-3^{\prime\prime}<x<-1^{\prime\prime}$ and $-0.2^{\prime\prime}<y<0.2^{\prime\prime}$) for reasons related to the S/N. We also considered both disk 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 disk flux.

We consider the model providing the minimum $\chi^{2}_{\nu,\,{\rm min}}$ as a best-fit model. Additionally, we determine the parameter sets of all well-fitting models bounded by a threshold of $\chi^{2}_{\nu}\leqslant\chi^{2}_{\nu,\,{\rm min}}+\Delta\,\chi^{2}_{\nu}$ (Thalmann et al., 2014) to take the systematics into account. The threshold is derived with $\Delta\,\chi^{2}_{\nu}=\sqrt{2/\nu}$ corresponding to a 1$\sigma$ deviation of the $\chi^{2}$ distribution with $\nu$ 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 $\chi^{2}_{\nu}$ 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 $\alpha_{\rm in}=3$ and $\alpha_{\rm out}=-4$ (H band data) or $\alpha_{\rm out}=-3$ (J band data), whereas the east side is best represented with $g=0.5$ and a much steeper index $\alpha_{\rm in}=8$ (H band data) and $\alpha_{\rm out}=-6$. In the J band data, there is no large difference for the inner radial index $\alpha_{\rm in}$ between the east and the west side. The aspect ratio of 0.03 between disk radius and width seems to give a good fit to both disk sides, in particular when modeling J band data.

As a final remark on the modeling results, we note that the $\chi^{2}$-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 disk could be explained by the presence of an inner belt of planetesimals with a radius of $R<2^{\prime\prime}$. 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 disk. The phase function in total intensity being different than in polarimetry, an inner component will dominate the disk image at very small phase angles precisely where the ADI post-processing strongly attenuates the disk signal. To further test the disk geometry we consider both one- and two-belts model 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 ($\chi_{\nu}^{2}=1.81$ for the west side and $\chi_{\nu}^{2}=1.40$ for the east side). The model is generated with the following parameters: $R_{0}=2^{\prime\prime}$, $i=86^{\circ}$, $g=0.3$, $\alpha_{\rm in}=5$, $\alpha_{\rm out}=-4$, $H_{0}/R_{0}=0.01$ and position angle $\theta_{\rm out\,belt}=279^{\circ}$. Figure 6 shows the appearance of the modeled 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 modeling 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 $R_{\rm inn\,belt}=1.3^{\prime\prime}$ in accordance with the location of the surface brightness (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 $\alpha_{\rm in}=5$ and $\alpha_{\rm 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 $R^{2}_{\rm out\,belt}/R^{2}_{\rm inn\,belt}$ 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 $\theta_{\rm out\,belt}=279^{\circ}$) has an inclination of $i=86^{\circ}$; the inner belt has a PA of $\theta_{\rm inn\,belt}=276^{\circ}$ as measured from the polarimetric data (see Sect. 4.1) and an inclination of $i=80^{\circ}$. Such inclination was estimated by Schneider et al. (2014) for the inner “half-ring” component of the HD 15115 disk (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 7c. Because the relative PA and inclination of belts are non-zero, we expect the maximum SB of the disk 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 disk 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 disk. 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 $\chi^{2}_{\nu}$ for the models presented in Figs. 7d and 8b setting the intensity ratio between the inner and outer belt equal to 0.5, 1, 2, 3 and 4. A simultaneous fit of both disk sides yields the $\chi^{2}_{\nu}$ in the range 3.38 and 3.67 for the non-coplanar model and in the range 3.68 and 3.89 for the coplanar model ($\nu=1060$, $N_{\rm par}=10$). The minimum $\chi^{2}_{\nu}$ is obtained for the 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-belt models, 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 using the 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 $\mathcal{L}$ of the model $\mathcal{M}$. We define the $BIC$ as follows (Schwarz, 1978):

[TABLE]

where $\hat{\@vec{p}}$ is a set of the best-fitting parameters for the model $\mathcal{M}$ providing the maximized value of the likelihood function $\mathcal{L}$. Under the assumption that the model errors are normally distributed, this function can be estimated as:

[TABLE]

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 disk flux was removed in the same way as was done with the $Q_{\varphi}$ image (see Sect. 4.1). Because all models are rotationally symmetric, their intensity is adjusted to match the intensity of the west disk side for an easier comparison with the asymmetric HD 15115 disk. 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 disk, for instance. On the contrary, the two-belts coplanar model (Fig. 9d) does not show the disk opening (back side is invisible).

We conclude that the two-belts model with the inner ring parameters $R_{\rm inn\,belt}=1.3^{\prime\prime}$, $i=80^{\circ}$ and PA $\theta_{\rm inn\,belt}=276^{\circ}$ provides a better fit to the H band data than the coplanar model or single-belt model according to the $\chi^{2}_{\nu}$ method and the $BIC$. However, we have to keep in mind that this conclusion is based on the $\chi^{2}$ 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.

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 disk. Due to the different angular dependence of the scattered phase function and polarized phase function on scattering angle, the distribution of the disk SB and therefore the appearance of the disk are different in the images of scattered and scattered polarized light. For a nearly edge-on disk, the polarized disk flux has a maximum at scattering angles close to $\theta=90^{\circ}$ (e.g. Engler et al., 2017) while the scattered flux peaks at $\theta=0^{\circ}$ due to diffraction of light. The polarimetric measurement traces the extent of the planetesimal belt because we measure a peak polarized surface brightness of the disk at the location of the belt ansae. Therefore, in the case of the edge-on disk with multiple planetesimal belts, we should measure several SB peaks along the disk axis. These peaks indicate the presence of the inner belts and their extent.

In order to demonstrate this behavior, we simulated images of polarized light from an edge-on disk consisting of either one or two belts (Figs. 10a and 10b). 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_{\rm sca}$ has an angular dependence with the scattering angle $\theta$ such as that for Rayleigh scattering:

[TABLE]

where $p_{m}$ is the maximum fractional polarization achieved at a scattering angle of $\theta=90^{\circ}$.

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 Fig. 10a to Fig. 10b 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 disk sides. In this case, the radial SB profile measured along the disk major axis reflects the location and the number of the planetesimal belts (see insets in Figs. 10a and 10b).

Compared to the disk 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 disk. 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 disk structure based on the modeling results, a higher S/N of the polarimetric data is needed. Therefore we do not perform a fitting of the disk models to the $Q_{\varphi}$ image to explore the whole parameter space but compare qualitatively (visually) and quantitatively ($\chi^{2}_{\nu}$-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_{\varphi}$ image which contains no disk 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 disk in comparison to the data and to the two-belts model (Fig. 10e) and therefore the latter matches the $Q_{\varphi}$ 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 $\chi^{2}_{\nu}=3.49$ contrary to the one-belt model with $\chi^{2}_{\nu}=3.78$.

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 disk 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 $\alpha_{\rm in},\alpha_{\rm out}$, the HG asymmetry parameter $g$, and aspect ratio $H_{0}/R_{0}$ (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 $\alpha_{\rm in}$, $\alpha_{\rm out}$, $g$, $H_{0}/R_{0}$ or inclination varies. For this analysis, we use the KLIP processed grid models with the fixed radius of the belt ($R_{0}=2^{\prime\prime}$) as described in Sect. 5.

Figure 11 illustrates the degeneracy between $\alpha_{\rm in}$, $\alpha_{\rm out}$, $g$, $H_{0}/R_{0}$, and disk 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 $\alpha_{\rm in}$), the higher the spine elevation above the midplane inside the parent belt, whereas the steeper fall-off of the grain number density (decreasing $\alpha_{\rm 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 disk aspect ratio $H_{0}/R_{0}$ on the position of the spine is similar to that of the asymmetry parameter - the larger the disk scale height, the larger the spine elevation (Fig. 11d). The vertical offset of the spine from the disk axis can also be affected by the disk 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 disk (Fig. 11e).

Interestingly, we note that the curvature of the spine in the outer region (beyond $r=2^{\prime\prime}$) changes with decreasing $\alpha_{\rm out}$: the spine is bowed for $\alpha_{\rm out}\leqslant-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^{\prime\prime}$ but the S/N of the disk 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 disk radius is fixed. Therefore, as an alternative to calculating a 2D $\chi^{2}$ on scattered light images, which implies a large number of degrees of freedom, we used the spine itself and derived $\chi^{2}_{\rm spine}$ for angular separations in the range $1-3^{\prime\prime}$ on the west side of the disk as already done in Boccaletti et al. (2012) for the HD 32297 debris disk. We define the $\chi^{2}_{\rm spine}$ according to:

[TABLE]

where $N_{\rm data}$ is the number of data points between $1^{\prime\prime}$ and $3^{\prime\prime}$, $y_{\rm{i\,data}}$ are the vertical offsets of the spine measured from the data with uncertainties $\sigma_{\rm{i\,data}}$ (see Sect. 4.3) and $y_{\rm{i\,model}}$ are the vertical offsets of the model spine.

The minimal value $\chi^{2}_{\rm spine}=2.47$ is obtained for $\alpha_{\rm in}=4$, $\alpha_{\rm out}=-7$, $g=0.6$, $H_{0}/R_{0}=0.03$ and $incl=86^{\circ}$. However, there are other combinations of parameters, in particular the model with $\alpha_{\rm in}=3$, $\alpha_{\rm out}=-7$, $g=0.5$ with $\chi^{2}_{\rm spine}=2.53$, 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 ($\chi^{2}_{\rm spine}=3.10$), where the inner belt is inclined at $80^{\circ}$ (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.

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 centered on the PSF with a radius of $r=0.6^{\prime\prime}$ ($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\pm 0.80)\cdot 10^{7}$ ADU/s in J band and $(6.52\pm 0.78)\cdot 10^{7}$ ADU/s in H band. To measure the disk 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 disk flux is corrected for the ADI-induced flux reduction. The total disk flux is integrated in two rectangular areas containing pixels with ($x,\,y$) coordinates lying within the bounds $1^{\prime\prime}<|x|<4^{\prime\prime}$ and $-0.14^{\prime\prime}<y<0.26^{\prime\prime}$. We obtain 2900 $\pm$ 300 ADU/s and 4100 $\pm$ 500 ADU/s in J and H, respectively. This yields a brightness ratio between the disk and the star of $(4.0\,\pm\,0.6)\cdot 10^{-5}$ in the J band and $(6.3\,\pm\,1.0)\cdot 10^{-5}$ in the H band for scattered light.

The polarized flux measured in the same areas in J band is 240 $\pm$ 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 disk amounts to 840 $\pm$ 240 ADU/s. With this value, the polarized light contrast of the disk in J band is $(1.2\pm 0.6)\cdot 10^{-5}$. The ratio between the intensity of the polarized light and total intensity of the scattered light is $0.29\pm 0.09$.

In the total intensity we measured that the east side flux is $\sim$40% of the west side flux in both J and H band data (more precise $\sim$38% and $\sim$44%, respectively) while in polarimetry $F_{\rm east}=0.8F_{\rm west}$.

6.1 Disk color

Total intensity and polarimetric SB profiles are displayed in Fig. 12 for comparison. Here, the disk flux is measured in square $0.25^{\prime\prime}\times 0.25^{\prime\prime}$ (20 $\times$ 20 pixels) apertures centered on the disk axis. This aperture size was chosen for a better comparison with the radial SB profiles of the HD 15115 disk in V and H bands presented by Kalas et al. (2007) (their Fig. 3) and in 1.1 $\mu$m band published by Debes et al. (2008) (their Fig. 3).

For stellocentric distances ranging between $1.5^{\prime\prime}$ and $3.0^{\prime\prime}$ we find that the west side is $\sim$ 0.2 - 0.5 mag$\cdot$arcsec*-2* brighter than the east side in the J band while this difference is as large as $\sim$ 0.7 - 1 mag$\cdot$arcsec*-2* in the H band.

The steeper fall-off of the SB in the H band ($\propto$ $r^{-3.3}$ in the range $1.2-4^{\prime\prime}$ on the west side and $\propto$ $r^{-5.9}$ in the range $1.2-3^{\prime\prime}$ on the east side) compared to the J band ($\propto$ $r^{-2.6}$ in the range $1.2-4^{\prime\prime}$ on the west side and $\propto$ $r^{-3.5}$ in the range $1.2-3^{\prime\prime}$ on the east side) produces a change of color at a particular separation to the star. We observe a red or grey color (J-H $\geq$ 0) at short separations, while the color turns to blue (J-H$<$0) at large separations. This color change occurs between 2.5 and $3.5^{\prime\prime}$ 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$\cdot$arcsec*-2*.

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 $\propto$ $r^{-3.4}$ for $r>2^{\prime\prime}$. 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^{\prime\prime}\times 0.25^{\prime\prime}$ apertures along the disk 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 disk sides show a similar behavior increasing from PF $\approx$ 0.03 at $r=1.25^{\prime\prime}$ up to $\sim$15-20% between $r=2^{\prime\prime}$ and $r=3^{\prime\prime}$, 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 $\beta$ Pictoris debris disk. A similar fractional polarization ($\sim$10-25%) is obtained by Asensio-Torres et al. (2016) for another edge-on disk HD 32297 in the H band.

At radial separations of $r\gtrsim 3^{\prime\prime}$, 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 disk around M star AU Mic (Graham et al., 2007).

6.3 Surface brightness profiles

The aperture size $0.25^{\prime\prime}\times 0.25^{\prime\prime}$ ($20\times 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^{\prime\prime}$. In Fig. 14, we used a finer sampling of $4\times 4$ pixels ($0.049^{\prime\prime}\times 0.049^{\prime\prime}$) together with a $\sigma=1$ pixel and $\sigma=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 centered 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\approx 2^{\prime\prime}$ from the star on both disk 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_{\rm east}=2.00^{\prime\prime}\pm 0.05^{\prime\prime}$ and for the west side $R_{\rm west}=1.91^{\prime\prime}\pm 0.05^{\prime\prime}$. The small difference in the belt radius between both sides indicates an offset of $\sim 0.045^{\prime\prime}$ of the belt towards the east with respect to the star.

Further inside the ring, the SB profiles exhibit another peak at $r\approx 1.2-1.3^{\prime\prime}$ also on both disk 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 disk (’F3’ in Fig. 14c) but at a lower significance level.

The increased surface brightness at the locations ’F1’, ’F2’, and ’F3’ is also visible in the disk 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.

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 disk. We assume that grains scatter stellar radiation preferentially in the forward direction, this means the northern part of the disk is the near side. We also consider that the circumstellar material is confined to a narrow ring with a radius of $R_{0}=1.95^{\prime\prime}$ (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 disk side in circular apertures with a radius of 6 pixels. The disk area within each aperture was attributed a unique scattering angle which was calculated for the near side of the disk according to:

[TABLE]

and for the back side of the disk

[TABLE]

where $(x_{c},\,y_{c})$ are coordinates of the aperture center in the coordinate system shown in Fig. 3a. The apertures were centered on the spine of the front and back side of the disk defined as the brightest points of the vertical cross-sections perpendicular to the disk 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 modeling 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 $\theta<50^{\circ}$ and gives rise to a relatively broad diffraction peak. The measured curve is rather flat between $60^{\circ}$ and $90^{\circ}$, whereas for larger angles corresponding to the fainter back side of the disk ($\theta>90^{\circ}$) 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 used the $\chi^{2}$-minimization algorithm to find a HG function (Eq. 2) 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 $\chi_{\nu}^{2}$-criterion, the better fit $f(\theta)$ is given by a linear combination of two HG functions (Engler et al., 2017):

[TABLE]

with $g_{1}=0.94\pm 0.06$, $g_{2}=0.21\pm 0.01$ and scaling factor $w=0.89\pm 0.10$ (blue curve in Fig. 16). For this combination, the reduced $\chi_{\nu}^{2}$ metric is 1.22. This is smaller than $\chi_{\nu}^{2}$ of 1.43 which is obtained for the best-fit of the HG function with the single asymmetry parameter $g=0.57\pm 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 significantly improve 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 $\mu$m which are covered by small regolith particles.333In the original paper by Min et al. (2010), grain diameters are given in the notation instead of grain radii. 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 $\theta=83^{\circ}$.

At large scattering angles ($\theta>90^{\circ}$) the measured phase function of the HD 15115 disk does not follow the phase functions obtained by Min et al. (2010). When comparing the diffraction parts of the phase functions ($\theta<80^{\circ}$), the scattering function for grains with a size of 5 $\mu$m (solid black curve) fits best the measured intensity profile. Assuming that this grain size is the average size $\langle a\rangle$ of the grain size distribution, we derive the minimum grain size $a_{\rm min}=3$ $\mu$m from the relation $a_{\rm min}=3/5\,\langle a\rangle$ 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 $\sim$1–3 $\mu$m assuming a stellar luminosity of 3.5$L_{\odot}$, a mass of 1.3$M_{\odot}$, and a grain average density in the range from 1 to 3 g cm*-3* (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 disks around HR 4796 A and Fomalhaut shows the compliance with the curves for grains with radii of $a=15$ $\mu$m (Milli et al., 2017) and $a=50$ $\mu$m (Min et al., 2010), respectively. This could suggest that the grains in the HR 4796 A and Fomalhaut disks are larger on average than the dust grains we observe on the west side of the HD 15115 debris disk.

7 Discussion

7.1 Unpolarized versus polarized light

As mentioned in Sect. 6, the polarized intensity (Stokes $Q_{\phi}$) measured for the east side of the disk 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 along the line of sight on the west side in front of the planetesimal belt. With polarimetry, we trace predominantly $90^{\circ}$ 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_{\phi}$ signal.

An alternative interpretation of the disk brightness asymmetry was recently proposed by Lee & Chiang (2016). They presented a model explaining the “needle”-like morphologies observed in the edge-on disks by a planet perturbation. Their model requires an eccentric parent belt with a brighter disk 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 disk is caused by planets, then the perturbation mechanism should be more complex than the model suggested by Lee & Chiang (2016).

7.2 Color

The photometric analysis (Figs. 12 and 14) suggests that inside the planetesimal belt ($r\lesssim 2^{\prime\prime}$) the disk color is red and becomes bluer with increasing distance from the star ($r>2^{\prime\prime}$). A similar color behavior has previously been observed for the HD 15115 disk by Kalas et al. (2007) and Debes et al. (2008). They argued that the change in the dust scattering efficiency beyond $r\approx 2^{\prime\prime}$ could be due to the absorption features of dust material at 1.6 $\mu$m and inside $r\approx 2^{\prime\prime}$ due to absorption of olivine at $\sim$1 $\mu$m, for instance. Another possibility to explain the complex color of the disk 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_{\rm 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)\propto 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 disk color bluer, with distance from the star.

7.3 Evidence for multiple belts?

Our images show an axisymmetrical disk 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 KS filters, the disk has a symmetric ring shape with a radius of $R_{0}\approx 2^{\prime\prime}$ and an inclination of $\sim$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 disk 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_{\varphi}$ image (Fig. 2b) and ripples of the SB profile along the disk axis, in particular on the east side of the disk (Fig. 14) are similar to the ripples in the SB profile of the polarized light detected by Tamura et al. (2006) in the $\beta$ Pic debris disk. The 2 $\mu$m polarimetry of the $\beta$ Pic disk has revealed many details in the inner part of the disk, including the previously observed warp. In analysing the disk 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\approx 2^{\prime\prime}$) of the HD 15115 disk. Therefore the other SB peaks, which we measure further inside the disk, may also indicate inner rings.

The modeling 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 $\chi^{2}$ criterion for the 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 disk sides simultaneously (see Table 2).

In Appendix C, we show that the two-belts model is also consistent with the spectroscopic and photometric data of HD 15115, if the modified blackbody function (Backman & Paresce, 1993) is used to fit the disk 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 disk SED can be 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 disk side with decreasing distance to the star for $r<1.6^{\prime\prime}$ (see their Figs. 2a and 4a, respectively) which is in contradiction to a bow-like disk shape reported by Rodigas et al. (2012) and Schneider et al. (2014). Both authors suggested a nearly edge-on disk model with a high HG scattering asymmetry parameter to explain the ”bowing” of the brighter (northern) disk side. In accordance with this model, the PA of the east disk side should decrease following the course of the disk spine when approaching the star. This means that the vertical offset of the spine from the disk 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 $\mu$m (Debes et al., 2008, their Fig. 2a) and in KS band for $r<1.6^{\prime\prime}$ (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 disk structure. However, our data also show a small decrease of the spine offset between $r=1.5^{\prime\prime}$ and $r\approx$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\approx 1.4^{\prime\prime}$ on the east side of the disk: from a power law index of -1.1 for $r>1.4^{\prime\prime}$ to a power law index of -5.27 for $1.2^{\prime\prime}<r<1.4^{\prime\prime}$ (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\approx 1.5^{\prime\prime}$ which could be a planetesimal ring. In Fig. 17, we reproduce Fig. 1a from Rodigas et al. (2012) presenting a KS band image of the debris disk; it shows the increasing surface brightness below the disk 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 disk 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\sim$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 $\sim$0.2*′′* (clump width). Such a position shift on a circular orbit of the radius $r=2^{\prime\prime}$ for a 1.3$M_{\odot}$ star would be observable in about 18.5 years.

8 Summary

In this paper, we present classical imaging and polarimetric observations of the HD 15115 debris disk known as the “blue needle” with SPHERE-IRDIS in the broad bands 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_{\varphi}$, Fig. 2) in the broad band J were obtained using the double difference method.

Our results from the analysis of the disk geometry and photometry can be summarized as follows:

• •

The overall disk 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 disk towards the east is $R_{\rm east}=2.00^{\prime\prime}\pm 0.05^{\prime\prime}$ ($\sim$96 au) and towards the west it is $R_{\rm west}=1.91^{\prime\prime}\pm 0.05^{\prime\prime}$ ($\sim$92 au). By fitting an ellipse to the spine curve we derive a disk inclination of $85.8^{\circ}\pm 0.7^{\circ}$.

• •

The PA of the disk major axis measured in the total intensity and polarized intensity data is $\theta_{\rm disk}=278.9^{\circ}\pm 0.1^{\circ}$.

• •

The disk photometry and radial SB profiles along the disk major axis show the east-west asymmetry. In two rectangular areas bounded by the coordinates $1^{\prime\prime}<|x|<4^{\prime\prime}$ and $-0.14^{\prime\prime}<y<0.26^{\prime\prime}$ in the total intensity data of J and H bands, we measure $\sim$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 $\sim$80% of the flux measured on the west side (Sect. 6). This might indicate a different average grain size in the disk to the east compared to the disk to the west due to the loss of small dust grains towards the west.

• •

Based on the forward modeling results, we measure a change in the disk color between 2*′′* and 3.5*′′* from red (for $r<2^{\prime\prime}$) to predominantly blue (for $r<3.5^{\prime\prime}$). This color 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 disk.

• •

The image of the polarized light from the disk exhibits some SB peaks which can also be seen in the radial SB profile along the disk axis (Fig. 14c). Similar SB peaks can be recovered in the total intensity data (Figs. 14a and 14b). 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^{\prime\prime}$ and $1.5^{\prime\prime}$ 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 $\chi^{2}$ 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 disk with ZIMPOL/SPHERE polarimetry. An observation at shorter wavelengths would provide a required high resolution of $\sim$25-30 mas when observing in Very Broad Band (VBB, $\lambda_{c}=735$ nm) or I band ($\lambda_{c}=790$ nm). The Strehl ratio is higher in the I band but the count rate is $\sim$40% lower compared to the VBB. Under good observing conditions (seeing $\leqslant 0.7$, coherence time $\geqslant 3.5$ ms) a high contrast of $10^{-6}$ can be achieved. Polarimetry as a differential technique has also a decisive advantage over ADI in that it does not suffer from the disk self-subtraction, which could be crucial if we want to study the close-in disk structure.

Acknowledgements.

We would like to thank the referee for many thoughtful comments which helped to improve this paper. This work is supported 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 19 demonstrates the S/N of the total intensity images in J and H bands.

Appendix B Position angle of the disk

See Fig. 20.

Appendix C Consistency of the two-belts configuration with the disk SED

To check the consistency of the two-belts model with the available spectroscopic and photometric data of the HD 15115 disk, we fitted the disk SED with a single temperature blackbody (BB) emission $B_{\nu}(\lambda,T_{\rm dust})$ and with the modified Planck function (MBB, Backman & Paresce 1993):

[TABLE]

where $B_{\nu}(\lambda,T_{\rm dust})$ is the Planck function defined by the wavelength $\lambda$ and the dust temperature $T_{\rm dust}$, $H$ is the Heaviside step function, $\lambda_{0}$ is the characteristical wavelength, above which the absorption efficiency of dust grains scales like $(\lambda_{0}/\lambda)^{\beta}$ with $\beta$ being the opacity index. The free parameters $\lambda_{0}$ and $\beta$ allow to account for the steeper falloff of the debris disk 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_{\rm dust}$, the opacity index $\beta$, the characteristic wavelength $\lambda_{0}$ 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_{\rm dust}=39$ K for the one-belt model, and additionally a radius of the second belt of 1.3*′′* (61 au) with $T_{\rm 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 disk SED is constructed using the sampled 5–35 $\mu$m Spitzer/IRS spectrum (Moór et al. 2011) and fluxes compiled from the literature (Table  4). The stellar photosphere emission of a star is calculated with $L_{\ast}=3.5L_{\odot}$ and $T_{\rm eff}=6780K$ (Moór et al. 2011).

The derived best fits for each model are plotted in Fig. 21. 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 $\mu$m (Chen et al. 2014) which also disagrees with the Spitzer/IRS spectrum (AOR 10885632) retrieved from the archive.