Protoplanetary disks in K s -band total intensity and polarized light ⋆

Context. Diverse morphology in protoplanetary disks can result from planet-disk interaction, suggesting the presence of forming planets. Characterizing disks can inform the formation environments of planets. To date, most imaging campaigns have probed the polarized light from disks, which is only a fraction of the total scattered light and not very sensitive to planetary emission. Aims. We aim to observe and characterize protoplanetary disk systems in the near-infrared in both polarized and total intensity light, to carry out an unprecedented study of the dust scattering properties of disks, as well as of any possible planetary companions. Methods. Using the star-hopping mode of the SPHERE instrument at the Very Large Telescope, we observed 29 young stars hosting protoplanetary disks and their reference stars in the K s -band polarized light. We extracted disk signals in total intensity by removing stellar light using the corresponding reference star observations, by adopting the data imputation concept with sequential non-negative matrix factorization (DI-sNMF). For well-recovered disks in both polarized and total intensity light, we parameterized the polarization fraction phase functions using scaled beta distribution. We investigated the empirical DI-sNMF detectability of disks using logistic regression. For systems with SPHERE data in Y - / J - / H -band, we summarized their polarized color at ≈ 90 ◦ scattering angle. Results. We obtained high-quality disk images in total intensity for 15 systems and in polarized light for 23 systems. Total intensity detectability of disks primarily depends on host star brightness, which determines adaptive-optics control ring imagery and thus stellar signals capture using DI-sNMF. The peak of polarization fraction tentatively correlates with the peak scattering angle, which could be reproduced using certain composition for compact dust, yet more detailed modeling studies are needed. Most of disks are blue in polarized J − K s color, and the fact that they are relatively redder as stellar luminosity increases indicates larger scatterers. Conclusions. High-quality disk imagery in both total intensity and polarized light allows for disk characterization in polarization fraction. The combination of them reduces the confusion between disk and planetary signals.


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. 2015Benisty et al. , 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 planetdisk 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. 2018Pinte et al. , 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. 2019Wagner et al. , 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 highcontrast 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. 2022Juillard et al. , 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 protoplanetary disks in total intensity from the ground.As many as 29 young stars are surveyed in the K s -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.
A114, page 2 of 26 Ren, B. B., et al.: A&A, 680, A114 (2023)   A114, page 3 of 26 Ren, B. B., et al.: A&A, 680, A114 (2023) 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 T eff from VizieR (Ochsenbein et al. 2000) and the photometry using the VizieR photometry tool1 .We calculated the stellar luminosity L star through a PHOENIX model (Hauschildt et al. 1999) scaled to the de-reddened brightness in V-band.From T eff and L star , 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 T eff and L star , and therefore of M star 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).
By observing these systems in K s -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 referencewhich 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.

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).

Preprocessing
In the star-hopping observation of a host-reference pair, the host has dedicated sky background observations for empirical K s -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.

Polarimetric differential imaging (PDI): Q ϕ 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-polarizationsubtracted Q ϕ 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 K s -band data from this study.

Reference diffraction imaging (RDI): I tot images
For total intensity data that can be generated using the polarized observations, we performed RDI data reduction to obtain total intensity (I tot ) 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 I tot 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 K s -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.
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 I tot 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-K s bands, we additionally convolved the data to reach same spatial resolution as K s -band, based on Rayleigh diffraction limits of 1.22λ/D where λ is the central wavelength 3 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.

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, I tot detections from RDI data reduction using DI-sNMF, we present the images in Fig. 2.

PDI Q ϕ 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 K s -band in Fig. 1 do not have significant morphological variations from them.The K s -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 U ϕ 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ϕ in Ma et al. 2023) Ren, B. B., et al.: A&A, 680, A114 (2023) simulation).To further reduce the limitation on U ϕ 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.

RDI I tot 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 K s -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 K s -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 K s -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.

RDI I tot 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 R p − K color and R p -band (or K-band) magnitude.The R p -band and K-band magnitudes are from Gaia DR3 (Gaia Collaboration 2023) and 2MASS (Skrutskie et al. 2006), respectively.
Beyond an empirical threshold of R p ≳ 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.

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 R p -and K s -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 K s -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 R p -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 R p ≳ 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 R p -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 K s -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).

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 R p -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 starhopping 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 K s -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 R p -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 increase in the S/N for the pixels containing control ring signals.

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 K s -band polarization fraction of 1.0 +1.5 −1.0 %, which is measured here using the PDI total polarized intensity and RDI total intensity results.The polarized HD 100453B suggests that there exists circumsecondary 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 starhopping RDI datasets.We also compare our K s -band results with previous claims and discuss potential reasons of their non-detection in K s -band here.

Detection limits
To quantify the detection capability of point sources, existing surveys (Nielsen et al. 2019;Vigan et al. 2021)   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.

Contrast calculation
Using the preprocessed and postprocessed data, we calculated the RDI contrast curves from DI-sNMF for the disks with highquality 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 Centre 9 and we did not observe any significant differences.

Contrast results
For RDI with DI-sNMF, we reached for K s -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 K 1 -or K 2 -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 3 , the filters used in Wahhaj et al. (2021) have narrower wavelength coverage (K 1 -and K 2 -band) than the K s -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 R p -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 K s -band data here with other observations, we used AMES-Cond models (Allard et al. 2012) to 9 https://sphere.osug.fr/spip.php?rubrique16&lang=en 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 fieldtracking, 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.

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 2020-12-08  LkHa 330 2020-12-08 MWC 758 2020-12-19   2020-12-22   MWC 758 2020-12-23   V1247 Ori 2020-12-24 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 K s -band.Second, even if existing claimed planetary objects are visible in K s -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 K s -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.

Polarization fraction and color
With the PDI Q ϕ and RDI I tot 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.

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 A114, page 11 of 26 Ren, B. B., et al.: A&A, 680, A114 (2023) between the incident light vector and the scattered light vector, the polarization fraction follows 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, f max pol , 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 diskmap in cylindrical coordinates follows where h 0 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 K s -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 f max pol , together with geometrical parameters (i.e., scale height h 0 , 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 loglikelihood function, ln  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.
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 X obs and X model 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 < h 0 < 0.2, 0 < γ < 2, 0 ≤ f max pol ≤ 1, 1 < α ≤ 5 and 1 < β ≤ 5. Using the polarization fraction data from Fig. 6, we present the bestfit 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 Table C.1.
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  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. 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.

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 selfshadowing 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 A114, page 13 of 26 Ren, B. B., et al.: A&A, 680, A114 (2023) 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.

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 K s -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 K s -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 ν = −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 K s -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 ∝ a 1.5 c with a c being the characteristic radius of an aggregate) with monomer size a mon = 0.2 µm in the IM Lup disk surface when observed in H-band.In comparison, their best-fitting aggregate model would suggest f max pol = 0.83 and θ max = 89 • at the K s -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 longlived 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.

Disk color in polarized light
From the color images, calculated by comparing Y-, J-, or Hband data with K s 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 K s bands for stars that are less luminous than ∼10 L ⊙ .In H and K s bands, the colors of the disks can vary between blue and red.When stellar luminosity increases, we observe that both the J − K s and the H − K s color in polarized light are redder, and that the H − K s polarized color changes slower than that of J − K s .The slow H − K s 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 K s -band in polarized light.

Summary
We obtained K s -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 K s -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 K s -band detection limits in comparison to the exploration in K 1 -/K 2 -bands in Wahhaj et al. (2021).Nevertheless, an actual detection is a trade-off between contrast and bandintegrated companion luminosity, and thus narrower bands do not necessarily always provide better detections.Given that starhopping 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) dust models 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.
A114, page 15 of 26 Ren, B. B., et al.: A&A, 680, A114 (2023) 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 J pol − K s pol and H pol − K s 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 − K s 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 R p ≳ 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.

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 ∈ R 1×N pix

≥0
with N pix ∈ N pixels, with an NMF component basis vector

≥0
, we can denote the corresponding coefficient with ω i ∈ R ≥0 .When a fraction of the data in T is missing (or artificially ignored here), the corresponding coefficient is ω ′ i ∈ R ≥0 .With these notations, Theorem 2 in Ren et al. (2020) states that 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: 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-order deviation for the cross-talk terms of H in the second multiplier without loss of generality.Similarly to Eq. (B.1), we have where 1 T ∈ B 1×N pix is an indicator matrix (in which 1 T j = 0 when the corresponding element in T j is missing, and 1 T j = 1 otherwise) which matches the dimension of T and H i .
Comparing the target modeling terms in Eq. 33 of Ren et al. (2020), 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 for future 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 (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: In terms of statistics, the probability density function (PDF) of a beta distribution follows: for x ∈ [0, 1], and Γ(•) is the gamma function with Γ(x) = ∞ 0 t x e −t dt for ∀x ∈ R + .We can normalize Eq. ( 1) to have θ scat π follow a beta distribution form, To enable observational description of polarization fraction, we can set the maximum polarization fraction to be f max pol ∈ [0, 1].We now have a polarization fraction curve of 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 (i.e., where the polarization fraction peaks), and thus it is used to regulate the maximum polarization fraction to be f max pol 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) , 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: where r ∈ R + is the stellocentric distance in the midplane of the disk, h 0 ∈ R + is the disk scale height at 1 au, and γ ∈ R + 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): 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.

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  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 A114, page 23 of 26 Ren, B. B., et al.: A&A, 680, A114 (2023) 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.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-dimension 11 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 K s -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 fivedegree 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. 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 A114, page 24 of 26 Ren, B. B., et al.: A&A, 680, A114 (2023) 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 I tot 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.
due to convolution effects.Observationally, studies including Engler et al. (2023) demonstrated that the U ϕ 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_PRI M band (∼0.790 µm), with the latter having higher angular resolution.For the K s -band (∼2.182 µm) observations in this study, we thus expect the U ϕ maps contain more leakage from Q ϕ data due to a finite angular resolution.We observe that the U ϕ 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 A114, page 6 of 26

Fig. 3 .
Fig. 3. Disks with confident recovery in K s -band via NMF data imputation, as a function of (a) Gaia DR3 R p magnitude or (b) 2MASS K s magnitude with the R p − 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.

Fig. 4 .
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 R p − K band color potentially contribute to disk detection.(c) Brighter reference stars in K s -band might tentatively contribute to disk detection.(d) Reference stars with fainter R p magnitude of target stars might tentatively contribute to disk detection, however that is likely biased by the non-detections with R p ≳ 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.

Fig. 5 .
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.

Fig. 6 .
Fig. 6.K s -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 I tot maps in Fig. 2. The data used to create this figure are available at the CDS.andMWC 758(Reggiani et al. 2018;Wagner et al. 2019Wagner et al. , 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 I tot 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 K s -band.Second, even if existing claimed planetary objects are visible in K s -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 inFranson 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 K s -band than in J-/H-band and thus worse spatial resolution (i.e., 1.22λ/D, see Fig. A.2 for the IRDIS filters),

Fig. 7 .
Fig. 7. K s -band polarization fraction curves, see Appendix C.1 and TableC.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.

Fig. 8 .
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 inBenisty et al. 2018.)  were excluded from Fig.7to illustrate this relationship (see Sect. 5.2.1 for a detailed discussion).(b) Gaussian random sphere (GRS) dust models fromTazaki & Dominik (2022) andTazaki 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.

Fig. 9 .
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 K s -band data in polarized light, respectively.In J − K s comparison, increase in stellar luminosity leads to more neutral color; while the correlation is less evidence in H − K s 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.

Fig. A. 1 .
Fig. A.1.K s -band U ϕ 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.
Figure C.1(a) for the results for both morphologies and Figure C.2(a) for the best-fit polarization curves.A114, page 22 of 26

Fig
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; B 1 : SZ Cha, on 2020 December 29; B 2 : 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.
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 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., B 1 and B 2 ).(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.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 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.

Table 1 .
System parameters and star-hopping observation log.

Table A .
1. Archival SPHERE observations of K s -band counterparts in the Y, J, or H bands in polarized light Letter identifiers of the targets in this paper, the • • • symbols are used for systems with no existing polarized observations in other bands or without confident detection in K s -band for polarized color extraction.Column (2): Target name.Columns (3), (6), and (9): UTC observation dates.Columns (4), (7), and (10): Total on-source exposure time for the target.Columns (5), (8), and (11): ESO Program ID.Different observation rows on the same observation night indicate different observation setups.For some observations, there were no PSF frames for relative flux measurement.
A114, page 19 of 26 Fig. A.2. Transmission profiles, as well as the central wavelengths, for SPHERE/IRDIS in the Y, J, H, and K s bands.