© The Authors 2023

1 Introduction

Thanks to numerous mid- and far-infrared surveys conducted over the last two decades, debris disks have been discovered around hundreds of young stars (age<100Myr; Hughes et al. 2018, and references therein). These discoveries have brought about unprecedented possibilities for us to study the native environments of stars with newly born and still rapidly evolving planetary systems.

Based on these studies, we expect nearly each star to be surrounded by some amount of dust and gas, depending on its mass and evolutionary stage. The formation pathway of a star from a molecular cloud to a main-sequence star is accompanied by a transition of primordial gas and dust, which build a protoplanetary disk in the first five million years, to a processed (secondary) disk referred to as a debris disk. Such a disk contains dust particles and sometimes gas generated by collisions of planetesimals or evaporated from the surfaces of large icy bodies. This secondary material mostly consists of silicates, carbonaceous minerals, and H₂O, CO, and CO₂ ices and gases and provides building blocks for the growing cores of young terrestrial planets (Wyatt 2018).

While at stellar ages of ten million years the formation of gas giants is almost complete, this process takes longer for terrestrial planets (e.g., Rice & Armitage 2008; Morbidelli 2007). They are expected to grow over the next few hundred million years, collecting material on their orbits around a star, and can experience violent collisions with other planetary bodies from time to time.

The SpHere INfrared survey for Exoplanets (SHINE; Vigan et al. 2021; Desidera et al. 2021; Langlois et al. 2021), conducted with the Spectro-Polarimetric High-contrast Exoplanet REsearch (SPHERE) imager (Beuzit et al. 2019), was designed to search for and study young planets at early stages of their evolution. Many debris disks have been observed in the course of the SHINE program (Engler et al. in prep.) since they are the integral components of planetary systems. One of them is HD 114082, a debris disk that was resolved for the first time in the course of the SPHERE High Angular Resolution Debris Disk Survey (SHARDDS; Wahhaj et al. 2016; Milli et al. 2017a; Choquet et al. 2018, PI: J. Milli), a survey that targeted cold nearby debris disks at less than 100 pc from the Sun and with high infrared (IR) excess (L_ir/L_* > 10⁻⁴).

HD 114082 is an F3V star (Houk & Cowley 1975) located at a distance of 95.1 ± 0.2 parsecs (Gaia Collaboration 2021) in the Scorpius-Centaurus (Sco-Cen) association and a member of the Lower Centaurus Crux (LCC) subgroup. The age of the young stars in the LCC subgroup is estimated to be ~17 Myr (Mamajek et al. 2002).

The circumstellar dust around HD 114082 was first detected at far-IR wavelengths (Chen et al. 2011) from observations with the Multiband Imaging Photometer for Spitzer (MIPS; Rieke et al. 2004). In these observations, the star showed a very large excess flux over the predicted stellar flux at 24 and 70 μm, yielding a fractional luminosity of L_ir/L_* of 3.01 × 10⁻³ (Jang-Condell et al. 2015).

HD 114082 was also targeted by the Gemini Planet Imager (GPI; Esposito et al. 2020) and several times by the Ata-cama Large Millimeter/submillimeter Array (ALMA; e.g., Lieman-Sifry et al. 2016). In the GPI images of total intensity, the disk’s western extension appears to be 0.2″ longer and up to 80% brighter than the eastern extension. In observations at millimeter wavelengths, the debris disk was marginally resolved (Kral et al. 2017), suggesting a spatial extent of the disk of less than 1.69″ and a disk flux of 0.72 ± 0.05 mJy at 1.27 mm (Kral et al. 2020). CO gas was not detected.

This work presents new optical and near-IR observations of the HD 114082 debris disk performed with the SPHERE InfraRed Dual-beam Imager and Spectrograph (IRDIS; Dohlen et al. 2008), the Integral Field Spectrograph (IFS; Claudi et al. 2008), and the Zurich Imaging POLarimeter (ZIMPOL; Schmid et al. 2018), along with an analysis of the disk structure and the scattering properties of its dust particles. We also report the detection of a transiting planet in the cavity of this disk in the time series data taken with the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2014).

In Sects. 2 and 3, we describe the acquisition of the SPHERE data and their reduction process. Section 4 discusses the disk morphology observed in the images of total and polarized intensities and presents the results of the disk modeling. Section 5 addresses the scattering phase functions (SPFs) and the polarization fraction (PF) of scattered light. In Sect. 6, we discuss scattering properties of dust particles by analyzing the disk reflectance spectrum and comparing it to the spectrum of the HD 117214 debris disk. HD 117214 is another young F star in the LCC subgroup that harbors a bright debris disk (Engler et al. 2020) with a fractional luminosity (2.53 × 10⁻³; Jang-Condell et al. 2015), and physical disk size that are close to those of the HD 114082 disk. The optical properties of dust particles in both disks could therefore be similar. So far, HD 117214 is the only F star from the same sky region that has been observed under excellent observing conditions, making a comparison of disk reflectance spectra possible.

Finally, in Sect. 7 we report on the discovery of a transiting planet in the disk and on companion candidates detected in the HD 114082 images, and present the sensitivity limits on bound companions for a combination of direct imaging and radial velocity (RV) data with Gaia astrometry. Our results are summarized in Sect. 8.

2 Observations

We present in this paper new SHINE and SHARDDS observations of HD 114082 with SPHERE in three different observing modes: (i) IRDIFS observations, (ii), IRDIS polarimetric observations, and (iii) ZIMPOL polarimetric observations. The polarimetric observations were part of the SHARDDS survey (Dahlqvist et al. 2022) as a follow-up and characterization strategy of the disk. They complement archival data of the disk acquired with the classical imaging (CI) mode of IRDIS, using the broadband H filter and the apodized Lyot coronagraph N_ALC_YJH_S (diameter of 185 mas, Carbillet et al. 2011; Guerri et al. 2011). The log and details of all observations are presented in Table 1.

2.1 IRDIFS observations

The IRDIFS observations of HD 114082 were performed on 2017 May 17 in the IRDIFS-EXT mode, which provides a simultaneous data acquisition with IRDIS and the IFS. We used the pupil-stabilized mode to record a sequence of images consisting of 16 × 6 × 64 s IRDIS exposures and 16 × 5 × 64 s IFS exposures while the field of view was rotating. The observation started shortly before the meridian passage and covered a total range of parallactic angle variation of 39° in the IRDIS data set (33° in the IFS data because of shorter duration of the IFS observation). The detector integration time (DIT) per frame, the total exposure time, and the observing conditions can be found in Table 1.

The IRDIS data were taken in the dual-band imaging (DBI) mode (Vigan et al. 2010) using the filter pair K1K2 (λ_K1 = 2.110 μm, Δλ_K1 = 0.102 μm; λ_K2 = 2.251 μm, Δλ_K2 = 0.109 μm). The IFS was operated in the Y-H mode (0.95–1.66 μm, with a spectral resolution R_λ = 35). The field of view (FOV) of the IRDIS detector is approximately 11” × 12.5”, and that of the IFS is 1.73” × 1.73”. The apodized Lyot coronagraph N_ALC_Ks (diameter of 240 mas, Carbillet et al. 2011; Guerri et al. 2011) was used to block the stellar light.

To measure the stellar flux, several short exposures, where the star was offset from the coronagraphic mask, were taken at the beginning of the observation using a neutral density (ND) filter ND1.0 with a transmission of about 15% (see the SPHERE User Manual¹). The DIT of these flux calibration frames is 2 s (DIT = 4sforthe IFS).

Additionally, the “center frames” were taken before and after the science sequence using the deformable mirror waffle mode. They provided a measurement of the star position behind the coronagraph.

2.2 IRDIS polarimetrie observations

To extract the near-IR linearly polarized intensity of scattered light, we observed the disk with the dual-beam polarimetric imaging (DPI; de Boer et al. 2020; van Holstein et al. 2020) mode of IRDIS in the broadband H filter (λ = 1.625 μm, Δλ = 0.290 μm) using the apodized Lyot coronagraph N_ALC_YJH_S. These polarimetric observations were field-stabilized, because pupil-stabilized observations were not yet offered in DPI at the epoch of our observations.

A raw IRDIS frame has format of 2048 × 1024 pixels and contains two side by side star images provided by a beam splitter in front of the instrument detector. In DPI mode, two polarizers in the optical paths downstream from the beam splitter ensure that the images of two opposite polarization states are simultaneously recorded on the left and right halves of detector. The Stokes parameters +Q, -Q, +U, and -U are measured using a half-wave plate (HWP) inserted upstream from the beam splitter and rotated by 0°, 45°, 22.5°, and 67.5°, respectively.

We performed seven polarimetric cycles of consecutive measurements of the Stokes parameters +Q, -Q, +U, and -U consisted of five exposures with a DIT of 32 s for each parameter. The stellar flux was obtained with an off-axis image of the star using a ND filter, ND2.0, with a 1.5% transmissivity in the broadband H filter (see the SPHERE User Manual) and a DIT of 4s.

2.3 ZIMPOL polarimetric observations

The SPHERE instrument also has an optical arm, equipped with the ZIMPOL polarimeter. We present here polarimetric observations obtained in the I_PRIME filter (λ = 0.790 μm, Δλ = 0.153 μm) in the field-stabilized polarimetric mode (P2 mode) of the instrument using the slow polarimetry mode of detector. The FOV of ZIMPOL is 3.6” × 3.6”.

The principle of polarimetric data recording with ZIMPOL is different from that using IRDIS in DPI mode. In ZIMPOL, the incoming signal is modulated with a ferro-electric liquid crystal retarder and a polarization beam splitter downstream from the HWP, and sent to two demodulating CCD cameras of ZIMPOL. The demodulation technique provides quasi-simultaneous recording of intensities of opposite polarization states with the same detector pixels (Thalmann et al. 2008). The modulation/demodulation cycle frequency can be set to 967.5 Hz (fast polarimetry mode) or 26.97 Hz (slow polarimetry mode). A high detector gain of 10.5 e⁻ ADU⁻¹ in the fast polarimetry mode allows for the measurements of bright sources with short integration times without detector saturation. The slow polarimetry mode has a lower detector gain of 1.5 e⁻ ADU⁻¹ and a much lower read-out noise level providing a higher sensitivity for longer integration times compared to the fast polarimetry mode.

During three observing nights (see Table 1) we recorded 20 polarimetric cycles arranged in six blocks. For each cycle, four measurements of the Stokes parameters +Q, -Q, + U and -U were performed by rotation of the HWP by 0°, 45°, 22.5°, and 67.5°, respectively. All frames were taken without a coronagraph and with a DIT of 10 s. Such a long exposure saturated the science frames inside a circular region with a radius of ~10 pixels (~36mas).

In order to determine the position of the star in the saturated frames and correct for the beam shift in frames of opposite polarization states (Schmid et al. 2018), two additional unsaturated polarimetric cycles were recorded at the beginning and end of each block (on 2017 February 17 only at the block beginning) using a ND filter, ND2.0, with a transmission of 0.95% (see the SPHERE User Manual). These frames were also used for the measurement of the stellar flux to calibrate the photometry.

3 Data reduction

3.1 IRDIFS observations

The IRDIS and IFS data were calibrated with the SPHERE data reduction and handling pipeline esorex (Pavlov et al. 2008) and processed at the SPHERE Data Center (Delorme et al. 2017). The calibration of raw data consisted of background subtraction, bad pixel correction, flat fielding, correction of the pixel distortion (Maire et al. 2016), and extraction of the IFS spectral data cube.

The final calibrated data sets include two IRDIS temporal data cubes (K1 and K2 filters) and 39 IFS temporal data cubes (39 wavelength channels) with 48 frames each. The pixel scale of the IRDIS detector is 12.27 mas (K band), and a science frame is 1024 × 1024 pixels. In the IFS data the pixel scale is 7.46 mas and a science frame has a format of 290 × 290 pixels.

To subtract the stellar light, we used the SpeCal pipeline (Galicher et al. 2018), which provides several algorithms based on angular differential imaging (ADI; Marois et al. 2006) such as classical ADI (cADI) and template locally optimized combination of images (TLOCI; Marois et al. 2014), which were used to process both the IRDIS and IFS data sets. Figure 1 shows the cADI images of the spectrally combined IFS (left panel) and IRDIS (right panel) data.

Fig. 1 Total intensity images of the HD 114082 debris disk obtained with the cADI data reduction of the spectrally combined IFS data (left panel) and IRDIS data (rightpanel). The position of the star is marked with a yellow cross. Color bars show the flux values in counts per second (cts s−1).

Fig. 2 Polarized intensity images of the HD 114082 debris disk. The Qφ images are shown in the left panels and Uφ images in the right panels. IRDIS data taken in the H band are displayed in the top row. ZIMPOL data recorded with the I_PRIME filter are shown in the bottom row. In the IRDIS Qφ image, a strong residual signal at the coronagraph edge is masked out for clarity. The positions of the star are marked with yellow crosses. All images are smoothed via a Gaussian kernel with σkemei = 1 pixel and shown in the same spatial scale, indicated in the top-left image. The color bar shows the flux values in cts s−1.

3.2 IRDIS polarimetrie observations

For the data reduction of IRDIS polarimetrie data set we used the IRDIS Data reduction for Accurate Polarimetry (IRDAP) pipeline (van Holstein et al. 2020). This pipeline provides the basic data calibration steps such as sky-subtraction, flat-fielding and bad pixel masking. The processed frames are centered using center calibration frame with the satellite spots.

IRDAP then corrects for the instrumental polarization and polarization crosstalk using a Mueller matrix model describing the telescope and instrument. The pipeline also offers an additional correction for the linear polarization of stellar light if needed. We measured no significant polarization from the halo of starlight in the images of HD 114082.

The calibrated images of the Stokes parameters Q and U are finally obtained with a double difference method. They are used to calculate images of the azimuthal Stokes parameters Q_φ and U_φ presented in Fig. 2: $$Q_{\phi }=-Q\cos 2\phi -U\sin 2\phi$$(1) $$U_{\phi }=-Qsin2\phi -U\cos 2\phi ,$$(2)

where φ is the polar angle measured east of north in the coordinate system centered on the star and the sign convention for the Q_φ = −Q_r and U_φ = −U_r parameters defined in Schmid et al. (2006) is adopted.

The azimuthal Stokes parameter Q_φ image contains the linearly polarized intensity of the disk-scattered light as expected for an optically thin debris disk (e.g., Engler et al. 2017). In this work, we refer to the Stokes parameter Q_φ as simply polarized intensity because the component of circular polarization is much smaller in our observation and can be neglected.

3.3 ZIMPOL polarimetrie observations

The ZIMPOL polarimetric data set was reduced with the ZIMPOL data reduction pipeline developed at ETH Zurich. The pipeline includes preprocessing and calibration of the raw data: subtraction of the bias and dark frames, flat-fielding and correction for the modulation/demodulation efficiency. The instrumental polarization is corrected through the forced equalization of fluxes measured in the frames of perpendicular polarization states within an annulus of inner radius 110 pixels (400mas) and outer radius 200 pixels (730 mas).

To determine the shift between positions of the star in frames of different polarization states (Schmid et al. 2018) and position of the star in the combined intensity image, we fitted a 2D Gaussian function to the unsaturated point spread function (PSF) profiles in the frames taken with the ND filter. The position of the star on detector and the beam offset in the saturated science frames were interpolated as a function of the local sidereal time using the measurements with the ND filter recorded immediately before and after the respective block of science frames.

The format of the reduced images is 1024 × 1024 pixels with a pixel size of 3.6 × 3.6 mas. The spatial resolution of our data is about 25 mas.

For the analysis of ZIMPOL data we selected only the best quality frames where the star position could be well determined from the unsaturated frames: all cycles from 2017 March 14 and 15 and the first block of science data taken on 2017 February 17.

4 Disk morphology

4.1 Total intensity images

The IFS and IRDIS images of HD 114082 (Fig. 1) show a nearly symmetric debris ring with a radius larger than 0.3″ inclined at ~80° (~10° away from an edge-on configuration). Both images reveal a cavity inside the ring from ~0.3″ down to at least the inner working angle of the coronagraph at 0.11″.

The southeast side of the disk, close to the edge of the coronagraph, appears to be brighter than the southwest side. This could be a result of a higher particle number density on the southeast side or an offset of the ring with respect to the star toward the northwest. However, our data do not show any significant or measurable ring offset. A similar surface brightness (SB) asymmetry is observed in the image of the HD 117214 debris disk (Engler et al. 2020), which has no noticeable ring offset either.

A possible reason for such an asymmetric brightness, seen in images of both targets after data post-processing, could be instrumental effects, for instance the wind-driven halo, the low-wind effect (LWE), or stellar leakages outside the coronagraph (see Cantalloube et al. 2019, for a review of these effects). This explanation is also supported by the images of polarized intensity (Fig. 2, left panels), which show a rather symmetric east-west SB distribution. On the contrary, the north-south SB asymmetry of the debris ring clearly results from the asymmetric scattering of stellar light by dust particles, which are commonly assumed to scatter more light in the forward direction.

4.2 Polarized intensity images

Figure 2 displays the polarized intensity data taken in the H band with IRDIS (top row) and in the I_PRIME filter with ZIMPOL (bottom row). Both Q_φ images (left panels) reveal the same disk morphology as seen in the images of total intensity (Fig. 1). The back side of the disk is particularly visible in the ZIMPOL image.

The U_φ images (right panels in Fig. 2) contain only stellar residuals and, in the IRDIS image, a faint residual disk signal that comes from the effect of convolution of the PSF with the disk image (see the appendix in Engler et al. 2018, for description of the effect). This spurious U_φ signal is not seen in the ZIMPOL image (right lower panel in Figs. 2 and A.2) because the spatial resolution of ZIMPOL is ~3 times higher than the resolution provided by IRDIS.

To support this statement, we modeled the images of the Stokes parameters Q and U to calculate the parameter U_φ for the HD 114082 disk. Figure A.1 shows the modeled U_φ image resulting from the convolution of parameters Q and U with the IRDIS PSF. This image reproduces the pattern of the IRDIS U_φ image (right upper panel in Fig. 2) well.

4.3 Model for scattered light

To estimate the geometrical parameters of the disk we used a 3D model for single scattering of stellar photons by dust particles. The detailed description of models for scattered and polarized intensity images can be found in Appendix B. All model parameters are listed in the first column of Table 2.

In order to find a set of parameters that best fits the total intensity image in the K1 band (Fig. 3a), we ran a custom Markov chain Monte Carlo (MCMC) code. At each iteration in this code, a new model with a parameter set drawn from the uniform prior distributions (Col. 2 in Table 2) is created, convolved with the instrument PSF and inserted into an empty data cube at different position angles (PAs) to mimic the rotation of the sky field during the observation. A cADI forward-modeling is then performed and the final model image is compared to the imaging data.

We evaluated two families of models with two different SPFs: the first SPF with a single Henyey–Greenstein (HG) function (hereafter we refer to this model as model 1), and the second SPF with a two-component HG function referred to as model 2. The detailed description of the HG functions and the applied MCMC fitting method is given in Appendix B.

The posterior distributions of fitted parameters for both families of models are shown in Figs. E.1 and E.2. We consider the median values of these distributions as the best-fit parameters and list them in Cols. 3 and 4 of Table 2. Both models yield similar best-fit parameter sets with values for the radius of the planetesimal belt R_belt of 0.33 ± 0.03 (model 1) and 0.35 ± 0.03 arcsec (model 2), and disk inclinations of 82-83°. The index of the inner radial power law for the particle number density α_in is the only parameter that remains unconstrained. The MCMC sample for this parameter shows a clear tendency to values larger than 40 indicating a sharp inner edge of the disk.

Wahhaj et al. (2016) fitted a different model for the spatial distribution of dust particles to the scattered light images of the HD 114082 disk obtained in the H band with SPHERE. The authors modeled the disk surface density within the range of radial distances r_in and r_in + Δr, with r_in being the inner radius and Δr the width of the disk, using three images from various data reductions: principal component analysis (PCA), matched locally optimized combination of images (MLOCI; Wahhaj et al. 2015), ADI, and MLOCI+reference-star differential imaging (RDI). Their results related to the radial location of the belt are comparable with the radius for the planetesimal belt R_belt we derived for our models 1 and 2: r_in + Δr = 0.31″ + 0.05″ for the PCA image. However, the disk inclination of 84.9° in the PCA model is slightly outside of the 68% credible intervals for the inclination provided by our models. In contrast, the other fitting procedures (MLOCI+ADI and MLOCI+RDI) preferred also a lower inclination of the disk down to 81°. As a SPF, Wahhaj et al. (2016) used a single HG function with one asymmetry parameter g. Their best-fit values for g are between 0.07 and 0.23, which are much lower than our result (g = 0.56 in model 1). These differences can be explained by a lower quality of the images due to unfavorable observing conditions² and more aggressive data reduction techniques.

With the best-fit parameters (Cols. 3 and 4 of Table 2) we created two images for a visual and statistical comparison of model 1 and model 2 to the data. Figure 3 displays forward-model images of both models (Figs. 3b and d). It shows also the residual images (Figs. 3c and e) obtained after subtraction of forward-model images from the data (Fig. 3a). The visual comparison of the residual images suggests that model 2 fits the data slightly better.

The χ² analysis, which we performed in order to assess the goodness-of-fit of both models, also supports this conclusion. We calculated the reduced ${\chi }_v^2$ metric according to $${\chi }_v^2=\frac{1}{v}{\displaystyle \sum _{i=1}^{N_{\text{data}}}\frac{{\left[F_{i\text{,data}}-F_{i,\text{model}}\left(p\right)\right]}^2}{{\sigma }_{i\text{,data}}^2}},$$(3)

where N_data is a number of modeled resolution elements, F_i,data is the flux measured in element i with an uncertainty σ_i,data and F_i,model(p) (Eq. (B.6)) is the modeled flux of the i element. The degree of freedom of the fit is denoted by v = N_data − N_par(ℳ), where N_par(ℳ) represents the number of variable parameters $p=\left(p_1,p_2,\dots ,p_{N_{\text{par}}}\right)$ of model ℳ. In total, model 1 has nine parameters, and model 2 eleven (Col. 1 of Table 2). The full width at half maximum of the stellar profile in the IRDIS flux images is ~3 pixels. Therefore, we consider an image area of 3 × 3 binned pixels as one resolution element. The uncertainties σ_i,data are computed as the standard deviation of the flux distribution in concentric annuli of the total intensity image excluding regions with the disk flux.

We obtained a reduced ${\chi }_v^2=1.79$ for model 1 (v = 458) and ${\chi }_v^2=1.67$ for model 2 (v = 456). Although both ${\chi }_v^2$ values are relatively close to each other, the Bayesian information criterion (BIC) indicates that model 2 provides a better fitting result. We calculated the BIC according to Schwarz (1978): $$\text{BIC}\left(ℳ\right)={\chi }_v^{2}v+N_{\text{par}}\left(ℳ\right)\ln \left(N_{\text{data}}\right)+\text{const.}$$(4)

For models 1 and 2, the difference BIC(model 1) – BIC(model 2) is positive and equal to 46. The lower BIC value of model 2 implies that this model should be preferred over model 1 even if the latter has fewer fitting parameters. We have to note, however, that a correlation between the adjacent resolution elements is not taken into account in Eqs. (3) and (4), and, thus, the actual BIC difference value could deviate from the adopted one.

Fig. 3 Comparison of the total intensity image (panel a) with the forward models of scattered light with a single HG function (panel b) and a two-component HG function (panel d). Panels c and e: residual images obtained after subtracting the forward model images from the total intensity image. The color bar shows the flux values in cts s−1.

5 Phase functions and polarization fraction

The SPF and the polarized SPF (pSPF) of a debris disk are respectively defined as a function that describes how the total and polarized flux scattered by a population of dust particles is spatially distributed as a function of scattering angle θ (0 < θ < 180°). The polarized flux represents only a fraction of the total scattered flux. This fraction depends on the scattering angle as well and can be defined as a PF phase function p(θ): $$p\left(\theta \right)=\frac{p\text{SPF}\left(\theta \right)}{\text{SPF}\left(\theta \right)}.$$(5)

In this work, when we refer to the PF, we always mean a degree of linear polarization of scattered light because the circular polarization is negligibly small in the observations of debris disks and, therefore, not taken into account.

All three phase functions together with the value of maximum PF of scattered light p_max are unique for each debris disk and, as a tool, offer an opportunity to study the reflectivity of small dust particles in the disks surrounding distant stars. Light scattering theories, like the Mie theory, predict, and laboratory measurements (e.g., Muñoz et al. 2017; Frattin et al. 2019) show that the phase functions contain an information about dust optical properties. Therefore, many previous debris disk studies at optical and near-IR wavelengths used the phase function analysis to infer the characteristics of observed dust populations, for instance particle effective size or composition (e.g., Graham et al. 2007; Milli et al. 2017b, 2019; Esposito et al. 2018; Ren et al. 2019; Engler et al. 2020).

For this kind of analysis, the phase functions can be measured from the images of debris disks in scattered light. However, this is a challenging task that requires data with a high signal-to-noise ratio (S/N), additional data calibrations and corrections for different effects such as disk flux losses due to the data post-processing technique or convolution with the instrument PSF. The correction factors for these effects can be estimated only by modeling the disk geometry and forward-modeling of the scattered (polarized) light image.

Furthermore, the calculation of the PF requires two images in the same wavelength range: a total intensity image and a polarized intensity image of the disk. HD 114082 was observed in field-stabilized polarimetric mode with both IRDIS and ZIMPOL. Therefore, two images showing the disk polarized intensity are available (in the H band and in the I_PRIME filter, Fig. 2) but there are no suitable images of the disk in total intensity from the same observations. The IRDIS broadband H image of the disk from the previous SHARDDS observation suffers a lot from the LWE in the regions of interest, as mentioned in Sect. 4.3. Therefore, in order to measure the SPF and PF of the HD 114082 disk at near-IR wavelengths, we created a combined total intensity image corresponding to the H band by summing up the total intensity images in the IFS channels from 1.45 to 1.64 μm. This wavelength range covers approximately the range of the broadband H IRDIS filter used in the DPI observations of HD 114082. Thus, the obtained IFS image for the total intensity is directly comparable to the polarized intensity image in the broadband H filter shown in Fig. 2 (left panel).

5.1 Scattering phase function

To derive the SPF of the HD 114082 disk shown in Fig. 4, we measured the total intensity of scattered light within circular apertures with a radius of 19 mas (0.5 λ/D; i.e., the half width at half maximum of the PSF) or an area of 1.13 × 10⁻³ arcsec² placed along the disk rim in the combined IFS image. Before applying aperture photometry, we corrected the IFS ADI data for the disk flux losses resulting from the ADI post-processing. For this correction, we multiplied the IFS image by a scaling factor map that was computed as a flux ratio in each pixel between the image of the best-fit model convolved with the instrument PSF and the image obtained after the ADI post-processing of the same model. The post-processing was applied to a data cube consisting of 80 frames, where the model image (convolved with the instrument PSF) was injected at different PAs corresponding to the parallactic angle change during observation.

The flux values shown in Fig. 4 are normalized to a sum of stellar fluxes in the IFS channels corresponding to the H band (wavelengths 1.45–1.64 μm). Because of the high inclination and small angular size of the disk, we were able to perform aperture photometry only within the range of scattering angles 25° < θ < 155°.

By using this measurement method, we assume that the disk is thin in the vertical direction so that the derived aperture flux can be associated with a single scattering angle. The vertical thickness of the disk, however, affects the shape of the measured SPF at scattering angles close to 90° for highly inclined disks (Olofsson et al. 2020). This effect, sometimes referred to as limb brightening, due to the increased column density of dust particles around the disk major axis (θ = 90°) is also observed in our measurements: the SPF curve (black solid line in Fig. 4) shows a slight bump between θ = 60° and θ = 120° compared to the smoothly decreasing model SPFs. As seen in Fig. 4, the SPF of model 2 seems to match the measurement better.

We compare the measured SPF of HD 114082 disk in Fig. 5 with the empirically derived phase functions of other debris dust systems as compiled in Engler et al. (2020). This figure includes the SPFs of several debris disks: HD 117214 (Engler et al. 2020), HR 4796 A (Milli et al. 2017b), HD 61005 (Olofsson et al. 2016), HD 191089 (Ren et al. 2019), HD 35841 (Esposito et al. 2018) and HD 15115 (Engler et al. 2019). The phase functions obtained for the zodiacal dust in the Solar System (Leinert et al. 1976), ring D68 of Saturn (Hedman & Stark 2015) and comet 67P/Churyumov–Gerasimenko (data set MTP020 taken on 2015 August 28; Bertini et al. 2017) are also plotted in this figure, as well as the theoretical SPFs derived using the Fraunhofer diffraction model for the particles with radii of 1, 5, 15, and 50 μm (Min et al. 2010; Engler et al. 2020) and Hapke reflectance theory (Hapke 1981) for the dust grains covered by small regolith particles (microasteroid SPFs in Fig. 5; Min et al. 2010). For easier comparison, all displayed SPFs are normalized to their values at scattering angle of 90°.

Between θ = 40° and θ = 90° the SPF of the HD 114082 disk has a negative slope like the SPFs of other debris disks except for HR 4796 A. For scattering angles larger than 90°, the HD 114082 SPF seems to decrease in a similar way as the SPFs of the solar zodiacal dust and Saturn’s ring. In the whole range of measured scattering angles (25° < θ < 155°), the trend of the HD 114082 phase function is most similar to that of the debris disk surrounding HD 191089. HD 191089 is a young F star (F5V) in the β Pictoris moving group that harbors an exo-Kuiper belt with a radius of ~45 au (Ren et al. 2019).

As seen in Fig. 5, the SPF curve of the HD 114082 disk rises for scattering angles smaller than 90° and is contained between the microasteroid curves of 5 and 15 μm. If the forward-scattering peak of phase function is governed by Fraunhofer diffraction (see Sect. 6 in Engler et al. 2020), the average effective size of particles probed with this near-IR observation (H band) is expected to be between 5 and 15 μm. The size of 5 μm would be consistent with the blow-out size of 1–3 μm for the dust particles, which are mainly removed from the system by the radiation pressure of the host star. Adopting stellar luminosity L_* = 3.86 L_⊙ and mass M_* = 1.4 M_⊙ (Pecaut et al. 2012), and the average bulk density of dust particles ρ of 2–3 g cm⁻³, we derive this blow-out size from relation between radiation pressure force and gravitational force acting on dust particles (Burns et al. 1979): $$a_{\min }=\frac{3L_{\ast }{Q}_{\text{pr}}}{16\pi G{M}_{\ast }c\rho },$$(6)

where Q_pr is the radiation pressure coupling coefficient, G is the gravitational constant and c is the speed of light. For a F star, most of the stellar energy is concentrated at wavelengths near 0.5 μm. Therefore, for dust particles with sizes larger than 1 μm, the size parameter x = 2πα/λ ≫ 1 and the geometrical optics approximation holds. This means that the effective cross section of particles can be approximated by their geometric cross section resulting in Q_pr ≈ 1, and thus the particle minimum size a_min is around 3 μm for HD 114082.

If we assume that dust is produced through a collisional cascade, then its size distribution should follow a power-law of index q = −3.5, for which most of the geometrical cross section should be in the smallest (bound) particles with sizes a > 3 μm, so that the scattered flux should be dominated by these particles, at least at wavelengths smaller than 2π a_min (e.g., Thébault & Augereau 2007). Thus, we can expect to trace dust particle population with the effective particle size of ~5 μm in this observation.

Fig. 4 Scattering phase function of the HD 114082 debris disk normalized by the stellar flux. The gray-shaded area indicates the uncertainty on the mean flux (black line) of the east and west disk sides.

5.2 Polarized scattering phase function

To retrieve the pSPFs (dashed lines in Fig. 6), we applied the same aperture photometry to the Q_φ images in the H band and I_PRIME filter (Fig. 2, left panels): the circular apertures with a radius of 19 mas were placed at the same positions in the disk as in the IFS H band image of the total intensity. The measured aperture fluxes F_aper plotted in Fig. 6 are normalized by the stellar flux F_star in the corresponding filter and therefore represent the disk contrast in polarized light at measured position.

Both retrieved pSPFs have a similar shape and exhibit a bow between θ = 80° and θ = 110°. The highest disk contrast is measured at the disk ansae (around θ = 90°) and at the smallest scattering angle θ = 25° we are able to measure. This contrast reaches (4.0 ± 1.0) × 10⁻³ arcsec⁻² in the H band and (3.9 ± 0.8) × 10⁻³ arcsec⁻² in the I_PRIME filter. A higher polarized contrast at the disk ansae can be explained by a higher column density of dust particles along the line of sight (LOS) at the position of the planetesimal belt (limb brightening effect) and by the highest PF of scattered light at angles close to 90°. At smaller scattering angles, the increasing polarized flux results from the rapidly growing amount of scattered light due to the forward-scattering behavior of dust particles. For θ > 100° the polarized flux is continuously decreasing and apparently diminishing toward zero at θ > 140°. The shapes of both measured pSPFs thus are quite similar to that of HR 4796 A (Milli et al. 2019).

The impact of limb brightening on the shape of measured pSPFs can be estimated with the model of a disk that scatters light isotropically. The result of such an estimation is model-dependent but, if corrected, the measured phase function should reproduce the shape of the actual phase function better. Therefore, we applied the aperture photometry to the image of model 1 generated with the HG asymmetry parameter g = 0 (isotropic scattering) to calculate the correction factor for the pSPFs in dependence on scattering angle θ. The correction factor was derived as a ratio between the flux value measured within a circular aperture at a given position and the smallest flux value we obtained for θ = 25°. Multiplying the measured pSPFs by this factor, we derived the corrected phase functions, which are flat between θ = 60° and θ = 110° (solid lines in Fig. 6).

Fig. 5 Scattering phase functions measured for different debris disks and dust populations in the Solar System in comparison with the theoretical microasteroid phase functions and diffraction model for micron-sized particles.

Fig. 6 Polarized scattering phase functions of the HD 114082 debris disk in the H band and I_PRIME filter normalized by the stellar fluxes. The pSPFs corrected for the limb brightening effect are shown with the solid lines. The light red-shaded area indicates the uncertainty on the measured disk contrast in the H band. The light blue-shaded area indicates the uncertainty on the measured disk contrast in the I_PRIME filter.

5.3 Degree of linear polarization of scattered light

5.3.1 Polarization fraction in the H band

To measure the PF of the scattered flux in the H band, we used the aperture photometry results presented in Sects. 5.1 and 5.2. As mentioned above, we preferred to use the IFS H band image for this analysis instead of the archival IRDIS H broadband image of total intensity. Nevertheless, we utilized the latter to put constraints on the PF at the disk ansae, where the disk flux is least disturbed by the LWE.

The PF p(θ) in Fig. 7 is derived by dividing polarized intensity by the total intensity, corrected for the disk flux losses caused by the data post-processing (see Sect. 5.1), which are both normalized to the median value of the stellar flux measured in the IRDIS broadband H (F_* (IRDIS)) and in the IFS H band image (F_*(IFS)): $$p\left(\theta \right)=\frac{F_{{\text{aper,Q}}_{\varphi }}\left(\theta \right)}{F_{\text{aper,IFS}}\left(\theta \right)}\cdot \frac{F_{\star }\left(\text{IFS}\right)}{F_{\star }\left(\text{IRDIS}\right)}\cdot$$(7)

The purpose of this flux normalization is to take into account the disk flux variations due to different observing techniques with two different instruments as well as observing conditions and slightly different wavelength ranges, which we compare with each other. All these factors contribute to the error of the PF measurement but their effects are essentially canceled out when applying the flux normalization.

Figure 7 shows the PF derived from the aperture photometry for both disk sides and their mean value (black line) that represents the measured PF before the corrections for the dilution of the polarized light. Within the range of investigated scattering angles (25° < θ < 155°), the PF has a bell shape with a maximum value of p_max ≈ 11.5% between θ = 70° and θ = 115° and a minimum value of p_min < 2.5% beyond θ = 140°.

The PF derived from archival H broadband image of total intensity reaches similar values at disk ansae (θ = 90°) although with a larger uncertainty. For the west ansa of the disk we obtain p_max = (11.2 ± 2.5)% and for the east ansa p_max = (10.8 ± 2.5)%, which are in good agreement with the PF estimated using the IFS data.

The gray-shaded area in Fig. 7 shows the measurement uncertainty that takes into account only statistical noise of the flux measurement when applying aperture photometry to the Q_φ and IFS H band images. This uncertainty does not include the uncertainty resulting from the model-dependent correction of the disk flux losses caused by the ADI data post-processing and the methodological errors mentioned above, which are difficult to quantify. In particular, the observing conditions, for example seeing, coherence time, or atmospheric transparency, can strongly decrease disk flux and have a significant impact on observational results. For instance, if the coherence times are short (τ < 2−3 ms), the position of the star behind the coronagraph is unstable. In this case, a very accurate registering of frames is not possible, since the position of the star is not precisely known and deviates from the stellar position in the center frames. For a disk with a small angular size, such as the HD 114082 disk, the location of the maximum value of the PF (scattering angle at which the PF maximum occurs) depends on position of the star in the image. Thus, the frame centering should be performed as accurately as possible.

At this point, we had only considered the disk flux losses in the IFS total intensity data. However, the disk flux is reduced not only in the total intensity image but also in the polarized intensity image. The polarized flux we measure in the Q_φ image, and hence the estimated PF shown by black line in Fig. 7, is lower compared to the actual disk polarized flux. The reason for this flux reduction is the extended instrument PSF, which leads to the polarization cancelation effect in the Stokes Q and U signals (Schmid et al. 2006).

There is another important effect that is responsible for the underestimation of the polarized flux amount in the range of observed scattering angles. This phenomenon is not related to the instrument but to the disk morphology itself. When performing the aperture photometry, we often measure a smaller PF than the intrinsic PF of scattered light in the disk midplane as a result of polarized flux dilution along the LOS. We call such an attenuation of polarized flux the LOS effect in this work and discuss it in more details in Sect. 5.3.3.

In order to derive the “true” PF of scattered light, the measurement needs to be corrected for both effects by estimating their magnitudes with forward-modeling (e.g., Engler et al. 2018). For this purpose, we created several images of disk polarized light using the best-fit parameters of both models and iterating the shape of PF phase function each time.

For the first iteration, the measurement from the data was used as a first guess for the PF. To model the image we adopted a pSPF given by a product of the model SPF and a polynomial fit (solid green line in Fig. 8) to the measured mean PF of the disk (dotted black line in Fig. 8).

The forward-model Q_φ image contains less polarized flux compared to the initial guess because the LOS and PSF convolution effects reduce the polarized flux as discussed above. Therefore, the PF that we measure from the forward-model Q_φ image is lower than that given by the initial PF. As an example, the dashed blue line in Fig. 8 shows the PF measured from the forward-model images of total and polarized intensities of model 2³. As expected, it is lower than the polynomial fit used as the initial guess for the PF (solid green line). It is also flatter.

Therefore, for the second iteration, we compute a new PF for our model multiplying the measured PF (dotted black line) by a ratio between the first polynomial fit (solid green line) and the PF measured from the forward model after first iteration (dashed blue line). The dash-dotted blue line in Fig. 8 shows a polynomial fit to the data corrected by this ratio. This fit is used as the PF in the second iteration.

The PF measured from the forward-model images of second iteration (solid blue line in Fig. 8) matches the data measurement very well. At smaller (θ = 30°) and larger (θ = 130°) scattering angles there is no large discrepancy in values between the PF used to create model images (dash-dotted blue line) and the PF measured from the processed model. However, the actual PF maximum value p_max used to model the data is much higher than 11.5% measured from the data and reaches 17.4% in model 2 (16.6% in model 1). We conclude that the actual maximum PF should be about 17% for the dust particles in HD 114082 disk observed at near-IR wavelengths, and estimate the modeling uncertainty on this value to be 3–4% including modeling uncertainty from the ADI post-processing effect.

Fig. 7 Polarization fraction of the debris disk HD 114082 derived as a mean value of the measurements on the west and east disk sides in the H band (black line). The gray-shaded area indicates the uncertainty on this mean PF from the aperture photometry. The blue arrows show a lower limit for the PF in the I_PRIME filter.

Fig. 8 Polarization fraction of scattered light in debris disk HD 114082. The dotted black line shows a mean value of the PF derived for the west and east disk sides. The green line shows a polynomial fit to it.

Fig. 9 Comparison of the Qφ image (panel a) to the model images of polarized light obtained with parameters of model 1 (panel b) and model 2 (panel c). The color bar shows flux in cts s−1.

5.3.2 Model validation

For a visual examination of the goodness of fit, we compare the images of models 1 and 2 showing the disk in scattered polarized light with the Q_φ image in the H band in Fig. 9. The model images shown in Figs. 9b and 9c are created using the pSPF given by the product of the SPFs of both models and the intrinsic PF (blue dash-dotted line in Fig. 8) found with the forward-modeling approach in Sect. 5.3.1.

Both models provide similar images that fit the Q_φ image well, except for the regions very close to the coronagraph edge and where the coronagraph throughput is higher than 50%. At these radial separations (0.10″ < r < 0.15″) model 2 seems to match the SB distribution of the debris ring better because the SPF of model 2 has a much higher forward-scattering peak (g₁ = 0.83) compared to model 1 (g = 0.56). The higher scattered flux means that the polarized scattered flux is also higher in model 2 at smaller scattering angles because the same PF was adopted in both models. The stronger forward-scattering becomes apparent through the higher polarized SB close to the coronagraph edge in the image of model 2 (Fig. 9c). Similar SB peaks are observable in the Q_φ image (Fig. 9a), which would support the preference for model 2 over model 1. However, we should note that this region in the Q_φ image is relatively noisy because of the coronagraph and telescope spider effects.

5.3.3 Line-of-sight effect

Our modeling results confirm that not only does the effect of the convolution with the instrument PSF significantly influence the outcome of the PF measurement, but the disk morphology also has a large impact.

As mentioned in Sect. 5.1, we attributed the disk flux measured within the circular apertures to a single scattering angle corresponding to the dust location in the midplane of a narrow planetesimal belt. By doing so, we neglected the effect of disk thickness on the amount of polarized flux measured along the LOS direction. However, the nonzero disk extension in both the radial and vertical directions means that scattering from many locations with different scattering angles contributes to the disk flux that we measure within the circular apertures in the images of both total and polarized intensities. Through the integration along the LOS, in particular the polarized light can become diluted because of contributions from the scattering angles with lower degree of polarization compared to the disk midplane. This leads to the underestimation of the disk PF, for instance the maximum PF value at θ = 90°.

We can probe this LOS effect by measuring the PF from the non-convolved model images (total and polarized intensities) using the same aperture photometry method. This measurement is indicated by a dotted blue line in Fig. 8 (for model 2). It shows the impact of disk morphology on the measured PF. Even without considering the effect of the instrument PSF, we derive a lower degree of disk polarization and a lower PF maximum value of 13.9% compared to that used to create the model images.

Not only is the magnitude of the measured maximum PF p_max influenced by the LOS effect but its position is also impacted. As seen in Fig. 8, the intrinsic PF of the model has a maximum at θ = 85° (dash-dotted blue line) whereas the measured PF maximum is shifted to θ = 95° for both convolved and non-convolved models (blue solid and dotted line, respectively).

In Appendix B, we discuss how the inclination and scale height of the disk affect the PF measurement. Depending on these two parameters the measured maximum PF p_max can be reduced by 30% compared to the intrinsic maximum PF of the disk. Since the inclination of the HD 114082 disk is constrained, we estimated the impact of the scale height on the measured p_max. We used model 2 to create images of the disk with an inclination of 83° and different scale heights: H₀ = 0.004″ (best-fit parameter of model 1), 0.007″ (best-fit parameter of model 2), 0.01″ and 0.1″. Our results show that the scale height has a rather small effect in this particular case. The measured maximum PF p_max is 13.9% for models with the scale heights of 0.007″, 0.01″ and 0.1″. It is 14.2% in case the scale height equals to 0.004″.

5.3.4 Polarization fraction along the disk major axis

Using the same aperture photometry, we measured the PF along the disk major axis. The derived PF mean value between the east and west disk sides is shown by the dotted black line in Fig. 10. It continuously increases when approaching the radial position of the planetesimal belt, where the peak SB is located (r ≈ 0.36″). The highest mean PF value of ~11.2% is measured between 0.35″ and 0.37″, thus at the location of the belt.

In our ideal model, shown in Fig. 10, the PF is constant and equal to p_max independently of the separation, because all dust particles located in the disk along its major axis scatter stellar radiation at a scattering angle of 90° toward the observer. In our best-fit model this scattering angle corresponds to a maximum PF of 17.4% (see the dash-dotted blue line in Fig. 8). This value does not depend on the distance from the star because we assume an uniform population of dust particles throughout the disk in the model. This assumption is certainly valid for the range of particle sizes that we probe in the observations at near-IR wavelengths.

However, when applying this measurement to the image of the non-convolved best-fit model, we obtain the maximum PF value p_max of 14.2% for r > 0.37″ (dotted blue line in Fig. 10). As expected, the PF measured from the model image is lower than the actual model value of 17.4% due to the LOS effect discussed in Sect. 5.3.3. The PF also varies with the radial distance (0.14″ < r < 0.37″) because the aperture encloses dust grains emitting at different scattering angles, in particular decreasing the PF at smaller separations.

The effect of convolution with the instrument PSF additionally decreases the observed PF. To demonstrate this, we perform the measurement with the convolved best-fit model and show the result with the solid blue line in Fig. 10. The PF measured from the convolved model image matches the data very well reaching ~11.5% at a radial separation of 0.35″. Beyond r = 0.35″, it slowly increases up to 12.6% and remains at this level following the PF trend derived from the non-convolved model image. This contrasts with decreasing PF measured from the data. The apparent discrepancy between the model and the data could be explained by the low polarized flux measured in the Q_φ image at larger distances.

In general, we expect the PF of the scattered light measured along the disk major axis to be slowly growing beyond the planetesimal belt because the effects of LOS and PSF convolution are getting smaller with the distance to the star. In our Q_φ image, the polarized flux at larger separations might be diminished through the applied correction for the instrumental polarization. Though, the other reasons, for instance different optical properties of dust particles at larger distances because of different shape or composition, cannot be completely ruled out.

Fig. 10 Degree of linear polarization of scattered light for the debris disk HD 114082 calculated along the disk major axis. The black line shows a mean value of the west and east disk sides. The radial dependence of the PF is significantly moderated by the difference in scattering angles, such as in Fig. 8 due to the LOS and aperture size effects.

5.3.5 Constraints on the polarization fraction from ZIMPOL data

As mentioned above, only the polarized intensity image of the disk is available in the ZIMPOL I_PRIME filter. The ZIMPOL data were taken in the field-stabilized mode, which does not allow the disk to be detected in total intensity because the region of interest is contaminated with the stellar PSF. For this reason, it is not possible to derive the PF of scattered light using the ZIMPOL data set alone.

However, we can put constraints on PF in the I_PRIME filter comparing the magnitude of the polarized flux measured in the ZIMPOL Q_φ image (Fig. 2, bottom left panel) with the total intensity of the disk measured at shortest IFS wavelengths. For this purpose, we applied the aperture photometry to the Q_φ image in the I_PRIME filter (λ = 0.790 μm, Δλ = 0.153 μm) and the total intensity IFS image combining wavelength channels from 0.95 to 0.99 μm as described in Sect. 5.3.1. The measured fluxes were also normalized according to Eq. (7).

The PF estimated in this way is indicated by the blue arrows in Fig. 7. It is slightly lower than the PF obtained for the H band and represents a lower limit on the PF in the I_PRIME filter because the corrections for the attenuation of polarized flux (LOS effect and convolution with the instrument PSF effect) are not applied to this measurement.

5.3.6 Comparison with other debris disks and dust populations

Based on modeling results (Sect. 5.3.1), we conclude that the PF maximum value in the H band is 17 ± 4% at the radial position of the HD 114082 planetesimal belt.

Similar PFs were measured for a number of other exo-Kuiper belts, although the results of previous studies are most probably underestimated because they do not take into account the convolution of polarized signals with the instrument PSF and LOS effects decreasing the measured PF. Using the High-Contrast Coronographic Imager for Adaptive Optics (HiCIAO) at Subaru Telescope, Tamura et al. (2006) detected ~10% polarization in the K band (2.2 μm) in the β Pictoris debris disk at radial separations 50 < r < 120 au: at distances, where the edges of one or two planetesimal belts are expected to be (Heap et al. 2000; Golimowski et al. 2006). The optical (BVRI) imaging polarimetry of the β Pictoris disk showed a higher level of polarization of 15–20% but at larger radial distances (r > 150 au; Gledhill et al. 1991; Wolstencroft et al. 1995). Asensio-Torres et al. (2016) also used the HiCIAO instrument to measure ~15% of polarization for the outer belt in the HD 32297 disk. Bhowmik et al. (2019) obtained a polarized flux fraction of ~15% in the J band (1.2 μm, SPHERE/IRDIS) for the HD 32297 disk as well. A recent study of another double-belt system HD 15115 reported a maximum PF of 10–15% at radial location of the outer belt measured in the J band with SPHERE/IRDIS (Engler et al. 2019).

The shape of the measured HD 114082 PF is also very similar to that derived for the dust in debris belt HR 4796 A (Milli et al. 2019). It is interesting to note that the HD 114082 PF seems to approach a zero value at scattering angle close to 160° like the HR 4796 A phase function does. Many studied samples of cometary dust in the Solar System (Hadamcik et al. 2016; Shestopalov & Golubeva 2017) and their analogs (Frattin et al. 2019) behave in the same way: they exhibit a negative polarization branch for scattering angles larger than 160°.

6 Disk near-IR spectrum

The IFS and IRDIS data allow us to investigate the scattering efficiency of dust particles in the HD 114082 disk at near-IR wavelengths. We perform this analysis comparing the dust reflectivity in two debris disks: HD 114082 and HD 117214. This comparison was motivated by the idea to show that the HD 114082 disk is much brighter than a debris disk with a similar IR excess, temperature of dust grains and host star characteristics.

The HD 117214 is a young F6V star (Houk & Cowley 1975) located at the distance of 107.4 ± 0.3 parsecs (Gaia Collaboration 2021) in the LCC subgroup as well. It possesses a bright debris disk (Engler et al. 2020) with a radius of ~0.45″ (48 au) and fractional luminosity of 2.53 × 10⁻³ (Jang-Condell et al. 2015), which are similar to the parameters of the HD 114082 disk. HD 117214 is thus far the only F star to have been observed under very good observing conditions with SPHERE using the same instrument configuration (Engler et al. 2020), making such a comparison possible. Since both targets are young F stars in the LCC subgroup, their circumstellar dust might potentially exhibit similar scattering properties.

6.1 Spectral analysis

In order to derive the spectra of HD 114082 and HD 117214 debris disks, we applied the aperture photometry to the IFS and IRDIS images as described in Sect. 5.1. The disk contrast presented in Fig. 11 was measured along the disk major axes at positions of the peak SB in the r²-scaled disk images, which indicate the location of the planetesimal belt. These positions were selected because the disk flux is less affected (reduced) by the ADI data processing at larger radial distances. We placed the apertures at radial separations of 0.37″ (35 au) from the star in the image of HD 114082 disk, and at radial separations of 0.46″ (49 au) in the HD 117214 disk image.

The displayed spectra (Fig. 11) are corrected for the disk flux losses (see Sect. 5.1) and normalized by the value of stellar flux obtained from the flux measurement frames in the corresponding IFS wavelength channel or IRDIS band. The spectra show thus the disk SB in terms of contrast to the host star. This contrast was measured in the circular apertures with an area of 4.25 × 10⁻³ arcsec². The normalization by the stellar flux allows us to remove disk flux fluctuations that reflect the variation of the stellar spectrum, among other reasons due to spectrally dependent transmission of both instruments and the Earth’s atmosphere. As practice shows, especially the correction for the telluric water absorption bands close to 1.4 μm is challenging, because the absorption of water molecules is very strong in this wavelength range. The error bars in Fig. 11 indicate the uncertainty of the aperture photometry and do not include uncertainties from the modeling of flux losses or any calibration issues.

Except for the specific regions of telluric water vapor absorption (see Fig. 11), the derived reflectance spectra should represent the underlying disk spectra well, in particular because the spectrum of HD 117214 disk shows some dips at locations where the stellar spectrum is featureless, for instance at 1.2 μm and in the region between 1.5 and 1.6 μm. In these spectral ranges, the CO₂ ice exhibits relatively strong absorption bands: at 1.21,1.44, and 1.54–1.58 μm (Hansen 2005). Therefore, the measured spectrum of HD 117214 disk maybe indicates the presence of CO₂ ice as a constituent of debris material at the investigated radial distance from the star.

In fact, it would not be surprising to identify this ice in a young debris system. The CO₂ is one of the most abundant ices after H₂O ice. It was found in many lines of sight toward molecular clouds and young stellar objects (Öberg et al. 2007). The CO₂ gas is also present in the Earth atmosphere but in a very small amount. The absorption features, which we possibly observe in the HD 117214 spectrum, should have another origin. Otherwise, we would clearly see them in the stellar spectrum as well.

The CO₂ gas freezes out over a narrow range of dust temperatures between 40 and 80 K (Harsono et al. 2015). Under the assumption of local temperature equilibrium, the dust temperature at location of planetesimal belt can be estimated using the canonical expression from Backman & Paresce (1993): $$T_{\text{\hspace{0.17em}dust}}=\left(278\text{\hspace{0.17em}K}\right){\left(\frac{L_{*}}{L_{\odot }}\right)}^{0.25}{\left(\frac{1}{R_{\text{belt}}}\right)}^{0.5},$$(8)

where the L_* is the stellar luminosity, and R_belt is the radius of the planetesimal belt in au. The luminosity of F6V star HD 117214 is found to be 5.64 L_⊙ (Jang-Condell et al. 2015) whereas the luminosity of F3V star HD 114082 is 3.86 L_⊙ (Pecaut et al. 2012). The estimated dust temperature in the HD 117214 disk at the distance of 49 au from the star (disk ansae) is 61 K. For the HD 114082 disk we obtain a dust temperature of 66 K at the distance of planetesimal belt at 35 au. Thus, both debris belts could be located at the CO₂ snow line.

The spectrum of HD 114082 disk shows no obvious absorption features except at 1.4 μm, which most likely comes from the insufficient telluric correction of spectrum. Nevertheless, the presence of CO₂ cannot be ruled out for HD 114082.

In the wavelength range between 0.96 and 1.66 μm both disk reflectance spectra are rather flat or have a gray color at measured radial position in the disk. However, they seem to show a red color at longer wavelengths than 1.66 μm because the contrasts in K1 and K2 bands are higher than the contrasts measured in the IFS channels with the K2 band contrast being the highest. In this respect, the HD 114082 and HD 117214 disks are similar to the debris disks HR 4796 A (Milli et al. 2017b) and HD 115600 (Gibbs et al. 2019), which also show a red color in their reflectance spectrum.

As seen in Fig. 11, the measured scattered flux is significantly higher in the HD 114082 disk compared to the HD 117214 disk. When comparing these two disks, we have to keep in mind that planetesimal belts are located at different radial distances from their host stars and have different inclinations. For a comparison of measured fluxes, these geometric parameters have to be taken into account. Therefore, the HD 117214 disk fluxes plotted in Fig. 11 are multiplied by a factor of 1.55. This factor is given by the ratio of squared radial distances of the belts used in our measurement: $R_{\text{HD117214}}^2/R_{\text{HD114082}}^2={\left(0.46/0.37\right)}^2=1.55$. The correction factor compensates for the different illumination of disks because of dilution of stellar light with a distance that scales as 1/r².

Likewise, the amount of light scattered by dust particles into the LOS of an observer strongly depends on the disk morphology, in particular the disk inclination. With a larger inclination, more flux should be measured along the disk rim because more dust particles will be in the LOS of the observer due to the projection effect. Using the best-fit model for the HD 117214 disk (Engler et al. 2020), we estimated this effect by varying the disk inclination between 73° (inclination of HD 117214 disk) and 83° (inclination of HD 114082 disk). In the disk model with a higher inclination, the scattered flux measured within the same circular apertures placed at the disk ansae is 1.75 times higher than in the disk with a lower inclination. Multiplying the HD 117214 data (Fig. 11) by this correction factor for the disk inclination, in addition to the previous correction factor of 1.55 for the flux dilution, we still obtain a ~2.5 lower flux for HD 117214 compared to HD 114082.

Fig. 11 Near-IR reflectance spectra of debris disks HD 114082 and HD 117214 measured with the IFS (solid lines) and IRDIS K1 and K2 bands (dots). The HD 117214 spectrum has been additionally scaled by a factor of 1.55 for comparison purposes. The spectra are normalized by the stellar fluxes shown by the dashed lines for the IFS data. Stellar flux values are given on the right y-axis of the plot. Black dotted lines show the location of the CO2 ice absorption bands at 1.21, 1.44 and 1.54–1.58 μm. The blue-shaded regions schematically indicate the absorption bands of the atmospheric water vapor.

6.2 Constraints on the disk albedo

We define a spectral albedo of a debris disk ω_{λ disk} as a relation between the amount of attenuation of the stellar flux F_{* λ} by its scattering off dust particles in the whole debris disk (including reflection, refraction and diffraction) and the total amount of stellar flux attenuation by dust scattering and absorption. The disk spectral albedo can be expressed in terms of the averaged contrast of the disk representing the disk scattering cross section, σ_{λ disk}, at a wavelength λ and the fractional IR luminosity of the disk, L_{IR disk}/L_*, which can be used as a proxy for the averaged absorption cross section of the disk κ_disk: $${\omega }_{\lambda \text{\hspace{0.17em}disk}}=\frac{{\sigma }_{\text{sca}}{}_{\lambda \text{\hspace{0.17em}disk}}}{{\sigma }_{{{}_{\text{sca}}}_{\lambda \text{\hspace{0.17em}disk}}}+{\kappa }_{\text{disk}}}=\frac{F_{\text{sca\hspace{0.17em}}\lambda }/F_{*\lambda }}{F_{\text{sca\hspace{0.17em}}\lambda }/F_{*\lambda }+L_{\text{IR\hspace{0.17em}disk}}/L_{*}}.$$(9)

The scattering albedo ω_disk for the HD 114082 disk at near-IR wavelengths can be evaluated based on the best-fitting model 1 and model 2. With their help we are able to estimate how much flux is scattered into the LOS of the observer and the entire space.

For this purpose, we measured the scattered flux in the K1 band image of total intensity (Fig. 3a) within two rectangular boxes with a height of 0.16” centered on the disk major axis and covering the image area between radial distances r = 0.12″ and r = 0.54″ from the star. For the east side of the disk we obtained a total flux of 410 ± 20 cts s⁻¹, for the west side 340 ± 18 cts s⁻¹. This flux is significantly reduced compared to the actual disk flux scattered into the LOS because of the ADI-processing of the data. Using the models 1 and 2 we derive the total intrinsic flux for the left and right boxes together of 4400 ± 500 cts s⁻¹.

The stellar flux in the K1 band sums up to (1.13 ± 0.00) × 10⁶ cts s⁻¹. This gives a contrast for the measured total intrinsic flux of (3.9 ± 0.5) × 10⁻³. With this contrast we can estimate the effective spectral disk albedo, meaning the albedo derived with the contrast of the probed regions only. Since one part of the disk is behind the coronagraph, the fractional IR luminosity in Eq. (9) needs to be reduced by a factor of 0.93 taking into account only visible part of the disk. With this roughly estimated factor we obtain the effective spectral albedo for the HD 114082 disk of 0.58: $$\begin{array}{l}{\omega }_{\text{eff\hspace{0.17em}}\lambda \text{\hspace{0.17em}HD\hspace{0.17em}}114082}=\frac{\left(3.90\pm 0.50\right)\times {10}^{-3}}{\left(3.90\pm 0.50\right)\times {10}^{-3}+0.93\cdot 3.01\times {10}^{-3}}\hfill \\ =0.58\pm \mathrm{0.08.}\hfill \end{array}$$(10)

We note that this value does not include scattered flux from the region behind the coronagraphic mask where the peak of forward-scattering occurs for the highly inclined disks such as HD 114082 disk. Moreover, the effective albedo derived with the flux scattered into LOS depends on disk inclination. This means that if the HD 114082 disk would have a different inclination, then the effective albedo would be also different. Therefore, in the following paragraphs, we estimate the average spectral albedo for HD 114082 disk based on the modeled disk luminosity.

The total scattered flux in the images of non-convolved and unprocessed models for the K1 band data is (2.0 ± 0.4) × 10⁴ cts s⁻¹. The scattered flux averaged over the full sphere solid angle, which we define as the scattering luminosity of the disk divided by 4π, is F_{sca λ} = (6.32 ± 1.25) × 10³ cts s⁻¹. The uncertainty on these values mostly results from the unknown amount of scattered flux from the disk part behind the coronagraph. Models with a large HG asymmetry parameter ɡ describing the diffraction peak of forward-scattering, like for instance model 2 (ɡ₁ = 0.83), predict much higher unseen flux compared to the models with a smaller parameter g, like model 1 (ɡ = 0.56).

This yields an averaged contrast of the disk F_{sca λ}/F_{* λ} = (5.61 ± 1.10) × 10⁻³ in scattered light at λ = 2.110 μm. Based on the estimated averaged flux value, we derive a relatively high mean albedo of 0.65 ± 0.15 for the HD 114082 disk: $${\omega }_{\lambda \text{\hspace{0.17em}HD\hspace{0.17em}}114082}=\frac{\left(5.61\pm 1.10\right)\times {10}^{-3}}{\left(5.61\pm 1.10\right)\times {10}^{-3}+3.01\times {10}^{-3}}=0.65\pm \mathrm{0.15.}$$(11)

This value is, however, close to the derived effective spectral albedo of 0.58.

Our analysis of the SPF curve in Sect. 5.1 constrains the effective particle size in the disk observed at near-IR wavelengths to be within a range of 5–15 μm. This relatively large effective particle size and the low maximum PF of 17% (Sect. 5.3) agree well with the predictions of light scattering theory for high-albedo particles. Aggregates of debris particles with relatively large sizes tend to produce multiple scattering of stellar photons and, thus, generate a higher scattered flux similar to high-albedo surfaces. Due to multiple scattering stellar light becomes less polarized compared to the single scattering event. This depolarization effect, the so-called Umov effect (Umov 1905), could explain the low PF of scattered light in the HD 114082 disk.

The scattering albedo of the HD 117214 disk might be lower than the value we estimated for the HD 114082 disk. According to the results of the spectral analysis (Sect. 6.1), the micron-sized particles from HD 117214 scatter ~2.5 times less stellar flux than the particles from the HD 114082 disk at scattering angle θ = 90°. This means that, if the scattering albedo of the dust particles in both disks were similar, the mass of micron-sized particles we trace in both observations should be ~2.5 times lower in HD 117214 compared to HD 114082.

Assuming the same particle size distribution in both systems, we can roughly estimate the actual relation between these masses by comparing the disk IR luminosities L_{IR disk}. Since the particle temperatures we derived for both disks (see Eq. (8)) are approximately the same, the IR luminosities of disks should scale with the dust masses of the micron-sized particles M_dust as $$M_{\text{dust}}\propto L_{\text{IR\hspace{0.17em}disk}}=f_{\text{IR\hspace{0.17em}disk}}\cdot L_{*},$$(12)

where f_{IR disk} is the fractional IR luminosity of the disk, which is equal to 3.01 × 10⁻³ for the HD 114082 disk and 2.53 × 10⁻³ for the HD 117214 (Jang-Condell et al. 2015). The ratio of the disk masses should scale as $$\frac{M_{\text{disk}}\left(\text{HD\hspace{0.17em}}117214\right)}{M_{\text{disk}}\left(\text{HD\hspace{0.17em}}114082\right)}\sim \frac{2.53\times {10}^{-3}\cdot 5.64L_{\odot }}{3.01\times {10}^{-3}\cdot 3.86L_{\odot }}=\mathrm{1.23.}$$(13)

If our assumptions are valid, the HD 117214 disk should possess a higher mass than the HD 114082 disk. This conclusion is consistent with the estimated stellar mass being 1 .14 times higher for HD 117214 compared to HD 114082 (Jang-Condell et al. 2015; Pecaut et al. 2012). This indicates that the mass of the original molecular cloud of the former was higher, and thus the mass of debris disk surrounding it is most likely also higher. Higher dust mass and lower scattered flux mean that the albedo of micron-sized particles in HD 117214 disk maybe 2–2.5 times lower at near-IR wavelengths compared to the HD 114082 disk.

Different albedo of the disks can be explained by different optical properties of dust: slightly different effective sizes, shapes, or composition of dust particles. If the particle size is the most important parameter affecting the amount of scattered flux, then, from a comparison of the particle reflectivity between the HD 114082 and HD 117214 disks, we conclude that the effective particle size in the HD 117214 disk is smaller than 5 μm. If the particle composition is responsible for the lower reflectivity of dust in HD 117214 disk, we can think of particles composed of “dirty” CO₂ and H₂O ice particles with inclusions of dark material, for example complex organic material (Levasseur-Regourd et al. 2018). The reflectance of the CO₂ ice particles could be lower than 50% in the near-IR, as the spectroscopic studies of the Martian polar regions covered by a layer of the CO₂ ice report (Appéré et al. 2011).

Different composition of dust particles, in particular the ice volume fraction, in HD 114082 and HD 117214 disks could be a direct consequence of different evolutionary stages at which both stars might be. In Appendix C, we discuss in more details the possibility for HD 117214 to be much younger than the median age of 17 Myr assumed for the F stars in the LCC. At younger ages, the circumstellar disks possess more gas and ices. Their amounts decrease over time due to gas release, when icy particles collide with other debris, evaporate or undergo photo-desorption with subsequent gas photodissociation caused by the stellar and interstellar radiation (e.g., Grigorieva et al. 2007) or due to accretion onto growing giant planets.

The upcoming studies with the James Webb Space Telescope (JWST) can shed light on the content of various ice species, such as H₂O, CO, CO₂, CH₃OH, NH₃, and CH₄, in the dust particles of debris disks surrounding stars at different ages. In particular, the JWST Near Infrared Camera (NIRCam) will provide sensitive imaging and spectroscopy at near-IR wavelengths from 0.6 to 5.0 μm which includes H₂O and CO₂ frost features. The future JWST observations can reveal the locations of ices, and thus locations of various snow lines in the disks and possibly clarify how much primordial gases remain bound in the debris around young stars entering the main sequence.

7 Planetary system of HD 114082

The planetary system of HD 114082 consists of at least two spatially separated planetesimal belts: a bright exo-Kuiper belt with a radius R_belt of 0.37″ (35 au) discussed in this work and a warm dust belt. The latter was identified from the disk’s spectral energy distribution, which is better matched with a thermal emission from two separate dust populations. The radial location of the warm dust belt is however unclear: different studies placed the second dust component at different radial positions between and 18 au from the star (e.g., Mittal et al. 2015; Jang-Condell et al. 2015), and even a dust population with a temperature of 500 K at 0.8 au provided a good fit to the spectrum measured with the Spitzer InfraRed Spectrograph (Chen et al. 2014). Therefore, the HD 114082 disk might harbor multiple components with inner debris rings that cannot be resolved with current instruments because of small angular sizes (the radius of the coronagraphic mask is 0.12″ or 11 au in our observations) and high inclination of the bright outer belt.

No planetary companions have been found in the cavity of the outer belt with direct imaging so far. The analysis of recent TESS data shows that a Jupiter-sized planet resides within the cavity and is hidden by the coronagraph in our SPHERE observations (Fig. 12).

7.1 Discovery of a transiting planet

HD 114082 was observed by TESS in sector 38 (2021 April 28 to 2021 May 26) at a 2-min cadence and additional observations are planned in sectors 64 and 65 (2023 April 06 to 2023 June 02). A single transit event is detected and announced as a community TESS object of interest (CTOI) in July 2021. This single transit event is not concurrent with any asteroid crossing or centroid shifts, which could indicate that the signal is on a background star. The transiting candidate is on target and passed the first photometric vetting stage.

We downloaded the TESS light curve available on the Mikulski Archive for Space Telescopes⁴ and we selected the Pre-search Data Conditioning Simple Aperture Flux (PDCSAP) and associated errors to use in our modeling. The PDCSAP flux is a product of the data reduction done at the Science Processing Operation Center (SPOC; Jenkins et al. 2016). The PDCSAP flux is corrected for systematic trends and the photometric dilution due to contamination by other stars in the aperture.

We perform the modeling of the photometric data with the software package Juliet (Espinoza et al. 2019). The transit model is provided by batman (Kreidberg 2015) and we model the stellar activity with a Gaussian process. The fit is done using the nested sampler dynesty implemented in Juliet. We set a log-uniform prior for the orbital period between 21.5 and 1000 days, the lower period bound being the minimum period allowed by the TESS data. We chose to fit for the stellar density with a normal period based on the values from the TICv8 catalog (Stassun et al. 2019). The stellar density information combined with the transit shape (e.g., transit duration, impact parameter) constrain the speed at which the planet is crossing the stellar disk and thus give us an estimate of its orbital period. Figure 13 shows the detrended TESS light curve together with the median transit model and the associated residuals.

We find that the planetary candidate corresponds to a ${0.977}_{-0.025}^{+0.026}$ R_J with an orbital Period of ${197}_{-109}^{+171}$ days, which corresponds to a semimajor axis of ${0.7}_{-0.3}^{+0.4}$ au. The total transit duration equals to ${14.58}_{-0.06}^{+0.06}$ h, and the impact parameter is ${0.42}_{-0.18}^{+0.23}$ that translates into inclination of the orbit of ${89.78}_{-0.25}^{+0.10}$ degrees.

In parallel with our work another study describing the discovery of this transiting planetary candidate appeared on arXiv (Zakhozhay et al. 2022b). This study combines the constraints on the orbital period and planet mass from the TESS transit event in 2021 and RV data obtained from the 4-yr monitoring campaign of HD 114082 with the Fibre-fed Extended Range Échelle Spectrograph (FEROS) and publicly available spectra from the High Accuracy Radial velocity Planet Searcher (HARPS; Mayor et al. 2003) archive. The authors confirm the presence of a close-in giant planet in the system, which was assigned the identification name HD 114082 b. According to the analysis by Zakhozhay et al. (2022b), the planet has a mass of 8.0 ± 1.0 M_Jup and moves around the star on an eccentric orbit (e = 0.40 ± 0.04) with a semimajor axis of a = 0.51 ± 0.01 au and orbital period of 109.8 ± 0.4 days. The last two parameters are within 1σ-interval of our parameter values from the fitting of TESS light curve.

Since the orbital plane of HD 114082 b is almost edge-on, it is non-coplanar with the disk midplane and has an offset of more than six degrees from it. The possible ellipticity of the planet orbit is not apparent in the disk images presented in this work: the disk morphology is consistent with the rotationally symmetric model for the dust distribution.

The non-coplanarity of planet and disk planes might be a result of interaction of HD 114082 b with the other planets residing inside the outer planetesimal belt. If such planets exist and move in the disk plane, they do not transit the star and can be currently detected only with RV measurements. We conclude thus that further photometric and RV monitoring of HD 114082 is needed to confirm the parameters of HD 114082 b and its orbit and investigate the presence of other planets in this system.

Fig. 12 Illustration of the transiting planet orbit. The orbit and planet sizes are shown schematically and not in the scale indicated in the upper left corner: the radius of the disk is 35 au, whereas the radius of the planet orbit is 0.7 au.

Fig. 13 Detrended TESS light curve with the median transit model (black line) and the associated residuals. The standard deviation of the TESS residuals is 380 ppm.

7.2 Companion candidates in the SPHERE image

The SHINE observations of HD 114082 revealed 23 companion candidates in the IRDIS FOV, with the S/N ranging from ~3 up to ~70 and projected semimajor axes from ~60 au up to more than 500 au (Fig. 14). To assess their status, the SHARDDS data from 2016 February 14 were reprocessed using the open source vlt-sphere⁵ python pipeline (Vigan 2020). A PCA analysis with 5 subtracted modes was applied to the reprocessed data cubes and all the SHINE candidates are clearly re-detected in the SHARDDS data. Due to the time baseline of more than 1 yr between the SHARDDS and SHINE observations, and to the relatively high proper motion of HD 114082 (44.4 mas yr⁻¹), the motion of the candidates between the two epochs is obvious. We simply use a Gaussian fit of their PSF in the SHARDDS data to obtain their astrometry with an accuracy estimated to approximately half an IRDIS pixel (i.e., 6.1 mas).

A plot of the relative motion of all the candidates with respect to the central star between the two epochs is presented in Fig. 15. The detailed astrometry is also reported in Table D.1. The relative motion of all the candidates clearly follow the track expected for stationary background star, so we conclude that all candidates are background stars. Such a relatively large number of background stars in the IRDIS FOV is explained by the proximity of HD 114082 to the galactic plane: its galactic latitude is 2.5° (e poch = J2000).

Fig. 14 Companion candidates for HD 114082 detected with SPHERE/IRDIS. The color bar shows flux units normalized by the stellar flux.

Fig. 15 Relative motion of the candidate companions with respect to the central star between the SHARDDS and SHINE observations, on 2016 January 14 and 2017 May 17, respectively. The blue curve represents the track expected for stationary background objects. From this analysis we conclude that all candidates detected in the SHINE data are background stars.

7.3 Companion mass sensitivity limits combining direct imaging, radial velocity, and Gaia

To bring constraints on the presence of additional companions in the system, we combined the constraints from direct imaging with those from RV (HARPS and FEROS) and from the Gaia astrometry. We used the Multi-epochs multi-purposes Exoplanet Simulation System 3 (MESS3) tool (Kiefer et al., in prep.), an extension of the MESS2 tool (Lannier et al. 2017) that combines only direct imaging and RV data.

In direct imaging, we used as input the sensitivity maps deduced from the IFS image of the night of 2017 May 16 reduced in spectral imaging and ADI, and the IRDIS broadband H image of the night of 2016 February 13 reduced with PCA. Two mass sensitivity maps were therefore created assuming the AMES-COND evolutionary model (Baraffe et al. 2003), and used as input for MESS3. We also retrieved the archival HARPS and FEROS RV data of the target as described in Zakhozhay et al. (2022b). We complemented that with the excess noise and proper motion anomaly from Gaia Data Release 3 (Gaia Collaboration 2021), with the technique described in Kiefer et al. (in prep.). For the generation of orbits of planetary companions, we assumed inclinations distributed between 0° and 90° and eccentricities between 0 and 0.9.

The detection limits obtained when combining those three techniques are displayed in Fig. 16 where we plot the isocontours of probability of detection, in a mass versus semimajor axis diagram. At a semimajor axis of 0.7 au (the most likely semi-major axis of the transiting planet based on our analysis), the 1 − σ sensitivity (68% probability of detection) corresponds to companions more massive than 30 M_Jup. We can exclude bound companions more massive than 40 M_Jup on any orbits with the same confidence interval. Exterior to the disk, at 50 au we reach a sensitivity down to ~5 M_Jup, and interior to the disk at 30 au we are sensitive to planets more massive than ~10 M_Jup.

8 Summary

In this work we present the new optical and near-IR images of scattered light (total and polarized intensity) from the debris disk surrounding the young F star HD 114082, a member of the Sco-Cen association. The images were obtained using the SPHERE instruments IRDIS, IFS, and ZDVIPOL. They show a very bright axisymmetric planetesimal belt with a radius of ~0.37″ (35 au) and an inclination of ~83°. The disk has a PA of 106° ± 2° (east of north).

Our modeling results show that the distribution of disk brightness in the image of total intensity can be well reproduced with a SPF that combines two HG functions with asymmetry parameters ɡ₁ = 0.83 (weight 0.73), describing the forward-scattering peak, and ɡ₂ = 0.06 (weight 0.27), describing a flat middle region of the SPF. When compared to the diffraction part of theoretical SPFs calculated for the microasteroids (Fig. 5), the measured SPF of the HD 114082 disk might indicate that the effective size of dust particles traced in the SPHERE observation at near-IR wavelengths is slightly larger than or close to 5 μm.

We used polarimetric data of HD 114082 to compare the disk morphologies seen in total and polarized intensity images (Figs. 1 and 2) and to estimate the degree of linear polarization of scattered light. Comparing IRDIS DPI observations in the H band to the IFS data that cover the same wavelength range, we obtained a maximum PF of 17 ± 4% reached between scattering angles of 70° and 115° and a minimum PF of less than 2% between 140° and 155°. The decrease in PF at scattering angles larger than 115° is similar to the behavior of cometary dust in the Solar System, which shows zero PF at θ ≈ 160°. The overall shape of the disk PF resembles the bell-shaped PF measured for the cometary dust analogs.

A comparison of the near-IR spectra of two debris disks, HD 114082 and HD 117214, reveals a higher reflectivity for HD 114082 dust particles, with an albedo of 0.65 ± 0.15. Although both host stars are young F-type stars in the LCC subgroup, their debris disks most likely contain particles with different optical properties in the near-IR. This can be explained by a different composition or structure of debris material that can change with the stellar age. The nearly twice higher luminosity of HD 117214 compared to the luminosity of HD 114082 suggests that HD 117214 is either a close binary or much younger than the adopted median age of 17 Myr for the LCC members. In the latter case, HD 117214 would be at an earlier evolutionary stage than HD 114082 and could harbor a more gas- and ice-rich disk.

The reflectance spectra of both disks, measured along the disk major axes at the radial position of planetesimal belts, show a rather flat flux density contrast within the wavelength range 0.96 to 1.66 μm. It is lower than the flux density contrast derived for the IRDIS K1 and K2 filters (λ_K1 = 2.110 μm λ_K2 = 2.251 μm), where the reflectance spectra show a red color. Apart from that, the HD 114082 spectrum appears to be featureless, whereas the HD 117214 spectrum may indicate absorption features. They could be interpreted as the absorption bands of CO₂ ice. However, to conclude about the actual content of CO₂ and other ices in the dust particles of HD 117214 or other debris disks, space-based observations with high spatial and angular resolutions in the near-IR are needed. For some debris disks, detection of the ice absorption bands will be possible with the instruments on board the JWST.

The TESS photometric data show a dip in the light curve of HD 114082, indicating that a Jupiter-sized planet revolves around the star within the cavity of the outer planetesimal belt in an orbit with a semimajor axis of 0.7 ± 0.4 au. Thus far, the transit of the planet has only been observed in 2021, and its parameters remain to be confirmed via further photometric and RV monitoring.

Within the IRDIS FOV, we detected 23 companion candidates in the HD 114082 data. Based on a comparison with available SHARDDS images, we conclude that all these candidates are background stars. By combining the SPHERE data with archival RV data and astrometry, we reach a high sensitivity down to ~5 M_Jup at 50 au, and ~ 10 M_Jup at 30 au when using the AMES-COND evolutionary models to convert the derived contrast limits into the companion mass.

Fig. 16 Sensitivity limits on bound companions, obtained by combining the constraints from IRDIS/IFS direct imaging, HARPS/FEROS RV and Gaia astrometry. The darkest area shows the parameter space where the probability of detecting a planet is less than 10%, while the white area corresponds to a probability higher than 99%.

Acknowledgements

We thank the anonymous referee for many thoughtful comments which helped to improve this paper. This work has made use of the SPHERE Data Centre, jointly operated by OSUG/IPAG (Grenoble), PYTHEAS/LAM/CeSAM (Marseille), OCA/Lagrange (Nice), Observatoire de Paris/LESIA (Paris), and Observatoire de Lyon/CRAL, and supported by a grant from Labex OSUG@2020 (Investissements d’avenir - ANR10 LABX56). This research has made use of the Direct Imaging Virtual Archive (DIVA), operated at CeSAM/LAM, Marseille, France. This project has received funding from PSL-OCAV and from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (COBREX; grant agreement no. 885593). AV acknowledges funding from the ERC under the European Union’s Horizon 2020 research and innovation program (grant agreement no. 757561). JPM acknowledges research support by the Ministry of Science and Technology of Taiwan under grants MOST107-2119-M-001-031-MY3 and MOST109-2112-M-001-036-MY3, and Academia Sinica under grant AS-IA-106-M03. CdB acknowledges support by Mexican CONACYT research grant FOP16-2021-01-320608. AZ acknowledges support from the FONDECYT Iniciación en investigación project number 11190837 and ANID – Millennium Science Initiative Program – Center Code NCN2021_080. SPHERE is an instrument designed and built by a consortium consisting of IPAG (Grenoble, France), MPIA (Heidelberg, Germany), LAM (Marseille, France), LESIA (Paris, France), Laboratoire Lagrange (Nice, France), INAF — Observatorio di Padova (Italy), Observatoire de Genève (Switzerland), ETH Zurich (Switzerland), NOVA (Netherlands), ONERA (France) and ASTRON (Netherlands), in collaboration with ESO. SPHERE was funded by ESO, with additional contributions from CNRS (France), MPIA (Germany), INAF (Italy), FINES (Switzerland) and NOVA (Netherlands). SPHERE also received funding from the European Commission Sixth and Seventh Framework Programmes as part of the Optical Infrared Coordination Network for Astronomy (OPTICON) under grant number RII3-Ct-2004-001566 for FP6 (2004-2008), grant number 226604 for FP7 (2009-2012) and grant number 312430 for FP7 (2013-2016). This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

Appendix A U_φ images

Figure A.1 shows the modeled image of the IRDIS U_φ signal after convolution of the modeled Stokes parameters Q and U with the IRDIS PSF from the night of observation. We used model 2 to create the images. The convolution of the measured Stokes parameters Q and U with the instrument PSF always generates an artificial U_φ signal and results in the reduction of the polarized flux measured in the Q_φ image. This effect is stronger for the disks with smaller angular sizes, such as the HD 114082 disk (Engler et al. 2018).

Figure A.2 shows the ZIMPOL U_φ image, the same as in Fig. 2 but here with different scale range and colors. This image does not contain a signal from the disk.

Fig. A.1 IRDIS Uφ model image.

Fig. A.2 ZIMPOL Uφ image.

Appendix B Modeling the scattered (polarized) light from a debris disk

Model for the image of scattered (polarized) light

To create a synthetic image of light scattered by small dust particles, we use a 3D rotationally symmetric model for the spatial distribution of particle number density in a debris disk. We describe such a distribution by the product of the radial power law R(r) (Eq. B.1, Augereau et al. 1999) and the Lorentzian profile Z(r, h) for the vertical distribution of particle density (Eq. B.3). Compared to the exponential profile adopted by Augereau et al. (1999) for the study of HR 4796 A disk, the Lorentzian profile offers a model with a smaller number of parameters. This profile seems to fit the vertical cuts of debris disks well, and was used for the disk modeling in previous studies (e.g., Krist et al. 2005; Engler et al. 2018).

Our model is based on the theory of the parent body belt that consists of massive planetesimals moving on a circular orbit with a radius R_belt. Mutual collisions between the planetesimals produce a large amount of micron-sized dust particles that are radially redistributed with time. The maximum density of particles occurs at the site of their production, namely at the radial distance R_belt. Therefore, the particle number density in the disk midplane can be described by the radial power laws with the exponents α_in > 0 inside the belt and α_out < 0 outside the belt: $$R\left(r\right)\propto {\left({\left(\frac{r}{r_0}\right)}^{-2{\alpha }_{in}}+{\left(\frac{r}{r_0}\right)}^{-2{\alpha }_{out}}\right)}^{-1/2},$$(B.1)

where r₀ is the reference radius that is related to the belt radius R_belt as $$R_{\text{belt}}=r_0{\left(-\frac{{\alpha }_{in}}{{\alpha }_{out}}\right)}^{1/\left(2{\alpha }_{in}-2{\alpha }_{out}\right)},$$(B.2)

and thus R_belt = r₀ if α_in = −α_out.

We define the vertical profile for the particle density distribution in disk cross sections perpendicular to the disk midplane according to $$\text{Z}\left(r,h\right)\propto {\left[1+{\left(\frac{h}{H\left(r\right)}\right)}^2\right]}^{-1},$$(B.3)

where h is the height above the disk midplane. The scale height of the disk H(r) is defined as a half-width at half-maximum of the vertical profile at radial distance r and scales like H(r) = H₀ (r/r₀)^β, where H₀ = H(r₀) and β is the disk flare index.

How much stellar light is scattered into LOS of observer depends on the scattering angle θ. This dependence is given by the SPF, which is often approximated by the HG function (Henyey & Greenstein 1941): $$f\left(\theta ,g\right)=\frac{1-g^2}{4\pi {\left(1+g^2-2gcos\left(\theta \right)\right)}^{3/2}},$$(B.4)

where ɡ is the HG scattering asymmetry parameter $\left(-1⩽g⩽1\right)$, and θ is the scattering angle measured between direction of the incoming stellar radiation and the direction of the scattered radiation. The latter is the LOS of the observer.

We assume that dust particles scatter more radiation in the forward direction (at smaller scattering angles). This means that the asymmetry parameter ɡ has a positive value and the brighter side of the disk is closer to the observer. We also consider that dust particles have the same properties everywhere in the disk and scattered-light images preferentially trace dust particles with a size comparable to the wavelength of observation.

As mentioned in Sect. 4.3, in this work we tested two models with different SPFs. The phase function of model 1 is the HG function (Eq. B.4). The phase function of model 2 is a linear combination of two HG functions (Engler et al. 2017): $$f\left(\theta ,g_1,g_2\right)=w\cdot f\left(\theta ,g_1\right)+\left(1-w\right)\cdot f\left(\theta ,g_2\right),$$(B.5)

where the first parameter ɡ₁ describes a forward-scattering peak of the SPF, the second parameter ɡ₂ represents the more isotropic part of the SPF, and ω is the scaling parameter, $0⩽\omega ⩽1$.

Fig. B.1 Comparison of the PF derived with models 1 and 2.

For the modeling of polarized-scattered flux, the pSPF is used instead of SPF. The pSPF is given by a product of SPF and PF (Eq. 5). The PF can be approximated with an analytic expression for the small dust particles (e.g., Graham et al. 2007; Engler et al. 2017) or measured from the data, if the data quality is good.

Using this model, we create a 2D image of scattered or polarized light from a debris disk by summing up all contributions to the scattered (polarized) flux F_i along the LOS of the observer at each pixel (resolution element) i of the image: $$F_i=A{\displaystyle \underset{LOS}{\int }\frac{f\left(\theta ,g\right)R\left(r\right)Z\left(r,h\right)}{r^2}},$$(B.6)

where A is a scaling factor, and the SPF f(θ, ɡ) is replaced by the pSPF for the polarized flux image.

The obtained model image can be compared to the science image using the reduced ${\chi }_v^2$ metric (Eq. 3).

To determine the best-fit parameter set for the total intensity image in the K1 band (Fig. 1, left panel), we used a custom MCMC code employing the Python package emcee (Foreman-Mackey et al. 2013). We ran the MCMC sampler with 1000 walkers for each model family separately using the uniform priors as specified in Col. 2 of Table 2.

At each walker step, a synthetic image of scattered light for a model with a given parameter set $p=\left(p_1,p_2,\dots ,p_{N_{\text{par}}}\right)$ is created. This image is then inserted in a data cube at different PAs to reproduce the rotation of the sky field during the observation and perform the cADI forward-modeling. The forward model is used to estimate the likelihood ℒ of parameter set p by calculating the χ² metric according to Eq. 3 and the In ℒ according to Hogg et al. (2010): $$\ln ℒ=-\frac{1}{2}\left({\chi }_v^{2}v+{\displaystyle \sum _{i=1}^{N_{\text{data}}}\text{ln}\left(2\pi {\sigma }_{i,\text{data}}^2\right)}\right).$$(B.7)

Fig. B.2 Polarization fraction of scattered light derived for model 2 with different inclinations and scale heights of the disk.

The MCMC sampler converges to the parameter values with the highest likelihood after two thousand steps providing the posterior distributions of the fitted parameters, which we obtain after removing the burn-in phase consisting of first 400 walker steps. The derived posterior distributions are shown in Figs. E.1 and E.2. Their median values are given in Cols. 3 and 4 of Table 2. The lower and upper uncertainties on these values are defined as the 16th and 84th percentiles of the posterior distributions.

Polarization fraction phase functions of models 1 and 2

Figure B.1 shows a comparison between forward-modeled PF for model 1 and model 2. After two iterations a very good match between the measured PF from the data and measured PF from both models was obtained.

Impact of disk inclination and scale height on the measured polarization fraction of scattered light

The outcome of PF measurement from the images of debris disks in scattered and polarized light is significantly affected by the disk morphology, in particular by its radial and vertical extent. We referred to this effect as the LOS effect in Sect. 5.3.3. Since this is a geometrical effect, we expected the inclination and scale height to be the disk parameters that most influence the measured PF value. To investigate this, we used model 2 to obtain the measured PFs for disks with different inclinations and scale heights. The inclinations were varied from 10° (nearly pole-on disk) to 90° (edge-on disk). The scale heights were set to 0.004″ (H₀/r₀ = 0.01) and 0.04″ (H₀/r₀ = 0.1).

The PFs shown in Fig. Β.2 were derived applying aperture photometry to the non-convolved model images of the total and polarized intensities. They demonstrate that the measured PF decreases with the growing disk inclination and scale height and that the position of maximum value is shifted to the larger scattering angles for the disks with inclination larger than 80°.

Appendix C On the age of HD 117214

In Sect. 6 we discussed different reflectivity of dust particles in debris disks HD 114082 and HD 117214. The difference in optical properties between these two disks could arise from different properties of their host stars. It is assumed that HD 117214 and HD 114082 are approximately 17 Myr old (Pecaut et al. 2012). Since they are both F-type stars with a similar surface temperature, the much higher brightness of HD 117214 (1.5-2 times brighter than HD 114082) raises the question of why this star is much more luminous than expected for its estimated age. There are two possible answers to this question: (1) the star is a close (spectroscopic) binary with two stars of nearly equal masses, or (2) the HD 117214 is significantly younger than 17 Myr old, the age that is commonly adopted for the LCC members.

Indeed, HD 117214 is not a single star. According to the Gaia Early Data Release 3 catalog (Gaia Collaboration 2021), it is accompanied by at least one M dwarf (2MASS J13300763-5829026), which has a G magnitude of 13.4^m and a contrast to the central star of dG = 5.44^m. This is certainly a bound companion to the HD 117214 since both stars have the same parallax of about 9.3 mas within the error bars and a similar proper motion of 38-39 mas/year. However, this late-type companion has a projected separation of 10.72″ or 1150 au (Gaia Collaboration 2021) and therefore cannot be responsible for the over-luminosity of HD 117214. If this companion is on a circular orbit with the same inclination of ~70° as the debris disk (Engler et al. 2020), then its actual separation from the primary star is ~3500 au.

Fig. C.1 TESS photometric time series (top panel), the ACF (middle panel), and periodogram of HD 117214 (bottom panel).

Fig. C.2 Full width at half maximum time variations for HD 117214.

Fig. C.3 Color-magnitude diagram for the HD 117214 primary and secondary stars. The extinction vector just indicates the direction of reddening.

However, HD 117214 might be itself an unresolved spectroscopic binary. If so, the two components of star should be separated by at most a few astronomical units. This follows from our analysis of the SHINE survey data (Engler et al. 2020) and proper motion analysis of this system based on Hipparcos and Gaia Data Release 2 data by Kervella et al. (2019). In such a case, the binarity of HD 117214 can be tested by means of an analysis of the photometric variability or RV data.

In fact, we know that unresolved close binaries consisting of components of similar brightness, as assumed in our hypothesis, and of late spectral type (equal to or later than the F type), generally host activity phenomena, like starspots, which can induce quasi-periodic variation of their flux. A periodogram analysis of photometric time series can reveal the rotation period of each component, thus unveiling the binary nature of the system (see, e.g., Rebull et al. 2018; Messina 2019). Indeed, our target was observed by TESS, and we analyzed its photometric time series (top panel of Fig C.1) using both the auto-correlation function (ACF; McQuillan et al. 2013) and the Lomb-Scargle (Horne & Baliunas 1986) periodogram analysis.

The ACF (second panel of Fig C.1) displays the typical pattern with one or two lower-amplitude local maxima followed by the highest amplitude maximum. This pattern is generally produced when the spots responsible of the photometric variability are unevenly distributed across the stellar surface. The highest local maximum at P = 1.644 d is assumed to represent the stellar rotation period.

This interpretation is supported by the results of the Lomb-Scargle periodogram (bottom panel of Fig C.1) in which we observe two major peaks, at P = 0.551± 0.006 d and P = 0.844 ± 0.012 d. The latter period is half the rotation period inferred from ACF, confirming the presence of two spot groups on opposite hemispheres. The P = 0.551 d is a beat period of P = 0.844 d as well as the other lower-power peaks at P = 2.94±0.13 d and P = 2.56±0.10 d are all aliases. Therefore, no evidence of binarity is inferred from our periodogram analysis.

The RV measurements of this star do not show significant variations on the time baseline of several years. Chen et al. (2011) reported a RV of 9.7 ± 0.3 km/s whereas Gaia’s measurements resulted in 10.2 ± 0.2 km/s. The IR spectra of HD 117214 obtained with the HARPS have been analyzed by Grandjean et al. (in prep.), who found no indication of RV variation over a period of 682 days (measured amplitude of RV variations is 96 m/s). This makes it quite difficult but not completely impossible for the star to be an equal mass binary. In fact, if two stellar components have very similar brightness, with a RV difference not larger than a half of v sin i (40 km/s, Chen et al. 2011), lines of two components in spectrum are blended. In this case the barycenter of the observed line does not change along the orbit but the width of the line should change. However, the full width at half maximum is stable around 56.4 km/s over time in HARPS data (Fig. C.2), and so the variation of the line width is not observable.

Since available RV measurements do not suggest that HD 117214 is a close binary, the analysis of photometric variability and RV data seems to rule out the possibility for HD 117214 to be a close binary. We suggest thus that the star may be much younger than the median age of 17 ± 0 − 8 Myr found for the F-type members of the LCC (Pecaut et al. 2012).

Both the central star and the late-type companion are much redder than expected for the age of 17 Myr (see the compilation of available photometry for the primary and secondary stars in Table C.1). If corrected for the distance, the absolute G magnitude of companion is 8.258^m and its bp-rp color is 2.864. This is inconsistent with the color-magnitude diagram isochrones from the BT-Settl CIFIST evolutionary models by Baraffe et al. (2015). Such an inconsistency can be explained either by the much younger age of the system or by large interstellar reddening toward the system, which makes stars much redder compared to the stars on the main sequence. If the last case were true, the stars would be intrinsically much brighter and massive, but this would be inconsistent with their spectral classification.

From comparison of the M-type companion position with the isochrones in the color-magnitude diagram (Fig. C.3), the age of this star would be between 2 and 3 Myr. However, extensive comparisons for stars in Sco-Cen show that isochrones underestimate the ages for M-stars (Feiden 2016), likely because they do not properly take into account the impact of stellar activity. Once corrected for this effect, a more appropriate age for the secondary star should be about 5 Myr, if it is not an equal mass binary and does not possess a circumstellar disk itself. Such a young age is not an exception for the LCC members: Murphy et al. (2015) identified two accreting M-dwarfs as young as 5 - 10 Myr in this subgroup.

The primary star is too massive for comparison with the BT-Settle CIFIST isochrones. Using the PARSEC models (Bressan et al. 2012), we find a good match of stellar photometry for stellar ages between 7 and 9 Myr, and masses between 1.63 and 1.53 M_⊙, respectively. Based on the analysis of the primary and secondary stars, we conclude that the age of the HD 117214 might be 7 ± 2 Myr. Our conclusion is supported by the results of the RV Survey for Planets around Young stars (RVSPY; Zakhozhay et al. 2022a), which classifies this target as a pre-main-sequence star. Given the large intrinsic age spread for the LCC members, the HD 117214 might be one of the last formed stellar systems across the region.

If the star is younger than 10 Myr, the upper limits to the masses of companions possibly responsible for the disk gap presented in Engler et al. (2020) (their Fig. 7) should be slightly reduced. In this case the configuration of planetary system with several small planets rather instead of a single giant planet is even more compelling.

Appendix D Companion candidates’ astrometry

Table D.1 presents the astrometry of all detected companion candidates around HD 114082. The 2016-02-14 points correspond to the SHARDDS data, while the 2017-05-17 points correspond to the SHINE data. All candidates are classified as background.

Appendix E Posterior distributions of the fitted parameters

Figures E.1 and E.2 show the posterior distributions of the fitted parameters for models 1 and 2.

Fig. E.1 Posterior distributions of the fitted parameters for a model with one HG parameter (model 1).

Fig. E.2 Posterior distributions of the fitted parameters for a model with two HG parameters (model 2).

References

All Tables

All Figures

	Fig. 1 Total intensity images of the HD 114082 debris disk obtained with the cADI data reduction of the spectrally combined IFS data (left panel) and IRDIS data (rightpanel). The position of the star is marked with a yellow cross. Color bars show the flux values in counts per second (cts s−1).
In the text

Fig. 2 Polarized intensity images of the HD 114082 debris disk. The Qφ images are shown in the left panels and Uφ images in the right panels. IRDIS data taken in the H band are displayed in the top row. ZIMPOL data recorded with the I_PRIME filter are shown in the bottom row. In the IRDIS Qφ image, a strong residual signal at the coronagraph edge is masked out for clarity. The positions of the star are marked with yellow crosses. All images are smoothed via a Gaussian kernel with σkemei = 1 pixel and shown in the same spatial scale, indicated in the top-left image. The color bar shows the flux values in cts s−1.

In the text

	Fig. 3 Comparison of the total intensity image (panel a) with the forward models of scattered light with a single HG function (panel b) and a two-component HG function (panel d). Panels c and e: residual images obtained after subtracting the forward model images from the total intensity image. The color bar shows the flux values in cts s−1.
In the text

	Fig. 4 Scattering phase function of the HD 114082 debris disk normalized by the stellar flux. The gray-shaded area indicates the uncertainty on the mean flux (black line) of the east and west disk sides.
In the text

	Fig. 5 Scattering phase functions measured for different debris disks and dust populations in the Solar System in comparison with the theoretical microasteroid phase functions and diffraction model for micron-sized particles.
In the text

Fig. 6 Polarized scattering phase functions of the HD 114082 debris disk in the H band and I_PRIME filter normalized by the stellar fluxes. The pSPFs corrected for the limb brightening effect are shown with the solid lines. The light red-shaded area indicates the uncertainty on the measured disk contrast in the H band. The light blue-shaded area indicates the uncertainty on the measured disk contrast in the I_PRIME filter.

In the text

	Fig. 7 Polarization fraction of the debris disk HD 114082 derived as a mean value of the measurements on the west and east disk sides in the H band (black line). The gray-shaded area indicates the uncertainty on this mean PF from the aperture photometry. The blue arrows show a lower limit for the PF in the I_PRIME filter.
In the text

	Fig. 8 Polarization fraction of scattered light in debris disk HD 114082. The dotted black line shows a mean value of the PF derived for the west and east disk sides. The green line shows a polynomial fit to it.
In the text

	Fig. 9 Comparison of the Qφ image (panel a) to the model images of polarized light obtained with parameters of model 1 (panel b) and model 2 (panel c). The color bar shows flux in cts s−1.
In the text

	Fig. 10 Degree of linear polarization of scattered light for the debris disk HD 114082 calculated along the disk major axis. The black line shows a mean value of the west and east disk sides. The radial dependence of the PF is significantly moderated by the difference in scattering angles, such as in Fig. 8 due to the LOS and aperture size effects.
In the text

Fig. 11 Near-IR reflectance spectra of debris disks HD 114082 and HD 117214 measured with the IFS (solid lines) and IRDIS K1 and K2 bands (dots). The HD 117214 spectrum has been additionally scaled by a factor of 1.55 for comparison purposes. The spectra are normalized by the stellar fluxes shown by the dashed lines for the IFS data. Stellar flux values are given on the right y-axis of the plot. Black dotted lines show the location of the CO2 ice absorption bands at 1.21, 1.44 and 1.54–1.58 μm. The blue-shaded regions schematically indicate the absorption bands of the atmospheric water vapor.

In the text

	Fig. 12 Illustration of the transiting planet orbit. The orbit and planet sizes are shown schematically and not in the scale indicated in the upper left corner: the radius of the disk is 35 au, whereas the radius of the planet orbit is 0.7 au.
In the text

	Fig. 13 Detrended TESS light curve with the median transit model (black line) and the associated residuals. The standard deviation of the TESS residuals is 380 ppm.
In the text

	Fig. 14 Companion candidates for HD 114082 detected with SPHERE/IRDIS. The color bar shows flux units normalized by the stellar flux.
In the text

	Fig. 15 Relative motion of the candidate companions with respect to the central star between the SHARDDS and SHINE observations, on 2016 January 14 and 2017 May 17, respectively. The blue curve represents the track expected for stationary background objects. From this analysis we conclude that all candidates detected in the SHINE data are background stars.
In the text

	Fig. 16 Sensitivity limits on bound companions, obtained by combining the constraints from IRDIS/IFS direct imaging, HARPS/FEROS RV and Gaia astrometry. The darkest area shows the parameter space where the probability of detecting a planet is less than 10%, while the white area corresponds to a probability higher than 99%.
In the text

	Fig. A.1 IRDIS Uφ model image.
In the text

	Fig. A.2 ZIMPOL Uφ image.
In the text

	Fig. B.1 Comparison of the PF derived with models 1 and 2.
In the text

	Fig. B.2 Polarization fraction of scattered light derived for model 2 with different inclinations and scale heights of the disk.
In the text

	Fig. C.1 TESS photometric time series (top panel), the ACF (middle panel), and periodogram of HD 117214 (bottom panel).
In the text

	Fig. C.2 Full width at half maximum time variations for HD 117214.
In the text

	Fig. C.3 Color-magnitude diagram for the HD 117214 primary and secondary stars. The extinction vector just indicates the direction of reddening.
In the text

	Fig. E.1 Posterior distributions of the fitted parameters for a model with one HG parameter (model 1).
In the text

	Fig. E.2 Posterior distributions of the fitted parameters for a model with two HG parameters (model 2).
In the text