© The Authors 2025

1 Introduction

The field of exoplanet research, which has rapidly evolved in recent decades, has uncovered an immense diversity in planetary structure and composition, ranging from small rocky worlds to massive gas giants, orbiting their stars in periods spanning days to hundreds of years (Jontof-Hutter 2019; Winn & Fabrycky 2015; Dawson & Johnson 2018; Zhu & Dong 2021; Wordsworth & Kreidberg 2022, and references therein). The vast variety of exoplanets might arise from the distinct environments of circumstellar gas and dust in which planets form and evolve over millions, or even billions, of years. These environments undergo continuous transformations: beginning with the collapse of a molecular cloud that gives rise to a new star, progressing through a protoplanetary disk (PPD), where planets are born, and eventually becoming a debris disk as the star enters the main sequence after several million years. Studying these environments is essential to answering fundamental questions in exoplanet science.

The circumstellar material that provides the building blocks for future planets has different origins and properties in protoplanetary and debris disks. In PPDs, gas and dust are pristine, originating directly from the initial molecular cloud. In contrast, the primary mechanism of dust production in debris disks is the collisions between kilometer-sized rocky bodies. These collisions supply the disk and the forming planets with substantial amounts of dust grains of various sizes and small amounts of gas (e.g., Wyatt 2018; Hughes et al. 2018).

The evolution of dust particle properties from the protoplanetary to the debris disk phase can be studied using two distinct ranges of the electromagnetic spectrum, as dust particles interact with stellar light in two primary ways. Some stellar photons are scattered by dust grains in all directions, particularly at optical and near-infrared (IR) wavelengths. Meanwhile, other stellar photons are absorbed by the dust grains and re-emitted as thermal radiation predominantly in the IR to millimeter wavelength range.

In the past decades, space-based mid-IR and far-IR observations with the Spitzer, IRAS, and Herschel Space Observatory played a crucial role in advancing our understanding of debris disks. In particular, the Disc Emission via a Bias-free Reconnaissance in the IR/Submillimeter (DEBRIS; Sibthorpe et al. 2018) and DUst around NEarby Stars (DUNES; Eiroa et al. 2013) surveys provided comprehensive statistical studies of debris disks around nearby main-sequence stars. These surveys enabled precise measurements of IR excess and dust temperatures, revealing trends with stellar type and age. They also established a framework for estimating dust luminosity distributions and the incidence rate of debris disks, especially around solar-type and early-type stars.

In addition to mid-IR and far-IR observations, the Hubble Space Telescope (HST) made a groundbreaking contribution to the imaging of debris disks in scattered light. Using its high-contrast imaging (HCI) capabilities, HST provided the first resolved views of numerous debris disks, revealing their morphology and fine structures such as rings, warps, and asymmetries (e.g., Golimowski et al. 2006; Kalas et al. 2007). Systematic surveys of circumstellar environments led by Schneider et al. (2014, 2016) offered valuable complementary insights to thermal emission data, helping to constrain disk geometries and the scattering properties of dust grains.

The detailed studies of both scattered and thermal light from circumstellar disks using the ground-based telescopes became possible with the start of operation of high-contrast and high-resolution instruments such as the Spectro-Polarimetic High contrast imager for Exoplanets REsearch (SPHERE; Beuzit et al. 2019) at VLT, the Gemini Planet Imager (GPI; Macintosh et al. 2014) or the Subaru Coronagraphic Extreme Adaptive Optics (SCExAO; Jovanovic et al. 2015), along with interferometric facilities like the Atacama Large (sub)Millimeter Array (ALMA) which delivered unprecedented images of many protoplanetary and debris disks around young stars (e.g., Perrot et al. 2016; Andrews et al. 2018; Avenhaus et al. 2018; Boccaletti et al. 2020; Columba et al. 2024). These observations have targeted both individual disks (e.g., Garufi et al. 2016; Milli et al. 2017b; Olofsson et al. 2018; Ménard et al. 2020) and large disk samples (e.g., Ansdell et al. 2017; Ginski et al. 2024; Garufi et al. 2024; Matrà et al. 2025), facilitating the first demographic studies that address the morphology of these objects.

Most optical, near-IR and (sub)millimeter imaging campaigns have focused on studying PPD evolution and searching for forming planets within them (Benisty et al. 2023, and references therein). In contrast, only a few comparable studies have investigated direct imaging (DI) data for a large sample of debris disks (e.g., Schneider et al. 2014; Esposito et al. 2020; Crotts et al. 2024). One key reason is that debris disks are significantly older (≳7 Myr) and contain roughly three orders of magnitude less dust than PPDs (Wyatt 2008). As a result, they are much fainter and more challenging to image directly.

This study focuses on SPHERE observations of debris disks. SPHERE is an extreme adaptive optics (AO) instrument optimized for observing circumstellar environments. Since its commissioning in 2014, SPHERE has been extensively utilized and has proven to be one of the most productive HCI instruments. As part of the SPHERE guaranteed time observation (GTO) program, numerous debris disks have been observed and detected in the course of the dedicated disk program, and sometimes as a by-product of the SpHere INfrared survey for Exoplanets (SHINE; Chauvin et al. 2017; Desidera et al. 2021; Vigan et al. 2021; Langlois et al. 2021). To perform a comprehensive analysis of these observations, we compiled a sample of targets known to host debris disks from the archival datasets of GTO and various open-time programs, including all targets from the SPHERE High Angular Resolution Debris Disk Survey (SHARDDS, PI: J. Milli; Milli et al. 2017b; Dahlqvist et al. 2022).

This study aims to consistently characterize the structural and compositional properties of debris disks, focusing on both their radial and vertical extents, as well as the nature of their constituent dust. A primary goal is to explore how the architecture of debris disks, particularly the radial locations of planetesimal belts and their dust masses, which are key to understanding the evolution of debris disks and their interaction with planets (Krivov & Wyatt 2021), relates to the fundamental properties of their host stars. By analyzing a large sample of spatially resolved debris disks, we investigated how the belt radii and disk dust masses scale with stellar luminosity and mass, and how these relationships evolve over time. Additionally, we assessed the conditions that distinguish detected from non-detected disks in scattered light, accounting for both intrinsic disk properties and observational biases. We further investigated the architecture of planetary systems within the sample, analyzing how the presence, absence, or configuration of planetary companions correlates with the structure and detectability of debris disks.

The paper is organized as follows. In Sect. 2, we present the stellar parameters of the targets included in our sample. Section 3 describes the SPHERE observing modes and the data reduction techniques used throughout the study. In Sect. 4, we analyze the morphological structure of the detected disks, emphasizing the radial locations of the planetesimal belts. These constraints are then incorporated into spectral energy distribution (SED) modeling in Sect. 5, allowing us to derive the parameters of the dust grain size distributions (SDs) and to estimate plausible ranges for the scattering albedo.

A key part of our investigation, detailed in Sect. 6, addresses the question of why the majority of young debris disks remain undetected in scattered light imaging. In particular, we examine the role of the dust’s optical properties and the observational biases associated with viewing geometry and disk structure. Special focus is placed on the evaluation methods for the dust albedo and polarization efficiency based on polarimetric imaging (Sect. 6.3).

In Sect. 7, we shift focus to planetary system architectures. We analyze systems where both exoplanets and debris disks are detected, as well as those with no detected planets, by estimating the locations and masses of planets that could dynamically shape the observed disk structures. This includes modeling scenarios in which unseen planets are responsible for clearing gaps or truncating the inner edges of planetesimal belts. The key findings of the study are summarized in Sect. 8.

2 Sample description

We compiled a sample of debris disks from archival SPHERE observations, selecting main-sequence stars with an IR excess above 10⁻⁶, based on data from the Jena debris disk database¹. This sample comprises 161 stars spanning a broad range of spectral types and ages. The majority are young F-type (35%) and A-type (29%) stars (Fig. 1). These targets usually exhibit strong IR excesses, indicating significant amounts of dust, and are bright enough to serve as reference stars themselves for the AO system. For instance, ZIMPOL observations require a reference star with a G magnitude brighter than 9.5^m for a good AO correction in the optical. Additionally, the two surveys that contributed most to our sample, SHINE and SHARDDS, primarily focus on A-type and solar-type stars. The presence of an IR excess or debris disk was not a part of the selection criteria for the SHINE statistical sample (Desidera et al. 2021). However, in specific cases, the known presence of exoplanets or a disk, suggesting a higher likelihood of harboring young, directly imageable planets (e.g., Meshkat et al. 2017), led to a target being classified as a special object, thereby increasing its observational priority. In contrast, SHARDDS was a dedicated debris disk survey, with targets selected based on the predicted brightness of their disks (f_disk > 10⁻⁴)². The SHARDDS survey included 55 main-sequence stars observable from the Southern hemisphere, covering spectral types A through M and stellar ages ranging from 10 Myr to 6 Gyr. Its aim was to provide a comprehensive overview of planetary system properties and their temporal evolution.

In addition to A- and solar-type stars, our sample includes eight B-type and eight M-type stars, with the latter group exhibiting the highest detection rate among all spectral types in our sample. However, this high detection rate is in part due to the unexpected discovery of a debris disk around the M1Ve star GSC 7396-0759 (Sissa et al. 2018), which was not previously known to exhibit an IR excess and was routinely observed within the SHINE program. The median stellar mass of the full sample is 1.43 M_⊙, while for the subsample of targets with detected disks it is slightly lower at 1.38 M_⊙ (Fig. 1).

The sample comprises 18 binary systems, including spectroscopic, visual, and astrometric binaries, as well as spectroscopic binary candidates (SBCs), four triple systems, three quadruple systems and two systems with higher-order multiplicity (N > 4) according to the Washington Double Star catalog (WDS; Mason et al. 2001) as of January 15, 2024. In two of the triple systems, debris disks are known around two different components, which were observed individually and listed separately in Table G.2. These include HD 216956 (Fomalhaut A) and GSC 06964-1226 (Fomalhaut C), as well as HD 181296 (A component), which shares a common proper motion with HD 181327 (B component). Among the quadruple systems, HD 20320 consists of a spectroscopic binary (SB) as its A component and an astrometric binary as its B component, while HD 98800 features a pair of SBs orbiting each other (Kennedy et al. 2019). Another quadruple star system in the sample is HD 102647 (Denebola). The components of multiple systems that host debris disks and were observed with SPHERE are specified in Col. 6 of Table G.2.

Our sample includes five chemically peculiar stars classified as Lambda Boo stars: HD 30422, HD 31295, HD 110411, HD 183324, and HD 218396. These stars exhibit surface deficiencies in iron-peak elements while maintaining nearly solar abundances of carbon, nitrogen, oxygen, and sulfur (e.g., Paunzen 2001; Gray et al. 2017). This anomaly may be explained by preferential gas accretion over dust from a dynamically evolving debris disk, possibly influenced by migrating planets or accretion from the atmospheres of hot Jupiters (Murphy & Paunzen 2017). The debris disk hypothesis is further supported by the high fraction (up to 77%) of Lambda Boo stars exhibiting IR excess, which is often linked to the presence of a debris disk. (Draper et al. 2016b). Some of these disks have been imaged with the Herschel Space Observatory at 70,100 and 160 μm, including several debris disks analyzed in this study (Su et al. 2009; Draper et al. 2016b). In Sect. 4, we present a scattered light image of the inner belt surrounding a Lambda Boo star HD 218396 (HR 8799).

Stellar ages were compiled from the literature, with their lower and upper boundaries listed in Col. 11 of Table G.2. For some targets, particularly field stars, there are significant discrepancies, up to 3000 Myr, between ages reported in different studies. This large scatter arises from the use of diverse age-dating techniques, such as isochrone fitting, kinematic group membership, and indicators of stellar activity (e.g., Ca II H and K line strength or X-ray luminosity). For instance, published age estimates for HD 15115 include ${12}_{-4}^{+8}$ Myr (Moór et al. 2006), 100 Myr (Zuckerman & Song 2004), or ${500}_{-500}^{+1500}$ Myr (Holmberg et al. 2009).

In such cases, we adopted an age range that covers the full span of results derived from various methods. For bona fide members of moving groups (MGs) and targets lacking literature age estimates, upper and lower age limits were assigned based on the most probable MG membership. In Column 12 of Table G.2, we list the MG with the highest probability of association for each star, along with the corresponding probability percentage (in parentheses), as determined using the BANYAN ∑ tool (Gagné et al. 2018). According to this analysis, the majority of our sample consists of field stars (50%), followed by members of the β Pictoris MG (βPMG, 9%).

The median age of our sample is 100 Myr, with approximately half of the targets estimated to be between 10 and 100 Myr old (Fig. 1). The debris disks around the youngest stars (<10 Myr) such as Herbig Ae/Be stars HD 141569 (e.g., Perrot et al. 2016) and HD 156623, as well as T Tauri stars such as TWA 7 (e.g., Olofsson et al. 2018; Ren et al. 2021), exhibit structures with multiple rings and spiral arms (Figs. D.1 and B.1), features typically associated with PPDs. The fractional IR luminosities of these young stellar objects are generally below 0.1, leading to their classification as debris disks. However, these systems may represent an intermediate stage between the protoplanetary and debris disk phases. We categorize such disks as transition disks due to their evolutionary status.

Furthermore, transition disks often contain high CO masses, comparable to those found in PPDs (Lieman-Sifry et al. 2016; Moór et al. 2017, 2019). Our sample includes several debris disk systems with a significant gas reservoir, commonly referred to as hybrid disks: HD 9672 (Moór et al. 2011; Choquet et al. 2017; Pawellek et al. 2019), HD 21997 (Kóspál et al. 2013), HD 121617 (Perrot et al. 2023), HD 131488 (Pawellek et al. 2024), HD 131835 (Hung et al. 2015; Feldt et al. 2017), and HD 141569 (Dent et al. 2005). In these hybrid systems, dust evolution may have progressed more rapidly than gas dissipation (Péricaud et al. 2017).

Similar to stellar ages, a wide range of metallicity values for the same star can be found in the literature. Depending on the method used to determine metallicity, discrepancies of up to 0.5 dex can arise between different studies. To ensure consistency, we opted to use the median metallicity value from all studies recorded in the SIMBAD database³.

The target distances, listed in Col. 7 of Table G.2, were derived using stellar parallaxes from the Gaia DR3 catalog (Gaia Collaboration 2023b). The closest object in our sample is HD 22049 (ϵ Eridani), located at 3.6 pc from the Sun, while the most distant target is HD 149914 at 154.3 pc. Two-thirds of all targets in the sample are located within 80 pc of the Sun, with a minor peak at ~100 pc, corresponding to the distance of the Scorpius-Centaurus OB association, a region rich in young stars (Fig. 1).

3 SPHERE observing modes and data reduction

The disk observations presented in this work were performed with different SPHERE subsystems (see Table 1 and the SPHERE User Manual⁴): the InfraRed Dual-beam Imager and Spectrograph (IRDIS; Dohlen et al. 2008), the Integral Field Spectrograph (IFS; Claudi et al. 2008) and the Zurich Imaging POLarimeter (ZIMPOL; Schmid et al. 2018). A variety of instrument modes were used, including pupil- and field-stabilized configurations, along with different filters ranging from optical to near-IR. Observations were conducted with classical or apodized pupil Lyot coronagraphs (Boccaletti et al. 2008; Carbillet et al. 2011; Guerri et al. 2011), or in some cases, without a coronagraph.

Many disks were observed only once using either classical or polarimetric imaging (de Boer et al. 2020; van Holstein et al. 2020) modes of IRDIS with the broadband H filter (λ_c = 1.625 μm, Δλ = 0.290 μm), or polarimetric imaging mode of ZIMPOL (Schmid et al. 2018), often employing the Very Broad Band filter (VBB or RI; λ_c = 0.735 μm, Δλ = 0.290 μm). The observations of all SHINE targets were performed in either the IRDIFS or IRDIFS_EXT modes (Langlois et al. 2021), which provide a simultaneous data acquisition with both IRDIS and IFS. With these instrument setups, the IRDIS is operated in the dual-band imaging mode (Vigan et al. 2010) with the filter pair H2H3 (λ_H2 = 1.593 μm, Δλ_H2 = 0.052 μm; λ_H3 = 1.667 μm, Δλ_H3 = 0.053 μm) for the IRDIFS mode, or with the filter pair K1K2 (λ_K1 = 2.110 μm, Δλ_K1 = 0.102 μm; λ_K2 = 2.251 μm, Δλ_K2 = 0.109 μm) for the IRDIFS_EXT mode, whereas the IFS is operated in the IRDIFS Y-J mode (0.95–1.35 μm, with a spectral resolution of R_λ = 50), or IRDIFS_EXT Y-H mode (0.95–1.65 μm, with a spectral resolution of R_λ = 35).

The IRDIS and IFS datasets were processed at the High-Contrast Data Center⁵ (HC-DC, Delorme et al. 2017, formerly known as the SPHERE Data Center). For both instruments, the pre-processing steps are based on the SPHERE Data Reduction and Handling pipeline (Pavlov et al. 2008) to correct for bad pixels, flat-field non-uniformity, optical distortions, and telescope or sky background. In addition, for the IFS, the pre-processing includes a wavelength calibration and a correction for cross-talks between spectral channels. Coronagraphic images are centered via four satellite spots used to determine the accurate position of the star hidden behind the coronagraphic mask.

Pre-processed IRDIS and IFS datasets form spectral and temporal cubes of centered images, to which dedicated stellar subtraction algorithms can be applied. Such algorithms include classical Angular Differential Imaging (ADI; Marois et al. 2006), Principal Component Analysis (PCA; Soummer et al. 2012; Amara & Quanz 2012) or the Locally Optimized Combination of Images (LOCI; Lafrenière et al. 2007), implemented in the HC-DC in a template-oriented version (T-LOCI; Marois et al. 2014; Galicher et al. 2018). For several datasets, post-processing employing the reference-star differential imaging (RDI) technique was also applied, as described in Xie et al. (2022).

The IRDIS polarimetric datasets were processed using the IRDAP pipeline (van Holstein et al. 2020), while the ZIMPOL polarimetric datasets were reduced with a pipeline developed at ETH Zürich, as described in Engler et al. (2017) and Hunziker et al. (2020). Both pipelines are currently implemented in the HC-DC and include, as part of the pre-processing steps, subtraction of bias and dark frames, flat-fielding, and correction for instrumental polarization. Additionally, the ZIMPOL frames are corrected for modulation and demodulation efficiency (Schmid et al. 2018).

In both pipelines, the Stokes parameter Q and U images are computed from the calibrated and centered polarimetric frames using the double-difference method. These Q and U images are then transformed into the azimuthal Stokes parameter Q_φ and U_φ: $$\begin{array}{rl}& Q_{\phi }=-Q\cos 2\phi -U\sin 2\phi \\ & U_{\phi }=Q\sin 2\phi -U\cos 2\phi ,\end{array}$$

where φ is the polar angle measured east of north in a coordinate system centered on the star, and the sign convention Q_φ = −Q_r and U_φ = −U_r is adopted from Schmid et al. (2006) (see also Monnier et al. 2019).

4 Morphology of resolved debris belts

Out of 161 targets, 51 debris disks were successfully detected with SPHERE in total intensity of scattered light (Fig. 2, all images are presented in Fig. D.1), linearly polarized intensity⁶ (Fig. 3, all images are presented in Fig. F.1), or both. Four of these debris disks, BD-20 951 (Perrot et al., in prep.), HD 36968 and the inner belts of HD 218396 (HR 8799) and HD 36546 systems had not been imaged with any instrument before. Debris disks HD 38206, HD 36546, HD 38397 (Perrot et al., in prep.), HD 98800, HD 182681 are resolved in scattered light for the first time. Scattered light images of the HD 105, HD 377, TWA 25 (Langlois et al., in prep.), HD 30447, HD 92945, HD 145560, HD 192758 and HD 202917 debris disks have previously been obtained with instruments such as HST or GPI. However, the SPHERE images of these disks have not yet been published. The majority of detections were around F-type stars, with 23 discoveries, corresponding to a 45% detection rate in our sample (Fig. 1 upper left panel).

We determined the radii of the resolved debris belts by analyzing the r²−scaled images of total or polarized intensities. The radial position of the peak surface brightness (SB) along the disk’s major axis was measured, and the resulting disk radius, referred to as $R_{\text{belt}}^{\text{mes}}$, is listed in Col. 2 of Table 2. The inclination and position angle (PA) of each disk, provided in Table 2, were derived by fitting ellipses to the visible contours of the disk rims. Additionally, disk images with a higher signal-to-noise ratio (S/N) were modeled in more detail to obtain the fundamental geometrical and scattering parameters necessary for a comprehensive characterization of disk properties (Sect. 4.5.1).

In all images presented in Figs. 2, 3, D.1, and F.1, sky north is up and east to the left. The PA defines the orientation of the disk’s major axis on the sky and is measured counterclockwise from sky north to east. The PA values of the eastern disk extensions are listed in Col. 4 of Table 2 (0° ⩽ PA ⩽ 180°). The inclination of a debris disk is conventionally defined as the angle between the sky plane and the disk’s minor axis, where a pole-on disk has an inclination of 0°, and an edge-on disk has an inclination of 90°. In this study, disk inclinations (Col. 3 of Table 2) follow the convention that an inclination is less than 90° when the brighter side of the disk is oriented southward, whereas an inclination is greater than 90° when the brighter side is oriented northward. In many debris disk studies, the inclination is reported as less than 90°, regardless of the orientation of the disk’s brighter side. To enable comparison with such studies, the inclinations of disks with the brighter side oriented northward (i > 90° in Table 2) can be converted using i_<90° = 180° − i_>90°.

4.1 Disk radii in SPHERE versus ALMA observations

We detect planetesimal belts in 33 debris disks which have also been resolved with ALMA and SMA at wavelengths of 0.856–1.34 mm as part of the REASONS survey (Matrà et al. 2025). The nature of dust emission observed in scattered light images (from optical to near-IR wavelengths) and thermal emission images (from mid-IR to millimeter wavelengths) is fundamentally different. In SPHERE images (both total and polarized intensities), we observe stellar photons scattered off dust grains into our line of sight. In contrast, the thermal emission detected by ALMA and SMA originates from the absorption of stellar photons by dust particles, which raises their temperature and leads to re-emission at longer wavelengths. Additionally, scattered light images trace predominantly dust particles with sizes smaller than a few microns, whereas (sub)millimeter imaging is more sensitive to (sub)millimeter particles. As a result, the disk morphology, particularly the radial position and extent of the belt, can differ between scattered-light and thermal-emission images.

Since the spatial resolution of many millimeter observations is sufficiently high to examine the relationship between belt radii measured in both near-IR and millimeter wavelengths, we analyze this correlation and present our results in Fig. 4. For this comparison, we used disk radii measured from r²−scaled SPHERE images (Col. 2 in Table 2), while the belt radii observed in thermal emission were obtained from the REASONS survey (Tables 1 and A.1 in Matrà et al. 2025). In that study, all targets are fitted with a single planetesimal belt model, where the radial surface density of particles is described by a Gaussian distribution. Consequently, the derived belt radii represent the centroid radii of this distribution.

To ensure consistency, we excluded from this comparison the REASONS targets that were only marginally resolved in millimeter observations or exhibited more than one planetesimal belt, with two exceptions:

• HD 15115: the radial locations of its two cold belts were taken from the two-belt model fit of the ALMA image presented by MacGregor et al. (2019).

• HD 92945: the SB profile of the disk, as shown in Fig. 2 by Marino et al. (2019), was used to determine the radial positions of its two belts in ALMA images.

Figure 4 demonstrates a good agreement between the belt radii measured from SPHERE and ALMA images, indicating a near 1:1 relationship between the locations of the radial SB peaks in near-IR scattered light and thermal emission images. A linear fit to the data (black solid line in Fig. 4) yields a slope of 1.05 ± 0.04, representing the ratio $R_{\text{belt}}^{\text{mes}}(near-IR)/R_{\text{belt}}(\mathrm{m}\mathrm{m})$. This value is lower than the average ratio of 1.39 reported by Esposito et al. (2020) in a similar comparison. This finding highlights the need for higher-sensitivity and higher-resolution observations to better understand the connection between disk structures observed in scattered light and thermal emission.

4.2 Ratio of radii in multiple belt systems

Observations across multiple wavelengths, from optical to millimeter, suggest that many young debris disks likely consist of multiple planetesimal belts (e.g., Golimowski et al. 2006; Bonnefoy et al. 2017; Marino et al. 2019). This is further supported by the fact that many disk SEDs are better modeled using two blackbody (BB) components with distinct equilibrium temperatures (T_bb), requiring dust populations at different radial distances from the host star⁷ (e.g., Chen et al. 2014).

Multiple belt configurations are rarely detected in DI, as disks with high inclinations are more easily resolved, as discussed in Sect. 6.2. Among the debris disks imaged with SPHERE, excluding transition disks such as HD 141569, seven systems exhibit spatially resolved double-belt structures, while HD 131835 (Feldt et al. 2017) and TWA 7 images reveal three distinct planetesimal belts. In these systems, both the inner and outer belts belong to the category of cold exo-Kuiper belts (Sect. 4.3) and are listed separately in Table 2. Interestingly, the ratio between the outer and inner belt radii is consistently around 1.5 or 2, with a median value of 1.69 across the nine resolved systems (Fig. 5). These similar belt spacing ratios may hint at similar evolutionary pathways of debris systems or could indicate the presence of mean-motion resonances, possibly due to unseen planets shaping these structures.

4.3 Empirical correlation between belt radius and star luminosity

As mentioned in Sect. 4.2, the SEDs of many debris disks require a two-BB model fit. In such two-temperature debris disks, the dust populations are typically classified into warm dust belts (belts with BB temperature between ~100 and 200 K) and cold dust belts (exo-Kuiper belts with BB temperature below 100 K). This bi-modal temperature distribution has been investigated in previous studies (e.g., Kennedy & Wyatt 2014) and is often explained by the preferential formation of planetesimals at the ice lines of volatile compounds such as water, ammonia, carbon dioxide, or carbon monoxide (Morales et al. 2011).

The ice line, also referred to as the frost or snow line, of a volatile compound marks the minimum radial distance from a star at which the temperature is sufficiently low for the compound to condense. Beyond this distance, gas condensation promotes the formation and growth of icy dust particles, thereby facilitating the development of planetesimals. Consequently, debris belts may be more likely to form just beyond the ice lines of common volatile compounds such as H₂O, CO or CO₂.

The positions of ice lines within a disk are not fixed throughout a star’s lifetime, as they are influenced by the evolving stellar luminosity and the opacity of surrounding material. As a result, the disk’s radial temperature profile and the condensation thresholds of different volatile substances change over time. This implies that the range of radial distances at which a specific volatile compound may condense into ice can be relatively broad. For example, in the solar nebula, the water snow line has been predicted to lie at 2.7–3.2 au, with grain temperatures between 170 and 143 K, depending on the model (Hayashi 1981; Podolak & Zucker 2004), In contrast, the current water snow line in the Solar System is estimated to be at ~5 au from the Sun (Jewitt et al. 2007). Moreover, the condensation temperature of volatile compounds is influenced by the properties of debris particles onto which the gases freeze. Kim et al. (2019), for instance, found that in the case of the young A-type star β Pic (HD 39060), the water snow line could be located anywhere between 4.4 and 28.3 au, depending on the dust grain composition, grain size and ice phase (amorphous or crystalline).

To explore the correlation between the locations of ice lines and planetesimal belts within the subsample of debris disks spatially resolved with SPHERE, we estimated the temperature of their BB grains (Col. 5 in Table 2). These grains, being significantly larger than the peak wavelength of the emitted disk spectrum, allow their temperature to be determined using the following expression (Backman & Paresce 1993): $$T_{\mathrm{b}\mathrm{b}}=(278\mathrm{K}){\left(\frac{L_{\star }}{L_{\odot }}\right)}^{0.25}{\left(\frac{1\mathrm{a}\mathrm{u}}{R_{\mathrm{b}\mathrm{e}\mathrm{l}\mathrm{t}}^{\mathrm{m}\mathrm{e}\mathrm{s}}}\right)}^{0.5},$$

where L_⋆ is the stellar luminosity, and $R_{\text{belt}}^{\text{mes}}$ is the measured belt radius in au (Col. 2 in Table 2).

According to this estimation, all resolved debris belts fall into the category of exo-Kuiper belts containing cold dust (T_bb < 100 K), with the exception of the warm dust belt around HD 172555, which was detected with ZIMPOL (Engler et al. 2018). In Fig. 6, we present the derived BB temperatures of the belts as a function of their measured radii. The shaded regions in the plot indicate the upper temperature limits at which H₂O, CO₂ and CO may condense in young disks, depending on gas pressure and dust temperature (Harsono et al. 2015). For comparison, the plot also includes the locations of the Edgeworth–Kuiper belt at 40 au (Stern & Colwell 1997), with an estimated BB temperature of T_bb = 44 K, and the main asteroid belt in the Solar System at 3.5 au (Wyatt 2008), with T_bb = 150 K.

As shown in Fig. 6, the majority of planetesimal belts are located within the CO₂ and CO ice formation zones. This trend is also evident in seven double-belt systems, where both components reside within the same ice-species region. Notably, all three resolved planetesimal belts of HD 131835 lie within the CO₂ ice zone, suggesting that they originate from a common, broad debris disk in which gaps may have been sculpted by planetary bodies. The disk around HD 172555 is located in a region where water molecules can accumulate on grain surfaces. This analysis supports the hypothesis that planetesimals preferentially form beyond the ice lines of various gas species.

If this statement holds true, the radial distance of a debris belt should correlate with the luminosity of its host star. This relationship has recently been examined in samples of debris disks resolved at millimeter and far-IR wavelengths (Matrà et al. 2018; Marshall et al. 2021). The subsample of debris disks spatially resolved with SPHERE (Table 2) provides an opportunity to explore the correlation further. To quantify this relationship, we applied a power-law fit $$R_{\text{belt\hspace{0.17em}}}=R_{L_{\odot }}{\left(\frac{L_{\star }}{L_{\odot }}\right)}^{\alpha }$$(1)

to the data points in Fig. 7, where we show the distribution of the exo-Kuiper belts in our subsample in the parameter space [$R_{\text{belt}}^{\text{mes}}$, L_⋆]. The scaling factor R_L⊙ is in au and represents the expected radial position of a planetesimal belt around a star with solar luminosity.

We obtained a relatively shallow linear dependence in logarithmic space log(R_belt) = α log(L_⋆/L_⊙) + log(R_L⊙) (magenta dash-dotted line in Fig. 7) with R_L⊙ = 74 ± 7 au and α = 0.11 ± 0.05. These values remain within the 1σ uncertainties of similar parameters reported in studies at millimeter and far-IR wavelengths (Matrà et al. 2018; Marshall et al. 2021).

The HD 172555 disk was excluded from this analysis, as it is the only system in our subsample that contains warm dust. However, it is likely that the disks included in our subsample formed in connection with the ice lines of various volatile species, such as CO₂ and CO gases. If this is the case, analyzing the relationship between the radial distance of a belt and stellar luminosity requires categorizing the sample based on disk BB temperature, which may correspond to the freeze-out temperature of a specific gas specie.

Therefore, we divided our subsample into two groups of disks based on their temperatures T_bb. Taking into account the uncertainties in the estimated temperature, we set T_bb = 35 K as the upper limit for the coldest disks in the subsample, where the CO gas may freeze out (CO subsample). By fitting a power-law function (Eq. (1)) to this group of disks, we obtained best-fit parameters of R_L⊙ = 96 ± 7 au and α = 0.30 ± 0.07. This fit is represented by the orange dashed line in Fig. 7. For the group of disks with a local equilibrium temperature above 35 K (CO₂ subsample), the best-fit parameters are R_L⊙ = 43 ± 8 au and α = 0.30 ± 0.08, as indicated by the blue solid line in Fig. 7.

As expected, the power-law functions for both disk groups are steeper than the function fitted to the entire sample. Notably, the parameter R_L⊙ for the CO₂ subsample is found to be 43 au. This radial distance closely corresponds to the location of the Edgeworth–Kuiper belt in the Solar System.

As previously discussed in this section, the radial locations of the ice lines for volatile compounds vary over the course of stellar evolution. To investigate whether a corresponding evolution in the radial positions of planetesimal belts is observable, we plotted the measured belt radii as a function of stellar age for the systems in the CO₂ subsample. This subsample provides a relatively larger number of systems with stars of similar spectral type but different ages, allowing for a more meaningful comparison.

Given that stellar luminosity is a key parameter in this context, we examined three groups of stars categorized by spectral type and mass: (1) A-type stars with 2.2 M_⊙ < M_⋆ < 2.6 M_⊙, (2) A-type stars with 1.9 M_⊙ < M_⋆ < 2.1 M_⊙, and (3) F-type stars with 1.2 M_⊙ < M_⋆ < 1.6 M_⊙. For stars with multiple resolved belts, we adopted the mean belt radius, as all detected belts in these systems have T_bb >35 K.

We note, that the mean estimated age (Col. 11 in Table G.2) of most stars with detected disks (90% of detections) is below 50 Myr. According to stellar evolutionary models (e.g., Palla & Stahler 1999; Baraffe et al. 2015), such young objects are likely located either on the pre-main-sequence (PMS) or on the zero-age main sequence (ZAMS) in the Hertzsprung-Russell diagram.

Intermediate-mass stars (1 M_⊙ < M_⋆ < 3 M_⊙) exhibit the most pronounced luminosity evolution during their PMS phase. These stars begin as fully convective objects with large radii and high luminosities, which decrease as the stars contract. This phase is followed by an increase in both temperature and luminosity as a radiative core begins to develop because it becomes hotter and denser during the contraction phase. This process leads to an increase in the rate of nuclear fusion, ultimately stabilizing the star and placing it on the main sequence.

The luminosity evolution during the PMS phase depends sensitively on the stellar mass and chemical composition (e.g., Tognelli et al. 2011). To illustrate this, Fig. 8 shows the evolution of the radial location of the CO₂ freeze-out zones (red, orange and blue shaded regions), corresponding to the CO₂ zone presented in Fig. 6, for stars with metallicity Z = 0.028 and helium abundance Y = 0.304, calculated for stellar masses of 1.5 M_⊙, 2 M_⊙ and 2.3 M_⊙ (Tognelli et al. 2011).

As shown in Fig. 8, stars in all three mass groups show a trend of increasing belt radius with stellar age within the region corresponding to the CO₂ freeze-out zone. This behavior may reflect the outward migration of the CO₂ ice line due to increasing luminosity as the star evolves. To better quantify this result, we performed a multiple regression analysis. We considered log t_age (in Myr) and log M_⋆ (in solar masses) as independent variables, and log R_belt (in au) as dependent variable. The best-fit relation is the following: $$\begin{array}{rl}\log \left(R_{\text{belt\hspace{0.17em}}}/\mathrm{a}\mathrm{u}\right)=& (0.37\pm 0.11)\log \left(t_{\text{age\hspace{0.17em}}}/\mathrm{M}\mathrm{y}\mathrm{r}\right)\\ & +(0.59\pm 0.23)\log \left(M/M_{\odot }\right)+(1.27\pm 0.15).\end{array}$$

Both coefficients are significant (with a 0.002 probability of being a chance result for the dependence on the age, and of 0.02 for the dependence on the mass). The dependence on age is then highly significant. This relation predicts log R_belt for stars in this range of ages and masses with an accuracy of 0.12 dex.

4.4 Morphology of selected targets

In this section, we discuss debris systems that have been imaged in scattered light for the first time, as well as debris disks whose morphology exhibits notable features, such as multiple belts, that warrant further examination.

HD 9672 (49 Ceti). The A1V star HD 9672 is one of the youngest and brightest stars in the sample. It is likely a member of the ~40 Myr old Argus MG (99% membership probability; Zuckerman 2018). The debris disk surrounding HD 9672 is gas-rich, with an estimated CO mass exceeding 2.2 × 10⁻⁴ M_⊕. The spatial distribution of CO gas closely resembles the structure of the outer debris disk, whereas no molecular gas has been detected within ~90 au of the star (Hughes et al. 2008).

The thermal emission of the HD 9672 disk is well characterized by two distinct dust populations: a warm component at 136–160 K and a cold component at 47–60 K (e.g., Wahhaj et al. 2007; Roberge et al. 2013; Chen et al. 2014, this work). Models reproducing the disk’s emission at λ = 12.5 μm and λ = 17.9 μm suggest that the warm dust grains are located within 60 au (Wahhaj et al. 2007).

The debris disk has been observed multiple times with SPHERE using different filters (e.g., Pawellek et al. 2019), including the first detection of its scattered light (Choquet et al. 2017). The PCA-reduced data taken with the IRDIS B_Y filter, shown in Fig. 9a, reveal a broad debris ring extending to the image edges (~6″). This structure may consist of multiple narrow rings. In the r²−scaled image, the radial SB peak is located between 144 and 156 au. Additionally, image residuals hint at an inner planetesimal ring with a radius of ~105–110 au (see Fig. 9a). If this second cold debris ring is confirmed, it would contain dust grains with a blackbody temperature of T_BB = 54 K and would not account for the warm dust excess observed in the HD 9672 SED.

Polarized scattered light observations of the disk, conducted with IRDIS in polarimetric mode using the B_Y filter, led to a detection, albeit with a relatively low S/N (Fig. F.1). This may be attributed to an intrinsically low polarization fraction of the dust in this disk.

HD 16743. The ALMA image of the debris disk around the F0 star HD 16743 was recently published by Marshall et al. (2023). While the authors report only a marginal detection of the disk in the IRDIS H-band data, our image, obtained with the IRDIS H23 filter (Fig. D.1), provides the first clear detection of this debris disk in scattered light. This detection allows us to constrain the disk’s geometrical parameters (Table 2). We determined that the PA of the disk is ~169°, which closely matches the value reported by Marshall et al. (2023). Additionally, our image reveals an extended feature at a PA of approx. 17° (see Fig. B.2). The origin of this feature remains unclear; it is most likely a residual PSF artifact, possibly caused by the telescope spider. However, the possibility of a scattered light signal cannot be entirely ruled out.

HD 36546. HD 36546 is another A-type star in our sample (A0V-A2V; Lisse et al. 2017; Currie et al. 2017), located at a distance of 100.1 pc. It is a probable member of the Mamajek 17 group, with an estimated age between 3 and 10 Myr. For the first time, its debris disk has been spatially resolved using the Subaru/HiCIAO camera in the H band (Currie et al. 2017).

Our observations reveal a well-defined debris ring in the IFS data at a radial separation of 55 ± 10 au, as shown in Fig. 9c. This ring is also visible in the IRDIS image in the K band (Fig. 9b) albeit with a lower S/N. Additionally, the IRDIS image, along with some IFS data, reveals a more extended debris belt with a radius of 110 ± 30 au (outer ring in Fig. 9b). This belt appears both wider and brighter than the inner ring and may consist of multiple components. The presence of two cold debris rings at approx. 55 and 110 au aligns with the modeling results of Currie et al. (2017) and Lawson et al. (2021), who determined that the debris disk of HD 36546 extends between 60 and 115 au. However, due to the system’s relatively high inclination (~80°), precisely determining the number of rings present remains challenging.

The residual pattern observed within the inner ring (Fig. 9b) resembles the structure of a smaller, distinct ring, particularly in its alignment along the major axis of the outer rings. If this feature represents a genuine ring rather than PSF residuals, an alternative that cannot be ruled out, its radial separation from the star would be ~30 au. Interestingly, Lisse et al. (2017) previously reported a debris belt at ~135 K, which is expected to be located between 20 and 40 au from HD 36546. Moreover, when fitting the disk spectrum and photometric data from multiple instruments, Lisse et al. (2017) also predicted the existence of an additional inner belt with a temperature of ~570 K, corresponding to hotter dust located between 1.2 and 2.2 au within the system.

HD 36968. HD 36968 is a young (~20 Myr) F2V star located at 149 pc in Octans Association (Moór et al. 2011; Murphy & Lawson 2014). The debris disk surrounding this star exhibits a high IR excess of 1.34 × 10⁻³ and has been detected for the first time in both total scattered intensity and polarized intensity (Figs. D.1 and F.1). The fundamental geometrical parameters of the disk are provided in Table 2.

HD 92945. HD 92945 is a nearby K1V star located at a distance of 21.51 pc from the Sun. We resolved the inner debris belt with a radius of ~56 au and a part of the outer belt at ~119 au (Fig. D.1), both of which were previously imaged with ALMA (Marino et al. 2019). The PCA reduction of the H-band data shows residual structures that suggest the presence of a third dust ring with a possible radius of ~38 au.

HD 98800. HD 98800 is an intriguing quadruple system located 42.1 pc from the Sun (Gaia Collaboration 2021; Boden et al. 2005). As a member of the TW Hydrae Association, it has an estimated age between 7 and 10 Myr (Ducourant et al. 2014). The system consists of two pairs of spectroscopic binaries which orbit each other with a semimajor axis of ~50 au and period of 246 years (Kennedy et al. 2019; Zúñiga-Fernández et al. 2021). A double-lined SB BaBb (M_Ba = 0.70 M_⊙, M_Bb = 0.58 M_⊙, P = 315 days; Boden et al. 2005) is surrounded by a bright transitional debris disk, previously imaged at 1.3 mm with ALMA (Kennedy et al. 2019) and at 8.8 mm and 5 cm with the Very Large Array (VLA) (Ribas et al. 2018).

Imaging this disk in scattered light presents a significant challenge due to its small angular size and the presence of a close central binary. The disk has a radius of less than 0.1″ and an inclination of less than 45°, making non-coronagraphic polarimetric differential imaging (PDI) with ZIMPOL the most suitable observing strategy for resolving it in scattered polarized light. Such observations were conducted on April 14, 2016 (ESO 097.C-0344, PI: Kennedy), and our data reduction successfully detected scattered light from the disk. Figure 3 presents the Q_ϕ image obtained using the ZIMPOL R_PRIME filter. The HD 98800 transitional disk is the smallest detected with SPHERE/ZIMPOL in both angular and physical size (R_belt ≈ 0.07″ or 3 au). In Fig. 3, the white bar in the HD 98800 panel showing the Q_ϕ image represents 0.5″ (20 au), whereas in all other panels it corresponds to 1″.

The Stokes Q and U signals, used to compute the Q_ϕ image, are partially reduced due to the significant PSF convolution effect caused by the small angular size of the disk (Engler et al. 2018). Additionally, the disk signal may be affected by residual flux from the central binary, as complete removal of stellar flux in the center of image is not possible even with PDI. These stellar residuals are always present at the image center and originate either from nonzero polarization of the star(s) or from slight mismatches in the PSF shapes of the two orthogonal polarization states, which do not perfectly align. These mismatches arise due to short coherence times and, particularly relevant for a binary system, differences between the PSF of a binary star and that of a single point source. For the HD 98800 disk, we estimated the extent of the stellar residuals from the BaBb binary in the Q_ϕ image by analyzing the residuals from the AaAb binary, which is located at ~0.7″ from BaBb, just outside of the frame in Fig. 3. These stellar residuals have been masked in the central region of the presented image.

The scattered polarized light from the disk is detected between 0.06″ (2.52 au) and 0.12″ (5.06 au). Two dips in SB are visible on the eastern (PA = 109°) and western (PA = 280°) sides of the disk, which may be attributed to stellar PSF effects. Based on the SB distribution in the Q_ϕ image, we derived the geometrical parameters of the disk, as listed in Table 2. Within uncertainties, these parameters are in good agreement with those obtained from VLA and ALMA images (Ribas et al. 2018; Kennedy et al. 2019).

HD 111520. The strong brightness asymmetry in the scattered light of the nearly edge-on disk around HD 111520 (F5/6V star at d = 108 pc) has been previously observed with HST (Padgett & Stapelfeldt 2016) and GPI (Draper et al. 2016a). In the SPHERE/IFS image, the northern extension of the debris disk appears significantly brighter than the southern extension, where a dip in SB is observed at approx. 0.5″. Within 0.8″, the disk morphology resembles that of AU Mic disk. The SB variations along the major axis in both systems may be explained by the presence of a spiral disk structure or a set of non-coplanar debris rings. Indeed, the SED of HD 111520 is best fitted with multiple dust populations at different temperatures, suggesting the existence of radially separated debris belts containing both warm and cold dust.

HD 120326. The two distinct cold dust belts around the F0V star HD 120326 were first resolved in scattered light with SPHERE (Bonnefoy et al. 2017). We measure a radial distance of ~119 au (1.05") for the larger planetesimal belt and ~50 au (0.44″) for the smaller one. The polarimetric data of HD 120326 reveal polarized light between 0.25″ and 0.7″ with a tentative SB peak at ~0.5″, potentially indicating the presence of an additional inner debris belt in this system.

HD 129590. HD 129590 is a G3V star with one of the highest IR excesses in our sample (f_disk = (6.3 ± 1.8) × 10⁻³). The star is surrounded by two planetesimal belts, forming a structure reminiscent of a “moth” shape (Matthews et al. 2017; Olofsson et al. 2023) similar to that observable in the HD 61005 disk (Buenzli et al. 2010). The inner belt, located at ~49 au, is bright and exhibits an extended halo of small dust particles. In contrast, the outer planetesimal ring is significantly fainter but remains clearly visible in the PCA-reduced total intensity images (Fig. 9d). The region between the two belts does not appear to be completely cleared. The polarized intensity data show that although the outer ring is less pronounced in the halo, it remains detectable (Fig. 9e). This ring likely extends between 80 and 92 au, with a peak SB measured at ~82 au in the r²−scaled polarized intensity image. Recently, CO gas was detected in the system (Kral et al. 2020), supporting the possibility of gas pileup as a contributing factor to the observed disk structure (Olofsson et al. 2023).

HD 145560. We resolve the debris disk around HD 145560 (F5V star at 121.23 pc) in total intensity using the RDI technique (Xie et al. 2022) applied to the H2H3 dataset taken with IRDIS (Fig. D.1), as well as in polarized intensity using the VBB filter of ZIMPOL (Fig. 9h). This disk has also been observed with ALMA (Lieman-Sifry et al. 2016; Matrà et al. 2025) and GPI (Esposito et al. 2020). Among all available data for this target, the ZIMPOL image provides the highest spatial resolution and appears to reveal a spiral-like structure on the southern side of the disk, as well as a point-source-like residual (denoted as “PS?” in Fig. 9h) on the western side. However, this image is affected by low-wind effects and a short coherence time during the observation, which lowered the S/N of the polarimetric data, making the detection of this structure uncertain. The apparent point source could, in reality, be a bright part of the disk.

The scattered light in both SPHERE images exhibits an elliptical structure, with a major axis PA = 39 ± 5.0° and a radius of r = 87 ± 5 au, as measured from the r²−scaled image, and an inclination of 48.7 ± 7.0°. These geometrical parameters are within 1σ in good agreement with the GPI measurement (Esposito et al. 2020). However, there is a noticeable offset when comparing the disk’s orientation on the sky as measured with ALMA: 20 ± 7.0° at 1.24 mm (Lieman-Sifry et al. 2016) and 28 ± 8.0° at 1.3 mm (Matrà et al. 2025). This discrepancy may be due to the different spatial and angular resolutions between the scattered light images (SPHERE, GPI) and the thermal emission images (ALMA). Since the HD 145560 disk is the only resolved disk in our sample that shows such a PA deviation compared to ALMA data, this offset may suggest a more complex disk structure than a simple ring. Possible explanations include a spiral structure or the presence of multiple planetesimal rings at different PAs and inclinations.

HD 157587. The debris disk around F5V star HD 157587 is resolved with SPHERE instruments in both total and polarized scattered intensities. In the Q_ϕ image taken in the broadband H (Fig. 9f), the residual pattern inside the disk resembles the morphology of a smaller ring with a radius of ~50 au. These suspicious residuals are particular visible within the southeast extension of the disk and are also present in the Stokes U image.

HD 182681. HD 182681 is a B8.5V star located at 70.69 pc and a member of the βPMG. The debris disk surrounding this star was recently resolved with ALMA at 1.27 mm (Matrà et al. 2025). In the IRDIS H-band image, we detected extended scattered light emission from the debris belt, which has a radius of ~2.27″ (160 au). Additionally, there is evidence of a possible second belt at a radial distance of ~2.94″ (208 au).

HD 218396 (HR 8799). HD 218396, better known as HR 8799, is classified as an F0+VkA5mA5 Lambda Boo star (Gray et al. 2006) and is located at a distance of ~41 pc. It is surrounded by an extended exo-Kuiper belt, previously imaged at far-IR, sub-millimeter, and millimeter wavelengths using various facilities (e.g., Hughes et al. 2011; Faramaz et al. 2021, and references therein). Within the large disk cavity (r ~ 100 au), four giant planets with masses below 10 M_Jup have been discovered (Marois et al. 2008, 2010) and extensively studied (e.g., Esposito et al. 2013; Zurlo et al. 2016; Wang et al. 2018, see also Sect. 7). The relatively low inclination (~30°) of the debris belt facilitates planet detection but makes the belt itself challenging to observe using DI with the ADI technique.

HD 218396 was observed in various modes with all SPHERE instruments. In the imaging modes, the cold debris belt, extending between ~80 and 350 au (2″−7.5″) and peaking in SB at 180–200 au (~4.5″, Faramaz et al. 2021), was not detected. Similarly, in the H-band polarized intensity image, the disk remains either undetectable or barely visible, likely due to its low SB in polarized light. However, the image reveals a bright ring-like structure near the coronagraph, at a radial distance of ~0.4″ or ~15 au (Fig. 9i).

The possibility that this structure results from stellar PSF residuals cannot be entirely excluded. HD 218396 was observed with IRDIS in polarimetric mode on two nights: October 11 and October 13, 2016. The ring-like structure is detected in the data from October 11, when the observing conditions were significantly better (seeing between 0.44″ and 0.78″, and coherence times between 3.5 and 6.1 ms) compared to those on October 13 (seeing between 0.92″ and 2.82″, and coherence times between 1.7 and 4.0 ms). Poor observing conditions, such as those during the second night, can completely prevent the detection of a debris disk (see discussion in Sect. 6). Therefore, the non-detection of the ring in the data from the second night does not rule out the presence of a warm planetesimal belt at the considered radial position. Additionally, PSF residuals of this kind, especially at locations farther from the coronagraph and AO ring, are uncommon in IRDIS polarimetric data. Thus, the imaged ring likely traces polarized scattered light from dust particles in a second, inner debris belt.

The idea that this feature originates from a warm dust belt is supported by flux measurements of HD 218396 obtained with IRAS, ISO, and Spitzer space telescopes. Based on these data, Su et al. (2009) modeled the disk SED with three distinct dust components: a warm belt, a cold belt, and an extended halo. Stronger evidence for the presence of a warm belt at ~15 au comes from recent JWST/MIRI observations of HD 218396 at mid-IR wavelengths (Boccaletti et al. 2024) which provided spatially resolved signatures of the inner disk component.

If these interpretations are correct, the IRDIS H-band image resolves the inner warm dust belt in scattered polarized light for the first time. This belt has a radial distance of r = 15.5 ± 1.8 au, an inclination of i = (32 ± 7)°, and a PA of (39 ± 22)°, and it may have an offset from the central star.

TWA 7. TWA 7 is a 4.4 ± 1.4 Myr old (Herczeg & Hillenbrand 2014) M2Ve star in the TW Hydrae association. The IRDIS polarimetric H-band image (Fig. 9g) reveals a nearly pole-on system consisting of three rings at approx. 27, 52 and 93 au (Olofsson et al. 2018; Ren et al. 2021). The disk in the region between the inner and middle rings exhibits a clumpy structure, with arc-like streamers particularly evident in the area extending from the middle ring to the outer boundary of the detected scattered-light emission at ~96 au. These structural features are most clearly visible on the southern side of the disk, which is inclined toward the observer and exhibits enhanced SB due to the forward scattering of stellar light by dust grains. In the Q_ϕ image (Fig. 9g), we highlight three such features, labeled “F1”, “F2”, and “F3”, all of which have been detected in at least three separate epochs of IRDIS polarimetric observations, albeit with varying S/N (see also Appendix B and Fig. B.1). The most pronounced of these, “F2”, was also identified in HST/STIS and HST/NICMOS data (Ren et al. 2021).

These arc-like features may share a similar origin with the fast-moving clumps observed in the edge-on disk of AU Mic (Boccaletti et al. 2018). In both systems, sub-micron dust grains may be expelled by strong stellar winds from their active M-type host stars (Strubbe & Chiang 2006; Schüppler et al. 2015), potentially triggered by collisions in a secondary belt or in the vicinity of a planetary companion (e.g., Chiang & Fung 2017; Sezestre et al. 2017). Notably, a candidate Saturn-mass planet located at a projected distance of ~52 au or 1.5″, coincident with the position of TWA 7’s second planetesimal ring, has recently been detected with JWST/MIRI (Lagrange et al. 2025). This ring is both very narrow and flanked by two gaps, appearing underluminous at the planet’s location relative to other azimuths (Fig. 9g). Such a morphology supports the scenario of a resonant planetesimal ring sculpted by the planet, which may be carving the adjacent gaps and generating a local void.

If confirmed, this planetary companion could be responsible for gravitational perturbations that locally enhance dust production. Once released, small grains are redistributed by interactions with stellar wind and radiation pressure, giving rise to asymmetric structures such as arcs, streamers, or clumps, depending on the disk inclination and viewing geometry. The nearly pole-on orientation of the TWA 7 disk may thus offer a complementary view of the dynamic processes that shape AU Mic’s edge-on disk.

BD-20951. The highly inclined circumbinary debris disk around the SB2 BD-20951 (Torres et al. 2008) is resolved for the first time in both total and polarized scattered light with SPHERE in the H band (Perrot et al., in prep). The primary is a K1V(e) star, and the binary components have an estimated flux ratio of ~0.25 (Elliott et al. 2014). The system may be a member of the Carina MG (28±11 Myr; Gratton et al. 2024) or TucanaHorologium association (37±11 Myr; Gratton et al. 2024) as proposed by Torres et al. (2008), although the BANYAN Σ tool classifies it as a field star.

The IR excess was identified by Moór et al. (2016) who noted that the colder component may significantly contribute to the total near-IR flux of the system. The residuals in the PCA-reduced image suggest that the disk possesses sweap-back wings (Fig. 2). The geometrical parameters of belt are specified in Table 2.

4.5 Modeling of selected planetesimal belts

To investigate the morphology of the detected debris disks and potential correlations between disk parameters and host star properties, we fitted images of several debris belts (listed in Table G.1) using a grid of models that simulate single scattering of stellar photons by dust particles in an optically thin disk. A key advantage of our approach, compared to studies focused on individual disks, is the use of a uniform modeling framework for all systems in our sample. This consistency allows for a more direct and meaningful comparison of the derived disk parameters.

4.5.1 Model for scattered light

To estimate the fundamental geometric parameters of the debris belts, we generated synthetic images of scattered (or polarized) light and compared them with the observed disk images. To create these synthetic images, we employed a 3D, rotationally symmetric model to describe the spatial distribution of grain number density, n_gr(r, h), within the disk. Following the approach of Augereau et al. (1999), we characterize this distribution as the product of a radial profile R(r) and a Gaussian function Z(r, h). The profile R(r) defines the variation of grain number density in the disk midplane as a function of radial distance from the star r. Meanwhile, the Gaussian function determines the vertical profile of n_gr(r, h), shaping its distribution in the direction perpendicular to the disk midplane, as described by the height coordinate h: $$\begin{array}{rl}& n_{\mathrm{g}\mathrm{r}}(r,h)\sim R(r)\times Z(r,h)\\ & ={({\left(\frac{r}{r_0}\right)}^{-2{\alpha }_{in}}+{\left(\frac{r}{r_0}\right)}^{-2{\alpha }_{out}})}^{-1/2}\times \exp [-\ln 2{\left(\frac{|h|}{H(r)}\right)}^2],\end{array}$$

where r₀ is the reference radius of the debris belt, α_in > 0 and α_out < 0 are the exponents of the radial power laws for the dust distributions inside and outside of the belt, respectively. The scale height of disk the H(r) is defined as a half-width at halfmaximum (HWHM) of the Gaussian profile at radial distance r. In this work, the scale height is assumed to increase linearly with radial distance, following the relation H(r) = H₀(r/r₀)^β, where H₀ = H(r₀) and the disk flaring index β is fixed at β = 1.

In the model, each dust grain location within the disk is associated with a scattering angle θ, defined as the angle between the incident stellar ray striking a dust particle and the observer’s line of sight, where θ = 0 corresponds to forward scattering. The scattering angle is a crucial model parameter, as it governs the fraction of incident stellar light scattered in a particular direction, described by the so-called scattering phase function (SPF). To generate synthetic disk images, we employed the Henyey–Greenstein (HG) function as the SPF (Henyey & Greenstein 1941), which provides a convenient parametrization of anisotropic scattering by dust grains: $$SPF(\theta ,g)=\frac{1-g^2}{4\pi {(1+g^2-2g\cos \theta )}^{3/2}},$$

where g represents the HG scattering asymmetry parameter. This parameter quantifies the preferential direction of scattering, ranging from g = −1 (backward scattering) to g = 1 (forward scattering), with g = 0 corresponding to isotropic scattering.

Based on observations of scattered light from zodiacal and cometary dust in the Solar System (e.g., Leinert et al. 1976; Bertini et al. 2017), as well as laboratory experiments with dust analogs (e.g., Frattin et al. 2019; Muñoz et al. 2017), interplanetary dust grains are expected to preferentially scatter radiation in the forward direction. As a result, in disk images, the side of the disk that is closer to the observer appears brighter.

The forward-scattering behavior of dust particles can be described using a HG function with a positive asymmetry parameter (g > 0). However, real SPFs derived from observational data often exhibit a more complex structure: a pronounced diffraction peak at small scattering angles (g >> 0), a relatively flat mid-range (g ~ 0), and an enhanced backscattering component (g < 0). Consequently, a more accurate representation of an actual SPF would require a combination of three HG functions.

Nevertheless, we opted to model the data using a single HG function. This choice is justified by the fact that, in most cases, the inclination of resolved debris disks does not permit the measurement of scattering intensity across the full range of scattering angles (0° to 180°). Instead, we can only fit the portion of the SPF that is accessible in the data, which can be adequately approximated by a single HG function.

Another simplification adopted in our modeling approach is the assumption that the optical characteristics of dust grains, defined by their composition, shape, and size, are spatially uniform across the disk and can be represented by a single SPF. While this assumption simplifies the modeling process, it remains a coarse approximation, as the SPF is inherently dependent on dust properties that are expected to vary with radial distance. For instance, at the radial location of the peak grain number density, R_belt, debris spanning a wide range of sizes is typically present, from submicron grains to kilometer-sized planetesimals (e.g., Wyatt 2008). In contrast, the outer disk is expected to form a halo of small grains that are collisionally produced within the main ring and subsequently placed on high-eccentricity orbits by stellar radiation pressure. This halo is anticipated to exhibit strong size segregation (Thebault et al. 2014), with the dominant grain size decreasing with increasing radial distance.

Depending on the quality of the available data, we fitted either an image of total intensity or polarized intensity. To generate a polarized intensity image (Q_ϕ image), we employed a polarized scattering phase function (pSPF) given by a functional form (e.g., Engler et al. 2017) $$pSPF(\theta ,g)=p_{max}\frac{1-{\cos }^2\theta }{1+{\cos }^2\theta }\frac{1-g^2}{4\pi {(1+g^2-2g\cos \theta )}^{3/2}},$$(2)

where p_max is the maximum polarization fraction of dust particles.

In this equation, we employed a polarization fraction function p(θ) characteristic of Rayleigh scattering, in which the maximum polarization fraction is attained at a scattering angle of θ = 90° (Bohren & Huffman 1983): $$p(\theta )=\frac{pSPF(\theta )}{SPF(\theta )}=p_{max}\frac{1-{\cos }^2\theta }{1+{\cos }^2\theta }.$$(3)

Rayleigh scattering describes the scattering of light by particles whose sizes are at least an order of magnitude smaller than the wavelength of the incident radiation. In near-IR observations of debris disks, this condition is generally satisfied in the outer disk regions, where the abundance of larger grains declines and the contribution of small particles to the scattered light becomes increasingly significant. The polarization fraction function of this small-particle population as well as that of micron-sized grains, typically traced in near-IR scattered-light images, can be reasonably approximated by the Rayleigh scattering function p(θ), as defined in Eq. (3) (see Appendix C).

For each disk specified in Table G.1, we generated a large set of models using a grid of fitting parameters. The parameter ranges were individually defined based on the findings of previous studies on these targets. We deliberately chose to use a model grid approach rather than a Markov Chain Monte Carlo (MCMC) algorithm, which, although widely employed in disk modeling, often yields unrealistically small uncertainties on the best-fitting parameters. Mazoyer et al. (2020) showed indeed that MCMC often leads to underestimated uncertainties, an effect probably due to the non-Gaussian statistics of the residual noise in coronagraphic images (Pairet et al. 2019).

To identify a family of models that adequately reproduce the observed disk images, we therefore adopted a more conservative threshold for the χ² value than the one derived from the probability distribution of parameters, which assumes normally distributed errors: $${\chi }^2<{\chi }_{min}^2+\mathrm{\Delta }{\chi }^2,$$

where ${\chi }_{\text{min}}^2$ is the minimum chi-square value obtained from the fits, and $\mathrm{\Delta }{\chi }^2=\sqrt{2\nu }$ with ν denoting the number of degrees of freedom⁸ (Thalmann et al. 2013). The mean values and standard deviations of the parameter distributions from the family of well-fitting models are adopted as the best-fitting parameters and their corresponding uncertainties, as reported in Table G.1. Note that these estimates can be affected by the parameter step size, as it determines the sample size and can influence the precision of the evaluated sample mean and standard deviation⁹. Also, the best-fitting parameters do not correspond to the single model with the minimum χ², but rather reflect the statistical properties of the family of acceptable models.

HD 39060. For the HD 39060 (β Pic) debris disk, the model consists of outer and inner planetesimal belts and has, therefore, a double number of parameters. We chose a double-belt model for this particular target because IRDIS images in both total and polarized intensity reveal an inner disk located within the main outer belt (Fig. 10). While the inner belt shares a similar inclination with the outer belt, it has a slightly different PA. Consequently, its extensions become visible, producing the characteristic “butterfly pattern” in the scattered light distribution observed in HD 39060 disk images (e.g., Golimowski et al. 2006; Ahmic et al. 2009).

In the polarized intensity image of HD 39060 (Q_ϕ image), the inner disk becomes more distinct when the polarized flux from the outer belt is removed. To achieve this, we rotated the Q_ϕ image (Fig. 10b) counterclockwise by 60°, aligning the major axis of the outer belt horizontally. We then subtracted the left half of the Q_ϕ image from the right half and vice versa. The resulting image, shown in Fig. 10b, reveals the near side of the inner disk, which becomes visible in the lower left and upper right quadrants. In this image, the polarized flux from the outer belt is largely eliminated due to its symmetrical distribution with respect to the vertical axis. The polarized flux from the inner belt is partially reduced, particularly near the image center and in the upper left quadrant, due to the asymmetrical distribution of flux from the inner belt relative to the vertical axis of the image.

HD 39060 is the only target for which we applied a double-belt model. Other targets, such as HD 15115 (Engler et al. 2019; MacGregor et al. 2019), HD 120326 and HD 129590, also exhibit indications of a second planetesimal belt, though it is only marginally resolved (Sect. 4.4). The quality and spatial resolution of the available data do not allow for a reliable fit using a double-belt model to obtain robust constraints on the parameters of the secondary component. We modeled these systems using a single-belt approach despite their multiple-belt structure. In such cases, the fitted parameter values may be influenced by the presence of the second component, particularly affecting the derived belt radius or scale height, as discussed in the next section.

4.5.2 Discussion of modeling results

Comparison between measured and modeled radial distances of the planetesimal belts

The reference radius r₀, in combination with the model parameters α_in and α_out obtained from disk image modeling, determines the modeled radial position of the debris belt $R_{\text{belt}}^{mod}$: $$R_{\mathrm{b}\mathrm{e}\mathrm{l}\mathrm{t}}^{mod}=r_0{(-\frac{{\alpha }_{in}}{{\alpha }_{out}})}^{1/(2{\alpha }_{in}-2{\alpha }_{out})}.$$(4)

The modeled radius $R_{\text{belt}}^{\text{mod}}$ defines the location of the peak grain volume density in the radial profile of the disk midplane and determines the region where collisions between larger debris fragments or planetesimals generate dust particles. In the ADI-processed scattered-light images, the radial position of the SB peak measured along the disk’s major axis may slightly deviate from the actual location of the planetesimal belt. This discrepancy arises from a combination of factors, including stellar illumination, spatial resolution of instruments, geometrical projection effects and the asymmetry of the SPF. For pole-on disks, for instance, the observed SB peak of the radial profile in the r²−scaled images is expected to be at a radial position of the modeled peak surface density of grains, $R_{max(\sigma )}^{mod}$, which can be evaluated through the integration over the whole disk height in the vertical direction (Augereau et al. 1999): $$R_{max(\sigma )}^{mod}=r_0{(-\frac{{\alpha }_{in}+\beta }{{\alpha }_{out}+\beta })}^{1/(2{\alpha }_{in}-2{\alpha }_{out})},$$

where β = 1 is adopted in our model (see Sect. 4.5.1). With this value, the difference between the derived $R_{\text{belt}}^{mod}$ and $R_{max(\sigma )}^{mod}$ remains within 4% for all modeled debris disks listed in Table G.1 except for the HD 61005 disk, where the discrepancy reaches 10%.

In Fig. 11a, we compare the belt radii derived from disk image modeling (Col. 2 of Table G.1) with the radial locations of the SB peaks measured directly in the r²−scaled images (Col. 2 of Table 2). The black line in Fig. 11a represents an empirical fit and coincides with the 1:1 relation, indicating that the modeled belt radii are in good agreement with the directly measured values. For most targets, the difference between $R_{\text{belt}}^{\text{mes}}$ and $R_{\text{belt}}^{\text{mod}}$ is within 10%, and within 12% between $R_{\text{belt}}^{\text{mes}}$ and $R_{max(\sigma )}^{mod}$. There are two noticeable exceptions: HD 120326 and HD 129590. Both of these targets are likely to host at least two distinct planetesimal belts (Sect. 4.4), which may explain the observed deviations.

Asymmetry parameter

The derived values for the HG asymmetry parameter g range from 0.79 (for HD 106906 disk) to 0.5 (for HD 172555 disk). We observe a slight trend toward higher asymmetry parameters when modeling disks with higher inclinations (Fig. 11b). This trend can be attributed to the broader range of scattering angles accessible in highly inclined disks, which allows for a more pronounced forward-scattering peak to be observed. Consequently, this suggests that low-inclination disks may intrinsically exhibit higher asymmetry parameters, but their forward-scattering component remains less apparent due to the limited range of observable scattering angles. If these disks were viewed at higher inclinations, their asymmetry parameters might appear larger.

Exponents of radial power law for the radial distribution of grain number density

The best-fit values for the exponent of the radial power law α_in span a relatively wide range, from 2.3 for the HD 115600 disk to 25.9 for the HD 114082 disk. However, it is important to note that three of the four highest exponents (>20), shown as open circles in Fig. 11c, are not well constrained. Their distributions lack a clear peak indicating an optimal fit within the tested parameter space. Instead, these values consistently trend toward the upper limit of the parameter range, even when the maximum tested value is set as high as α_in = 80. Notably, these unconstrained values were obtained from the fitting of total intensity images of the HD 114082, HD 117214 and HD 131488 disks, which may suggest a limitation in accurately determining this parameter using ADI forward-modeling, particularly for disks with high inclinations and small angular sizes.

Alternatively, these disks may indeed possess extremely sharp inner edges, a feature often interpreted as evidence of unseen planets clearing the space at the edges of the planetesimal belts. Milli et al. (2017b) tested methods to constrain large values of the α_in and α_out parameters and concluded that the modeling of these parameters is limited by the intrinsic steepness of the PSF. Specifically, belt edges that are steeper than the PSF wings are inherently blurred by convolution with the PSF, making it impossible to constrain α_in and α_out values steeper than approximately 30 in the IRDIS H band.

If the aforementioned targets are excluded, a trend emerges in which the α_in parameter increases with stellar luminosity (Fig. 11c). Regarding α_out (Fig. 11f), the derived values for most systems are consistent with the archival data presented in Table 1 of Thebault et al. (2023). As discussed in that study, the expected α_out value for a typical belt-like system with an outer halo composed of small grains placed on high-eccentricity orbits by radiation pressure is approx. −2.5. When considering error bars, we find that roughly half of the systems in Table G.1 have α_out estimates that are compatible with this reference value. For the remaining systems, we obtain steeper outer profiles, with α_out reaching values as low as −8 or even −12. In these cases, additional dynamical processes, such as perturbations by external planets or stellar companions, are likely responsible for clearing out the outer disk regions (Thebault et al. 2023).

Vertical disk structure

As discussed in Sect. 4.5.1, we adopted a Gaussian function to describe the vertical profile of each disk, accounting for its nonzero vertical width. In this study, we define the HWHM of the Gaussian profile as the scale height H(r) of the disk at a given radial distance r. The ratio of the scale height to the radial distance, A_disk = H(r)/r, is referred to as the disk aspect ratio, also known as the disk half-opening angle. Since, in our model, the scale height H(r) varies linearly with distance r, the aspect ratio remains constant throughout the disk and can be expressed as A_disk = H₀/r₀, where H₀ is the scale height at the reference radius r₀.

The scale height and aspect ratio of a disk serve as key indicators of the dynamical excitation within a debris system. Thebault (2009) numerically estimated a minimum aspect ratio of 0.04 ± 0.02 for a collisionally evolving disk observed at visible to mid-IR wavelengths. If the disk experiences additional dynamical perturbations from massive bodies such as giant exoplanets, its aspect ratio is expected to exceed 0.06.

In Figures 11d and 11e, we show the results of our analysis of the aspect ratios obtained for the modeled disks. We consider the disk to be resolved in direction perpendicular to the disk midplane, if the fitted FWHM of the vertical profile is larger than the FWHM of the stellar PSF. For four modeled disks (HD 114082, HD 117214, HD 120326, HD 131488) which are shown in Fig. 11d with open markers this condition is not fulfilled. We note that the parameters of these four disks are derived from the fitting of the total intensity images applying ADI forward-modeling approach. Therefore, the value of the scale height might be underestimated. The other two disks (HD 109573, HD 181327) are only marginally resolved in vertical direction. Both disks have an inclination lower than 80°. In particular the HD 181327 disk has an inclination of ~30° which is the lowest one we modeled. This might reflect the challenge of modeling the disk scale height when using images obtained with the ADI technique or images of low inclined disks.

Most of the modeled disks have a scale height between 0.02 and 0.06 (Fig. 11d). This relatively small vertical extent can be explained by the combined effects of radiation pressure acting on small dust particles and their mutual collisions, which naturally regulate the disk’s thickness (Thebault 2009). In contrast, the HD 172555 disk exhibits a significant larger scale height of 0.1, which may indicate the influence of additional massive perturbers, although this result could be affected by the lower S/N of the image. This system is particularly intriguing, as it contains detected gas and a notable abundance of hot, small dust grains, features that may be the aftermath of a recent, violent collision between planetary bodies (Lisse et al. 2017; Riviere-Marichalar et al. 2012; Kiefer et al. 2014; Engler et al. 2018; Schneiderman et al. 2021; Samland et al. 2025).

Three targets (HD 106906, HD 115600, and HD 129590) with modeled scale heights ranging between 0.06 and 0.07 are enclosed by the red ellipse in Fig. 11d. All three disks are highly inclined and are suspected to host at least two cold belts. If so, one possible explanation for the relatively large scale heights inferred from the models is that the inner belts remain unresolved but lie in close projected proximity to the outer belts along the minor axis in the scattered-light images. This configuration could mimic the appearance of a geometrically thicker planetesimal belt. An alternative explanation, at least for HD 106906, is dynamical excitation of the disk by a massive substellar companion known to be present in the system (see Sect. 7). It is also worth noting that the HD 129590 disk contains small amounts of CO gas (M_CO = 10⁻⁵−10⁻⁴M_⊕; Kral et al. 2020). In contrast, no gas has been detected so far in HD 106906 and HD 115600, leaving the influence of gas on the observed disk scale heights in these systems uncertain.

However, there are four gas-rich systems with the estimated CO masses exceeding 10⁻² M_⊕, namely HD 9672, HD 32297, HD 121617 (Moór et al. 2019), and HD 131488 (Moór et al. 2017; Pawellek et al. 2024). These systems are marked by triangles in Figs. 11d and 11e, and their aspect ratios span a range between 0.026 and 0.011. In these disks, the CO emission has been observed to be axisymmetric and co-located with the millimeter-sized dust particles (Hughes et al. 2017; Moór et al. 2017; MacGregor et al. 2018). The relatively low aspect ratios derived for these gas-rich systems may suggest that micron-sized dust grains are dynamically coupled to the gas, leading to their settling toward the disk midplane.

There appears to be a trend of decreasing disk aspect ratio with increasing stellar luminosity, as observed in Fig. 11e. This trend could also account for the low aspect ratios of gas-rich disks, as three of these systems are associated with high-luminosity A-type stars (L_⋆ > 13 L_⊙), and, especially as we derive an aspect ratio of A_disk = 0.017 for the debris belt around the AOV star HD 109573 (L_⋆ = 25.2 L_⊙), which is not known to be gas-rich, further supporting this possible correlation. However, confirming this trend would require a significantly larger sample of measured aspect ratios.

5 SED modeling

We applied SED modeling to characterize the thermal emission of debris disks in our sample, using two different approaches to fit the photometric data. The modified BB (MBB) approach (Backman & Paresce 1993) provides a uniform fitting method for all targets, making it particularly suitable for a statistical analysis of the sample. In contrast, the particle size distribution (SD) approach (e.g., Müller et al. 2010; Pawellek et al. 2021) allows for a more detailed characterization of dust grain properties, including the dominant grain size, the SD steepness, and the bulk optical properties of the dust. The SD model is only applicable to spatially resolved debris disks, as determining the disk radius is necessary to break the degeneracy between the location of dust particles and their sizes (Pawellek et al. 2014).

5.1 Modeling procedure

We utilized photometric data for our sample from published catalogs, such as 2MASS (Cutri et al. 2003), the WISE All-Sky Release Catalog (Wright et al. 2010), the AKARI All-Sky Catalog (Ishihara et al. 2010), the Spitzer Heritage Archive (Carpenter et al. 2008; Lebouteiller et al. 2011; Chen et al. 2014; Sierchio et al. 2014), and the Herschel Point Source Catalog (Marton et al. 2015; Marshall et al. 2021). In addition, we used data published in the literature (e.g., Chen et al. 2014; Matrà et al. 2017; Marshall et al. 2021). These data allow us to analyze the SEDs and assess the presence of IR excess emission beyond what is expected from the stellar photosphere.

To find excess emission, we fitted an SED model consisting of a star and a disk. Firstly, we fitted PHOENIX stellar photosphere models (Brott & Hauschildt 2005) for each target using the stellar luminosity and the stellar temperature as model parameters, and photometric data in the VIS/NIR where the stellar emission is supposed to dominate the SED and the disk emission is negligible. The resulting stellar luminosities and temperatures are listed in Table G.2. Secondly, knowing the stellar contribution to the mid- and far-IR data, we derived the excess emission in the appropriate wavelength bands between ~5 and ~1000 μm taking into account the uncertainties of the photometry and the photospheric model.

We followed the four criteria given in Ballering et al. (2013) and Pawellek et al. (2014) to check for the presence of a warm disk component in addition to the cold Kuiper belt analog. Firstly, the number of photometric data points must be large enough so that the data are not over-fitted. If that was the case in a second step we considered a warm component to be present, if there is a significant excess (≥3σ) in either the WISE/22 or MIPS/24 in excess of that which could originate in a single ring fitted to longer wavelength data. Thirdly, the fit of the two-component SED has to be much better than the one-component fit according to the Bayesian information criterion (BIC) $$\mathrm{B}\mathrm{I}\mathrm{C}={\chi }^2+J{\log }_e\left(N_{\mathrm{d}\mathrm{a}\mathrm{t}\mathrm{a}}\right),$$(5)

where J represents the number of free parameters and N_data the number of data points. We used the classification given in Kass & Raftery (1995) to infer whether a one- or a two-component model is more likely (Pawellek et al. 2021). As a fourth criterion, we required the inferred ring containing the warm dust to be located outside the sublimation radius (assuming 1300 K as the sublimation temperature for astrosilicate). If all four criteria were met, we assumed the SED to consist of a two-component model.

The uncertainties of the fit parameters were inferred in the following way. We started at the position of the minimum χ² in parameter space, e.g., from the best-fitting fractional luminosity, f_disk, and BB radius, R_MBB, in case of an MBB model (Sect. 5.2). A set of new parameter values was randomly generated from which we calculated the SED. This led to a new χ² value that was compared to the previous minimum value. The χ² parameter estimates how likely the set of parameter values fits the SED. If the probability was larger than a certain threshold value, the set was saved. In the end, it was counted how often the code reached a certain set of f_disk and R_MBB. The closer the parameters get to the best-fitting values, the higher the probability. The resulting distribution in parameter space represents an estimate for the probability distribution of the parameters and thus, allows us to calculate the confidence levels for the parameters assuming that the values follow a normal distribution in parameter space (simulated annealing; e.g., Pawellek 2017).

5.2 Modified blackbody model

Every disk in the sample was fitted with a MBB model for which the thermal emission of the dust is described as $$F_{\lambda }^{th}\sim B_{\lambda }(\lambda ,T_{\mathrm{M}\mathrm{B}\mathrm{B}})[u({\lambda }_0-\lambda )+u(\lambda -{\lambda }_0){\left(\frac{\lambda }{{\lambda }_0}\right)}^{-{\beta }_{\mathrm{o}\mathrm{p}}}],$$(6)

where $F_{\lambda }^{\mathrm{th}}$ is the spectral flux density of thermal emission, B_λ is the Planck function and u the Heaviside step function. The parameter λ₀ represents the characteristic wavelength, while β_op is the spectral opacity index.

From this model we derived the dust temperature, T_MBB, and the resulting BB radius of the disk, R_MBB, as well as the fractional luminosity, f_disk. Here, R_MBB is the distance from the star that the temperature implies if the dust acted like BB in equilibrium with the stellar radiation. If a warm component was present (Sect. 5.1), we modeled it also with a MBB model, but assumed that the β_op parameter of the warm and cold component are similar to keep the number of free parameters as low as possible, and to avoid degeneracies between the component parameters.

For a single ring model the free modeling parameters are the fractional luminosity f_disk, the BB radius R_MBB, the characteristic wavelength λ₀, and the opacity index β_op. Hence, at least five data points are needed to not over-fit the photometric data. In case of a two-component model we added the BB temperature, the fractional luminosity, and the characteristic wavelength of the inner ring as free parameters which required at least eight data points.

We considered the SED model to be unreliable, i.e. it cannot be fitted to the data, if one of the following conditions is not fulfilled: (1) The number of available photometric data points is higher than a minimum number of data points needed to fit an SED, the number of needed data points varies with the modeling approach (MBB or SD); (2) There is a clear detection of the IR excess emission compared to the stellar photosphere, i.e. the total flux density exceeds the stellar photosphere by at least 3σ; (3) The photometric data cover the peak of the thermal emission to constrain the model.

Following these criteria, the SEDs are counted as unreliable in case for the following targets: GSC 7396-0759, HIP 63942, HD 35114, HD 36968, HD 53842, HD 69830, HD 122705, HD 135379, HD 141011, HD 141943, HD 181869, and HD 274255.

Two targets (HD 17390 and TWA 25) fulfill points (2) and (3), but not point (1), as they have only three data points. For these targets we assumed a pure BB model where we only have two free parameters and thus, only needed three data points to achieve a proper fit.

5.3 Results of MBB modeling

The results of the one-component MBB fitting are specified in Table G.3. Based on the criteria outlined in Sect. 5.1, fourteen SEDs were fitted using a two-component model, incorporating both warm and cold dust components. The best-fit parameters from this modeling are provided in Table G.4. Figure E.1 (top row) shows examples of SEDs fitted with MBB models consisting of one or two components. Using the results of these fits, we investigated the evolution of the disk IR excess, the correlation between the estimated dust mass in the belt and its radial distance from the star for A-type and F-type stars in our sample, as well as the correlation between dust mass and host star mass for all targets.

The disk IR excess is one of the key factors influencing disk detectability (Sect. 6). Younger debris disks exhibit higher fractional luminosities, as their luminosity scales with dust mass, which is at its peak in early evolutionary stages. Various dust removal mechanisms, including the blowout of the smallest dust particles due to radiation pressure, the Poynting-Robertson drag, and photoevaporation, play a crucial role in shaping the evolution of disk mass and luminosity (e.g., Dominik & Decin 2003). Previous studies (e.g., Kalas et al. 2000; Rieke et al. 2005) have shown that the fractional luminosities of debris disks, f_disk, and consequently their dust content, follow a power-law dependence on time $$f_{\text{disk\hspace{0.17em}}}\propto t_{\text{age\hspace{0.17em}}}^{{\alpha }_t}.$$

Figure 12 shows the evolution of disk fractional luminosity for the debris disks around A- and F-type stars in our sample. By fitting the data points in the [f_disk, t_age] parameter space, we derived the power-law indices α_t = −1.18 ± 0.14 for A-type stars and α_t = −0.81 ± 0.12 for F-type stars. These exponents are in general agreement, though formally distinct within 1σ uncertainties with the steady-state collisional evolution theory of planetesimal belts, which predicts a disk luminosity decline proportional to $t_{\text{age}}^{-1}$ (Dominik & Decin 2003). We note that the α_t = −1 slope (black dash-dotted line in Fig. 12) is only valid for collisional systems in steady-state, where the disk age exceeds the collisional lifetimes of the largest planetesimals it contains. As pointed out by Löhne et al. (2008), this condition might be met only in very old systems (>1 Gyr), or in very massive disks and/or disks at short radial distances from their host stars. Given that the median age of our sample is ~100 Myr, and the median belt distances are around 70–80 au, our results do not align with these expectations.

Our results show however that debris disks around A-type stars tend to decline more rapidly in brightness than those around F-type stars. This behavior can be interpreted in the context of differences in initial disk masses and their influence on collisional evolution timescales. Debris disks with higher masses produce dust more efficiently due to the increased frequency of planetesimal collisions (e.g., Wyatt 2008). However, these systems also deplete more rapidly, since both large planetesimals and the dust generated by their collisions are removed more quickly, either through further collisional grinding or via radiative forces acting on small grains (Löhne et al. 2008; Krivov 2010). A-type stars, being more massive and luminous than F-type stars, are expected to host more massive planetesimal belts, which evolve faster due to the enhanced dynamical excitation and stronger radiation pressure that efficiently clears small grains from the system. In contrast, the more gradual decay of IR excess observed in F-type systems suggests slower collisional processing, consistent with initially less massive disks and weaker dynamical stirring.

To further investigate this trend and quantify the effect of stellar mass on debris disk evolution, we estimated dust masses for our targets using the MBB best-fit parameters and equation for the disk dust mass following Wyatt (2008): $$M_{\mathrm{d}\mathrm{u}\mathrm{s}\mathrm{t}}=12.6f_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{k}}{R}_{\mathrm{M}\mathrm{B}\mathrm{B}}^2{k}_{850}^{-1}{\left(\frac{850\mu \mathrm{m}}{{\lambda }_0}\right)}^{-{\beta }_{\mathrm{o}\mathrm{p}}}\text{,\hspace{0.17em}}$$(7)

where k₈₅₀ = 45 au²M_⊕⁻¹, and λ₀ and β_op are the characteristic wavelength and spectral opacity index, respectively, obtained from the MBB fit (see Eq. (6)).

The derived dust masses of the planetesimal belts, plotted as a function of stellar mass, are shown in Fig. 13. Since the majority of the debris disks detected in our sample are younger than 50 Myr (90% of all detections), we divided the sample into two groups based on stellar age: systems aged 10–50 Myr and older systems. These are represented by orange dots and red diamonds in Fig. 13, respectively. The detected debris disks are indicated by filled markers.

To investigate the M_dust − M_⋆ relation and its potential evolution over time, we fitted a power-law function in log-log space log(M_dust/M_⊕) = α_mass log (M_⋆/M_⊙) + β_mass to three datasets: disks in the 10–50 Myr range, disks older than 50 Myr, and all disks older than 10 Myr. The best-fit parameters α_mass and β_mass for the three subsamples are reported in Table 3 (rows 3–5) and the fits are visualized in Fig. 13, along with 68% and 95% confidence intervals for the full sample. In all cases, we find a steeper than a linear relation (α_mass > 1) between the estimated disk dust mass and the mass of the host star. Notably, the scaling relation for the younger group of debris disks (see row 3 of Table 3) is steeper than that of the older group (row 4 of Table 3). This decrease in steepness with age may indicate a more rapid dust mass depletion in initially more massive disks, a trend that is also apparent in Fig. 12.

It is well established that PPD masses correlate with stellar mass as well (e.g., Andrews et al. 2013; Barenfeld et al. 2016; Ansdell et al. 2016). For example, Pascucci et al. (2016) analyzed the scaling of disk masses with stellar mass in star-forming regions, using various assumptions for the dust temperature–luminosity relation, and similarly found super-linear trends. To place our findings in context, we included in Table 3 (rows 1–2) the best-fit parameters reported by Pascucci et al. (2016) for the Chamaeleon I star-forming region and overplot the corresponding fits in Fig. 13 as black solid lines. As evident from this figure, the PPD relation is very similar to that we found for the debris disks and supports the idea that debris disks reflect the initial conditions set during the protoplanetary phase. This is particularly significant, as PPD masses represent the material reservoir available for planet formation. A super-linear scaling of disk mass with stellar mass therefore implies that the resulting planet populations, both in terms of typical masses and occurrence rates, are also expected to have a positive correlation with stellar mass.

Comparing the offset of the linear relations (Col. 4 in Table 3) between PPD and debris disk dust masses (Fig. 13) we find that the average dust mass decreases by ~3 dex within the first 50 Myr, and by 3.7 dex for the older stars in our sample. However, the dust masses exhibit a relatively large spread, which can be partially attributed to the variation in belt radii among the systems with the same mass of the host star, and scaling the disk mass with its radius.

A correlation between dust mass and radial location is consistent with previous studies suggesting that belts at larger distances from the star are likely to be more massive (e.g., Andrews et al. 2013; Matrà et al. 2025). This trend can be explained by the fact that a wider belt spans a larger volume, thereby allowing for a greater population of dust-producing planetesimals, assuming a roughly constant surface density or collisional activity per unit area. Moreover, the collisional timescales in outer regions of debris disks are longer due to lower orbital velocities and decreased dynamical stirring, allowing dust to persist for extended periods and accumulate to a higher levels (e.g., Wyatt 2008).

In Fig. 14, we show the derived belt dust masses as a function of the BB belt radius R_MBB for A-type stars (panel a) and F-type stars (panel b). The contour plots in this figure represent the bivariate probability density function (PDF) for subsamples of stars with estimated ages between 10 and 50 Myr (red contours) and stars older than 50 Myr (violet contours). A comparison of the PDF peak positions across different stellar age bins indicates that the dust mass declines by 1–1.5 dex on average for older stars, particularly among A-type stars, in line with our results presented in Fig. 13.

We also observe a trend of increasing belt dust mass with growing radial distance from the star across all age bins. This trend follows the fitted radial power law $M_{\text{dust\hspace{0.17em}}}\propto R_{\mathrm{M}\mathrm{B}\mathrm{B}}^{{\alpha }_{\mathrm{R}}}$, with α_R = 2.3 ± 0.4 for A-type stars in the 10–50 Myr age bin and α_R = 1.9 ± 0.4 for A-type stars older than 50 Myr. Similarly, for F-type stars, we find α_R = 2.6 ± 0.4 for the 10–50 Myr bin and α_R = 2.1 ± 0.5 for older stars. This result is expected, as it is consistent with the power-law form used in Eq. (7), within the given uncertainties. Therefore, we further investigate this relationship using belt dust masses derived from SD modeling (Sect. 5.5).

5.4 Size distribution model

In case of spatially resolved disks, we used the SONATA code (Pawellek et al. 2014; Pawellek & Krivov 2015) to model the SEDs with a dust SD. While for the MBB model we simply fitted a dust temperature and a fractional luminosity without consideration of dust properties, the SONATA code calculates the temperature and the thermal emission of dust particles at different distances to the star following $$r\left(T_{\text{grain\hspace{0.17em}}}\right)=\frac{R_{\star }}{2}{\left[\frac{{\int }_0^{\mathrm{\infty }}{Q}_{\lambda }^{\mathrm{a}\mathrm{b}\mathrm{s}}(a)B_{\lambda }\left(T_{\star }\right)\mathrm{d}\lambda }{{\int }_0^{\mathrm{\infty }}{Q}_{\lambda }^{\mathrm{a}\mathrm{b}\mathrm{s}}(a)B_{\lambda }(T_{\text{grain}}(r))\mathrm{d}\lambda }\right]}^{1/2},$$(8)

where B_λ is the Planck function, T_⋆ and R_⋆ the stellar temperature and radius, and T_grain the grain temperature. The parameter $Q_{\lambda }^{\text{abs}}(a)$ gives the absorption efficiency dependent on wavelength λ and grain radius a. As seen in Eq. (8), this equation must be solved iteratively to determine the particle temperature as a function of distance from the star.

We assumed compact spherical grains and used Mie theory to compute the absorption efficiencies (Bohren & Huffman 1983). The dust composition was set to pure astronomical silicate (Draine 2003a), with a bulk density of ϱ = 3.3 g/cm³. The SONATA code integrates the emission from particles over a range of sizes to generate the SED. As mentioned before, flux densities at wavelengths shorter than 5 μm were excluded from the dust disk fitting, as the stellar photosphere dominates the emission in this regime.

We applied a power law for the SD of the dust and assumed a Gaussian radial distribution for the spatially resolved ring using the surface number density N_SED(r, a) similar to Pawellek et al. (2021), $$N_{\mathrm{S}\mathrm{E}\mathrm{D}}(r,a)\sim a^{-q}\frac{1}{\sqrt{2\pi }\mathrm{\Delta }{R}_{\mathrm{b}\mathrm{e}\mathrm{l}\mathrm{t}}}\exp [-\frac{1}{2}{\left(\frac{r-R_{\mathrm{b}\mathrm{e}\mathrm{l}\mathrm{t}}^{\mathrm{m}\mathrm{e}\mathrm{s}}}{\mathrm{\Delta }{R}_{\mathrm{b}\mathrm{e}\mathrm{l}\mathrm{t}}}\right)}^2],$$

where r represents the distance to the star, $R_{\text{belt}}^{\text{mes}}$ the belt radius measured from the r²−scaled scattered-light images along the disk’s major axis, and ΔR_belt is taken to be $0.1R_{\text{belt}}^{\text{mes}}$, assuming a radial dust distribution width of approx. 20% of the belt radius. This assumption aligns reasonably well with the results from well-resolved debris belts exhibiting a wide range of eccentricities, such as HD 22049 (Booth et al. 2017), HD 109085 (Marino et al. 2017), HD 109573 (Milli et al. 2019), HD 181327 (Marino et al. 2016), HD 202628 (Faramaz et al. 2019), and HD 216956 (MacGregor et al. 2017; Kennedy 2020). The parameter a represents the grain radius, while q is the SD power-law index. The surface number density, N_SED(r, a), is directly related to the surface density, Σ(r, a), by the equation $$\mathrm{\Sigma }(r,a)da=\pi a^2{N}_{\mathrm{S}\mathrm{E}\mathrm{D}}(r,a)da.$$

We considered grain sizes ranging from a minimum value, a_min, to a fixed maximum value of a_max = 5000 μm. Grains larger than this were assumed to contribute negligibly to the SED, as particles absorb and emit radiation efficiently only at wavelengths shorter than their size. For instance, the efficiency of interaction with radiation drops for weakly absorbing materials at wavelengths longer than a/2π (Backman & Paresce 1993), while for moderately absorbing materials such as “dirty ice”, this critical wavelength is approximately equal to the particle size (Greenberg 1978). Therefore, the adopted maximum size of 5 mm is sufficiently large to encompass all grain sizes that significantly contribute to the SED flux within the wavelength range covered by the available observations.

A one-ring model has three free parameters: the minimum grain radius, a_min, the SD index, q, and the amount of dust, $M_{\text{dust}}^{\mathrm{S}\mathrm{D}}$, for particles between a_min and a_max assuming a bulk density ϱ. Hence, at least four photometric data points (and the disk radius) are needed to fit a SD model to the data.

If the scattered light images provided evidence for both an outer and an inner belt, and both rings were spatially resolved, we applied a two-component SD model to fit the data. To address the degeneracy in dust mass estimates between the two components, we determined the mass ratio as follows. First, we fitted the SED with a single-component model, considering a dust distribution spanning from the central radius of the inner ring to the central radius of the outer ring, allowing us to estimate the total dust mass of the entire disk. For simplicity, we assumed that the dust mass of a belt scales with the square of its radial distance from the star, such that it is given by $$M_{\text{outer\hspace{0.17em}}}=\frac{M_{\text{total\hspace{0.17em}}}\times R_{\text{outer\hspace{0.17em}}}^2}{R_{\text{inner\hspace{0.17em}}}^2+R_{\text{outer\hspace{0.17em}}}^2}{M}_{\text{inner\hspace{0.17em}}}=\frac{M_{\text{total\hspace{0.17em}}}\times R_{\text{inner\hspace{0.17em}}}^2}{R_{\text{inner\hspace{0.17em}}}^2+R_{\text{outer\hspace{0.17em}}}^2},$$

where M_inner and M_outer are the dust masses of the inner and outer belt, M_total the estimated total dust mass, and R_inner and R_outer the central radii of the inner and outer belt.

We then fixed the masses of the individual belts and fitted the SD parameters a_min and q, assuming that both cold dust rings share the same SD. In this approach, only the dominant grain size and the SD index are the only free parameters, as the dust masses are fixed. Consequently, at least three photometric data points are required to fit the SED.

With our chosen approach, we focused on spatially resolved data from scattered-light images. However, warm components that remain undetected in the images or are spatially unresolved are not included in the modeling process. As a result, the NIR/MIR data of such SEDs are not well-fitted, and part of the emission that could originate from a warm dust component is instead incorporated into the fit of the cold dust belt. In such cases, the fitted minimum grain size is underestimated.

To mitigate this issue, for targets with a sufficiently large number of available photometric points, we tested various fitting approaches. These included a fit with a warm component modeled as a BB, which provided an upper limit on a_min, and a fit excluding the warm component, which yielded a lower limit on a_min. These approaches allowed us to derive a mean value for the minimum grain size, along with its upper and lower boundaries. For targets, where a warm component fit was not feasible, only lower limits on the minimum grain size were obtained.

5.5 Results of SD modeling

Figure E.1 (bottom row) presents examples of SEDs fitted with SD models, featuring either one- or two-component configurations. The results of SD modeling are presented in Table 4, which provides the best-fitting parameters for the targets with reliable fits. These include the minimum grain radius, a_min, the SD power-law index, q, the BB temperature corresponding to the peak of the fitted emission, $T_{\text{dust}}^{SD}$, the IR excess, $f_{\text{disk}}^{SD}$, of the resolved cold belts, and their dust masses, $M_{\text{dust}}^{SD}$, integrated over grain sizes up to 5000 μm and assuming a pure astrosilicate composition. Additionally, the table lists the radiation pressure blowout size a_blow, which represents the grain size threshold below which particles are expected to be expelled from the disk by stellar radiation pressure. Dust grains become unbound when the radiation pressure force, F_rad, acting on them exceeds half of the gravitational force, F_gr (e.g., Krivov 2010). Based on the condition F_rad/F_gr = 0.5, the blowout grain size was calculated for a star with mass M_⋆ and luminosity L_⋆ following the formulation by Burns et al. (1979): $$a_{\text{blow\hspace{0.17em}}}=\frac{3L_{\star }⟨Q_{\mathrm{p}\mathrm{r}}⟩}{8\pi G{M}_{\star }c\varrho },$$

where ⟨Q_pr⟩ is the mean radiation pressure coupling coefficient averaged over the stellar flux (e.g., Augereau et al. 1999), G is the gravitational constant, and c is the speed of light. We assume ⟨Q_pr⟩ = 1 which is close to the values of 1.1 and 1.4 for stars with luminosities of L_⋆ = 10 L_⊙ and L_⋆ = 1 L_⊙, respectively, for astrosilicate composition (Pawellek et al. 2014). For low-luminosity stars (L_⋆ < 1 L_⊙), the blowout grain size is not calculated, as the radiation pressure is too weak to serve as an effective dust removal mechanism in such systems.

The derived minimum grain sizes, or their lower limits in cases where the warm component could not be fitted, are presented in Fig. 15, while the ratios of these values to the blowout grain sizes are shown in Fig. 16 as a function of stellar luminosity. Notably, nearly all derived a_min values are close to, but consistently larger than, the corresponding blowout sizes a_blow. This is expected, as grains smaller than a_blow are efficiently removed from the systems on short timescales due to radiation pressure. Figure 16 reveals that the ratio a_min/a_blow tends to increase with decreasing stellar luminosity, reaching values of up to ~9 in the case of HD 105. This trend corroborates earlier findings (see Pawellek et al. 2014; Pawellek & Krivov 2015) and could have several, not mutually exclusive, explanations: limitations in grain surface energy that inhibit the production of small collisional fragments (Krijt & Kama 2014; Thebault 2016) or lower dynamical excitation levels in cold disks, leading to a depletion of small grains (Thébault & Wu 2008).

For nearly all systems, the derived q values (Col. 4 in Table 4) are close to the canonical value of 3.5, which is expected for an idealized infinite, self-similar collisional cascade (Dohnanyi 1969). The mean value of q = 3.62 is slightly higher (Fig. 17), which is consistent with expectations for more realistic collisional systems, where the critical specific energy required for fragmentation increases with decreasing particle size in the strength-dominated regime (O’Brien & Greenberg 2003).

Figure 18 presents the belt dust masses derived from the SD model as a function of the measured disk radius, in order to examine the correlation between these two quantities. The contour plot in the figure represents the PDF for a subsample of resolved single belts around A-, F- and G-type stars with estimated ages between 10 and 200 Myr. The double-belt data points (red open circles in Fig. 18) are excluded from the PDF calculation due to the assumed relationship between the masses of two components in double-belt systems (Sect. 5.4).

We observe a tendency for increasing belt dust mass with growing radial distance from the star, following the fitted radial power law $M_{\text{dust}}^{\mathrm{S}\mathrm{D}}\propto R_{\text{belt\hspace{0.17em}}}^{{\alpha }_{\mathrm{R}}}$ with α_R = 2.1 ± 0.4. Although subject to small-number statistics bias, this trend is consistent with our findings for a larger sample of A- and F-type stars, based on the results of MBB modeling. The actual radial dependence may be less steep than our result suggests if single belts consist of multiple components that remain undetected due to insufficient spatial resolution. In such a case, the total cold dust mass would be distributed across multiple components, thereby reducing the mass assigned to a single belt and leading to a shallower mass distribution.

HD 107146. To investigate whether the non-detection of a disk could be attributed to the low reflectivity of its dust grains, we applied the SD model to fit the SED of the HD 107146 debris disk. This nearly pole-on (i = 19°) disk consists of two broad cold planetesimal belts located at ~50 and 120 au from its G-type host star, as previously observed with ALMA (Marino et al. 2018). The debris belts were not clearly detected in the IRDIS H-band polarimetric observations. To compare the dust optical properties inferred from SD modeling for this disk with those of a detected debris disk (see Sect. 6.3), we derived the best-fit parameters for HD 107146 as well. The results are listed in the last rows of Table 4.

6 Detections versus non-detections

To comprehensively understand the diversity of debris disk architectures, it is essential not only to analyze systems where belts are detected in scattered light, but also to interpret cases of non-detections, which provide valuable constraints. They may indicate disks that are intrinsically fainter, narrower, more poleon, or more evolved, and thus consistent with lower dust masses or smaller planetesimal belts located closer to the star. Non-detections in scattered light may still host massive disks visible in the IR, and their lack of scattered-light visibility could be explained by viewing geometry or dust properties. Additionally, the observing conditions during the runs, such as atmospheric seeing, coherence time, or instrument stability, can significantly influence detection sensitivity, particularly for faint disks.

Two-thirds of the debris disks in our sample were not detected in SPHERE observations. Therefore, in the following sections, we use the results of disk and SED modeling to explore the potential causes of these non-detections in more detail. Particular attention is given to the optical properties of dust grains, focusing on two complementary approaches for deriving dust scattering characteristics, such as albedo and maximum polarization fraction: one based on theoretical predictions using Mie theory (Mie 1908), and the other on measurements of total and polarized scattered fluxes.

6.1 Disk luminosity and dust mass

Debris disks with an IR excess below 10⁻⁴ are very faint, making them challenging to image in scattered light with current instrumentation. The faintest debris disk successfully imaged with SPHERE is that of HD 141943 (Boccaletti et al. 2019), with a disk fractional luminosity of f_disk = 1.2 × 10⁻⁴.

As illustrated in Fig. 12, where targets with detected disks are represented by filled circles, stellar age appears to be one of the most critical factors limiting the detectability of scattered light from debris disks around distant stars. Among all disks detected with SPHERE, 90% of the targets have a mean estimated age below 50 Myr. The stellar age histogram in Fig. 1 (lower left panel) also shows that debris disks are detected in more than 50% of systems younger than 100 Myr, whereas for older stars, the detection rate drops to just 5%. This declining detection rate with increasing stellar age can be attributed to the progressive reduction in dust mass over time, leading to a decrease in both scattered and thermal emission from dust particles (e.g., Wyatt 2008; Krivov 2010).

The relative small number of detected debris disks older than 50 Myr in our sample, which spans system ages from a few Myr to several Gyr, may provide constraints on the size of the largest planetesimals formed by the end of the PPD phase. To address the issue of overly high inferred debris disk masses, Krivov & Wyatt (2021) proposed that young debris disks may contain relatively small largest planetesimals (on the order of 1 km in size), suggesting that “planetesimals are born small”. The detection statistics in our sample support this hypothesis, as disks formed with small planetesimals are expected to appear bright at young ages (tens of Myr) but fade rapidly within a few hundred Myr, whereas disks formed with larger planetesimals would maintain their brightness over several Gyr.

6.2 Disk geometry and observing techniques

In addition to system age, the viewing geometry of the disk and the observing technique significantly influence the detectability of debris disks through DI. In particular, disk inclination often plays a decisive role, as highly inclined (nearly edge-on) systems are generally easier to detect (e.g., Esposito et al. 2020). Such disks exhibit increased SB along the edges of the planetesimal belt, where the column density of dust particles is highest, and on the disk’s front side, where forward-scattering enhances the intensity of scattered light. These effects make inclined disks more readily observable, whereas pole-on systems, which lack strong forward-scattering features and appear more diffuse, are inherently more challenging to detect.

To illustrate this effect, we generated model images of a typical debris disk observed at inclinations 0° (pole-on), 45° and 90° (edge-on). For this purpose, we used the model described in Sect. 4.5.1 adopting parameters commonly found in disk studies. Figure 19 presents the corresponding scattered light images, with total intensity shown in the top row and polarized intensity in the bottom row. When viewed edge-on (top left panel), the disk exhibits the highest SB, significantly enhancing its detectability.

Additionally, the observing technique most commonly used for imaging debris disks in scattered light is the ADI (Marois et al. 2006). This method is particularly sensitive to edge-on systems, as they generate a signal that differs from the stellar PSF when the sky field rotates. For disks with lower inclinations (i ≲ 70°), ADI is less effective, and for rotationally symmetric pole-on disks, it is inapplicable. As a result, despite their high IR excesses (f_disk > 10⁻³), the debris disks around HD 107146 and HD 95086 remained undetected in the SPHERE ADI datasets. Nonetheless, a large sky rotation angle during observations under good conditions can improve the detectability of low-inclination disks. This is demonstrated in the case of HD 105 debris disk (Fig. D.1), which was successfully detected, even though it has a relatively low inclination of 50.5°.

For imaging debris disks at low inclinations, PDI is more suitable than ADI. PDI yields images of the polarized intensity of scattered light, referred to as polarimetric images, which complement total intensity images of the same disk. However, the SB distribution in polarimetric images differs from that in total intensity images, as the two are governed by distinct phase functions: the SPF and the pSPF, respectively. This distinction is illustrated in Fig. 19, which shows simulated images of total and polarized intensities for a debris disk generated using our model with a parametric representation of the phase functions: a HG function with an asymmetry parameter of g = 0.6 for the SPF, and the function given by Eq. (2) with a maximum polarization fraction p_max = 0.3 for the pSPF.

The total flux, or integrated intensity, defined as the sum over all image pixels containing disk emission, also differs between total scattered and polarized intensity images and varies systematically with disk inclination. This effect is demonstrated in Fig. 20, which is based on the same model as in Fig. 19. The left panel of Fig. 20 shows the total scattered flux, F_sca, while the middle panel displays the total polarized flux, F_pol, both plotted as functions of the scattering asymmetry parameter g and disk inclination. Each flux is normalized by a factor of 4π/L_sca, where L_sca denotes the scattered luminosity of the disk, or total intensity integrated over the full solid angle.

As expected, the total scattered flux reaches a maximum for an edge-on disk with g = 0.9, corresponding to strongly forward-scattering grains (Fig. 20 left panel). The maximum polarized flux is obtained for a pole-on disk with an isotropic scattering parameter of g = 0 (Fig. 20 middle panel). In this configuration, most scattering occurs at θ = 90°, where, according to our model, the degree of linear polarization reaches its maximum (Eq. (3)). However, the polarized scattered intensity represents only a fraction of the total scattered intensity (Fig. 20 right panel), and this fraction can decrease rapidly with lower disk inclination, particularly if the dust particles exhibit strong forward-scattering behavior. A comparison of the model images in Figs. 19 and 20 demonstrate that debris disks appear faintest in polarized intensity when viewed pole-on, and that they consistently exhibit lower brightness in polarized light than in total scattered light, irrespective of the dust’s optical properties.

Despite this, thirty six debris disks in our sample were successfully detected using the PDI modes of IRDIS and ZIMPOL (Fig. F.1). The PDI data processing methodology enables a more effective subtraction of the stellar PSF from the science frames compared to the ADI technique. As a result, polarimetric images can achieve higher contrast levels, reaching up to 10⁻⁸−10⁻⁷ with ZIMPOL (Hunziker et al. 2020; Tschudi et al. 2024). This allows for the easy detection of young, low-inclination disks with a high polarization fraction of scattered light, where dust particles are confined to a narrow, bright ring such as HD 181327 planetesimal belt (Fig. 3; Milli et al. 2024). Conversely, older low-inclination disks, which exhibit a broad radial dust distribution in some ALMA images and may consist of multiple faint planetesimal rings, are more challenging to resolve using PDI. Broad pole-on disks generally have a lower surface density resulting in lower SB in polarized light. This could explain the marginal detection of debris belts in HD 107146 using IRDIS H-band polarimetry, despite relatively favorable observing conditions.

6.3 Influence and derivation of dust albedo and polarization characteristics

The optical properties of dust particles, specifically their scattering and absorption efficiencies, can also be responsible for a non-detection of a debris disk. If the dust scattering efficiency is low at a particular wavelength, the SB of the debris disk is correspondingly low at that wavelength, decreasing the probability of disk detection.

The optical characteristics of dust grains are fundamental parameters in the study of circumstellar environments. They are intrinsically linked to the grains’ composition, structure, and SD, and thus provide indirect constraints on the primordial solid-phase reservoir from which exoplanets may have formed. Understanding these properties is therefore a key objective in debris disk research.

One approach to constraining dust composition is through the evaluation of the dust albedo, which quantifies the relative efficiency of scattering versus total extinction (scattering plus absorption). A higher albedo than 0.5 indicates that scattering dominates over absorption, while a lower albedo suggests that absorption is more significant. If the dust grain SD within a disk is known, often inferred from modeling the SED (e.g., Pawellek et al. 2019) or scattered-light imaging (e.g., Olofsson et al. 2016), then theoretical predictions of albedo can be made for various dust compositions using Mie theory (Mie 1908) or more complex models accounting for grain porosity and non-sphericity (Draine & Flatau 1994; Min et al. 2005).

These theoretical predictions can then be compared to observational constraints, derived from combined measurements of scattered stellar light and thermal re-emission (SED). By jointly analyzing these datasets, the range of plausible dust compositions can be narrowed down. For example, grains composed primarily of astronomical silicates, carbonaceous materials, or ices each exhibit distinct scattering and absorption efficiencies, and thus different albedo values. This comparison allows to exclude certain grain compositions and structures, providing a more refined picture of the physical nature of dust in debris disks.

In the following section, we present the methodology used to calculate the albedo of dust grains in debris disks, employing Mie theory to derive the scattering and absorption efficiencies for particles of specified composition and size. We then describe how these theoretical predictions are compared with observational constraints, obtained from the analysis of scattered-light images of spatially resolved debris disks in our sample and their IR excess measurements.

6.3.1 Albedo

The amount of stellar photons with a wavelength λ which are scattered by a spherical dust particle of radius a is determined by its spectral cross section for scattering, $C_{\lambda }^{\text{sca}}$. This cross section quantifies the relationship between the intensity of the incident radiation, I_{⋆ λ}, and the scattered power or spectral radiant flux, $F_{\lambda }^{\text{sca}}$ (Bohren & Huffman 1983): $$C_{\lambda }^{\mathrm{s}\mathrm{c}\mathrm{a}}=\frac{F_{\lambda }^{\mathrm{s}\mathrm{c}\mathrm{a}}}{I_{\star \lambda }}=Q_{\lambda }^{\mathrm{s}\mathrm{c}\mathrm{a}}A,$$

where $Q_{\lambda }^{\mathrm{s}\mathrm{c}\mathrm{a}}$ is the scattering efficiency and A = πa² is the geometrical cross section of particle.

The fraction of stellar light attenuation caused by scattering, and thus the role of scattering in the overall extinction process, is characterized by the single scattering albedo ω_λ, which describes the proportion of extinction resulting from scattering rather than absorption: $${\omega }_{\lambda }=\frac{C_{\lambda }^{\mathrm{s}\mathrm{c}\mathrm{a}}}{C_{\lambda }^{\mathrm{e}\mathrm{x}\mathrm{t}}}=\frac{C_{\lambda }^{\mathrm{s}\mathrm{c}\mathrm{a}}}{C_{\lambda }^{\mathrm{s}\mathrm{c}\mathrm{a}}+C_{\lambda }^{\mathrm{a}\mathrm{b}\mathrm{s}}}=\frac{Q_{\lambda }^{\mathrm{s}\mathrm{c}\mathrm{a}}}{Q_{\lambda }^{\mathrm{s}\mathrm{c}\mathrm{a}}+Q_{\lambda }^{\mathrm{a}\mathrm{b}\mathrm{s}}}=\frac{F_{\lambda }^{\mathrm{s}\mathrm{c}\mathrm{a}}}{F_{\lambda }^{\mathrm{s}\mathrm{c}\mathrm{a}}+F_{\lambda }^{\mathrm{a}\mathrm{b}\mathrm{s}}},$$

where $C_{\lambda }^{\mathrm{e}\mathrm{x}\mathrm{t}}=Q_{\lambda }^{\mathrm{e}\mathrm{x}\mathrm{t}}A=(F_{\lambda }^{\mathrm{s}\mathrm{c}\mathrm{a}}+F_{\lambda }^{\mathrm{a}\mathrm{b}\mathrm{s}})/I_{\star \lambda }$ and $C_{\lambda }^{\mathrm{a}\mathrm{b}\mathrm{s}}=Q_{\lambda }^{\mathrm{a}\mathrm{b}\mathrm{s}}A=F_{\lambda }^{\text{abs\hspace{0.17em}}}/I_{\star \lambda }$ are the spectral extinction and absorption cross sections, respectively, and $F_{\lambda }^{\text{abs}}$ is the power of light absorbed by dust particle.

For a given dust composition characterized by a complex refractive index, Mie theory predicts the scattering efficiency $Q_{\lambda }^{\text{sca}}$ and extinction efficiency $Q_{\lambda }^{\text{ext}}$ as functions of the size parameter x (Bohren & Huffman 1983): $$x=\frac{2\pi a}{\lambda }.$$(9)

For dust particles following a SD characterized by a differential grain number density n(a)¹⁰, the effective size parameter can be defined as the ratio of the third to the second moment of the distribution¹¹ (Hansen & Travis 1974): $$x_{\mathrm{e}\mathrm{f}\mathrm{f}}=\frac{2\pi }{\lambda }{a}_{\mathrm{e}\mathrm{f}\mathrm{f}}=\frac{2\pi }{\lambda }\frac{\int a^{3}n(a)da}{\int a^{2}n(a)da}.$$

The total cross sections for scattering ${\sigma }_{\lambda }^{\text{sca}}$ and extinction ${\sigma }_{\lambda }^{\text{ext}}$ are obtained by averaging over the distribution: $${\sigma }_{\lambda }^{\mathrm{s}\mathrm{c}\mathrm{a}}=\frac{\int C_{\lambda }^{\mathrm{s}\mathrm{c}\mathrm{a}}(a)n(a)da}{\int n(a)da}{\sigma }_{\lambda }^{\mathrm{e}\mathrm{x}\mathrm{t}}=\frac{\int C_{\lambda }^{\mathrm{e}\mathrm{x}\mathrm{t}}(a)n(a)da}{\int n(a)da}.$$

In this case the spectral albedo is given by $${\omega }_{\lambda }=\frac{\int C_{\lambda }^{\mathrm{s}\mathrm{c}\mathrm{a}}(a)n(a)da}{\int C_{\lambda }^{\mathrm{e}\mathrm{x}\mathrm{t}}(a)n(a)da}.$$(10)

6.3.2 Scattering albedo of various dust compositions

Using Eq. (10), we modeled the spectral albedo for a range of dust compositions, including astrosilicates, amorphous carbon, and silicate grains coated with either dirty or water ice, motivated by both observational evidence and theoretical considerations. Astrosilicates are widely used to represent the dominant silicate emission features observed in mid-IR spectra of circumstellar environments and reflect the mineralogical composition inferred from both debris disks and Solar System dust populations (e.g., Draine 2003a; Dorschner et al. 1995). Amorphous carbon is included to represent more absorbing, featureless materials, commonly invoked to model the continuum emission in disks and supported by the presence of carbonaceous material in interplanetary dust and meteorites (e.g., Zubko et al. 1996; Li & Greenberg 1997). To account for conditions in the outer, colder regions of disks, we also considered silicate grains coated with water ice, which are expected beyond the ice line and significantly modify scattering properties due to their high albedo and distinct optical constants (e.g., Donaldson et al. 2013; Xie et al. 2025). Additionally, dirty ice grains, incorporating refractory inclusions, provide a more realistic representation of ice mantles processed by collisions and irradiation (e.g., Preibisch et al. 1993; Li & Greenberg 1998). This set of compositions spans a physically plausible range and enables the exploration of how material properties influence key observables such as scattered light brightness.

For these four dust compositions, we computed the scattering albedo ω_λ and present the results as 2D maps in Fig. 21. The calculations are based on a grain SD model n(a)da ∝ a^−qda, with grain sizes spanning the range a_min ⩽ a ⩽ a_max. We varied a_min between 0.9 and 5 μm while keeping a_max fixed at 5 mm, following our SED fitting procedure (Sect. 5.4). We considered three values for the SD power-law exponent, q = 3.0, 3.5 and 4.0, assuming that for most debris systems, the exponent falls within this range (see Col. 4 in Table 4).

The optical data for astrosilicates were taken from Draine (2003a), for water ice from Warren & Brandt (2008), and for amorphous carbon and dirty ice from Preibisch et al. (1993). The refractive indices of silicate grains coated with water or dirty ice were calculated assuming a volume fraction of 50% for each component. This corresponds to a mass fraction of 79% for silicates and 21% for water ice, based on their material densities of ϱ_sil = 3.5 g cm⁻³ and ϱ_ice = 0.92 g cm⁻³, respectively (Pollack et al. 1994). The optical constants of the dirty ice coating were derived for a mixture of H₂O- and NH₃-ices with a volume ratio of 3:1 polluted by amorphous carbon with a volume fraction of 10%. Adopting the material densities for carbon ϱ_car = 2.3 g cm⁻³ and for NH₃-ice ϱ_{NH3 ice} = 0.85 g cm⁻³, we obtained a mass fraction of 23% for the dirty ice mantle.

The left column in Fig. 21 presents the single scattering albedo for particles composed of pure astrosilicates. For q = 3.0 the albedo remains nearly constant (~0.56) for all SDs and wavelengths, corresponding to a gray disk, meaning that the reflectance spectrum of the disk does not vary with wavelength. For q = 4.0, the spectral variation of ω_λ is more significant, with the albedo increasing towards the lower right corner of the plot as the effective size parameter decreases x_eff. This trend is due to the decreasing effective grain size a_eff and the increasing wavelength λ. In the considered range of x_eff the scattering cross section of astrosilicates is larger for smaller values of effective size parameter. Within the wavelength range covered by SPHERE filters (0.5–2.25 μm), the scattering cross section of dust grains increases with wavelength, resulting in a higher albedo and thereby an increase in the relative scattered flux (i.e., the disk flux normalized by the stellar flux). Consequently, a disk composed of such dust particles is expected to exhibit a red reflectance spectrum¹², which becomes more pronounced for dust compositions with a higher fraction of small particles (i.e., smaller x_eff or larger q).

For dust particles composed of amorphous carbon and following a SD n(a) ∝ a⁻³ (top panel of second column in Fig. 21) the reflectance spectrum color remains essentially unchanged regardless of variations in the minimum grain size or wavelength. In this case, the bulk albedo is ~0.63. In contrast to astrosilicates, the albedo decreases for steeper SDs, reaching a value of ~0.52 for a distribution with a minimum grain size of a_min = 0.9 μm (bottom panel of second column in Fig. 21). However, even for steeper SDs, the bulk albedo exhibits minimal spectral variation, meaning that a disk composed of such dust particles would appear gray to an observer using SPHERE, irrespective of the specific SD parameters.

The albedo maps of astrosilicate particles coated with either dirty or pure water ice reveal distinct water ice absorption features (third and fourth columns in Fig. 21), with the most prominent one at 3.1 μm attributed to O-H stretching vibrations of water ice. The spectral albedo of a mixture containing pure water ice is significantly higher, reaching up to 0.88, compared to that of pure astrosilicates. When observed using the broadband H filter with IRDIS, a disk composed of such icy grains would exhibit up to twice the scattered flux of a similar disk with identical viewing geometry and stellar irradiance but composed of amorphous carbon or even astrosilicate particles.

Based on the results of the SED fitting with the SD model (Cols. 3 and 4 in Table 4), we calculated the range of possible spectral albedo values for disks observed at λ = 1.6 μm (central wavelength of the IRDIS broadband H filter). For all resolved exo-Kuiper belts listed in Table 4, the albedo values lie within the range [0.54, 0.68], assuming an astrosilicate dust composition. Higher albedo values are obtained for SDs with smaller minimum grain sizes. Given the narrow range of derived albedo values, it is unlikely that dust albedo is the primary factor behind the disk non-detections.

6.3.3 Detection of the HD 181327 debris belt versus non-detection of the inner belt around HD 107146

The difference in scattering efficiency of the dust material may explain the detection of the HD 181327 debris belt and the non-detection of the inner belt around HD 107146 in the H-band polarimetric observations with IRDIS. Both stars were observed using the same instrumental setup, under comparable observing conditions, and with nearly identical total exposure times.

The F6V star HD 181327 (L_⋆ = 2.88 L_⊙, 18–23 Myr, f_disk = (2.6 ± 0.7) × 10⁻³) hosts a debris belt at a radial distance of 82 au, inclined at 30° (Table 2). As mentioned in Sect. 5.5, the G2V star HD 107146 (L_⋆ = 1.04 L_⊙, 50–2400 Myr, f_disk = (1.1 ± 0.3) × 10⁻³) possesses two low-inclination (i = 19°) planetesimal belts located at ~50 and 120 au. The stellar illumination of the HD 181327 belt is comparable to that of the inner belt of HD 107146, as it is governed by the ratio $L_{\star }/R_{\text{belt}}^2$, which yields a value of 0.55 W m⁻² in both cases. This is two orders of magnitude lower than the solar flux received by Jupiter in the Solar System. Additionally, the small difference in inclination between the two planetesimal belts is not expected to significantly impact the amount of observed polarized scattered light, as seen in the middle panel of Fig. 20. Thus, the stellar illumination and disk viewing geometry are not the primary factors responsible for the non-detection of the 50 au belt in the HD 107146 system, or at least, they do not play a decisive role.

Using the results of SD modeling, we estimated the scattering efficiencies of dust grains in both systems assuming they are composed of astrosilicates. Interestingly, for both systems we obtained the same SD power-law index of 3.42 but different minimum grain sizes: a_min = 1.00 ± 0.18 μm for the HD 181327 belt and a_min = 2.79 ± 1.17 μm for the HD 107146 belt. Such a difference in minimum grain sizes would lead to different bulk albedo values in the H band (λ_c = 1.6 μm): ω_H = 0.61 for HD 181327 and ω_H = 0.56 for HD 107146. Consequently, the polarized scattered flux from the HD 181327 belt could be 1.1 times higher than that from the HD 107146 inner belt, assuming a comparable number of scattering particles in both belts.

However, it is entirely possible that the dust particles in the HD 107146 belt have a different composition, with lower scattering efficiency or a lower maximum polarization fraction of scattered light compared to the dust around HD 181327. Moreover, the dust spatial distribution and total mass may be the key factors contributing to the non-detection. Despite both systems exhibiting a high IR excess, the majority of dust in the HD 181327 system is confined to one relatively narrow debris belt, whereas in HD 107146, the total dust mass is distributed across at least two cold belts. According to our SD modeling results, the dust mass of the inner belt in HD 107146, and therefore the number of scattering particles, is estimated to be approx. 5.5 times lower than that of the HD 181327 debris belt (Col. 7 in Table 4). ALMA observations of both disks (Marino et al. 2016, 2018) show that the belt area of the HD 181327 disk (R_belt × ΔR_belt = 86 × 23.2 = 2 × 10³ au²) is larger than the area of the inner belt in HD 107146 (R_belt × ΔR_belt ≈ 50 × 30 = 1.5 × 10³ au²). This implies that the dust surface density, and consequently the disk’s SB in scattered light, may be up to four times higher for the HD 181327 belt, making its detection in polarized light with SPHERE significantly more likely.

6.3.4 Parametric approach for deriving optical properties of dust grains

In the following sections, we introduce a new diagnostic approach for deriving the bulk albedo and maximum polarization fraction of dust grains based on scattered-light and polarized-light images combined with parametric modeling. This method provides independent albedo estimates that can be directly compared with the values obtained from Mie theory discussed in Sects. 6.3.1 and 6.3.2.

In most cases, the single scattering albedo of disk material cannot be directly retrieved from disk images. This is because the observer measures only a fraction of the total scattered flux, which is governed by the material’s SPF or, equivalently, its differential scattering cross section $d{\sigma }_{\lambda }^{\mathrm{s}\mathrm{c}\mathrm{a}}/d\mathrm{\Omega }$. This quantity describes the fraction of incident light scattered into a specific direction per unit solid angle Ω, so that $${\sigma }_{\lambda }^{\mathrm{s}\mathrm{c}\mathrm{a}}=\int \frac{d{\sigma }_{\lambda }^{\mathrm{s}\mathrm{c}\mathrm{a}}}{d\mathrm{\Omega }}d\mathrm{\Omega }=\int {\sigma }_{\lambda }^{\mathrm{s}\mathrm{c}\mathrm{a}}\times SPFd\mathrm{\Omega }.$$

Additionally, the total scattered flux, F_{sca λ}, measured in disk images¹³ is an integrated flux from disk regions with varying scattering angles θ, and is further influenced by the disk’s viewing geometry (Sect. 6.2). This geometry is characterized by parameters such as the disk inclination, i, the radius of the planetesimal belt, R_belt, the disk opening angle, H₀/r₀, and the exponents of the radial power laws, α_in and α_out (see Sect. 4.5.1): $$F_{\mathrm{s}\mathrm{c}\mathrm{a}\lambda }=f(\frac{d{\sigma }_{\lambda }^{\mathrm{s}\mathrm{c}\mathrm{a}}}{d\mathrm{\Omega }},i,R_{\mathrm{b}\mathrm{e}\mathrm{l}\mathrm{t}},H_0/r_0,{\alpha }_{\mathrm{i}\mathrm{n}},{\alpha }_{\mathrm{o}\mathrm{u}\mathrm{t}})$$

Nevertheless, optical properties of the dust in debris disks can be constrained with a comparison between the amount of observed scattered radiation and the IR excess, which represents very roughly the relation between dust scattering and absorption. To examine this relation, several debris disk studies have estimated the so-called disk single scattering albedo ω_{disk λ}. This was done by computing the ratio of the total scattered flux derived from the disk image, F_{sca λ}, to the disk’s IR excess, L_{IR disk}/L_⋆, (e.g., Schneider et al. 2014; Choquet et al. 2018; Engler et al. 2023) according to equation $${\omega }_{\mathrm{disk}\lambda }=\frac{F_{\mathrm{s}\mathrm{c}\mathrm{a}\lambda }/F_{\star \lambda }}{F_{\mathrm{s}\mathrm{c}\mathrm{a}\lambda }/F_{\star \lambda }+L_{\mathrm{I}\mathrm{R}\mathrm{disk}}/L_{\star }},$$(11)

where F_{⋆ λ} denotes the stellar flux at wavelength λ, and L_{IR disk}/L_⋆ is used as a proxy for the absorbed flux at that wavelength. The latter represents a rough approximation, as the IR excess reflects the total emission integrated over the entire IR spectrum, and the dust generating the thermal flux may be located not only at the position of planetesimal belt resolved in scattered light.

When calculated in this manner, the disk scattering albedo is proportional to the dust albedo ω_λ (Eq. (10)). However, it also depends on the spatial distribution of dust particles in the disk, as discussed above: $${\omega }_{\text{disk\hspace{0.17em}}\lambda }\propto {\omega }_{\lambda }\cdot f(\frac{d{\sigma }_{\lambda }^{\text{sca\hspace{0.17em}}}}{d\mathrm{\Omega }},i,R_{\text{belt\hspace{0.17em}}},H_0/r_0,{\alpha }_{\text{in\hspace{0.17em}}},{\alpha }_{\text{out\hspace{0.17em}}}).$$

Therefore, this type of disk albedo can be useful for comparing the scattering properties of two debris disks with the same viewing geometry and stellar irradiance and, ideally, the same SPF. In any case, the measured scattered flux should always be corrected for flux losses introduced by the post-processing of ADI data to ensure accurate comparisons.

The single scattering albedo of dust material can be determined by considering the full angular distribution of scattered light. This requires disk modeling to estimate the scattered spectral luminosity of the disk L_{sca λ}. Once this value is obtained, the dust albedo can be derived using the following relation: $${\omega }_{\lambda }\approx \frac{L_{\mathrm{s}\mathrm{c}\mathrm{a}\lambda }/L_{\star \lambda }}{L_{\mathrm{s}\mathrm{c}\mathrm{a}\lambda }/L_{\star \lambda }+L_{\mathrm{I}\mathrm{R}\text{\hspace{0.17em}disk\hspace{0.17em}}}/L_{\star }}=\frac{⟨F_{\mathrm{s}\mathrm{c}\mathrm{a}\lambda }⟩/F_{\star \lambda }}{⟨F_{\mathrm{s}\mathrm{c}\mathrm{a}\lambda }⟩/F_{\star \lambda }+L_{\mathrm{I}\mathrm{R}\text{\hspace{0.17em}disk\hspace{0.17em}}}/L_{\star }},$$(12)

where ⟨F_{sca λ}⟩ = L_{sca λ}/4π is the scattered disk flux averaged over the full solid angle.

The ratio between the observed disk flux F_{sca λ} and the averaged flux ⟨F_{sca λ}⟩ can be determined if the shape of the SPF is known, for example, from the disk image model. To demonstrate it, we plot in Fig. 22a the ratio F_{sca λ}/⟨F_{sca λ}⟩ for different HG functions and disk inclinations. By applying a correction factor for the disk inclination, hereafter referred to as the view factor for scattered flux and defined as f_{sca λ} = F_{sca λ}/⟨F_{sca λ}⟩, we obtain $${\omega }_{\lambda }\approx \frac{F_{\mathrm{s}\mathrm{c}\mathrm{a}\lambda }/F_{\star \lambda }}{F_{\mathrm{s}\mathrm{c}\mathrm{a}\lambda }/F_{\star \lambda }+f_{\mathrm{s}\mathrm{c}\mathrm{c}\lambda }\cdot L_{\mathrm{I}\mathrm{R}\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{k}}/L_{\star }}.$$(13)

The averaged scattered flux ⟨F_{sca λ}⟩ is equal to the measured disk flux for all inclinations only in the case of isotropic scattering (g = 0). As shown in Fig. 22a (see also Schmid 2021), for each HG parameter g, there is a specific disk inclination where f_{sca λ} = 1. In such a case, e.g., for a disk inclined at 60° with a HG parameter g = 0.6, the measured scattered flux does not require any correction, and the single scattering dust albedo can be approximated as ω_λ ≈ ω_{disk λ}.

The view factor f_{sca λ} depends only on the shape of the SPF and the disk inclination, making it independent of other disk geometrical parameters. If the scattered-light image is modeled using a single HG function (i.e., without a combination of multiple HG functions), the view factor f_{sca λ} can be directly obtained from Fig. 22a.

6.3.5 Disk albedo with polarized flux. The Λ parameter

In Equation (13), the scattered flux F_{sca λ} can be replaced by the measured polarized flux F_{pol λ} using another multiplicative factor f_{pol λ} = F_{pol λ}/⟨F_{sca λ}⟩ (see Sect. 6.3.6) to derive the scattering albedo of the disk: $${\omega }_{\lambda }\approx \frac{F_{\mathrm{p}\mathrm{o}\mathrm{l}\lambda }/F_{\star \lambda }}{F_{\mathrm{p}\mathrm{o}\mathrm{l}\lambda }/F_{\star \lambda }+f_{\mathrm{p}\mathrm{o}\mathrm{l}\lambda }\cdot L_{\mathrm{I}\mathrm{R}\text{\hspace{0.17em}disk\hspace{0.17em}}}/L_{\star }}.$$(14)

Polarized flux is often easier to measure, particularly for debris disks with lower inclinations, and can be corrected for polarimetric signal losses more reliably (Engler et al. 2018).

Equation (14) can be reformulated in terms of the Λ parameter (see Eq. (17)). This parameter quantitatively characterizes the relationship between the measured polarized flux and the fractional IR luminosity of a debris disk (Engler et al. 2017): $${\mathrm{\Lambda }}_{\mathrm{p}\mathrm{o}\mathrm{l}\lambda }=\frac{F_{\mathrm{p}\mathrm{o}\mathrm{l}\lambda }/F_{\star \lambda }}{L_{\mathrm{I}\mathrm{R}}/L_{\star }}.$$(15)

Similar to the disk single scattering albedo, ω_{disk λ}, this observational parameter is proportional to the ratio of dust cross sections for scattering and absorption. Additionally, it depends on two key factors: the polarization fraction function, p_λ(θ), and the maximum polarization fraction of dust material, p_{max λ}.

Analogously to Eq. (15), we define the Λ parameters for the measured scattered flux F_{sca λ} and disk-averaged flux ⟨F_{sca λ}⟩ as $${\mathrm{\Lambda }}_{\mathrm{s}\mathrm{c}\mathrm{a}\lambda }=\frac{F_{\mathrm{s}\mathrm{c}\mathrm{a}\lambda }/F_{\star \lambda }}{L_{\mathrm{I}\mathrm{R}}/L_{\star }}\text{\hspace{0.17em}and\hspace{0.17em}}⟨{\mathrm{\Lambda }}_{\mathrm{s}\mathrm{c}\mathrm{a}\lambda }⟩=\frac{⟨F_{\mathrm{s}\mathrm{c}\mathrm{a}\lambda }⟩/F_{\star \lambda }}{L_{\mathrm{I}\mathrm{R}}/L_{\star }},$$(16)

respectively. Using these definitions, the Eqs. (13) and (12) can be re-expressed in terms of the Λ parameters as well.

6.3.6 Disk polarized contrast versus IR excess for the studied debris disk sample

To derive the Λ parameter for debris disks detected with SPHERE in polarimetric modes, we measured their polarized contrast relative to the star, given by F_{pol λ}/F_{⋆ λ}. Most of these disks were observed with IRDIS using broadband filters H and J, indicated in Fig. 23 by blue and red markers, respectively. This figure displays the measured polarized contrast as a function of disk fractional IR luminosity. The elliptical shape of each marker reflects the disk inclination. However, in the case of polarized contrast, the dependence on inclination is relatively weak, significantly less pronounced than for the total scattered flux (Figs. 22a, b).

Figure 23 reveals a positive correlation between polarized contrast and fractional IR luminosity. This correlation is expected and can be attributed to the dust scattering and dust absorption opacities which both depend mainly on the amount and spatial distribution of dust particles.

Although the number of disks with similar inclinations is limited, we fitted a simple linear relation between polarized contrast and fractional IR luminosity to the H-band data, focusing on disks with inclinations greater than 75° (which represent the majority of detected disks). For these disks, the polarized contrast follows the fractional IR luminosity, f_disk, according to the expression 0.074 f_disk −8.1 × 10⁻⁶, as shown by the blue dashed line in Fig. 23. This fit corresponds to a Λ_{disk H} parameter of 0.074 ± 0.007.

A few debris disks with inclinations lower than 75° observed in the H band, as well as most disks observed in the J band, exhibit a lower contrast and consequently a lower Λ parameter. The scatter seen in Fig. 23 may reflect variations in the scattering properties of dust grains among disks with similar inclinations or may result from differences in viewing geometry.

The dependence of the data on viewing geometry can be eliminated by dividing the disk’s polarized contrast by the view factor f_{pol λ} = F_{pol λ}/⟨F_{sca λ}⟩ (Fig. 22b). Similar to f_{sca λ}, the view factor for polarized flux, f_{pol λ}, is independent of the radial structure for axisymmetric dust distributions. However, its value is determined by the shape of the pSPF and scales linearly with the maximum polarization fraction of the dust material, p_{max λ}.

Once the Λ_{pol λ} parameter and the view factor f_{pol λ} are determined, the single scattering albedo of the dust material can be estimated as follows: $${\omega }_{\lambda }\approx \frac{{\mathrm{\Lambda }}_{\mathrm{p}\mathrm{o}\mathrm{l}\lambda }}{{\mathrm{\Lambda }}_{\mathrm{p}\mathrm{o}\mathrm{l}\lambda }+f_{\mathrm{p}\mathrm{o}\mathrm{l}\lambda }}.$$(17)

For example, for disks with inclinations greater than 75° observed in the H band, we derived a value of Λ_{pol H} = 0.074 ± 0.007. Assuming a maximum polarization fraction p_{max H} = 0.3 and a scattering asymmetry parameter g = 0.9 for these disks, we obtain a view factor of f_{pol λ} = 0.033 from Fig. 22b. Substituting these values into Eq. (17) yields a dust albedo of 0.69 ± 0.02. According to Fig. 21, such an albedo is consistent, for instance, with dust grains composed of astrosilicates, either pure or coated with a dirty ice mantle, and having a minimum size of a_min = 1 μm.

6.3.7 Constraining the maximum polarization fraction from combined scattered and polarized light observations

The combined analysis of total and polarized intensity images enables a more comprehensive assessment of dust scattering properties. Specifically, the ratio of polarized to total scattered flux provides an observational constraint on the maximum polarization fraction, p_max,λ, of the dust grains at a given wavelength.

This parameter represents the peak of the polarization fraction phase function p_λ(θ) and is sensitive to grain composition, porosity, and SD. As such, it plays a critical role, alongside the single scattering albedo ω_λ and the shape of the SPF, in constraining the optical constants and morphology of the dust population (e.g, Graham et al. 2007; Kirchschlager & Wolf 2014). Because different grain materials and structures (compact vs. aggregate particles) produce distinct polarization signatures, the combination of albedo and p_{max λ} allows for a significantly narrower range of viable dust models than can be achieved using albedo alone.

The analysis presented in this section requires accurately measured total scattered fluxes, which are best obtained from images processed with RDI, a technique that preserves the photometric fidelity of extended disk structures. Using this technique, we measured the total scattered flux of the HD 129590 debris disk (Olofsson et al. 2024) in the H band. This is the only IRDIS image in our sample from which the total scattered flux could be reliably extracted. To broaden the sample and enable a meaningful comparison, we complemented our analysis with scattered light measurements obtained from HST observations at optical wavelengths, where RDI processing is routinely applied and photometric calibration is robust. In the following, we compare total intensity contrasts derived from HST data (Schneider et al. 2014) to the polarized intensity contrasts measured with sphere.

Using broadband optical¹⁴ images of ten debris disks obtained with the Space Telescope Imaging Spectrograph (STIS), Schneider et al. (2014) investigated the correlation between IR excess and the optical scattering fraction, analogous to the relation shown in our Fig. 23. Their result is reproduced in Fig. 24 as the black solid line (see also their Fig. 7). They derived a proportionality factor between the fractional scattered flux and the IR excess of Λ_{sca opt} ≈ 1.05 (according to Eq. (16)), which is significantly higher than the Λ_{pol H} parameter derived from our polarized flux analysis in the near-IR (Sect. 6.3.6). This difference is expected, as the polarized flux represents only a fraction of the total scattered flux (Sect. 6.2).

We measured the polarized flux for five of the ten HST targets using SPHERE: in the J band for HD 15115, HD 32297, and HD 197481 (AU Mic), and in the H band for HD 61005 and HD 181327. The HST measurements for these targets are represented by black ellipses in Fig. 24. The SPHERE measurement of total scattered flux in the H band for the HD 129590 debris disk is marked with an orange ellipse in the same figure. In total, we thus have six targets for which both total and polarized scattered fluxes are available. However, for five of these, the measurements come from different instruments, HST for total intensity and SPHERE for polarized intensity, while HD 129590 is the only case where both measurements were obtained from the same dataset.

To provide a general comparison between HST/STIS and SPHERE/IRDIS measurements, we performed a linear fit to the STIS data points, yielding a best-fit slope of Λ_{sca opt} = 0.56 ± 0.24. This fit, shown as the black dash-dotted line in Fig. 24, is consistent within a 2σ confidence interval with the trend reported by Schneider et al. (2014). The result suggests that the fractional scattered fluxes measured in optical total intensity images are, on average, approx. 7.5 to 10 times higher than those derived from near-IR polarized intensity images.

A comparison of the five individual targets for which both optical scattered and near-IR polarized fluxes are available shows that the optical fluxes exceed the polarized fluxes by factors ranging from 7 (HD 32297) to 33 (HD 197481). In terms of magnitudes, this corresponds to differences of 2.18^m for HD 32297 (A0V, ~30 Myr) and 3.80^m for HD 197481 (M1V, 18-23 Myr), both of which are edge-on systems. For the other three debris disks, the magnitude difference falls within a comparable range: 2.61^m for HD 181327 (F6V, i = 30°, 18-23 Myr), 2.81^m for HD 15115 (F4V, i = 85°, 10-500 Myr), and 3.09^m for HD 61005 (G8V, i = 82°, 45-55 Myr). For the IRDIS measurement of the HD 129590 disk (G3V, i = 81°, 14-18 Myr), the ratio between the total scattered and polarized fluxes is 16.9 ± 2.4 (3.07^m ± 0.14^m).

Using this ratio, we can estimate the value of maximum polarization fraction of the HD 129590 debris particles, as illustrated in Fig. 22c. In this figure, the blue arrow and shaded area represent the measured flux ratio and its uncertainty, while the orange arrow and shaded area indicate the measured disk inclination and its associated uncertainty. The intersection point of these two indicators should lie on the curve corresponding to the HG asymmetry parameter derived from the modeling of the disk image. In the case of HD 129590, the intersection falls on the curve with g = 0.6, which is the lower bound of the modeling result, suggesting a maximum polarization fraction of approx. 0.3. To test the compatibility of modeling result with higher values of p_max, the plotted curves can be rescaled. As seen in Fig. 22c, the measured flux ratio and inclination are also consistent with a solution adopting g = 0.7 and p_max = 0.4 (black dotted line), which is closer to the best-fit value for g obtained for this target.

A similar conclusion can be drawn from Fig. 25, where we plot the ratio of the measured polarized to scattered flux as a function of the disk IR fractional luminosity for systems with available flux measurements from both HST/STIS and SPHERE/IRDIS images of total scattered intensity, as presented in Fig. 24. As discussed above, the data point for HD 129590 (orange marker in Fig. 25) is the only case where both the polarized and scattered fluxes were derived from the same IRDIS H-band dataset. For the remaining five data points (black markers in Fig. 25) the total scattered flux was measured from HST/STIS, while the polarized flux was obtained from IRDIS polarimetric observations in either the J or H band. Consequently, when estimating the polarization fraction in the near-IR, one should consider that the positions of the black markers in this figure may shift, if the contrast in scattered light differs between the optical and near-IR wavelengths.

According to results presented in Fig. 25, the ratios of polarized to scattered flux are below 15% for all six considered debris disks, regardless of inclination or IR excess. These values are broadly consistent with the inverse of the modeled ratio F_{sca λ} / F_{pol λ}, assuming a maximum polarization fraction of p_max = 0.3 (see right panel in Fig. 20). This means, that, if the fractional scattered flux, or contrast, in the near-IR is the same as at the optical wavelengths, the near-IR maximum polarization fraction of scattered light for these disks might be in the range 25–35%.

The largest difference to the modeling results adopting p_max = 0.3 is observed for HD 181327. Two possible explanations may account for this deviation. First, if the scattered-light contrast is similar across optical and near-IR wavelengths, the HD 181327 disk may intrinsically exhibit a lower maximum polarization fraction than p_max = 0.3. This scenario can be readily assessed, as the modeled ratio F_{sca λ} / F_{pol λ} varies only marginally with g for disks inclined at 30° (see Fig. 22c). By rescaling this ratio, we estimate the location of HD 181327 disk in Fig. 25 for assumed values of p_max = 0.1 and p_max = 0.2, indicated by annotated green markers. The marker for p_max = 0.1 lies closest to the observed position of HD 181327, suggesting a maximum polarization fraction of ~12% for this disk. This value is about half of that estimated by Milli et al. (2024), who used HST/NICMOS F110W filter data to calibrate the total scattered flux inferred from the IRDIS total intensity image obtained using RDI technique, along with the same polarimetric dataset analyzed in this work. Therefore, an alternative explanation maybe valid, namely, that the optical scattered-light contrast measured by Schneider et al. (2014) is significantly higher than the near-IR scattered-light contrast. This system features a prominent halo of small grains, which is well resolved in the HST data due to the disk’s low inclination. The contribution from this extended halo likely inflated the total scattered flux measured in the optical, resulting in a lower estimated flux ratio that does not accurately reflect the maximum polarization fraction in the near-IR.

An advantage of using low-inclination disks (i < 40°) in this diagnostic is that the ratio F_{sca λ}/F_{pol λ} remains nearly constant across all values of the HG asymmetry parameter (see Fig. 22c). Consequently, the positions of these disks in the diagnostic diagram are primarily sensitive to the assumed maximum polarization fraction p_max and are largely independent of g. This is not the case for higher-inclination systems, such as HD 129590, where the ratio F_{sca λ}/F_{pol λ} shows a stronger dependence on both g and p_max. In Figure 25, for instance, the calculated position of the disk corresponding to p_max = 0.4 and g = 0.7 lies closer to the observed value than the one with p_max = 0.3 and the same g. However, a similar good agreement as for p_max = 0.4 can also be obtained with p_max = 0.3 and a lower asymmetry parameter of g = 0.6, since the modeled positions shift upward with decreasing g or increasing p_max. Although values of p_max > 0.4 may appear to better match the data for HD 129590, given its best-fit asymmetry parameter of g = 0.78 (see Table G.1), the actual maximum polarization fraction in this system could be lower than 0.4. As discussed in Sect. 4.4, HD 129590 exhibits a double-belt structure, with the inner belt being significantly brighter than the outer one. For simplicity, we modeled this system using a single-belt model, resulting in a best-fit radius that lies between the measured radii of the two components. The dominance of the bright inner belt may have biased the modeling toward a higher apparent value of g. Therefore, the true asymmetry parameter could be lower than the derived value, which in turn could imply a maximum polarization fraction below 0.4.

It is important to note that results discussed in this section assume a simple HG function for the SPF and a Rayleigh-like function for the polarization fraction phase function. If alternative forms of the SPF and pSPF are used, the diagnostic relationships F_{sca λ}/F_{pol λ} would need to be recalculated using a disk model.

6.3.8 Averaged scattered flux for the studied debris disk sample

For each disk shown in Fig. 24, the averaged scattered flux ⟨F_{sca λ}⟩ can be derived by dividing the measured scattered and polarized fluxes by the corrections factors f_{sca λ} and f_{pol λ}, respectively, as described in Sect. 6.3.6. To illustrate this approach, we determined the view factors using Fig. 22 and assuming an asymmetry parameter of g = 0.7 (the average value from our disk modeling) and a maximum polarization fraction of p_max = 0.3 for all disks, regardless of the observation wavelength. Figure 26 displays the positions of the disks in the parameter space [⟨F_{sca λ}⟩/F_⋆, f_disk] following this flux correction.

If both the scattered and polarized fluxes for a given disk are measured at the same wavelength, as is the case for the HD 129590 disk (orange and blue markers on the far right in Figs. 24 and 26), and if the adopted values for g and p_max are appropriate for this disk, then the corrected scattered and polarized flux markers should overlap. This overlap indicates a consistent estimate of the disk’s averaged scattered flux ⟨F_{sca λ}⟩. As shown in Fig. 26, a good match using the adopted view factors is achieved for two HST targets with the lowest IR contrast. For the remaining four targets, the corrections yield different values of the average scattered flux, although the discrepancy for HD 129590 remains within the error bars.

For the adopted values for g and p_max, the location of the data points in the diagram lies close to the line ⟨F_{Sca λ}⟩/F_⋆ = L_IR/L_⋆, indicating that comparable fractions of stellar radiation interacting with dust particles are either scattered or absorbed. This is equivalent to the statement that the typical scattering albedo is ω ≈ 0.5. This rough estimate relies on the assumption of a constant albedo across wavelengths, which may not be valid. For large particles, one would expect ω > 0.5 because a value of ω ≈ 0.5 is already expected from diffraction. A scattering contribution from radiation interacting with the particle would then provide a total ω > 0.5. For small particles, such as interstellar dust, one typically expects ω ⪅ 0.5 in the near-IR, along with a relatively low asymmetry parameter, g ⪅ 0.5, but still a high scattering polarization p_max would be possible (Draine 2003b). This type of dust would not be consistent with the collected observational data of debris disks.

As a final remark, we note that adopting higher values of the asymmetry parameter g or lower values of the maximum polarization fraction p_max to derive the view factor for the polarized contrast shifts the data points upward in the diagram presented in Fig. 26, placing them above the line ⟨F_{sca λ}⟩/F_⋆ = L_IR/L_⋆. This would imply a typical scattering albedo greater than 0.5, which is more consistent with the results obtained from our SD modeling (Sect. 6.3.1).

6.4 Observing conditions

Due to the faint and extended nature of debris disks in scattered light, their detection is strongly influenced by both ambient conditions (e.g., coherence time, seeing, wind speed) and observational parameters such as airmass, detector integration time (DIT or exposure time per frame), total exposure time, and, in the case of ADI, the total rotation angle of the sky field during the observation. For successful DI of debris disks with SPHERE instruments, the coherence time should ideally exceed 3–4 ms, while the seeing conditions should remain below 0.8″. If the coherence time is short (<3 ms), the stellar position behind the coronagraph can deviate by more than one pixel from the intended center position. Misalignment in individual frames reduces the S/N in the final science image after frame combination. Similarly, high atmospheric turbulence (seeing >0.8″) or strong variability in the PSF adversely impacts detection.

This effect was notably observed during polarimetric observations of HD 117214 with ZIMPOL, where the degrading seeing conditions over the two-hour observing run significantly reduced the detectability of the debris disk (Engler et al. 2020). When the seeing exceeds 0.95–1″, the stellar PSF often becomes highly variable, making it unlikely to resolve a debris disk, even if the disk is intrinsically bright. However, in ADI observations, a sufficiently large sky rotation angle and a robust number of frames allow for selective frame combination, which can improve the S/N by applying thresholds on seeing (or equivalently PSF FWHM).

The quality of data obtained with SPHERE instruments is also affected by wind speed at Paranal Observatory. The optimal wind speed range for observations lies between 2 and 5 m s⁻¹. When wind speeds drop below 1 m s⁻¹, data quality is impacted by the “low wind effect” (LWE; Milli et al. 2018). To mitigate this issue, a special low-emissivity coating for the telescope spiders was introduced in August 2017. Before this modification, SPHERE observations conducted at wind speeds below 4 m s⁻¹ were often degraded by the LWE. Conversely, highaltitude turbulence with wind speeds exceeding 5 m s⁻¹ produces a wind-driven halo within the AO correction radius, thereby reducing achievable contrast (Cantalloube et al. 2018).

Furthermore, debris disks observed at high airmass (>1.7) are essentially undetectable. The angular size of the disk also plays a crucial role; it must fit within the FOV of the instrument or exceed its inner working angle. This requirement is not met for some targets in our sample, such as HD 3003, where the suspected debris disk around the primary star of this binary system likely has an angular size smaller than 0.1″.

From the analysis of detections and non-detections of debris disks, we conclude that numerous observational requirements must be fulfilled for successful imaging of debris disks in scattered light with SPHERE. Under optimal conditions, debris disks with low IR excess (~10⁻⁴) can be imaged, whereas even bright disks (f_disk > 10⁻³) may remain undetected under suboptimal conditions. These considerations are critical for interpreting debris disk brightness in scattered light and estimating their scattering albedos, particularly when comparing different systems.

7 Debris disks as integral components of stellar systems

Along with stellar, substellar and planetary companions, debris disks represent key components of stellar systems, each tracing different aspects of their formation and dynamical evolution. Bright debris disks around young stars serve as visible signposts of ongoing or past planet formation, and their morphology offer insights into the dynamical environment of a system, often shaped by gravitational interactions with nearby massive bodies. Understanding how debris disks relate spatially and dynamically to both exoplanets and stellar companions is essential for building a more complete picture of young planetary systems.

In this section, we explore these relationships by examining debris disks detected with SPHERE in the context of the companions present around the stars in our sample. Among these systems are 23 young stars with confirmed exoplanets, discovered through various detection techniques, including DI, astrometry, transits, and radial velocity (RV) measurements. Table H.1 lists the parameters of these exoplanets, including only those with masses below 13 M_Jup, the threshold below which companions are categorized as planets. All other types of companions, both confirmed and candidates, are described in Appendix H and listed in Table H.2.

7.1 Companions to the program stars

The debris disks detected with SPHERE occupy a specific region within stellar systems, clearly distinct from the regions where stellar and planetary companions are found. This distinction is likely due to both the intrinsic characteristics of debris disks and selection effects that lead to their detection. To explore these, we considered the separation - mass plane presented in Fig. 27. In this figure, we compared the distribution of debris disks with that of companions to the stars listed in Appendix H.

To interpret the locations of debris disks and compact companions around stars, we first briefly review the main selection effects at play. The FOV of SPHERE allows the detection of disks at projected separations of at least 0.1–0.15″ from the star. Disks located closer are obscured by the coronagraph and are difficult to detect due to the brightness of the stellar halo. Conversely, disks at separations beyond 5.5″ may fall outside the instrument’s FOV, though they may still be detected if observed at a high inclination. Given that the distances to the program stars range from a few to over a hundred parsecs, the number of stars effectively surveyed for debris disks of a given physical size varies. This variation is illustrated by the background grayscale in Fig. 27. We note that the distribution of disk radii is narrower than the region where the survey is complete across the program stars, suggesting that an underlying physical cause, rather than observational bias alone, may be responsible.

Regarding the distribution of stellar companions, very few are found at separations similar to those of the debris disks in this sample. However, we know that the range from a few to hundreds of au is populated by many stellar companions and is the peak of the distribution of stellar companions for solar and A-type stars (Duquennoy & Mayor 1991; Raghavan et al. 2010; De Rosa et al. 2014; Gratton et al. 2023b). The lack of such companions in our sample may reflect the SPHERE observation strategy, as targets with known companions in this separation range were excluded from the SHINE survey, and only shallow observations were performed when such systems were inadvertently included (see Bonavita et al. 2022). However, stellar companions at these separations can destabilize debris disks, so it is reasonable to expect that they are absent in most of the stars with disks considered in this study. In an unbiased search for multiples among stars hosting debris disks, Rodriguez et al. (2015) found that the properties of disks in binary systems are not statistically different from those around single stars. However, their sample contained very few stellar companions with semimajor axes in the 10–100 au range, where the debris disks in our sample are found.

7.2 Architectures of individual planetary systems

The spatial distribution of planetary companions relative to debris disks provides critical insights into the architecture and dynamical history of planetary systems. In particular, the location of planets within or near debris belts can indicate past migration processes, sculpting effects, and zones of dynamical stability or instability. In the sample of stars observed with SPHERE, the planetary companions (M_p < 13 M_Jup) are located closer to their host stars than the resolved debris disks (Fig. 27). This separation is expected from both formation models and observational biases, but it also carries important implications for the interpretation of system dynamics.

Among the systems in our sample, only five (HD 39060, HD 106906, HD 114082, HD 197481, and HD 218396) have both debris disks resolved with SPHERE and confirmed planetary companions. In these cases, the relative radial positions of disks and planets are well constrained, allowing for a detailed investigation of system architecture and the study of how planets shape and interact with circumstellar debris.

Figure 28 presents the configurations of these five systems in comparison with the Solar System, highlighting both similarities and diversity in planet–disk arrangements. This comparative approach helps to identify trends and anomalies that can guide future searches and refine theoretical models. By analyzing these architectures, we can also place constraints on the potential locations and masses of additional, as-yet-undiscovered planets that may reside between known companions and debris structures.

HD 39060 (β Pic). The HD 39060 disk is the most prominent debris disk among all targets, as it was the first to be directly imaged (Smith & Terrile 1984). It is highly asymmetric and extends beyond 1000 au (e.g., Janson et al. 2021). This vast disk consists of two cold exo-Kuiper belts, both observed edgeon but with slightly different PAs (Golimowski et al. 2006; Ahmic et al. 2009, this work). Additionally, several inner planetesimal belts are likely present (Okamoto et al. 2004; Wahhaj et al. 2003), though their exact radial positions and orientations remain uncertain; therefore they are not depicted in Fig. 28.

Two giant planets reside in the inner regions of the disk: HD 39060 b, a ~11 M_Jup planet located at ~9 au (Lagrange et al. 2010), and HD 39060 c, a ~10 M_Jup at ~3 au (Lagrange et al. 2019). The latter is positioned near the radial distance from the star where the stellar incident flux matches the solar incident flux on Earth, I_⊕. This location is indicated by the blue vertical line in the right panel of Fig. 28, which compares the radial positions of exoplanets in terms of the incident flux or irradiance I_⋆ from their host stars. For the two planets, a similar irradiance level suggests that their equilibrium temperatures, T_eq (or BB-equivalent temperatures), are comparable, given that $T_{\text{eq}}\sim I_{\star }^{1/4}$.

HD 218396 (HR 8799). The planetary system of HD 218396 exhibits an even stronger resemblance to the Solar System in terms of stellar irradiance levels experienced by its giant exoplanets compared to the Jovian planets (Fig. 28 right panel). Its architectural structure, featuring a warm belt (R_belt ~ 15 au), a cold belt (R_belt ~ 180 au), and four giant planets residing in the space between them, closely parallels the Solar System’s main asteroid belt, Edgeworth-Kuiper belt, and the orbits of the Jovian planets in between (Marois et al. 2010).

As discussed in Sect. 4, the polarized intensity image taken in the H band with IRDIS reveals the warm belt spatially resolved for the first time (Fig. 9i). Its radial position, indicated in left panel of Fig. 28, is close to that of planet HR 8799 e, which has a projected separation of 0.39″ ± 0.01″ (~16 au) in the same dataset. This spatial coincidence hints at a possible role for the planet in shaping the inner belt’s edge or maintaining its structure through dynamical shepherding.

HD 106906. The HD 106906 is another prominent planetary system hosting a debris disk around a spectroscopic binary, composed of two F5V stars with nearly equal masses (Rodet et al. 2017). The edge-on disk is oriented ~21° away from a planetary-mass companion, which is located at a large projected separation of 650 au from HD 106906 AB (e.g., Bailey et al. 2014; Lagrange et al. 2016; Kalas et al. 2015). The disk appears symmetric in SPHERE near-IR images (Lagrange et al. 2016), in contrast to optical-wavelengths images, which reveal a needle-like structure. A similar morphological difference is observed in the near-IR and optical images of the HD 15115 debris disk (Kalas et al. 2007; Engler et al. 2019) suggesting that the asymmetric, needle-like appearance of edge-on disks is predominantly shaped by submicron-sized dust particles.

HD 114082. HD 114082 is a young F3V star belonging to the Lower Centaurus Crux subgroup of the Sco-Cen association. Its debris disk, with a fractional luminosity of f_disk = 3.8 × 10⁻³, ranks among the brightest disks observed in scattered light to date (Wahhaj et al. 2016; Engler et al. 2023). A transiting super-Jovian planet was detected orbiting the star at a radial distance of 0.5–0.7 au using the RV and transit techniques (Zakhozhay et al. 2022b; Engler et al. 2023). If the predicted orbital parameters of this planet are accurate, its orbit is not co-planar with the midplane of the debris disk, which has an inclination of ~83°.

A second transiting candidate has been detected around HD 114082 in photometric data from the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015). The system was observed by TESS in Sectors 38, 64, and 65, with a second distinct transiting event clearly identified in Sector 64. However, this event does not correspond to HD 114082 b, as the observed transit duration and depth do not match those previously reported by Zakhozhay et al. (2022b). This second single-transit event is not associated with an asteroid crossing or centroid shifts, suggesting that it may be caused by a second planetary body.

To analyze this candidate, we download the TESS light curves produced by the Science Processing Operation Center (Jenkins et al. 2016) from the Mikulski Archive for Space Telescopes. For photometric modeling, we selected the pre-search data conditioning simple aperture flux (PDCSAP) and associated errors. Using the software package juliet (Espinoza et al. 2019) with the dynesty sampling method, we modeled the TESS light curve of Sector 64, incorporating a transiting planet model and a Gaussian process to account for stellar variability. A log-uniform prior ranging from 33 to 1000 days was set for the orbital period of the candidate planet. Additionally, we used the stellar parameters derived in Zakhozhay et al. (2022b) to impose a normal prior on the stellar density, which, when combined with the transit model, constrains the planetary transit speed across the stellar disk.

For this second planetary candidate, we obtain a large radius of ${1.29}_{-0.05}^{+0.05}$ R_Jup and orbital parameters listed in Table 5 along with the orbital parameters of planet HD 114082b for a comparison. The TESS light curve with the median transit model is shown in Fig. 29.

The system HD 114082 offers a compelling example of a planetary system where giant planets orbit close to the host star, while a spatially extended debris disk lies much farther out and appears misaligned with the inner planetary orbits, so that their planes are inclined by at least 6–7 degrees relative to the disk plane.

Such orbital misalignment between planets and debris disks have important implications for the dynamical history of the system. Several mechanisms have been proposed to explain such configurations. One possibility is planet–planet scattering, where gravitational interactions between giant planets lead to the ejection of one or more bodies and a reconfiguration of the survivors onto eccentric and inclined orbits. This process can disrupt the coplanar architecture established during the PPD phase and has been invoked to explain the high eccentricities and inclinations observed in many giant exoplanets (e.g., Chatterjee et al. 2008; Jurić & Tremaine 2008; Raymond et al. 2011).

Another mechanism is secular perturbations from additional, possibly undetected, massive companions. If a distant planet is present on an inclined orbit, it can induce long-term precession in the orbits of inner planets or in the disk itself, leading to gradual misalignment over time. This mechanism, particularly the Kozai-Lidov effect, has been shown to cause significant orbital inclination variations under specific conditions (e.g., Nagasawa et al. 2008; Naoz et al. 2013).

In the case of HD 114082, the wide separation between the planetary system and the debris belt, more than a factor of 30 in radius, creates a dynamical environment in which such longterm perturbations are plausible. Furthermore, the large radial cavity separating the planets from the debris disk may itself be sculpted by additional low-luminosity companions or represent the aftermath of dynamical clearing by now-ejected bodies.

The HD 114082 system offers a compelling case for investigating the early dynamical evolution of tightly packed planetary systems accompanied by outer debris structures. Although TESS has detected single transit events for both planets, their orbital geometries remain poorly constrained due to the limited transit coverage, and the PAs of their orbital planes are still unknown. Continued photometric and spectroscopic monitoring of HD 114082 is essential to refine the orbital parameters of the planets, including their mutual orientation and orbital eccentricities.

In particular, precisely determining the transit timing would enable targeted RV observations during the planet’s passage across the stellar disk. This could allow for the detection of the Rossiter–McLaughlin effect (Triaud 2018), a spectroscopic signature that occurs when a transiting planet temporarily blocks part of the rotating stellar surface. Measuring this effect provides the sky-projected angle between the planetary orbital axis and the stellar spin axis, offering insight into the mutual alignment (or misalignment) between the planetary orbital planes and the host star’s equatorial plane. This, in turn, helps to reconstruct the three-dimensional architecture of the system and test scenarios of planet migration or dynamical perturbation.

HD 197481 (AU Mic). HD 197481 is an active flaring M1Ve star belonging to the βPMG and hosts a highly dynamical debris system. A monitoring campaign of AU Mic conducted with SPHERE between 2015 and 2017 revealed multiple dust clumps appearing above and below the disk midplane (Boccaletti et al. 2018), seemingly moving in non-Keplerian orbits around the star. The origin of these fast-moving dust structures is still debated. One possibility is that they are generated in a collisional avalanche, where the main planetesimal belt intersects a debris stream resulting from the catastrophic disruption of a large asteroid-like body (Chiang & Fung 2017). Another explanation suggests that these dust clumps are expelled from the system by AU Mic’s strong stellar wind, having been emitted by a parent body, such as a planet or a disk substructure within the main debris belt (Sezestre et al. 2017). Indeed, AU Mic is known to host three confirmed planets that follow close orbits, possibly in a 4:6:9 mean-motion resonance (see Table H.1; Plavchan et al. 2020; Martioli et al. 2021; Wittrock et al. 2023). However, they are too close to the star to be responsible for the observed dust clumps (Boccaletti 2023). In Fig. 28, these planets are shown inside two dust belts, positioned at the locations of radial SB peaks as measured along the disk’s major axis in the r²−scaled H-band image from May 20, 2017. The AU Mic’s morphology, however, varies over time and resembles an edge-on spiral structure more than a stable belt system.

The comparative architectures shown in Fig. 28 indicate that our Solar System is not unusual but instead falls within the general distribution of planetary systems hosting both debris disks and planets. In most benchmark systems, the giant planets are confined to the region inside the exo-Kuiper belt, typically spanning orbital radii from a few to a few tens of au, while the planetesimal belts extend from several tens to more than one hundred au. An exception is HD 106906, where a massive (13 M_Jup) companion resides well outside the belt. When scaled by stellar irradiance (right panel in Fig. 28), the exo-Kuiper belt locations cluster around the equivalent solar-system value, underscoring the analogy. This alignment demonstrates that multi-planet arrangements inside wide Kuiper-belt analogs are common, and that the Solar System’s configuration, with its giant planets located interior to an outer cold belt, represents a recurring outcome of planetary system formation and evolution.

7.3 Theoretical predictions for the masses of undetected planets

To investigate the potential planetary architectures in other systems with detected debris disks in our sample, we applied a set of analytical prescriptions under the assumption that observed gaps and edges of planetesimal belts have a dynamical origin. The adopted method depends on whether the system hosts a single resolved belt or multiple belts. For all systems, we assumed a radial belt width equal to 20% of its resolved central radius in consistency with our SED modeling approach (Sect. 5.4).

For systems with only one resolved belt, we estimated the range of planet masses and semimajor axes for a single planet on a circular orbit capable of sculpting the inner edge of the belt. We assumed that such a planet dynamically clears a region interior to the belt by gravitationally ejecting dust particles. To relate the width of the cleared zone, Δa, to the planet’s mass, M_p, we adopted the expression derived by Morrison & Malhotra (2015): $$\mathrm{\Delta }a=1.7{\left(\frac{M_{\mathrm{p}}}{M_{\star }}\right)}^{0.31}{a}_{\mathrm{p}},$$

where a_p is the semimajor axis of the planet and M_⋆ is the mass of the host star. Planet masses were explored in the range of 0.1–13 M_Jup, with corresponding semimajor axes calculated to match the observed belt edge. The resulting minimum and maximum allowed semimajor axes, corresponding to the assumed mass bounds, are reported in Table 6.

For systems with two or more resolved belts, we examined the gap between each pair of belts, testing dynamical configurations involving one, two, or three planets on circular orbits capable of maintaining the observed separations. This analysis follows the framework of Lazzoni et al. (2018). In multi-planet scenarios, we assumed equal-mass planets arranged in a maximally packed configuration, thus with the minimum orbital spacing required for long-term dynamical stability, in order to reduce degeneracy in the possible solutions. For each configuration, the resulting planet masses and semimajor axes are reported in Table 7. For HD 141943, only the single-planet scenario is reported, as the resulting mass required to account for the observed belt structure is already extremely low (5 × 10⁻⁶ M_Jup), making the consideration of multi-planet solutions unnecessary.

The results of our dynamical analysis suggest that many debris disk systems may host planetary companions whose gravitational perturbations shape the observed disk morphology. This is particularly evident in the case of HD 218396 (HR 8799). In a three-planet configuration with equal-mass planets of 6.78 M_Jup, our model yields a configuration that closely approximates the known planetary architecture of this system (Sect. 7.2), albeit without reproducing it in full detail.

For three other systems, HD 36546, HD 92945 and HD 120326, we find that a single planet located in the gap between two resolved belts must possess a mass greater than 3.3 M_Jup, to account for the observed structure. This mass range lies within the detection capabilities of current HCI instruments. In the case of HD 92945, contrast limits achieved in IRDIS and IFS observations on 27 January 2018 were analyzed by Mesa et al. (2021). Using AMES-COND evolutionary models (Allard et al. 2012) to convert IRDIS contrast limits into mass constraints, they derived mass limits between 1 and 2 M_Jup at the gap radial position. This comparison indicates that the observed gap is more likely carved by multiple lower-mass planets, although Marino et al. (2019) and Mesa et al. (2021) found that a single planet with a mass of 0.3–0.6 M_Jup could reproduce the structure, assuming a narrower gap width of 20 au. In contrast, our analysis adopts a gap width of 45.6 au (see Table 7), which would require a more massive perturber. According to our model, the masses of planets in a multi-planet configuration for the HD 92945 system would fall below the SPHERE detection limit of 1 M_Jup.

The SPHERE data for HD 120326 were investigated by Bonnefoy et al. (2017), who discovered the double-belt structure around this star and identified ten candidate companions. Based on their positions in the color–magnitude diagram and comparison with earlier HST/STIS observations (Padgett & Stapelfeldt 2016), all candidates were classified as background objects. By converting the 5σ-detection limits into mass constraints, the authors ruled out the presence of giant planets with masses greater than 2 M_Jup at radial separations between 55 and 100 au, corresponding to the edges of the gap in our model (see Table 7). This excludes our single-planet model but still allows for a two-planet configuration with M_p = 1.49 M_Jup, which remains below the detection threshold of available data.

In all other cases investigated, the predicted range of planet masses capable of sculpting the inner edges of single belts or clearing the gaps between two belts extends to values below the current detection limits of HCI surveys. This implies that a significant fraction of planetary systems could contain sub-Jovian or Neptune-like planets in wide orbits, which remain elusive to existing instruments but play a central role in shaping circumstellar dust distributions.

The inferred masses and orbital distances of these hypothetical planets in our analysis align well with the population of exoplanets identified through RV and transit surveys, especially in the super-Earth to sub-Saturn mass range (Winn & Petigura 2020). These studies have shown that such planets are common around a wide range of stellar types and may often occur in multi-planet configurations. The orbital distances we derive, typically tens to over a hundred au, complement these findings by probing an otherwise underexplored region of radial separations. Our results further support the notion that debris disk morphology can serve as an indirect tracer of planetary companions (e.g., Raymond et al. 2011; Lee & Chiang 2016) and highlight the need for combined approaches incorporating disk modeling, DI, and indirect detection methods, to fully characterize the architectures of planetary systems.

8 Summary

In this work, we performed a demographics study of debris disks around young stars observed with SPHERE at optical and near-IR wavelengths. By analyzing a large sample of 161 targets with IR excess, we addressed morphological, photometric, and polarimetric properties of debris disks across different stellar types and ages. From this sample, compiled from archival GTO and open-time program observations, we resolve 40 disks in scattered light and 36 in polarized light, identifying seven systems with two planetesimal belts and two with three distinct belts (HD 131835 and TWA 7). Newly resolved structures include the disks around HD 36968 and BD-20 951, as well as the inner belts of HR 8799 and HD 36546.

Using SPHERE images, we measured geometrical parameters of detected disks through ellipse fitting, and applied a grid of models for higher-quality data to derive disk radii, aspect ratios, radial density slopes, and scattering asymmetry parameters. This uniform modeling enables a consistent comparison of structural disk properties and reveals systematic dependencies of geometrical parameters on stellar luminosity. The inner slopes of grain density distributions steepen with increasing stellar luminosity, while disk vertical aspect ratios tend to decrease. For most systems, aspect ratios are between 0.02 and 0.06, consistent with expectations for collisionally active belts, though gas-rich disks show unusually small values.

A direct comparison between the radii measured in SPHERE scattered-light images and those derived from ALMA and SMA thermal-emission observations demonstrate a close spatial agreement, with a mean radius ratio of 1.05 ± 0.04. In double-belt systems, the outer belts are typically 1.5–2 times larger than the inner ones.

Almost all resolved planetesimal belts contain cold dust with BB temperatures below 100 K, with HD 172555 being the only exception. A weak correlation between belt radius and stellar luminosity is found, following $R_{\text{belt\hspace{0.17em}}}\propto L_{\star }^{0.11\pm 0.05}$. When dividing belts according to dust temperatures associated with CO and CO₂ freeze-out, this correlation becomes steeper: with α = 0.30 ± 0.08 for CO subsample (T_bb < 35 K) and α = 0.30 ± 0.07 for CO₂ subsample (T_bb >35 K). We also investigated how the locations of debris belts evolve with stellar age in relation to the possible migration of ice lines. For disks in the CO₂ subsample, we find that belt radii increase systematically with stellar age, following $R_{\text{belt}}\propto t_{\text{age}}^{0.37\pm 0.11}$, indicating that debris architectures evolve alongside stellar luminosity growth during the PMS phase.

Complementary SED modeling with MBB and SD approaches allowed us to connect disk morphology with photometric properties across both single- and multi-belt systems estimating their luminosities, dust masses and grain SDs. From this analysis, we find that disk fractional luminosities evolve approximately as $t_{\text{age}}^{-1.18\pm 0.14}$ for A-type stars and $t_{\text{age}}^{-0.81\pm 0.12}$ for F-type stars, supporting collisional evolution as the main driver of long-term decay.

In addition, we examined the scaling of disk dust masses with stellar properties and compared them with relations established for PPD. The dust masses derived with MBB approach, scale super-linearly with stellar mass following $M_{\text{dust}}\propto M_{\star }^{{\alpha }_{\text{mass}}}$ with α_mass = 1.6 ± 1.0 for systems aged 10–50 Myr and α_mass = 1.4 ± 0.9 for older disks, similar to the relations observed in 2–3 Myr old star-forming regions (e.g., Pascucci et al. 2016). This continuity suggests that the initial conditions set during the protoplanetary phase strongly influence debris disk evolution. Furthermore, typical debris disk masses decrease by about three orders of magnitude within the first 50 Myr and by nearly four orders at later ages, consistent with collisional depletion. The M_dust − R_belt relation also follows a power law with index larger than 2, implying that more extended belts tend to host proportionally larger dust reservoirs. Together, these trends demonstrate that debris disks preserve imprints of their primordial PPD phase while revealing the efficiency of collisional evolution in regulating dust content over time.

To model the SEDs of the detected planetesimal belts with a SD approach, we adopted belt radii measured from SPHERE scattered-light images. The fits yield minimum grain sizes consistently larger than 0.8 μm and an average power-law slope of q = 3.62, slightly steeper than the canonical collisional cascade value of 3.5. The resulting dust masses of exo-Kuiper belts, integrated over particle sizes from a_min to 5 mm and assuming astrosilicate composition, lie in the range 0.01–1 M_⊕ and agree with those inferred from MBB modeling. Moreover, these dust masses scale approx. as $R_{\text{belt\hspace{0.17em}}}^{2.1}$ with radial distance in the sub-sample of A-, F-, and G-type stars with ages between 10 and 200 Myr.

Building on the SD modeling results, we estimated bulk dust albedo values using Mie theory for four different grain compositions. The derived values are consistently higher than 0.5, but the variation between compositions is relatively small, indicating that dust albedo is unlikely to be the primary factor behind disk non-detections. In addition, we introduced in this work a parametric approach based on image modeling and flux measurements from scattered-light and polarized-light images, demonstrating how both the dust albedo and the maximum polarization fraction can be derived with this method.

Analysis of polarized-light images reveals a correlation between polarized fluxes in H and J bands and modeled IR excesses. The slope of this relation is shallower than that found for total-intensity optical images from HST, consistent with the fact that polarized flux traces only a fraction of total scattered light, dependent on dust properties and disk inclination.

The analysis of non-detections further shows that 90% of the detected disks have estimated ages below 50 Myr, indicating that many undetected systems can be explained by intrinsically low dust masses resulting in faint scattered-light emission. Nevertheless, even bright disks (f_disk > 10⁻³) may remain undetected in cases of unfavorable viewing geometry or suboptimal observing conditions, whereas under optimal conditions disks with excesses as low as 10⁻⁴ are detectable. These findings establish practical detection thresholds and observational biases that are critical for interpreting the demographics of debris disks in current and future HCI surveys.

The connection between debris disks and planetary architectures was probed through dynamical modeling. In systems with resolved disks, we estimated the masses and semimajor axes of planets that could sculpt gaps or belt inner edges. For single-belt systems, planets capable of shaping the inner edges span masses from 0.1 to 13 M_Jup, with orbital radii consistent with observed disk edges. In double-belt systems, multiple planets with sub-Jovian masses can reproduce the gaps while remaining below current detection thresholds. These findings imply that Neptune- to sub-Saturn-mass planets at tens to hundreds of au may be common but remain undetected. The inferred planet populations align with those found by RV and transit surveys in the super-Earth to sub-Saturn range, though at larger orbital separations.

Among the systems in our sample, only five, HD 39060 (β Pic), HD 106906, HD 114082, HD 197481 (AU Mic), and HD 218396 (HR 8799), host both resolved debris disks and confirmed planets. In HD 114082, we identify a second transiting giant planet candidate in TESS data, with a radius of 1.29 ± 0.05 R_Jup and orbit near 1 au. Its misaligned orbit relative to the debris disk and the known planet, along with their close proximity to the host star, well inside the typical formation region for giant planets, suggests a history of dynamical evolution involving planet–planet scattering or planet migration. Continued monitoring of HD 114082 is required to refine the orbital parameters of both planets, particularly their eccentricities and the orientation of their orbital planes. These constraints are critical for assessing the system’s long-term dynamical stability and for distinguishing between potential migration and scattering scenarios.

Finally, we considered stellar companions. Few are found at separations overlapping the debris disks in our sample, consistent with selection biases excluding known multiples from SPHERE surveys. Nonetheless, stellar companions at 10–100 au are known to destabilize debris disks, likely explaining their scarcity in our resolved sample.

This study presents the largest homogeneous analysis of debris disks imaged in scattered light with SPHERE, complementing previous surveys conducted with HST, Herschel, GPI, and ALMA. It is important to note, however, that our results are based on a sample that is inherently biased, as it primarily includes targets selected for their high IR excess or previously known disk structures. This selection effect should be carefully considered when interpreting trends in disk properties and dust evolution. A more comprehensive and unbiased statistical assessment would require a broader sample, including fainter and lower-mass disks that remain undetected with current HCI techniques.

Future studies will benefit from advancements in observational capabilities, such as those provided by JWST, which can probe debris disks at mid-IR wavelengths (e.g., Boccaletti et al. 2024; Mâlin et al. 2024; Su et al. 2024), as well as next-generation ground-based instruments including ELT/METIS (Brandl et al. 2024). Additionally, combining high-resolution scattered-light imaging with ALMA millimeter observations will allow for a more complete picture of disk morphology, grain composition, and spatial segregation of dust populations. Expanding such multiwavelength approaches, together with advancements in disk modeling techniques, is essential for refining our understanding of debris disk evolution and the architectures of planetary systems.

Acknowledgements

We would like to thank the anonymous referee for many thoughtful comments which helped to improve this paper. SPHERE is an instrument designed and built by a consortium consisting of IPAG (Grenoble, France), MPIA (Heidelberg, Germany), LAM (Marseille, France), LESIA (Paris, France), Laboratoire Lagrange (Nice, France), INAF Osservatorio di Padova (Italy), Observatoire de Genève (Switzerland), ETH Zurich (Switzerland), NOVA (Netherlands), ONERA (France) and ASTRON (Netherlands) in collaboration with ESO. SPHERE was funded by ESO, with additional contributions from CNRS (France), MPIA (Germany), INAF (Italy), FINES (Switzerland) and NOVA (Netherlands). SPHERE also received funding from the European Commission Sixth and Seventh Framework Programmes as part of the Optical Infrared Coordination Network for Astronomy (OPTICON) under grant number RII3-Ct-2004-001566 for FP6 (2004-2008), grant number 226604 for FP7 (2009-2012) and grant number 312430 for FP7 (2013-2016). This work has made use of 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 work has been supported by the DDISK ANR contract number ANR-21-CE31-0015. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This paper includes data collected with the TESS mission, obtained from the MAST data archive at the Space Telescope Science Institute (STScI). Funding for the TESS mission is provided by the NASA Explorer Program. STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5–26555. JM and JCA acknowledge funding from the "Programme National de Planétologie" (PNP) of CNRS-INSU in France through the EPOPEE project (Etude des POussières Planétaires Et Exoplanétaires). A.Z. acknowledges support from ANID – Millennium Science Initiative Program – Center Code NCN2024_001 and Fondecyt Regular grant number 1250249.

References

---

Appendix A ESO program IDs

This study is based on data collected at the European Southern Observatory in Chile under programs 095.C-0192, 095.C-0273, 096.C-0388, 097.C-0394, 096.C-0640, 097.C-0344, 098.C-0505, 098.C-0686, 098.C-0790, 099.C-0147, 0100.C-0548, 0101.C-0016, 0101.C-0015, 0101.C-0128, 0101.C-0420, 0101.C-0422, 0101.C-0502, 0102.C-0916, 0102.C-0861, 0104.C-0436, 0104.C-0456, 105.20GP.001, 1104.C-0416, 0102.C-0453, 1100.C-0481, and 198.C-0209.

Appendix B SPHERE images of disks TWA 7 and HD 16743

To complement the findings discussed in Sect. 4.4, Figs B.1 and B.2 present supplementary SPHERE/IRDIS images of the TWA 7 and HD 16743 systems, highlighting additional structural features within their respective debris disks. Figure B.1 shows the polarized intensity image of the TWA 7 disk, combining H-band data from three separate observing epochs. In this composite image, all three planetesimal belts are discernible, with the outer belt particularly prominent. Also visible are arc-like features connecting the middle and outer belts, labeled “F1” and “F2” in Figs. 9g and B.1. The most extended of these structures, “F2”, reaches the outer edge of the SPHERE FOV and was previously detected in both HST and ALMA observations (Ren et al. 2021; Bayo et al. 2019).

Figure B.2 displays the total intensity image of the HD 16743 disk, showing a planetesimal belt oriented at a PA of ~170°, as well as an extended feature at PA = 17° marked with a question mark. The origin of this feature is uncertain; it may represent either a residual PSF artifact or genuine scattered light from disk material.

Appendix C Polarization fraction function of micron-sized dust particles

To model the disk images of polarized scattered light, we adopted the pSPF incorporating a polarization fraction function derived under the Rayleigh scattering assumption (Eq 3). Rayleigh scattering describes the interaction of light with particles significantly smaller than the wavelength of the incident radiation, and it exhibits a pronounced angular dependence, with maximum polarization occurring at scattering angle of 90°. This approximation provides a useful analytic form for the polarization behavior of small and compact grains in debris disks, although it becomes less accurate for larger or more complex dust particles.

To illustrate this, we present in Fig. C.1 the polarization fraction functions measured for Mg-rich olivine samples with three different particle SDs alongside the theoretical Rayleigh scattering polarization fraction function for comparison. We selected olivine as a representative material because its presence is commonly inferred from cometary and debris disk spectra (e.g., Kolokolova & Jockers 1997; Chen et al. 2006). The measurements for samples with effective grain radii of 2.6 μm and 3.8 μm were conducted using a laser source at 633 nm (Muñoz et al. 2000), while the data for the sample with an effective radius exceeding 20 μm were obtained using a white light source and the spectral response of imaging polarimeter cameras in the wavelength range of 1.5 – 1.6 μm (Renard et al. 2014)¹⁵.

The scaled Rayleigh scattering function p(θ) is plotted in Fig. C.1 with a maximum polarization fraction p_max of 9%, chosen to match the measurement for the sample with an effective grain radius of a_eff = 3.8 μm. As shown, the overall shape of the Rayleigh polarization function reproduces the observed trend reasonably well. The largest discrepancy arises at large scattering angles (θ > 150°), corresponding to small phase angles (phase angle = 180° – scattering angle), where the measured polarization fraction (<4%) becomes negative, forming a so-called negative polarization branch. This sign inversion indicates a change in the orientation of the polarization vector from azimuthal to radial in the image of polarized intensity, meaning that the scattered light becomes preferentially polarized parallel to the scattering plane¹⁶ rather than perpendicular to it.

The maximum polarization fraction of dust particles may also occur at scattering angles different from 90°, though still typically close to it. For example, in Fig. C.1, the sample with an effective grain size of a_eff = 2.6 μm exhibits a peak polarization fraction at a scattering angle of approx. 100°. In practice, the polarization fraction function of a debris disk, along with the SPF and pSPF, is expected to be smoothed due to averaging over a distribution of particle sizes and a mixture of dust populations across different disk regions. We therefore conclude that the Rayleigh scattering function p(θ) offers a reasonable, though simplified, approximation that can be applied to model polarimetric images of debris disks. If a more accurate representation is required, the location of the maximum polarization fraction can be parametrized using a scaled version of the beta distribution, as proposed by Ren et al. (2023). However, this approach introduces two additional degrees of freedom into the model.

Appendix D Images of the total intensity of scattered light

Figure D.1 presents images of the total intensity of scattered light from debris disks detected with IRDIS, IFS, or ZIMPOL.

Appendix E Examples of the SED fits using MBB and SD models

Figure E.1 shows several examples of the SED best-fitting models applying various approaches. Top row presents the MBB models consisting of one or two planetesimal rings. Bottom row shows the SD models based on the belt radii measured from the r²−scaled disk images.

Appendix F Images of the polarized intensity of scattered light

Figure F.1 presents images of the polarized intensity of scattered light from debris disks detected with IRDIS, IFS, or ZIMPOL.

Appendix G Sample stars and results from disk and MBB modeling

Table G.2 lists the sample stars along with their stellar parameters. Tables G.1, G.3, and G.4 present the results of disk modeling and the one- and two-component MBB modeling, respectively.

$F_{\lambda }^{th}$	spectral flux density of thermal emission
fdisk	disk fractional luminosity
$f_{\text{disk}}^{SD}$	disk fractional luminosity obtained with SD model
fpol λ	view factor for polarized flux
fsca λ	view factor for scattered flux
G	gravitational constant
g	HG asymmetry parameter
H	disk scale height
H0	disk scale height at reference radius
h	disk height coordinate
I⋆ λ	intensity of the incident radiation
I⊕	intensity of the solar radiation on Earth
J	number of model free parameters
L⋆	stellar luminosity
L⋆ λ	spectral stellar luminosity
L⊙	solar luminosity
LIR disk	disk IR luminosity
Lsca λ	spectral disk scattered luminosity
M⋆	stellar mass
M⊙	solar mass
M⊕	Earth mass
MB	mass of stellar companion B
MCO	carbon monoxide mass
Mdust	dust mass obtained with MBB model
$M_{\text{dust}}^{\mathrm{S}\mathrm{D}}$	dust mass obtained with SD model
Minner	dust mass of inner belt in SD model
Mouter	dust mass of outer belt in SD model
Mp	planet mass
Mtotal	total dust mass of inner and outer belts in SD model
Ndata	number of data points
Nobs	number of photometric points for wavelengths >22 μm
NSED	surface grain number density in the SD model
n(a)	differential grain number density
ngr	volume grain number density in the disk
P	orbital period
pmax	maximum polarization fraction of scattered light
Q	Stokes parameter Q
⟨Qpr⟩	mean radiation pressure coupling coefficient averaged over the stellar flux
Qr	radial Stokes parameter Q
$Q_{\lambda }^{\mathrm{a}\mathrm{b}\mathrm{s}}$	spectral absorption efficiency
$Q_{\lambda }^{\text{ext}}$	spectral extinction efficiency
$Q_{\lambda }^{\text{sca}}$	spectral scattering efficiency
Qφ	azimuthal Stokes parameter Q
q	SD index
R⋆	stellar radius
Rbelt	planetesimal belt radius
$R_{\text{belt}}^{\text{mes}}$	belt radius measured from the r2–scaled SPHERE images
$R_{\text{belt)}}^{\text{mod}}$	modeled peak volume density of grains
Rinner	radius of inner belt in SD model
RL⊙	scaling factor representing the expected radial position of a belt around a star with solar luminosity
RMBB	disk radius obtained with MBB model
$R_{\text{max}(\sigma )}^{\text{mod}}$	radial distance of the modeled peak surface density
Router	radius of outer belt in SD model
Rp	planet radius
r	radial distance from the star
r0	reference radius of planetesimal belt
T⋆	stellar temperature
Tbb	BB temperature of dust grains
$T_{\text{dust}}^{\text{SD}}$	dust temperature obtained with SD model
Tgrain	grain temperature
Teff	effective stellar temperature
Teq	equilibrium temperature
TMBB	BB temperature of dust grains obtained with MBB model
Ttr	transit duration
tage	stellar age
U	Stokes parameter U
Uφ	azimuthal Stokes parameter U
Ur	radial Stokes parameter U
u	Heaviside step function
x	size parameter
xeff	effective size parameter
α	index of power-law function for the relation between belt radius and stellar luminosity
αin	index of power-law function for the grain number density in the inner region of the disk
αmass	index of power-law function for the relation between belt dust mass and stellar mass
αout	index of power-law function for the grain number density in the outer region of the disk

Appendix H Companions to the program stars

Table H.1 presents the parameters of confirmed exoplanets with masses below 13 M_Jup in the stellar systems of our sample, as of August 1, 2024. The data were retrieved from the NASA Exoplanet archive¹⁷.

To identify companions to the program stars with masses greater than 13 M_Jup, we followed the methodology outlined in Gratton et al. (2023b,a, 2024, 2025). Our search incorporated both direct observations of visual companions, primarily from Gaia data and HCI, and indirect evidence based on photometry (eclipsing binaries), RV measurements (spectroscopic binaries), and astrometry (mainly from Gaia).

H.1 Visual binaries

Wide visual companions (separation >0.7″) can be detected either as separate entries in Gaia DR3 or by HCI. Both are available for all the program stars. Given the age of the stars, Gaia DR3 can reveal companions down to the star-brown dwarfs for separation larger than a few hundred au. HCI is sensitive to massive planetary companions, but only in the separation range from a few tens to about 1000 au.

We considered as physical companions all point sources that have similar parallax and proper motion, and are within 60″ (that is, about 1000-10000 au, depending on distance of the target) of each of the program stars. We also notice that the astrometric solution may be missing for faint sources very close to brighter objects in the Gaia DR3 catalogue. Given the low-density of the background fields, we additionally considered as physical companions the objects listed in Gaia and projected within 2" of each star, but lacking an astrometric solution.

H.2 Eclipsing binaries

We searched for eclipsing binaries (EBs) within the TESS dataset but did not identify any among the program stars.

H.3 Spectroscopic binaries

We inspected the S_B⁹ (Pourbaix et al. (2004), Tokovinin (2018)) and Gaia DR3 binary Gaia Collaboration (2023a) catalogs for spectroscopic and spectrophotometric binaries. In addition, we considered high-precision RV series from the literature and the low-precision one from Gaia DR3. Companions detected using RVs typically have separation less than a few au’s.

H.4 Astrometric binaries

We inspected various catalogs looking for astrometric binaries based on Gaia data: nss_two_body_orbit (Gaia Collaboration 2023a), nss_acceleration_astro (Gaia Collaboration 2023a, Holl et al. 2023).

We considered proper motion acceleration (PMa) from Kervella et al. (2022). The PMa is the difference between the proper motion in Gaia DR3 (baseline of 34 months) and that determined using the position at Hipparcos (1991.25) and Gaia EDR3 (2016.0) epochs. This quantity is available for a vast majority of the program stars. PMa is sensitive to binaries with a projected separation between 1 and 100 au. We consider any value of PMa with a S/N > 3 as an indication of the presence of companions.

We also considered the re-normalised unit weight error (RUWE) as an indication of binarity. This parameter is an indication of the goodness of the 5-parameter solution found by Gaia (Lindegren et al. 2018). Belokurov et al. (2020) showed that when this parameter is > 1.4 the star is likely a binary, at least for stars that are not too bright (G > 4) and saturated in the Gaia scans. This method is sensitive to systems that have periods from a few months to a decade (Penoyre et al. 2022). The RUWE parameter is available for the vast majority of the program stars.

H.5 Parameters for the components

The semimajor axis and the mass for the companions are listed in Table H.2. They are obtained following the methods considered in Gratton et al. (2023b, 2024, 2025), briefly summarised in the following.

For unresolved systems, the sum of the masses is made compatible with the apparent G magnitude of the system, using the mass-luminosity relation for the Gaia G band appropriate for the age of the system. We assumed that the semimajor axis is equal to the projected separation divided by the parallax. On average, this corresponds to the thermal eccentricity distribution considered by Ambartsumian (1937) of f(e) = 2e (see Brandeker et al. 2006). Uncertainties in the masses derived using these recipes are small (well below 10%), while those for the semimajor axes are about 40% (see Figure A.1 in Brandeker et al. 2006).

The indication of binarity for many objects comes from RUWE (>1.4) or PMa (S/N > 3 objects) or a combination of these techniques. The secondary of these stars is not imaged, and no period or semimajor axis is determined. Since RV variations, RUWE, and PMa have different dependence on the semimajor axis and the mass, we may better constrain the parameters of the companions by combining different methods using all available information rather than considering only a single technique. We did this by means of exploring the semimajor axis – mass ratio plane using a Monte Carlo code. We assumed eccentric orbits, with uniform priors between 0 and 1 in eccentricity (which is in agreement with Hwang et al. 2022 for this range of separations), 0 and 180 degrees in the angle of the ascending node Ω, and 0 and 360 degrees in the periastron angle ω, and left the inclination and phase to assume a random value. In addition, the period was used to fix the solution whenever it was available. Uncertainties are large and only give order-of-magnitude estimates.

αR	index of power-law function for the relation between belt dust mass and MBB radius
αt	index of power-law function for the evolution of disk fractional IR luminosity
β	disk flaring index
βmass	scaling factor of power-law function for the relation between belt dust mass and stellar mass
βop	dust spectral opacity index
Δa	width of the cleared zone
ΔRbelt	belt width in SD model
Δλ	filter wavelength range
θ	scattering angle
i	disk or planet orbit inclination
Λpol λ	Λ parameter for polarized flux
Λsca λ	Λ parameter for scattered flux
⟨Λsca λ⟩	Λ parameter for disk-averaged scattered flux
λ	wavelength
λ0	SED characteristic wavelength
λc	central wavelength of filter
ϱ	grain material density
Σ	surface density in the SD model
${\sigma }_{\lambda }^{\text{ext}}$	spectral total cross section for extinction
${\sigma }_{\lambda }^{\text{cca}}$	spectral total cross section for scattering
ν	model degrees of freedom
Ω	solid angle
ωλ	spectral single scattering albedo

Appendix I List of symbols

In this section, we provide a list of the symbols employed throughout this work.

A	geometrical cross section of particle
Adisk	vertical disk aspect ratio
a	particle radius
ablow	blowout grain size
amax	maximum grain radius in SD model
amin	minimum grain radius in SD model
ap	semimajor axis of planet orbit
ap max	maximum semimajor axis of planet orbit
ap min	minimum semimajor axis of planet orbit
Bλ	Planck function
b	impact parameter
$C_{\lambda }^{\mathrm{a}\mathrm{b}\mathrm{s}}$	particle spectral absorption cross section
$C_{\lambda }^{\text{ext\hspace{0.17em}}}$	particle spectral extinction cross section
$C_{\lambda }^{\mathrm{s}\mathrm{c}\mathrm{a}}$	particle spectral scattering cross section
c	speed of light
d	distance between the Earth and a star
F⋆ λ	spectral stellar flux
Fgr	gravitational force acting on particle
Frad	radiation pressure force acting on particle
$F_{\lambda }^{\mathrm{a}\mathrm{b}\mathrm{s}}$	spectral absorbed power
$F_{\lambda }^{\mathrm{s}\mathrm{c}\mathrm{a}}$	spectral scattered power
Fpol λ	polarized flux measured from disk image
Fsca λ	scattered flux measured from disk image
⟨Fsca λ⟩	scattered disk flux averaged over the full solid angle

All Tables

All Figures