© The Authors 2025

1 Introduction

Recent observations have confirmed that protoplanetary disks are the birthplaces of planets (Keppler et al. 2018; Pinte et al. 2018; Teague et al. 2018). Beyond these fascinating discoveries, instrumental improvements have led to the detection of a wide morphological diversity in disks, with spiral, gaps, and rings being prominent and reflecting the complex planetary formation process (Benisty et al. 2023). Indeed, submillimeter-sized dust grains must agglomerate to form planetesimals in a very short time to explain the birth of planets (Drazkowska et al. 2022). The size and structure of dust aggregates influence their ability to grow (Johansen et al. 2014), and understanding these properties is important to better track the evolution of dust in disks.

The interaction of an electromagnetic wave with a dust particle depends on the wavelength, on the size and shape of the particle, and on the refractive index of the material. Therefore, measuring the scattering phase function (SPF), describing the intensity of scattered light as a function of the scattering angle (e.g., Bohren & Huffman 1983), can be a powerful tool with which to study the dust’s structural properties (e.g. Ginski et al. 2016; Stolker et al. 2016; Tobon Valencia et al. 2022; Ginski et al. 2023; Tazaki et al. 2023; Tobon Valencia et al. 2024; Chen et al. 2024).

In optically thin debris disks, for which only single scattering is relevant, the interpretation of extracted SPFs is direct (Milli et al. 2017). The potential for significant absorption as well as multiple scattering and geometric complexity in protoplanetary disks makes the interpretation of their SPF more challenging. Yet, the image quality (angular resolution, contrast) provided by modern-adaptive-optics-enhanced imaging instruments has led to a dramatic increase in the number of star and disk systems for which SPFs can be reliably extracted (Ginski et al. 2023). This opens the way to a deeper exploration of dust grain properties in the early stages of their evolution.

To date, a few tools have been developed to extract the SPF from protoplanetary disk surfaces, such as Diskmap (Stolker et al. 2016). However, these tools require a prior knowledge of the disk geometry, which is generally obtained by a customized ring fitting (Boer et al. 2016; Avenhaus et al. 2018), a process that we generalize here, or by full radiative transfer model fitting. In the latter, the time-consuming model-fitting process is an obstacle to large-scale studies of dust grains in protoplanetary disks using SPF. To overcome this constraint, we developed a new comprehensive tool, DRAGyS, to analyze the structure and scattered light properties of protoplanetary disks without prior information. Our approach stands out by its ability to estimate the disk geometry directly from an image, and directly extract the SPF. This method offers a fast and efficient alternative to existing methods, making it particularly well suited to the study of large samples of protoplanetary disks, and a better understanding of the properties of dust grains they contain.

In Section 2, we present the method of estimating the geometrical parameters of ring-shaped disks; that is, their inclination, position angle, and aspect ratio. Then, we outline the steps involved in the extraction of the SPF from the surface of the protoplanetary disks. In Section 3, we first apply our new tool to synthetic images of disks, then evaluate the impact of limb brightening and show how it can be corrected for. In Section 4, we test DRAGyS on a sample of six ring-shaped disks observed in the near-infrared. These disks were chosen because their SPF have already been extracted, and therefore serve as useful benchmarks (Ginski et al. 2023). In Section 5, we describe the advantages and limitations of this new tool. Our conclusions are summarized in Section 6.

Fig. 1 Illustration of the DRAGyS method, illustrating the ellipse estimation of the disk (red), and the geometrical method of recovering the different structure parameters. Left: image seen by the observer, with the ellipse in red and different estimated quantities, such as the semiminor and semimajor axes, m and M, the position angle, PA, and the projected distance between the star and the major axis, D, used to compute the aspect ratio. Right: sketch of protoplanetary disk observation. The point of view of observer 1 highlights the geometrical process to recover the disk geometry, namely the scattering surface height, h, at a radius, rmid, and the inclination, i. The point of view of observer 2 illustrates the definition of the scattering angle.

2 Methodology

DRAGyS is designed to avoid the time-consuming step of model fitting to determine the disk structure necessary to extract the SPF. Instead, it estimates the disk geometry by assuming the disk is circular and by exploiting visible structures such as rings that are used to delineate the position of the disk surface where scattering occurs. For geometrically thick disks, these elliptical rings are offset from the center, and the offset contains the information about the thickness and flaring of the disk surface (Lagage et al. 2006; Avenhaus et al. 2018).

DRAGyS provides a comprehensive analysis of the disks, from their geometric parameter estimation to the SPF extraction. It follows three steps: (i) the detection of the brightness peaks on the image, (ii) fitting by ellipses to estimate the disk geometry, and finally (iii) extraction of the SPF over the full disk or a restricted sector only, using the estimated geometric parameters.

We characterized the disk structure by the aspect ratio of the scattering surface, defined as the disk scattering surface height, h, over the midplane radius, r, from the star. Throughout our study, the power law exponent that governs the flaring of the disk scattering surface is χ, and we refer to it as the scattering surface flaring exponent. We defined the local height of the protoplanetary disk as $h(r)=h_{\text{ref}}{\left(\frac{r}{r_{\text{ref}}}\right)}^{\chi },$(1)

where h_ref is the height of the disk at a distance, r_ref, from the star. The orientation of the disk is important to determine the scattering angles necessary to construct the SPF. So, DRAGyS also fits for the inclination and the position angle of the disks. Together, these three parameters (respectively. h(r), i, and PA) define the geometry of the disk surface (see Fig. 1, right).

The method developed here is suitable for any disk observed in scattered light, in total and polarized intensity, at any wave-length. It is not limited by the geometric thickness of the disk, and can also be applied to optically thin disks. This approach requires the presence of well-defined concentric structures, such as rings. Most disks around T Tauri stars have a ring structure (Rapson et al. 2015; Thalmann et al. 2016; Monnier et al. 2017; Bertrang et al. 2018; Muro-Arena et al. 2018; Avenhaus et al. 2018; Marel et al. 2019; Benisty et al. 2023). Incident starlight is attenuated by dust close to the central star (Boer et al. 2016; Ginski et al. 2016). Thus, we assume that maximum flux comes from the surface of the disk, illuminated by the star, and not at all or only partially affected by the inner rim of the ring. It is also necessary that instrumental effects (e.g., coronagraphic mask, imperfect PSF subtraction, low throughput) do not significantly alter the measured surface brightness.

2.1 Step 1: Identification of the ring structure

The first step relied on image analysis and searches for maxima along radial profiles spaced by 1° in azimuth (see Fig. 2). That allows determining the location of the elliptical shape formed by the luminosity peaks within the ring from which the SPF is studied. Spatial and intensity filtering can be applied to optimize this detection, avoid confusion with other rings in the system, and eliminate spurious peaks due to weak local maxima in simulated images, or noise in observations (see Section 4). In addition, we performed a manual selection of peaks associated exclusively with the selected ring. Then, we down-selected peaks that were spaced by a resolution element to ensure that the fitting was not biased by a cluster of correlated points, for example, at the ellipse ansae. For this identification process and for SPF extraction, we did not use an r²−scaled image, contrary to Stolker et al. (2016); Ginski et al. (2023) studies, as we did not have prior knowledge of the disk geometry.

Fig. 2 Illustration of the elliptical fit process on a Qphi simulated disk image with multiple rings. Left: disk image with a radial cut at one azimuthal angle in white. Points used for ellipse fitting are shown in green and unwanted deleted points in smaller orange dots. Middle: radial profile extracted along the white line, where the purple curve corresponds to the smoothed profile. The maxima corresponding to the selected ring and the undesired maxima from other rings are shown in green and red, respectively. Right: disk image with the fit ellipse shown in white. A set of random points used to estimate uncertainties is plotted in blue.

2.2 Step 2: Estimation of geometric parameters

The positions of the previously detected brightness peaks define an elliptical shape, whose parameters such as its center (X_c, Y_c), semiminor (m) and semimajor (M) axes, and position angle (PA) can be estimated (see Fig. 1, left). For this purpose, we used the least-square ellipse fitting method described in Halíř & Flusser (1998).

Then, the estimated ellipse parameters were transformed into disk geometrical parameters, taking into account the effects of projection onto the image plane (see Fig. 2). The inclination was determined by the ratio of the semiminor and semimajor axis. The position angle, PA, was measured in the usual way, from the horizontal left axis to the ellipse semiminor axis corresponding to the disk near-side and counterclockwise. Finally, the scattering surface height, h, at a given midplane radius, r_mid = M, was defined as the distance, D, between the estimated ellipse center and the central star, taking into account projection effects due to disk inclination¹. We defined the aspect ratio of the scattering disk surface as h/r = D/(r_mid sin i).

Uncertainties associated with the disk geometric parameters were estimated by propagating the ellipse fitting error. Concretely, we computed the standard deviation, σ, of the relative distances between each maxima point and the fit ellipse. From this standard deviation, we generated N_p sets (typically 100 sets) of points randomly distributed around the ellipse, following a normal distribution in radius using r_mid for the mean and σ for the standard deviation. The ellipse fitting process was applied to each of these sets to obtain a statistical distribution of geometric parameters. The uncertainty of each parameter was finally defined as the standard deviation of the estimated parameter values across the N_P sets, providing a robust measure of the errors.

2.3 Step 3: SPF extraction

The SPF must be extracted from a restricted zone bordered by two ellipses located at the estimated disk surface. They were computed using estimated geometrical parameters, and followed the formulation from Stolker et al. (2016). To accurately capture the flux from a ring, the extraction zone was limited by minimum and maximum midplane radii. This was necessary to exclude unwanted areas, such as gaps or other rings, which could alter the SPF, but also to focus the study on a well-defined part (radius) of the disk. This method therefore also allows one to treat each ring separately, without assuming the same SPF throughout. This is particularly interesting for disks with multiple rings that could host different dust populations. It enables the analysis of the dust properties as a function of radius. In the same way, it enables one to observe the radial dependence of light scattered on the surface of the disk by probing zones that are more or less distant from the center. Consequently, the extraction zone was chosen to be as wide as possible, within the ring observed. The influence of its position and width is discussed in Appendix B.

To construct the SPF, the flux was extracted for each pixel in this extraction zone, and the corresponding scattering angle, Ψ, was computed as $\cos (\text{Ψ})=\cos (\gamma )\cos (\varphi )\sin (i)+\sin (\gamma )\cos (i),$(2)

where γ is the disk opening angle, which corresponds to the angle between the midplane and the disk scattering surface as seen from the star.

To incorporate uncertainties from estimated geometric parameters, the SPFs were extracted over a fixed sampling in scattering angle, using N_s different sets of parameters randomly taken at 1σ around their estimated value. This yields N_s SPFs and their error given by the mean value and standard deviation of flux, respectively, in each bin. We averaged these SPFs to obtain the final curve, and the standard deviation of the mean gives the uncertainty linked to the geometric parameters. Lastly, we quadratically combined this uncertainty with the one from the flux extraction to obtain the total uncertainty of the extracted SPF.

2.4 Features and requirements

DRAGyS is primarily designed to provide a complete analysis of ringed protoplanetary disks. It determines the disk geometry and extracts the SPF, covering all azimuths as well as each side separately to reveal eventual asymmetries. However, DRAGyS can extract the SPF from any disks, including smooth ones, provided the disk surface geometry is known a priori (or modeled). It further includes several advanced options, relevant in specific cases. For example, it handles situations in which the actual location of the star is offset from the center of the disk, and automatically corrects the extracted SPF from the limb brightening effect, whose details are discussed in Section 3.3. The tool deals with both total and polarized intensity equally well. To simplify the rest of the paper, we focus on polarized intensity results only.

3 Benchmark tests with simulated data

3.1 Making a mock observational image

As a first test, we applied DRAGyS to simulated images of protoplanetary disks. These simulations were performed using the 3D radiative transfer code MCFOST (Pinte et al. 2006, 2009). This software follows the Stokes formalism (Bohren & Huffman 1983) to generate I, Q, and U images of the disk. The Stokes formalism allows one to follow the evolution of a light beam through successive scattering events. On top of the disk images, MCFOST also outputs the theoretical SPF expected from the dust distribution used in the model. This allows for a direct comparison between the SPF extracted from a simulated disk image and the theoretical SPF computed by MCFOST (see Fig. 4, right).

Synthetic disks are located at a distance of 100 pc and composed of three rings at 60–90 au, 110–140 au, and 160–190 au. They are composed of silicate grains (Draine & Lee 1984) with a power-law size distribution, n(a)da = a^−pda, ranging from 0.005 to 3 μm with a power-law index, p = 3.5 (Dohnanyi 1969). Although these grains may not be representative of those expected in protoplanetary disks, they produce a complex SPF with multiple minima and maxima, placing maximum stress on the method. This is well suited for testing DRAGyS and evaluating its strengths, biases, and weaknesses. We varied parameters such as the position angle (PA = 70°, 200°, 270°, 300°) and inclination (i = 20°, 50°, 70°) to test the robustness of this method. The gas pressure height was defined as h_gas = 5 au(r/100 au)^1.125. The density profile follows a power-law, Σ(r) ∝ r^−p, with surface density exponent, p = 0.5. The dust mass was fixed to M_d = 10⁻⁴ M_⊙ in the fiducial model. All computed Q_phi images (Schmid et al. 2006) were obtained at 1.6 μm, and convolved using the VLT/SPHERE instrument diameter. We display all simulated Q_phi images in appendix A.

3.2 Geometric parameters fitting

We directly compared the input model values with the estimated inclination, position angle, and scattering height values extracted with DRAGyS. The models are composed of three rings, which also enable an analysis of the potential dependence on the position of the ellipse used for estimation. The estimated parameters for each ring are listed in Table 1.

An overall comparison shows that parameter estimation tends to improve with increasing radius. Analysis of the inner ring shows larger discrepancies with reference values, both for the position angle and inclination, mainly due to the smaller number of pixels available for fitting when looking at the inner ring. For the same reason, measurement errors decrease with increasing radius. Specifically, our method recovers the inclination within ∼ 2° of the theoretical value for all cases. The worstcase error is for the inner ring estimation of the 70° inclined-disk. We also observe an underestimation of inclination of slightly less than 1° on average. In terms of position angle, the estimated value nearly always includes the theoretical value within its confidence range, with a maximum offset of ∼ 1.5°, when analyzing a weakly inclined inner ring. Regarding the measured scattering height, we compare it with the τ_s = 1 surface computed from the star. This surface corresponds to the last scattering event that emerges from the disk (see Fig 3). For very optically thick and flared disks, as is usually the case for protoplanetary disk in the optical and near-infrared, the photons reaching the observers come mainly from the disk surface and this surface is well estimated by this τ_s = 1 surface (see Appendix D for a comparison of the measured surface versus the τ_s = 1 surface for various disk masses). For the inner ring, the estimations are all lower in the disk than the τ_s = 1 surface, and more spread out. This is due to the direct illumination of the inner rim, combined with the decreasing number of points available for ellipse fitting when approaching the central star. Our method reaches its limits as the distance between the projected stellar position and the ellipse center becomes smaller, and thus harder to measure precisely. Comparing the estimated scattering heights with the pressure scale height used by MCFOST gives a mean ratio of ∼ 3.7, within the same range as was found in previous studies on different targets (Chiang et al. 2001; Avenhaus et al. 2018; Ginski et al. 2016). We discuss this comparison for several disk dust masses in Appendix D.

For these tests, the uncertainties are based on noiseless simulated data. They represent a minimum uncertainty intrinsic to the parameter estimation method. In other words, the estimated geometric parameter uncertainties are optimistic relative to the case of a real dataset.

Fig. 3 Comparison of estimated scattering height estimation with the τs = 1 surface integrated from the star (red) for each ring (respectively between 60–90 au, 110–140 au, and 160–190 au). We overplot the scattering height surface estimated using DRAGyS, for each inclination and position angle, in different colors and markers. We have added the dust density map in the background and, we also plot the gas pressure scale height with a dashed gray curve, using the input MCFOST parameters.

Fig. 4 Typical case of SPF extraction for simulated data, using a disk with PA = 270° and i = 70°. Left: disk image in polarized intensity, with an extraction zone defined by the green-shaded annulus. Right: polarized SPF normalized at 90° directly extracted from the disk surface (gray) and corrected from the limb brightening effect (black). The theoretical intrinsic SPF is plotted with a dashed purple line.

3.3 Biases in extracted SPF

We compared the “true” (intrinsic) SPF from the model, with the SPF directly extracted from the disk surface (see Fig. 4). The extracted SPF has the correct qualitative behavior, with troughs at the same scattering angles and peaks roughly at the correct position. The shift in the maximum of the extracted SPF is due to two major effects that modify the shape of the SPF.

The first one is multiple scattering, which is characterized by photons that undergo several scattering events on the dust grains inside the protoplanetary disk before reaching the surface and leaving the disk. Tazaki et al. (2023) show that multiple scattering affects the polarized SPF of protoplanetary disks, particularly when the single scattering albedo of grains is large. There is no simple prescription to correct for this effect, since it depends on the disk structure and optical depth.

The second one is the limb brightening effect, which results in the near side of the disk appearing brighter than the far side. In the context of optically thick disk observations, Tazaki et al. (2023) show that this purely geometric effect is caused by the optical depth, which varies according to the location on the inclined disk. In other words, the optical depth measured from the surface to the observer does not reach unity at the same height in the disk over all azimuths, which significantly influences the intensity received from different surface parts, and then, the measurement of the SPF. The method developed here estimates, and can correct the SPF for, the limb brightening effect, following the formulation from Tazaki et al. (2023): $\text{LB}=\frac{\cos (\gamma )\sin \left({\gamma }^{\prime }-\gamma \right)}{\cos (\gamma )\sin \left({\gamma }^{\prime }-\gamma \right)+\cos (i)\cos \left({\gamma }^{\prime }-\gamma \right)-\sin \left({\gamma }^{\prime }\right)\cos (\text{Ψ})},$(3)

where γ = h/r and γ^′ = dh/dr = χh/r.

The limb brightening effect plotted in Figure 5 shows that there is only a slight dependence on the disk flaring exponent.

On the contrary, limb brightening is dominated by the aspect ratio. We therefore assumed a fixed flaring exponent value close to unity for the DRAGyS SPF extraction process and for the limb brightening effect estimation. Uncertainties of geometrical parameters were propagated through the limb brightening correction process. Figure 4 shows that the limb brightening-corrected SPF is a much better match to the expected SPF from the model.

Fig. 5 Limb brightening effect for different disk geometries. The forward scattering part is more affected than the backward scattering.

4 Application to observational data

Several disks have clearly visible rings for which the SPF has already been estimated in the past (Ginski et al. 2023), which is useful to test the DRAGyS extraction process. Here we focus on six protoplanetary disks observed in scattered light with the Spectro-Polarimetric High-contract Exoplanet REsearch facility (SPHERE; Beuzit et al. (2019)) at VLT, using the IRDIS near-infrared camera (Dohlen et al. 2008). To obtain a polarized intensity image of the protoplanetary disk, the instrument is operated in dual-beam polarization (DPI) mode and observations are made in the J and/or H bands (see Fig. 6). It should be noted that although the data are obtained in the classic Stokes Q and U formalism, we used the azimuthal Stokes Q images, Q_ϕ, to carry the SPF extraction from each observation (Schmid et al. 2006).

Geometric fitting of protoplanetary disks from real observations presents a few challenges. The presence of noise complicates the identification of true intensity maxima, so additional steps must be applied to the fitting process compared with the method presented in Section 2. These consist essentially of a spatial averaging of the image and a more aggressive filtering on the intensity maxima detection.

We list the parameters estimated with DRAGyS in Table 2 and maps of points detected for ellipse fitting are shown in Figure 6. The comparison between the reference geometrical parameters and the ones estimated with DRAGyS yields consistent results in most cases. Larger deviations in the PA are sometimes found, especially when the disk inclination is low; for example, in the PDS 66 (J-band), V4046 Sgr (H-band), and RXJ 1852 at worst. For the aspect ratio, some estimations are consistent between our method and the reference data, but there can be large discrepancies, as in the case of RX J1852 or LkCa 15. Aspect ratio estimations also vary from band to band. Therefore, the best way to compare our estimation is to observe the extracted SPF.

The geometric parameters estimated in Table 2 were applied to extract the SPF from an extraction zone on the disk surface using DRAGyS (see Fig. 7). We compared our results with reference SPFs using Diskmap and reference geometries from Ginski et al. (2023), keeping the same width and position of the extraction zone for both processes. The difference between the two SPFs was quantified using the relative error over the common interval in the scattering angle for both SPFs.

Both methods cover very similar scattering angle ranges. As a result, we typically obtain a good agreement in terms of the SPF, with differences in the range of 2% to 18%. As with the geometrical parameter estimation, the most problematic case remains V4046 Sgr, where the scattering angle ranges and shapes are not consistent with Diskmap. This system presents the largest difference (∼ 17% in both band), in line with the fact that its geometric estimation was the worst in the sample. If the defined geometry is different, SPFs will usually be different. In contrast, a much better match is observed in the case in which the geometries are consistent for both methods, such as for RX J1615, HD 163296 (H-band), and PDS 66 (J-band), with a percentage difference of 2.57%, 4.25%, and 3.25%, respectively. These two extremes demonstrate the significant influence of the estimated geometry on the SPF, but also the consistent SPF produced by the two methods for similar geometries.

Now, consider the SPF corrected for limb brightening, which is estimated with DRAGyS (see Fig. 7). The general shape of the curves does not significantly change, with troughs and peaks still present at the same scattering angles, but the slope becomes shallower in all the cases considered here. This deviation may vary according to the disk geometry: it can be negligible in some cases, such as for V4046 Sgr and RX J1852 in the J and H bands, and strong enough to significantly modify the SPF in other cases, as occurs for LkCa 15, RX J1615, and HD 163296. Moreover, we also computed the percentage difference between SPF with and without correction. Limb brightening typically impacts the SPF by about 10% on average, and can even reach ∼30% in some cases. It is therefore important to correct for it before interpreting SPFs.

In summary, the method described here not only provides a relatively accurate estimation of SPF, but also corrects an effect that can have a significant amplitude depending on the observations, which was not corrected in previous results. We discuss this point in Section 5.

Fig. 6 Qϕ images of six protoplanetary disks observed in J or H bands using DPI with the VLT/SPHERE IRDIS instrument. The maximum pixel positions for parameter estimation are shown as white dots, and the fit ellipse with a dashed black line.

5 Discussion

5.1 Limiting factors

The key factor for SPF extraction is the positioning of the extraction zone. A correct positioning ensures accurate coverage of the disk scattering surface (here, the ring). Incorrectly estimated geometric parameters lead to a mismatch between the extraction zone and the ring. In the worst case, this can lead to the inclusion of zones outside the disk surface, or that come from another ring, leading to an incorrect SPF.

For simulated disk images, with known inclinations and PAs, DRAGyS finds a systematic underestimation of ∼ 1° for the inclination, but the correct PA. The exact value depends on the disk model, but also on the manual filtering of the intensity maxima (see section 3.2). Obviously, it is harder to estimate the position angle for a weakly inclined disk, whereas a highly inclined disk makes inclination estimation more difficult, and even more if combined with a large flaring. For VLT/SPHERE images, the true inclination values are unknown. DRAGyS recovers parameters similar to the ones found by Ginski et al. (2016).

The extracted SPF depends on the intrinsic SPF, limb brightening, and multiple scattering effects. The limb brightening effect can amplify the flux received at low scattering angles and must be corrected for. However, multiple scattering remains an uncalibrated factor. It can introduce a significant bias into the polarized SPF, especially since multiple scattering leads to a reduction in the net polarization fraction due to interactions between light and dust grains, especially at small scattering angles (see Fig. 9 in Tazaki et al. 2023). This effect can be evaluated a posteriori using radiative transfer models. While desirable, the development of a method to quantify this effect and obtain a closer-to-reality SPF is beyond the scope of this paper. For the simulation described here (Fig. 4), the extracted SPF is almost fully consistent with the intrinsic SPF of dust grains. However, there remains the multiple scattering effect, which depends on dust properties that have only been partially explored in our simulations. To further assess the impact of this effect, we extract SPFs from disks with different dust compositions and grain size distributions in Appendix E. Overall, we find that multiple scattering only has a modest influence on SPFs across the tested conditions.

As these two effects of limb brightening and multiple scattering have opposite influences on the polarized SPF, they may counterbalance each other. These findings underscore the importance of considering both geometric and physical factors when interpreting SPFs.

5.2 Achievements and highlights

Despite these limitations, DRAGyS accurately retrieves the disk geometric parameters on both simulated and telescopic data. It allows for the fast and reliable extraction of SPFs without prior knowledge of disk geometry.

The accuracy of the SPF recovery depends on the correct positioning of the extraction zone over the disk surface. It can be affected if the zone is not well matched to the ring position and width (see Appendix B).

High-inclination disks are particularly sensitive to errors in the extraction zone position. However, a modest error on the inclination (± 2°), or on the aspect ratio (typically 0.025), does not significantly alter the extracted SPF on this type of disk (see Appendix C). The simulations performed here are focused on disks with narrow, closely spaced rings, an extreme scenario that demonstrates the robustness of the tool on complex structures.

Concerning the advantages of the method, DRAGyS allows one to extract the SPF in total and polarized intensity from the protoplanetary disk surface, avoiding long model-fitting computations for each disk, thanks to a geometric fitting approach directly on the disk images. In addition, the method presented here enables the complete study of a wide range of disks from geometric fitting on the disk image to SPF extraction and limb brightening correction. The only constraint is that a structure such as a ring is required for the geometrical fitting process only. But it is also possible to manually define geometrical parameters, bypass that part of the fitting process, and apply the SPF extraction directly for disks with already known geometry.

A capability of DRAGyS that we have not used here, but that will prove powerful in the future, is the possibility to separate the extracted SPF into two zones, typically on either side of the minor axis of the disks. Indeed, observations often reveal disks with brightness asymmetries on their surface. These could be due to an asymmetric dust density distribution around the disks, but also because of shadows cast by a misaligned inner disk, for example. These asymmetries can therefore bias the inferred dust properties and DRAGyS provides an automated process to identify systems with significant asymmetry, and allows the user to extract the SPF from a subset of the data that is deemed uncontaminated by intrinsic asymmetries.

We explored the SPFs extracted for all disks considering each half separately (see Appendix F). This is particularly interesting for sources HD 163296 and RX J1615, where the brightness ratio between the disk ansae is around 2.75 and 1.75, respectively. This ratio is close to unity for the other sources, which does not have a noticeable impact on the SPFs.

Fig. 7 Polarized SPFs extracted from SPHERE Qphi observations (Ginski et al. 2023) of protoplanetary disks. The left and right panels are for H and J band observations, respectively. Solid and dashed colored lines are, respectively, SPFs uncorrected and corrected for the limb brightening effect, extracted using DRAGyS with fit geometric parameters. Dotted black lines show the reference SPFs from Diskmap (Stolker et al. 2016) following the geometric parameters from Ginski et al. (2023). Limb brightening uncorrected and Diskmap SPFs are normalized to 90°. We normalized the corrected SPF to the minimum scattering angle in common with the uncorrected one, emphasizing the slope change induced by the limb brightening correction. An offset has been added between each curve for better rendering. The percentages of the mean offset between SPF before and after removing the limb brightening are shown in color, and the ones between DRAGyS (without LB correction) and Diskmap are in black.

6 Conclusion

Here, we present DRAGyS, a new method of extracting the SPF by analyzing a disk from geometric parameters estimation to SPF extraction. DRAGyS represents a major advance in the study of the SPF extracted from ring-shaped protoplanetary disks. Unlike previous methods requiring prior modeling of the disk, this tool directly uses images to estimate the geometry and extract the SPF, while correcting for geometric effects such as limb brightening, which are often neglected but which can significantly influence the shape of the SPF.

Tests on synthetic images have shown that the geometric parameters of the disk, such as the inclination and position angle, can be estimated within 1° accuracy, although minor biases may occur, particularly in the case of disks with a high inclination or more complex structures. These estimations are made by detecting an elliptical shape created by brightness peaks on the disk image. The application to real data from polarized nearinfrared observations demonstrates the robustness of the tool, which provides consistent results with existing techniques such as Diskmap while offering advantages in terms of speed and execution simplicity. Therefore, DRAGyS is particularly suitable for studies of large disk samples.

Finally, this tool also has the potential to improve the characterization of asymmetric or complex disks, by enabling the separate extraction of SPF for different regions of the disk. The correction of limb brightening, combined with a detailed analysis of surface asymmetries, could open up new perspectives in the study of dust grain properties and planet formation in protoplanetary disks.

Data availability

The DRAGyS tool is available in GitHub repository: https://github.com/mroumesy/DRAGyS

Acknowledgements

We would like to thank the anonymous referee whose intuition led to a better paper. 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).

Appendix A MCFOST disk models for testing DRAGyS

In this section, we display images of the twelve disks simulated with MCFOST (Fig. A.1) used to test the robustness of our new tool, both for geometric disk estimation and for SPF extraction.

Fig. A.1 Panel of the twelve simulated disk from MCFOST. All disks have the same structure with 3 rings, but for three inclinations (columns) and four position angle (rows).

Appendix B Extraction zone impact on SPF extraction

In this section, we investigate the impact of the width and position of the extraction zone on the accuracy of the extracted SPF. The aim is to optimize the positioning of this zone to obtain the most accurate SPF possible. We simulate a protoplanetary disk with three rings at 30-50 au, 70-150 au, and 170-190 au, and we focus on the intermediate one. The disk is inclined at 70° to cover a wide range of scattering angles, facilitating in-depth comparisons. We first extract the SPF from an extraction zone centered at 110 au, midway across the ring, and we increase zone width from 2 to 60 au. Then we fix the width to 2 au and move its position from 81 to 134 au.

First, SPFs extracted for different zone widths are all very similar (see Fig. B.1, top), showing that the width does not influence the extraction method. However, it should be noted that we are dealing with simulations, and the dust grain properties are the same for the whole ring, which may not be true for a real disk, where a narrower zone will be preferred.

Considering now various positions, the SPFs are very similar once again, although there is a noticeable difference when the zone is close to the inner ring edge (see Fig. B.1, bottom). This suggests that the extraction zone may include areas outside the ring. Surprisingly, this is not the case for the SPF on the outer ring edge. An extraction zone far enough from the ring edges is the optimal choice for good extraction of the disk intrinsic SPF.

When extracting the SPF from the disk surface, we assume that surface brightness only depends on scattering angles, neglecting the radial effect linked to the natural brightness decrease proportional to 1/r². In an extended disk, a larger radius leads to a lower surface brightness. However, the resulting error is uniform in each direction and can be globally canceled out, especially for symmetrical disks.

Fig. B.1 Tests of SPF extraction from a 70° inclined simulated disk surface with different extraction zones. Top: The extraction zone is centered on the ring and its width increases. Bottom: The extraction zone is 2 au wide and its radial position increases. Extraction zones are shown as sets of two colored ellipses on the disk image, and corresponding SPF are plotted with the same color. We show the limb brightening corrected SPF in full color, and the uncorrected one in faint color. The theoretical intrinsic dust SPF is displayed in a dotted black line.

Appendix C Impact of parameters estimation error on extracted SPF

Errors in the disk geometry parameters affect the accuracy of SPF extraction from the disk surface. In this appendix, we investigate the impact of estimation errors on the extracted SPF, for inclination, position angle, and scattering surface height. Note that the SPFs presented here are all limb brightening corrected, also dependent on geometric estimation. So the study generalizes the impact of estimation error on SPF extraction and correction. We use the same standard model as in section 2.2, inclined at 70°and with a position angle of 270°. The SPFs are extracted by defining extraction zones with different inclinations, position angles, and aspect ratios. We measure the relative error with the reference SPF.

All three panels in figure C.1 show that the SPFs shape is always well recovered, even for errors of 4° in inclination, ± 4° in position angle, and 0.04 in aspect ratio, especially the 120° trough. In particular, an error in inclination results in a shift of the low scattering angle part of SPF towards the higher scattering angles, although the percentage error does not exceed 20% for an error of 2° in inclination. The PA primarily alters the amplitude of the SPF peak, and a 2° error results in an 11% error on the SPF. Finally, the aspect ratio is probably the most sensitive parameter, judging by the SPF extracted. Indeed, the small scattering angle part is strongly impacted by an error on h/r, as well as peak amplitude. An error of a few hundredths gives a percentage error of around 15% on the SPF. The largest percentage error is induced by inclination errors, but it should be noted that the computation of the relative error is biased by the shift of the low scattering angle part towards the higher angles. In addition, the shape of the SPF extracted for different inclinations remains very similar to each other.

Fig. C.1 Tests of the SPF extraction from a 70° inclined simulated disk surface with different offsets of ±4° on inclination (left), ±4° on position angle (middle), and ±0.04 aspect ratio (right). SPFs for underestimated and overestimated parameters are shown in blue, and red, respectively.

Appendix D Impact of dust mass on scattering surface height and pressure scale height

In this section, we investigate the influence of the disk mass on the position of the measured surface with respect to the optical depth τ_s = 1. (see Fig. D.1). These tests show that the surface measured by DRAGyS is very close to the τ_s = 1 surface for high disk dust masses (> 10⁻⁴ M_⊙ per ring). We note that the surface estimated by DRAGyS also moves progressively below the τ_s = 1 for the lower disk masses (here 10⁻⁶ and 10⁻⁷ M_⊙ in each ring). As expected, the τ_s = 1 surface moves closer to the disk midplane for decreasing disk masses and, accordingly, the ratio between the τ_s = 1 surface and the pressure scale height also decreases with decreasing disk mass. This ratio also decreases with radius because the pressure scale height increases (power-law index > 1) but the volume density in each ring decreases (because each ring has the same mass but increasing volumes with radius).

Fig. D.1 Scattering surface estimation hs using DRAGyS (orange cross), compared to the τs = 1 surface from the star (red curves) for different dust masses from Md = 10−2, 10−4, 10−6 and 10−7 M⊙ for each ring (respectively between 60-90 au, 110-140 au, and 160-190 au). As in section 3.2, we add the dust density map in the background, and we overplot the input gas pressure scale height as a gray dashed curve. The ratio (τs = 1 surface / pressure scale height) is indicated in white for each ring.

Appendix E SPF extraction for different dust properties

In this section, we investigate the impact of dust composition and grain size distribution on the extracted and limb brightening-corrected SPF. In section 3.3, we showed that applying limb brightening correction yields a remarkable agreement between the extracted SPF and the intrinsic theoretical SPF. However, multiple scattering has the potential to affect the extracted SPF even though this was not the case of silicates with a size distribution ranging from 0.005 to 3 μm (see figure 4). We therefore performed a study of SPFs extracted from the surface of disks composed of silicate, amorphous carbon, graphite and ice, with maximum grain sizes of 1, 3, 5 and 10 μm (figure E.1).

As we found with our standard model, the correction of the limb brightening effect is a necessity to get closer to the intrinsic SPF. However, it also appears that some SPFs are harder to recover, particularly at small scattering angles, such as for Amorphous carbon, with a_max = 10μm. This is likely due to the effect of multiple scattering, which tends to depolarize the scattered flux, especially at small scattering angles (Tazaki et al. 2023).

In conclusion, even if the SPF extracted from our standard model (see Section 3) is in particularly good agreement with intrinsic one, the method we applied generally reproduces well the intrinsic SPFs, although some effects due to multiple scattering can introduce biases at small scattering angles.

Fig. E.1 SPFs extracted from simulated disk images using MCFOST, changing the composition (rows) and the size distribution (columns). For each case, we display the intrinsic SPF as the dotted purple curve, and the extracted SPF before and after limb brightening correction as gray and black lines, respectively.

Appendix F 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. Figures F.1, … F.9 correspond to the datasets discussed in section 4.

Fig. F.1 Qphi image (left) with extraction zone (middle) and extracted polarized SPF (right) for HD 163296 disk in H-Band. We plot SPF for the entire disk (solid black), east side (dotted red), west side (dotted blue), and without correcting for limb brightening effect (dotted grey).

Fig. F.2 Same as Fig F.1 for the LkCa 15 disk in the J band

Fig. F.3 Same as Fig F.1 for the PDS 66 disk in the H band

Fig. F.4 Same as Fig F.1 for the PDS 66 disk in the J band

Fig. F.5 Same as Fig F.1 for the RX J1615 disk in the H band

Fig. F.6 Same as Fig F.1 for the RX J1615 disk in the J band

Fig. F.7 Same as Fig F.1 for the RX J1852 disk in the H band

Fig. F.8 Same as Fig F.1 for the V4046 Sgr disk in the H band

Fig. F.9 Same as Fig F.1 for the V4046 Sgr disk in the J band

References

All Tables

All Figures

Fig. 1 Illustration of the DRAGyS method, illustrating the ellipse estimation of the disk (red), and the geometrical method of recovering the different structure parameters. Left: image seen by the observer, with the ellipse in red and different estimated quantities, such as the semiminor and semimajor axes, m and M, the position angle, PA, and the projected distance between the star and the major axis, D, used to compute the aspect ratio. Right: sketch of protoplanetary disk observation. The point of view of observer 1 highlights the geometrical process to recover the disk geometry, namely the scattering surface height, h, at a radius, rmid, and the inclination, i. The point of view of observer 2 illustrates the definition of the scattering angle.

In the text

Fig. 2 Illustration of the elliptical fit process on a Qphi simulated disk image with multiple rings. Left: disk image with a radial cut at one azimuthal angle in white. Points used for ellipse fitting are shown in green and unwanted deleted points in smaller orange dots. Middle: radial profile extracted along the white line, where the purple curve corresponds to the smoothed profile. The maxima corresponding to the selected ring and the undesired maxima from other rings are shown in green and red, respectively. Right: disk image with the fit ellipse shown in white. A set of random points used to estimate uncertainties is plotted in blue.

In the text

Fig. 3 Comparison of estimated scattering height estimation with the τs = 1 surface integrated from the star (red) for each ring (respectively between 60–90 au, 110–140 au, and 160–190 au). We overplot the scattering height surface estimated using DRAGyS, for each inclination and position angle, in different colors and markers. We have added the dust density map in the background and, we also plot the gas pressure scale height with a dashed gray curve, using the input MCFOST parameters.

In the text

Fig. 4 Typical case of SPF extraction for simulated data, using a disk with PA = 270° and i = 70°. Left: disk image in polarized intensity, with an extraction zone defined by the green-shaded annulus. Right: polarized SPF normalized at 90° directly extracted from the disk surface (gray) and corrected from the limb brightening effect (black). The theoretical intrinsic SPF is plotted with a dashed purple line.

In the text

	Fig. 5 Limb brightening effect for different disk geometries. The forward scattering part is more affected than the backward scattering.
In the text

	Fig. 6 Qϕ images of six protoplanetary disks observed in J or H bands using DPI with the VLT/SPHERE IRDIS instrument. The maximum pixel positions for parameter estimation are shown as white dots, and the fit ellipse with a dashed black line.
In the text

Fig. 7 Polarized SPFs extracted from SPHERE Qphi observations (Ginski et al. 2023) of protoplanetary disks. The left and right panels are for H and J band observations, respectively. Solid and dashed colored lines are, respectively, SPFs uncorrected and corrected for the limb brightening effect, extracted using DRAGyS with fit geometric parameters. Dotted black lines show the reference SPFs from Diskmap (Stolker et al. 2016) following the geometric parameters from Ginski et al. (2023). Limb brightening uncorrected and Diskmap SPFs are normalized to 90°. We normalized the corrected SPF to the minimum scattering angle in common with the uncorrected one, emphasizing the slope change induced by the limb brightening correction. An offset has been added between each curve for better rendering. The percentages of the mean offset between SPF before and after removing the limb brightening are shown in color, and the ones between DRAGyS (without LB correction) and Diskmap are in black.

In the text

	Fig. A.1 Panel of the twelve simulated disk from MCFOST. All disks have the same structure with 3 rings, but for three inclinations (columns) and four position angle (rows).
In the text

Fig. B.1 Tests of SPF extraction from a 70° inclined simulated disk surface with different extraction zones. Top: The extraction zone is centered on the ring and its width increases. Bottom: The extraction zone is 2 au wide and its radial position increases. Extraction zones are shown as sets of two colored ellipses on the disk image, and corresponding SPF are plotted with the same color. We show the limb brightening corrected SPF in full color, and the uncorrected one in faint color. The theoretical intrinsic dust SPF is displayed in a dotted black line.

In the text

	Fig. C.1 Tests of the SPF extraction from a 70° inclined simulated disk surface with different offsets of ±4° on inclination (left), ±4° on position angle (middle), and ±0.04 aspect ratio (right). SPFs for underestimated and overestimated parameters are shown in blue, and red, respectively.
In the text

Fig. D.1 Scattering surface estimation hs using DRAGyS (orange cross), compared to the τs = 1 surface from the star (red curves) for different dust masses from Md = 10−2, 10−4, 10−6 and 10−7 M⊙ for each ring (respectively between 60-90 au, 110-140 au, and 160-190 au). As in section 3.2, we add the dust density map in the background, and we overplot the input gas pressure scale height as a gray dashed curve. The ratio (τs = 1 surface / pressure scale height) is indicated in white for each ring.

In the text

	Fig. E.1 SPFs extracted from simulated disk images using MCFOST, changing the composition (rows) and the size distribution (columns). For each case, we display the intrinsic SPF as the dotted purple curve, and the extracted SPF before and after limb brightening correction as gray and black lines, respectively.
In the text

	Fig. F.1 Qphi image (left) with extraction zone (middle) and extracted polarized SPF (right) for HD 163296 disk in H-Band. We plot SPF for the entire disk (solid black), east side (dotted red), west side (dotted blue), and without correcting for limb brightening effect (dotted grey).
In the text

	Fig. F.2 Same as Fig F.1 for the LkCa 15 disk in the J band
In the text

	Fig. F.3 Same as Fig F.1 for the PDS 66 disk in the H band
In the text

	Fig. F.4 Same as Fig F.1 for the PDS 66 disk in the J band
In the text

	Fig. F.5 Same as Fig F.1 for the RX J1615 disk in the H band
In the text

	Fig. F.6 Same as Fig F.1 for the RX J1615 disk in the J band
In the text

	Fig. F.7 Same as Fig F.1 for the RX J1852 disk in the H band
In the text

	Fig. F.8 Same as Fig F.1 for the V4046 Sgr disk in the H band
In the text

	Fig. F.9 Same as Fig F.1 for the V4046 Sgr disk in the J band
In the text