© The Authors 2026

1 Introduction

Circumstellar disks are the birthing grounds of planets. Observations with a high spatial resolution taken in the past decade have revealed complex disk morphologies, such as multi-ringed disks, spirals, and gaps, which are evidence of ongoing dynamical processes. Understanding the physical geometry of these disks is essential for constraining the disk evolution and the processes that lead to planet formation. However, it still remains unclear how the current observed planet population relates to the disk properties. In this context, understanding the vertical height structure of the dust is a necessary step towards linking observed disk features to the mechanisms that shape planetary systems. While several studies have measured the dust vertical height in individual disks or small samples (de Boer et al. 2016; Ginski et al. 2016; Avenhaus et al. 2018; Rich et al. 2021; Ginski et al. 2024), a large-scale analysis is needed to identify broader trends.

The first set of imaged protoplanetary disks mostly consisted of Herbig Ae/Be stars (e.g. Grady et al. 1999; Fukagawa et al. 2006), which are traditionally divided into two groups: Group I sources, whose spectral energy distributions show strong far-infrared excesses that are interpreted as flared disks with large inner cavities, and Group II sources, which display weaker far-infrared emission that corresponds to flatter disks with small or no inner cavities (e.g. Garufi et al. 2014; Honda et al. 2015). These two groups led to the long-standing study of the correlation between a disk cavity and its effect on the disk height. Group II disks have been found to be inherently faint in scattered light (Grady et al. 2005), and thus, any large sample is systematically biased towards their exclusion, making a comparison in scattered light difficult. However, by investigating the dust vertical height structure for different categories of disk morphology, it is possible to identify trends that may indicate the presence or absence of a cavity.

Further, the disk height may be affected by other properties of the disk. The degree of disk flaring is expected to correlate with the total disk mass. The disk becomes more optically thick with an increase in total dust mass, and observations of the optically thick scattering layer therefore trace higher layers within the disk. Therefore, a relation between dust mass and the measured aspect ratio is physically motivated. Over time, the disk aspect ratio is expected to decrease as the system evolves. This reflects the physical processes that cause the disk to become geometrically thinner, such as radiative cooling, mass depletion due to accretion or expulsion from the system, and reduced vertical support from a lack of gas pressure, resulting in a smaller dust vertical height relative to the radius. If this is correct, the aspect ratio should decline with increasing stellar age. Alternatively, the clearing of the cavity observed in older (transition) disks might enhance the irradiation of the outer disk, heating it more strongly and thereby increasing its flaring. Lastly, the vertical height of the dust may correlate with stellar mass (see Appendix G for the dependence derivation), since the classical mass-luminosity relation (MLR) established for main-sequence stars is not expected to hold during the pre-main-sequence phase (e.g. Baraffe et al. 2002) characteristic of protoplanetary disks. Further, Eker et al. (2018) proposed that the MLR for main-sequence stars might be described by $L\propto M_{*}^{\beta }$, with differing values for ß depending on stellar mass, reflecting transitions in the dominant energy transport mechanism. Therefore, vertical height and stellar mass might be correlated.

We present a large-sample analysis of the vertical height structure of the optically thick dust-scattering surface (τ = 1) within 92 protoplanetary disks (examples shown in Fig. 1), comprised of polarimetrie near-infrared (NIR) scattered-light observations. This sample represents a substantial fraction of all NIR scattered-light images of protoplanetary disks currently available.

This paper is organised as follows. In Sect. details of the dataset used throughout this work are presented. Section 3 outlines the methods used to extract the disk vertical height structure. Finally, Sect. 4 contains the results and discussion of the full dataset, two example individual systems, multi-ringed systems, an investigation into system properties that may affect the dust vertical height structure, and the calculation of planet masses opening gaps, a use case for the constraint of the dust vertical height structure.

Fig. 1 Gallery of representative disks from the sample, grouped by morphological geometry. The ellipse fits we made to the scattered-light disk structures are overplotted. The scale bar in each panel denotes 1 arcsecond.

Fig. 2 Distribution of year and filter within the initial dataset (totalling 294 images), showing a greater tendency towards H-band images.

2 Observations

The dataset we used comprises the full set of protoplanetary disk observations available in the European Southern Observatory (ESO) archive from the Very Large Telescope using the Spectro-Polarimetric High-contrast Exoplanet REsearch instrument (VLT/SPHERE; Beuzit et al. 2019) as of September, 2024, representing a substantial fraction of all existing polarised scattered-light images of protoplanetary disks to date. The dataset consists of 186 unique disks with 294 total images, observed with VLT/SPHERE using the Infra-Red Dual-beam Imager and Spectrograph (IRDIS; Dohlen et al. 2008) subsystem. The dual-beam polarimetric imaging (DPI; Langlois et al. 2014) mode was used in each case, measuring linear polarised light corresponding to the scattered light from the dust-scattering surface. The dataset comprised J- (1.1–1.4 μm), H- (1.5–1.8 μm), and Ks- (2.0–2.4 μm; hereafter K-band) band images taken between 2015 and 2024 (see Fig. 2). All data were reduced using IRDIS Data reduction for Accurate Polarimetry (IRDAP), which is described in detail in van Holstein et al. (2020). The reduction included a model-based correction for instrumental polarisation and crosstalk between optical components.

To differentiate signal strength in the images and to evaluate the disks' suitability for the structure extraction/ellipse fitting (SEEF) algorithm (see Sect. 3), they were classified as high-confidence, moderate-confidence or low-confidence cases. High-confidence disks required little to no user intervention, moderate-confidence disks presented challenges such as asymmetries or spiral features, and low-confidence disks were dominated by noise or exhibited insufficient disk signal for reliable detection. 71 disks were categorised as high-confidence, 21 as moderate-confidence, and 92 as low-confidence. Thus, the SEEF algorithm was run for 92 disks total (high- and moderate-confidence disks), using 1 image per disk (as many disks had observations in two or more filters). The image with the least distortions within the field-of-view was chosen.

Fig. 3 Schematic illustration of the projected geometry of an inclined multi-ringed circum-stellar disk (this also holds for continuous disks). When viewed at an inclination angle i, the disk's intrinsically (and assumed) circular structure appears elliptical due to projection. The offset u between the star position and the centre of the ellipse arises from the projected height of the disk surface along the semiminor axis. This geometric offset becomes observable in scattered-light images and is used to infer the disk height structure.

3 Methods

Following de Boer et al. (2016), to extract the vertical height profile of the optically thick dust-scattering surface of the disk, ellipses are fit to the observed structures to measure the offset, u, between the centre of the fitted ellipse and the stellar position along the minor axis. Shown in Fig. 3, the height of the scattering surface of the dust, where the disk becomes optically thick to the wavelength of light (h_τ=1, hereafter simply h), is calculated as h = u/ sin(i), where i is the inclination of the disk. The separation, inclination, and position angle (PA) are also derived from the ellipse fit.

The algorithm¹ created to fit ellipses and measure the vertical height structure consists of two main stages: structure extraction and ellipse fitting (SEEF). In the first stage, positional information of the disk structure is extracted from the image (Sect. 3.2). In the second stage, elliptical models are fitted to this information to characterise the geometry (Sect. 3.3).

3.1 Pre-processing

Pre-processing of the images was often required to increase the algorithm capability of structure extraction. Firstly, the image is rotated to an optimal starting position for outlier removal (see Appendix C). Optionally, a high pass filter, specifically unsharp masking (Malin 1977), may be applied to enhance structures. A smoothed representation was generated by convolving with a Gaussian kernel (σ = 5 pixels). The high-frequency component was obtained by subtracting the smoothed image from the original. This residual was then added back to the original image with a gain factor of 2. This step was used only for disks with high S/N and relatively smooth edges (e.g. high-confidence images with edges of S/N > 5). For disks containing faint features, unsharp masking can alter the apparent structure due to the low contrast between signal and noise. The high-pass filter was therefore used primarily to aid edge detection, but not for Gaussian fitting, as it often disrupted the intrinsic shape of the signal profile.

The second optional pre-processing tool is an interactive ellipse-based masking tool, developed in Python² to isolate specific regions of interest, specifically, for use during structure extraction in this work. The tool allows the user to define elliptical regions directly onto images using mouse-driven selection, with support for three masking modes: masking the interior of a single ellipse (inner), the exterior (outer), or the annular region between two ellipses (ring). Users can toggle between modes, zoom in and out, interactively adjust the position and size of ellipses, and apply or reset masks in real time. Once finalised, the tool outputs the masked image to the SEEF algorithm, effectively enabling precise control over which regions of the disk are included or excluded in subsequent measurements. The tool was used for structures within multi-ringed disks where signal from adjacent rings interfered with extracting the target structure or for images containing artifacts that required masking.

Fig. 4 Left: idealistic example of the edge detection method using a synthetic circumstellar disk with added Gaussian background noise. The radial cuts (dashed white lines) sample intensity profiles iteratively. Magenta dots mark detected disk edges where intensity first (measured inwards towards the star) exceeds a 5σ noise primary threshold and passes the step function criteria. Edge labels are shown only on the right side of the disk (angles between 0° and 180° clockwise) for clarity. Background noise level (σ) and primary threshold (5σ) are annotated in the figure. Right: radial intensity profile extracted at 180° from the centre of GM Aur in H-band, comparing the real high-S/N case (blue) and a synthetically degraded low-S/N case (orange). Horizontal dashed lines indicate the 5σ and 3σ primary thresholds used for edge detection. Step points used to confirm edge detection are marked with red crosses, and the final detected edge locations are indicated by vertical dashed-dotted lines and filled circles. The low-S/N profile shows degraded detection performance and premature truncation due to the outer edge being buried in noise.

3.2 Structure extraction

Two structure-extraction methods are available for use within the SEEF algorithm, edge detection and Gaussian fitting, following and building on methods described by Ginski et al. (2024). Edge detection is performed by placing a radial cut through the image to detect the outer disk edge. The most outward point along the radial cut with intensity above a user-defined threshold (hereafter called the primary threshold) is stored as a signal point (see Fig. 4). The direction of the radial cut is then rotated through the centre and the process repeated azimuthally in 1° steps. Binning is applied based on the radial extent of the feature, where features at larger radii (r > 75 px/0.92") are left unbinned to avoid undersampling. For features within approximately 100 px, Nyquist sampling is achieved, whereas extended features (r > 100 px/1.22") are undersampled, resulting in the loss of minor fine structure along the disk edge. Nonetheless, this level of sampling is sufficient for accurately extracting the disk edge location. The primary threshold is defined relative to the background standard deviation (cr), with 5σ used as a standard primary threshold for the majority of disks, with lower σ primary thresholds used for faint structures when required.

To address the issue of single-pixel high-intensity noise, a step-function was implemented. In this technique, the algorithm checks the next n pixels along the radial direction, and if their intensities also exceed the set primary threshold, the point is marked as a valid signal point. The choice of n is a user-defined parameter and dependent on the target structure: if set too high, narrow features may be overlooked or true signal points discarded due to individual noisy pixels, whereas setting it too low makes the method more susceptible to positive noise. Following Ginski et al. (2024), the radial uncertainty (Δr) is calculated as $$\Delta r=\frac{3\sigma }{m_{\text{profile}}},$$(1)

where m_profile is the slope of the step function n pixels calculated along the radial cut.

The effect of the scattering phase function (i.e. a brighter forward-scattering and fainter back-scattering side of the disk) must be considered when using the edge detection technique. Using a fixed, user-defined azimuthal threshold on inclined disks can cause the algorithm to trace different physical regions of the disk surface due to strong azimuthal intensity variations. This represents a significant methodological caveat of the approach. The impact of the scattering phase function on this method is examined in detail by Ginski et al. (2024), who apply the technique to two sets of radiative transfer models with differing degrees of flaring (α = 1.09 and α = 1.30), with both sets spanning a range of inclinations. By comparing the recovered inclination, position angle, and aspect ratio to the known model inputs, they show that the measured parameters remain consistent with the true values within the uncertainties. Notable caveats include a systematic offset of approximately 0.05 of the aspect ratio for the higher flaring case (α = 1.30), where the measured values were systematically lower than the true model values. Disks with inclinations(i ≥ 65°) were not included in this analysis and for lowly inclined disks (i ≤ 10° ), the uncertainties on PA and inclination were large, making the measurements unreliable. Further, as discussed in Appendix C, the SEEF algorithm permits the manual exclusion of azimuthal bins; for inclined disks (i ≥ 45°), bins on the back-scattering side were frequently excluded to mitigate the effects of azimuthal intensity variations introduced by the scattering phase function.

The edge detection method is well-suited for identifying the boundaries of continuous disks and faint extended structures, where the target is resolved and exhibits a distinct edge as it detects the most outward signal and the user-defined threshold can be adjusted depending on the S/N of the target, with a higher threshold level for higher S/N disks. However, for inclined disks, the method is sensitive to azimuthal intensity variations introduced by the scattering phase function, which can cause the detected edge to correspond to different physical regions of the disk surface in forward- and back-scattering directions. It is also not effective for detecting radially thin ring-like features.

Gaussian fitting was implemented as a complementary technique (shown in Fig. 5). This technique leverages the signal strength and the shape of the intensity profile as criteria for extracting the ring location. The implementation involved sampling intensity profiles along a window of a radial cut and fitting a Gaussian function of the form $$I\left(x\right)=A\exp \left(-\frac{{\left(x-\mu \right)}^2}{2{\Sigma }^2}\right),$$(2)

where A is the amplitude, μ is the centre (peak position), and Σ is the standard deviation (related to the profile width). A sliding-window technique was used to scan through each profile in fixed-size segments.

To filter meaningful detections, the algorithm enforced several constraints: the peak amplitude A had to exceed a user-defined minimum (generally a more liberal 3σ of the background), and the fitted width Σ had to fall within reasonable physical bounds, based on the width of the ring structure. For a ring structure approximately 10 pixels wide, this corresponds to Σ ≈ 4.25 pixels, since the full width at half maximum (FWHM) is related to Σ by FWHM = 2.355 • Σ. Therefore, the bounds for acceptable fits were chosen to bracket this value, typically between 2 and 7 pixels, to allow for moderate variation while excluding unrealistically narrow or broad features. This helped suppress noise and false, narrow peaks. On fulfilling this criteria, the peak of the Gaussian is stored as the detected point. Following Ginski et al. (2025), an uncertainty estimate was derived as the FWHM divided by the S/N, where the S/N was defined as A divided by the standard deviation of the background.

Unlike edge detection, which identifies the location of sharp intensity changes to define the outer edge of a ring, Gaussian fitting targets the centre of the intensity profile, effectively detecting the brightest point in the emission rather than the boundary. While this makes it conceptually different from edge detection, the resulting measurements of vertical height (see Table A.1) are functionally equivalent. Gaussian fitting also mitigates the impact of the scattering phase function compared to edge detection by providing a centroid estimate rather than relying on raw intensity. However, the stricter criteria and the additional complexity of fitting often result in fewer valid detections overall, especially in low S/N regions or when the profile deviates from a clean Gaussian shape.

Following the extraction of data points, a series of outlier removal steps were performed on the data and are explained in detail in Appendix C.

Fig. 5 Idealistic example of the Gaussian fitting. We show a synthetic circumstellar ring with added noise (left) showing the sampled radial cut (white) and the Gaussian fit window (cyan). The right panel displays the extracted intensity profile along the window and the Gaussian fit (dashed red line) modelling the ring edge, with the detected peak position indicated by the vertical green line. Note: here, σ = Σ in the main text.

3.3 Ellipse fitting

Two forms of ellipse fitting were employed: a constrained and a free version of the least-squares method.

The free ellipse fitting algorithm employed is the direct least-squares (DLS) method, introduced by Fitzgibbon et al. (1999). The method minimises the algebraic distance between the set of points derived from structure extraction and a fitted conic, subject to a quadratic constraint (4ac – b² = 1) to ensure the fitted conic is an ellipse, avoiding degenerate solutions. It reduces to solving a generalised eigenvalue problem, making it fast and numerically stable.

It is important to note that the DLS fitting algorithm inherently biases the solution toward smaller ellipses due to its use of algebraic distance as the objective function, rather than the true geometric (orthogonal) distance between data points and the ellipse. While algebraic distance offers computational efficiency and enables a closed-form solution under certain constraints, it does not uniformly penalise deviations across the ellipse and disproportionately favours solutions with smaller semi-axes. This issue has been discussed in detail by Kanatani (1994), who demonstrated that minimising algebraic distance introduces a systematic bias that is not easily corrected.

The constrained fit is built upon the SciPy least_squares function (Virtanen et al. 2020) which utilises the Levenberg-Marquardt (LM) algorithm (Moré 1978). The constrained version of the algorithm incorporates prior measurements of PA and inclination by initialising the solution with values obtained from previous ALMA observations, thereby anchoring the fit in physically motivated bounds and improving convergence in cases with noisy edge data, particularly when only a limited arc of the structure is visible. This technique enforces the conic constraint to ensure that the solution represents a valid ellipse that realistically defines the unseen structure. Although PA and inclination could be fixed, they were instead provided as initial guesses for the least-squares function rather than hard constraints. This allows the fit to remain flexible and account for uncertainties in those parameters. This form of the ellipse fit is suitable for cases with limited edge coverage (see Appendix D for further details on the coverage required for each ellipse fit form) and to address the free fit tendency towards the minimum solution. Therefore, for systems with prior available measurements of inclination and PA from ALMA-based observations, the constrained ellipse fit was always used (Table A.1 indicates which form was used for each system).

Uncertainty estimation for the fitted ellipse parameters differs depending on whether constraints are applied. In the constrained case, errors are derived analytically from the Jacobian matrix returned by the LM optimisation. For the free fit, a non-parametric bootstrap method is used. This involves repeatedly resampling the edge point data with replacement and re-fitting the ellipse in each iteration (we used 100 000 iterations). The distribution of recovered parameters across bootstrap trials provides empirical estimates of the standard deviation, which are taken as the uncertainty on each fitted parameter. Most disks are adequately sampled near the Nyquist limit; however, extended or low-contrast disk structures are often undersampled in specific azimuthal directions. As a result, uncertainty estimates from the Jacobian-based least-squares fitting and non-parametric bootstrap resampling may be dominated by well-sampled regions of the disk and may underestimate errors on global shape parameters.

The offset (used to calculate the height, h), inclination, radial separation (semi-major axis, r) and PA (see Appendix D for the PA convention we used) are measured directly from the ellipse fit. The aspect ratio (h/r) is then calculated as a measure of the disk flaring, which can be described by a power law, $$\frac{h}{r}=\frac{h_0}{r_0}{\left(\frac{r}{r_0}\right)}^{\alpha -1},$$(3)

where h₀ describes the height at radius r₀, and α is the flaring index, which quantifies the degree of flaring found within the disk.

Fig. 6 Aspect ratio (h/r) vs. separation (r) for the structures within the 92 disks, with 133 measurements (rose). A power-law model of the form $h/r=A\cdot {\left(\frac{r}{1\text{au}}\right)}^{\alpha -1}$ was fit to the data using least squares (black). (R2) is shown in the bottom right corner. For comparison, power-law fits from Avenhaus et al. (2018) and Ginski et al. (2016) are overlaid along with their corresponding data points in grey and light blue, respectively.

4 Results and discussion

The properties measured for each of the 92 disks can be found in Table A.1 and the images with ellipses overlaid in Appendix B.

A series of plots are presented within this section, showing h versus r and h/r versus r, illustrating the vertical structure of disks. The coefficient of determination, R², is used to evaluate how well the fitted models reproduce the observed data. A value of R² = 1 indicates a perfect fit, while R² = 0 implies that the model explains none of the variance in the data. Negative values of R² can occur when the model fits the data worse than a horizontal line at the mean, meaning the sum of squared residuals from the model exceeds the variance of the data. For vertical height, a power-law of the form $h=A\cdot {\left(\frac{r}{1\text{au}}\right)}^{\alpha }$ is fit to the dataset, where α is the flaring index and A is the scaling parameter, with R² used to quantify how well this curve captures the observed spread. For the aspect ratio, R² is computed using a fitted power-law model of the form $h/r=A\cdot {\left(\frac{r}{1\text{au}}\right)}^{\alpha -1}$. In both cases, higher R² values indicate stronger agreement between model predictions and data, serving as a comparative metric.

The R² value for aspect ratio versus radius is lower than that for height versus radius because dividing by r suppresses large-scale radial trends and reduces the variance of the dependent variable, inherently limiting the fraction of variance that can be explained by a power-law model in radius. However, by removing global radial scaling, h/r suppresses the dominance of large-radius points in unnormalised fits. Since this analysis compares many different disks, parameters A and α from h/r fits are adopted as less biased, more consistent measures of disk flaring. Accordingly, the R² from aspect ratio space is reported as the meaningful measure of conformity.

4.1 The full sample

To evaluate trends across the dataset, the full sample of disks is analysed collectively. By combining the measured vertical height structures from all systems, patterns arising from disk properties or morphology can be identified.

The plot of aspect ratio vs. the separation is shown in Fig. 6. The resulting value of R² = 0.158 indicates that a single power-law model explains only a small fraction of the variance in the observed aspect ratios. This suggests a weak correlation between h/r and radius across the sample, and implies that disk flaring may not follow a homogenous power-law relation for all disk morphologies. Given the low overall correlation, reflected in the weak R² value of 0.158, a single flaring (this plot yields α = 1.31) relation for the full-disk sample may not adequately capture the diversity in disk structures. The disk flaring is compared against two previous studies with height measurements; Ginski et al. (2016) measured the rings and gaps of a single multi-ringed disk HD 97048, while Avenhaus et al. (2018) measured the heights of structures within five disks around T Tauri stars. Our measured trend closely follows that of Avenhaus et al. (2018) who find a flaring index α = 1.22, but with a lower overall mean height. Compared to Ginski et al. (2016), we found a stronger flaring at small separations (r ≲ 60 au), while their larger flaring index (α = 1.73) leads to greater flaring at larger separations.

To investigate whether underlying trends exist within subsets of the population, the disks were classified into distinct morphological categories based on observed structure. The categorisation morphological criteria was as follows: those with a total radius exceeding 150 au were classified as extended; disks with an inclination greater than 70° were labelled as inclined; disks with an inclination lower than 15° were labelled as face-on; systems exhibiting visible spiral arms were identified as spiral disks; and those with a radius smaller than 50 au were designated as compact. Disks not meeting the criteria for the specific groups were assigned to the uniform category. These disks exhibit a mostly smooth, axisymmetric morphology with no distinct features that would allow further classification.

Three groups that posed a challenge to the fitting routine due to their morphology and the effect of the scattering phase-function as detailed in Sect. 3.2 were identified: spiral, inclined and face-on. This is primarily due to their distinct structural characteristics and the limitations of the fitting procedure. Spiral disks often exhibit non-axisymmetric features, such as arms or warps, that deviate from the smooth elliptical morphology of a typical protoplanetary disk, introducing challenges that can compromise structure detection. Highly and lowly inclined disks both show a common sensitivity issue: any slight deviation in the fitted ellipse parameters leads to a disproportionately large change in the inferred offset, owing to the strong dependence of the offset term on the inclination angle. For highly inclined systems, this effect is amplified because projection effects stretch the disk along the major axis, making the semi-minor axis, and thus the height measurement, extremely sensitive to even small fitting errors. For lowly inclined systems, the offset is intrinsically small because sin i ≈ 0. The offset is therefore nearly undetectable in the image, and any deviation in the measurement will cause a large relative change in the height measurement.

Excluding the aforementioned categories and plotting the remaining groups together as a restricted dataset (totalling 79 unique disks and 107 measurements), namely the extended, compact, and uniform disks, is shown in Fig. 7. The restriction results in a negligible increase in the coefficient of determination (R² = 0.166). This increase suggests that even when the difficult-to-fit disks are removed, the scatter is still present. An R² value of 0.166 implies that 16.6% of the variance in the observed aspect ratios can be explained by the model and while this is still a relatively weak correlation, it indicates that the power-law trend is not entirely absent but that a single flaring prescription is insufficient to capture the full diversity of vertical structures in the entire population.

When plotting the same graphs for extended disks only (see Fig. 8), R² increases significantly across both representations ($R_{\text{height}}^2=0.938$ and $R_{\text{h/r}}^2=0.732$). This stronger correlation indicates that extended disks more closely follow the expected power-law flaring behaviour, consistent with theoretical models (Chiang & Goldreich 1997; Kenyon & Hartmann 1987). This subset thus provides the clearest empirical support for a flared disk geometry which may be confidently described by a single power-law, among all morphological classes analysed.

The power-law fit to the extended disk category gives a flaring index of 1.70 with A = 0.0043, which is in close agreement with the flaring index of 1.73 reported by Ginski et al. (2016), though the scaling factor is slightly lower than their value of A = 0.0064. This suggests that while the overall radial dependence of the vertical height in the extended disk subsample is consistent with previous findings, the vertical heights are systematically smaller. This discrepancy in A could reflect differences in the sample selection (this work includes 11 unique disks, whereas Ginski et al. (2016) analysed only a single disk) or fitting method. Nonetheless, the agreement in the flaring exponent reinforces the presence of a broadly similar underlying disk geometry. Compared to the full sample shown in Fig. 6, which follows a trend similar to Avenhaus et al. (2018), the extended disk sample trend more closely resembles the results of Ginski et al. (2016), whose measurements were based a single extended disk. Considering Avenhaus et al. (2018) had a more varied dataset (containing disks that are classified as extended and uniform), this may have introduced a similar trend to that found in Fig. 6.

The vertical scale height of millimetre-sized dust grains has been measured in multiple ALMA studies (e.g. Villenave et al. 2020; Doi & Kataoka 2021; Jiang et al. 2025), consistently yielding values of a few au at radii equal to 100 au, corresponding to an aspect ratio of approximately 0.05 or lower. Independent JWST/MIRI observations (Villenave et al. 2024; Tazaki et al. 2025), which estimate the widths of dark lanes in edge-on disks, are consistent with these values and find similar aspect ratios. Equivalently, this work finds the aspect ratio to be approximately 0.1 at 100 au, suggesting the micron-sized dust is better coupled to the gas and are distributed at higher layers when compared to millimetre-sized dust which tends more towards the midplane of the disk. This further demonstrates that infrared observations predominantly trace the optically thick upper layers of the disk rather than the bulk dust mass.

Fig. 7 Two-panel visualisation of the disk vertical height (h) and aspect ratio (h/r) as a function of separation (r) for the restricted dataset, without spiral and highly or lowly inclined disks. Leaving extended, compact, and uniform disk types totalling 107 measurements. R2 is shown in the bottom right corner of each plot. Top: power-law fit to h vs. r shows an increasing vertical structure consistent with flaring. Bottom: h/r vs. r with power-law fits overlaid for this work (black), and comparison to Avenhaus et al. (2018) and Ginski et al. (2016) in grey and light blue respectively.

4.2 Multi-ringed systems and WIS PIT 2

This section explores the entire multi-ringed population and WISPIT 2 (TYC 5709-354-1; van Capelleveen et al. 2025), a multi-ringed disk recently resolved in the NIR for the first time.

From Fig. 10, a range of different flaring profiles is observed for multi-ringed disks, which all show an increase in h/r with separation (except for WRAY 15-788 (Bohn et al. 2019), which contains a significantly high error for the outermost ring). Excluding WRAY 15-788, IM LUP (Pinte et al. 2008) exhibits the lowest flaring index (α = 1.32) but the highest overall flaring, with HD 97048 having the lowest flaring at large separations. V4046 Sgr (Sissa et al. 2018) exhibits the steepest flaring index (α = 3.15), although this value should be interpreted with caution due to the disk's limited radial extent.

RXJ 1615 (α = 1.80; de Boer et al. 2016), TW HYA (α = 1.70; van Boekel et al. 2017), and WISPIT 2 (in H-band; α = 1.79) align with the flaring index found for extended disks (α = 1.72) and Ginski et al. (2016) for HD 97048 (α = 1.73).

An interesting comparison case for multi-ringed systems is the newly identified disk around the young (approximately 5 Myr) solar-mass star WISPIT $2\left(\text{mass=1}{.08}_{-0.17}^{+0.06}{M}_{\odot }\right)$, recently resolved for the first time as part of the WISPIT (Wide Separation Planets In Time) survey. A detailed description of the system and its discovery is presented in van Capelleveen et al. (2025). Unlike the rest of the sample, for which geometric measurements were obtained in a single filter for each disk, WISPIT 2 was fitted in the H and K bands. This allows a more direct comparison of wavelength-dependent structure (see Appendix E for the wavelength dependence analysis of the full sample), in which it is expected that longer wavelengths probe deeper scattering layers of dust within the disk. This makes WISPIT 2 a valuable case for detailed analysis, serving as a reference point for comparison against the general disk population.

The geometric fits and the resulting properties measured can be found in van Capelleveen et al. (2025). The rings exhibit a mean inclination of i = 43.99 ± 0.87° and a PA of 358.7 ± 1.1° in the H-band with i = 45.86 ± 1.27°, PA = 359.3 + 2.3° in K-band.

Since the disk was only recently discovered and no prior ALMA measurements are available, the free ellipse fit was adopted. All structures were fitted in the H-band. In the K-band image, Ring 0 (the outermost ring) was too faint to be detected.

Notably, Ring 3 exhibits a significant asymmetry, with excess signal concentrated along the minor axis in the western direction and within the forward scattering section of the disk (investigated in detail in van Capelleveen et al. 2025). This localised brightness enhancement does not appear in the other rings. A similar signal excess is observed in R1 of RXJ 1615 (see Fig. F.1) and Ring 2 of HD 97048 (Ginski et al. 2016), which are both also located along the disk minor axis. These signal excesses are shown and further explored in Appendix F. Whether these asymmetries arise from intrinsic structural variations within the disk or result from the effects of polarimetric imaging remains uncertain. In the case of WISPIT 2, the discovery of WISPIT 2b adjacent to the excess suggests that the excess may be associated with the presence of a companion. Excesses pose clear challenges for edge detection if not properly identified and addressed during fitting. The SEEF algorithm can address this issue by allowing for the exclusion of specific azimuthal sections during analysis or by increasing the primary threshold level used. Since the algorithm is not fully automated and occurs iteratively, presenting the geometric detection output directly to the user for each pass, such aberrations are readily identifiable and were manually accounted for in the fitting process of the main sample. However, for Ring 3, Gaussian fitting failed to produce an azimuthally stable fit. To address this, the edge detection method was modified to identify the peak signal within the ring region (isolated via masking) rather than the outermost edge, effectively circumventing the issue.

Figure 9 shows the heights and aspect ratios of each ring structure, plotted for both filters to assess their consistency/differences with the trend found for extended disks in Sect. 4.1. In the H band, the height and the aspect ratio show an exceptionally strong positive correlation to the extended disk trend, yielding coefficients of determination of $R_{\text{h}}^2=0.993$ and $R_{\text{h/r}}^2=0.975$, respectively. In contrast, the K band shows a moderate correlation for height ($R_{\text{h}}^2=0.501$), but the aspect ratio exhibits a weak negative correlation with radius ($R_{\text{h/r}}^2=-0.155$), suggesting that the K-band image is more affected by noise, as evident from the lower S/N and the reduced number of detected structures. Limited radial coverage weakens the reliability of any inferred geometric trends in the absence of a measurement of Ring 0. In addition, opacity effects play a significant role: the longer wavelength of the K-band allows it to probe deeper into the disk, closer to the midplane, where the vertical height is intrinsically smaller. As seen in Fig. 9, the measured heights and aspect ratios in K band are systematically lower compared to the H band, which traces higher scattering layers. This trend is clearly visible in the plots and aligns with expectations, supporting the interpretation that different wavelengths probe distinct vertical layers within the disk.

The wavelength dependence of disk height and aspect ratio is a critical factor in interpreting the results of this study. To assess this effect, the analysis was repeated using only the H-band observations (see Appendix E), which showed no significant impact on the overall sample. However, as demonstrated by the results for WISPIT 2, this dependence must be considered in individual studies that include multi-wavelength data when comparing measurements.

The evident flaring for multi-ringed disks and extended disks (with most multi-ringed disks also being extended disks) suggests that the outer rings remain illuminated and observable, rather than being shadowed by the inner disk regions. If there is a population of disks that are not flared, the outer rings could be shadowed by the inner regions and therefore might not be observable in scattered light with current limitations. This suggests that the presence of multiple visible rings provides indirect evidence for a flared geometry across these disks.

Fig. 8 Two-panel visualisation of the disk vertical height (h) and aspect ratio (h/r) as a function of separation (r) for extended disks only (with a total extent of >150 au). R2 is shown the in bottom right corner of each plot. Top: power-law fit to h vs. r shows an increasing vertical structure consistent with flaring. Bottom: h/r vs. r with power-law fits overlaid for this work (black), and comparison to Avenhaus et al. (2018) and Ginski et al. (2016) shown in grey and light blue respectively.

Fig. 9 Vertical structure of the WISPIT 2 disk. Top: h vs. r for H-band (blue) and K-band (orange). A power-law fit derived from the extended disk sample is overlaid (black). Bottom: h/r vs. r. The extended disk power law is again shown for comparison (black), and the trends found by Ginski et al. (2016), and Avenhaus et al. (2018) are also plotted in light blue and grey, respectively. The R2 values for both filter bands are shown to quantify the agreement with the extended disk trend for the height and aspect ratio. Note: in K band, Ring 0 was not detected.

Fig. 10 Aspect ratio (h/r) as a function of radius (r) for multi-ringed disks with three or more measurements only. Each disk is individually fitted with a power-law model $h/r=A{\left(\frac{r}{1\text{au}}\right)}^{\alpha -1}$, with the resulting flaring index α indicated in the legend. Overlaid dashed lines represent the trend for extended disks (black), as well as comparisons to Ginski et al. (2016) and Avenhaus et al. (2018), in light blue and grey respectively.

4.3 System properties affecting the vertical height structure

To explore whether disk flaring correlates with broader system properties, the aspect ratios derived for the full sample (see Appendix A for all measured values) were compared against dust mass, stellar age, and stellar mass. These stellar and disk properties were computed by Garufi et al. (in prep.), following the method described in Garufi et al. (2024). The spectral energy distribution and effective temperature was collected for each source through Vizier. Following this, the stellar luminosity was calculated from a de-reddened Pheonix model of the stellar photosphere (Hauschildt et al. 1999). The stellar age and mass were then constrained from isochrones (Siess et al. 2000; Bressan et al. 2012; Baraffe et al. 2015). The dust mass was directly from the observed millimetre photometry of each source with standard assumptions on optically thin emission and dust opacity (Beckwith et al. 1990). Not all disks within our dataset are represented in the analysis by Garufi et al. (in prep.). Stellar masses and ages were available for 61 out of the 92 unique disks, while dust masses were reported for 46 disks.

Ideally, all measurements should be taken at similar radial distances since h/r increases with radius, comparing measurements at different radii can introduce scatter unrelated to the investigated property. The plot shown in Fig. 11 explores the potential relation between aspect ratio and dust mass, using only measurements within a radius range between 40 and 60 au, as this bracket contained the most measurements (15 unique disks out of 46 in the total sample). However, for all brackets with a reasonable number of data points (n > 5), no correlation has been found.

This dataset includes too few data points at common radial distances (ideally within 5 au of each other) to construct a meaningful comparison. Another factor to consider is the high uncertainty in the calculation of the dust mass. A gas-to-dust mass is assumed to be the same as the interstellar medium (ISM) with a ratio of 100 (Bohlin et al. 1978). However, if the gas/dust ratio varies between disks, the aspect ratio/height profile is expected to vary also. Williams & Best (2014) found a wide range of inferred gas/dust ratios for nine class II disks with a mean ratio value of 16, with all ratios below the assumed ISM ratio. Further, the dust mass calculation assumes the dust is optically thin. In reality, a significant fraction of the disk can be optically thick at millimetre wavelengths (Garufi et al. 2025). This particularly affects compact disks as a larger fraction of the material along the line of sight is concentrated in a small area, leading to an underestimation of dust mass that is stronger for smaller disks. Thus, the results presented in Fig. 11 may contain a trend that is concealed by the effect of the opacity.

Figure 12 shows the aspect ratio vs. the stellar age for four different radial distance brackets of 10 au. Stellar age is assumed to represent the time since disk formation, with the underlying premise that the circumstellar disk forms simultaneously with its host star. The aspect ratio of protoplanetary disks is expected to decrease with stellar age as the disk cools and its dust and gas dissipate, leading to reduced pressure support and less flared geometries. Figure 12 reveals no significant correlation found between stellar age and aspect ratio for any single bracket. However, while understanding that this sample represents small number statistics and the correlation is weak, it is noted that the lowest aspect ratios are found for the youngest systems and the older systems typically have larger aspect ratios. A plausible explanation for this is the effect of self-shadowing.

Self-shadowing is a geometric effect in circumstellar disks where the inner regions block starlight from reaching the outer parts of the disk. This typically occurs when the disk has a puffed-up or vertically extended inner rim, often caused by direct irradiation, which casts a shadow over the outer disk (Natta et al. 2001; Dullemond et al. 2001). As a result, the shadowed region receives lower stellar flux, leading to lower temperatures and flatter structures. Self-shadowing can significantly alter the observed brightness, thermal emission, and vertical structure of a disk, and is a key factor in explaining non-uniform illumination patterns seen in scattered-light images. Assuming as the disk evolves, through planet formation carving gaps or other effects such as photoevaporation or accretion, the disk begins to carve an inner cavity (i.e. a transition disk), the disk would become less self-shadowed and the temperature of the outer regions of the disk would increase leading to a greater aspect ratio for older systems. Garufi et al. (2018) found that younger disks (≤5 Myr) are systematically fainter than older disks in scattered light and Garufi et al. (2022) did not find a significant trend but noted that young sources (≤2 Myr) were averagely fainter than older sources. Should this increase in brightness in older systems be caused by the reduction of self-shadowing, this would lead to greater aspect ratios in older disks in agreement with the weak correlation between aspect ratio and stellar age found here, despite the relation of a decreasing dust mass with age (Garufi et al. 2018).

Considering heights were measured at the most outward radial extent with the edge detection method and this method was applied to all continuous disks not containing gaps (i.e. the majority of disks with $\text{r\hspace{0.17em}}\underset{˜}{\underset{˜}{\text{<}}}150\text{\hspace{0.17em}au}$), if self-shadowing plays a critical role in the aspect ratio scatter observed in Sect. 4.1, then this could explain the large scatter found at these radial extents. As Garufi et al. (2024) note, the inner rim is expected to lie within 0.02–2 au of the star; therefore, the measured feature cannot correspond to the inner rim. However, if a moderate fraction of these disks are partially self-shadowed, this effect could explain the wide spread in measured aspect ratios at radii $\underset{˜}{\underset{˜}{\text{<}}}150\text{\hspace{0.17em}au}$. Extended disks do not exhibit this scatter, suggesting that if self-shadowing is indeed responsible, it is not present in these systems. This is consistent with the interpretation that such disks have cleared their inner cavities (i.e. Group I or transition disks), reducing shadowing and allowing for more uniform flaring, in agreement with theoretical expectations (Natta et al. 2001). These systems may have evolved beyond the self-shadowed phase (i.e. cleared the inner region), or never experienced it to a significant degree. Conversely, more compact disks may be preferentially affected by self-shadowing, suppressing/increasing the aspect ratio significantly because their intrinsic size is smaller, and this introduces the observed scatter at small separations (r < 50 au). Further investigation of the role of self-shadowing would confirm whether this is indeed a factor within the observed scatter in vertical height between the disks in this dataset.

The same analysis was repeated for stellar mass which is expected to correlate with the vertical height (see Appendix G), assuming that the classical MLR does not hold for pre main-sequence stars. However, similar to the stellar age, no correlation was found for aspect ratio vs. stellar mass.

For all of these previous analyses, a key limitation is the need to measure aspect ratios at a consistent radius across systems, which drastically reduces the sample size for each analysis, leading to small-number statistics. Additionally, if varying disk morphologies exhibit different aspect ratios, this would complicate any simple relation between the dust mass, stellar age or stellar mass and the disk geometry.

Fig. 11 Aspect ratio (h/r) vs. disk dust mass (in Earth masses) for all measurements within a separation bracket of 40–60 au.

Fig. 12 Aspect ratio (h/r) plotted against stellar age (Myr) for four separation brackets to examine how the disk vertical structure evolves with age.

4.4 Planets opening gaps

The formation of gaps in protoplanetary disks by embedded planets is a well-established outcome of planet-disk interactions, driven by the exchange of angular momentum between the planet and the disk gas (Lin & Papaloizou 1979; Goldreich & Tremaine 1980; Lin & Papaloizou 1993). The relation between planet mass and gap width has been extensively explored through hydrodynamic simulations (Kanagawa et al. 2016; Dong & Fung 2017; Zhang et al. 2018). Because the gap width and aspect ratio are measured through ellipse fitting and are key parameters in planet-disk interaction models, the results of the SEEF algorithm provide a natural foundation for exploring the theoretical masses of planets responsible for opening these gaps.

We applied the empirical model developed by Kanagawa et al. (2016) to estimate the planet masses from observed gap widths, assuming the gap is opened by a single planet located at the gap centre, $$\frac{M_{\text{p}}}{M_{*}}=2.1\times {10}^{-3}{\left(\frac{{\Delta }_{\text{gap}}}{R_{\text{p}}}\right)}^2{\left(\frac{h_p}{0.05R_{\text{p}}}\right)}^{\frac{3}{2}}{\left(\frac{{\alpha }_{\text{visc}}}{{10}^{-3}}\right)}^{\frac{1}{2}},$$(4)

where M_p,R_p, and h_p/R_p are the planet mass, orbital radius, and local disk aspect ratio of the gas pressure scale height, respectively; M_* is the stellar mass; and α_visc is the disk viscosity parameter (Shakura & Sunyaev 1973). Following Ginski et al. (2016); Ginski et al. (2025), the gas pressure scale height is estimated from the dust-scattering surface height by dividing by a factor of 4, under the assumption that gas and dust are well mixed (Chiang et al. 2001). We note that the conversion between the scattering-surface height and the gas pressure scale height (here assumed to be a factor of four) is uncertain and depends on the dust optical depth and vertical dust distribution, introducing a systematic uncertainty in the inferred gas scale height, and consequently in the estimated planet masses. For each gap, the table lists the gap width (Δ_gap), the corresponding planet's semi-major axis (R_p), and the local disk aspect ratio, derived from observations. The gap width was calculated as R₀ – R₁, where R₀ and R₁ are the radii of the adjacent rings, measured along the major-axis and R₀ is the most outward ring of the pair. This definition formally describes the gas gap; throughout this work, we assume that the dust is sufficiently well coupled to the gas that the measured dust ring locations trace the underlying gas structure, consistent with the assumptions adopted in Kanagawa et al. (2016). This assumption neglects the fact that gaps seen in scattered light may also be shaped by shadowing and by the geometry of the disk scattering surface. The planet is assumed to orbit at the midpoint of the gap. The aspect ratio was calculated using the scale parameter A and flaring index α derived from the power-law fits for individual systems (see Fig. 10). For systems with fewer than three height measurements, the extended disk trend was adopted instead (see Sect. 4.1). For WRAY 15-788 with three measurements, which had a large uncertainty for the outermost ring, the extended disk trend values were used. These values were then applied using the power-law relation $$\frac{h_{\text{p}}}{R_{\text{p}}}=A{\left(\frac{R_{\text{p}}}{1\text{au}}\right)}^{\alpha -1}.$$(5)

Planet masses are calculated under the assumption that the gaps are opened by single planets and are reported for three values of the disk viscosity parameter α_visc = 10⁻², 10⁻³ and 10⁻⁴ to reflect different plausible disk conditions. Uncertainties were estimated as follows: the gap width uncertainty reflects the propagated uncertainties in the positions of the adjacent rings; a 10% uncertainty was assigned to the stellar mass; and a fixed uncertainty of 0.05 was adopted for the aspect ratio, consistent with the scatter observed in the extended disk population (see Fig. 8). These estimates provide a comparative view of how inferred planet masses vary with assumed disk viscosity. Figure 13 shows the inferred planet masses for the sample of multi-ringed disks for the three values of disk viscosity.

The bottom panel of Fig. 13 presents the inferred planet masses and separations derived, overlaid on the broader population of detected exoplanets for reference. Among the models explored, the high viscosity α_visc = 10⁻² case aligns most closely with the current exoplanet detections in the mass-separation plane yet still below current detection limits. However, this apparent agreement should be interpreted with caution, as planet migration may shift planets inward from their original formation locations, meaning that the observed disk gaps may not correspond directly to current orbital positions (Goldreich & Tremaine 1980; Fernandez & Ip 1984). In the low-viscosity cases (α_visc = 10⁻³ and 10⁻⁴), the inferred planets occupy a largely unexplored region of parameter space, potentially due to significant inward migration after formation. It is noted that this region of the mass-separation plane aligns with the population of debris-inferred exoplanets using sculpting and stirring arguments (Pearce et al. 2022).

For disk viscosity parameters between 10⁻² and 10⁻⁴, the range of inferred planet is between approximately 2 M_J and 0.2 M_J for larger gaps or 0.02 M_J and 0.2 M_J for smaller gaps. Observationally, current detection limits within disk gaps are primarily constrained by achievable contrasts. Ginski et al. (2025) derived a limit of approximately 5 M_J for the gap in 2MASS J1612, while Asensio-Torres et al. (2021) reported a range of detection limits from 9 M_J at close separations (approximately 10 au) down to 4 M_Jup at wider separations (approximately 100 au), highlighting how sensitivity improves with increasing orbital radius. Considering the detection limits, then all the inferred planets would not be directly detectable by current limitations. However, the large uncertainties in the inferred planet masses (see Table H.1), often exceeding the nominal values, arise primarily from the ±0.05 uncertainty in the aspect ratio, so these results should be interpreted more as an investigation into the dependence of the inferred planet masses on the viscosity parameter than an estimate of the true planet masses.

Fig. 13 Top: planet mass (MJ) vs. separation of the inferred planet (Rp, au) for three different values of the disk viscosity parameter αvisc. The top left panel shows all αvisc values overlaid together. The remaining panels show individual plots for αvisc = 10−2, 10−3, and 10−4, respectively. The error bars represent uncertainties in the inferred planet mass and the measured gap radius. The plot highlights the dependence of inferred planet mass on αvisc, with lower viscosity requiring lower planet masses to carve equivalent gaps. Bottom: distribution of detected exo-planets plotted as mass (MJ) vs. separation (au). Blue points represent confirmed exoplanets from the exoplanet.eu database, where true masses are used when available and minimum masses (M × sin i) are used otherwise. Other coloured points correspond to the planet masses we inferred for different viscosities.

5 Summary and conclusions

We analysed the vertical height structure for a sample of 92 class II systems containing disks. This expands on the samples in the literature for vertical height measurements (de Boer et al. 2016; Ginski et al. 2016; Avenhaus et al. 2018; Rich et al. 2021; Ginski et al. 2024). We also presented a method that expands the methods introduced by Ginski et al. (2024) for robustly measuring the vertical heights and aspect ratios for differing morphologies of protoplanetary disks. The vertical heights, aspect ratios, separations, inclinations, and position angles of the structures with the 92 disks are listed in Table A.1.

For the full sample, including the investigation into removing spiral, face-on, and inclined disks and further, limiting the dataset to H-band filter images alone (see Appendix E), we found that the vertical height structure cannot be confidently described by a single power law. This is evidenced by the scatter and the moderate to low correlation observed in the h/r versus r distributions. The impact of observational filter used on the inferred flaring was found to be significant, as seen in the case of WISPIT 2, and it must be considered when comparing vertical height measurements made using different filters. However, in this work, where the majority of images were in H band, the choice of filter did not significantly affect the trends we presented. The low correlations suggest that either (1) different disk morphologies intrinsically follow distinct vertical structures, or (2) an additional, unidentified factor contributes to the scatter across the population. Notably, the categorical breakdowns and system properties (Sect. 4.3) we explored revealed no significant trends within any underlying factor or subgroup, except for the extended disks, which remain the only category that consistently follows a well-defined flaring power law.

Extended disks (r_outer ≥ 150au) were found to be well described by a single power law, $$h=0.0043{\left(\frac{r}{1\text{au}}\right)}^{1.70},$$(6)

suggesting that some underlying physical property of these systems is resilient to the factors driving the observed spread in flaring among other disk categories. This consistency implies thatextended disks mightoccupy a distinct evolutionary orstructural regime in which processes such as dust settling, turbulence, or self-shadowing are either stabilised or play a diminished role in shaping the vertical structure.

While no statistically significant correlation was found between aspect ratio and disk properties such as stellar mass, age, or dust mass, a trend was noted in which the youngest systems exhibit the lowest aspect ratios, with older systems showing higher values. If self-shadowing plays a dominant role in shaping vertical height and disk flaring, this trend would be consistent with older disks clearing their inner region, that is, Group I disks. However, the investigation of these dependences is strongly limited by small-number statistics. Disk self-shadowing might contribute to the observed scatter in vertical height, regardless of stellar mass, age, disk category, or dust mass.

Using the vertical height trends, we estimated the masses of gap-opening planets. The inferred planet masses span approximately 0.2–2 M_J for larger gaps and 0.02–0.2 M_J for smaller gaps, depending on the assumed disk viscosity parameter (α_visc). These planets remain below current direct-detection limits. Given the large uncertainties in these estimates, the derived values should therefore be regarded as tentative. With reference to the detected exoplanet population (Fig. 13) and assuming the traced disk gaps are predominantly opened by planets forming gaps, our analysis suggests two potential scenarios. If the planets remain at their current location in the mass-separation plane, then the disk structures trace an entirely new and lower-mass population of wide-separation planets compared to those previously detected. Alternatively, these planets might still be undergoing significant inward migration within the disk, in which case these may be precursors to the lower separation planet population detected with indirect methods.

Data availability

A copy of the reduced images is available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/710/A38.

Acknowledgements

We would like to thank the anonymous reviewer for their careful reading of the manuscript and for the many constructive comments that helped improve the quality of this work. This work is based on observations collected with the SPHERE instrument on the Very Large Telescope at the European Southern Observatory under multiple ESO programmes, retrieved from the ESO Science Archive Facility. The research has made use of the following software and Python libraries: NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), and Astropy (Astropy Collaboration 2013, 2018, 2022). The author acknowledges the support and supervision provided by the School of Natural Sciences and the Centre for Astronomy at the University of Galway. This paper is based on the author's Master's thesis carried out at the University of Galway.

References

Appendix A Table of results

Table A.1 details the height (au), inclination (degrees), separation (au), aspect ratio (h/r), PA (degrees), and whether the ellipse fit had prior estimates of the inclination and PA for the constrained ellipse fit (Constr?) for all measured disks. The region of the ellipse fit is noted to differentiate between substructures of the same disk, where inner (2) and (3) refer to inner rings at a greater radial extent than the most inner structure which is simply referred to as inner. (*) The PA shown is measured from the +y axis to the first instance of the major axis counter-clockwise, as the front scattering side could not be determined.

Appendix B Figure of ellipse fits for all disks

Shown below is the figure containing all the ellipse fits derived through this work, overlaid on their respective image. The ellipse fit shown is the constrained fit if available and the free fit otherwise. The disks are sorted into their respective disk morphologies defined in Sect. 4.1 for comparison.

Fig. B.1 Ellipse fits described in Sect. 3.3 overlaid on all 92 unique disks, categorised by the disk type defined in Sect. 4.1. Scattered-light images are displayed using either a logarithmic stretch or a linear scale with fixed υmin and υmax to optimise visibility of disk structure. The scale bar in each panel denotes 1 arcsecond.

Appendix C Outlier removal

Due to the faint and noisy nature of the images, the detected points are often irregular and inconsistent, with spurious detections occurring across closely spaced radial passes. Several stages of outlier removal are implemented within the algorithm to combat clear outliers, hereafter also referred to as false detections.

First, extreme outliers are removed by computing the radial distance of each detected point from the image centre (i.e. the coronagraph position) and discarding those exceeding the mean radial distance of all points by a user-defined threshold (typically ≥50 pixels), which changes depending on whether one is extracting information from a feature close (smaller threshold) or distant from the centre (larger threshold).

The points were then spatially binned based on azimuthal angle. High S/N disks with well-defined edges require little to no binning, while poorly defined or noisy edges benefit from binning to suppress spurious structure. Binning reduces the total number of data points and smooths out local irregularities, but it may interfere with other outlier removal processes, and their user-defined parameters must be adjusted to account for the decrease in points and for the increase in distance between them. Therefore, the bin angle selected varies between disks, with little/no binning for extended structures and greater degrees of binning for compact, rough-edged disks. It is noted that one must consider that using a large binning angle (≥ 5°) causes the averaged coordinates to systematically truncate inward, since the binned points approximate a curved structure, as illustrated in Fig. C.1.

Fig. C.1 Illustration of point binning with excessive angular bin size. When the bin angle is too large, true edge points (blue) are truncated, and the resulting averaged point (red) no longer accurately represents the underlying distribution.

An outlier removal step is now performed within each bin by discarding points exceeding the mean radial distance of the bin by a smaller threshold (typically 10 pixels), aiming to remove extreme outliers from their local bin. After averaging each bin, a final pass removes bins whose mean radius deviates from the total mean of all bins by more than 20-30 pixels. This helps eliminate bins comprised of entirely false detections but can be relaxed for highly eccentric structures, where the disk edge greatly varies in distance from the star azimuthally.

With global and binned outliers removed, to identify and remove local outliers that deviate from the natural curve of the structure along the extracted azimuthal profile, a two-stage filtering process was implemented, visualised entirely in Fig. C.3.

In the first stage (shown in Fig. C.2 and together with stage 2 in Fig. C.3), data points were scanned sequentially using a sliding-window technique. Each point was compared to the local mean radial distance, obtained from the trimmed mean of the previous ten points (central eight by value) to reduce outlier influence. Naturally, the first ten points were excluded from this removal step (see Fig. C.3); therefore the image must be rotated to an angle where the edge is smooth and well-defined (e.g. the forward-scattering side). Points deviating beyond the trimmed mean plus a fixed threshold (we used 10 pixels) were identified as false detections. To handle structural discontinuities (e.g. gaps or disk breaks due to shadowing or faint signal on the back-scattering side), a new section was initiated once the number of consecutive outliers exceeded a predefined 'bridge' threshold (typically 20 rejections), preventing over-rejection near features such as the noisy back-scattering side.

In the second stage of the local outlier removal (shown alongside stage 1 in Fig. C.3), each section was independently processed using 1.5σ clipping around its mean radial distance to remove residual outliers. This targets clusters of false detections that passed initial filtering, especially in low S/N regions where noise near the disk edge can bias binned averages outward. These biased bins inflate local standard deviations, reducing the effectiveness of earlier filtering. To prevent over-rejection, the first (typically largest) section on the high S/N front-scattering side is exempt, as strict clipping would wrongly exclude points near the curve ends due to natural curvature.

Outside of the previously mentioned, predominately automatic outlier removal steps, the algorithm allows for the manual exclusion of bins within specific angular ranges. These user-defined removals can be applied either before or after the outlier rejection process. When applied beforehand, they prevent problematic regions, such as regions of consistently poor signal (e.g. the back-scattering side), from influencing the initial threshold-based or rolling average calculations. When applied after outlier rejection, they allow for a targeted clean-up of residual artefacts, particularly in edge cases where structured noise or imaging artefacts remain.

Fig. C.2 Diagram illustrating the first step of local outlier removal, where points are identified as outliers if their radial distance exceeds the standard deviation of the previous ten points plus a defined threshold.

Appendix D Coverage and PA convention

To evaluate how azimuthal coverage of the extracted points from the structure extraction step of the algorithm impacts the accuracy of free ellipse fit, two configurations were tested on the same model disk, shown in Fig. D.1. In each case, edge points were incrementally added until the ellipse solution converged. For the first configuration, convergence occurred after accumulating 270° of continuous coverage on one side of the disk. In the second configuration, convergence was achieved with only 40° of total azimuthal coverage, where the points were distributed equally on opposing sides of the ring. This result demonstrates that accurate ellipse fitting depends more critically on azimuthal distribution than on total angular span. Even limited angular coverage can constrain the fit effectively, providing the points sample diverse regions of the disk geometry.

Fig. C.3 Diagram illustrating the two-step local outlier removal process. The elliptical disk is divided into angular sections based on series of outlier rejections from step 1. Section 1 (green) spans the bright, front-scattering side and undergoes step 1 filtering. Sections 2 (yellow) and 3 (purple) lie on the fainter, back-scattering side and are targeted by the second-stage rejection (black crosses). Red arcs (a bridge) mark angular gaps where the local standard deviation is inflated due to lower S/N, preventing the sliding window from identifying outliers. Defining new sections in these regions enables localised standard deviation recalculation, allowing for successful outlier detection across the back-scattering side.

The PA obtained directly from the ellipse fit is measured counter-clockwise from the positive x-axis to the first instance of the major axis. To make this consistent with the typical convention, which is illustrated in Fig. D.2 and used throughout this work. The typical convention for PA is measured from north through east to the first instance of the major axis such that the forward-scattering side lies between this axis and its opposite counterpart. To correct the PA measured from the ellipse fit, an offset of 90^° was added to convert the measurement to start from the +x axis to +y axis. If this angle did not contain the forward scattering side, then an additional offset of 180^° was added, giving offset options of 90° and 270°. The forward-scattering side was identified by comparing phase functions derived for each offset, inspecting the corresponding disk-region overlay fits, and assessing the global disk geometry and brightness distribution for each side of the disk to determine which side is the forward and back scattering region (further discussed in Coyne et al., in prep.).

Fig. D.1 Example of two free ellipse fits (cyan) applied to the same model disk, where edge points (red) were extracted at increasing azimuthal angles in one direction. Left: the detected edge points span 270° of a single side of the ring. Right: points span only 40° but are located on opposing sides of the ring. The two ellipses accurately map the underlying structure, showing the importance of extracting information from various azimuthal regions of the disk.

Fig. D.2 Diagram showing the PA convention. The PA is measured north to east, counter-clockwise, and to the major-axis containing the forward-scattering region of the disk.

When the forward-scattering side was ambiguous, a vector was drawn from the stellar position to the centre of the fitted ellipse; for an inclined disk this vector should point toward the back-scattering side, since the star appears offset toward the forward-scattering side in the image plane. In some cases, the vector lay approximately along the major axis and was therefore unreliable. For these systems, the reported PA is defined simply as the angle from north to the first occurrence of the major axis, independent of any forward-scattering determination.

Appendix E Wavelength dependence analysis

The dataset includes observations in the J (1.1 – 1.4 μm), H (1.5 – 1.8 μm), and K (2.0 – 2.4 μm) bands, each sampling slightly different heights in the disk surface. To assess the impact of this wavelength dependence on the derived aspect ratios, the analysis was repeated using only the H-band subset, with the previously identified difficult-to-fit groups removed (see Sect. 4.1), as shown in Fig. E.1 (totalling 40 unique disks and 58 measurements). Since H-band observations constitute the majority of the dataset, the resulting correlations remained mostly consistent, suggesting that band-dependent effects do not significantly influence the overall trend. The coefficient of determination R² decreased from 0.166 to 0.139, indicating a slight, but not significant, decrease in correlation when limiting the sample. This decrease contrasts with the expected improvement following dataset standardisation, although the decrease is minor.

Fig. E.1 h (Top) and h/r (Bottom) versus separation (r) for the restricted dataset, excluding spiral and lowly or highly inclined disks and now including only H-band observations. As the majority of the dataset consists of H-band images, the resulting distributions and fitted trends show minimal deviation from the restricted dataset containing all filters.

As seen from the analysis on WISPIT 2 (Sect. 4.2), wavelength dependence does indeed play a significant role for individual systems imaged in different Alters and would also influence datasets with a more evenly spread filter distribution. However, since the dataset contains predominately H-band images, the effect does not substantially contribute to the scatter observed in the h/r versus r relations presented here.

Appendix F Signal excesses found in multi-ringed systems

For three multi-ringed systems RXJ 1615, HD 97048, and WISPIT 2, an excess emission has been observed along the minor axis. It is found within the forward scattering region for RXJ 1615 and WISPIT 2 and the back scattering region for HD 97048. The excess signal for WISPIT 2 is discussed in Sect. 4.2 and in detail by van Capelleveen et al. (2025). To the best of our knowledge, there has been no previous discussion of the signal excess within RXJ 1615.

For HD 97048, Ginski et al. (2016) noted excess emission along the eastern (back-scattering) side of Ring 2, which they attributed to a partially illuminated inner wall of the structure. Their fitting technique, which targets the peak flux, is largely insensitive to this extended signal at the ring's edge. In contrast, this work employs an edge-detection method that traces the outer boundary of the emission (Gaussian fitting was not used for Ring 2 due to a lack of retained points during extraction), thereby including the excess flux in the fit. If this extended emission indeed originates from a vertical wall rather than the true disk signal, the resulting ellipse At, offset, and inferred vertical proAle may not reflect the true geometry of the ring. Rich et al. (2021) investigated this scenario by generating scattered-light radiative transfer models using TORUS (Harries et al. 2019) for single-ring and two-ring disk geometries, with radial extents comparable to the inner disk and Ring 2 of HD 97048. In the single-ring case, the ring is fully illuminated, producing a bright wall-like feature. In contrast, the two-ring conflguration (as observed in HD 97048) experiences shadowing from the inner disk, and no such wall is observed. Following this, Rich et al. (2021) Atted an ellipse to Ring 2 with h = 32.3 ± 1.5 au at r = 172 ± 4 au, which results in h/r = 0.188, a higher flaring than this work (h/r = 0.160) but signiflcantly lower than the aspect ratio reported by Ginski et al. (2016, h/r = 0.307).

Appendix G Deriving the stellar mass and vertical height correlation

Further to Sect. 4.3, we explicitly derived the expected dependence of the disk vertical scale height on stellar mass to demonstrate that if luminosity is not proportional to M⁴, the cancellation between the sound speed and Keplerian frequency no longer holds. In this case, the disk scale height acquires an explicit dependence on stellar mass through the mass-luminosity relation, such that variations in stellar luminosity at Axed mass, as expected for pre-main-sequence stars, may introduce trends in the disk aspect ratio.

The vertical scale height of a disk in hydrostatic equilibrium is given by $$h=\frac{c_s}{{\Omega }_K},$$(G.1)

where c_s is the isothermal sound speed and Ω_K is the Keplerian angular frequency. The Keplerian frequency scales as $${\Omega }_K=\sqrt{\frac{G{M}_{*}}{r^3}}\propto M_{*}^{1/2},$$(G.2)

at fixed radius. The sound speed is given by $$c_s=\sqrt{\frac{k_{B}T}{\mu m_H}}\propto T^{1/2}.$$(G.3)

Assuming the disk temperature is dominated by stellar irradiation, $$T\propto L_{*}^{1/4}.$$(G.4)

For a general mass-luminosity relation of the form $$L_{*}\propto M_{*}^{\beta },$$(G.5)

the sound speed scales as $$c{}_s\propto M_{*}^{\beta /8}.$$(G.6)

Combining these relations yields $$h\propto \frac{c_s}{{\Omega }_K}\propto M_{*}^{\left(\beta -4\right)/8}.$$(G.7)

For the classical main-sequence case where ß = 4, the stellar-mass dependence cancels exactly. For pre-main-sequence stars, where ß is not Axed, an explicit dependence of the disk scale height on stellar mass is expected.

Fig. F.1 WISPIT 2 (left) and RXJ 1615 (right) with the excess signal annotated for each image. Both signal excesses appear along the minor axis and on the forward scattering side of the disk.

Appendix H Table of planets opening gaps

Table H.1 summarises the inferred planet masses associated with observed gaps in multi-ringed disks within the sample, based on the gap-opening model of Kanagawa et al. (2016) and as discussed in Sect. 4.4. The gap width is measured as the midpoint between the two adjacent rings/disk, where the point is measured at the peak of the Gaussian or at the edge of a disk, depending on whether Gaussian fitting or edge detection was used. While the two methods measure different points within their respective structure, the signal (and distance) added or lost from the difference in the location of the points is negligible compared to the total gap width. R_p denotes the separation of the inferred planet (e.g. the midpoint of the gap), the aspect ratio was calculated using the scale parameter A and flaring index a derived from the power-law fits for individual systems in Sect. 4.2 and where measurements were not made, the extended disk trend parameters (see Sect. 4.1) were used instead. Finally, the planet masses are listed for the three considered viscosity parameters.

All Tables

All Figures

	Fig. 1 Gallery of representative disks from the sample, grouped by morphological geometry. The ellipse fits we made to the scattered-light disk structures are overplotted. The scale bar in each panel denotes 1 arcsecond.
In the text

	Fig. 2 Distribution of year and filter within the initial dataset (totalling 294 images), showing a greater tendency towards H-band images.
In the text

Fig. 3 Schematic illustration of the projected geometry of an inclined multi-ringed circum-stellar disk (this also holds for continuous disks). When viewed at an inclination angle i, the disk's intrinsically (and assumed) circular structure appears elliptical due to projection. The offset u between the star position and the centre of the ellipse arises from the projected height of the disk surface along the semiminor axis. This geometric offset becomes observable in scattered-light images and is used to infer the disk height structure.

In the text

Fig. 4 Left: idealistic example of the edge detection method using a synthetic circumstellar disk with added Gaussian background noise. The radial cuts (dashed white lines) sample intensity profiles iteratively. Magenta dots mark detected disk edges where intensity first (measured inwards towards the star) exceeds a 5σ noise primary threshold and passes the step function criteria. Edge labels are shown only on the right side of the disk (angles between 0° and 180° clockwise) for clarity. Background noise level (σ) and primary threshold (5σ) are annotated in the figure. Right: radial intensity profile extracted at 180° from the centre of GM Aur in H-band, comparing the real high-S/N case (blue) and a synthetically degraded low-S/N case (orange). Horizontal dashed lines indicate the 5σ and 3σ primary thresholds used for edge detection. Step points used to confirm edge detection are marked with red crosses, and the final detected edge locations are indicated by vertical dashed-dotted lines and filled circles. The low-S/N profile shows degraded detection performance and premature truncation due to the outer edge being buried in noise.

In the text

Fig. 5 Idealistic example of the Gaussian fitting. We show a synthetic circumstellar ring with added noise (left) showing the sampled radial cut (white) and the Gaussian fit window (cyan). The right panel displays the extracted intensity profile along the window and the Gaussian fit (dashed red line) modelling the ring edge, with the detected peak position indicated by the vertical green line. Note: here, σ = Σ in the main text.

In the text

Fig. 6 Aspect ratio (h/r) vs. separation (r) for the structures within the 92 disks, with 133 measurements (rose). A power-law model of the form $h/r=A\cdot {\left(\frac{r}{1\text{au}}\right)}^{\alpha -1}$ was fit to the data using least squares (black). (R2) is shown in the bottom right corner. For comparison, power-law fits from Avenhaus et al. (2018) and Ginski et al. (2016) are overlaid along with their corresponding data points in grey and light blue, respectively.

In the text

Fig. 7 Two-panel visualisation of the disk vertical height (h) and aspect ratio (h/r) as a function of separation (r) for the restricted dataset, without spiral and highly or lowly inclined disks. Leaving extended, compact, and uniform disk types totalling 107 measurements. R2 is shown in the bottom right corner of each plot. Top: power-law fit to h vs. r shows an increasing vertical structure consistent with flaring. Bottom: h/r vs. r with power-law fits overlaid for this work (black), and comparison to Avenhaus et al. (2018) and Ginski et al. (2016) in grey and light blue respectively.

In the text

Fig. 8 Two-panel visualisation of the disk vertical height (h) and aspect ratio (h/r) as a function of separation (r) for extended disks only (with a total extent of >150 au). R2 is shown the in bottom right corner of each plot. Top: power-law fit to h vs. r shows an increasing vertical structure consistent with flaring. Bottom: h/r vs. r with power-law fits overlaid for this work (black), and comparison to Avenhaus et al. (2018) and Ginski et al. (2016) shown in grey and light blue respectively.

In the text

Fig. 9 Vertical structure of the WISPIT 2 disk. Top: h vs. r for H-band (blue) and K-band (orange). A power-law fit derived from the extended disk sample is overlaid (black). Bottom: h/r vs. r. The extended disk power law is again shown for comparison (black), and the trends found by Ginski et al. (2016), and Avenhaus et al. (2018) are also plotted in light blue and grey, respectively. The R2 values for both filter bands are shown to quantify the agreement with the extended disk trend for the height and aspect ratio. Note: in K band, Ring 0 was not detected.

In the text

Fig. 10 Aspect ratio (h/r) as a function of radius (r) for multi-ringed disks with three or more measurements only. Each disk is individually fitted with a power-law model $h/r=A{\left(\frac{r}{1\text{au}}\right)}^{\alpha -1}$, with the resulting flaring index α indicated in the legend. Overlaid dashed lines represent the trend for extended disks (black), as well as comparisons to Ginski et al. (2016) and Avenhaus et al. (2018), in light blue and grey respectively.

In the text

	Fig. 11 Aspect ratio (h/r) vs. disk dust mass (in Earth masses) for all measurements within a separation bracket of 40–60 au.
In the text

	Fig. 12 Aspect ratio (h/r) plotted against stellar age (Myr) for four separation brackets to examine how the disk vertical structure evolves with age.
In the text

Fig. 13 Top: planet mass (MJ) vs. separation of the inferred planet (Rp, au) for three different values of the disk viscosity parameter αvisc. The top left panel shows all αvisc values overlaid together. The remaining panels show individual plots for αvisc = 10−2, 10−3, and 10−4, respectively. The error bars represent uncertainties in the inferred planet mass and the measured gap radius. The plot highlights the dependence of inferred planet mass on αvisc, with lower viscosity requiring lower planet masses to carve equivalent gaps. Bottom: distribution of detected exo-planets plotted as mass (MJ) vs. separation (au). Blue points represent confirmed exoplanets from the exoplanet.eu database, where true masses are used when available and minimum masses (M × sin i) are used otherwise. Other coloured points correspond to the planet masses we inferred for different viscosities.

In the text

	Fig. B.1 Ellipse fits described in Sect. 3.3 overlaid on all 92 unique disks, categorised by the disk type defined in Sect. 4.1. Scattered-light images are displayed using either a logarithmic stretch or a linear scale with fixed υmin and υmax to optimise visibility of disk structure. The scale bar in each panel denotes 1 arcsecond.
In the text

	Fig. C.1 Illustration of point binning with excessive angular bin size. When the bin angle is too large, true edge points (blue) are truncated, and the resulting averaged point (red) no longer accurately represents the underlying distribution.
In the text

	Fig. C.2 Diagram illustrating the first step of local outlier removal, where points are identified as outliers if their radial distance exceeds the standard deviation of the previous ten points plus a defined threshold.
In the text

Fig. C.3 Diagram illustrating the two-step local outlier removal process. The elliptical disk is divided into angular sections based on series of outlier rejections from step 1. Section 1 (green) spans the bright, front-scattering side and undergoes step 1 filtering. Sections 2 (yellow) and 3 (purple) lie on the fainter, back-scattering side and are targeted by the second-stage rejection (black crosses). Red arcs (a bridge) mark angular gaps where the local standard deviation is inflated due to lower S/N, preventing the sliding window from identifying outliers. Defining new sections in these regions enables localised standard deviation recalculation, allowing for successful outlier detection across the back-scattering side.

In the text

Fig. D.1 Example of two free ellipse fits (cyan) applied to the same model disk, where edge points (red) were extracted at increasing azimuthal angles in one direction. Left: the detected edge points span 270° of a single side of the ring. Right: points span only 40° but are located on opposing sides of the ring. The two ellipses accurately map the underlying structure, showing the importance of extracting information from various azimuthal regions of the disk.

In the text

	Fig. D.2 Diagram showing the PA convention. The PA is measured north to east, counter-clockwise, and to the major-axis containing the forward-scattering region of the disk.
In the text

	Fig. E.1 h (Top) and h/r (Bottom) versus separation (r) for the restricted dataset, excluding spiral and lowly or highly inclined disks and now including only H-band observations. As the majority of the dataset consists of H-band images, the resulting distributions and fitted trends show minimal deviation from the restricted dataset containing all filters.
In the text

	Fig. F.1 WISPIT 2 (left) and RXJ 1615 (right) with the excess signal annotated for each image. Both signal excesses appear along the minor axis and on the forward scattering side of the disk.
In the text