© The Authors 2026

1 Introduction

Debris disks offer a window into the final stages of planet formation, revealing details about the dynamical evolution of systems on timescales from tens of millions to billions of years (for reviews, see: Wyatt 2008; Matthews et al. 2014; Hughes et al. 2018). In particular, the vertical distribution of dust in debris disks contains information about the orbital distribution of planetesimals and collision rates within the disk (e.g., Thébault & Augereau 2007; Quillen et al. 2007; Moore et al. 2023). A measurement of the disk scale height at (sub)millimeter wavelengths is effectively a measurement of the inclination dispersion of the dust grains, which in turn enables constraints to be placed on the system dynamics (e.g., Quillen et al. 2007). At (sub)millimeter wavelengths, dust grains are largely unaffected by radiation pressure and allow direct probing of the mass distribution and level of dynamical excitation of large planetesimals in the system.

With high angular resolution (sub)millimeter observations from the Submillimeter Array (SMA) and the Atacama Large Millimeter/submillimeter Array (ALMA), dozens of debris disks have now been resolved, including some with marginally resolved or partially constrained vertical structures (e.g., Kennedy et al. 2018; Daley et al. 2019; Matrà et al. 2019; Vizgan et al. 2022; Hales et al. 2022; Marshall et al. 2023; Terrill et al. 2023; Matrà et al. 2025). Matrà et al. (2019) resolved the vertical structure in the disk around β Pictoris and found that a two-component vertical structure including both a dynamically cold and hot population was required to match the data – a structure that is reminiscent of the dual populations of planetesimals in the classical Kuiper belt caused by Neptune’s migration early in the history of the Solar System (e.g., Brown 2001; Morbidelli et al. 2008; Petit et al. 2011). Daley et al. (2019) resolved the vertical structure of the edge-on debris disk around AU Mic using 1.3 mm data and used grain-size dependent velocity dispersions in the collisional cascade models of Pan & Schlichting (2012) to infer a total mass of stirring bodies in the debris belt no greater than 1.8 M_⊕. Incorporating additional observations, Vizgan et al. (2022) examined the relationship between grain size and scale height by obtaining scale height values using both 1.3 mm and 450 µm data. Kennedy et al. (2018) modeled the HR 4796A (HD 109573) debris disk and determined that their modeling consistently returned a marginally resolved value for the scale height; however, this value was somewhat inconsistent with the indications from scattered light of a dynamically cold population of planetesimals (e.g., Kenyon et al. 1999). These observations revealed the promise of using vertical structure to interpret the dynamics of debris disk systems, although larger samples are needed to draw more robust conclusions.

While (sub)millimeter observations closely trace the masses and dynamical histories of debris disks, scattered light observations probe smaller grains that are sensitive to both radiation forces and collisions. At these wavelengths, debris disks are expected to have a minimum vertical aspect ratio full width at half maximum (FWHM) of h_min = 0.04 ± 0.02, which is driven by radiation pressure and collisions inflating inclinations and eccentricities of the small grains (Thébault 2009). The vertical structures of debris disks have already been characterized for several disks in scattered light, with most scale height measurements at or above the theoretically predicted minimum (e.g., Milli et al. 2014; Millar-Blanchaer et al. 2015, 2016; Engler et al. 2019; Boccaletti et al. 2019; Olofsson et al. 2020, 2022b).

However, a larger sample of scale height measurements at (sub)millimeter wavelengths is needed to complement this sample of scattered light scale heights and illuminate connections between grain size and system dynamics. For example, small grains may be stirred to higher inclinations by a giant planet in the disk, inflating the observed scale height (Quillen 2006; Pan & Schlichting 2012). Furthermore, comparing the vertical structures from millimeter and scattered light observations may enable constraints on the gas mass of debris disks as well as on the origins of the gas (Olofsson et al. 2022b). However, if grains are fully coupled to the gas or turbulent diffusion is significant, inclination damping could be reduced, preventing micron-sized grains from fully settling into the midplane and resulting in a more vertically extended grain distribution (Marino et al. 2022). More observational measurements of vertical disk structures at multiple wavelengths will reveal how these processes sculpt debris disks.

The ALMA survey to Resolve exoKuiper belt Substructures (ARKS) sample contains ALMA observations of 24 debris disks at a high enough resolution and sensitivity to carry out the detailed characterization of radial and vertical dust structures (Marino et al. 2026). In order to resolve the vertical distribution of dust, observations of highly inclined (i.e., close to edge-on) debris disks at high angular resolution are necessary. The most highly inclined debris disks in the ARKS sample provide the opportunity to begin building a larger sample of vertically well-resolved debris disks at millimeter wavelengths.

In this paper, we analyze these highly inclined ARKS debris disks in order to place new constraints on their vertical dust distributions. The observations and sample of debris disks are described in Sect. 2. Our modeling framework is detailed in Sect. 3, and our analyses are presented in Sect. 4. In Sect. 5 we relate our parametric modeling results to theoretical models of debris disk dynamics to infer the masses of stirring bodies and the dynamical history of the system, compare the parametric and nonparametric modeling results, and investigate the relationship between vertical structure at millimeter and scattered light wavelengths. Our main conclusions are summarized in Sect. 6. Some technical details are provided in the appendices.

2 Observations

ARKS (2022.1.00338.L, PI: S. Marino) aimed to resolve the FWHM of the vertical distribution of the most highly inclined (i > 75^◦) disks in the sample. The new Band 7 observations were configured to resolve the vertical FWHM of the disks with at least two resolution elements, assuming a vertical aspect ratio¹ of h_σ = 0.05. The observing strategy for ARKS targets with i < 75^◦was not optimized for vertical structure recovery, and observations of these targets often lack the resolution or S/N to firmly constrain vertical dust distributions (assuming the same vertical aspect ratio of h_σ = 0.05). However, parametric and/or nonparametric methods can recover vertical structure constraints for some disks at i < 75^◦. In addition to the methods described in this work, for example, azimuthal brightness variations in moderately inclined, optically thin dust rings can be exploited to estimate the disk scale height, provided that the belt is relatively vertically thick and radially narrow (Doi & Kataoka 2021; Villenave et al. 2025). In this paper, we include all disks with vertical structure constraints, though the focus is on the most highly inclined disks, which are amenable to detailed vertical characterization via parametric modeling.

The ARKS sample also includes six archival observations of well-resolved debris disks. The archival observations for HD 39060 (β Pic), HD 107146, HD 197481 (AU Mic), and HD 206893 were taken in Band 6. The vertical aspect ratio of a disk is known to vary with observing wavelength; this variation has been measured in AU Mic between 1.3 mm and 450 μm, with the scale height decreasing by a factor of ∼1.3 with a factor of about three decrease in wavelength (Vizgan et al. 2022). The difference between scale height measurements at Band 7 (0.89 mm) and Band 6 (1.3 mm) is likely more subtle given that the wavelength differs by less than a factor of two, but there are no Band 7 observations of these targets at sufficient resolution to confirm this. The full target selection and observing strategy is justified in Marino et al. (2026), as well as further details about the observations and data calibration. In this study, we use the corrected data produced and presented in Marino et al. (2026), which account for weights and phase center and flux offsets between execution blocks.

3 Modeling

3.1 Definitions of scale height

In general, a vertical scale height characterizes the thickness of the disk by quantifying the size of the region around the disk midplane from which the majority of disk emission originates. This can be defined in a number of ways, and can depend on the exact vertical distribution of disk material. Here, we parameterize the disk scale height, H, as $$H(r)=hr$$(1) such that the vertical aspect ratio, h, is constant at all radii, r.

It is common to express a disk scale height as the standard deviation of a Gaussian distribution, H_σ. In this work, we consider multiple vertical distributions, including Lorentzian and multi-component vertical dust distributions which do not have a clearly defined standard deviation. Therefore, we report scale heights in terms of the half width at half maximum (HWHM) of the vertical profile, H_HWHM. We use the HWHM rather than the full width at half maximum (FWHM) to be more directly comparable with legacy H_σ values. The different scale height formats are related by $${}_{h_{\text{HWHM}}=\frac{12}{}{h}_{\text{FWHM}}=\sqrt{2\ln 2}{h}_{\sigma }.}^{H_{\text{HWHM}}=\frac{12}{}{H}_{\text{FWHM}}=\sqrt{2\ln 2}{H}_{\sigma }}$$(2)

We use the above notation to express which scale height or aspect ratio definition is used. For a more direct comparison across different methodologies and archival studies, we preferentially use H_HWHM and h_HWHM, converting to this format when possible. A complete list of aspect ratio values in their original and converted formats is available in Table C.1.

3.2 Parametric modeling

We use a parametric debris disk modeling code adapted from the publicly available code described in Flaherty et al. (2015)², which was developed from earlier disk modeling works (Dartois et al. 2003; Rosenfeld et al. 2012, 2013). The modeling framework, optimized for modeling the continuum (dust) emission of axisymmetric debris disks³, is described in detail in Fehr (2023) and summarized below.

As debris disks are optically thin, their temperatures and densities are fully degenerate. Since we are modeling observations at a single wavelength, the geometric results are insensitive to the temperature structure and we assume a blackbody equilibrium temperature structure, allowing us to fit for the density structure. We computed the dust temperature by setting the flux received by the disk from its host star at some distance d from the star $d=\sqrt{r^2+z^2}$, where r and z are the radial and vertical distances in cylindrical coordinates, respectively) equal to the flux emitted by the dust grains at that distance. The bolometric flux emitted at the surface of a dust grain is $$F_{\text{dust}}(r)=\frac{L_{\ast }}{16\pi d^2},$$(3) where L_∗ is the bolometric luminosity of the host star. Assuming grains emit and absorb in a manner similar to black bodies, the dust temperature is then $$T_{\text{dust}}(r)={\left(\frac{L_{\ast }}{16\pi \sigma d^2}\right)}^{\frac{\mathrm{1/4}}{}},$$(4) where σ is the Stefan-Boltzmann constant.

The density structure ρ(r, z) of the disk is defined by two components: a radial profile form and a vertical profile form. Careful treatment of the radial structures is critical, as they can be degenerate with the vertical structure of the disk. Thus, the set of radial functional forms applied for the vertical parametric modeling presented here was informed by the ARKS radial structure analysis (Han et al. 2026). Most of our models use either a radial double power law or a radial double Gaussian; the functional forms and associated free parameters are shown in the top section of Table 1. For a more complete examination of the disk radial structures and functional forms, refer to Han et al. (2026).

For the disk vertical structure, we parameterize the disk scale height according to Eq. (1), such that the vertical aspect ratio h is constant throughout the disk in order to minimize additional free parameters. We discuss this assumption and its implications in Sect. 5.6. The vertical functional forms used in this study are shown in the bottom section of Table 1, and are motivated by the array of vertical functional forms used in other debris disk studies (e.g., Golimowski et al. 2006; Lagrange et al. 2012; Apai et al. 2015; Matrà et al. 2019).

3.2.1 Radiative transfer

Once the temperature and density structures of the disk have been defined, we compute the flux throughout the disk in order to obtain a model image of the disk with projected major axis x, projected minor axis y, and distance from disk center along the line of sight s. The intensity of the dust emission is $$I_v(x,y)={\displaystyle {\int }_0^{\infty }{S}_v(x,y,s)\exp \left[-{\tau }_v(x,y,s)\right]K_v(x,y,s)ds},$$(5) where S_ν(x, y, s) is the source function, τ_ν(x, y, s) is the optical depth, and K_ν(x, y, s) is the absorption coefficient. We approximated the source function with the Planck function: $$S_v(x,y,s)=\frac{2h_p{v}^3}{c^2}\frac{1}{\exp \left[\frac{h_{p}v}{k_{B}T(x,y,s)}\right]-1},$$(6) where h_p is Planck’s constant, c is the speed of light, and k_B is the Boltzmann constant.

The optical depth is $${\tau }_v(x,y,s)={\displaystyle {\int }_0^s{K}_v\left(x,y,s^{\prime }\right)d{s}^{\prime }}.$$(7)

We assumed $$K_v(x,y,s)={\kappa }_v\rho (x,y,r),$$(8) where κ_ν is the dust mass opacity. We adopted κ_ν =1.9 cm² g⁻¹ at 0.89 mm and κ_ν =1.3 cm² g⁻¹ at 1.3 mm as in Marino et al. (2026), who assumed a power law grain size distribution with a slope of –3.5, and a composition of astrosilicates (70%), crystalline water ice (15%), and amorphous carbon (15%) by mass (Draine 2003; Li & Greenberg 1998).

After generating the model disk image and multiplying it by a Gaussian approximation of the primary beam⁴, we used

galario

to obtain synthetic visibilities sampled at the same spatial frequencies as the data (Tazzari et al. 2018). The stellar flux is added to the image in the visibility domain at the disk center, which is set by free parameters that enable offsets in right ascension (dRA) and declination (dDec) from the phase center.

3.2.2 Model parameters

We fit the synthetic model visibilities to the observed visibilities using Markov Chain Monte Carlo (MCMC) methods. We used the affine-invariant MCMC implementation in

emcee

to obtain best-fit model parameters and their posterior distributions (Foreman-Mackey et al. 2013). For each disk model, the free parameters included the disk inclination (i), position angle (PA), surface density normalization (Σ_c), stellar flux (F_⋆) at the observing wavelength⁵ (F_∗), and offsets from the phase center (dRA and dDec). The additional free parameters for each model depend on the choice of radial and vertical density forms, presented in Table 1.

For each model, we initialized the MCMC with a 64 walker ensemble (more than double the number of free parameters for each model). To ensure that the parameter space was sufficiently sampled, we required at least 10⁵ samples after convergence, which corresponds to approximately 1500 steps for each of the 64 walkers. Convergence was determined by visually inspecting the evolution of all walkers in both parameter space and log-likelihood space. We identified convergence once the log-likelihood stabilized and the walkers fluctuated around consistent parameter ranges, discarding earlier steps as burn-in. We further discarded walkers with log-likelihood values at least 5σ lower than the maximum log-likelihood value in order to prevent individual unconverged walkers from affecting the results.

The log-likelihood is proportional to −0.5χ², where the model and data are compared using a χ² metric⁶: $${\chi }^2={\displaystyle \sum }_{j=1}^N{w}_j\left[{\left({\text{ReV}}_{{\text{obs}}_j}-{\text{ReV}}_{\mathrm{mod}j}\right)}^2+{\left({\text{ImV}}_{{\text{obs}}_j}-{\text{ImV}}_{\mathrm{mod}j}\right)}^2\right],$$(9) where w_j are the visibility weights and V_{obs j} and V _{mod j} are the observed and model visibilities, respectively.

3.2.3 Model comparison

We tested several model parameterizations for each source, typically beginning with two simple cases: a vertical Gaussian with a radial double power law, and a vertical Gaussian with a radial double Gaussian (Tables 1 and A.1). If radial parametric modeling results favored a different radial density distribution, we included additional models (see Han et al. 2026). After running the initial suite of models with a vertical Gaussian density distribution, we expanded the range of models to include additional vertical complexities.

We compared different model parameterizations using the Akaike and Bayesian Information Criteria (AIC and BIC, respectively), which quantify the trade-off between how well the different models fit the data and how many free parameters are needed (e.g., Akaike 1974; Schwarz 1978; Burnham et al. 2002; Liddle 2007). These criteria are defined as $$\begin{array}{c}\text{AIC}\equiv -2\ln {ℒ}_{\max }+2k,\\ \text{BIC}\equiv -2\ln {ℒ}_{\max }+k\ln N,\end{array}$$(10) where ℒ_max is the maximum likelihood achievable by the model, k is the number of parameters of the model, and N is the number of data points (in this application, the number of visibilities).

The set of model parameterizations tested for each source are shown in Table A.1, including the AIC and BIC comparisons. These criteria were computed such that $$\begin{array}{c}{\text{ΔAIC}}_{\text{i}}={\text{AIC}}_{\text{i}}-{\text{AIC}}_{\min },\\ {\text{ΔBIC}}_{\text{i}}={\text{BIC}}_{\text{i}}-{\text{BIC}}_{\min },\end{array}$$(11) where AIC_min and BIC_min are the minimum values of the criteria across the full suite of models, indicating the best fit to the data. Thus, for both ∆AIC and ∆BIC, a value of 0 indicates a best-fit model. AIC comparisons can be converted to a relative likelihood, which is defined as p_i = exp [−∆AIC_i/2]. Assuming a normal distribution with a mean of 0, we compute a symmetric interval which contains probability 1 − p_i. We present AIC comparisons using the upper bound of this interval, which gives the relative likelihood in terms of a confidence level σ. BIC results are directly presented as ∆BIC. Generally, ∆BIC > 10 is regarded as strong evidence against a particular model (>99% confidence, see Kass & Raftery 1995; Raftery 1995).

We use the AIC and BIC to identify a fiducial model for each source. When the ∆AIC and ∆BIC are minimized for the same model, this model is selected as the fiducial model. When the ∆AIC and ∆BIC are minimized for different models, the model which is least strongly ruled out by the AIC and BIC (AIC confidence < 3σ and ∆BIC < 10) is selected as the fiducial model.

3.3 Nonparametric modeling

Vertical information can also be extracted using nonparametric methods. We used both frank, which operates in the visibility domain, and rave, which operates in the image domain. Each is briefly summarized below. (For more detailed descriptions in the context of radial profile recovery, see Han et al. 2026.)

3.3.1 Frank

frank is an open-source code which reconstructs the 1D (azimuthally averaged) radial intensity profile of a disk by fitting the visibilities with a Fourier-Bessel series (Jennings et al. 2020). While frank is primarily a tool for extracting radial profiles, we use an extension which enables constraints to be placed on the vertical aspect ratio of the disk (originally implemented in Terrill et al. 2023). The specific implementation of frank used for analyzing the ARKS sample (arksia, Jennings et al. 2024) is publicly available⁷.

For all fits, we adopt the position angle, inclination, and RA and declination offsets presented in Marino et al. (2026). We assumed the disk to be vertically Gaussian, with the scale height H_σ proportional to the radius, thereby defining an aspect ratio of h_σ = H_σ(r)/r. Following the approach in Terrill et al. (2023), we repeated the radial profile fit over a range of scale height assumptions, computing the χ² value for each height and estimating its likelihood as exp(−χ²/2). The radial profile fits are presented in Han et al. (2026); here, we focus only on the vertical constraints obtained with frank.

Starting from a wide range in h_σ, we narrowed the height grid to more densely sample the region of h_σ in which the bulk of the probability mass is centered. The sample points were then interpolated quadratically to obtain a finer probability density distribution (PDF). For approximately half of the disks, the PDF exhibits a clear peak and approaches 0 or a value close to 0 on either side of the peak, and the best-fit height was estimated with the median and the uncertainties with the 16th and 84th percentiles. For the rest of the disks, the PDF does not approach 0 as h_σ approaches 0. In such cases, an upper limit was estimated with the 99.7th percentile if the PDF approaches 0 toward large values of h_σ. The fitted scale heights are listed in Table 3, and the PDFs are shown in Fig. B.2.

3.3.2 rave

rave is an open-source package for recovering the radial brightness profiles of edge-on disks by fitting models of concentric annuli in the image domain (Han et al. 2022). Like frank, rave also enables vertical height aspect ratio fitting. Assuming the disk to be vertically Gaussian, the radial profile was fit at a range of h_σ assumptions including 0, and nine additional points spaced uniformly in logarithmic space between 0.01 and 0.5. The χ² value for the model fit at each height assumption was calculated using the root mean square (rms) noise of the CLEAN image and assuming the number of independent elements to be equal to the number of beams in the region of the image containing disk flux (defined as a rectangle within r_max from the star along the major axis and y_max along the minor axis, as listed in the appendix of Han et al. 2026). The likelihood of each height was estimated to be proportional to exp(−χ²/2).

As with frank, the sample points were interpolated to obtain a finer PDF. For disks with PDFs that show a clear peak, the best-fit height was estimated with the median and the lower and upper uncertainties were estimated with the 16th and 84th percentiles. For the remainder of disks, the PDF does not approach 0 even when h_σ is 0. In such cases, if the PDF only drops to 0 at h_σ ⪆ 0.3, we take the 84th percentile to be an upper limit on h_σ; otherwise, the height was unconstrained. The scale heights fitted with rave are listed in Table 3, and the PDFs are shown in Fig. B.3.

4 Analysis

The fiducial parametric models are summarized briefly in Table 2. The fiducial images and residuals for each disk are shown in Figure 1 using the same imaging parameters as in Marino et al. (2026). For most sources the fiducial models closely match the observations, with no positive or negative features in the residuals at the 5σ level. Figure 2 shows the best-fit vertical profile for each disk at its reference radius (see Table 4), defined as the radial location of peak intensity in our models, which are axisymmetric. We find aspect ratios in the range between roughly 0.003 ≤ h_HWHM ≤ 0.19. We note that while we present h in terms of the HWHM of the vertical distribution, many previous studies (particularly those which assume a Gaussian vertical distribution) present h_σ, the standard deviation of the vertical profile. For comparison with other studies, we have converted the published h_σ values mentioned in this work to h_HWHM values.

For each source, the aspect ratio is meaningfully constrained, i.e., the posterior probability distribution of h_HWHM values has a nonzero peak (Figure 3). To determine whether the scale heights are fully or marginally resolved, we take the difference between the median and 16th percentile h values in order to approximate a 1σ confidence interval (this difference is equivalent to a 1σ confidence interval for a perfectly Gaussian posterior probability distribution, though not all of the posteriors are Gaussian). We use this interval to extrapolate to a 3σ confidence interval; for sources with a h_HWHM ≤ 0 at the 3σ lower bound, we consider h_HWHM to be marginally resolved, placing only an upper limit on h_HWHM. We apply this extrapolation because h_HWHM ≤ 0 is not allowed by our model. Therefore, the 0.15th percentile (corresponding to ∼3σ for a Gaussian distribution) will always be positive and non-zero.

We additionally examine how well-resolved h_HWHM is with respect to the resolution of the observations. Importantly, the beams shown in Figure 1 do not reflect the smallest scales resolvable in our analyses. Because our parametric modeling occurs in the visibility domain, we retain spatial information up to the nominal resolution of the interferometer, θ_min = λ_obs/b_max, where λ_obs is the observing wavelength and b_max is the longest baseline in the array during the observations. As a result, h_HWHM is more easily resolved in Fourier space than in image space, and some values of h_HWHM may be resolvable even if the vertical extent of the disk is on roughly the same scale as the CLEAN beam (Figure 2). Furthermore, it is not necessary that θ_min < h_HWHM in order to resolve the aspect ratio, as h_HWHM does not represent the full vertical extent of the disk. Rather, h_HWHM only describes the distance from the disk mid-plane to the height at which the disk has half of its peak density. Depending on the vertical parametrization, the disk could thus host a substantial amount of material at several times h_HWHM.

For example, a vertical Gaussian distribution contains about 95% of the material within ±2 h_σ from the midplane, equivalent to ∼3.4 h_HWHM. In this case, a scale height may not be resolved by the observations if θ_min > 3.4 h_HWHM. A vertical Lorentzian distribution has thicker tails, with more disk material further from the midplane: ∼95% within ∼25 h_HWHM centered on the mid-plane. In this case, a scale height may not be resolved by the observations if θ_min > 25 h_HWHM. Measurements of h_HWHM that are limited by the data resolution are diskussed in the source-specific sections, and values of θ_min are provided in Table 4 for reference.

4.1 Marginally resolved scale heights

4.1.1 HD 109573 (HR 4796A)

Parametric models in Han et al. (2026) suggest that the radial density distribution of this source is well-fit by a double Gaussian, assuming vertical Gaussian profile with a fixed vertical aspect ratio of h_HWHM = 0.015. Incorporating scale height as a free parameter, we examined double power law, double Gaussian, Gaussian, and asymmetric Gaussian radial profiles. We find that the double Gaussian continues to outperform the other radial density distributions even when incorporating vertical flexibility into the model. We then modeled HD 109573 with a radial double Gaussian density distribution and the four different vertical parameterizations presented in Table 1.

The AIC and BIC strongly favor the vertical Gaussian and vertical Lorentzian models over the other vertical forms, with the vertical Gaussian model preferred over the vertical Lorentzian only at low confidence (AIC σ = 0.09). However, the posterior probability distributions are broad, with a vertical aspect ratio of $h_{\text{HWHM}}={0.0082}_{-0.0038}^{+0.0033}$ (vertical Gaussian) or $h_{\text{HWHM}}={0.0049}_{-0.0014}^{+0.0016}$ (vertical Lorentzian). Both of these models exhibit a non-Gaussian posterior probability distribution for h_HWHM, truncating near zero. As a result, we only marginally resolve the scale height for HD 109573. We take the upper limit to be the 84th percentile of the posteriors, h_HWHM < 0.0116, using the fiducial vertical Gaussian model. Nonparametric fitting also constrains only upper limits, with frank finding h_HWHM < 0.012 and rave finding h_HWHM < 0.022.

The vertical structure of HD 109573 has previously been studied using Band 7 ALMA observations (∼0.″17 resolution, Kennedy et al. 2018). Kennedy et al. (2018) modeled the disk with a radial Gaussian and a vertical Gaussian, finding H_FWHM = 7 ± 1 au at the radial location of the ring (78.6 au), corresponding to an aspect ratio of $h_{\text{HWHM}}={0.0045}_{-0.006}^{+0.006}$. Matrà et al. (2025) reanalyzed these observations, also modeling HD 109573 with radial and vertical Gaussians. At the radial ring location of 77.8 au, they obtain H_HWHM = 5.3 au, or $h_{\text{HWHM}}={0.0068}_{-0.002}^{+0.002}$. Han et al. (2025) used fave, an extension of rave optimized for less inclined disks, to fit the aspect ratio of HD 109573, finding $h_{\text{HWHM}}={0.0057}_{-0.022}^{+0.016}$. With the new Band 7 ARKS data (0.″08), we find a similar ring location (77.6 au), but a much smaller vertical extent (H_HWHM < 0.9 au). This discrepancy is likely due to the lower resolution of previous observations combined with the different best-fit functional form found here (double rather than single Gaussian); the degeneracy between radial and vertical structure makes it difficult to place strong constraints on the scale height, and this issue is more pronounced with lower-resolution data. Despite the improved resolution of the ARKS observations, we may still be limited by the data resolution of the data (θ_min = 3.63 au). At a reference radius of 77.6 au (the location of peak brightness in the fiducial model), we obtain H_HWHM = 0.636 au. For a vertical Gaussian profile, 95% of the disk material would then lie within 2.16 au centered on the disk midplane, about 60% of the distance that could be resolved by the observations.

4.1.2 HD 131488

Parametric models in Han et al. (2026) suggest that the radial density distribution of this source is best fit by a double Gaussian, assuming a vertical Gaussian profile with a fixed vertical aspect ratio of h_HWHM = 0.005. We find that a double Gaussian radial density distribution is strongly preferred over a double power law or single Gaussian, assuming a vertical Gaussian profile with the vertical aspect ratio as a free parameter. We additionally modeled HD 131488 using a radial double Gaussian density distribution and both a vertical exponential and a vertical Lorentzian.

We find that the vertical Gaussian model performs best when considering both the AIC and BIC (AIC σ = 0.90 over the exponential, ∆BIC = 0), so we select this as the fiducial model for HD 131488. The radial double Gaussian with a vertical exponential is most strongly favored by the AIC, but disfavored by the BIC (∆BIC = 11.63). The radial double Gaussian with a vertical Lorentzian form is not strongly ruled out by either the AIC or BIC tests (AIC σ = 1.26, ∆BIC = 1.14), but because the fiducial Gaussian model provides a better fit with the same number of parameters, we do not increase model complexity any further (i.e., we do not run a vertical double Gaussian model, which could have a similar vertical profile to the Lorentzian model at the cost of two additional free parameters). The aspect ratio is well-constrained, with $h_{\text{HWHM}}={0.0059}_{-0.0010}^{+0.0011}$. The aspect ratio constraints from frank and rave are in agreement, finding h_HWHM < 0.008 and $h_{\text{HWHM}}={0.0016}_{-0.006}^{+0.005}$, respectively. These values are consistent with the uncertain constraint reported in Worthen et al. (2024), $h_{\text{HWHM}}={0.09}_{-0.05}^{+0.05}$, which was derived from Band 6 ALMA observations of HD 131488 (∼0.″5 resolution).

We note that while our measurements of h_HWHM are constrained according to the posterior probability distributions, the interpretation of our measured values of h_HWHM may be limited by the resolution of the data. The observations of HD 131488 have a nominal resolution of θ_min = 3.32 au. At a reference radius of 90.8 au, we obtain H_HWHM = 0.536 au. For a vertical Gaussian profile, 95% of the disk material would then lie within 1.82 au centered on the disk midplane, nearly half the distance that could be resolved by the observations. We take the upper limit of h_HWHM to be the 84th percentile of the posteriors for the fiducial model, h_HWHM < 0.0070, though follow-up observations at higher resolution are necessary to confirm this measurement.

4.2 Resolved Gaussian scale heights: HD 14055

The only source in our sample which is best fit with a vertical Gaussian distribution and has a fully resolved scale height is HD 14055. Parametric models in Han et al. (2026) suggest that the radial density distribution of this source is well-fit by a double power law, though a double Gaussian is not ruled out at high confidence. In the radial models, a vertical Gaussian profile was used with a fixed vertical aspect ratio of h_HWHM = 0.0475. Incorporating scale height as a free parameter, we examined double power law and double Gaussian radial density distributions with a vertical Gaussian. The AIC favored the double Gaussian model, while the BIC favored the double power law model, so we modeled both radial double power law and radial double Gaussian density distributions with the four different vertical parameterizations presented in Table 1. Of these eight models, we select the radial double power law with a vertical Gaussian as the fiducial model, as it is favored at high confidence with the BIC and not strongly ruled out with the AIC (σ = 0.56 over a radial double Gaussian with a vertical Gaussian).

We measure $h_{\text{HWHM}}={0.053}_{-0.008}^{+0.007}$ with the fiducial parametric model, using the 0.″23 ARKS observations in Band 7. The vertical aspect ratio is also resolved with frank and rave: frank finds $h_{\text{HWHM}}={0.037}_{-0.016}^{+0.010}$, while rave finds $h_{\text{HWHM}}={0.047}_{-0.031}^{+0.038}$. These values are in agreement with Terrill et al. (2023), Matrà et al. (2025), and Han et al. (2025), who used Band 6 observations (∼2″ resolution) of HD 14055 to find limits of h_HWHM < 0.12, $h_{\text{HWHM}}={0.04}_{-0.02}^{+0.02}$, and $h_{\text{HWHM}}={0.131}_{-0.075}^{+0.073}$, respectively.

4.3 Resolved complex structures

In this section we present sources that are consistent with having an extended or multi-component millimeter vertical emission profile (i.e., statistically preferred models use either a Lorentzian or a double Gaussian). We note that while we report a single h_HWHM for each vertical Lorentzian model, the broad tails of a Lorentzian can approximate a combination of broad and narrow components with different widths (see Sect. 5.1, where we also discuss the possible mechanisms that could generate these distributions).

4.3.1 HD 9672 (49 Ceti)

Parametric models in Han et al. (2026) suggest that the radial density distribution of this source is best fit by a double power law, assuming a vertical Gaussian profile with a fixed vertical aspect ratio of h_HWHM = 0.07. We modeled HD 9672 using a radial double power law density distribution and the four different vertical parametrizations presented in Table 1, as well as with a radial double Gaussian and vertical Gaussian. We find that the vertical Lorentzian model is most favored by the BIC, and not strongly ruled out by the AIC (σ = 0.27). The vertical double Gaussian model is most strongly favored by the AIC, and disfavored by the BIC (∆BIC = 25.81). We identify the vertical Lorentzian model as the fiducial model for HD 9672.

We obtained a vertical aspect ratio of $h_{\text{HWHM}}={0.0046}_{-0.008}^{+0.008}$ for the fiducial model. As the double Gaussian model is only disfavored by the BIC, which penalizes additional parameters more heavily than the AIC, we may tentatively resolve two distinct dynamical populations within the disk. The first has a measured aspect ratio of $h_{\text{HWHM}}={0.031}_{-0.013}^{+0.022}$ containing ${68}_{-15}^{+11}\%$ of the dust mass, and the other has an aspect ratio of $h_{\text{HWHM}}={0.105}_{-0.040}^{+0.021}$ containing the rest of the mass. The nonparametric aspect ratios fall in between, with $h_{\text{HWHM}}={0.061}_{-0.004}^{+0.004}$ from the frank fit and $h_{\text{HWHM}}={0.058}_{-0.024}^{+0.017}$ from the rave fit.

The aspect ratio measured for the dynamically cooler population of the parametric vertical double Gaussian model is in agreement with the values reported in Delamer (2021) from an analysis of Band 8 (614 µm, ∼0.″ 25 resolution) and Band 7 (868 µm, ∼0.″13 resolution) continuum emission of HD 9672 ($h_{\text{HWHM}}={0.088}_{-0.005}^{+0.005}$ and $h_{\text{HWHM}}={0.0081}_{-0.008}^{+0.006}$, respectively), assuming a Gaussian vertical profile.

With the same Band 8 observations, Terrill et al. (2023) found $h_{\text{HWHM}}={0.0059}_{-0.008}^{+0.008}$ and Han et al. (2025) found $h_{\text{HWHM}}={0.039}_{-0.021}^{+0.016}$, with both studies using a nonparametric model with an assumed Gaussian vertical density distribution. Matrà et al. (2025) reanalyzed the Band 7 observations with radial and vertical Gaussians, finding an aspect ratio of $h_{\text{HWHM}}={0.089}_{-0.011}^{+0.011}$. Finally, Delamer (2021) also analyzed Band 6 (1.33 mm, ∼0.″85 resolution) observations, finding a slightly larger aspect ratio of $h_{\text{HWHM}}={0.018}_{-0.02}^{+0.02}$.

These differences could be partially driven by the range of observations represented, as the dust scale height may vary with observing wavelength. However, the variation among aspect ratios derived from Band 7 observations is likely driven by the fact that a Lorentzian or two Gaussians have broader tails, while a single Gaussian drops off more quickly away from the midplane. When fitting the disk with only one vertical Gaussian, the model may settle on a larger value of h_HWHM as it attempts to capture the more dynamically active population in the disk. HD 9672 also hosts a likely warp asymmetry, discussed in detail in Lovell et al. (2026). While we do not see significant features in the residual map (Fig. 1), we note that a subtle warp feature could bias our axisymmetric models to larger scale heights.

4.3.2 HD 10647

Parametric models in Han et al. (2026) suggest that the radial density distribution of HD 10647 is well-fit by a double Gaussian, assuming a vertical Gaussian profile with a fixed vertical aspect ratio of h_HWHM = 0.055. We modeled the disk using radial double power law and double Gaussian density distributions with the four different vertical parametrizations presented in Table 1. The AIC strongly favors the model with both radial and vertical double Gaussian density distributions, while ∆BIC = 9.042. The BIC favors the radial double power law with a vertical Lorentzian, though the AIC strongly disfavors this model (AIC = 6.151σ). We thus identify the radial and vertical double Gaussian model as the fiducial model.

We resolve two distinct dynamical populations with the fiducial model of HD 10647, one with a measured scale height of $h_{\text{HWHM}}={0.053}_{-0.007}^{+0.006}$ containing ${88}_{-6}^{+6}\%$ of the dust mass, and the other with a scale height of $h_{\text{HWHM}}={0.225}_{-0.002}^{+0.002}$ containing the rest of the mass. This is in line with the constraint from rave, $h_{\text{HWHM}}={0.040}_{-0.021}^{+0.016}$. However, the upper limit derived from the frank fit is much smaller, h_HWHM < 0.008. While we do not see significant features in the residual map of HD 10647 (Fig. 1), this source has a major axis asymmetry that could bias the results of our axisymmetric modeling (Lovell et al. 2026).

Several studies have already made measurements which are in agreement with our measurement of the dynamically cooler population of the fiducial model (which contains most of the disk mass): Lovell et al. (2021) found the scale height to be consistent with $h_{\text{HWHM}}={0.05}_{-0.001}^{+0.001}$, while Han et al. (2025) found $h_{\text{HWHM}}={0.054}_{-0.012}^{+0.011}$. Terrill et al. (2023) and Matrà et al. (2025) found marginal measurements of $h_{\text{HWHM}}={0.044}_{-0.008}^{+0.009}$ and $h_{\text{HWHM}}={0.06}_{-0.001}^{+0.001}$, respectively.

4.3.3 HD 15115

The parametric modeling in Han et al. (2026) shows that HD 15115 is well-represented by a double Gaussian radial density distribution, assuming a vertical Gaussian profile with a fixed vertical aspect ratio of h_HWHM = 0.02. We modeled HD 15115 using a radial double Gaussian density distribution and the four different vertical parametrizations presented in Table 1. We identify the vertical Lorentzian model as the fiducial model for HD 15115, as it is strongly favored by both the AIC and BIC. The scale height is well-constrained, with a vertical aspect ratio of $h_{\text{HWHM}}={0.222}_{-0.0015}^{+0.0013}$ for the vertical Lorentzian model. The frank and rave values are comparable, $h_{\text{HWHM}}={0.019}_{-0.001}^{+0.001}$ and $h_{\text{HWHM}}={0.022}_{-0.002}^{+0.002}$, respectively.

MacGregor et al. (2019) found that HD 15115 is not vertically resolved in 0.″6 (29 au) Band 6 observations. Using the same Band 6 observations, Matrà et al. (2025) marginally detected the scale height ($h_{\text{HWHM}}={0.06}_{-0.02}^{+0.01}$) by modeling the disk with a radial and vertical Gaussian. This is in good agreement with results from Terrill et al. (2023) and Han et al. (2025), who used nonparametric methods with a vertical Gaussian assumption to measure $h_{\text{HWHM}}={0.56}_{-0.008}^{+0.008}$ and $h_{\text{HWHM}}={0.79}_{-0.040}^{+0.035}$, respectively. The differences between these values from previous observations and our measurement are likely due to the improved resolution of the ARKS data (0.″12), which is 5 times higher than the previous observations of HD 15115, combined with the different assumed radial and vertical distributions. In addition, there is evidence that HD 15115 is eccentric (Lovell et al. 2026), though this may be compensated for by the offset parameters dRA and dDec in our modeling (see Table A.1).

4.3.4 HD 32297

Parametric models in Han et al. (2026) suggest that the radial density distribution of HD 32297 is well-fit by a double power law, assuming a vertical Gaussian profile with a fixed vertical aspect ratio of h_hwhM = 0.012. We modeled the disk using a radial double power law density distribution and the four different vertical parametrizations presented in Table 1. The vertical Lorentzian model is favored by both the AIC and BIC, so we select it as the fiducial model. The aspect ratio of the fiducial model is well-constrained ($h_{\text{HWHM}}={0.092}_{-0.0012}^{+0.0012}$) and in line with fits from frank ($h_{\text{HWHM}}={0.010}_{-0.001}^{+0.001}$) and rave ($h_{\text{HWHM}}={0.015}_{-0.004}^{+0.003}$).

The AIC score of the vertical double Gaussian model indicates that it also reproduces the data well (σ = 0.97); however, due to the additional parameters in the vertical double Gaussian model, the BIC score is significantly worse compared to the Lorentzian model (∆ΒΙC = 27.95). Nevertheless, we marginally resolve two independent dynamical populations within HD 32297: one with a scale height of $h_{\text{HWHM}}={0.83}_{-0.0034}^{+0.0018}$ containing ${92}_{-3}^{+2}\%$ of the dust mass, and another with a scale height of $h_{\text{HWHM}}={0.146}_{-0.047}^{+0.027}$ containing the remainder.

There exist several previous scale height measurements for HD 32297 at ~millimeter wavelengths. Terrill et al. (2023) and Han et al. (2025) use nonparametric methods (assuming a Gaussian vertical distribution) to find $h_{\text{HWHM}}={0.09}_{-0.01}^{+0.01}$ and $h_{\text{HWHM}}={0.139}_{-0.052}^{+0.041}$, respectively. Worthen et al. (2024) find $h_{\text{HWHM}}={0.018}_{-0.002}^{+0.002}$ and Matrà et al. (2025) find $h_{\text{HWHM}}={0.09}_{-0.02}^{+0.01}$, both modeling the disk with Gaussian radial and vertical density distributions. Our measured h_hwhM is smaller than all of the above, likely due to several factors. First, the ARKS data have a higher angular resolution (0.″06). Second, the vertical and radial distributions of disk material can be degenerate with disk inclination; uncertainties in the disk inclination make it more difficult to constrain the scale height and belt width, and all of the above studies find/use smaller inclinations than the value we obtain from our fiducial model ($i=88.48{°}_{-0.03}^{+0.04}$, which does not appear to be affected by this degeneracy). We find that HD 32297 is best modeled with a vertical Lorentzian profile and a radial double power law profile; the common assumption of vertical and/or radial Gaussian profiles used for previous scale height measurements may be biasing the measurement of h_hwhM, resulting in artificially high values. Lastly, there is evidence that HD 32297 is eccentric (Lovell et al. 2026), though this may be compensated for by the offset parameters dRA and dDec in our modeling (see Table A.1).

4.3.5 HD 39060 (Beta Pic)

Parametric models in Han et al. (2026) suggest that the radial density distribution of this source is well-represented by a double power law, assuming a vertical Gaussian profile with a fixed vertical aspect ratio of h_hwhM = 0.0445. We modeled HD 39060 using a radial double power law density distribution and the four different vertical parametrizations presented in Table 1. The vertical Lorentzian (AIC σ = 1.55, ∆ΒΙC = 0) and vertical double Gaussian (AIC σ = 0, ∆ΒΙC = 22.85) models are both good fits to the data. For the fiducial vertical Lorentzian model, the aspect ratio is well-constrained with, $h_{\text{HWHM}}={0.1122}_{-0.0013}^{+0.0013}$. The frank fit yields $h_{\text{HWHM}}={0.073}_{-0.002}^{+0.002}$, while rave finds $h_{\text{HWHM}}={0.086}_{-0.005}^{+0.005}$.

We identified two distinct dynamical populations within the disk with the vertical double Gaussian model; the first has a measured aspect ratio of $h_{\text{HWHM}}={0.046}_{-0.003}^{+0.005}$ containing ${58}_{-3}^{+4}\%$ of the dust mass, and the other has a scale height of $h_{\text{HWHM}}={0.212}_{-0.003}^{+0.002}$ containing the rest of the mass. Matrà et al. (2019) also measured two distinct aspect ratios for this source, finding $h_{\text{HWHM}}={0.130}_{-0.007}^{+0.009}$ containing ${80}_{-5}^{+4}\%$ of the dust mass, and $h_{\text{HWHM}}={0.016}_{-0.007}^{+0.007}$ containing the rest of the mass.

The differences are likely partially driven by the choice in radial density distribution. Here we use a double power law, finding a relatively shallow slope in the inner disk ($\left({\alpha }_{\text{in}}={1.78}_{-0.10}^{+0.14}\right)$) with a sharper decrease in density in the outer disk ($\left({\alpha }_{\text{out}}={5.2}_{-0.5}^{+0.6}\right)$). Matrà et al. (2019) use a radially Gaussian density distribution with $\sigma ={36.7}_{-1.3}^{+1.0}$ au. Interestingly, Matrà et al. (2019) find that the broad component contains most of the disk mass, while we find that the narrow component contains slightly more mass. Matrà et al. (2019) also find marginal evidence of a radially varying aspect ratio in HD 39060, with h(r) = h₀(r/r₀)β⁻¹ and $\beta ={0.7}_{-0.2}^{+0.2}$, which our model specification does not account for.

Asymmetries in the disk could also contribute to differences in the measured scale heights. Though the models presented both here and in Matrà et al. (2019) are axisymmetric, there is evidence of asymmetric structures in the millimeter emission of HD 39060 (Lovell et al. 2026). An asymmetric feature at the ~3σ level can be seen in the northeast of the residual map in Figure 1.

4.3.6 HD 61005

The parametric modeling in Han et al. (2026) shows that HD 61005 is well-represented by a double Gaussian radial density distribution, assuming a vertical Gaussian profile with a fixed vertical aspect ratio of h_hwhM = 0.02. We modeled HD 61005 using a radial double Gaussian and the four different vertical parametrizations presented in Table 1. We find that the vertical Lorentzian model is favored by both the AIC and the BIC at high significance, corresponding to an aspect ratio of $h_{\text{HWHM}}={0.0143}_{-0.0021}^{+0.0020}$. Compared to the vertical double Gaussian model, the AIC only favors the Lorentzian model at 1.31σ. We thus find marginal evidence of two resolved populations in the disk; the first has a measured scale height of $h_{\text{HWHM}}={0.014}_{-0.006}^{+0.005}$ containing ${90}_{-12}^{+6}\%$ of the dust mass, and the other has a scale height of $h={0.117}_{-0.087}^{+0.043}$ containing the rest of the mass. Aspect ratio fits from frank and rave find $h_{\text{HWHM}}={0.019}_{-0.001}^{+0.001}$ and $h_{\text{HWHM}}={0.020}_{-0.007}^{+0.005}$, respectively.

Assuming Gaussian radial and vertical density distributions, Matrà et al. (2025) find $h_{\text{HWHM}}={0.052}_{-0.005}^{+0.004}$. Similarly, Terrill et al. (2023) and Han et al. (2025) find $h_{\text{HWHM}}={0.046}_{-0.005}^{+0.004}$ and $h_{\text{HWHM}}={0.034}_{-0.011}^{+0.008}$, respectively, using nonparametric methods and assuming a Gaussian vertical distribution. These values may be inflated due to the models fitting a more dynamically active disk to a simple vertical Gaussian, driving up the scale height measurement to account for the material further from the disk midplane.

It is also possible that asymmetries in the disk could bias our axisymmetric model toward a vertical structure indicative of multiple populations. While we do not identify structure in our model residuals (see Figure 1), Lovell et al. (2026) find evidence of asymmetry in HD 61005, which could affect the vertical structure measurements presented here.

4.3.7 HD 76582

Parametric models in Han et al. (2026) suggest that the radial density distribution of HD 76582 is well-fit by an asymmetric Gaussian (which is defined by some central radius and two different σ values inside and outside of the central radius), assuming a vertical Gaussian profile with a fixed vertical aspect ratio of h_HWHM = 0.015. We examined the vertical structure of HD 76582 by first modeling the vertical structure of the disk with a Gaussian, and the radial structure with a variety of different radial forms. The best-fitting radial forms were a double power law, asymmetric Gaussian, and error function. We modeled HD 76582 using each of these radial density structures and the four different vertical parametrizations presented in Table 1. The model with a radial asymmetric Gaussian and vertical Lorentzian is favored by both the AIC and BIC, yielding $h_{\text{HWHM}}={0.192}_{-0.009}^{+0.009}$; we thus selected this as the fiducial model for HD 76582. The frank aspect ratio fit ($h_{\text{HWHM}}={0.065}_{-0.012}^{+0.012}$) and the rave upper limit (h_HWHM ≤ 0.097) are smaller.

The rave and frank results are in agreement with the upper limits identified by Matrà et al. (2025) (h_HWHM < 0.1) and Han et al. (2025) (h_HWHM < 0.180). All of these values were obtained by assuming Gaussian vertical density distributions. Furthermore, it is likely that the radial Gaussian density distribution assumed by Matrà et al. (2025) resulted in a biased scale height value. We found that an asymmetric Gaussian yields a good fit to the data, with substantially different widths on each side of the Gaussian peak (σ_in = 45 ± 8 au and ${\alpha }_{\text{out}}={114}_{-15}^{+23}$ au). This corresponds to a FWHM of $\Delta R={187}_{-27}^{+37}$ au, while the Gaussian model of Matrà et al. (2025) has ∆R = 210 ± 20 au. Our asymmetric Gaussian model has its peak at R = 181 ± 11 au, in contrast to $R={219}_{-8}^{+9}$ au in Matrà et al. (2025). As the Gaussian model of Matrà et al. (2025) is both marginally wider and significantly farther out compared to the best-fit asymmetric Gaussian model presented here, it is likely that the scale height measurement was underestimated to compensate for the extra radial width in the model.

4.3.8 HD 131835

Parametric models in Han et al. (2026) suggest that the radial density distribution of this source is well-fit by a double Gaussian, assuming vertical Gaussian profile with a fixed vertical aspect ratio of h_HWHM = 0.015. We find that a radial double Gaussian continues to perform significantly better than a radial double power law or a radial triple Gaussian when the vertical aspect ratio is included as a free model parameter. We then modeled HD 131835 using a radial double Gaussian density distribution and the four different vertical parameterizations presented in Table 1.

The AIC and BIC strongly favor the vertical Gaussian and vertical Lorentzian models, with the Lorentzian model favored at low confidence over the Gaussian model (AIC σ = 0.548, ∆BIC = 1.077). Therefore, we select the Lorentzian model as the fiducial model for HD 131835. The posterior probability distributions are broad for both the fiducial and vertical Gaussian model, with an aspect ratio of $h_{\text{HWHM}}={0.022}_{-0.006}^{+0.007}$ for the fiducial Lorentzian model and $h_{\text{HWHM}}={0.030}_{-0.010}^{+0.010}$ for the vertical Gaussian model. frank and rave find $h_{\text{HWHM}}={0.039}_{-0.006}^{+0.006}$ and $h_{\text{HWHM}}={0.057}_{-0.028}^{+0.020}$, respectively.

4.3.9 HD 161868

Parametric models in Han et al. (2026) suggest that the radial density distribution of this source is well-fit by a double power law, assuming a vertical Gaussian profile with a fixed vertical aspect ratio of h_HWHM = 0.015. We modeled HD 161868 with a radial double power law and the four different vertical parametrizations presented in Table 1. We find that the vertical Lorentzian model is favored by the BIC (AIC σ = 1.255), while the double Gaussian model is favored by the AIC (∆BIC = 20.786). We select the vertical Lorentzian model as the fiducial model for HD 161868.

For the fiducial model, we measure a scale height of $h_{\text{HWHM}}={0.192}_{-0.016}^{+0.012}$. The double Gaussian model is likely dis-favored by the BIC more because of the additional parameters in the model rather than due to a poor fit to the data, as evidenced by the AIC. We thus find evidence that we resolve two distinct dynamical populations within the disk. The first has a measured aspect ratio of $h_{\text{HWHM}}={0.099}_{-0.045}^{+0.034}$ containing ${89}_{-2}^{+3}\%$ of the dust mass, and the other has an aspect ratio of $h_{\text{HWHM}}={0.380}_{-0.023}^{+0.037}$ containing the rest of the mass. Fits with frank and rave yield only upper limits, finding h_HWHM < 0.097 and h_HWHM < 0.135, respectively.

Matrà et al. (2025) marginally resolve the scale height of HD 161868, finding $h_{\text{HWHM}}={0.15}_{-0.06}^{+0.05}$ with a parametric fit assuming Gaussian radial and vertical density distributions. By nonparametrically fitting the visibilities and assuming a Gaussian vertical density distribution, Terrill et al. (2023) find a marginal measurement of $h_{\text{HWHM}}={0.18}_{-0.05}^{+0.04}$ and Han et al. (2025) find $h_{\text{HWHM}}={0.16}_{-0.05}^{+0.04}$. All three of these measurements are consistent with the aspect ratio of our fiducial model.

4.3.10 HD 197481 (AU Mic)

The parametric modeling in Han et al. (2026) shows that models of HD 197481 perform comparably well with a variety of radial density structures, assuming a vertical Gaussian profile with a fixed vertical aspect ratio of h_HWHM = 0.015. We modeled HD 197481 using both radial double power law and double Gaussian density distributions with the four different vertical parametrizations presented in Table 1. There is no single model which stands out as a true best fit, but we identify the radial double power law and vertical Lorentzian model as the fiducial model due to its moderate AIC confidence (σ = 2.491) and strong BIC (∆BIC = 0). For this model, we resolve an aspect ratio of $h_{\text{HWHM}}={0.028}_{-0.0002}^{+0.0009}$, the smallest in our sample.

The radial double Gaussian and vertical exponential model was favored by the AIC (∆BIC = 23.328), yielding a scale height measurement of $h_{\text{HWHM}}={0.0338}_{-0.0015}^{+0.0014}$ with a profile that drops off sharply ($\left(\gamma ={16}_{-3}^{+3}\right)$). frank and rave find similar aspect ratios: $h_{\text{HWHM}}={0.023}_{-0.002}^{+0.002}$ and $h_{\text{HWHM}}={0.026}_{-0.006}^{+0.005}$, respectively. Lastly, our model with a radial double power law and vertical double Gaussian density distribution also cannot be ruled out (AIC σ = 1.196, ∆BIC = 15.567), corresponding to two measured scale heights: one with $h_{\text{HWHM}}={0.0009}_{-0.0003}^{+0.0005}$ containing ${78}_{-14}^{+9}\%$ of the dust mass, and the other with $h_{\text{HWHM}}={0.058}_{-0.017}^{+0.022}$ containing the rest of the mass.

Many of these measurements are comparable to other aspect ratio measurements of HD 197481; for example, Daley et al. (2019) model the 1.35 mm observations with a radial power law and vertical Gaussian, finding $h_{\text{HWHM}}={0.037}_{-0.005}^{+0.006}$. Vizgan et al. (2022) repeated this analysis, with a resulting measurement of $h_{\text{HWHM}}={0.029}_{-0.006}^{+0.005}$. They also analyze 450 µm observations of HD 197481, yielding a ratio of h_1.35mm/h₄₅₀ µ_m = 1.35 (h_HWHM ≈ 0.022 at 450 µm). Recently, Matrà et al. (2025) reanalyzed the 1.35 mm data assuming a radially Gaussian belt with a vertical Gaussian density distribution, finding $h_{\text{HWHM}}={0.020}_{-0.004}^{+0.004}$. Terrill et al. (2023) used a nonparametric approach to fit the radial structure of the 1.35 mm observations; assuming a vertical Gaussian profile, they find $h_{\text{HWHM}}={0.024}_{-0.002}^{+0.002}$. Similarly, Han et al. (2025) find $h_{\text{HWHM}}={0.026}_{-0.004}^{+0.004}$.

The variety of measured h_HWHM for HD 197481 suggest that the characterization of its vertical structure is highly dependent on the assumed density structure. Furthermore, while the fiducial measurement of h_HWHM – the smallest aspect ratio measured for this source – is resolved according to the criteria outlined in Sect. 4, it remains unclear whether the data resolution is sufficient to fully resolve such a small aspect ratio. The Band 6 observations of HD 197481 have a nominal resolution of θ_min = 1.86 au. At a reference radius of 36.3 au, an aspect ratio of h_HWHM = 0.0028 yields a scale height of H_HWHM ≈ 0.1 au. A vertical Lorentzian profile with a HWHM of 0.1 au would contain 95% of the disk material within 2.5 au centered on the disk midplane. This is less than a factor of two larger than θ_min; while such small spatial scales could be supported by the longest baseline observations, they cannot be resolved with two resolution elements. Given the spread of measured h_HWHM values under different density distribution assumptions, higher resolution follow-up observations may be helpful in breaking this degeneracy.

5 Discussion

5.1 Non-Gaussian vertical distributions as the new norm

Previous debris disk studies commonly assume that the vertical density structure follows a Gaussian distribution (e.g., Kennedy et al. 2018; Daley et al. 2019; Terrill et al. 2023; Matrà et al. 2025). This common assumption arises because self-stirring processes can generate a Rayleigh inclination distribution (Ida & Makino 1992; Lissauer & Stewart 1993), which corresponds to a Gaussian vertical density distribution (Matrà et al. 2019). Furthermore, distinguishing between different vertical parameterizations may not have been possible in studies analyzing lower-resolution data, particularly since the radial dust distribution is also less certain in these conditions.

However, high resolution observations of debris disks have enabled a more detailed exploration of the precise dust distributions. We find that most sources analyzed here (10/13) are significantly better fit with a thick-tailed vertical density distribution (i.e., either a Lorentzian or a double Gaussian) compared to a single Gaussian distribution. For each of the three remaining disks, vertical Lorentzian models are not ruled out at high confidence (HD 14055 with a radial double power law and vertical Lorentzian, AIC σ = 1.616 and ∆BIC = 3.376; HD 109573 with a radial double Gaussian and vertical Lorentzian, AIC σ = 0.090 and ∆BIC = 0.149; and HD 131488 with a radial double Gaussian and vertical Lorentzian, AIC σ = 1.255 and ∆BIC = 1.141; see Table A.1).

Some studies have already used non-Gaussian vertical structures for debris disks, in particular for HD 39060 (β Pic). For example, scattered light observations of HD 39060 have been fit using two vertical Lorentzians (e.g., Golimowski et al. 2006; Lagrange et al. 2012) and two vertical Gaussians (e.g., Apai et al. 2015), though the authors emphasize that these parameterizations were not physically motivated, but rather chosen by qualitatively assessing the suitability of different profiles. Similarly, Matrà et al. (2019) find that modeling HD 39060 using two vertical Gaussians provides a significantly better fit to the 1.33 mm dust observations than a single Gaussian, and use this parameterization to find marginal evidence of a radially varying vertical aspect ratio.

Now equipped with a larger sample of highly resolved debris disk observations, we have begun teasing out trends in the vertical dust distributions of disks. In our detailed modeling of 13 disks, we find evidence that the millimeter-sized dust typically exhibits a more gradual decline from the midplane density compared to the commonly used, steeper Gaussian profile. Rather, the best models tend to have vertical dust distributions with the heavy-tailed behavior seen in a Lorentzian distribution, or in the combination of a narrow Gaussian with a broader Gaussian. While we do not resolve two distinct scale height values in many disks at high confidence (only HD 10647 is identified as having a fiducial model with a vertical double Gaussian), we note that a Lorentzian distribution can be qualitatively similar to 1) the sum of two Gaussians with different widths and 2) an exponential distribution with ζ < 2, both of which are shown in Figure 4. Thus, some disks that are best fit by a Lorentzian may have multiple dynamical populations even if we are unable to measure multiple distinct scale heights at high confidence.

When comparing a model with fewer free parameters (e.g., one with a Gaussian vertical density distribution) to a more complex model like the vertical double Gaussian models presented here, the choice of model selection criteria may have a large impact on which model is identified as the best fit. Here, four of the 13 sources (HD 9672, HD 10647, HD 39060, and HD 161868) have a vertical double Gaussian model which is most strongly favored by the AIC, while no sources have a vertical double Gaussian model which is most strongly favored by the BIC. As the BIC penalizes model complexity more aggressively than the AIC, particularly for large datasets like ALMA observations (Eq. (10)), it may be useful to employ multiple selection criteria (e.g., both the AIC and BIC) to robustly evaluate model fit.

Multi-component fits have long been used for the classical Kuiper belt objects (KBOs), including double Gaussians, a Gaussian and a Lorentzian, and two von Mises–Fisher functions (Brown 2001; Kavelaars et al. 2009; Gulbis et al. 2010; Van Laerhoven et al. 2019; Malhotra & Roy 2023). Single Lorentzians have also been used to fit the classical KBOs (Adams et al. 2014). The inclinations and eccentricities of KBOs were likely inflated by Neptune’s outward migration during the early history of the Solar System (e.g., Malhotra 1995; Hahn & Malhotra 2005; Levison et al. 2008; Nesvorný 2015). Therefore, one possible interpretation of the apparent preference for Lorentzian or multi-component distributions in debris disks is that Neptune-like migration histories are common in these systems.

However, such structures can also arise from secular (i.e., long-term) planet-disk interactions without invoking planetary migration. For instance, Farhat et al. (2023) find a two-component inclination distribution emerges for the case of HD 106906, accounting for the central binary star as well as the known external planetary companion, while assuming a massless disk. Sefilian et al. (2025) consider a scenario involving a massive, back-reacting (but not self-gravitating) debris disk with an internal planetary perturber, finding that non-Rayleigh inclinations and non-Gaussian vertical profiles emerge naturally regardless of the disk mass. The range of disk processes and conditions which could generate multiple vertical populations is still an active area of study. For example, most previous works on planet-debris disk interactions neglect at least one or more of: the disk’s gravitational back-reaction and/or self-gravity, self-stirring processes, the possibility of multiple perturbers, and collisions within the disk. Gas-bearing debris disks should also be treated with care, as gas effects can also have a large impact on grain dynamics depending on the uncertain gas densities (though this may not have an observable effect on the vertical profile at millimeter wavelengths, Olofsson et al. 2022b). Regardless of the mechanisms driving the vertical grain distribution, non-Gaussianity appears to be nearly ubiquitous.

Furthermore, if the vertical distribution of solids changes with observing wavelength, then even a Gaussian vertical distribution of solids could produce a non-Gaussian emission profile (e.g., Terrill et al. 2023). This degeneracy could be broken by modeling well-resolved observations at multiple wavelengths. Currently, most debris disks which have been observed at multiple wavelengths lack the resolution for detailed vertical characterization at all wavelengths (Appendix C).

5.2 Constraints on the disk mass

The vertical thickness of debris disks can be influenced by a number of mechanisms, including early stirring in the protoplanetary disk phase (e.g., Walmswell et al. 2013; Booth & Clarke 2016), stirring by one or more planetary perturbers (e.g., Mustill & Wyatt 2009; Pearce & Wyatt 2014; Jílková & Portegies Zwart 2015; Dong et al. 2020; Chai et al. 2024, also see Sect. 5.5), or stellar flybys (e.g., Kenyon & Bromley 2002a; Lestrade et al. 2011; Moore et al. 2020; Batygin et al. 2020). Here, we primarily consider self-stirring of the disk, which is commonly invoked as a key mechanism driving disk scale heights (e.g., Ida 1990; Ida & Makino 1993; Kenyon & Bromley 2008; Kennedy & Wyatt 2010).

Under the assumption that the vertical thickness of disks is set by only self-stirring, measured scale height values can be used to place an upper limit on the mass of bodies stirring the disk. Following Daley et al. (2019), we defined the rms eccentricities and inclinations as $\overline{e}=\sqrt{\langle e^2\rangle }$ and $\overline{i}=\sqrt{\langle i^2\rangle }$, respectively, where ⟨e²⟩ and ⟨i²⟩ are the eccentricity and inclination dispersions of the grain population. We assume that grains are gravitationally excited isotropically in space such that ī = ē/2 (Inaba et al. 2001). At millimeter wavelengths, the inclination of a self-stirred disk is related to the observed aspect ratio by $\overline{i}~\sqrt{2}{h}_{\text{HWHM}}$ (Quillen et al. 2007). The relative velocities of particles are then given by $$\langle v_{\text{rel}}\rangle \approx v_{\text{Kep}}(r)\sqrt{{\overline{i}}^2+1.25{\overline{e}}^2}\approx 2\sqrt{3}{v}_{\text{Kep}}(r)h_{\text{HWHM}},$$(12) where v_Kep(r) is the Keplerian velocity at radius r (Lissauer & Stewart 1993; Wetherill & Stewart 1993; Wyatt & Dent 2002).

In the absence of significant damping, ⟨v_rel⟩ will be excited to roughly the escape velocity, v_esc, of the largest bodies in the system (which dominate the disk dynamics and drive most of the stirring, e.g., Pan & Schlichting 2012; Schlichting 2014), at which point collisions begin to dominate, limiting the growth of ⟨v_rel⟩. The escape velocity is given by $$v_{\text{esc}}=\sqrt{\frac{2\pi G{D}^2{\rho }_d}{3}}$$(13) for a spherical object of diameter D and bulk density ρ_d. By equating ⟨v_rel⟩ and v_esc, we can obtain lower limits on the size and mass of the largest bodies in the disk. We assumed a bulk density of ρ_d = 2.0 g/cm³, corresponding to the bulk densities of trans-Neptunian objects (Carry 2012). The minimum radii R_stir and masses M_stir of the largest bodies stirring each disk are presented in Table 4.

We also estimated the total mass of each disk using the collisional cascade model of Pan & Schlichting (2012). The model assumes a steady-state disk in which grains are stirred by close encounters with large bodies and damped by dynamical friction and collisions. Mass is conserved in all catastrophic collisions, and we assumed efficient collisional damping such that collisions with small bodies contribute significantly to the velocity damping (see Sect. 2 of Pan & Schlichting 2012). We adopted an eccentricity of 0.03 for all perturbing bodies. For each disk, we used Eq. (12) to estimate the typical random velocity of millimeter-sized grains within the disk, using the fiducial scale height values presented in Table 2. We used the stellar masses shown in Table 4, which were derived in Marino et al. (2026). Solving for the steady-state sizes and velocity distributions enabled the dust characteristics obtained from our parametric modeling (e.g., the millimeter dust mass and scale height) to be extrapolated to estimate the total amount of mass in the disk.

The derived masses, M_disk, for each disk are presented in Table 4 and plotted in Figure 5. The main panel shows the connection between h_HWHM and M_disk (we note that the positive correlation occurs due to the model specification, as h_HWHM is used to obtain the typical grain velocity, which directly affects the derived disk mass), while the top panel shows a histogram of the disk masses. Under half of the disks (5/13) contain less than the mass of Neptune, while the more massive disks reach up to several times the mass of Neptune. The color bar corresponds to the parametrically derived fractional widths (∆R/R) presented in Han et al. (2026). We find that belts which are vertically thicker and radially wider tend to be more massive. Furthermore, we see that the gas-bearing disks in the sample tend to be more massive than gas-poor disks with similar scale heights; however, gas is not considered in the collisional cascade model used to obtain our estimates on M_disk.

To verify whether the disk masses are large enough to stir the disks to their observed morphologies within the age of the system, we compute the stirring timescale T_stir using Eq. (28) in Krivov & Booth (2018). This calculation requires assumptions about the maximum planetesimal radius, belt location and fractional width, disk mass, stellar mass, bulk density, and fragmentation velocity v_frag. We use the values of R_stir, R_ref, M_disk, and M_⋆ in Table 4, and the ∆R/R values from Han et al. (2026). As with our stirring mass estimates, we assume a bulk density of ρ_d = 2.0 g/cm³. Because v_frag is not well-constrained, we adopt the nominal value of v_frag = 30 m/s used in Krivov & Booth (2018). Lastly, the equation includes the parameter γ, which ranges from 1 ≤ γ ≤ 2 and bounds the upper and lower limits on T_stir; we choose an intermediate γ = 1.5. The stirring timescales are shown in Table 4. We find that all but three disks (HD 32297, HD 131488, and HD 197481) could have been self-stirred within the age of the system. The disks that cannot be explained solely by self-stirring require an alternative mechanism, such as the planet-disk stirring scenarios discussed in Sect. 5.5.

There exist other methods of estimating M_disk. For example, Krivov & Wyatt (2021, hereafter K21) compute a minimum disk mass $M_{\text{disk}}^{\min }$ contained within the collisional cascade by extrapolating the observed mm-sized dust, assuming a power law size distribution with q = 3.7 up to some maximum size at the top of the cascade (s_max, estimated using observables). The Pan & Schlichting (2012) model employed here, on the other hand, estimates the mass of the largest bodies driving the stirring within the disk, which may not necessarily contribute directly to the collisional cascade (as the collision timescale becomes larger than the age of the system at these sizes). In this case, we assume that the dust mass is negligible, such that M_disk ≈ M_large.

These different modeling frameworks result in different mass estimates; for many of the disks, our estimates of M_disk are incompatible with the minimum masses reported in K21. However, much of this difference may stem from the choice of assumed millimeter dust mass, as the $M_{\text{disk}}^{\min }$ estimates from K21 use different millimeter dust masses than found here. To mitigate this effect, we rescaled the K21 estimates of $M_{\text{disk}}^{\min }$ to use the millimeter dust masses derived from our fiducial models (Table 4), finding that the estimates of $M_{\text{disk}}^{\min }$ are compatible with the estimates of M_disk found using the Pan & Schlichting (2012) model for four disks, and incompatible for the remaining three disks. These results are presented and discussed in more detail in Appendix D.

5.3 Comparison to other vertical measurements

The REASONS survey uniformly analyzed observations of 74 resolved debris disks (Matrà et al. 2025), many of which are also present in the ARKS sample. The modeling and analysis resulted in scale height measurements or constraints for many of the disks, including 11 out of the 13 disks in the sample presented here. In Figure 6, we compare the REASONS aspect ratios (circles) with the aspect ratios derived from the ARKS parametric modeling (stars).

We find that some sources with well-resolved aspect ratios from REASONS have significantly smaller (but still well resolved) values derived from the ARKS parametric modeling. The disks with the largest discrepancies tend to be some of the thinnest disks in the sample, according to the parametric constraints presented in this work (e.g., HD 197481, HD 109573, HD 32297). Some theoretical predictions suggest an increase in inclination dispersion (or scale height) with age (e.g., without considering damping and following $\sqrt{\langle i^2\rangle }\propto t^{\frac{14}{}}$ as in Ida & Makino 1993, assuming that stirrers form early in the disk lifetime). While this trend was not resolved with the REASONS sample, we find some evidence of an increase in ⟨i²⟩ with disk age in the ARKS scale height measurements.

The large differences between the REASONS and ARKS aspect ratios likely stems from differences in the model parametrizations, data resolution, and observing wavelength. Matrà et al. (2025) modeled each of the REASONS disks with a simple model that assumes the disk is both radially and vertically Gaussian. This had the advantage of enabling direct comparison between disks across a large sample of debris disks, including characterizing basic characteristics of the dust belts. However, the ARKS radial analysis reveals a diverse set of belt morphologies, very few of which are consistent with a Gaussian.

The work presented here further suggests that vertically, a Gaussian dust distribution may not be representative of many debris disks.

High resolution was not the focus of the REASONS survey, while the ARKS observing strategy focused specifically on resolving the belts well enough to conduct a more detailed analysis of each source (see Marino et al. 2026). On average, the resolution of the REASONS observations is several times poorer than the ARKS observations of a given target. In some cases, the angular resolution of the ARKS observations is nearly an order of magnitude higher (e.g., HD 32297 is resolved 0.″06 with ARKS and 0.″5 with REASONS). This resolution may have resulted in biased scale height measurements, even if they appear to be well-constrained using the simple Gaussian parametrization. Some of the aspect ratios recovered with the ARKS data still appear to be limited by the data resolution to some extent (Figure 2); follow-up observations of the thinnest disks at higher angular resolution may be useful to fully understand the dependence of measured scale height on resolution.

Lastly, variations in scale height with observing wavelength have already been observed (e.g., for AU Mic in ALMA Bands 6 and 9, Vizgan et al. 2022). The REASONS sample includes observations at 0.855 ≤ λ [mm] ≤ 1.36, so differences in the observing wavelengths could contribute to differences in measured aspect ratios. However, this is likely only a minor contribution, as it is not clear that such a small difference in wavelength would result in a measurable difference in h. A list of ARKS targets with measured aspect ratios (both from this work and previous studies) is available in Table C.1, along with the corresponding data resolution and observing wavelength.

Many ARKS targets have also been detected in scattered light; a detailed comparison of the scattered light and thermal emission of these disks is presented in Milli et al. (2026). Here, we briefly discuss the vertical profiles and scale heights from the scattered light observations. The vertical profiles of the disks modeled in scattered light are parametrized with a Gaussian profile of standard deviation r tan ψ linearly increasing with the radius r, where ψ is the vertical opening angle of the disk. The HWHM of such a vertical profile is therefore $\sqrt{2\ln 2}\times r\tan \psi$ and the corresponding aspect ratio is $\sqrt{2\ln 2}\tan \psi$.

This quantity is plotted in Figure 7 together with the h_HWHM derived from the fiducial ARKS parametric models (Table 2). Systems are plotted with increasing millimeter scale height, and gas-bearing disks are highlighted with a light blue background. We find that most disks have a similar scale height or are thicker in scattered light compared to thermal imaging. HD 39060 has a substantially larger scale height in both thermal emission and scattered light compared to the rest of the targets.

Theoretical works predict the gas-bearing disks to have a smaller scale heights caused by dust vertical settling (Olofsson et al. 2022b). For second-generation gas present at a low level (≲ 0.1M_⊕, depending on the dust mass), this settling is expected to be most visible in scattered light while the larger grains seen in thermal images are less affected. On the other hand, we expect larger quantities of gas (≳ 0.1M_⊕) to cause significant damping in both millimeter and scattered light observations. In Figure 7, no clear trend emerges from the comparison of the gas-bearing disks shown shaded in blue with the disks containing no detectable amount of CO.

These comparisons come with several caveats. First, not all disks are well resolved vertically in scattered light. For instance, this is the case for HD 32297, as mentioned in Olofsson et al. (2022). Second, depending on the dataset used to model the disk vertical distribution in scattered light (the spectral filter for HD 15115 or the choice between total or polarized light for most other disks, see details in Milli et al. 2026), the scale height may be different, mostly because the signal-to-noise of the scattered light observations is not the same.

5.4 Theoretical implications

One way to probe which dynamical effects are at work in these disks is to consider the ratio ⟨e_rms⟩/⟨i_rms⟩. In self-stirred disks an energy equipartition between vertical and radial motion is expected, such that ⟨e_rms⟩/⟨i_rms⟩ ∼ 2 (Ida & Makino 1992). If a disk is stirred by a (nearly) co-planar planet, eccentricities will be excited to a greater degree than the inclinations, and we may expect that ⟨e_rms⟩/⟨i_rms⟩ > 1 (Pearce et al. 2024). Alternatively, if a disk is pre-stirred (e.g., if the dynamical excitation is set during the protoplanetary phase), and if the eccentricity and the inclination distributions subsequently evolve due to fragmentation-dominated collisional evolution, we may expect that the ratio decreases to ⟨e_rms⟩/⟨i_rms⟩ < 1 with inefficient collisional damping (Jankovic et al. 2024). The vertical and radial constraints from this work and Han et al. (2026) enable observational constraints to be placed on this ratio.

We used the upper limit estimates of ⟨e_rms⟩ presented in Han et al. (2026), which assume a self-stirred disk (Marino 2021; Rafikov 2023). For disks with no reported ⟨e_rms⟩, we estimated the upper limit of ⟨e_rms⟩ using the belt fractional width ∆R/R as in Marino (2021), $$\langle e_{\text{rms}}\rangle \le 0.6\Delta R/R.$$(14)

For these sources, we adopted the median fractional widths reported in Han et al. (2026). We estimated ⟨i_rms⟩ as described in Sect. 5.2. We find that the majority of disks have an upper limit of ⟨e_rms⟩/⟨i_rms⟩ > 2, while four (HD 15115, HD 39060, HD 76582, and HD 161868) have ⟨e_rms⟩/⟨i_rms⟩ ≲ 2 (Figure 8). While none of the disks have an upper limit of ⟨e_rms⟩/⟨i_rms⟩ < 1, the true ratio could be lower, so we cannot rule out that some disks in the sample may be dynamically dominated by catastrophic collisions.

We emphasize that the values of ⟨e_rms⟩ and ⟨i_rms⟩ shown in Figure 8 should be taken as estimates. Because ⟨e_rms⟩ is derived from the radial distribution of material in the disk, there is a degeneracy between the semi-major axis distribution and the eccentricity distribution in the disk. Thus, the true values of ⟨e_rms⟩ could be smaller. Similarly, while the errors on ⟨i_rms⟩ shown in Figure 8 stem only from uncertainties in the measured scale height, there is additional uncertainty stemming from the approximations and assumptions used to compute ⟨i_rms⟩. We calculate ⟨i_rms⟩ using the approximation from Quillen et al. (2007), such that $\langle i_{\text{rms}}\rangle ~\sqrt{2}{h}_{\text{HWHM}}$. Matrà et al. (2019), on the other hand, assume a Rayleigh distribution of inclinations in a self-stirred disk to derive $\langle i_{\text{rms}}\rangle ~\sqrt{2}{h}_{\sigma }$. We note that h_σ is approximately 85% of h_HWHM. However, the non-Gaussianity of most vertical profiles means that neither of these equations are strictly true for this application.

Another way to compare different theoretical models to the observations is to consider how h varies with wavelength. If a disk is fragmentation-dominated with no significant collisional damping, h is expected not to vary with wavelength, in viscously stirred disks (Kenyon & Bromley 2002b, 2004, 2008; Kennedy & Wyatt 2010) as well as in pre-stirred disks (Jankovic et al. 2024). On the other hand, if h increases with wavelength, it could be a sign of efficient collisional damping in the system (Pan & Schlichting 2012).

We compiled ALMA measurements of h_HWHM, shown in Table C.1 and summarized in Figure C.1. A few systems show a consistent increase in measured h_HWHM: HD 15115, HD 32297, HD 61005, and HD 131488. The remaining disks with multiwavelength constraints (HD 9672, HD 10647, HD 14055, HD 107146, HD 161868, and HD 197481⁸) appear more consistent with an aspect ratio that remains approximately constant with wavelength. The remaining disks do not have measurements in multiple ALMA bands.

The predictions come with the caveat that results are difficult to directly compare if the h measurements were obtained with different methods or with different data resolutions. Unfortunately, this makes it impossible to draw firm conclusions with the current data and analyses. Follow-up measurements at similar resolutions to the ARKS observations could help to confirm or rule out these tentative trends.

5.5 Constraints on planetary companions

In this work, we have assumed that disks have aspect ratios which are driven by internal processes. However, there is a rich literature exploring alternative scenarios in which planetary (or stellar) companions drive the required dynamical excitation. In these cases, the properties of the hypothetical planet(s) and the disk-to-planet mass ratio (M_disk/m_pl) set the disk particles’ eccentricities and inclinations. For example, low mass disks may have observable differences in their vertical structure depending on whether the planet’s orbit is inclined with respect to the disk: a coplanar planet would excite the disk radially rather than vertically (Pearce & Wyatt 2014). On the other hand, initially misaligned, low-inclination planets can generate warps that propagate through the disk over time, causing the disk to settle into a vertically inflated configuration. This can occur regardless of whether the planet lies interior (Mouillet et al. 1997; Dawson et al. 2011; Pearce & Wyatt 2014; Smallwood 2023) or exterior to the disk (Nesvold et al. 2017; Farhat et al. 2023), assuming the dynamics are not significantly affected by the disk gravity (typically requiring M_disk/m_pl ≪ 1; Sefilian et al. 2025). Assuming a given configuration, one can then use the measured vertical scale heights to infer the planet properties.

For relatively massive disks, the collective gravity of the disk may be significant enough to resist planetary perturbations, maintaining a low dispersion in both radial (Sefilian 2024) and vertical directions (Poblete et al. 2023; Sefilian et al. 2025). Alternatively, a massive disk could give rise to secular-inclination resonances and localized warps (Sefilian et al. 2025). This significantly complicates the interpretation of planetary parameters inferred from massless disk models, which typically result in artificially low planetary masses, particularly when the planet lies close to the disk (Sefilian et al. 2025). A full exploration of planet-driven dynamical excitation scenarios is outside the scope of this study. However, a future study will provide a detailed investigation on planetary inference in the ARKS sample, building on the vertical aspect ratio measurements presented in this work (Jankovic et al. in prep).

5.6 Limitations

We have presented the modeling and analysis of the vertical aspect ratios and dust distributions of the highly inclined ARKS targets. However, the analyses presented here have several limitations which could impact how the results are interpreted. First, we consider only axisymmetric disks with vertical distributions that peak at the disk midplane. The axisymmetric modeling could obscure more complex structures both radially and vertically. For example, the presence of a warp could bias the axisymmetric models to larger scale heights in order to account for the warped material sitting further from the midplane. An analysis of the asymmetric features in the ARKS sample is presented in Lovell et al. (2026).

Second, in this work we assume that aspect ratio is constant with radius. Using a constant aspect ratio is beneficial for enabling direct comparison with the majority of previous vertical structure measurements (see Table C.1), as well as for minimizing the number of free parameters in our models. However, recent work has called the validity of this assumption into question, as more disks with evidence of a variable h are identified. For example, Matrà et al. (2019) parametrized the aspect ratio of HD 39060 as h ∝ r^β−1. In this parametrization, β = 1 corresponds to a constant aspect ratio (as in our Eq. (1)), β > 1 to a flared disk, and β = 0 to a constant H. The model fit yielded β = 0.7 ± 0.2, which is consistent with a constant h(r), but also provides marginal evidence (at < 2σ) of a slowly increasing scale height resulting in a radially varying h(r). This, however, did not result in a significantly improved fit compared to models assuming a constant h. Similarly, Terrill et al. (2023) allowed for a varying aspect ratio with radius when fitting HD 9672 (49 Ceti), and found a loosely constrained value of $\beta ={0.79}_{-0.35}^{+0.29}$.

More recently, Han et al. (2025) fit H(r) nonparametrically, finding several disks that are not consistent with a constant aspect ratio. Some ARKS targets were included in this study: HD 14055 was found to be consistent with a constant aspect ratio, while HD 32297, HD 61005, and HD 197481 were found to have a steeper H(r) (scale height increasing more sharply with radius). Interestingly, HD 15115 was found to have a decreasing H(r), suggesting that the inner ring may be thicker than the outer ring.

From a theoretical perspective, Sefilian et al. (2025) find that in debris disks perturbed by inner planets, the aspect ratio can follow h(r) ∝ r^−7/2, with the extent of this behavior depending on the disk-to-planet mass ratio, M_disk/Mplanet. This effect arises from the disk’s ability to mitigate planetary perturbations (see also Sefilian 2024). When M_disk/M_planet ≳ 1, this trend holds across the entire disk; for M_disk/M_planet ≲ 1, it applies only in the disk’s outer regions, while the inner parts maintain h(r) ≈ const. This suggests that planet-disk interactions can produce aspect ratios that vary and decrease with radius. The work presented here only considers the case of a radially constant h, but follow-up studies with nonconstant h would be beneficial to fully understand how common scale height variability is in debris disks.

Lastly, the specific radial and vertical density parametrizations can impact the measured scale heights. A poor choice of radial profile can result in biased vertical results, (and vice versa). This may explain some of the differences between h_HWHM measurements derived from REASONS and ARKS (Figure 6). Similar discrepancies can be seen in Figure C.1, which shows all of the measured h_HWHM listed in Table C.1. We have attempted to address this issue by leveraging the extensive radial structure analysis presented in Han et al. (2026) to ensure that we are using the most appropriate parametrizations for the radial dust distributions before extending the exploration to the vertical dust structures. However, it is possible that the radial and/or vertical density distributions adopted in this work do not fully capture the morphology of certain sources.

6 Conclusions

In this paper we have presented the detailed vertical structure analysis of the ARKS debris disks, including constraints from parametric modeling, frank, and rave. We successfully measured the scale heights for the 13 most highly inclined targets in the ARKS sample using all three methods, and we recovered scale height constraints for a handful of moderately inclined ARKS targets with frank and rave (Table 3). Our key results are summarized here:

• Many disks are quite thin (h_HWHM ≲ 0.01), while others host millimeter-sized dust grains much farther from the midplane (h_HWHM ≳ 0.1);

• Lorentzian (or otherwise thick-tailed) vertical dust distributions are common, suggesting that systems may commonly host multiple dynamical populations within the disk. This could hint at Neptune-like migration histories in these systems or could be a remnant of long-term planet-disk interactions in general;

• We used our aspect ratio measurements to estimate the total mass of solids in the disk, and we placed limits on the mass of the bodies that could be gravitationally stirring the disks (Table 4, Figure 5). We find that fewer than half of the disks contain less mass than Neptune and that the most massive disks also have the broadest belts radially;

• We measure much smaller aspect ratios (thinner disks) than comparable measurements made with lower-resolution data at the same observing wavelength. As a result, we may for the first time be observing the theoretically predicted increase in scale height with disk age. However, because the disks are so thin, follow-up observations at higher resolution could help us better understand this trend;

• Aspect ratios obtained from scattered light observations tend to be similar to or larger than the corresponding ALMA millimeter aspect ratios (Figure 7). We do not observe a clear trend when comparing the gas-bearing debris disks against those with no detected gas, though several of the disks are not well resolved vertically in scattered light.

These measurements provide new empirical constraints on the vertical structures of debris disks and the range of dynamical processes that are likely to sculpt them. Future high resolution and multiwavelength observations will help to further refine these measurements, strengthening our population level understanding of the debris disk morphologies and dynamics.

Data availability

The ARKS data used in this paper and others can be found in the ARKS dataverse. The data products produced by this work can be found at https://doi.org/10.7910/DVN/6NNKLY. For more information, visit arkslp.org.

Acknowledgements

This paper makes use of the following ALMA data: ADS/JAO.ALMA# 2022.1.00338.L, 2012.1.00142.S, 2012.1.00198.S, 2015.1.01260.S, 2016.1.00104.S, 2016.1.00195.S, 2016.1.00907.S, 2017.1.00167.S, 2017.1.00825.S, 2018.1.01222.S and 2019.1.00189.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The project leading to this publication has received support from ORP, that is funded by the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 101004719 [ORP]. We are grateful for the help of the UK node of the European ARC in answering our questions and producing calibrated measurement sets. This research used the Canadian Advanced Network For Astronomy Research (CANFAR) operated in partnership by the Canadian Astronomy Data Centre and The Digital Research Alliance of Canada with support from the National Research Council of Canada the Canadian Space Agency, CANARIE and the Canadian Foundation for Innovation. Support for BZ was provided by The Brinson Foundation. This material is based upon work supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE 2140743. AMH, AJF, and EM acknowledges support from the National Science Foundation under Grant No. AST-2307920. EM acknowledges support from the NASA CT Space Grant. MP acknowledges support from NASA grant 80NSSC21K1334. JM acknowledges funding from the Agence Nationale de la Recherche through the DDISK project (grant No. ANR-21-CE31-0015) and from the PNP (French National Planetology Program) through the EPOPEE project. TDP is supported by a UKRI Stephen Hawking Fellowship and a Warwick Prize Fellowship, the latter made possible by a generous philanthropic donation. A.A.S. is supported by the Heising-Simons Foundation through a 51 Pegasi b Fellowship. MB acknowledges funding from the Agence Nationale de la Recherche through the DDISK project (grant No. ANR-21-CE31-0015). AB acknowledges research support by the Irish Research Council under grant GOIPG/2022/1895. CdB acknowledges support from the Spanish Ministerio de Ciencia, Innovación y Universidades (MICIU) and the European Regional Development Fund (ERDF) under reference PID2023-153342NB-I00/10.13039/501100011033, from the Beatriz Galindo Senior Fellowship BG22/00166 funded by the MICIU, and the support from the Universidad de La Laguna (ULL) and the Consejería de Economía, Conocimiento y Empleo of the Gobierno de Canarias. EC acknowledges support from NASA STScI grant HST-AR-16608.001-A and the Simons Foundation. S.E. is supported by the National Aeronautics and Space Administration through the Exoplanet Research Program (Grant No. 80NSSC23K0288, PI: Faramaz). MRJ acknowledges support from the European Union’s Horizon Europe Programme under the Marie Sklodowska-Curie grant agreement no. 101064124 and funding provided by the Institute of Physics Belgrade, through the grant by the Ministry of Science, Technological Development, and Innovations of the Republic of Serbia. This work was also supported by the NKFIH NKKP grant ADVANCED 149943 and the NKFIH excellence grant TKP2021-NKTA-64. Project no.149943 has been implemented with the support provided by the Ministry of Culture and Innovation of Hungary from the National Research, Development and Innovation Fund, financed under the NKKP ADVANCED funding scheme. JBL acknowledges the Smithsonian Institute for funding via a Submillimeter Array (SMA) Fellowship, and the North American ALMA Science Center (NAASC) for funding via an ALMA Ambassadorship. SMM acknowledges funding by the European Union through the E-BEANS ERC project (grant number 100117693), and by the Irish research Council (IRC) under grant number IRCLA- 2022-3788. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. SM acknowledges funding by the Royal Society through a Royal Society University Research Fellowship (URF-R1-221669) and the European Union through the FEED ERC project (grant number 101162711). JPM acknowledges research support by the National Science and Technology Council of Taiwan under grant NSTC 112-2112-M-001-032-MY3. LM acknowledges funding by the European Union through the E-BEANS ERC project (grant number 100117693), and by the Irish research Council (IRC) under grant number IRCLA- 2022-3788. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. SP acknowledges support from FONDECYT Regular 1231663 and ANID – Millennium Science Initiative Program – Center Code NCN2024_001. PW acknowledges support from FONDECYT grant 3220399 and ANID – Millennium Science Initiative Program – Center Code NCN2024_001.

References

---

Appendix A Summary of parametric modeling

Appendix B frank and rave vertical inference

Appendix C Summary of new and archival aspect ratio measurements

Appendix D Comparison of disk mass estimates

All Tables

All Figures