© The Authors 2026

1 Introduction

Recent observations provide significant insight into the early stages of the planetary formation process, especially through the direct detection of planets in protoplanetary disks, as illustrated by the direct imaging of PDS 70 (Keppler et al. 2018; Blakely et al. 2025). These discoveries further confirm that protoplanetary disks, rich in gas and dust, are the sites of ongoing planetary formation (Pinte et al. 2018; Teague et al. 2018). However, the growth of sub-micrometric grains into kilometer-sized planetesimals over a very short timescale remains a mystery in the process of planet formation (Drążkowska et al. 2023). Several studies have investigated the physical conditions and mechanisms governing this process, including the dynamical aspects of dust vertical settling (Dubrulle et al. 1995; Villenave 2025), radial drift (Weidenschilling 1977), streaming instabilities (Youdin & Goodman 2005; Johansen & Youdin 2007; Villenave 2025), and the influence of dust aerodynamic behavior and coagulation efficiency (Ormel & Cuzzi 2007).

Today, we are in a position to better understand the physical properties of dust grains evolving in protoplanetary disks, such as their size, shape, and composition, to carry out these investigations and move toward a deeper understanding of the dust growth process. In this context, the number of planetforming disks observed in scattered light has increased substantially over the past decade, driven by large-scale observational programs such as GTO (Beuzit et al. 2019) and DESTINYS (Ginski et al. 2021). These VLT/SPHERE observations revealed a wide diversity of substructures in the disks, including rings, gaps, and spirals (Benisty et al. 2023), possibly indicating close interactions between the disks and the forming planets. These near-infrared scattered-light observations are strongly dependent on the properties of dust lying on the disk’s surface layer (Mulders et al. 2013). It is therefore possible to analyze the ability of the disk’s dust to scatter the stellar radiation. This measurement of the polarized or total brightness is the scattering phase function (SPF; Bohren & Huffman 1983) and is defined by the properties of dust grains (Min et al. 2016; Tazaki et al. 2016). The SPF is extracted by measuring the brightness coming from different zones on the disk’s surface as a function of the scattering angle. The extraction of the SPF has been frequently performed on debris disks (Olofsson et al. 2020; Milli et al. 2019; Monnier et al. 2017), but its application to protoplanetary disks is relatively new due to their optically thick nature and their more complex geometry caused by the surface flaring. Research on young planet-forming disks started with specific cases such as IM Lup (Tazaki et al. 2023) to estimate dust properties. Meanwhile, Ginski et al. (2023) analyzed SPFs of polarized light on several disks to study their shape and try to identify general trends.

It is now established that dust grains are not perfectly spherical and that the Mie theory does not provide the most suitable method for describing their optical properties. Recent work has shown that the grains observed in comets such as 67P/Churyumov-Gerasimenko are fractal aggregates (Blum et al. 2017, 2022; Mannel et al. 2019). These results are in line with numerical simulations (Lombart & Laibe 2021; Lombart et al. 2022) and laboratory experiments (Blum 2018), in which dust growth by coagulation and fragmentation produces complex aggregates. In this context, the development of the AggScatVIR¹ database (Tazaki & Dominik 2022; Tazaki et al. 2023) represents a significant step forward, providing SPFs for irregular grains and aggregate particles - be they compact or fractal - composed of sub-micrometric monomers.

The complex geometry of protoplanetary disks has long limited statistical studies of SPFs. The new DRAGyS² tool (Roumesy et al. 2025) overcomes this barrier by deriving disk geometry from ring structures and extracting the SPF directly from the disk’s surface. This approach generalizes previous individual geometric analyses (de Boer et al. 2021; Avenhaus et al. 2018; Ginski et al. 2016) and removes the need for time-consuming radiative-transfer simulations.

In Sect. 2, we present the 30 ring-shaped disk-polarized-intensity Qφ images observed with the VLT/SPHERE instrument (Beuzit et al. 2019). Section 3 describes the methods and main results. We first estimate the disk geometry, including flaring in Sect. 3.1, then we describe the SPF extraction method using DRAGyS in Sect. 3.2. Finally, we analyze the dust properties with AggScatVIR in Sect. 3.3 to identify general trends in the properties of dust within protoplanetary disks. A discussion of these results is presented in Sect. 4.

2 Observations

The objects selected in this study are moderately inclined protoplanetary disks with one or multiple well-resolved rings in polarized scattered light in the near-infrared bands. Our dataset combines observations from the SPHERE/ESO archives with more recent observations from the Disk Evolution Study Through Imaging of Nearby Stars (DESTINYS) ESO large program (Ginski et al. 2020, 2021), mostly located in the Orion (Valegård et al. 2024), Chamaeleon I (Ginski et al. 2024), and Taurus (Garufi et al. 2024) star-forming regions. All of these data were obtained using the IRDIS (Dohlen et al. 2008) instrument of VLT/SPHERE (Beuzit et al. 2019) in dual-beam polarization mode to obtain Q_φ polarized images with a coronagraph to improve the contrast. The Q_φ images of selected objects are shown in Fig. 1.

As an additional source, the disks around TYC 6213-1122-1 and Haro 1-1, both coming from the Scorpius star-forming region, were observed in polarimetry for the first time during the DESTINYS ESO large program (PI:Ginski). We provide an estimate of inclination and position angle for these two systems based on their ring structure using the DRAGyS tool (Roumesy et al. 2025). Haro 1-1 lies at 151, 1 pc (Gaia Collaboration 2020).

The disk has an inclination of 54°, a position angle of 225.0° ± 3.8°, and seems to have an inner ring that is partially hidden by the coronagraph. Another ring is detectable, centered around 50 au from the central star and reaching out to ~80 au. TYC 6213-1122-1 lies at 154.6 pc (Gaia Collaboration 2020). Its inclination is around 57°, with a position angle of 31.8° ± 4.1°. It reveals a relatively broad ring extending from about 30 au to over 100 au. A bright annulus stands out in the inner part of this ring, at around 40 au from the central star (see Table A.1).

3 Methods and results

In this study, we aim to estimate the properties of dust grains in protoplanetary disks by analyzing the SPF. Extracting this measurement requires knowledge of the disk geometry-inclination, position angle, scattering surface height, and flaring-to define the scattering angle, i.e., the angle describing the change of propagation direction between the incident and scattered-light beams (0° meaning no change in direction or forward-scattering); and to constrain the analysis within a specific ring of the disk. The observations presented here are scattered-light data, so the flaring refers to the last scattering surface and not to the gas scale height. This study applies the DRAGyS tool (Roumesy et al. 2025), which enables a comprehensive and consistent scattered-light analysis, from the geometry to the SPF extraction. This section describes the methods used and the results obtained for the disk geometry and the SPFs and identifies the best dust model for each one.

3.1 Disk geometry

Determining the disk geometry is based on the identification of elliptical shapes corresponding to rings within disks to estimate the inclination, position angle, and scattering-height surface, assuming circular axisymmetric rings. This is respectively performed by measuring the minor and major axes of the fitted ellipse, the orientation of the disk, and the offset between the ellipse center and the position of the central star (see Roumesy et al. (2025) for a detailed description of the DRAGyS geometric fitting process). We provide an estimation of the geometry for each disk observation, based on this ring detection in Table A.1.

For disks observed at multiple wavelengths, the surface of last scattering may differ from one band to another, although the observation bands are relatively close to each other. For example, slight deviations of a few degrees are observed in the inclination estimation for the disks around V4046 Sgr, PDS 66, and HD 34282, depending on the wavelength of observation, but they remain within the confidence range. However, there is a significant deviation in the position-angle estimation for UX Tau, which is measured at 339.3° ± 1.8° in the J band and 355.0° ± 2.8° in the H band. This offset of more than 15° is unexpected considering the proximity of the wavelength, but UX Tau is a specific case as it is affected by a perturber (Ménard et al. 2020). To ensure consistency between the geometry and the SPFs, we decided to estimate the disk geometry and extract the corresponding SPFs separately for each band.

We estimated the scattering height at a given radius for a ring for each disk, but it is not possible to estimate the disk flaring with only one measurement of height. However, our sample includes several multiple-ring systems, such as HD 97048, HD 34282, RX J1615.3-3255, and SZ Cha. We assumed that the scattering surface height as a function of radius is governed by a power law, $$h_s(r)=h_0\times (r/r_0{)}^{\alpha },$$(1)

where h₀ is the reference height at a reference radius of r₀, and α is the scattering surface flaring exponent. Therefore, we estimated the scattering height surface for each ring of multiple ring systems individually, and we adjusted this h_s(r) function, fixing the reference radius, r₀, at the mean radius of each target (see Fig. 2). For the rest of our study, we used the scattering scale height and flaring estimations using this power-law fitting for these multiple-ring systems.

An interesting remark is that the measurement of the scattering height in different bands is very similar, reinforcing the accuracy with which DRAGyS determines the disk’s surface height. This is the case for the disk around HD 34282, with h_s/r ≃ 0.28 and α ≃ 1.22 in the J and H bands, α ≃ 1.17 in the K band, and around RX J1615.3-3255 with four rings, with h_s/r ≃ 0.16 and α ≃ 1.3. In contrast, the single-ringed disk around PDS 66 exhibits a decreasing scattering aspect ratio (h_s/r) with increasing wavelength, which may indicate that different vertical layers of the disk are being probed. However, given the image noise (particularly in the K-band observations) and the absence of a similar trend in other disks observed across multiple wavelengths, it remains difficult to draw a firm conclusion about the layering of the scattering surface. In the case of the disk around HD 97048, previous studies reported a surprisingly high flaring of 1.73 (Ginski et al. 2016) and suggested that the flux from the disk’s back side was mainly dominated by the inner rim of the disk rather than its surface. However, the new observation provides better contrast and reveals a greater part of the disk’s back side, leading to a more robust estimation of the disk’s geometry, especially to a lower flaring exponent of 1.23 (see Fig. 3), which is in good agreement with the previous estimation by Lagage et al. (2006).

In order to have a global view of the scattering height surface tendencies of protoplanetary disks, we considered all the scattering surface height measurements for all disks simultaneously to adjust the power-law function defined by Eq. (1), setting the reference radius at r₀ = 100 au. We estimate an average scattering scale height of h₀/r₀ = 0.176 ± 0.011 and an average flaring exponent of α = 1.36 ± 0.11 (see Fig. 4), consistent with a previous estimation of average scattering scale height (h₀/r₀ = 0.1617), but with a stronger disk flaring (α = 1.219 in Avenhaus et al. 2018). Our sample is larger than in the previous study, resulting in greater dispersion in the surface heights measured for protoplanetary disks. This highlights the wide diversity of disk geometry, which is expected given that the disk’s vertical structure depends on many physical parameters (e.g., stellar mass, luminosity, and dust mass) that differ significantly from one system to another. For the rest of our study, we used the individual flaring estimate for disks with multiple rings, and we took this global estimated flaring exponent of α = 1.357 for the others.

The overall disk-geometry estimate is sufficiently precise to ensure that the SPF extraction zone is properly positioned on the ring surface. As the scattering surface height estimation was performed directly on the target ring, the influence of flaring on the extraction zone remains minor, especially in the case of narrow rings where scattering height variation is smaller. However, for wider structures, this approximation becomes less reliable, such as in the case of the PDS 70 protoplanetary disk. With a strongly flared ring, measuring the scattering surface height on the inner part of the ring with DRAGyS is no longer applicable to the entire ring surface. Therefore, to better encompass the ring surface within the extraction zone in the case of PDS 70, we manually increase the flaring value to 1.6 and the scattering scale height to 0.175 rather than the estimated value of 0.135. We can note the same feature for TYC 6213-1122-1. It also has a relatively wide ring, but with a very bright region in the inner part.

For this reason, we decided to focus the analysis on this bright region rather than the entire ring.

As an additional point, the disk flaring estimation based on the estimated scattering height for each ring does not necessarily follow the theoretical flaring expected for a smooth, continuous disk. Indeed, this flaring is well established for pro-toplanetary disks without substructures (Chiang & Goldreich 1997; Dullemond et al. 2001), but the formation of rings and gaps can significantly impact the vertical structure of the disk, particularly through variations of dust mass within each ring, which affects the scattering height from ring to ring. Consequently, the optical-depth surface, τ = 1, observed on each ring does not necessarily reproduce the global flaring of the disk. All of our flaring estimations are greater than 1, which is consistent with the geometric requirements for the observation of such rings. Curiously, the flaring values derived from the scattering surface height of each ring, as in our study, are consistent with the case of a continuous disk-flaring exponent, suggesting an underlying similarity between the surface properties of disks with and without marked substructures.

Fig. 1 VLT/SPHERE Qφ observations in the J-, H-, and K-bands. All images have dimensions of 3.4″ × 3.4" and are shown using an individual log-scale color map.

Fig. 2 Scattering-surface-height fitting for multiple ring-shaped disks. We plot adjustments as dashed colored curves, and we display each estimated value of scattering height ratio, h0/r0, and flaring exponent in Table 1.

Fig. 3 Ring fitting on HD 97048. We compare the estimation presented in Ginski et al. (2016) (in orange) and our estimation (in cyan). The image is on R2-scaled, and we only plot half of the ellipse to bring up substructures. Features to the northwest are artifacts from data reduction.

Fig. 4 Scattering surface height fitting considering all scattering height estimations. We display fitting in red, and the red zone corresponds to the range of the flaring index at one σα. We display the scattering height estimation from Avenhaus et al. (2018): hs(r) = 16.17 (r/100au)1.219 measured on a different sample as a dashed gray curve for comparison.

3.2 Phase-function extraction

The next step is the extraction of the SPFs, which was also performed with the DRAGyS tool (Roumesy et al. 2025). We focused our study on a specific area of the disk’s surface, usually a ring, by defining an extraction zone delimited by two concentric circles with radii of R_in and R_out; these were projected according to the previously estimated disk geometry (see Table A.1). All SPFs were normalized to a 90° scattering angle, and the uncertainties include the errors on the estimated geometry, combined with the standard deviation within each bin. For the rest of our study, especially for adjusting the SPFs using the dust models presented later in Sect. 3.3, we decided to add a floor uncertainty of 0.05 to all SPF points as an estimate of the systematic uncertainty of the method.

For multiple-ring systems, we present the SPFs extracted for each ring in Appendix C. The SPF shape varies from one ring to another, which could reflect a radial grain-size distribution. However, this point is beyond the scope of the paper and will be explored in a future study. We focused the study on the best-defined ring for the following statistical study. In the case where the target is observed at several wavelengths, we keep the same extraction zone limits R_in and R_out for all wavelengths. The extraction zone covers the entire azimuth of the ring. However, it could overlap the coronagraph in some cases, or the transmission of the coronagraph could contaminate the surface brightness due to large inclination or strong flaring of the disk. This is the case for HD 34282, PDS 201, and V 599 Ori, where the SPFs are particularly affected for small scattering angles. To avoid such contamination, we decided to exclude the lower scattering angle for these SPFs. In some cases, other parts of the disk needed to be excluded to preserve the SPF. As an example, PDS 66 observed in the K band shows that the farthest part of the disk appears strongly affected by background noise. It may be useful to exclude the extracted SPF from these scattering angles. However, to keep the SPFs as consistent as possible, we decided to ignore this last point. This should be kept in mind when analyzing SPFs.

As a last step, the protoplanetary disk’s surface brightness is affected by the limb brightening effect, making the forward side of the disk brighter than the back side (Tazaki et al. 2023; Ginski et al. 2023; Roumesy et al. 2025). This purely geometric effect alters our extracted SPFs and must be corrected for in our measurements to better evaluate the intrinsic dust SPF lying on the disk’s surface. All corrected SPFs are displayed in Fig. 5, although only the one from the well-defined ring is shown for the multiple ring-shaped systems. More detailed SPFs are presented in Appendix C. A clear distinction is identified between two shapes of SPFs in Ginski et al. (2023). They define Category I as presenting a monotonic decrease in SPF with the scattering angle, whereas a bump characterizes Category II. However, when considering a larger sample of disks and after correcting for limb brightening that tends to flatten the SPFs, it appears that the observed SPF shapes follow a more continuous distribution. For simplicity (and direct comparison), and because we aim to identify (broad) families of dust shapes in protoplanetary disks, we retained this division into two categories of SPFs at this point. As a consequence, in some cases, this classification remains ambiguous. For example, PDS 66 shows a small bump in the J band, so we placed it in Category II, although it was classified as Category I by Ginski et al. (2023). Conversely, RX J1615.3-3255 exhibits an almost monotonic shape in the J band, but a marked bump in the H band. So, its first Category II classification is maintained (Ginski et al. 2023). The disks around PDS 70 and UX Tau present relatively flat profiles in the H band, but reveal a bump in the J band for UX Tau and in the K band for PDS 70. We therefore classified them as Category II. Lastly, targets such as HD 97048, HD 31648, and V1094 Sco also have relatively flat shapes, but their detailed SPFs (see Appendix C) indicate a monotonic decrease, justifying their inclusion in Category I.

3.3 Dust properties

To determine which dust model better explains the SPFs from category I and category II, we used the AggScatVIR (Aggregate Scattering for Visible and Infrared wavelengths) database, developed by Tazaki & Dominik (2022); Tazaki et al. (2023), where they computed the theoretical SPFs for 360 different dust models. The database is divided into three types of particles. The first type is compact aggregates formed by at least eight and at most 4096 equally sized monomers (100 nm, 200 nm, or 400 nm), providing a maximum grain size ranging from 0.22 to 3.14 μm and an equivalent volume radius between 0.2 and 1.6 μm. Compact aggregates are divided into three categories: low porosity (CALP), medium porosity (CAMP), and high porosity (CAHP). The second type is fractal grains, with a fractal dimension of Df = 1.1, 1.3, 1.5, and 1.9, composed similarly to compact aggregates, but with additional monomers of 150 nm and 300 nm in size. The fractal grains offer a wider range of maximum grain size (from 3.6 to 12 μm), and the corresponding equivalent volume radius ranges from 0.2 to 2.0 μm. Finally, irregular grains consisting of a Gaussian random sphere (GRS) are defined by their equivalent radius, which ranges from 0.2 to 1.6 μm.

AggScatVIR allows optical properties to be averaged over a particle size distribution n(a_v)da_v = a_v^−3.5da_v, where a_v is the volume-equivalent radius of particles. Hence, the maximum grain size of particles is larger than the size specified by this value for higher porosity (Tazaki et al. 2023). This distribution is a discrete sampling based on the number of monomers but not their size and on the volume-equivalent radius for the Gauss. Particle size distribution is applied throughout the rest of the study to better approximate the dust population in protoplanetary disks, and we also performed the same analysis when turning off this size distribution in Appendix B. AggScatVIR provides two different compositions. Both are a mixture of water ice (20% in mass), silicate (32.9%), and troilite (7.4%), to which either amorphous carbon (amc, Zubko et al. 1996) or organic matter (org, Henning & Stognienko 1996) is added.

For each category, the SPFs extracted from the disk presented in Sect. 3.2 are compared to the theoretical SPFs from the AggScatVIR dust models. The compatibility between observations and models is assessed using a χ² adjustment, and only the five models with the lowest values are selected for each case. This fitting process is shown for the disk around LkCa 15 in Fig. 6. This approach aims to identify families of dust models to perform a more representative statistical analysis. Thus, for each category, the retained dust models are combined, and the global distribution of results is presented in the histograms in Fig. 7.

The statistics obtained on the SPFs in Category I reveal that fractal aggregates with D_f = 1.1 and very small monomers of 100 nm are favored, with a strong preference for organic matter composition. CAMPs are also a reasonable candidate, with large monomers of 400 nm offering a maximum grain size comparable to the observational wavelengths. The two candidates also differ in their porosity, as fractals have a high porosity greater than 90% while CAMPs are more compact, with a porosity around 70%. For category II, SPFs fitting shows a predominance of irregular grains or small compact aggregates, both consistent with a maximum grain size smaller than 0.5μm. However, organic and amorphous carbon compositions are almost equally probable. The study of SPF shapes in about thirty observations of protoplanetary disks extends the first analysis proposed in Ginski et al. (2023). Although a continuous distribution of SPF shapes seems to emerge, it appears that the distinction between the two categories provides a fairly clear identification of two families of dust grains in protoplanetary disks: porous aggregates and small compact grains.

Fig. 5 Extracted SPFs from our datasets for J, H, and K bands. All SPFs are shifted vertically for clarity. The colors correspond to observed systems.

Fig. 6 SPF fitting for LkCa 15 in the K band using the three families of dust models from AggScatVIR: irregular grains constructed by a GRS and compact and fractal aggregates (from left to right). The five best models are shown as colors, and all others are displayed as faint dark curves.

Fig. 7 Statistical panel of best models estimated from the Category I (top) and Category II (bottom) SPFs. The colors are linked to the type of particles from AggScatVIR and displayed in the left column. We kept the same colors for the other histogram to track their distribution in monomer size, number of monomers, maximum grain size, porosity, and composition (from left to right). The gray bars represent the total number of all models considered to show the range of possibilities. Note that irregular grains (in red) are not composed of monomers, so they are not visible in the histogram related to aggregates. Furthermore, as they are dense particles, we fixed their porosity to 0.

4 Discussion

4.1 SPF shapes: Category I

To describe the first category of SPF, the best dust models clearly show a predominance of aggregate-shaped grains, especially with medium porosity (CAMP). To describe the first category of SPF, the best dust models clearly show a predominance of fractals. The best model is defined by a fractal dimension of D_f = 1.1, composed of organic matter with a monomer size of 100 nm, and very high porosity (greater than 90%). Furthermore, no general trend emerges for the maximum grain size (i.e., number of monomers), suggesting that the SPF provides more constraints on the monomer size within the aggregate than on the overall size of the dust particle. This global analysis to explain category I SPFs generalizes the results obtained by Tazaki et al. (2023), whose study of the SPF in IM Lup also suggests a fractal dust model. However, this study concluded that a highly absorbing composition (amorphous carbon) is needed to explain the SPF of IM Lup. Because the study of Tazaki et al. (2023) is for an individual case, it is perhaps not surprising that some differences arise compared to our more global trend.

In addition to fractals, compact aggregates may also be reasonable candidates to explain the shape of SPFs from Category I. Especially compact aggregate with medium porosity, made of large monomers measuring approximately 400 nm, with a total number ranging from 8 to 64 monomers. This results in a preferred maximum grain size ranging from 0.2 μm to nearly 3 μm. This wide range of grain sizes suggests that the SPF of aggregates depends more on the size of the monomers than on that of the overall particle. In other words, the shape of the particle, defined by its fractal dimension or porosity, strongly influences the dependence of the scattered light on the specific properties of the dust grains. Specifically, the SPF of a spherical dust particle or an irregular grain depends mainly on the overall size of the particle. In contrast, for compact porous aggregates or more complex fractal structures, the SPF is more influenced by the size of the monomers that compose the aggregate than by its overall size. A more thorough exploration of the influence of particle shape on SPF is therefore required to better understand these dependencies.

4.2 SPF shapes: Category II

For the second category of bell-shaped SPFs, the most promising dust models are irregular grains (GRS) and compact aggregates with low porosity (CALP), both with a maximum grain size of less than 0.5 μm. These two models do not exhibit a strong preference for any composition. Although there seems to be a slight tendency toward amorphous carbon (amc), the shape and porosity of the grains seem to govern the SPF more than their composition. Overall, sub-micron-sized dust particles that are highly compact generally explain SPFs with a bump.

An additional finding is the bimodal distribution of monomer sizes observed in CALP-type dust models, for which no solution was found at 200 nm in our sample. This is due to the lack of a bell-shaped profile in the SPF of the CALP models at 200 nm, unlike at 100 nm and 400 nm, where a bump is clearly visible (see Fig. 8). This behavior goes beyond the distinction of SPF in categories; a similar trend is observed in Category I, where the two preferred dust models are clearly split between 100 nm for FA11 and 400 nm for CAMP due to the characteristic slope of their respective SPFs. The distinction observed in monomer sizes is the result of an observational bias. Current instruments and the disk geometry do not provide access to lower scattering angles. As shown in Fig. 8, this part of the SPF could exhibit a bump, as is the case for 200 nm CALP, and access to these smaller angles would provide a much more robust constraint on the dust properties derived from SPFs.

Fig. 8 Theoretical SPFs of compact aggregates with medium (left) and low (right) porosity for different monomer sizes, shifted for clarity. Solid lines represent the typical range of scattering angles available in protoplanetary-disk observations, while dashed lines indicate the inaccessible range. For aggregates with medium porosity (a dust model that may explain Category I scattering angles), a gray area indicates the range of SPFs for Category I.

4.3 SPF asymmetry

Some SPFs from Category II show strong asymmetry coming from the disk’s surface brightness. A typical example in our sample is HD 163296, in which a significant brightness enhancement is visible on the north-west side. If we separate the contribution of SPFs from the east and west sides along the disk’s minor axis, we clearly identify that one side has a strong bump at around 70° and belongs to Category II, whereas the other one has a decreasing shape and belongs to Category I. In our study, we assume that the surface brightness is dominated by scattering by dust grains well mixed in the disk. However, this asymmetry in SPFs, also visible in the disks around HD 97048, HD 34282, PDS 70, PDS 201, RX J1615.3-3255, and V599 Ori, breaks our assumption for these disks, as the surface brightness is affected by physical processes beyond the properties of dust and raises the question of which disk portion can be used to estimate the intrinsic dust SPFs. These asymmetries may arise from the presence of planetary companions within the disks, which induce geometric perturbations such as warps, shadowing, or misalignment of the inner disk. This alters the distribution of scattered light on the disk surface, leading to asymmetric SPFs unsuitable for estimating the intrinsic properties of dust grains. Within our sample, several disks have been studied to detect or predict such companions, such as PDS 70 (Keppler et al. 2018; Blakely et al. 2025), RX J1615.3-3255 (de Boer et al. 2016; Avenhaus et al. 2018; Asensio-Torres et al. 2021), HD 97048 Pinte et al. (2019), and HD 163296 (Pinte et al. 2018; Izquierdo et al. 2026). To estimate the properties of dust in such cases, a radiative-transfer model is required to disentangle the intrinsic dust SPF from other contributions. We repeated the statistical analysis described in Sect. 3.3 to identify the best dust model for each category, excluding SPFs with strong asymmetry. This selection reduced our Category II sample but does not change the main results obtained before. It even reinforces the dominance of very compact sub-micrometric grain models such as GRS or CALPs, explaining the shape of Category II SPFs, as fractal grain models (D_f = 1.1 and 1.3) present in the initial statistics (see Fig. 7) disappear once the asymmetric SPFs are excluded.

4.4 Individual source: Degeneracy

This study highlights general trends regarding the best dust model to reproduce the shape of SPFs in each category, but individual interpretation may remain highly degenerate for some systems. This is the case for the disk around LkCa 15, for example (see Fig. 6); while fractal aggregates are the favored candidates, some compact aggregate models can also provide a credible fit. Nevertheless, the shape of the SPF still imposes useful constraints on some properties within each category, such as trends in maximum grain size, porosity, or composition. In future studies, it will be necessary to break down this degeneracy to examine the dust within the disks individually. One effective strategy is to use multiwavelength observations to extract SPFs and measure the color of the disk, which is directly related to the optical properties of dust grains and independent of the disk geometry. This is the approach adopted for IM Lup in Tazaki et al. (2023), where combining several bands allows us to link the shape of the SPFs and the disk color to a more restricted set of grain properties. Moreover, the wavelength dependence of the SPF and color is inherent to the dust model. So, comparing the same model to multi-band SPFs simultaneously adds constraints consistent across wavelengths and can break the degeneracy encountered when using a single SPF.

4.5 Dust growth and origin of aggregates

The growth of dust in protoplanetary disks begins with so-called first-generation sub-micronic grains originating from the interstellar medium, whose low relative collision velocities in the protoplanetary disk promote “hit-and-stick” coagulation, forming fractal aggregates with D_f ≤ 2 (Blum & Wurm 2008). As the aggregates become larger, collision velocity increases, leading to hit-and-stick collisions with deformation, while maintaining high porosity (Blum & Wurm 2008; Testi et al. 2014). When velocities increase further, compaction, fragmentation, or erosion phenomena occur, still producing porous aggregates but not fractals (Blum & Wurm 2008; Birnstiel et al. 2016). The most suitable dust models for SPF categories suggest that Category II (compact sub-micronic grains) corresponds to the initial growth phase with small compact grains from the ISM, while Category I (porous fractal or compact aggregates) may correspond to a more advanced growth stage involving compaction and/or fragmentation (Blum & Wurm 2008; Testi et al. 2014). However, since our sample consists of disks with a ring structure, often interpreted as more advanced systems where planet formation might be occurring, it is plausible that the two categories of SPFs, thus the two associated dust families, reflect distinct stages of the planet formation process (Turrini et al. 2019). Bernabò et al. (2022) shows that the amount of dust within the disks either remains stable or increases during the first 1 to 3 Myr. This is consistent with a scenario of collisions between planetesimals excited by the planet formation process that increases the dust population and generates a second generation of sub-μm-sized dust particles (Bernabò et al. 2022; Birnstiel 2024). Given current state-of-the-art methods, we cannot distinguish between the first generation of dust originating from the ISM - suggesting a very young disk without a planetary core - and second-generation dust, which suggests the presence of planetesimals or even a planet, depending on the rate of planetesimal formation (Birnstiel et al. 2016; Bernabò et al. 2022; Tatsuuma et al. 2024). This explains the lack of correlation between the age of the system and the SPF category in our sample, as already mentioned in Ginski et al. (2023), despite the influence of turbulence on collision velocities, which affects dust growth (Birnstiel 2024).

As an additional step, we tried to find a connection between the geometry and the category of dust we found using SPFs. As expected, we do not observe any connection between the dust families and the inclination of the disk, as well as for the position angle and the scattering scale height of the disk. We also tried it with the stellar spectral type, and again, no correlation appeared. The last test we carried out was on the ring structure of the disk. We estimated the ring width and its radius to the star by fitting the disk’s radial profile with a Gaussian for some disks in our sample where the ring structure is well defined (see Fig. 9). As a general view, the ratio between ring width and the radius seems to be smaller for the small compact grains (i.e., Category II) than for the aggregate particles (i.e., Category I). We applied the Mann-Whitney and Kolmogorov-Smirnov statistical tests to confirm the differences in ring widths, with confidence levels of 94% and 97%, respectively.

Assuming that rings in protoplanetary disks indicate a relatively advanced stage of planetary formation, this remark suggests that the narrow rings could host a population of second-generation small compact grains formed after an initial phase of planet or planetesimal formation (Tatsuuma et al. 2024). Conversely, broader rings would contain more porous grains such as aggregates, as they may still be in a dust-growth phase induced by coagulation. However, the observed trend does not constitute a strict correlation, and it remains difficult to draw definitive conclusions about the ring geometry and the dust population.

Fig. 9 Ratio between ring width and radial distance from the central star for 18 out of 30 disk observations. Measurements are separated according to the SPF-based classification: red markers indicate Category I, and blue markers indicate Category II. The dashed colored lines represent the mean value of ∆r/r for each category, where ∆r is the FWHM of the fitted Gaussian.

5 Summary

In this paper, we present a statistical analysis of a large sample of thirty Q_phi polarized-intensity images of protoplanetary disks, observed in the near infrared by the VLT/SPHERE instrument. We used the DRAGyS tool to measure the geometry of all these ring-shaped disks and extract the polarized SPFs corrected for the limb-brightening effect. The key results are listed below:

• The SPF categorization, defined in Ginski et al. (2023), was validated and completed according to the SPF shape. Category I shows decreasing SPFs, while Category II is distinguished by a bell shape of around 60-80°;

• Combining all scattering surface height measurements provides a global estimation of the disk flaring, with a powerlaw exponent of α = 1.357 ± 0.11, consistent with the value reported in a similar previous study (Avenhaus et al. 2018);

• The disks around HD 97048, RX J1615.3-3255, SZ Cha, and HD 34282 have several rings, allowing for the individual estimation of their flaring, ranging from 1.08 to 1.38. For disks observed at multiple wavelengths, the scattering surface height does not significantly change from one band to another.

This study proceeded with an estimation of the best dust models for each SPF category, using the AggScatVIR database (Tazaki & Dominik 2022; Tazaki et al. 2023), containing a variety of dust-particle shapes. A χ² adjustment between these models and the extracted SPFs enabled us to identify some general trends:

• Category I SPFs are dominated by fractal dust models composed of small monomers of 100 nm, with very high porosity (greater than 90%) and a composition clearly dominated by organic matter (org). However, compact aggregates (CAMP) seem to be a reasonable candidate too, with large monomers (400 nm) and lower porosity (70%). Although the size of compact aggregates appears to be similar to the observation wavelength (from ~1 to 3 μm), Category I SPFs are preferentially affected by monomer size, composition, and particle shape rather than grain size;

• Category II SPFs are clearly dominated by small submicrometric dust particles with irregular shapes or low-porosity dust aggregates. This category shows significant dependence on particle size and porosity, but not so much on dust composition;

• Furthermore, strong asymmetry is often observed in Category II SPFs. For example, the disk around HD 163296 shows SPFs of different categories on each side of the disk. These asymmetries could result from various mechanisms such as dust traps in pressure bumps, vortices, and alignment issues between the inner and outer disks, or even interactions with planetary companions forming within the disk. However, current observations are not sufficient to conclude which mechanism is correct or to infer the presence of planets based on these asymmetric SPFs;

• Dust models for each category can probe different stages of the grain growth scenario, as porous aggregates are formed after a hit-and-stick collision between small particles and can explain the Category II SPFs. However, two scenarios can explain the small compact grains for Category I SPFs. The first one is that sub-micron grains may originate from the ISM and be directly incorporated during the early stage of grain growth. Alternatively, the collision between planetesimals could fragment them, releasing these dense sub-micron grains (Tatsuuma et al. 2024). The ring-shaped structure of the disk may be caused by various mechanisms. One of these involves an advanced stage of planetary formation and could lend more credibility to the second scenario.

We emphasize that this statistical study provides a general trend on the dust properties in disks and will guide future studies on dust populations existing within disks. However, a similar study on a single disk is not enough to accurately evaluate the properties of the dust grains it contains, due to the high degeneracy in the results. To overcome this degeneracy, a multiwavelength study is required, not only to obtain multiple SPFs but also to access information about the disk’s color, which is directly related to the optical properties of the dust and independent of the disk geometry. Among our sample, some polarized intensity color was already estimated from NIR observations of disks such as V4046 Sgr, PDS 70 (despite its notable temporal variability Ma et al. 2024a), PDS 66, LkCa 15, and RXJ 1615.3-3255 (Ma et al. 2024b; Avenhaus et al. 2018).

The key to fully exploiting these additional constraints will be to understand precisely how the disk’s color varies as a function of different dust models. This is particularly important for the aggregate models simulated with AggScatVIR (Tazaki et al. 2023; Tazaki & Dominik 2022), which remain the most plausible candidates for populating protoplanetary disks.

Acknowledgements

This project has received funding from the European Research Council (ERC) under the European Union’s Horizon Europe research and innovation program (grant agreement No. 101053020, project Dust2Planets, PI: F. Ménard). RT was supported by JSPS KAKENHI grant Number JP25K07351.

References

Appendix A Geometric parameters of the disks

In this section we present the geometric parameters estimated using DRAGyS (Roumesy et al. 2025) for all of the protoplanetary disk observations presented here. These results are listed in Table A.1.

Appendix B Statistical study using single-size dust particle

In this section we suggest estimates for the best dust particles for each SPF category, without applying grain size distribution, and compare this with the case where this distribution is applied (see Fig. 7). The first observation is the strong emergence of CAHPs to explain category I SPFs. This dust model becomes the best model, with a monomer size of 200nm, offering a dust size of ~ 2 μm and a relatively high porosity of over 80%. In this case, and unlike the fractals with a dimension D_f = 1.1 that dominate when the size distribution is taken into account, amorphous carbon becomes the preferred composition. As for category II, there is also an increase in the number of CAHP considered to be reasonable candidates, but small aggregates with low porosity and irregular grains still dominate, with submicrometric grain size and extremely low porosity, without any major constraints on composition.

Fig. B.1 Same as Fig. 7 but without grain size distribution.

Appendix C Detailed SPF for all observations

In this section we present the polarized SPF extracted for each system, including when considering each side of the disk separately. Figure C.1 corresponds to the datasets.

Fig. C.1 Panel of SPF extraction from each side of the ring. In each case, we display the Qφ image with and without the extraction zone on the left, and the extracted polarized SPF on the right. We plot SPF for the entire disk (solid black), east side (dotted red), west side (dotted blue), and without correcting for the limb-brightening effect (dotted gray).

All Tables

All Figures

	Fig. 1 VLT/SPHERE Qφ observations in the J-, H-, and K-bands. All images have dimensions of 3.4″ × 3.4" and are shown using an individual log-scale color map.
In the text

	Fig. 2 Scattering-surface-height fitting for multiple ring-shaped disks. We plot adjustments as dashed colored curves, and we display each estimated value of scattering height ratio, h0/r0, and flaring exponent in Table 1.
In the text

	Fig. 3 Ring fitting on HD 97048. We compare the estimation presented in Ginski et al. (2016) (in orange) and our estimation (in cyan). The image is on R2-scaled, and we only plot half of the ellipse to bring up substructures. Features to the northwest are artifacts from data reduction.
In the text

	Fig. 4 Scattering surface height fitting considering all scattering height estimations. We display fitting in red, and the red zone corresponds to the range of the flaring index at one σα. We display the scattering height estimation from Avenhaus et al. (2018): hs(r) = 16.17 (r/100au)1.219 measured on a different sample as a dashed gray curve for comparison.
In the text

	Fig. 5 Extracted SPFs from our datasets for J, H, and K bands. All SPFs are shifted vertically for clarity. The colors correspond to observed systems.
In the text

	Fig. 6 SPF fitting for LkCa 15 in the K band using the three families of dust models from AggScatVIR: irregular grains constructed by a GRS and compact and fractal aggregates (from left to right). The five best models are shown as colors, and all others are displayed as faint dark curves.
In the text

Fig. 7 Statistical panel of best models estimated from the Category I (top) and Category II (bottom) SPFs. The colors are linked to the type of particles from AggScatVIR and displayed in the left column. We kept the same colors for the other histogram to track their distribution in monomer size, number of monomers, maximum grain size, porosity, and composition (from left to right). The gray bars represent the total number of all models considered to show the range of possibilities. Note that irregular grains (in red) are not composed of monomers, so they are not visible in the histogram related to aggregates. Furthermore, as they are dense particles, we fixed their porosity to 0.

In the text

Fig. 8 Theoretical SPFs of compact aggregates with medium (left) and low (right) porosity for different monomer sizes, shifted for clarity. Solid lines represent the typical range of scattering angles available in protoplanetary-disk observations, while dashed lines indicate the inaccessible range. For aggregates with medium porosity (a dust model that may explain Category I scattering angles), a gray area indicates the range of SPFs for Category I.

In the text

	Fig. 9 Ratio between ring width and radial distance from the central star for 18 out of 30 disk observations. Measurements are separated according to the SPF-based classification: red markers indicate Category I, and blue markers indicate Category II. The dashed colored lines represent the mean value of ∆r/r for each category, where ∆r is the FWHM of the fitted Gaussian.
In the text

	Fig. B.1 Same as Fig. 7 but without grain size distribution.
In the text

	Fig. C.1 Panel of SPF extraction from each side of the ring. In each case, we display the Qφ image with and without the extraction zone on the left, and the extracted polarized SPF on the right. We plot SPF for the entire disk (solid black), east side (dotted red), west side (dotted blue), and without correcting for the limb-brightening effect (dotted gray).
In the text