EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 680 (December 2023) A&A, 680 (2023) A114 Full HTML
Open Access
Issue
A&A
Volume 680, December 2023
Article Number A114
Number of page(s) 26
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/202347353
Published online 19 December 2023

© The Authors 2023

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1 Introduction

In the past 10 yr, the advent of high angular resolution facilities enabled the detection of numerous disk substructures, such as rings, spirals, and dust-depleted cavities in the near-infrared scattered light (e.g., Benisty et al. 2015,2023; Wagner et al. 2018; Shuai et al. 2022) and in the (sub-)millimeter/mm regime (e.g., Francis & van der Marel 2020; Long et al. 2022), indicating the ubiquity of substructures in large, bright disks (Bae et al. 2023). These substructures can be interpreted as evidence of planet-disk interactions, suggesting the presence of an underlying, yet-undetected population of young exoplanets (e.g., Dong et al. 2012). Additional support for this interpretation has recently come from the detection of local velocity deviations in the gaseous outer disk velocity field probed with ALMA (e.g., Pinte et al. 2018, 2020; Teague et al. 2018; Wölfer et al. 2023; Stadler et al. 2023). Scattered light surveys have also pointed out a large fraction of infrared-faint disks, which appear more compact and featureless in scattered light because of self-shadowing effects (e.g., Garufi et al. 2022). However, these disks often host substructures in the sub-millimeter (e.g., Long et al. 2018) that could be attributed to planets.

The presence of massive planets inside cavities was also suggested for transition disks (disks with depleted inner cavities; Bae et al. 2019) and confirmed in the case of at least one system, PDS 70, with the detection of two protoplanets (Keppler et al. 2018; Haffert et al. 2019). The range of plausible mass for the companion(s) in these disks is quite large, as an eccentric stellar companion could be sculpting the cavity (e.g., Calcino et al. 2019) as found in the HD 142527 system (Balmer et al. 2022). In that specific instance, the companion is also leading to a misaligned inner disk, which casts a shadow on the outer disk (Price et al. 2018). Such misalignments were found in at least six transition disks (Bohn et al. 2022). Whether these features are of planetary or stellar nature, the search for perturbers, which are responsible for all the observed disk substructures (e.g., Asensio-Torres et al. 2021; Cugno et al. 2023), is of prime importance to understanding the formation and evolution of planetary systems. The detection of these perturbers would offer crucial observational evidence to test planet-disk interaction theories (e.g., Dong et al. 2015) and constrain the overall evolution of a planetary system (Bae et al. 2019). However, directly imaging planets embedded in bright and highly structured disks is very challenging with current instruments. Until now, all claims but PDS 70 still require confirmation (e.g., Kraus & Ireland 2012; Sallum et al. 2015; Quanz et al. 2015; Reggiani et al. 2018; Wagner et al. 2019, 2023; Boccaletti et al. 2020; Uyama et al. 2020; Currie et al. 2022; Hammond et al. 2023; Law et al. 2023).

To observe exoplanetary systems with high-contrast imaging, observation strategies including angular differential imaging (ADI; Marois et al. 2006, where the parallactic angle diversity of observations is used to remove star light) have enabled the detection of prototypical planetary systems (e.g., HR 8799; Marois et al. 2008). Nevertheless, ADI detections are still limited by self-subtraction at close-in regions from the stars (e.g., Milli et al. 2012; Wahhaj et al. 2021), yet these regions are where giant planets are expected to have the most occurrence (1-10 au; from a combination of radial velocity and high-contrast imaging surveys, e.g., Nielsen et al. 2019; Fulton et al. 2021). To overcome this limitation, on the one hand, better optimized post-processing methods for ADI datasets have been developed (e.g., Pairet et al. 2021; Flasseur et al. 2021; Juillard et al. 2022, 2023). On the other hand, the diversity in archival observational data can enable the usage of other stars as the templates to remove star light and speckles with the reference differential imaging (RDI) data reduction strategy (e.g., Ruane et al. 2019; Xie et al. 2022). Moving forward along the direction of RDI, the Spectro-Polarimetic High contrast imager for Exoplanets REsearch (SPHERE; Beuzit et al. 2019) at the Very Large Telescope (VLT) from European Southern Observatory (ESO) initiated the star-hopping mode (Wahhaj et al. 2021), which offers quasi-simultaneous observations of a science star and its reference star, unleashing the full potential in exoplanet imaging in close-in regions for SPHERE.

Determining dust properties is of fundamental importance for the early stage of grain growth and planetesimal formation, as these parameters will determine the efficiency of grain sticking and fragmentation (Birnstiel et al. 2012). In addition to the planet imaging capabilities with SPHERE, the star-hopping mode enables the optimized extraction of disks in scattered light in terms of total intensity. This goes beyond the polarimetric surveys that have been routinely carried out in the near-infrared (e.g., Avenhaus et al. 2018; Garufi et al. 2020; Ginski et al. 2020) and allows us to better study spatial distribution and properties of dust in the disk (e.g., Olofsson et al. 2023). With the observations taken in dual-polarimetry imaging (DPI; de Boer et al. 2020; van Holstein et al. 2020) mode, which probes polarized signals in the scattered light, star-hopping can also offer total intensity imaging from RDI. The combination of both can yield an estimate of the polarization fraction and, thus, to better constrain dust properties (e.g., shape, composition; Ginski et al. 2023; Tazaki et al. 2023).

In this study, we present the first large survey of proto-planetary disks in total intensity from the ground. As many as 29 young stars are surveyed in the Ks-band with VLT/SPHERE in the star-hopping mode. Our target sample consists of both transition disk systems to search for protoplanets that can potentially reside in the close-in regions with star-hopping, which are otherwise unachievable (Wahhaj et al. 2021), and non-transition disk sample of faint disks in the infrared to search for planets in their outer disk regions. We also aim to derive the polarization fraction whenever possible. The paper is structured as follows: Sect. 2 provides the description of the observations and data reduction procedure, while Sect. 3 presents the polarized light and total intensity maps. Section 4 shows the detection limits of companions and in Sect. 5 we present the polarization fraction maps. We present our summary and conclusions in Sect. 6.

2 Observations and data reduction

2.1 Sample of protoplanetary disks

The sample analyzed in this work includes 29 young stars from both Herbig AeBe and T Tau stars. These are 13 sources from the Taurus star-forming region (CI Tau, CQ Tau, CY Tau, DL Tau, DM Tau, DN Tau, DS Tau, GM Aur, HD 31648, IP Tau, IQ Tau, LkCa 15, and MWC 758), 5 from the Scorpius-Centaurus association (HD 100453, HD 100546, HD 143006, HD 169142, SAO 206462), 3 from Chamaeleon (HD 97048, SZ Cha, SY Cha), 3 from Orion (HD 34282, PDS 201, V1247 Ori), and 1 from each of the following regions: Lupus (V1094 Sco), ϵ Cha (PDS 66), ρ Ophiuchus (EM* SR 20), and Perseus (LkHa 330). There is also one isolated source (HD 163296). To ensure the identifiability of these targets, which often have different names from various database, we also list the SIMBAD identifiers (Wenger et al. 2000) for these targets in Table 1.

To carry out this survey of protoplanetary disks in total intensity, we selected transition disks with known substructures, previously observed in scattered light, as well as non-transition disks in Taurus, with R-band magnitude within the SPHERE limits. Using the Gaia DR3 distances (Gaia Collaboration 2023) and following the approach in Garufi et al. (2018), we uniformly calculated the stellar properties for the entire sample. In particular, we retrieved the effective temperature Teff from VizieR (Ochsenbein et al. 2000) and the photometry using the VizieR photometry tool1. We calculated the stellar luminosity Lstar through a PHOENIX model (Hauschildt et al. 1999) scaled to the de-reddened brightness in V-band. From Teff and Lstar, we derived the interval of values for the stellar age identified by different sets of pre-main-sequence tracks (Siess et al. 2000; Bressan et al. 2012; Baraffe et al. 2015; Choi et al. 2016). We summarize all the derived stellar properties2 in Table 1. Our sample has a wide, uniform coverage of both Teff and Lstar, and therefore of Mstar and the system age, with a stellar mass from 0.5 M⊙ (IQ Tau) to 2.7 M⊙ (LkHa 330) and age from ~1 Myr (a few Taurus sources) to ≳ 10 Myr (e.g., MWC 758).

Table 1

System parameters and star-hopping observation log.

2.2 Observations

We observed 29 disks between November 27 and July 7, 2020, using VLT/SPHERE in the star-hopping mode in Ks -band with the Infra-Red Dual-band Imager and Spectrograph (IRDIS; Dohlen et al. 2008). The data were obtained through six programs: 0103.C-0470, 105.209E, 105.20HV, 106.21HJ, and 108.22EE (PI: M. Benisty), and 105.20JB (PI: M. Keppler).

By observing these systems in Ks-band, we can image the selected protoplanetary disk hosts in the longest wavelength offered by IRDIS3. In this way, we can reach a contrast regime that is more suitable for giant exoplanet imaging (e.g., Spiegel & Burrows 2012; Currie et al. 2023), while simultaneously imaging circumstellar disks and complementing existing IRDIS studies in shorter wavelengths.

The observations were conducted in the DPI mode with pupil tracking, so that we could simultaneously obtain both polarized light and total intensity observations for these systems. In the star-hopping mode, we expect to obtain quasi-simultaneous capture of wavefront variations for an observation pair of a diskhosting star (hereafter “host”) and its corresponding well-chosen point spread function star (which does not host a disk or companion, hereafter “reference”). With star-hopping, we can better capture the stellar speckles for a disk host using the reference–which has a similar color, magnitude, and proximity to the corresponding host – to better reveal circumstellar structures and exoplanets (Wahhaj et al. 2021) than existing archival studies (e.g., Xie et al. 2022). To find the reference stars in our program, we used the SearchCal tool of the JMMC4. See Table 1 for the reference stars, which have identical observational setup as their corresponding host stars, used in our observations.

2.3 Data reduction

We retrieved the raw data in fits format (Pence et al. 2010) for our programs from the ESO archive for SPHERE5. To reduce the IRDIS data, we proceeded as detailed in the following subsections: we first used the IRDAP (van Holstein et al. 2017, van Holstein et al. 2020) package6 for polarimetric differential imaging, which includes both the pre-processing of the raw data in Sect. 2.3.1 and post-processing in polarized light (see Sect. 2.3.2). We then performed reference differential imaging in total intensity light using the IRDAP output files (see Sect. 2.3.3).

2.3.1 Preprocessing

In the star-hopping observation of a host-reference pair, the host has dedicated sky background observations for empirical Ks-band background removal. To remove the stellar signals and speckles in host images, reference stars serve as empirical point spread functions (PSFs; they are occulted by the coronagraph here), yet they do not have dedicated sky backgrounds. To enable sky background removal for a reference star, we used the sky frames of its corresponding host.

We customized the preprocessing procedure in IRDAP to reduce the star-hopping datasets. In addition to sky background removal for the reference, we recalculated the parallactic angle for the target exposures from the PynPoint7 (Stolker et al. 2019) pipeline for SPHERE. We also rescaled the preprocessed rectangular IRDIS pixels to square pixels by streching the pixels by 1.006 along the vertical direction of the detector (e.g., Schmid et al. 2018). The pre-processed images from IRDAP are 1024 × 1024 pixel, where 1 IRDIS pixel is 12.25 mas (Maire et al. 2016), with the stars located at the centers of the images.

2.3.2 Polarimetrie differential imaging (PDI): 𝒬ϕ images

For the polarized data, we performed the PDI data reduction using IRDAP. We ran IRDAP with the default set of parameters to perform the PDI reduction. We used the star-polarization-subtracted 𝒬ϕ data for further analysis (see Fig. 1). We also reduced available archival observations with DPI in other IRDIS bands (i.e., Y, J, or H), as described in Appendix A and Table A.1, to allow for a comparison with the Ks-band data from this study.

2.3.3 Reference diffraction imaging (RDI): Itot images

For total intensity data that can be generated using the polarized observations, we performed RDI data reduction to obtain total intensity (Itot) images of the systems using the data imputation technique described in Ren et al. (2020). For this purpose, we added the IRDAP preprocessed images from the left and the right IRDIS channels, then used the central 350 × 350 pixel for RDI post-processing.

For the reference cube, we obtained their exposure features using non-negative matrix factorization (NMF; Ren et al. 2018) with five sequentially constructed NMF components (i.e., sNMF). In this work, we experimented using 10 or more sNMF components for data interpretation, yet we did not observe clear difference or improvement in the quality of RDI outputs. Thus, we adopted five sNMF components for the rest of the study for computational efficiency. We note that this is due to the high speckle similarity of star-hopping observations, and thus for non-star-hopping observations (e.g., archival data analysis), more sNMF components should be used (e.g., Xie et al. 2023; Sai Krishanth et al. 2023).

To obtain the Itot disk images from the disk host cubes, we imputed stellar signals in a circular region that is within a radius of 80 pixel from the center using the sNMF components. Specifically, we used an annular region that is between 80 pixels and 175 pixels from the stars to model the entire field of view using the NMF data imputation (DI-sNMF) approach described in Ren et al. (2020). We chose the annular region since it both covers the control ring of the adaptive optics (AO) system for SPHERE in Ks-band (the regions can change in different wavelengths), which contains sufficient information to infer the light distribution across the entire field of interest and does not contain disk signals that generally lie within 1″ (i.e., ~82 IRDIS pixels). We then removed the DI-sNMF models for each image in the disk host cube, derotated them to north-up and east-left using previously calculated parallactic angles in Sect. 2.3.1, and adopted the element-wise median as the final disk image in total intensity. We present the total intensity images using RDI from DI-sNMF in Fig. 2.

thumbnail Fig. 1

Ks-band Qϕ maps with dimensions of 2″ × 2″ with different color bars in log scale. The letter identifiers are from Table 1. The rulers correspond to 50 au. The regions interior to 0”. 1 are not accessible with coronagraph usage. We note that the data used to create this figure are available at the CDS.

2.3.4 Polarization fraction and color

With both polarized light and total intensity images, we computed the linear polarization fraction maps. We divided the linearly polarized Qϕ data with PDI from IRDAP by the Itot data with RDI from DI-sNMF. With relatively negligible uncertainty in the Qϕ data from IRDAP, the corresponding uncertainty for the polarization fraction maps are propagated using the element-wise standard deviation map of the Itot results.

To obtain the relative reflectance of the disks, we divided the Qϕ images by the IRDAP measurement of the star. We removed shot noises in the Qϕ data with two-dimensional (2D) Gaussian profile smoothing with standard deviations of σ = 2 pixel (e.g., Olofsson et al. 2018). For observations in non-Ks bands, we additionally convolved the data to reach same spatial resolution as Ks-band, based on Rayleigh diffraction limits of 1.22λ/D where λ is the central wavelength3 and D is the telescope pupil size of 8.0 m. To obtain disk colors in Qϕ, we converted the relative reflectance to magnitudes, then subtracted shorter wavelength magnitudes from longer wavelength ones.

thumbnail Fig. 2

Ks-band Itot maps with dimensions of 2″ × 2″ with different color bars in log scale, for disks identifiable with DI-sNMF star-hopping RDI data reduction. Due to varying observation conditions, a well-chosen reference star may not lead to optimum reduction results, e.g., SZ Cha in panels x and x′. We note that the data used to create this figure are available at the CDS.

3 Disk imaging

We gathered the star-polarization-subtracted Qϕ results from IRDAP, with the images presented in Fig. 1. To our knowledge, the maps of ten objects had not been published before, namely: CY Tau, DL Tau, DM Tau, DN Tau, DS Tau, IP Tau, and IQ Tau, SR 20, SY Cha, and SZ Cha. For systems with a high quality in terms of the total intensity, Itot detections from RDI data reduction using DI-sNMF, we present the images in Fig. 2.

3.1 PDI 𝒬ϕ maps

We obtained disk detection for nearly all of the targets in Fig. 1. For faint disks such as those of CI Tau, CY Tau, DL Tau, DN Tau, DS Tau, IP Tau, IQ Tau, and SR 20, some of the observational data did not produce high-quality detections (e.g., IP Tau) and some are even closer to a non-detection (e.g., SR 20). In comparison with existing PDI data of the systems at shorter wavelengths (i.e., Y, J, or H bands), the Qϕ maps in Ks-band in Fig. 1 do not have significant morphological variations from them. The Ks-band data are less resolved due to an increase in observation wavelength. From the gallery, it is apparent that the M-star companion HD 100453B is polarized in Fig. 1j, and the polarization is also detected in IRDAP-reduced archival data in J and H bands. A polarized HD 100453B indicates that it hosts dust, similarly to CS Cha b (Ginski et al. 2018).

Qϕ signals are expected to primarily trace single scattering events on the surfaces of optically thick disks. The signals are expected to be positive when the polarization vectors are perpendicular to the direction of the incident light on the scatters (e.g., Monnier et al. 2019). Nevertheless, multiple scattering naturally occurs in protoplanetary disks that are optically thick, reducing the polarization fraction of disks (e.g., Tazaki et al. 2019, Fig. 4 therein). In the observations, multiple scattering signals can be revealed in the 𝒰ϕ images (e.g., Canovas et al. 2015; Monnier et al. 2019), which traces the light that is ±45° from the incident light. In addition, due to the finite angular resolution with VLT/SPHERE, the IRDAP Qϕ maps for IRDIS in our study could be lower limits of the actual Qϕ light (i.e., Q^ϕ${{\hat Q}_\phi }$ in Ma et al. 2023) due to convolution effects. Observationally, studies including Engler et al. (2023) demonstrated that the 𝒰ϕ maps (which, in principle, should not contain signals for single-scattering systems such as the HD 114082 debris disk) do contain spurious signals in SPHERE/IRDIS in H-band (central wavelength: ~1.625 μm) but not SPHERE/ZIMPOL in I_PRIM band (~0.790 μm), with the latter having higher angular resolution. For the Ks-band (~2.182 μm) observations in this study, we thus expect the 𝒰ϕ maps contain more leakage from Qϕ data due to a finite angular resolution. We observe that the 𝒰ϕ signals in Fig. A.1 are ≤5% of the Qϕ signals in absolute values, and thus the Qϕ images in this work are not severely impacted (see, e.g., Canovas et al. 2015, for non-negligible impacts in their simulation). To further reduce the limitation on 𝒰ϕ signal leakage, which are albeit of minor contribution in absolute values here, a forward modeling of the convolution effects is necessary (e.g., Engler et al. 2018; Tschudi & Schmid 2021; Ma et al. 2023) and such a solution (e.g., Ma et al. 2023, Sect. 4.2.1 therein) is beyond the scope of this study.

3.2 RDI Itot maps and their fidelity

Using the SPHERE control ring at a region of 80–110 pixel (0″.98-1″.35 from the center) and the region exterior to it in Ks-band (i.e., 80–175 pixel), we were able to recover disks in the RDI total intensity in Fig. 2 for 15 systems using DI-sNMF. While these disks were detected in total intensity, the rest of the disks in our sample were detected with high fidelity primarily only in polarized light (e.g., CI Tau, HD 163296).

To obtain the RDI images, we only used the signals outside the inner edge of the IRDIS control ring of a host (i.e., 80-175 pixel) to infer the speckles interior to the control ring. In other words, we used ~80% of the image to infer the entire image, yet thus such an imputation might fall into the regime of under-fitting. In fact, we also used using only the control ring in Ks-band (i.e., 80–110 pixel) for DI-sNMF reduction and the results do not have significant change: this further supports the importance of the control ring in inferring the PSF signals interior to it using data imputation. In comparison, Ren et al. (2020) showed that the region used for imputation has a second-order influence on the recovery quality of speckle signals. Therefore, with ~20–50% of the regions being masked out in the study here, we would have expected a ~25% change according to Ren et al. (2020), which should have resulted into inferior qualify for the imputed stellar signals and, thus, the recovered disks. To investigate the mathematical reason for the high-quality images in Fig. 2, we advanced the mathematical investigation by presenting a corresponding derivation for ideal imputation cases (i.e., when the Wahhaj et al. 2021 requirements for reference stars are fully satisfied) in Appendix B. In fact, when the individual matrix elements in the matrices are weighted equally in DI-sNMF, the contribution of “missing data” can only have a theoretical fourth-order impact.

With the new derivation in Appendix B showing an expected fourth-order deviation from the missing data, we now can further establish the mathematical background for the authenticity of DI-sNMF reduction, especially when data quality (i.e., speckle stability, speckle resemblance across disk host and reference stars) is ideal. Based on the quality of disk recovery in Fig. 2, we now expect that the masked out region have a ~6% difference even when ~50% of the regions are masked out: this supports the observed high-fidelity morphological similarity between the PDI and RDI results. Moving forward, this indicates that for extended structures that overlap with the control ring or for disk observations in shorter wavelengths where the control ring is angularly closer-in than the Ks-band data here, we do not have to use the entire control ring to recover the extended structures. We leave such an investigation to a future work.

3.3 RDI Itot detectability with stellar parameters

We show in Fig. 3 the detectability of disks in total intensity (the confirmed detections are from Fig. 2), as a function of their stellar Rp – K color and Rp-band (or K-band) magnitude. The Rp-band and K-band magnitudes are from Gaia DR3 (Gaia Collaboration 2023) and 2MASS (Skrutskie et al. 2006), respectively.

Beyond an empirical threshold of Rp ≳ 11 or K ≳ 8, the disks are not detected in total intensity with DI-sNMF even when they are detected in Qϕ. Given that DI-sNMF depends on the adaptive optics’ control ring signals for data reduction in Sect. 3.2, this illustrates that the importance of the adaptive optics performance (e.g., Jones et al. 2022) in producing the control ring for DI-sNMF data reduction.

3.3.1 Logistic regression detectability

To quantify disk detectability using other system parameters, we used R (R Core Team 2022) to perform logistic regression. We explored independent regression in Fig. 4 using the host star magnitude, host-reference color difference, host-reference magnitude difference in Rp- and Ks-band, or angular separation between the target and the reference star. When we fit these parameters independently, we observe that brighter stars can yield better detection, see Fig. 4a. While redder reference and brighter Ks -band reference are more likely to yield better detections (see Figs. 4b,c, respectively), the relations are marginal.

We originally observed that fainter reference in Rp-band or angularly closer-in references might negatively impact detections (see Figs. 4d,e, respectively). However, we argue that they result from selection bias. On the one hand, most of the non-detections with Rp ≳ 11 have brighter references, which bias the fit in Fig. 4d. On the other hand, specifically, the HD 97048 disk is detected in Fig. 2i, while the reference is 3°.80 from it in Fig. 3 (note: that reference star was chosen since there are no other good nearby references for HD 97048). As a result, to select reference stars for star-hopping observations, we do not recommend selecting fainter Rp-band references or angularly distant references. Instead, we recommend focusing on other parameters in Figs. 4a–c: brighter host stars, redder references, and slightly brighter references in observational wavelengths.

When we jointly fit the disk detectability with DI-sNMF in total intensity with these parameters, only the relationship on target star magnitude persist and that trend does not change when we use the Gaia G-band8 or 2MASS Ks-band magnitudes. After all, the reference stars were already selected based on their magnitude and color match, as well as on-sky proximity, the logistic regression here could suffer from severe selection bias that is less severe were the references chosen randomly. However, a random selection of reference stars is dissuaded in star-hopping observations due to observation efficiency.

Due to the high selection bias among the sample in this work, we did not explore the detectability of disks as a function of disk property such as mass and inclination. To perform this study, we first need a proper removal of dominating effects such as host star brightness in Fig. 4a. However, due to the limited RDI disk detections here with DI-sNMF, we do not have enough targets for the exploration on disk properties (see Sect. 4.3 of Ren et al. 2023 for a similar discussion).

3.3.2 Implications for star-hopping reference selection

Assuming the relationships in Figs. 4a–c are trustworthy, we explain the RDI disk detectability as follows. First, brighter disk-hosting stars in apparent light in Fig. 4a can offer more light for scattering in disk particles, increasing disk detectability. Second, the AO system of SPHERE operates in the visible wavelength that overlaps with the Rp- and G-band of Gaia (e.g., Jones et al. 2022, Fig. 5, therein), thus AO operations could potentially offer similar deformable mirror corrections in a star-hopping observation, especially when the reference stars is of similar magnitude as the target star (or slightly brighter to ensure comparable AO performance). Third, assuming there is a pair of host and reference stars with similar magnitudes in AO operation wavelengths, then the redder the reference is (i.e., brighter in Ks-band), the higher the signal-to-noise ratio (S/N) for the reference exposures, and thus the better DI-sNMF data reduction quality. In comparison, for broader-band filters such as the Hubble Space Telescope’s Space Telescope Imaging Spectrograph, which operates at ~0.2–~1.2 μm, using color-matching and slightly brighter references can also yield better host signal recovery (Debes et al. 2019). Following these arguments, the detection of HD 163296 with low fidelity can be explained in Fig. 3: one reference in Rp-band is 1.2 in magnitude fainter than HD 163296 and thus yielded different AO operational status, the other reference (which has similar magnitude and redder) did not produce high-fidelity results, since it was resolved by IRDIS as a binary system during the observation.

There is no clear evidence that on-sky proximity would enhance RDI disk detections. Therefore, star-hopping users should attribute a low priority to it in their reference star selection and instead focus on the above-mentioned parameters. Nonetheless, it is not clear if planet detection capability is impacted by on-sky proximity to the host in star-hopping, especially when there is no circumstellar disk.

For systems with control rings with high signal-to-noise ratios, the fact that the two IRDIS channels pass through different optical paths may yield slight difference in data reduction quality. A potential increase of the disk quality in this work is to perform DI-sNMF reduction for the two channels separately. However, this approach is beyond the current scope of this study, since the purpose of this section is to demonstrate the usage of the control ring in DI-sNMF data reduction. For this purpose, we added the data from the two IRDIS channels, to obtain a factor of ∼2$ \sim \sqrt 2 $ increase in the S/N for the pixels containing control ring signals.

thumbnail Fig. 3

Disks with confident recovery in Ks-band via NMF data imputation, as a function of (a) Gaia DR3 Rp magnitude or (b) 2MASS Ks magnitude with the Rp – K color. Each connected pair is a host-reference pair, with their on-sky angular separation from Gaia DR3 in degrees. Notes: (1) Certain systems with marginal detections in NMF data imputation are marked as non-detection (e.g., HD 163296 or m). (2) The size of reference star symbols reflects typical uncertainties in color-magnitude measurements, zoom in the figure for actual error bars colored yellow for all systems.

4 Companion search

From the RDI results, we did not identify companions except for HD 100453. HD 100453 hosts two prominent spirals and has an M star companion (i.e., HD 100453B, e.g., Wagner et al. 2015; Benisty et al. 2017) and the companion has a Ks-band polarization fraction of 1.0−1.0+1.5%$1.0_{ - 1.0}^{ + 1.5}{\rm{\% }}$, which is measured here using the PDI total polarized intensity and RDI total intensity results. The polarized HD 100453B suggests that there exists circum-secondary disk. See van Holstein et al. (2021) for a study on polarized companions using IRDIS.

In this section, we obtain the detection limits for our star-hopping RDI datasets. We also compare our Ks-band results with previous claims and discuss potential reasons of their non-detection in Ks -band here.

4.1 Detection limits

To quantify the detection capability of point sources, existing surveys (Nielsen et al. 2019; Vigan et al. 2021) adopted primarily the principal-component-analysis-based algorithms (e.g., KLIP; Amara & Quanz 2012; Soummer et al. 2012). However, the results from these methods are prone to contamination, especially when disk signals exist, which further prevent the detection of companions from a combination of algorithm choice (overfitting or over-subtraction) and the ADI observational strategy (self-subtraction). As a result, aggressively post-processed disk signals can resemble point sources (e.g., Rameau et al. 2017; Currie et al. 2019) and extensive optimization is needed to increase the significance level of detections (e.g., Adams Redai et al. 2023).

With the RDI results from DI-sNMF on star-hopping data, we can now better preserve disk signals to avoid them being regarded as of planetary origin. For the RDI data from DI-sNMF here, we do not see significant improvement when more than five sNMF components are used, and thus we use them to calculate the detection limits of companions as a function of angular separation from the star.

thumbnail Fig. 4

Individual logistic regression on RDI detectability of disks in DI-sNMF. The solid dots are for systems in Fig. 2. (a) Brighter host stars can likely yield disk detection. (b) Redder Rp – K band color potentially contribute to disk detection. (c) Brighter reference stars in Ks-band might tentatively contribute to disk detection. (d) Reference stars with fainter Rp magnitude of target stars might tentatively contribute to disk detection, however that is likely biased by the non-detections with Rp ≳ 11, see Fig. 3a. (e) Although reference at farther separation might yield better disk detection, this stems from a selection bias when close-in references are not available. Notes: (1) A joint multinomial logistic regression only preserves the trend for (a). (2) The trends here are subject to change due to the highly selected samples in this study. (3) From solid to dashed lines, the statistical significance decreases; dotted lines are biased trends that are not trustworthy.

4.1.1 Contrast calculation

Using the preprocessed and postprocessed data, we calculated the RDI contrast curves from DI-sNMF for the disks with high-quality detection from Fig. 2. Using the RDI images, for one specific angular separation from the star, we calculated the radial profiles for the median and standard deviation within a 3 pixel annulus. We then rescaled the standard deviation profile by taking into account the small sample statistics in Mawet et al. (2014). Given that these radial profiles are calculated in detector units, we obtained the ratio needed for contrast conversion using stellar counts as follows: using the preprocessed output from IRDAP, we obtain the peak-to-total ratio between the peak and the total stellar photometry from the fits files in the calibration/flux folder (i.e., the flux files and regions used by IRDAP for stellar flux measurement and background removal). We then converted the measured radial profiles to a ratio with the star by dividing them by the peak-to-total ratio.

Noticing the existence of disks, we present in Fig. 5 the 5σ contrast curves by multiplying the rescaled standard deviation profile in units of ratio by 5. For one 5σ contrast curve, we did not add its corresponding median profile which reflects the disk signal to it, since any point source with a brightness of five times the rescaled standard deviation will be super imposed onto the disk, thereby reaching a 5σ contrast. It should be noted that the RDI results from DI-sNMF do not have a mean of zero for each reduction image as KLIP. Instead, the RDI contrast with DI-sNMF is calculated only using the standard deviation due to the existence of disks. To ensure the reliability of the DI-sNMF contrast, we compared our results with the TLOCI (Marois et al. 2014) ones from the High Contrast Data Centre9 and we did not observe any significant differences.

thumbnail Fig. 5

RDI contrast curves (left) and AMES-Cond mass limits (right) of DI-sNMF reductions for systems in Fig. 2 (see Sect. 4 for the detailed calculation procedure). Within ~0″.2, the measured contrast is prone to transmission reduction near the coronagraphic edge. The contrast peak at −1″ for j (HD 100453) is the binary companion. The contrast curve for t (PDS 66) is the deepest due to the faintness of the disk. The data used to create this figure are available at the CDS.

4.1.2 Contrast results

For RDI with DI-sNMF, we reached for Ks-band a 5σ contrast of Δmag ≈ 6–10 at ~0″.2, and Δmag ≈ 11–13 at ~1″. In comparison with Wahhaj et al. (2021), where the detection limits are Δmag ≈ 12 and ~14 correspondingly in K1- or K2-band with SPHERE/IRDIS, the detection limits here are ~2 mag or more brighter. We note that at least two factors have resulted into the difference. First, while the wavelength coverage of the two studies are similar, the filters used in Wahhaj et al. (2021) have narrower wavelength coverage (K1 – and K2-band) than the Ks-band here. With a broader wavelength coverage, stellar speckles on detectors are more extended, thereby limiting companion detection in reaching lower contrasts (e.g., Groff et al. 2016; Desai et al. 2022). Nevertheless, this does not suggest that narrower bands can provide better detections, since a detection is a trade off between the contrast and the planetary luminosity integrated in a band. Second, the existence of bright disks does limit companion detection (e.g., Quiroz et al. 2022) and the asymmetric distribution of disk signals additionally increases the standard deviation in a given annulus for a contrast calculation. In fact, PDS 66 offers the best detection limit in Fig. 5, with Δmag ≈ 10 at −0″.2, which is closer to the Wahhaj et al. (2021) values. This is because that the PDS 66 disk is fainter and more symmetric than other detected disks in total intensity, especially in comparison with other systems with similar Rp- and K-band stellar magnitudes (e.g., h, l, and z in Figs. 2 and 3) and, thus, the impact of disk is smaller than other systems.

To compare the Ks-band data here with other observations, we used AMES-Cond models (Allard et al. 2012) to convert the contrast curves to 5σ upper limits of mass following Wallack et al. (2023). We used the mean and standard deviation of the two age estimates in Table 1 to calculate the mass limit for companions in Fig. 5. We note that we did not account for the extinction from the circumstellar disk, and thus the sensitivity can decrease with the existence of extinction (e.g., Cugno et al. 2023). However, the determination of system ages, together with the evolutionary models and their assumptions (e.g., Allard et al. 2012; Baraffe et al. 2003), could significantly change the detection limits (e.g., Asensio-Torres et al. 2021; Wallack et al. 2023). We note that the limits that we derived are consistent with the ranges of estimates of previous studies (e.g., Asensio-Torres et al. 2021, see Fig. 7 therein). Nevertheless, the star-hopping mode should be preferred for future planet-hunting efforts, since it does not have a requirement on sky rotation for ADI data reduction.

We note that the star-hopping observations in this study were executed in the pupil-tracking mode of SPHERE/IRDIS, for which both pupil-tracking and field-tracking modes are available (e.g., Maire et al. 2021). In pupil-tracking, the diffraction spikes that are evident in field-tracking are suppressed. In field-tracking, while we would expect obtain a sufficient sky rotation for ADI data reduction, the rotation of the diffraction spikes in an observation sequence effectively limit the data reduction quality in reducing star-hopping observations with RDI. This is because the diffraction spikes from a host image cannot be removed using a rotated set of diffraction spikes from a reference image. For RDI data reduction, we thus only recommend star-hopping under pupil-tracking.

4.2 Comparison with existing claims

Several targets in this study were reported to have exoplanet candidates or claims from high-contrast imaging, such as HD 100546 (Quanz et al. 2015), HD 169142 (Hammond et al. 2023), LkCa 15 (Kraus & Ireland 2012; Sallum et al. 2015), and MWC 758 (Reggiani et al. 2018; Wagner et al. 2019, 2023). However, some of them were later identified to be likely disk signals (e.g., HD 100546: Rameau et al. 2017; LkCa 15: Currie et al. 2019) or non-confirmation (e.g., MWC 758; Boccaletti et al. 2021). For the planets that are embedded in disks (e.g., PDS 70 c; Haffert et al. 2019), proper separation of the signal between the planet and the disk is needed to confirm the planetary existence (Wang et al. 2020; Zhou et al. 2023). Such planets that do not host circumplanetary disks are normally expected to be only visible in total intensity Itot images, but not in polarized Qϕ images. While circumplanetary disks might be potentially detectable in polarized light observations (e.g., simulations in Szulágyi & Garufi 2021), there are no confirmed detections with current instruments yet. A detection only in total intensity and not in polarized light would be a direct evidence of a planet that is embedded in a circumstellar disk (e.g., Currie et al. 2022). With the two observational modes in this work, we can investigate the existence of such planets: when a planet is embedded, the polarization fraction value of the region where it resides should be smaller. What is more, we can use the polarized observations to trace the leading and trailing spirals of forming planets (e.g., Hammond et al. 2023).

With the polarization fraction maps in Fig. 6, however, we cannot directly recover any of the existing claims in the data here. This could be due to several factors: first, existing claims are reported in different wavelengths, making them not necessarily visible in Ks-band. Second, even if existing claimed planetary objects are visible in Ks-band, they can be fainter during our observational epochs due to variability (e.g., Sutlieff et al. 2023). Third, even if they are visible, they might have moved behind the coronagraph during our observation (e.g., AF Lep b detected in Franson et al. 2023; De Rosa et al. 2023; Mesa et al. 2023, yet not in Nielsen et al. 2019). Last but not least, with longer wavelengths in Ks-band than in J-/H-band and thus worse spatial resolution (i.e., 1.22 λ/D, see Fig. A.2 for the IRDIS filters), the planetary signal would be apparently more spread out and embedded onto disks due to broadening PSFs, even if the former is not physically co-located with the latter. This PSF-broadening effect further prevents a proper separation of planetary and disk signals, since it spreads planetary signals and make them appear less evident in polarization fraction maps. In principle, using only total intensity observations, here we could use high-pass filters to recover the disk-embedded planets. However, such an approach may require extensive tuning of the filtering parameters, and the highly structured disk morphology in this study further prevents a proper categorization between planetary and localized disk signals.

thumbnail Fig. 6

Ks-band polarization fraction maps with dimensions of 2″ × 2″ with identical color bars in linear scale, obtained from dividing the Qϕ maps in Fig. 1 by the corresponding Itot maps in Fig. 2. The data used to create this figure are available at the CDS.

5 Polarization fraction and color

With the PDI Qϕ and RDI Itot images, we calculated the polarization fraction maps by dividing them. The polarization fraction maps, also known as degree of linear polarization, are presented in Fig. 6. To explore the ensemble properties of the scatterers, we study both their polarization fraction curves, which depict the polarization fraction dependence on the scattering angle, and their colors in polarized light.

5.1 Parametric description

Polarization fraction curves can peak at different scattering angles and be asymmetric about the peak in both observations (e.g., Muñoz et al. 2021, Fig. 6 therein; Kiselev et al. 2022, Fig. 3 therein; Frattin et al. 2022) and theoretical studies (e.g., Tazaki et al. 2019, Figs. 2 and 4 therein; Chen et al. 2020, Fig. 17 therein). Motivated by the expectation that the polarization fraction is zero when the scattering angle is 0° or 180°, we can describe a polarization fraction curve using a scaled beta distribution. For a scattering angle θscat ∊ [0, π], which is the angle between the incident light vector and the scattered light vector, the polarization fraction follows fpol (θscat )∝θscat α−1(π−θscat )β−1,$ {f_{{\rm{pol}}}}\left( {{\theta _{{\rm{scat}}}}} \right) \propto \theta _{{\rm{scat}}}^{\alpha - 1}{\left( {\pi - {\theta _{{\rm{scat}}}}} \right)^{\beta - 1}}, $(1)

where α > 1 and β > 1. See Eq. (C.2) for the full expression using scaled beta distribution in statistics, where we also used the peak polarization fraction parameter, fpolmax$f_{{\rm{pol}}}^{\max }$, to control the maximum polarization fraction.

We implement the analytical polarization fraction curve for analysis. To depict the observed polarization fraction maps, we use the three-dimensional (3D) symmetric and flared disk geometry in diskmap10 (Stolker et al. 2016) by customizing its polarization fraction curve function. Specifically, diskmap can generate a 2D image from a 3D parameterized disk model. Even when some pixels have the same radial separation from the star, they can have different scattering phase angles in a flared disk. Using diskmap, we can convert between a total intensity image and a polarized image using a polarization fraction phase function. To implement this in our analysis, the vertical height of a disk in di skmap in cylindrical coordinates follows h(r)=h0⋅(r1  au)γ,$ h\left( r \right) = {h_0} \cdot {\left( {{r \over {1\,\,{\rm{au}}}}} \right)^\gamma }, $(2)

where h0 is the scale height and γ is the flaring index (see Appendix C.2 for details). For the 3D setup of the disk, we adopted the Bohn et al. (2022) results of the outer disks for both their inclination and position angle of the systems when available, and from Wagner et al. (2020) for PDS 201. In fact, we also fit the two angles independent of the values in previous publications and we did not obtain any significant deviation for the position angle, nor did we have better constraints for the inclination of the systems. Thus, in order to reduce the computational cost with no loss of information, we adopted their published values.

To model the observed polarization fraction maps, we used the Savitzky–Golay filter in two-dimension11 to minimize the random noise in Qϕ data, see Appendix C.5. Specifically, based on the resolution in Ks-band, we used the Savitzky-Golay filter to remove the random noise for the observational data by fitting five-degree polynomials with an 11-pixel window. In this way, we can fit smooth polarization fraction models to them without resolution degradation. Without the Savitzky–Golay smoothing, the best-fit models are prone to shot noise and do not produce consistent results across multiple observations even for a same system (e.g., SZ Cha in Figs. C.1 and C.2).

We explored the parametric dust scattering parameters in Eq. (1) and the peak polarization fraction fpolmax$f_{{\rm{pol}}}^{\max }$, together with geometrical parameters (i.e., scale height h0, flaring index γ) in Eq. (2), using emcee12 (Foreman-Mackey et al. 2013) which performs a Markov chain Monte Carlo (MCMC) exploration. We minimize the residuals by maximizing the following log-likelihood function, lnℒ(Θ∣Xobs)=−12∑i=1N(Xobs,i−Xmodel, iσobs,i)2−∑i=1Nlnσobs, i−N2ln(2π),$ \matrix{{\ln {\cal L}\left( {{\bf{\Theta }}\mid {X_{{\rm{obs}}}}} \right)} = { - {1 \over 2}\sum\limits_{i = 1}^N {{{\left( {{{{X_{{\rm{obs}},i}} - {X_{{\rm{model}},\,i}}} \over {{\sigma _{{\rm{obs}},i}}}}} \right)}^2}} } \cr {}{}{ - \sum\limits_{i = 1}^N {\ln {\sigma _{{\rm{obs}},\,i}}} - {N \over 2}\ln \left( {2\pi } \right),} \cr } $(3)

where Θ denotes the set of parametric scattering and geometrical parameters in Eqs. (1) and (2), while σobs is the uncertainty map for the polarization fraction map. In the above function, we also assume that the pixels, i, follow independent normal distributions. The Xobs and Xmodel parameters denote the observation and model datasets for the polarization fraction maps, respectively.

To obtain the MCMC results, we limited the disk parameters Θ in Eqs. (1) and (2) using uniform priors, with 0 < h0 < 0.2, 0 < γ < 2, 0 ≤ 1, 1 < α ≤ 5 and 1 < β ≤ 5. Using the polarization fraction data from Fig. 6, we present the best-fit polarization fraction curves in Fig. 7 and the corresponding parameters to generate the curves in Table C.1. In Appendix D, we present the images for the corresponding models and residuals. For the α and β parameters of the beta distribution, we set the upper limit to 5 since we do not observe significant curve change when we change the upper limit to larger values, since the internal data variation across different observations (see MWC 758 in Fig. 7) dominates the statistical variations. In addition, the emcee posteriors have extremely narrow ranges for the investigated parameters, and thus we do not present the corresponding ranges in Fig. 7, nor do we present the credible intervals in Tabled.

From an alternative approach, we also experimented forward modeling to obtain the polarization fraction curves. We first generated a total intensity disk model using Qϕ data using diskmap, then removed it from the original observations, and performed RDI reduction using KLIP (Soummer et al. 2012) to minimize the residuals, see Appendix C.3. To reach high the computational efficiency and data quality, we do not adopt the KLIP forward modeling results, and instead present the results from direct polarization fraction modeling in Fig. 7.

With the best-fit 2D polarization fraction models, we notice that a one-component scaled beta distribution polarization fraction model in Eq. (1) cannot fully describe multi-component observations, especially, CQ Tau, HD 100546, HD 143006, PDS 201, for instance (see Appendix D). In fact, one-component models are unable to capture local variations in smaller spatial scales (e.g., MWC 758, SAO 206462, V1247 Ori). Nevertheless, these models are still able to depict the large-scale variation even for spiral systems, especially in reproducing the regions with less polarization (i.e., north-west of MWC 758, south-east of SAO 206462) that might be otherwise categorized as other effects such as shadows.

Detailed inspection of the modeling residuals (e.g., Fig. D.2) suggests that the polarization fraction curves could vary within individual systems. For example, the multiple rings of SZ Cha can have different polarization fraction curves. We focused only on the largest-scale structures for both ring and spiral systems here, to obtain a general understanding for the polarization fraction of protoplanetary disks. Future studies, including modeling the ring components separately for multi-ringed systems, and focused work on spiral systems, are necessary to quantify the difference in scattering properties (e.g., polarization fraction curve) within individual systems.

thumbnail Fig. 7

Ks-band polarization fraction curves, see Appendix C.1 and Table C.1 for the analytical expressions using scaled beta distribution. Best-fit curves with gray segments inaccessible from observation, based on system inclination assuming infinitely thin disks. The error bar on the top right, which could be used to infer the systematic uncertainty from the fitting method, is from the standard deviation of three sets of MWC 758 best-fit results.
thumbnail Fig. 8

Scattering peak values in angle and fraction (SPAF) plot. (a) Tentative correlation between peak polarization fraction value and peak scattering angle, with the gray lines randomly fitting 6 to 12 systems for correlation exploration. The best-fit expression is performed on all data points here with a 0.06 uncertainty for the peak polarization fraction, and certain systems (For b or CQ Tau and k or HD 100546, marked with hollow symbols; for l or HD 143006, a scaled beta distribution cannot describe its polarization fraction curve, potentially due to large scale shadowing in Benisty et al. 2018.) were excluded from Fig. 7 to illustrate this relationship (see Sect. 5.2.1 for a detailed discussion). (b) Gaussian random sphere (GRS) dust models from Tazaki & Dominik (2022) and Tazaki et al. (2023) by varying the minimum dust size for irregular compact grains, overlaid on the 1σ and 2σ ranges from the resampling results in (a). We observe a similar trend for absorptive material (“amc”) but not for less absorptive materials (“org”), with the latter neither providing peak scattering angles with the observations that are consistent beyond ~100° (see Sect. 5.2.2). We note that the dust models do not necessarily reproduce the polarization fraction curves in Fig. 7.

5.2 Interpretation

5.2.1 Empirical trend

From the polarization fraction modeling results in Fig. 7, we observe that the peak polarization fraction is ≲0.6 for the systems in this study. In addition, the peak polarization fraction occurs at ≲90° scattering angle for nearly all systems. For certain systems, we could observe a tentative trend: the peak polarization fraction positively correlates with the peak scattering angle, see the scattering peak values in angle and fraction (SPAF) plot in Fig. 8a, where we assigned an uncertainty of 0.06 in peak polarization fraction for all samples (from a combination of a statistical uncertainty in the residual maps from modeling, with a systematic uncertainty from MWC 758 observations, see Fig. D.2). We did not observe a dependence of peak polarization fraction or its corresponding scattering angle, on system inclination: Fig. 7 is therefore not prone to modeling bias due to system inclination angles.

To generate the SPAF trend in Fig. 8a, we did not include CQ Tau, HD 100546, or HD 143006. The three excluded systems have more than two asymmetric disk components that are mutually superimposed in their polarization fraction maps and the superimposition will lead to non-credible results if they are fitted using a one-component scaled beta distribution. Specifically, CQ Tau is spatially more elongated along the north-east region than the south-west region and, thus, a single scaled beta distribution might not be able to model such a polarization fraction map, yielding high polarization fraction at small scattering angles. Second, HD 100546 has a “bright wedge” in its southwest region (Garufi et al. 2016) and it is also observed in Fig. 6, yielding high polarization fraction at large scattering angles. Third, HD 143006 has two ring components and large-scale self-shadowing effects (Benisty et al. 2018) and we observe that its polarization fraction along the north is significantly higher than that along the south, thus a one-component scaled beta distribution cannot describe its polarization fraction map. After all, if these three systems were included, the trend in Fig. 8a is less evident, and the inclusion of CQ Tau would even make the trend a negative correlation.

With 13 SPAF samples in Fig. 8a, we notice a tentatively positive correlation between peak scattering angle and peak polarization fraction. To investigate the robustness of this SPAF correlation, we first randomly resampled 6–12 systems, then drew samples from normal distribution for the peak polarization fraction values, and performed linear fit. In this way, we can conservatively investigate the correlation for the entire population that are currently inaccessible. With dominating positive SPAF correlations from the resampling, we obtained negative correlation in ~10% of this investigation.

Given that ~10% of the resampled trends have negative correlation, and that we have excluded systems for SPAF trend analysis (i.e., CQ Tau, HD 100546, HD 143006), it is thus possible that positive SPAF trend in Fig. 8a could revert when more samples are available in the future. What is more, although the strong residuals in Fig. D.2 (e.g., MWC 758, PDS 201) are included in the resampling exploration of Fig. 8a, parameterizing the polarization fraction curves with scaled beta distributions could be limited. Moving forward, detailed dust model with different composition and geometry are needed to reproduce the samples in Fig. 8a and the observed polarization fraction curves in Fig. 7, especially when more observations are available, in the future.

5.2.2 Dust model comparison

Polarization fraction could be indicative of dust properties, such as the porosity of aggregates and the size of their constituent grains (i.e., monomers; Tazaki et al. 2023). With the extracted polarization fraction curves, on the one hand, the observed maximum polarization fraction of ≲0.6 for the Ks-band in indicates that the monomers should be bigger than 100 nm (Tazaki et al. 2023, Fig. 6 therein). On the other hand, we can compare the observed polarization fractions with the predictions from radiative transfer numerical studies on different dust morphology.

To compare with numerical predictions, we first used the dust model results in Ks-band from the AggScatVIR database13 (Tazaki & Dominik 2022; Tazaki et al. 2023) and focused on irregular compact grains. We noticed that some of the predicted polarization fraction curves follow skewed bell-shaped profiles resembling the extracted curves in Fig. 7. Following the dust model categorization in Tazaki & Dominik (2022), we adopted the Gaussian random sphere (GRS) results therein, where there are two families of dust particles based on their composition. In each dust family, the constituting dust has a fixed mixture of four compositions (pyroxene silicate, water ice, carbonaceous material, and troilite) from Tazaki & Dominik (2022), see Table 1 therein for the optical constants. From the predictions, we extracted the maximum polarization fraction and its corresponding peak scattering angle for comparison with Fig. 8a.

The GRS model, namely, the irregular compact grain, has a shape characterized by a power-law autocorrelation function with an index v = −3.4 and the relative standard deviation of the radius σ = 0.2 (see Nousiainen et al. 2003, for more detailed descriptions). We note that the GRS model is a single solid grain and, therefore, does not have a porosity or aggregate structure. The grains are assumed to obey a power-law size distribution with an index of −3.5 and the maximum grain radius of 1.6 μm. The minimum grain radius is a parameter of this study. For the two different compositions, on the one hand, in the Fig. 8b predictions, we observe that for absorptive materials (“amc”), when a maximum grain size is fixed to be 1.6 μm, the peak polarization fraction decreases as the minimum dust size increases. This SPAF trend might be consistent with Fig. 8a, with a caveat that the observed data are not representative of the SPAF population. On the other hand, however, such a SPAF trend cannot be reproduced using less absorptive materials (“org”). What is more, there could exist dust that are larger than 1.6 μm that are accessible in Ks-band. After all, the GRS model might not be a representative description for dust in protoplanetary disks and more careful modeling is needed to explain the observed polarization fraction for the systems in this study.

To explore beyond the GRS models, using the entire AggScatVIR database, which also includes different dust properties (e.g., fractal and compact aggregates) at various porosity levels, we noticed that the peak value decreases with either increasing the dust radius or decreasing porosity. To explain the relatively low polarization fractions, the dust radius and porosity can be degenerate in reproducing certain polarization fraction curves, indicating the potential diversity of scatterers in these systems. However, we emphasize that these models roughly only reproduce the peak polarization fraction dependence on peak scattering angle, but do not match the extracted individual polarization fraction curves. Specifically, the predictions could have multiple local maxima in the polarization fraction curves, yet the beta distribution can only allow one: this mismatch is a limitation for our parametric description, yet it is hidden in the large uncertainties and more complicated parameterization is needed to describe the polarization fraction curves from the AggScatVIR database. To explain the observed polarization fraction maps, first, detailed dust model with different composition and geometry, as well as modeling multiple scattering effects, are needed. Second, observationally, separating the contributions between dust surface density and dust scattering properties could reduce the degeneracy. Third but not least, adopting parametric polarization fraction curves beyond beta distribution would allow multiple local maxima that are suggested in numerical models.

The dust polarization inferred from our observations may differ from those seen in the IM Lup disk surface. Tazaki et al. (2023) found that fractal aggregates having a fractal dimension of 1.5 (i.e., dust mass m∝ac1.5$m \propto a_c^{1.5}$ with ac being the characteristic radius of an aggregate) with monomer size amon = 0.2 μm in the IM Lup disk surface when observed in H-band. In comparison, their best-fitting aggregate model would suggest fpolmax=0.83$f_{{\rm{pol}}}^{\max } = 0.83$ and θmax = 89° at the Ks-band. However, none of our disk samples show such a high level of polarization fraction. Fractal aggregates are naturally formed through hit-and-stick coagulation, which is expected to occur during the early phases of dust coagulation. However, these fractal aggregates may not be long-lived as they quickly grow into larger aggregates and settle into the midplane (Dullemond & Dominik 2005; Tanaka et al. 2005). Eventually, the surface of the disk is dominated by particles replenished through collisional fragmentation of dust particles. What is more, the IM Lup disk is a young Class II disk with an estimated age of ~1.1 Myr (Avenhaus et al. 2018), whereas our samples used to derive the polarization fraction are generally older than IM Lup. The observed differences may reflect different stages of collisional dust evolution in the disks. In fact, using physical and chemical modeling, Cleeves et al. (2016) found that even mm-sized dust particles in IM Lup are lofted on to disk surface. Therefore, the distribution of dust particles on disk surface in IM Lup is likely not representative of those in our targets here.

Studying individual polarization fraction curves may provide dust information from an experimental approach. While the scaled beta distributions given here cannot describe the negative polarization fraction for large scattering angles (small phase angles: e.g., Muñoz et al. 2021; Frattin et al. 2022) in experimental studies, these angles are less accessible due to the inclination of the disks (shown in Fig. 7). With laboratory measurements showing diverse polarization fraction curves (e.g., Muñoz et al. 2021; Frattin et al. 2022), we will have access to high-quality phase curves for comparison with the observed ones.

thumbnail Fig. 9

Polarized color at ≈90° scattering angle and stellar luminosity in Sect. 5.3. (a), (b), and (c) are the Y, J, and H-band data in polarized light in comparison with Ks-band data in polarized light, respectively. In J – Ks comparison, increase in stellar luminosity leads to more neutral color; while the correlation is less evidence in H – Ks comparison likely due to adjacent wavelengths (e.g., Fig. A.2). Note: the bands are 1σ, 2σ, and 3σ confidence intervals from bootstrapping fit. The data used to create this figure are available at the CDS.

5.3 Disk color in polarized light

From the color images, calculated by comparing Y -, J-, or H-band data with Ks band data, we measured the colors at ~90° scattering angle following Ren et al. (2023). We present in Fig. 9 the color dependence on stellar luminosity. Such a dependence can reflect the dust properties of the scatterers in these systems as well as limitations in the observations.

We observe in Fig. 9b that the observed protoplanetary disks are blue in J and Ks bands for stars that are less luminous than ~10 L⊙. In H and Ks bands, the colors of the disks can vary between blue and red. When stellar luminosity increases, we observe that both the J – Ks and the H – Ks color in polarized light are redder, and that the H – Ks polarized color changes slower than that of J – Ks. The slow H – Ks polarized color change is possibly caused by the adjacency of the two wavelengths in SPHERE/IRDIS, see Fig. A.2. In comparison, Crotts et al. (2023) showed that debris disks in polarized light do also transition to redder color using the Gemini Planet Imager. The relatively redder color when stellar luminosity increases indicates that the scatterers are larger.

Measured at a ~90° scattering angle, existing color studies in scattered light by Ren et al. (2023) for debris disks showed a ubiquitously blue color in total intensity light by comparing between visible (~0.6 μm) and the near-IR (~1.1 μm and ~ 1.6 μm). With both studies having blue colors, we however note that the two studies are not comparable. The blue debris disks are between visible and near-infrared wavelengths (the latter is close to J/H bands) in total intensity in Ren et al. (2023), while the protoplanetary disks included here are between the J- or H-bands and the Ks-band in polarized light.

6 Summary

We obtained Ks -band imaging of protoplanetary disks in scattered light using SPHERE/IRDIS on VLT for 29 systems in star-hopping mode. In the DPI setup of IRDIS imaging, we can obtain both polarized light observations and total intensity observations simultaneously.

By modeling the interior regions of the IRDIS Ks-band control ring using the information on the control ring with DI-sNMF, we have identified 15 systems in total intensity light with unprecedented data quality. For the RDI results from DI-sNMF, we calculated the companion detection limits for these observations with high-quality disk recovery: the existence of disks does raise the Ks-band detection limits in comparison to the exploration in K1-/K2-bands in Wahhaj et al. (2021). Nevertheless, an actual detection is a trade-off between contrast and band-integrated companion luminosity, and thus narrower bands do not necessarily always provide better detections. Given that star-hopping observation has no dependence on sky rotation in the pupil-tracking mode and that it can reach similar mass detection limits as ADI observations, it should be preferred to ADI observations in terms of observational schedulability.

Together with the IRDIS Qϕ data, we obtained the polarization fraction maps for these systems. With these polarization fraction maps, we can reduce the confusion by blob structures resembling planetary signals, since signals from giant protoplanets are not expected to be polarized. For the polarization fraction maps, we described the polarization fraction curves using analytical beta distributions. The polarization fractions peak between ~20% and ~50%, yet they could be smaller than the actual values due to convolution effects from instrumentation. Assuming these polarization fraction curves are a credible representation of the actual polarization fractions or they otherwise undergo similar convolution effects, then we observe a tentative trend: the peak polarization fraction increases with the peak scattering angle. Using the Tazaki & Dominik (2022) and Tazaki et al. (2023) dustmodels from the AggScatVIR database, we could reproduce such a trend using absorptive materials for GRS dust; nevertheless, such models do not produce the individual polarization fraction curves. In addition, there can be alternative explanations with different dust parameters and future analyses and dust modeling are needed to interpret the observed polarization fraction curves. Moving forward, a more comprehensive extraction of the polarization fraction curves (including a modeling of the disk components separately) can aid in comparing the scattering properties within each disk. In addition, lab measurements (e.g., Muñoz et al. 2021; Frattin et al. 2022) may provide important dust information for the observed polarization fraction curves.

For the 26 systems that have existing IRDIS observations in shorter wavelengths (Y-, J-, or H-band), we obtained the color of these systems at ~90° scattering angle in polarized light. For Jpol – Ks pol and Hpol – Ks pol color in polarized light, we observe trends that the color is relatively redder when stellar luminosity increases. Such a trend indicates that the scatterers are larger for more luminous stars (e.g., Ren et al. 2023; Crotts et al. 2023). In addition, while the polarized H – Ks color here has a marginal trend of being relatively redder as stellar luminosity increases, the color ranges from red to blue for systems similar stellar luminosity, demonstrating the diversity of scatterers in different systems. In order to obtain the properties of the scatterers (e.g., mineralogy, morphology, porosity, size), detailed radiative transfer modeling efforts based on realistic models (e.g., Tazaki & Dominik 2022; Tazaki et al. 2023) are needed.

Using the SPHERE/IRDIS control ring for RDI data reduction with DI-sNMF, we cannot yet recover the disks in total intensity for systems with Gaia DR3 Rp ≳ 11 or 2MASS K ≳ 8. For the sample with high selection bias here, our logistic regression results indicate that brighter hosts, redder references, and brighter references in observational wavelengths could aid in detecting disks. Given that there is no clear evidence that closer-in references can provide better RDI imagery for the hosts, star-hopping users can attribute a lower priority to on-sky proximity in the reference selection.

Acknowledgements

We thank the anonymous referee for their prompt and constructive comments. We thank Valentin Christiaens for comments on the manuscript. B.B.R. thanks Yinzi Xin for discussions on wavefront sensing in high-contrast imaging, Jie Ma on convolution effects, and Laurent Pueyo for support. Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programs 0103.C-0470, 105.209E, 105.20HV, 105.20JB, 106.21HJ, and 108.22EE. For the archival data in Sect. 5.3, based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programs 60.A-9389, 60.A-9800, 095.C-0273, 096.C-0248, 096.C-0523, 097.C-0523, 097.C-0702, 097.C-0902, 297.C-5023, 198.C-0209, 098.C-0486, 098.C-0760, 099.C-0147, 0100.C-0452, 0100.C-0647, 0101.C-0464, 0101.C-0867, 0102.C-0162, 0102.C-0453, 0102.C-0778, 1104.C-0415, 0104.C-0472, 0104.C-0850, 109.23BC, and 111.24GG. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (PROTOPLANETS, grant agreement no. 101002188). This project has received funding from the European Union’s Horizon Europe research and innovation programme under the Marie Sklodowska-Curie grant agreement no. 101103114. This work has made use of the High Contrast 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 SIMBAD database (Wenger et al. 2000), operated at CDS, Strasbourg, France. This research has made use of the VizieR catalogue access tool, CDS, Strasbourg, France (DOI: 10.26093/cds/vizier). The original description of the VizieR service was published in Ochsenbein et al. (2000). The VizieR photometry tool is developed by Anne-Camille Simon and Thomas Boch. This research has made use of the Jean-Marie Mariotti Center SearchCal service (http://www.jmmc.fr/searchcal) co-developed by LAGRANGE and IPAG.

Appendix A Auxiliary IRDIS imaging data

We present the stellar-signal-removed 𝒰ϕ images from IRDAP in Ks-band in Fig. A.1. The absolute values of the 𝒰ϕ signals are ≲5% of the Qϕ signals and, thus, the Qϕ images dominate the polarization signals for the protoplanetary disks in this study.

To study the Qϕ polarized color for the protoplanetary disks in this work, we summarize available SPHERE/IRDIS Y-, J-, and H-band observations in broadband polarized light (in Table A.1) and compared them with the Ks-band data form this study. For HD 100546 and HD 163296 in the J-band, we obtained the data from program 111.24GG and 109.23BC, respectively. We show in Fig. A.2 the transmission profiles and central wavelengths3 for the IRDIS filters with the data analyzed in this study.

Table A.1

Archival SPHERE observations of Ks-band counterparts in the Y, J, or H bands in polarized light

thumbnail Fig. A.1

Ks-band 𝒰ϕ maps with dimensions of 2″ × 2″ with different color bars in linear scale. The field of view of each panel corresponds to those in Fig. 1. The data used to create this figure are available at the CDS.
thumbnail Fig. A.2

Transmission profiles, as well as the central wavelengths, for SPHERE/IRDIS in the Y, J, H, and Ks bands.

Appendix B Data deviation from NMF imputation

Missing data can impact the minimization of the cost function for matrix decomposition and dimensionality reduction methods. For NMF, Ren et al. (2020) showed that the expected deviation due to missing data could follow a second-order form in their Eq. (33). Specifically, given a target image T∈ℝ≥01×Npix $T \in _{ \ge 0}^{1 \times {N_{{\rm{pix}}}}}$ with Npix∈ℕ${N_{{\rm{pix}}}} \in $ pixels, with an NMF component basis vector Hi∈ℝ≥01×Npix ${H_i} \in _{ \ge 0}^{1 \times {N_{{\rm{pix}}}}}$, we can denote the corresponding coefficient with ωi∈ℝ≥0${\omega _i} \in {_{ \ge 0}}$. When a fraction of the data in T is missing (or artificially ignored here), the corresponding coefficient is ωi′∈ℝ≥0$\omega _i^\prime \in {_{ \ge 0}}$. With these notations, Theorem 2 in Ren et al. (2020) states that | ωi−ωi′ |=o2(ωi)=ωi⋅o2(1),$ \left| {{\omega _i} - \omega _i^\prime } \right| = {o^2}\left( {{\omega _i}} \right) = {\omega _i} \cdot {o^2}(1), $(B.1)

where o is the little-o notation, meaning |o(x)| ≪ |x|. Here we provide a derivation of the deviation under a more ideal assumption.

Theorem 3 (ideal imputation): In the case of missing data, if the cross-talk among NMF components is of the same order as the target modeling procedure, the influence of the missing data can reach a fourth-order deviation.

Proof : The second order deviation for Theorem 2 in Ren et al. (2020) showed in their Eq. 27 that: ωi=THiTHiHiT⋅(1+∑j=1,j≠inωjωi⋅HjHiTHiHiT)−1$ {\omega _i} = {{TH_i^T} \over {{H_i}H_i^T}} \cdot {\left( {1 + \sum\limits_{j = 1,j \ne i}^n {{{{\omega _j}} \over {{\omega _i}}} \cdot {{{H_j}H_i^T} \over {{H_i}H_i^T}}} } \right)^{ - 1}} $(B.2)

displays a second-order deviation for the multiplicand and the multiplier on the right-hand side. Therefore, the multiplication of two second-order terms would result into a second-order deviation.

Focusing on the second term in the summand above, and following the same derivation process as Eq. (33) in Ren et al. (2020), we can obtain a second-orderdeviation forthe cross-talk terms of H in the second multiplier without loss of generality. Similarly to Eq. (B.1), we have | HjHiTHiHiT−(Hj∘𝟙T)(Hi∘𝟙T)(Hi∘𝟙T)(Hi∘𝟙T)T |=o2(HjHiTHiHiT),$ \left| {{{{H_j}H_i^T} \over {{H_i}H_i^T}} - {{\left( {{H_j} \circ {{\rm{\alpha }}_T}} \right)\left( {{H_i} \circ {{\rm{\alpha }}_T}} \right)} \over {\left( {{H_i} \circ {{\rm{\alpha }}_T}} \right){{\left( {{H_i} \circ {{\rm{\alpha }}_T}} \right)}^T}}}} \right| = {o^2}\left( {{{{H_j}H_i^T} \over {{H_i}H_i^T}}} \right), $(B.3)

where 𝟙T∈B1×Npix${1_T} \in {^{1 \times {N_{{\rm{pix}}}}}}$ is an indicator matrix (in which 𝟙Tj = 0 when the corresponding element in Tj is missing, and 𝟙Tj = 1 otherwise) which matches the dimension of T and Hi.

Comparing the target modeling terms in Eq 33 of Ren et al (2020), | THiTHiHiT−(T∘𝟙T)HiT(Hi∘𝟙T)(Hi∘𝟙T)T |=o2(THiTHiHiT),$ \left| {{{TH_i^T} \over {{H_i}H_i^T}} - {{\left( {T \circ {1_T}} \right)H_i^T} \over {\left( {{H_i} \circ {1_T}} \right){{\left( {{H_i} \circ {1_T}} \right)}^T}}}} \right| = {o^2}\left( {{{TH_i^T} \over {{H_i}H_i^T}}} \right), $(B.4)

with the component cross-talk term in Eq (B.3), if the target modeling term, THiTHiHiT${{TH_i^T} \over {{H_i}H_i^T}}$, and the cross-talk term among the components, HjHiTHiHiT∀j≠i${{{H_j}H_i^T} \over {{H_i}H_i^T}}\forall j \ne i$, have identical orders of magnitude, we can rewrite Eq. (B.2) with the missing data: ωi′=THiTHiHiT[ 1−o2(1) ]⋅(1+∑j=1,j≠inωjωi{ HjHiTHiHiT[ 1−o2(1) ] })−1,$ \omega _i^\prime = {{TH_i^T} \over {{H_i}H_i^T}}\left[ {1 - {o^2}\left( 1 \right)} \right] \cdot {\left( {1 + \sum\limits_{j = 1,j \ne i}^n {{{{\omega _j}} \over {{\omega _i}}}\left\{ {{{{H_j}H_i^T} \over {{H_i}H_i^T}}\left[ {1 - {o^2}\left( 1 \right)} \right]} \right\}} } \right)^{ - 1}}, $(B.5) =THiTHiHiT[ 1−o2(1) ]⋅(1+∑j=1,j≠inωjωiHjHiTHiHiT)−1[ 1+o2(1) ],$ = {{TH_i^T} \over {{H_i}H_i^T}}\left[ {1 - {o^2}\left( 1 \right)} \right] \cdot {\left( {1 + \sum\limits_{j = 1,j \ne i}^n {{{{\omega _j}} \over {{\omega _i}}}{{{H_j}H_i^T} \over {{H_i}H_i^T}}} } \right)^{ - 1}}\left[ {1 + {o^2}\left( 1 \right)} \right], $(B.6) =THiTHiHiT⋅(1+∑j=1,j≠inωjωiHjHiTHiHiT)−1⋅[ 1−o4(1) ]=ωi⋅[ 1−o4(1) ],$ \matrix{ { = {{TH_i^T} \over {{H_i}H_i^T}} \cdot {{\left( {1 + \sum\limits_{j = 1,j \ne i}^n {{{{\omega _j}} \over {{\omega _i}}}{{{H_j}H_i^T} \over {{H_i}H_i^T}}} } \right)}^{ - 1}} \cdot \left[ {1 - {o^4}\left( 1 \right)} \right]} \hfill \cr { = {\omega _i} \cdot \left[ {1 - {o^4}\left( 1 \right)} \right],} \hfill \cr } $(B.7)

which is a fourth-order deviation under small number approximation. To reach Eq. (B.7), we applied the assumption of similar orders (i.e., ideal imputation requirement) between Eq. (B.5) and Eq. (B.6).

In reality, the ideal imputation condition is not always guaranteed since the matrix elements are not equally contributing to the calculation. Therefore, following the same argument as Eq. 23 of Ren et al. (2020), the deviation introduced by missing data is between second and fourth order. The fourth-order deviation in Eq. (B.7) can be observationally approached in reality, as has been supported in this study in Fig. 2. Furthermore, the recovery of the protoplanetary disks in total intensity from our study demonstrates that the control ring of SPHERE’s adaptive optics system in an exposure is well-correlated with the interior PSF (e.g., Guyon et al. 2021, Fig. 2 therein). This demonstrates that we can use the control ring to infer the PSF interior to it, thus strategically avoiding the problem of overfitting that has been plaguing the post-processing of high-contrast imaging observations in total intensity.

In this study, the usage of the SPHERE control ring has yielded beyond state-of-the-art results for the majority of the detected systems. Nevertheless, it is still limited by not only the existence of control rings, but also the objects of interest not superimposed on the control rings. We leave the handling of these limitations forfuture engineering (e.g., Guyon et al. 2021) and methodological studies for potential joint works.

Appendix C Parametric polarization fraction

To convert polarized light observations of disk-only signals to total intensity, existing studies have adopted a bell-shaped polarization curve (e.g., Engler et al. 2017; Olofsson et al. 2018; Lawson et al. 2022), with the curve physically motivated under the Rayleigh polarization regime. Specifically, by dividing the stellar-signal-removed local Stokes Qϕ data from IRDAP (van Holstein et al. 2017, van Holstein et al. 2020) by a polarization fraction map, it is possible to convert polarized data to expected total intensity data. Combining this with a physically flared 3D disk geometry (e.g., diskmap: Stolker et al. 2016), we should, in principle, obtain a well-described total intensity image from polarized light observations.

To explore beyond the limitations from Rayleigh scattering (e.g., Ren et al. 2021), which is nevertheless valid only when dust particles are smaller than observation wavelength by more than one order of magnitude and the peak polarization is at 90° scattering angle, here we adopted a parametric approach for extracting the best-fit polarization fraction curve. We obtained the polarization fraction by comparing IRDAP Qϕ data with total intensity data using diskmap while adopting an axisymmetric geometry for a flared disk.

C.1 Polarization fraction curve: scaled beta distribution

To address the fact that polarization fraction curves do not have to be symmetric around or peak at a scattering angle of π2${\pi \over 2}$ (e.g., Figure 17 of Chen et al. 2020), here we adopted a parametric description of polarization fraction. For a scattering angle θscat ∈ [0, π], the polarization fraction in Eq. (1) is: fpol (θscat)∝θscat α−1(π−θscat )β−1.${f_{{\rm{pol}}}}\left( {{\theta _{{\rm{scat}}}}} \right) \propto \theta _{{\rm{scat}}}^{\alpha - 1}{\left( {\pi - {\theta _{{\rm{scat}}}}} \right)^{\beta - 1}}.$

In terms of statistics, the probability density function (PDF) of a beta distribution follows: B(x∣α,β)=Γ(α+β)Γ(α)Γ(β)xα−1(1−x)β−1,$ B\left( {x\mid \alpha ,\beta } \right) = {{{\rm{\Gamma }}\left( {\alpha + \beta } \right)} \over {{\rm{\Gamma }}\left( \alpha \right){\rm{\Gamma }}\left( \beta \right)}}{x^{\alpha - 1}}{\left( {1 - x} \right)^{\beta - 1}}, $(C.1)

for x ∈ [0,1], and Γ(·) is the gamma function with Γ(x)=∫0∞txe−tdt$\Gamma \left( x \right) = \int_0^\infty {{t^x}{e^{ - t}}} dt$ for ∀x ∈ ℝ+. We can normalize Eq. (1) to have θScat π${{{\theta _{{\rm{Scat}}}}} \over \pi }$ follow a beta distribution form, fpol (θscat|α,β)=1πα+β−2Γ(α+β)Γ(α)Γ(β)θscat α−1(π−θscat)β−1.${f_{{\rm{pol}}}}\left( {\left. {{\theta _{{\rm{scat}}}}} \right|\alpha ,\beta } \right) = {1 \over {{\pi ^{\alpha + \beta - 2}}}}{{{\rm{\Gamma }}\left( {\alpha + \beta } \right)} \over {{\rm{\Gamma }}\left( \alpha \right){\rm{\Gamma }}\left( \beta \right)}}\theta _{{\rm{scat}}}^{\alpha - 1}{\left( {\pi - {\theta _{{\rm{scat}}}}} \right)^{\beta - 1}}.$

To enable observational description of polarization fraction, we can set the maximum polarization fraction to be fpolmax∈[0,1]$f_{{\rm{pol}}}^{\max } \in [0,1]$. We now have a polarization fraction curve of fpol(θscat∣α,β,fpolmax)=fpolmax⋅1B(α−1α+β−2|α,β)⋅B(θscatπ|α,β),$ {f_{{\rm{pol}}}}\left( {{\theta _{{\rm{scat}}}}\mid \alpha ,\beta ,f_{{\rm{pol}}}^{\max }} \right) = f_{{\rm{pol}}}^{\max } \cdot {1 \over {B\left( {\left. {{{\alpha - 1} \over {\alpha + \beta - 2}}} \right|\alpha ,\beta } \right)}} \cdot B\left( {\left. {{{{\theta _{{\rm{scat}}}}} \over \pi }} \right|\alpha ,\beta } \right), $(C.2)

where B(x | α,β) is the original beta distribution PDF evaluated at x using Eq. (C.1). We use this parametric description of the polarzation fraction curves in this study.

In the polarization fraction curve in Eq. (C.2), its second multiplier is the inverse of the beta distribution PDF evaluated at its mode of α−1α+β−2${{\alpha - 1} \over {\alpha + \beta - 2}}$ (i.e., where the polarization fraction peaks), and thus it is used to regulate the maximum polarization fraction to be fpolmax$f_{{\rm{pol}}}^{\max } $ together the first multiplier. The α and β parameters also control the spread of the phase function, in the sense that the variance of a beta PDF is αβ(α+β)2(α+β+1)${{\alpha \beta } \over {{{(\alpha + \beta )}^2}(\alpha + \beta + 1)}}$, or (8α + 4)−1 when β = α. In this study, we have α > 1 and β > 1 to avoid mathematical divergence of the polarization fraction at θscat ∈ {0, π}.

C.2 Implementation: 3D geometry

To generate a polarization fraction map, diskmap needs the specification of the scale height, maximum polarization fraction, position angle, and inclination angle of the disk. For the disk scale height, we have: h(r)=h0⋅(r1au)γ,$ h\left( r \right) = {h_0} \cdot {\left( {{r \over {1{\rm{au}}}}} \right)^\gamma }, $(C.3)

where r ∈ ℝ+ is the stellocentric distance in the midplane of the disk, h0 ∈ ℝ+ is the disk scale height at 1 au, and γ ∈ ℝ+ describes the flaring of the disk.

For the position angle and inclination angle values of the systems, we adopted the outer disk information from Bohn et al. (2022) when these information are available therein. For the polarization fraction function, we use the parametric description in Eq. (C.2). To extract the polarization fraction curves, we used emcee (Foreman-Mackey et al. 2013) to explore the parameters in Eqs. (C.2) and (C.3): Θ={ h0,γ,α,β,fpolmax },$ {\rm{\Theta }} = \left\{ {{h_0},\gamma ,\alpha ,\beta ,f_{{\rm{pol}}}^{\max }} \right\}, $(C.4)

which are related with disk polarization to generate total intensity disk images using polarized images.

We can obtain the polarization fraction curve in direct polarization fraction map comparison or forward modeling. On the one hand, from a direct measurement approach, we use the data imputation results and directly compare them with the polarization fraction map models. On the other hand, from a forward modeling approach, for a given set of parameters, we subtracted the corresponding total intensity model from the preprocessed data, then performed a Karhunen-Loève image projection (KLIP; Soummer et al. 2012; Amara & Quanz 2012) data reduction. For both approaches, we distribute the calculations using the DebrisDiskFM (Ren et al. 2019) framework to reduce real-time cost of parameter exploration on a computer cluster. We minimized the residuals to obtain the best-fit parameters while assuming the pixels are independent from each other. We present the best-fit profiles in Fig. 7 from the direct measurement approach, with the corresponding values in Table C.1.

Table C.1

Best-fit scaled beta distribution description for polarization fraction curve for Fig. 6

C.3 Experiment: KLIP forward modeling

While we adopted the direct polarization fraction map modeling using the DI-sNMF results, KLIP has been the classical postprocessing method in the high-contrast imaging of circumstellar structures. To study the application of scaled beta distribution polarization curve to KLIP, we also investigated the forward modeling approach to extract polarization fraction. Given that relatively simple geometry including ring structures can inform the 3D structures of protoplanetary disks in a straightforward way (e.g., Ginski et al. 2016; de Boer et al. 2016), we first applied the approach to ring systems. We then explored the applicability of the approach to spirals, see Figure C.1(a) for the results for both morphologies and Figure C.2(a) for the best-fit polarization curves.

thumbnail Fig. C.1

Qϕ to total intensity conversion using polarization fraction with scaled beta distribution, (a): direct polarization map division, (b): Savitzky–Golay filtered Qϕ data. From top to bottom, A: LkCa 15, on 2020 December 8; B1: SZ Cha, on 2020 December 29; B2: SZ Cha, on 2020 December 30; C: V1247 Ori, on 2020 December 24. From left to right, 1: Qϕ data; 2: total intensity data from data imputation; 3: total intensity data converted from Qϕ data using polarization fraction maps; 4: RDI KLIP residuals using total intensity models; 5: polarization maps used to model total intensity data with RDI KLIP. A comparison between (a) and (b) shows an improvement of retrieving quality with the Savitzky–Golay filter. Nevertheless, the patterned residuals show that a single profile is limited in describing polarization fraction maps, especially when multiple disk components exist: modeling for separate disk components (e.g., rings in LkCa 15 and SZ Cha) is needed for authentic description of polarization fraction curves.
thumbnail Fig. C.2

Best-fit polarization curves assuming beta distribution from RDI KLIP forward modeling in Figure C.1, the peak polarization locations are marked with circles. The light gray curves denote scattering angle ranges that are not accessible for inclined thin disks with no flaring. (a) Direct conversion: the extracted profiles are different even for the same system (i.e., B1 and B2). (b) Conversion with Savitzky-Golay filtered Qϕ data: the extracted curves are less distinct, since such reductions are less prone to pixel-wise discrete noise.

C.3.1 Rings: LkCa 15 and SZ Cha

Using KLIP forward modeling, we present the results of LkCa 15 observed on UT 2020-12-08, SZ Cha on UT 2020-12-29, and SZ Cha on 2020-12-30. To produce Fig. C.1, we removed a total intensity model from the observations, then performed KLIP data reduction to compare the results.

LkCa 15 hosts a two-ringed structure in scattered light (e.g., Thalmann et al. 2016). By directly dividing the Qϕ data by a polarization fraction map, the best-fit total intensity model resembles qualitatively the total intensity observation with data imputation. Nevertheless, the residual map shows pixelated islands, since the original Qϕ data has shot noise which does not fade away even after field rotation. In addition, although the outermost ring has a moderate level of residuals in Figure C.1(a), strong residuals in the innermost ring suggests that the two rings have distinct polarization curves.

SZ Cha hosts a three-ringed structure in scattered light. Similar as the two-ringed LkCa 15 disk, there exists pixelated residuals. With three rings in the system, the residuals are more evident, suggesting that the KLIP forward modeling approach is not able to depict the system with a simple flared disk while assuming an identical polarization fraction curve. What is more, the recovered polarization fraction maps are at different brightness levels with different peak scattering angles in Figure C.2(a). The pixelated residuals and excess disk residuals for the classical KLIP forward modeling approach suggests that the unmodified Qϕ data cannot be directly applied to multi-ring systems with moderate inclinations.

C.3.2 Spirals: V1247 Ori

For spirals, we applied the KLIP forward modeling approach to the V1247 Ori data on 2020 December 24. Hosting at least a pair of spirals in scattered light (e.g., Ohta et al. 2016), the residual disk signals previous seen in ring systems are stronger for V1247 Ori. With the unideal experiments on ring systems, the unideal performance for the spiral system have been anticipated for the KLIP forward modeling approach.

C.4 Limitations in polarization to total intensity conversion

With the experiments above, in addition to existing evidences including the non-detection of polarized light in the most backward scattering regions for the HR 4796A system (e.g., polarized light: Perrin et al. 2015; total intensity: Milli et al. 2017; Ren et al. 2020), we observe that the joint effect from scattering phase function and polarization fraction can result in a non-detection of polarized signals for certain regions that host total intensity signals. However, the non-detection of such polarized signals does not specifically exclude their existence, namely: such signals might exist, yet they are beyond the sensitivity limits of the existing instruments for given exposure times.

Given the facts that dust properties can vary as a function of stellocentric radius, which less efficient backward scattering could redistribute less light in observation and that polarization fraction decrease can happen concurrently with a backward scattering decrease, we would not expect to succeed in converting polarized images with no modifications to them to obtain perfect total intensity data for any disk system.

As a potentially practical application, the method can be potentially applied to single-ring systems, and/or systems that have low inclinations. Even if this approach works, however, we note potential limitations including that the pixelated (not smooth as the radiative transfer or simple geometric models generated in existing disk modeling work) disk images in polarized light, when converted to total intensity light, can yield multi-modal distributions in the retrieved posteriors of the disk parameters in Eq. (C.4). While ignoring the existence of complex structures such as multiple rings or spirals, it is still necessary to remove the pixelated noise for the approach to work.

C.5 Savitzky–Golay filter: advancing the conversion from polarization to total intensity

The pixelated residuals in Figure C.1(a) should have originated from the rotation of the non-smooth pixelated data from Qϕ observations and, thus, a smoothed version of the Qϕ data should be minimally adopted to minimize these residuals. In addition, a smoothing should not disperse the disk signals. Otherwise, smoothing would make it not directly comparable with the actual total intensity data, since the two observation modes are conducted at the same wavelengths on the same telescope instrument.

The Savitzky-Golay filter was originally used to remove noise for 1D data in Savitzky & Golay (1964). To smooth the Qϕ data for total intensity modeling, we use the Savitzky-Golay filter in two-dimension15 to minimize the random noise in Qϕ data. In comparison with other classical convolution-based methods where signals are dispersed to remove noise, the Savitzky-Golay filter instead fit p-degree polynomials to data in moving windows. Motivated by the fact that the Rayleigh resolution for Ks-band with VLT/SPHRE is 69 mas (i.e., 2.182 μm for a 8.0 m the telescope pupil seen by VLT/SPHERE), we performed a five-degree polynomial ft for a window of an 11-pixel (134.75 mas) width for the smoothing with a Savitzky–Golay filter.

Using the Savitzky-Golay filtered Qϕ data, we repeated the test described in Section C.3 and present the images in Figure C.1(b). In comparison with the results with the original Qϕ data in Figure C.1(a), the Savitzky-Golay filtered data show smoother residuals, with the total intensity modeling being qualitatively more compatible with the total intensity observation using data imputation. Although there are still strong residuals for multi-ring or spiral systems, the Savitzky-Golay filter has removed the pixelated residuals in the original Qϕ data.

We present the extracted best-fit polarization fraction curves with KLIP forward modeling in Figure C.2. Without the Savitzky-Golay filter, the best-fit polarization fraction curves in Figure C.2(a) showed different profiles, yet such profiles are not valid since the extracted profiles are even distinct for the SZ Cha data observed at different nights. The fact that the best-fit polarization fraction curves in Figure C.2(a) are different for SZ Cha suggests that the original approach cannot be used to interpret or compare the profiles. With the Savitzky-Golay filter, the polarization fraction curves in Figure C.2(b) are less distinct from each other, and the similar profiles for SZ Cha suggest that its distinct profiles in Figure C.2(a) are algorithmic effects.

For the similar profiles in different systems in Figure C.2(b), however, it is possible that the scattering angle and intensity corresponding to peak polarization fraction can still vary across different systems. Due to the assumption of independent pixels in Eq. (3), the retrieved uncertainties are extremely small in Table C.1 and thus were not presented. To further quantify the similarity or difference for the extracted profiles, proper uncertainty estimates should be performed (e.g., Wolff et al. 2017), and such estimates are beyond the scope of this study.

For the polarization fraction modeling of high-contrast imaging observations, our experiment here shows that the usage of smoothed data is necessary to reduce bias for further comparison of different systems. To fully extract the profiles for comparison, more careful treatment of the correlated uncertainties are needed. In the main text of this study in Sect. 3, we obtain the polarization fraction profiles by directly modeling the polarization fraction maps generated from PDI Qϕ and DI-sNMF Itot data. In addition, we used the results from three high-quality observations of MWC 758 in different nights, for an empirical estimation of the potential uncertainty of the polarization fraction curve modeling.

Appendix D Polarization fraction models

Here, we present the models and residuals for modeling the polarization fraction maps described in Sect. 5 here. In Fig. D.1, we show the regions and best-fit models using scaled beta distribution and Fig. D.2 contains the residuals by subtracting the models from Savitzky–Golay-smoothed observation in Fig. 5.

In the residual maps, we witness an ≈0.1 residual in polarization fraction and certain levels of patterned residuals. The patterned residuals are due to the limitations in describing the polarization fraction curves using a single scaled beta distribution. In fact, multi-ringed systems, as well as non-ring systems such as spirals, likely have different scattering properties at different locations. This effect can be also seen in the residual maps in Fig. C.1.

thumbnail Fig. D.1

Models of Ks-band polarization fraction with dimensions of 2″ × 2″ with identical color bars in linear scale (see Fig. 6 for the observation). The non-masked areas are regions used for polarization fraction modeling with scaled beta distribution in Sect. 5. The data used to create this figure are available at the CDS.
thumbnail Fig. D.2

Residuals of Ks-band polarization fraction maps in Fig. 6 (smoothed by Savitzky-Golay filter in Sect. C.5) subtracted by the models in Fig. D.1, with dimensions of 2″ × 2″ with identical color bars in linear scale. Combining the standard deviation of each residual map of ≈0.05 statistically, with an uncertainty of ≈0.03 systematically from the three MWC 758 observations in Fig. 7, we assign a total uncertainty of 0.06 for the maximum polarization fraction values.

Appendix E Contrast curves from ADI

We compared the RDI contrast curves from DI-sNMF with the RDI contrast curves from the High Contrast Data Centre and obtained consistent results. Using the High Contrast Data Centre products, we present the ADI contrast curves for the datasets presented in this work here. Specifically, we present the ADI contrast curves for both the target stars and their reference stars here; when the field rotation is not sufficient for ADI reduction, we present the non-ADI contrast curves. We compared the ANDROMEDA (Mugnier et al. 2009; Cantalloube et al. 2015), TLOCI (Marois et al. 2014), and KLIP (Soummer et al. 2012; Amara & Quanz 2012) contrast curves, and present the deepest contrasts in Fig. E.1.

thumbnail Fig. E.1

ADI contrast curves obtained from High Contrast Data Centre (left) and AMES-cond mass limits (right) for systems in this study, the line colors are consistent with Fig. 5. Disk hosts with only PDI detections (or marginal PDI detections) are displayed with light gray color with solid lines, and annotated with gray symbol; reference stars are in dotted light gray lines. For HD 163296, the ADI contrast curves are displayed and annotated in blue. The data used to create this figure are available at the CDS.

References

  1. Adams Redai, J. I., Follette, K. B., Wang, J., et al. 2023, AJ, 165, 57 [NASA ADS] [CrossRef] [Google Scholar]
  2. Allard, F., Homeier, D., Freytag, B., & Sharp, C. M. 2012, EAS Publ. Ser., 57, 3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  3. Amara, A., & Quanz, S. P. 2012, MNRAS, 427, 948 [Google Scholar]
  4. Asensio-Torres, R., Henning, T., Cantalloube, F., et al. 2021, A&A, 652, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Avenhaus, H., Quanz, S. P., Garufi, A., et al. 2018, ApJ, 863, 44 [NASA ADS] [CrossRef] [Google Scholar]
  6. Bae, J., Zhu, Z., Baruteau, C., et al. 2019, ApJ, 884, L41 [NASA ADS] [CrossRef] [Google Scholar]
  7. Bae, J., Isella, A., Zhu, Z., et al. 2023, ASP Conf. Ser., 534, 423 [NASA ADS] [Google Scholar]
  8. Balmer, W. O., Follette, K. B., Close, L. M., et al. 2022, AJ, 164, 29 [NASA ADS] [CrossRef] [Google Scholar]
  9. Baraffe, I., Chabrier, G., Barman, T. S., Allard, F., & Hauschildt, P. H. 2003, A&A, 402, 701 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Baraffe, I., Homeier, D., Allard, F., & Chabrier, G. 2015, A&A, 577, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Benisty, M., Juhasz, A., Boccaletti, A., et al. 2015, A&A, 578, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Benisty, M., Stolker, T., Pohl, A., et al. 2017, A&A, 597, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Benisty, M., Juhasz, A., Facchini, S., et al. 2018, A&A, 619, A171 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Benisty, M., Dominik, C., Follette, K., et al. 2023, ASP Conf. Ser., 534, 605 [NASA ADS] [Google Scholar]
  15. Beuzit, J. L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Boccaletti, A., Di Folco, E., Pantin, E., et al. 2020, A&A, 637, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Boccaletti, A., Pantin, E., Ménard, F., et al. 2021, A&A, 652, A8 [Google Scholar]
  19. Bohn, A. J., Benisty, M., Perraut, K., et al. 2022, A&A, 658, A183 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127 [Google Scholar]
  21. Calcino, J., Price, D. J., Pinte, C., et al. 2019, MNRAS, 490, 2579 [Google Scholar]
  22. Canovas, H., Ménard, F., de Boer, J., et al. 2015, A&A, 582, A7 [Google Scholar]
  23. Cantalloube, F., Mouillet, D., Mugnier, L. M., et al. 2015, A&A, 582, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Chen, C., Mazoyer, J., Poteet, C. A., et al. 2020, ApJ, 898, 55 [CrossRef] [Google Scholar]
  25. Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102 [Google Scholar]
  26. Cleeves, L. I., Oberg, K. I., Wilner, D. J., et al. 2016, ApJ, 832, 110 [NASA ADS] [CrossRef] [Google Scholar]
  27. Coulson, I. M., & Walther, D. M. 1995, MNRAS, 274, 977 [NASA ADS] [Google Scholar]
  28. Crotts, K. A., Matthews, B. C., Duchêne, G., et al. 2023, ApJ, submitted [arXiv:2311.14599] [Google Scholar]
  29. Cugno, G., Pearce, T. D., Launhardt, R., et al. 2023, A&A, 669, A145 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Currie, T., Marois, C., Cieza, L., et al. 2019, ApJ, 877, L3 [Google Scholar]
  31. Currie, T., Lawson, K., Schneider, G., et al. 2022, Nat. Astron., 6, 751 [NASA ADS] [CrossRef] [Google Scholar]
  32. Currie, T., Biller, B., Lagrange, A., et al. 2023, ASP Conf. Ser., 534, 799 [NASA ADS] [Google Scholar]
  33. Debes, J. H., Ren, B., & Schneider, G. 2019, J. Astron. Telescopes Instrum. Syst., 5, 035003 [NASA ADS] [Google Scholar]
  34. de Boer, J., Salter, G., Benisty, M., et al. 2016, A&A, 595, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. de Boer, J., Langlois, M., van Holstein, R. G., et al. 2020, A&A, 633, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. De Rosa, R. J., Nielsen, E. L., Wahhaj, Z., et al. 2023, A&A, 672, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Desai, N., Llop-Sayson, J., Bertrou-Cantou, A., et al. 2022, SPIE Conf. Ser., 12180, 121805H [NASA ADS] [Google Scholar]
  38. Dohlen, K., Langlois, M., Saisse, M., et al. 2008, SPIE Conf. Ser., 7014, 70143L [Google Scholar]
  39. Dong, R., Hashimoto, J., Rafikov, R., et al. 2012, ApJ, 760, 111 [NASA ADS] [CrossRef] [Google Scholar]
  40. Dong, R., Zhu, Z., Rafikov, R. R., & Stone, J. M. 2015, ApJ, 809, L5 [Google Scholar]
  41. Dullemond, C. P., & Dominik, C. 2005, A&A, 434, 971 [CrossRef] [EDP Sciences] [Google Scholar]
  42. Engler, N., Schmid, H. M., Thalmann, C., et al. 2017, A&A, 607, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. Engler, N., Schmid, H. M., Quanz, S. P., Avenhaus, H., & Bazzon, A. 2018, A&A, 618, A151 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Engler, N., Milli, J., Gratton, R., et al. 2023, A&A, 672, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  45. Flasseur, O., Thé, S., Denis, L., Thiébaut, É., & Langlois, M. 2021, A&A, 651, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
  47. Francis, L., & van der Marel, N. 2020, ApJ, 892, 111 [Google Scholar]
  48. Franson, K., Bowler, B. P., Zhou, Y., et al. 2023, ApJ, 950, L19 [NASA ADS] [CrossRef] [Google Scholar]
  49. Frasca, A., Biazzo, K., Lanzafame, A. C., et al. 2015, A&A, 575, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  50. Frattin, E., Martikainen, J., Munoz, O., et al. 2022, MNRAS, 517, 5463 [NASA ADS] [CrossRef] [Google Scholar]
  51. Fulton, B. J., Rosenthal, L. J., Hirsch, L. A., et al. 2021, ApJS, 255, 14 [NASA ADS] [CrossRef] [Google Scholar]
  52. Gaia Collaboration (Vallenari, A., et al.) 2023, A&A, 674, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. Garufi, A., Quanz, S. P., Schmid, H. M., et al. 2016, A&A, 588, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. Garufi, A., Benisty, M., Pinilla, P., et al. 2018, A&A, 620, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Garufi, A., Avenhaus, H., Pérez, S., et al. 2020, A&A, 633, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Garufi, A., Dominik, C., Ginski, C., et al. 2022, A&A, 658, A137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Ginski, C., Stolker, T., Pinilla, P., et al. 2016, A&A, 595, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  58. Ginski, C., Benisty, M., van Holstein, R. G., et al. 2018, A&A, 616, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  59. Ginski, C., Ménard, F., Rab, C., et al. 2020, A&A, 642, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  60. Ginski, C., Tazaki, R., Dominik, C., & Stolker, T. 2023, ApJ, 953, 92 [NASA ADS] [CrossRef] [Google Scholar]
  61. Gray, R. O., Riggs, Q. S., Koen, C., et al. 2017, AJ, 154, 31 [NASA ADS] [CrossRef] [Google Scholar]
  62. Groff, T. D., Riggs, A. J. E., Kern, B., & Kasdin, N. J. 2016, J. Astron. Telescopes Instrum. Syst., 2, 011009 [Google Scholar]
  63. Guyon, O., Norris, B., Martinod, M.-A., et al. 2021, SPIE Conf. Ser., 11823, 1182318 [NASA ADS] [Google Scholar]
  64. Haffert, S. Y., Bohn, A. J., de Boer, J., et al. 2019, Nat. Astron., 3, 749 [Google Scholar]
  65. Hammond, I., Christiaens, V., Price, D. J., et al. 2023, MNRAS, 522, L51 [NASA ADS] [CrossRef] [Google Scholar]
  66. Hauschildt, P. H., Allard, F., Ferguson, J., Baron, E., & Alexander, D. R. 1999, ApJ, 525, 871 [Google Scholar]
  67. Herbig, G. H. 1977, ApJ, 214, 747 [Google Scholar]
  68. Herbig, G. H., Vrba, F. J., & Rydgren, A. E. 1986, AJ, 91, 575 [NASA ADS] [CrossRef] [Google Scholar]
  69. Herczeg, G. J., & Hillenbrand, L. A. 2014, ApJ, 786, 97 [Google Scholar]
  70. Houk, N., & Swift, C. 1999, Michigan Spectral Survey, 5, 0 [Google Scholar]
  71. Irvine, N. J., & Houk, N. 1977, PASP, 89, 347 [NASA ADS] [CrossRef] [Google Scholar]
  72. Jones, M. I., Milli, J., Blanchard, I., et al. 2022, A&A, 667, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  73. Joy, A. H. 1949, ApJ, 110, 424 [NASA ADS] [CrossRef] [Google Scholar]
  74. Juillard, S., Christiaens, V., & Absil, O. 2022, A&A, 668, A125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  75. Juillard, S., Christiaens, V., & Absil, O. 2023, A&A, 679, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  76. Keppler, M., Benisty, M., Müller, A., et al. 2018, A&A, 617, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  77. Kiselev, N. N., Rosenbush, V. K., Petrov, D., et al. 2022, MNRAS, 514, 4861 [CrossRef] [Google Scholar]
  78. Kraus, A. L., & Ireland, M. J. 2012, ApJ, 745, 5 [Google Scholar]
  79. Law, C. J., Booth, A. S., & Oberg, K. I. 2023, ApJ, 952, L19 [NASA ADS] [CrossRef] [Google Scholar]
  80. Lawson, K., Currie, T., Wisniewski, J. P., et al. 2022, ApJ, 935, L25 [NASA ADS] [CrossRef] [Google Scholar]
  81. Long, F., Herczeg, G. J., Pascucci, I., et al. 2018, ApJ, 863, 61 [NASA ADS] [CrossRef] [Google Scholar]
  82. Long, F., Andrews, S. M., Zhang, S., et al. 2022, ApJ, 937, L1 [NASA ADS] [CrossRef] [Google Scholar]
  83. Ma, J., Schmid, H. M., & Tschudi, C. 2023, A&A, 676, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  84. Maire, A.-L., Langlois, M., Dohlen, K., et al. 2016, SPIE Conf. Ser., 9908, 990834 [Google Scholar]
  85. Maire, A.-L., Langlois, M., Delorme, P., et al. 2021, J. Astron. Telescopes Instrum. Syst., 7, 035004 [NASA ADS] [Google Scholar]
  86. Marois, C., Lafrenière, D., Doyon, R., Macintosh, B., & Nadeau, D. 2006, ApJ, 641, 556 [Google Scholar]
  87. Marois, C., Macintosh, B., Barman, T., et al. 2008, Science, 322, 1348 [Google Scholar]
  88. Marois, C., Correia, C., Galicher, R., et al. 2014, SPIE Conf. Ser., 9148, 91480U [Google Scholar]
  89. Mawet, D., Milli, J., Wahhaj, Z., et al. 2014, ApJ, 792, 97 [Google Scholar]
  90. Mesa, D., Gratton, R., Kervella, P., et al. 2023, A&A, 672, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  91. Milli, J., Mouillet, D., Lagrange, A. M., et al. 2012, A&A, 545, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  92. Milli, J., Vigan, A., Mouillet, D., et al. 2017, A&A, 599, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  93. Monnier, J. D., Harries, T. J., Bae, J., et al. 2019, ApJ, 872, 122 [Google Scholar]
  94. Mora, A., Merin, B., Solano, E., et al. 2001, A&A, 378, 116 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  95. Mugnier, L. M., Cornia, A., Sauvage, J.-F., et al. 2009, J. Opt. Soc. Am. A., 26, 1326 [NASA ADS] [CrossRef] [Google Scholar]
  96. Mulders, G. D., Pascucci, I., Manara, C. F., et al. 2017, ApJ, 847, 31 [NASA ADS] [CrossRef] [Google Scholar]
  97. Muñoz, O., Frattin, E., Jardiel, T., et al. 2021, ApJS, 256, 17 [CrossRef] [Google Scholar]
  98. Nielsen, E. L., De Rosa, R. J., Macintosh, B., et al. 2019, AJ, 158, 13 [Google Scholar]
  99. Nousiainen, T., Muinonen, K., & RäIsäNen, P. 2003, J. Geophys. Res. (Atmos.), 108, 4025 [NASA ADS] [CrossRef] [Google Scholar]
  100. Ochsenbein, F., Bauer, P., & Marcout, J. 2000, A&AS, 143, 23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  101. Ohta, Y., Fukagawa, M., Sitko, M. L., et al. 2016, PASJ, 68, 53 [Google Scholar]
  102. Olofsson, J., van Holstein, R. G., Boccaletti, A., et al. 2018, A&A, 617, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  103. Olofsson, J., Thébault, P., Bayo, A., et al. 2023, A&A, 674, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  104. Pairet, B., Cantalloube, F., & Jacques, L. 2021, MNRAS, 503, 3724 [Google Scholar]
  105. Pecaut, M. J., & Mamajek, E. E. 2016, MNRAS, 461, 794 [Google Scholar]
  106. Pence, W. D., Chiappetti, L., Page, C. G., Shaw, R. A., & Stobie, E. 2010, A&A, 524, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  107. Perrin, M. D., Duchene, G., Millar-Blanchaer, M., et al. 2015, ApJ, 799, 182 [Google Scholar]
  108. Pinte, C., Price, D. J., Ménard, F., et al. 2018, ApJ, 860, L13 [Google Scholar]
  109. Pinte, C., Price, D. J., Ménard, F., et al. 2020, ApJ, 890, L9 [CrossRef] [Google Scholar]
  110. Price, D. J., Cuello, N., Pinte, C., et al. 2018, MNRAS, 477, 1270 [Google Scholar]
  111. Quanz, S. P., Amara, A., Meyer, M. R., et al. 2015, ApJ, 807, 64 [Google Scholar]
  112. Quiroz, J., Wallack, N. L., Ren, B., et al. 2022, ApJ, 924, L4 [NASA ADS] [CrossRef] [Google Scholar]
  113. Rameau, J., Follette, K. B., Pueyo, L., et al. 2017, AJ, 153, 244 [Google Scholar]
  114. R Core Team 2022, R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria [Google Scholar]
  115. Reggiani, M., Christiaens, V., Absil, O., et al. 2018, A&A, 611, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  116. Ren, B., Pueyo, L., Zhu, G. B., Debes, J., & Duchêne, G. 2018, ApJ, 852, 104 [Google Scholar]
  117. Ren, B., Choquet, É., Perrin, M. D., et al. 2019, ApJ, 882, 64 [Google Scholar]
  118. Ren, B., Pueyo, L., Chen, C., et al. 2020, ApJ, 892, 74 [Google Scholar]
  119. Ren, B., Choquet, É., Perrin, M. D., et al. 2021, ApJ, 914, 95 [CrossRef] [Google Scholar]
  120. Ren, B. B., Rebollido, I., Choquet, É., et al. 2023, A&A, 672, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  121. Riviere-Marichalar, P., Ménard, F., Thi, W. F., et al. 2012, A&A, 538, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  122. Ruane, G., Ngo, H., Mawet, D., et al. 2019, AJ, 157, 118 [NASA ADS] [CrossRef] [Google Scholar]
  123. Rydgren, A. E. 1980, AJ, 85, 444 [Google Scholar]
  124. Sai Krishanth, P. M., Douglas, E. S., Hom, J., et al. 2023, Proc. SPIE, 12680, 1268021 [NASA ADS] [Google Scholar]
  125. Sallum, S., Follette, K. B., Eisner, J. A., et al. 2015, Nature, 527, 342 [Google Scholar]
  126. Savitzky, A., & Golay, M. J. E. 1964, Analyt. Chem., 36, 1627 [NASA ADS] [CrossRef] [Google Scholar]
  127. Schmid, H. M., Bazzon, A., Roelfsema, R., et al. 2018, A&A, 619, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  128. Shuai, L., Ren, B. B., Dong, R., et al. 2022, ApJS, 263, 31 [NASA ADS] [CrossRef] [Google Scholar]
  129. Siess, L., Dufour, E., & Forestini, M. 2000, A&A, 358, 593 [Google Scholar]
  130. Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163 [Google Scholar]
  131. Soummer, R., Pueyo, L., & Larkin, J. 2012, ApJ, 755, L28 [Google Scholar]
  132. Spiegel, D. S., & Burrows, A. 2012, ApJ, 745, 174 [Google Scholar]
  133. Stadler, J., Benisty, M., Izquierdo, A., et al. 2023, A&A, 670, A1 [Google Scholar]
  134. Stolker, T., Dominik, C., Min, M., et al. 2016, A&A, 596, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  135. Stolker, T., Bonse, M. J., Quanz, S. P., et al. 2019, A&A, 621, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  136. Sutlieff, B. J., Birkby, J. L., Stone, J. M., et al. 2023, MNRAS, 520, 4235 [NASA ADS] [CrossRef] [Google Scholar]
  137. Szulágyi, J., & Garufi, A. 2021, MNRAS, 506, 73 [CrossRef] [Google Scholar]
  138. Tanaka, H., Himeno, Y., & Ida, S. 2005, ApJ, 625, 414 [NASA ADS] [CrossRef] [Google Scholar]
  139. Tazaki, R., & Dominik, C. 2022, A&A, 663, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  140. Tazaki, R., Tanaka, H., Muto, T., Kataoka, A., & Okuzumi, S. 2019, MNRAS, 485, 4951 [NASA ADS] [CrossRef] [Google Scholar]
  141. Tazaki, R., Ginski, C., & Dominik, C. 2023, ApJ, 944, L43 [NASA ADS] [CrossRef] [Google Scholar]
  142. Teague, R., Bae, J., Bergin, E. A., Birnstiel, T., & Foreman-Mackey, D. 2018, ApJ, 860, L12 [NASA ADS] [CrossRef] [Google Scholar]
  143. Thalmann, C., Janson, M., Garufi, A., et al. 2016, ApJ, 828, L17 [Google Scholar]
  144. Torres, C. A. O., Quast, G. R., da Silva, L., et al. 2006, A&A, 460, 695 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  145. Tschudi, C., & Schmid, H. M. 2021, A&A, 655, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  146. Uyama, T., Ren, B., Mawet, D., et al. 2020, AJ, 160, 283 [NASA ADS] [CrossRef] [Google Scholar]
  147. van Holstein, R. G., Snik, F., Girard, J. H., et al. 2017, SPIE Conf. Ser., 10400, 1040015 [Google Scholar]
  148. van Holstein, R. G., Girard, J. H., de Boer, J., et al. 2020, A&A, 633, A64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  149. van Holstein, R. G., Stolker, T., Jensen-Clem, R., et al. 2021, A&A, 647, A21 [EDP Sciences] [Google Scholar]
  150. Vieira, S. L. A., Corradi, W. J. B., Alencar, S. H. P., et al. 2003, AJ, 126, 2971 [Google Scholar]
  151. Vigan, A., Fontanive, C., Meyer, M., et al. 2021, A&A, 651, A72 [EDP Sciences] [Google Scholar]
  152. Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
  153. Wagner, K., Apai, D., Kasper, M., & Robberto, M. 2015, ApJ, 813, L2 [Google Scholar]
  154. Wagner, K., Dong, R., Sheehan, P., et al. 2018, ApJ, 854, 130 [Google Scholar]
  155. Wagner, K., Stone, J. M., Spalding, E., et al. 2019, ApJ, 882, 20 [Google Scholar]
  156. Wagner, K., Stone, J., Dong, R., et al. 2020, AJ, 159, 252 [NASA ADS] [CrossRef] [Google Scholar]
  157. Wagner, K., Stone, J., Skemer, A., et al. 2023, Nat. Astron., 7, 1208 [CrossRef] [Google Scholar]
  158. Wahhaj, Z., Milli, J., Romero, C., et al. 2021, A&A, 648, A26 [EDP Sciences] [Google Scholar]
  159. Wallack, N. L., Ruffio, J.-B., Ruane, G., et al. 2023, AJ, submitted [Google Scholar]
  160. Wang, J. J., Ginzburg, S., Ren, B., et al. 2020, AJ, 159, 263 [Google Scholar]
  161. Wenger, M., Ochsenbein, F., Egret, D., et al. 2000, A&AS, 143, 9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  162. Wilking, B. A., Meyer, M. R., Robinson, J. G., & Greene, T. P. 2005, AJ, 130, 1733 [NASA ADS] [CrossRef] [Google Scholar]
  163. Wölfer, L., Facchini, S., van der Marel, N., et al. 2023, A&A, 670, A154 [CrossRef] [EDP Sciences] [Google Scholar]
  164. Wolff, S. G., Perrin, M. D., Stapelfeldt, K., et al. 2017, ApJ, 851, 56 [NASA ADS] [CrossRef] [Google Scholar]
  165. Xie, C., Choquet, E., Vigan, A., et al. 2022, A&A, 666, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  166. Xie, C., Ren, B. B., Dong, R., et al. 2023, A&A, 675, A1 [Google Scholar]
  167. Zhou, Y., Bowler, B. P., Yang, H., et al. 2023, AJ, 166, 220 [NASA ADS] [CrossRef] [Google Scholar]
2

The derived uncertainties in this paper are 1σ (frequentist) or 68% credible intervals (Bayesian) unless otherwise specified.

8

We adopt Gaia Rp-band here, since Wahhaj et al. (2021) referenced R-band for star-hopping observations.

All Tables

Table 1

System parameters and star-hopping observation log.

In the text
Table A.1

Archival SPHERE observations of Ks-band counterparts in the Y, J, or H bands in polarized light

In the text
Table C.1

Best-fit scaled beta distribution description for polarization fraction curve for Fig. 6

In the text

All Figures

thumbnail Fig. 1

Ks-band Qϕ maps with dimensions of 2″ × 2″ with different color bars in log scale. The letter identifiers are from Table 1. The rulers correspond to 50 au. The regions interior to 0”. 1 are not accessible with coronagraph usage. We note that the data used to create this figure are available at the CDS.
In the text
thumbnail Fig. 2

Ks-band Itot maps with dimensions of 2″ × 2″ with different color bars in log scale, for disks identifiable with DI-sNMF star-hopping RDI data reduction. Due to varying observation conditions, a well-chosen reference star may not lead to optimum reduction results, e.g., SZ Cha in panels x and x′. We note that the data used to create this figure are available at the CDS.
In the text
thumbnail Fig. 3

Disks with confident recovery in Ks-band via NMF data imputation, as a function of (a) Gaia DR3 Rp magnitude or (b) 2MASS Ks magnitude with the Rp – K color. Each connected pair is a host-reference pair, with their on-sky angular separation from Gaia DR3 in degrees. Notes: (1) Certain systems with marginal detections in NMF data imputation are marked as non-detection (e.g., HD 163296 or m). (2) The size of reference star symbols reflects typical uncertainties in color-magnitude measurements, zoom in the figure for actual error bars colored yellow for all systems.
In the text
thumbnail Fig. 4

Individual logistic regression on RDI detectability of disks in DI-sNMF. The solid dots are for systems in Fig. 2. (a) Brighter host stars can likely yield disk detection. (b) Redder Rp – K band color potentially contribute to disk detection. (c) Brighter reference stars in Ks-band might tentatively contribute to disk detection. (d) Reference stars with fainter Rp magnitude of target stars might tentatively contribute to disk detection, however that is likely biased by the non-detections with Rp ≳ 11, see Fig. 3a. (e) Although reference at farther separation might yield better disk detection, this stems from a selection bias when close-in references are not available. Notes: (1) A joint multinomial logistic regression only preserves the trend for (a). (2) The trends here are subject to change due to the highly selected samples in this study. (3) From solid to dashed lines, the statistical significance decreases; dotted lines are biased trends that are not trustworthy.
In the text
thumbnail Fig. 5

RDI contrast curves (left) and AMES-Cond mass limits (right) of DI-sNMF reductions for systems in Fig. 2 (see Sect. 4 for the detailed calculation procedure). Within ~0″.2, the measured contrast is prone to transmission reduction near the coronagraphic edge. The contrast peak at −1″ for j (HD 100453) is the binary companion. The contrast curve for t (PDS 66) is the deepest due to the faintness of the disk. The data used to create this figure are available at the CDS.
In the text
thumbnail Fig. 6

Ks-band polarization fraction maps with dimensions of 2″ × 2″ with identical color bars in linear scale, obtained from dividing the Qϕ maps in Fig. 1 by the corresponding Itot maps in Fig. 2. The data used to create this figure are available at the CDS.
In the text
thumbnail Fig. 7

Ks-band polarization fraction curves, see Appendix C.1 and Table C.1 for the analytical expressions using scaled beta distribution. Best-fit curves with gray segments inaccessible from observation, based on system inclination assuming infinitely thin disks. The error bar on the top right, which could be used to infer the systematic uncertainty from the fitting method, is from the standard deviation of three sets of MWC 758 best-fit results.
In the text
thumbnail Fig. 8

Scattering peak values in angle and fraction (SPAF) plot. (a) Tentative correlation between peak polarization fraction value and peak scattering angle, with the gray lines randomly fitting 6 to 12 systems for correlation exploration. The best-fit expression is performed on all data points here with a 0.06 uncertainty for the peak polarization fraction, and certain systems (For b or CQ Tau and k or HD 100546, marked with hollow symbols; for l or HD 143006, a scaled beta distribution cannot describe its polarization fraction curve, potentially due to large scale shadowing in Benisty et al. 2018.) were excluded from Fig. 7 to illustrate this relationship (see Sect. 5.2.1 for a detailed discussion). (b) Gaussian random sphere (GRS) dust models from Tazaki & Dominik (2022) and Tazaki et al. (2023) by varying the minimum dust size for irregular compact grains, overlaid on the 1σ and 2σ ranges from the resampling results in (a). We observe a similar trend for absorptive material (“amc”) but not for less absorptive materials (“org”), with the latter neither providing peak scattering angles with the observations that are consistent beyond ~100° (see Sect. 5.2.2). We note that the dust models do not necessarily reproduce the polarization fraction curves in Fig. 7.
In the text
thumbnail Fig. 9

Polarized color at ≈90° scattering angle and stellar luminosity in Sect. 5.3. (a), (b), and (c) are the Y, J, and H-band data in polarized light in comparison with Ks-band data in polarized light, respectively. In J – Ks comparison, increase in stellar luminosity leads to more neutral color; while the correlation is less evidence in H – Ks comparison likely due to adjacent wavelengths (e.g., Fig. A.2). Note: the bands are 1σ, 2σ, and 3σ confidence intervals from bootstrapping fit. The data used to create this figure are available at the CDS.
In the text
thumbnail Fig. A.1

Ks-band 𝒰ϕ maps with dimensions of 2″ × 2″ with different color bars in linear scale. The field of view of each panel corresponds to those in Fig. 1. The data used to create this figure are available at the CDS.
In the text
thumbnail Fig. A.2

Transmission profiles, as well as the central wavelengths, for SPHERE/IRDIS in the Y, J, H, and Ks bands.
In the text
thumbnail Fig. C.1

Qϕ to total intensity conversion using polarization fraction with scaled beta distribution, (a): direct polarization map division, (b): Savitzky–Golay filtered Qϕ data. From top to bottom, A: LkCa 15, on 2020 December 8; B1: SZ Cha, on 2020 December 29; B2: SZ Cha, on 2020 December 30; C: V1247 Ori, on 2020 December 24. From left to right, 1: Qϕ data; 2: total intensity data from data imputation; 3: total intensity data converted from Qϕ data using polarization fraction maps; 4: RDI KLIP residuals using total intensity models; 5: polarization maps used to model total intensity data with RDI KLIP. A comparison between (a) and (b) shows an improvement of retrieving quality with the Savitzky–Golay filter. Nevertheless, the patterned residuals show that a single profile is limited in describing polarization fraction maps, especially when multiple disk components exist: modeling for separate disk components (e.g., rings in LkCa 15 and SZ Cha) is needed for authentic description of polarization fraction curves.
In the text
thumbnail Fig. C.2

Best-fit polarization curves assuming beta distribution from RDI KLIP forward modeling in Figure C.1, the peak polarization locations are marked with circles. The light gray curves denote scattering angle ranges that are not accessible for inclined thin disks with no flaring. (a) Direct conversion: the extracted profiles are different even for the same system (i.e., B1 and B2). (b) Conversion with Savitzky-Golay filtered Qϕ data: the extracted curves are less distinct, since such reductions are less prone to pixel-wise discrete noise.
In the text
thumbnail Fig. D.1

Models of Ks-band polarization fraction with dimensions of 2″ × 2″ with identical color bars in linear scale (see Fig. 6 for the observation). The non-masked areas are regions used for polarization fraction modeling with scaled beta distribution in Sect. 5. The data used to create this figure are available at the CDS.
In the text
thumbnail Fig. D.2

Residuals of Ks-band polarization fraction maps in Fig. 6 (smoothed by Savitzky-Golay filter in Sect. C.5) subtracted by the models in Fig. D.1, with dimensions of 2″ × 2″ with identical color bars in linear scale. Combining the standard deviation of each residual map of ≈0.05 statistically, with an uncertainty of ≈0.03 systematically from the three MWC 758 observations in Fig. 7, we assign a total uncertainty of 0.06 for the maximum polarization fraction values.
In the text
thumbnail Fig. E.1

ADI contrast curves obtained from High Contrast Data Centre (left) and AMES-cond mass limits (right) for systems in this study, the line colors are consistent with Fig. 5. Disk hosts with only PDI detections (or marginal PDI detections) are displayed with light gray color with solid lines, and annotated with gray symbol; reference stars are in dotted light gray lines. For HD 163296, the ADI contrast curves are displayed and annotated in blue. The data used to create this figure are available at the CDS.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.