A Uniform Analysis of Debris Disks with the Gemini Planet Imager II: Constraints on Dust Density Distribution Using Empirically-Informed Scattering Phase Functions

Spatially-resolved images of debris disks are necessary to determine disk morphological properties and the scattering phase function (SPF) which quantifies the brightness of scattered light as a function of phase angle. Current high-contrast imaging instruments have successfully resolved several dozens of debris disks around other stars, but few studies have investigated trends in the scattered-light, resolved population of debris disks in a uniform and consistent manner. We have combined Karhunen-Loeve Image Projection (KLIP) with radiative-transfer disk forward modeling in order to obtain the highest quality image reductions and constrain disk morphological properties of eight debris disks imaged by the Gemini Planet Imager at H-band with a consistent and uniformly-applied approach. In describing the scattering properties of our models, we assume a common SPF informed from solar system dust scattering measurements and apply it to all systems. We identify a diverse range of dust density properties among the sample, including critical radius, radial width, and vertical width. We also identify radially narrow and vertically extended disks that may have resulted from substellar companion perturbations, along with a tentative positive trend in disk eccentricity with relative disk width. We also find that using a common SPF can achieve reasonable model fits for disks that are axisymmetric and asymmetric when fitting models to each side of the disk independently, suggesting that scattering behavior from debris disks may be similar to Solar System dust.


INTRODUCTION
Since the first indirect detection of a circumstellar debris disk around Vega (Aumann et al. 1984), numerous studies have investigated the properties and structures of these dusty systems across a wide range of wavelength regimes and angular resolutions.Resolved disks have been observed in sizes ranging from tens to hundreds of AU in diameter (e.g.Schneider et al. 1999).Due to the short timescale processes of Poynting-Robertson drag and radiation pressure, debris disks must continually produce dust to sustain their large and extended structures.These dust-replenishing properties may indicate the presence of ongoing planet formation and/or dynamical influence, such as the collisional grinding of planetesimals (Backman & Paresce 1993) or collisions of planets (Cameron 1997).
Direct imaging provides significant insight into the architecture of a debris disk system.While spectral energy distributions (SEDs) can place some constraints on the radial extent, dust mass, and composition of system dust, debris disk images can more directly constrain overall properties of disk morphology such as the radial extents of dust and planetesimal belts (Esposito et al. 2020).Substellar companions can directly influence the shapes of these belts, inducing features such as gaps, warps, and clumps that can be identified from resolved imaging.In addition to informing studies of disk dynamics and evolution (e.g. Lee & Chiang 2016), resolved debris disk images can also constrain the gravitational interactions between planets and disks (e.g., Liou & Zook 1999;Kuchner & Holman 2003;Quillen & Faber 2006;Wyatt 2006).
Direct imaging studies at the scales of 0. ′′ 01-1 ′′ have been explored at optical and near-IR wavelengths with instruments such as the Space Telscope Imaging Spectrograph (STIS; e.g.Schneider et al. 2014Schneider et al. , 2016) ) and ground-based adaptive optics (e.g.Nasmyth Adaptive Optics System-COude Near Infrared CamerA; NACO; Lenzen et al. 2003;Coronagraphic High Angular Resolution Imaging Spectrograph;CHARIS;Groff et al. 2015).At near-IR wavelengths, (sub)micron-sized dust grains are expected to scatter light from the host stars they surround.These observations are technically challenging, as the disk brightness from dust scattering is typically ∼ 10 6 times fainter than the host star brightness.Current generation direct imaging instruments such as the Gemini Planet Imager (GPI; ★ E-mail: jrhom@asu.eduMacintosh et al. 2014) and the Spectro-Polarimetric High-contrast Exoplanet REsearch (SPHERE; Beuzit et al. 2019) utilize adaptive optics (AO) and coronagraphy to increase the contrast by minimizing residual atmospheric turbulence (e.g.Poyneer et al. 2014) and blocking light from the host star (e.g.Soummer et al. 2009).Increased sensitivity and spatial resolution allow for easier identification and characterization of subtle asymmetric features such as eccentricity, warps, and clumps that could be attributed to the presence of substellar companions or ongoing planet formation.
Since its first light in 2014, GPI has spatially resolved ∼ 15 debris disks for the first time (e.g.Esposito et al. 2018;Kalas et al. 2015;Hom et al. 2020;Hung et al. 2015;Millar-Blanchaer et al. 2016;Currie et al. 2015) at inner working angles and resolutions not accessible by previous generation instruments.Combined with an additional sample of newly-resolved scattered-light debris disks from SPHERE and other instruments (e.g.Thalmann et al. 2013;Wahhaj et al. 2016;Engler et al. 2018Engler et al. , 2020;;Bonnefoy et al. 2021) and previously resolved systems (e.g.Schneider et al. 1999Schneider et al. , 2005;;Padgett & Stapelfeldt 2016;Liu 2004;Choquet et al. 2016;Soummer et al. 2014;Hines et al. 2007;Kalas et al. 2007), the ensemble of scattered-light imaged disks provides a rich and diverse sample to investigate overall trends in system architectures of young planetary systems and the properties of dust grains around stars of different spectral types and ages.
Despite this new large sample of scattered-light resolved debris disks, few group/population studies have been conducted.Ren et al. (2023) investigated debris disk color through HST observations, identifying a predominantly blue color suggesting higher scattering efficiency at shorter wavelengths.Esposito et al. (2020) first reported on the statistics and properties of the 29 circumstellar disks imaged as a part of the Gemini Planet Imager Exoplanet Survey (GS-2015A-Q-500, PI B. Macintosh).While the study reported and characterized the geometries of a few newly-resolved disks, the majority of properties of the sample were collated from previous investigations which applied different analysis approaches of individual systems, including studies utilizing observations from other instruments.These individual investigations often utilize unique approaches to data reduction and disk characterization, making direct comparisons between observations of the same system with different instruments difficult.Data reduction, particularly stellar point spread function (PSF) subtraction, can substantially affect the apparent structural appearance of a disk system (Milli et al. 2012), biasing morphological characterization efforts.Different approaches to disk modeling lead to results that are not directly comparable and sometimes biased, depending on the assumptions.For example, some studies employ ellipse fitting to assess inclination and position angle, which assumes the disk is an infinitely narrow ring (e.g.Crotts et al. 2024).The treatment of disk scattering properties is also approached with different methods.Henyey-Greenstein (HG) functions (Henyey & Greenstein 1941) have long been used to described the SPF of dust populations in the solar system (Hong 1985).Although they are not based on any physical scattering theories, a linear combination of a few HG functions can be used to approximate a wide range of scattering patterns using only a few parameters.Other studies that opt to model grain properties more robustly are often limited in terms of complexity, having to assume a uniform shape for all dust grains (e.g.Mie theory, Mie 1908) and a limited number of grain species.Furthermore, the parameterization of grain properties has been observed to provide unrealistic results (e.g.Duchêne et al. 2020), where constrained dust properties may be unphysical or contradictory.
In the first publication of this series (Crotts et al. 2024), empirical measurements of a set of disk morphological and brightness characteristics were conducted in a uniform manner to GPI polarized intensity data.In this study, we applied a uniform data reduction and radiative-transfer modeling approach to facilitate more direct comparisons between morphological properties of a sample of debris disk systems in total intensity light, investigating a separate regime of scattered-light resolved structure from Crotts et al. (2024).This approach applied a commonly-used dust density distribution function and explores consistent parameter spaces for all debris disk targets.Rather than modeling grain properties, which can be highly biased depending on the underlying assumptions and free parameters explored, we chose to utilize two empirical scattering phase functions that are not reliant on an HG formalism or underlying assumptions of grain properties.In this work, we present our modeling results of eight debris disk targets imaged by GPI.In §2, we describe the samples from which our observations are derived and the additional criteria employed to construct our final sample.In §3, we describe the properties of the observations of our target sample.In §4, we describe the data reduction approach for all datasets in our sample.In §5, we describe our forward modeling setup.In §6, we present the images and constrained parameters of our model analysis.In §7, we describe the interpretations of our results in modeling and understanding debris disk systems.In §8, we summarize our findings and discuss future implications for the results of our analysis.

TARGET SAMPLE
Selected targets in this analysis originate from three distinct observational programs.The majority of the data originate from the Gemini Planet Imager Exoplanet Survey (GPIES) campaign (GS-2015A-Q-500, PI B. Macintosh), a 600-star direct imaging survey utilizing the Gemini Planet Imager (GPI; Macintosh et al. 2014) in the spectral integral field unit (IFU) mode at H-band.A subset of this campaign was dedicated to observing both previously known scattered-light and/or thermal emission-resolved circumstellar disks as well as stars with notable IR-excesses exceeding 10 −5 at H-band wavelengths.Two other programs, "Debris Characterization in Exoplanetary Systems" (PI C. Chen; GS-2016A-LP-6) and "Does the HR 4796 Debris Disk Contain Icy Grains?" (PI C. Chen; GS-2015A-Q-27) were used for the remaining observational sequences of this sample and constitute a Gemini Large and Long Program (LLP).These two programs added both J-and K1-band wavelengths to the sample.All of these programs utilized the capabilities of GPI to image debris disks at higher contrast (∼ 10 −6 ) and smaller angular separations (tens of mas) than were possible with previous generation instruments.
From the extensive GPIES Disk and GPI-LLP Disk samples, eight targets were selected.These targets satisfied three main criteria: (1) data from a GPI spectral mode dataset, (2) detection in total intensity light, and (3) average signal-to-noise ratio (SNR) must exceed ∼5 per pixel.As an IFU, GPI only operates in two observational modes: spectral and polarimetric.As described in §5.2, a simulated PSF core unique to each dataset is necessary for our analysis, as PSF structure and intensity can vary with observing conditions and target brightness.In GPI datasets, this PSF core is modeled from photometric measurements of satellite spots created by fiducial images of the host star superimposed on the GPI pupil apodizer in each science frame.In spectral mode datasets, these spots have well-defined shapes and flux ratios.Although GPI has resolved an extensive sample of disks in polarized intensity (Esposito et al. 2020;Crotts et al. 2024), they are not considered in this analysis.In polarimetric datasets, the satellite spots necessary for our PSF core modeling routine are elongated and cannot be used to accurately generate a model PSF core for disk forward modeling.Finally, although polarized intensity phase functions contain notable features among scattered-light resolved disks, no common trends have been identified to warrant a model analysis of multiple systems with a common polarized intensity phase function, and studies of cometary dust have revealed different shapes to polarized intensity phase functions depending on the composition of grains (Frattin et al. 2019).Additionally, polarizability curves which describe the nature of polarization throughout a resolved scatteredlight disk can vary distinctly between systems and would introduce an additional parameter space of investigation that is not relevant to studying morphological properties of disks (e.g.Hadamcik et al. 2007;Hadamcik & Levasseur-Regourd 2003).Therefore, our model methodology approach cannot be applied to polarized intensity images.
All of the targets considered for this sample have H-band images taken as a part of the GPIES campaign.From the GPI-LLP Disk Sample, only K1-band observations of HD 32297 and HR 4796A were considered due to the SNR threshold we applied.Summary information regarding our sample targets is shown in Table 1.

OBSERVATIONS
Observations in H-and K1-band were conducted over several semesters at Gemini-South with GPI as a part of the GPIES campaign and the GPI-LLP Disk program, summarized in Table 2.All observations were conducted in GPI's spectral mode and covered high field rotation (ΔPA∼ 17 − 80 • ) for increased effectiveness in utilizing angular differential imaging (ADI; Lafrenière et al. 2007) for PSF subtraction.Integration times were selected such that detector readout noise did not exceed signal from the disk while also avoiding angular smearing from rotation and saturation of speckles.K1-band observations can tolerate longer exposure sequences (∼90 s) due to PSF speckles being fainter at longer wavelengths.The GPI target prioritization scheme also emphasized observing targets near transit to further achieve high field rotation.

DATA REDUCTION
The GPI Data Reduction pipeline (DRP, Perrin et al. 2014;Wang et al. 2018) was used for reducing all science data from GPI.The pipeline performed dark subtraction, correlated noise cleaning, and bad pixel correction for all raw data.For spectral mode data in particular, flexure correction for satellite spots was also performed.Before a spectral sequence was observed, an Ar lamp exposure was collected for wavelength calibration.Geometric distortion was also corrected for in all data cubes and smoothed with a Gaussian kernel ( = 1 pixel).The DRP was responsible for assembling the integral field spectrograph (IFS) spectral data cube and determining the location of the satellite spots (necessary for performing the DiskFM analysis, see §5).For K1-band observations only, thermal/sky background images were subtracted.The science images were also destriped, and some spectral slices within the IFS spectral datacube were not considered for the analysis if the SNR of the satellite spots was low due to high thermal noise at K-band wavelengths.Finally, following the approach in Wang et al. (2014), the location of the occulted primary star in each science frame was determined by performing a least-squares fit to all visible satellite spot positions to a precision of 0.7 mas.Science frames were then shifted to align with the star center located at the center of every image.The science frames were also all rotated so that north pointed upward in the image and east pointed to the left.
Spectral data cubes were further PSF-subtracted with the pyKLIP (Wang et al. 2015) implementation of the Karhunen-Loève Image Projection (KLIP) algorithm (Soummer et al. 2012) in combination with ADI.ADI takes advantage of the fact that instrumental PSF artifacts do not rotate with the FOV in pupil-stabilized observations.Any astrophysical object within the FOV will rotate with respect to the center of the image throughout the observational sequence.A model PSF halo pattern can then be generated from the rotating frames without including astrophysical signal (except in the case of extended structures such as disks) and subsequently subtracted from images throughout the sequence.KLIP expands upon this approach by performing principal component analysis of PSF features, achieving more robust results than ADI alone.For this analysis, 3-7 Karhunen-Loève modes and 4 − 6 • of minimum rotation (  in Lafrenière et al. 2007) between science frames were used for PSF reconstruction.The choice in the number of KL-modes and minimum rotation angles is informed by the inclination of the system, maximizing the average disk SNR, and minimizing overlap of the disk between science frames to mitigate self-subtraction (Milli et al. 2012).The analysis was performed globally across a whole GPI image and not subdivided into concentric annuli or subsections, facilitating the creation of continuous and smooth reduced images.The full details of the data reduction parameters are shown in Table 3.The Karhunen-Loève basis vectors were saved and projected onto disk forward models for later analysis, as described in §5.

MODEL METHODOLOGY
As high contrast images of debris disks reduced with ADI often suffer from severe self-subtraction (Milli et al. 2012), forward-modeling is often necessary to constrain the properties of a system.The key steps in the forward modeling methodology employed in this study are shown in Figure 1.To assess morphological properties, we adopt a forward modeling approach using DiskFM (Mazoyer et al. 2020) for comparison with the science data.Disk parameters are constrained using an iterative MCMC process with the affine-invariant sampler emcee Python package (Foreman-Mackey et al. 2013), generating hundreds of thousands of disk models per target.

Disk Model
The disk model utilized is based on a built-in module of DiskFM and is originally described in Millar-Blanchaer et al. (2015) and Millar-Blanchaer et al. (2016), with two major modifications made, related to the dust density distribution and dust scattering properties.
The updated surface density profile (, ) follows a smoothly connected two power-law structure described in more detail in Augereau et al. (1999) and given in Equation 1:  where () is given as: where  is the radial distance from the host star and  C is the critical radius where a transition between two power law density regimes occurs.The indices of these power law density regimes are given as  in > 0 and  out < 0 for the inner and outer regions of the disk respectively. (, ) is given as where  is the distance from the disk midplane,  vert dictates the shape of the vertical density distribution, and ℎ() is the height above the disk midplane as a function of .The Augereau et al. (1999) dust density profile was selected as it is more commonly used in dust density distribution analyses and therefore facilitates more consistent comparisons of our results to previous studies.For our analysis, we set  vert = 2 for a Gaussian vertical profile and ℎ() is given by   × where   is the constant aspect ratio ℎ/, assuming a "bow-tie" shape for every disk.analysis, as we only seek to understand the overall geometry of each system and not the dust mass; a constant brightness scaling factor is used as a free parameter to achieve low residual model fits to the data.
The second significant modification is the treatment of light scattering properties of the model.In the original model framework, a choice of a one, two, or three component HG function is used for representing all scattering properties.In this study, scattering properties of models are given by an empirically-informed scattering phase function interpolated with a 3rd order spline function.
One of the premises of this study is the observation by Hughes et al. ( 2018) that scattering in most debris disks appears to follow a similar SPF, which itself matches qualitatively that observed in a number of solar system objects.Based on the apparent similarities of the few measured SPFs, we generated a generic scattering phase function in the following manner.First, we gathered the SPF from Saturn's D68 and G rings (Hedman & Stark 2015), Jupiter's ring (Throop et al. 2004), and multiple comets (Hanner & Newburn 1989;Schleicher et al. 1998;Moreno et al. 2012;Hui 2013).To avoid one particular dataset biasing the final SPF due to more intense and/or wider sampling of scattering angles, we rebinned each individual SPF to a common 3 • sampling.We then renormalized each SPF to the most completely sampled SPF-that of Saturn's rings-using the median of their ratio over all overlapping bins, and averaged all resulting SPFs.Incidentally, the SPF of the Saturn D68 has been shown to have a similar shape to the measured SPF of HD 114082 from Engler et al. (2023).Finally, we fitted a 9th-degree polynomial (the lowest order that avoids overfitting) to these datapoints and the result is taken as the "generic SPF".The generic SPF used in this study is shown in Figure 2, overplotted with measured SPFs from solar system dust environments and the markedly distinct SPF measured from SPHERE images of HR 4796A by Milli et al. (2017).
We note that most of the SPFs used in this process were observed in the optical (the main exception being the  band compilation from Hanner & Newburn 1989).Furthermore, the various SPFs are not statistically consistent with each other, even though their overall shapes are qualitatively similar.In this context, it is intriguing that this SPF matches observed near-infrared SPFs of debris disks.
The second empirical SPF used in this study was measured directly from SPHERE images of HR 4796A (Milli et al. 2017).The HR 4796A SPF is distinct from the SPF trends seen in Hughes et al. (2018) and was used as a test case for the HD 114082 and HR 4796A datasets to demonstrate the differences in best-fit model appearances and noise-scaled residuals (hereafter called residuals) when using different SPFs.Although SPFs have been measured for HD 32297 (Duchêne et al. 2020), HD 35841 (Esposito et al. 2018), and HD 114082 (Engler et al. 2023), they are not considered for this analysis as their shapes were found to be largely similar to the trends observed  (Milli et al. 2017).A 9th-degree polynomial fit is used to account for the similar shapes observed in many dust environment SPFs.
in Hughes et al. (2018).SPFs for the other targets in this sample have not been measured.
Physical disk models are axisymmetric and generated in a threedimensional space and rotated according to constraints on position angle and inclination with respect to the observer to match the viewing geometry of each system.After the disk model is properly rotated, the three-dimensional disk model cube is collapsed along the lineof-sight direction, creating a two-dimensional image of the disk.

Forward Modeling Process
DiskFM (Mazoyer et al. 2020) was used to perform forward modeling on all disk models.For every GPI integration, a grid pattern imprinted on the pupil plane mask diffracts on-axis light from the target star to create a pattern of four fainter satellite spots in the image plane.Each disk model, as described in the previous section, is convolved with a PSF core generated from all satellite spots in all wavelength channels of every science frame where the disk does not overlap and the satellite spot has a SNR > 3. Finally, the KL basis vectors from the pyKLIP-ADI reduction of the original science dataset are projected onto the convolved model.This allows for a PSF subtraction to be applied to the forward model in an identical manner to the science data, without having to inject the model into an "empty" observational sequence and perform KLIP-ADI again to recalculate KL basis vectors, which would have different patterns of self-subtraction artifacts than the reduced science dataset.
To generate dataset-specific noise maps for calculating forward modeling likelihood, we rotate science frames in randomly-generated orientations to create a time-collapsed datacube that medians out a PSF halo pattern while preserving radial noise properties.This approach is similar to the approach described in Gerard & Marois (2016), although we choose random orientations as opposed to opposite orientations from North-up to decrease the likelihood that disk signal overlaps between individual science frames.PyKLIP-ADI is then utilized on the randomly-rotated science dataset to eliminate traces of the disk signal entirely and create a noise map that has been reduced with the same KLIP parameters as the non-rotated science dataset.From the randomly-rotated reduced image, the standard deviations of concentric rings 3 pixels in width are calculated.The final noise map consists of 3-pixel wide concentric rings with each ring containing the standard deviation of that same region in the "randomly-rotated" noise map.The noise map is then multiplied by a scalar of 3 for all model analyses, as described in Chen et al. (2020) and Mazoyer et al. (2020).We introduce this factor as we expect our noise maps to be underestimated, due to correlated noise features that our approach cannot reproduce.This method is also used to recover accurate error bars for planet photometry (e.g.Galicher et al. 2018).Mazoyer et al. (2020) showed that in almost all cases this factor was enough to accurately recover error bars for disk parameters, although they only tested this approach for a small sample of model systems.As shown in Table 5, this factor may not be the most ideal choice for all model analyses, in some cases leading to overestimations of the noise and therefore very low  2 red .This noise scaling factor is likely not consistent between different datasets, but we retain the factor of 3 to be consistent to previous high contrast imaging studies and so that all disk analyses are uniform.A synopsis of the steps and components of DiskFM is shown in the gold and purple boxes of Figure 1.

Likelihood Calculation
The forward model is compared to the data by measuring where  describes the region over which model likelihood was calculated.DiskFM utilizes the emcee package (Foreman-Mackey et al. 2013) for MCMC iterative analysis.Likelihood calculation is performed with an MCMC wrapper that maximizes  −  2 /2 in a masked region unique to each disk in the sample for both the disk forward model and science data.This calculation is performed until chain convergence.For all disk target analyses, 120 walkers were used, calculated for at least 5000 iterations.More iterations were added as necessary to confirm converged behavior within the MCMC chains.
After at least 300 iterations of converged behavior, all iterations up to that point were excluded as burn-in when generating the posterior distribution functions of our parameters.
Model likelihood was only calculated over specific regions where disk flux is most apparent and some background with no apparent disk signal is present.For systems where we tentatively resolve the back sides of disks-HD 35841, HD 106906, and HR 4796A-concentric rings enclosing the spine of the disk are created, with some space given between the apparent edge of the disk and the edge of the mask for the inclusion of the noise background in the likelihood calculation.If the front side of the disk passed in front of the FPM, the edge of the FPM was used for the inner boundary of the likelihood calculation mask.Regions interior to these concentric ellipses were further excluded if they contained significant amounts of noise.For close to edge-on systems HD 32297, HD 110058, and HD 111520, a flared bar shape was used as the likelihood calculation mask, with the inner radial boundary set slightly outside of the FPM.Interior to this inner boundary, noise is expected to be high due to the proximity to the host star and FPM.
The uncertainty in Equation 4 is taken from the generated noise map described in §5.2.Although spatially and spectrally-correlated noise among pixels is a concern for direct imaging IFS datasets (Greco & Brandt 2016), the scalar factor applied to the generated noise map as described in Mazoyer et al. (2020) was found to be sufficient in describing noise properties of an image.
Disk models have either 8 or 9 physical free parameters depending on disk inclination and apparent disk thickness.The 8 free parameters common to all disks are the inner cutoff radius  in , critical radius  C , surface dust density inner power law index   , surface dust density outer power law index  out , aspect ratio   , inclination , position angle  , stellocentric offset in the disk plane along the projected major axis  (where positive  corresponds to a stellocentric offset in au in the disk plane as defined in Table 4).A flux normalization scaling is used as an additional free parameter but is only used for scaling values in the model images to match the data; the value of flux normalization does not have a physical interpretation; this parameter is marginalized in the corner plots for our model analyses in Appendix A. For disks where we resolve the front and back sides of the disk, an additional free parameter  is parameterized (where positive  corresponds to a stellocentric offset in au in the disk plane pointing toward the observer along the minor axis), describing the stellocentric offset in the disk plane along the projected minor axis.Although we do not resolve the back side of HD 106906, we also parameterize  due to its evidence of asymmetric structure (Kalas et al. 2015;Crotts et al. 2021).The other targets in this sample do not contain enough spatial information to provide a meaningful constraint on .
MCMC initial parameters for  in ,  C , , and   are set as the constrained parameters of these disk targets as described in previous studies of scattered-light resolved imaging of disks, many of which were summarized in Esposito et al. (2020).Initial  in and  out parameters are chosen to represent median cases of power law indices, while the initial   parameters for HD 35841, HD 114082, HD 146897, and HR 4796A are derived to be in similar orders of magnitudes of previous studies (Esposito et al. 2018;Engler et al. 2023;Goebel et al. 2018;Milli et al. 2015;Olofsson et al. 2022).For HD 32297, HD 106906, HD 110058, and HD 111520 we opt to choose initial   parameters from the apparent vertical structure evident in our reduced GPI images.MCMC prior ranges are designed to be broad for radii, power law index, aspect ratio, and flux normalization parameters and narrow for inclination, position angle, and stellocentric offsets to match gross overall morphology.The vertical density distribution power law index  vert , the disk flaring index, and the outer radial cutoff  out are fixed.A summary of initial MCMC parameters and prior ranges are shown for each disk in Table 4.A summary of all DiskFM inputs (reduced science dataset, free parameters, fixed parameters, and the SPF) are shown in the blue boxes of Figure 1.The outputs of DiskFM, shown in green boxes, include the maximum likelihood best-fit model and posterior probability density distributions of free parameters.
Three disk systems -HD 106906, HD 110058, and HD 111520have significant asymmetrical structure in their scattered-light images (Kalas et al. 2015;Kasper et al. 2015;Draper et al. 2016;Crotts et al. 2024).The disk model utilized in this analysis does not have any parameterization that accounts for the strong disk asymmetries seen in these systems.Although a stellocentric offset can account for some brightness asymmetries, a reasonable offset is not sufficient to explain the significant brightness asymmetry observed in HD 111520.As seen in §6.1, strong residuals are present in the best model fits.As a result, we choose to conduct independent disk forward-modeling analyses on each side of these asymmetric disks (half mask models).Although we cannot fully constrain certain morphological properties because of the nature of this process, we seek to understand the value of applying the generic scattering phase function described in §5.1 to debris disk systems in general.

Uniform Disk Model Results
The best-fit models for each disk based on the uniform approach with the generic SPF, a continuous likelihood mask, and their residuals are shown in Figures 3-5.Median likelihood constrained parameters with 1 error bars and 3 upper and lower limits are reported in Table 5, except for HR 4796A 1 .Posterior distributions for all model analyses are reported in Appendix A. Figure 3 shows best-fit models of HD 32297, HD 35841, and HD 106906.Figure 4 shows the best-fit models of HD 110058, HD 111520, and HD 114082.Figure 5 shows the best-fit models of HD 146897 and HR 4796A.
Tight constraints (1 < 0.3 • ) on  were achieved for all systems except for HD 110058, for which we only find a lower limit, and HD 114082, where we find a double-peaked solution with 1 > 0.3 • .In both of these cases, the compact nature of the system may hinder the ability of DiskFM to identify a clear  solution.  was also tightly constrained (1 < 0.3 • ) for all systems.In general,  and   are typically the easiest parameters to constrain in a debris disk model, mostly due to clearly resolved positions.
Our constraints on the radial dust density distribution suggest two families of results: disks where the radial dust density distribution can be described by one power law ( C <  in or  C >  out ) and disks where the radial dust density distribution can be described by two power laws ( in <  C ,  C and  in and/or  out well-defined).Inner disk properties ( in and  in ) were difficult to constrain for most systems, likely due to the lack of line-of-sight resolution for these regions in close-to-edge-on systems.Outer disk properties ( out ) were comparatively easier to constrain, as a consequence of these regions being more clearly resolved in most systems.In the case of HR 4796A, the outer radial profile prefers a steep power law index outside of our prior range, suggesting a very sharp outer edge.From our posterior distributions, we conclude that the dust density distribution of disks around HD 35841 and HD 111520 can be modeled with two power laws, while the dust density distribution of disks around HD 106906, HD 110058, HD 114082, HD 146897, and HR 4796A can be modeled with a single power law. 1 We report median likelihood parameters separately for the generic SPF analysis of HR 4796A, as the poor model fit compared to the Milli et al. (2017) SPF analysis (see §6.3) does not provide meaningful constraints on disk morphology:  in = 74.17+0.08  −0.07 AU,  C = 70.14+0.22 −0.10 AU,  in > 1.34, 1 constraints were also found for the aspect ratios of all systems except for HD 146897 and HR 4796A, for which we achieve upper limits.The less-inclined and compact nature of HD 146897 as well as the bright PSF halo feature overlapping the front side of the disk can bias estimations of inclination and vertical extent.As HR 4796A is also relatively bright, nonlinearity in the KLIP reduction may also similarly bias estimates of some structural properties of the disk such as the vertical density distribution.
Narrow constraints (1 < 0.3 au) on  were identified for all systems except for HD 111520 for which we find a tightly constraining lower limit.Four systems appear to have significant offsets at or beyond the 3 confidence level-HD 106906, HD 110058, HD 111520, and HR 4796A-consistent with previous analyses of these systems.In scenarios where  was a parameter, 1 constraints were found for all systems (1 < 1.1 au) except for HD 106906, for which we find a lower limit.The lower limit suggests strong asymmetry, which has already been inferred from previous investigations of the system (Kalas et al. 2015;Lagrange et al. 2009;Crotts et al. 2021;Crotts et al. 2024).In our analyses, we identify well-fitting models that attempt to account for this asymmetry by introducing offsets along both the major () and minor () axes of the disk.However, residuals are still present, suggesting that stellocentric offsets alone may not be enough to account for such asymmetries.
We also calculate  2 red for all of our maximum likelihood models.While this value is typically used as a metric for quality of the model fit, we emphasize that this quantity is highly sensitive to the shape and size of the likelihood calculation mask and should be treated as more of a guideline.The likelihood calculation regions contain pixels with no disk signal in either model or data, naturally reducing  2 red by including additional degrees of freedom with little informational value.Additionally, this value is highly sensitive to the choice in multiplicative scaling factor applied to the noise map, with the ideal choice in scaling factor not likely to be consistent between different datasets.For some of our analyses, the noise is potentially overestimated, but reducing this factor runs the risk of underestimating error bars on our constrained parameters, a problem encountered in Mazoyer et al. (2020).We retain this scaling factor of 3 to be consistent with previous literature analyses and for uniformity in our analysis, in addition to providing conservative estimations of our constrained parameters.
Compared to other targets in this analysis, HD 35841 and HD 114082 display the smoothest and overall lowest residual maps, with the strongest residuals (≳ 2) appearing close to the FPM, likely due to deviations from Gaussian noise in this region.No strong residuals appear to overlap with disk signal.
For the asymmetric systems HD 106906, HD 110058, and HD 111520, structured, positive residuals (≳ 1) appear in regions that overlap with disk signal.In the case of HD 106906, DiskFM prefers solutions favors a best-fit model that places the center of the ring SE of the FPM.For HD 111520, the best-fit model appears to prefer a solution that attempts to leverage the brightnesses of both sides, strongly preferring solutions with offsets that push against the prior boundary of 5 au NW of the star.Finally, for HD 110058, the bestfit model itself appears to have an offset in   compared to the perceived location of the disk.This offset may be the best attempt of the model to match the observed "S"-shape seen in the reduced science image.The results for all three of these systems highlight the limitation of our model setup to constrain properties of highly asymmetric systems.
For HD 32297, strong positive residuals (Res.≳ 2) are seen in the region within 0. ′′ 2 of the inner edge of the likelihood calculation mask, centered at the midplane of the disk.This region overlaps with the highest SNR values in the reduced data, and are likely inducing nonlinearity in the KLIP-ADI reduction (Pueyo 2016), breaking down the linear expansion assumed in the forward modeling process.Regions where the absolute value of the residuals exceed 1 are also present far from the star and throughout most of the mask.
For HD 146897, significant positive and negative residuals are structured on the front and back sides of the disk respectively.The inability to achieve a low residual (|Res.|≲ 1.5) forward model in this case can likely be attributed to the overall noise properties present in regions near the front side, with residual PSF halo structure present that cannot be disentangled from disk signal.Because of the halo structure biasing the brightness of the front side of the disk, DiskFM cannot effectively identify models with a bright enough front side without making the back side too bright compared to the data using the generic SPF.
One target, HR 4796A, exhibits significant and azimuthallystructured residuals (|Res.|≳ 2) within its likelihood calculation regions.The distribution of residuals in this case appears to be bimodal rather than normally distributed, with high positive residuals present from the ansae to the back side of the disk and high negative residuals present on the front side of the disk.This behavior is expected, given our choice in using the generic SPF.In comparing the measured SPF of HR 4796A from Milli et al. (2017), the phase function is comparatively brighter at > 75 • than the generic SPF that we initially chose for the system.The generic scattering phase function is extremely forward-scattering, with less emission at the ansae of the disk and little back-scattering comparatively.Strong positive residuals (Res.≳ 2) are shown on the back side and ansae of the HR 4796A disk.

Half Mask Model Results
Four of the eight targets (HD 106906, HD 110058, HD 111520, and HR 4796A) have demonstrated asymmetric features in previous studies (Kalas et al. 2015;Kasper et al. 2015;Draper et al. 2016;Chen et al. 2020;Crotts et al. 2024).Median likelihood parameters related to asymmetry ( and ) for these four systems were also found to be greater than zero within 3.
As test cases, we perform additional analyses on three of these asymmetric disks (HD 106906, HD 110058, and HD 111520), treating each each side of the disks separately, as seen in Figures 6-8.These test cases assume that the disk is axisymmetric on either side of the star but are independent from each other.The large sizes of these systems allow for a sufficient region over which to calculate the likelihood, even when half of the disk is not being considered.For example, disk spine curvature is typically well-defined even when only considering half of the disk.For these analyses, the initial MCMC parameters and priors are consistent when modeling both sides of each disk, with the exception of the NW side of HD 110058, where initial investigations of the system preferred model solutions in the   outside of the initially chosen prior range.Therefore, for the NW side of HD 110058, we extend the upper bound of the prior range to 165 • .We do not expand the prior range for the SE side, as model solutions preferred a   well within the initially chosen prior range.The likelihood masks for these systems are split in half, allowing the likelihood calculation to ignore one side of the disk that may appear to have different morphological properties compared to the other side.We chose not to conduct a half-mask modeling analysis for HR 4796A despite its asymmetric nature, as the median likelihood value of  corresponds to an offset less than one pixel in length.
When determining best-fit models for each side separately, we are able to achieve lower residuals without systematic patterns of disk signal in the regions where likelihood was calculated.Regardless, radial and vertical dust density parameters are still found to be consistent with the initial analyses for HD 110058 and HD 111520 where the entire likelihood mask is used.These two systems do not appear to have morphological asymmetries as significant as HD 106906, with HD 110058 containing diffuse asymmetrical features around the ansae and HD 111520 containing a strong brightness asymmetry and tentative warp on the SE side (Draper et al. 2016).
In modeling the two sides of HD 106906 separately, we are able to constrain positional and outer disk properties and achieve upper/lower limits for stellocentric offsets.The most significant differences are evident in the median likelihood values of ,  ,   , and .While  may be somewhat degenerate with   , the distinct differences further highlight the asymmetric nature of the system, suggesting that the NW side is less inclined but vertically thicker than the SE side.Both of these models still appear to suggest large offsets, potentially suggesting that stellocentric offsets alone are still not sufficient in finding best-fit models for the disk even when analyzing both sides separately.To more robustly characterize the morphology of this system, a disk model that can create truly eccentric disks may be needed.
For HD 110058, 1 constraints for both sides are achieved for  out .The compact nature of the system is still the most limiting factor in modeling this system, as the lack of spatial resolution makes all parameters more difficult to constrain, particularly radii and  in .Otherwise, parameters for both sides appear consistent with each other within 3 except for  , likely due to the asymmetric warps on both sides of the disk (Stasevic et al. 2023;López et al. 2023).Appendix B3 contains a more detailed discussion of the properties of the HD 110058 disk and comparisons to previous literature.
In the two-sided model analysis of the HD 111520 disk, residuals no longer contain structure and are less (|Res.|≲ 1) compared to modeling the whole disk.1 constraints in  out ,   , and   were found for both sides.Constrained parameters for both analyses are consistent with each other except for  , likely due to the tentative warp observed on the SE side.
While constraints on some density properties of modeling either side of a disk are consistent with each other, namely  out and   , the greatest differences are highlighted in the constraints of   and  (and   in the case of HD 106906).The inconsistencies further highlight the asymmetric nature of these systems, as it suggests that a single model approach cannot unify both sides of a strongly asymmetric system; additional model mechanisms should be explored to treat such systems.

Different Scattering Phase Function Approach Results
Given that the previously measured SPF for HR 4796A is distinct from other systems, we performed an additional model analysis with the measured SPF from Milli et al. (2017).The initial parameters and priors are kept consistent from the previous analysis, and our best-fit model is shown in Figure 9.
The improvement in the best-fit model for HR 4796A is evident from the significantly lower residuals of the best-fit model, particularly at the back side and ansae of the disk (|Res.|≲ 2), in addition to the relatively lower  2 red (0.41 for the Milli et al. 2017 SPF com-pared to 1.90 for the generic SPF).While this result supports that the generic SPF is not applicable to all dust systems, it bolsters the methodology of this approach, namely that well-constrained morphological properties of disks can still be determined with a fixed SPF.Additionally, all median likelihood parameters and limits except for the aspect ratio were consistent between both modeling analyses, suggesting that using an inaccurate SPF may still be able to provide a first-order estimate of some morphological properties for a well-resolved and bright disk system.Our findings suggest much steeper radial profiles than our prior limits, which is consistent with the sharply defined edges observed in the disk.The vertical structure profile of the disk is also difficult to constrain, likely due to either a lack of resolution along the vertical direction and/or the removal of astrophysical signal from self-subtraction.To further assess the impact of SPF choice in constraining morphological parameters, we also conduct an analysis of the fainter HD 114082 ring with the HR 4796A SPF, even though the generic SPF analysis of this system provided a good model fit.The prior ranges, initial model parameters, and likelihood mask are kept consistent from the generic SPF analysis (except for the prior range and initial  , as discussed in this section).Given the high degree of backscattering in the HR 4796A SPF compared to our generic SPF, the faint and marginally-resolved back side of the disk, and the forward scattering peak of the disk being located behind the FPM, a reasonable model for this system would require "flipping" the orientation of the disk such that the apparent front side of the disk aligns with the back-scattering regime of the HR 4796A SPF.To account for this, we added 180 • to both the prior range and initial  , simulating placing the front side of the disk on the North side.While we achieve a reasonable maximum likelihood model (Figure 9), the fit is still poor compared to our initial analysis with the generic SPF and the front side of the disk on the South side, with a  2 red of 1.99 compared to 0.90 with the generic SPF.Furthermore, to achieve this model, we had to place the front side of the disk on the Northern side of the image, inconsistent with previous studies of the system (Crotts et al. 2024;Engler et al. 2023).Additionally, we find that most constrained parameters are consistent between the generic and HR 4796A SPF analyses, aside from ,  in , and   .
Both of the HD 114082 and HR 4796A analyses with the Milli et al. (2017) SPF suggest that choosing inaccurate scattering properties will lead to poorer model fits.We stress however that these are only two test cases, and we leave a deeper investigation of SPF choice for multiple systems to a future study.Given the significant differences between the HR 4796A SPF and our generic SPF, we cannot make a  3 but for two more targets in the sample.The structured and high residuals present in the HD 146897 residual map are likely due to diffuse PSF halo structure biasing the forward modeling process towards a lower inclination solution.Another possible explanation could be the choice in SPF; at a higher inclination, the model would produce too much forward scattering, and the curvature around the ansae would not be resolved.The highly structured residuals in HR 4796A are likely due to an incorrect choice of SPF, which cannot reproduce the distribution of forward and back scattering present in the data.conclusion on how much the model fits are affected by more subtle changes in the fixed SPF used; this investigation is also left as a future study.

K1-band Modeling Results
For the two systems with K1-band datasets meeting our SNR threshold, best-fit model results are shown in Figure 10.For HD 32297, we apply the generic scattering phase function as the proxy for grain properties.As we have already demonstrated that the generic scattering phase function does not provide a low residual model for HR 4796A, we apply the measured scattering phase function for HR 4796A from Milli et al. (2017).
Compared to the H-band reduced image of HD 32297 and bestfit model, the K1-band reduced image of HD 32297 and best-fit model appear less affected by the effects of self-subtraction from PSF subtraction and nonlinearity in forward modeling, although high residuals (|Res.|≳ 2) are still present.As the emission from the disk at K1-band is relatively less than emission at H-band, self-subtraction and nonlinearity in the forward modeling process is likely to be less severe.Additionally, similarities in relative errors among parameters with a lower disk SNR could suggest lower residual amplitudes.
The broader PSF at K1-band may also contribute to the smoothing of finer structural features in the disk, allowing the model to more easily match the data.Regardless, our K1-band model fit still has significant residuals within 0. ′′ 1 of the inner edge of our likelihood mask.
The HR 4796A best-fit K1-band model achieves a low residual (|Res.|≲ 1) model fit to the K1-band reduced science image.Some positive residuals persist on the NE side of the disk, particularly around the ansae and the front side of the disk.A brighter NE side has been observed in previous studies (e.g.Milli et al. 2017) suggesting that HR 4796A is an eccentric system, and while our disk model can produce modest eccentricities by inducing stellocentric offsets, we are unable to reproduce the observed asymmetry, suggesting that a dust density enhancement may be present that cannot be created by our model setup.We are able to tightly constrain most of our free parameters, with the exception of  in and the aspect ratio. in can only be constrained while  C >  in , and the lack of distinction between  in and  C suggests that the inner radial profile of the disk is sharply defined.Similar to the H-band result, our posteriors suggest steep radial density profiles.2017) SPF was used for HR 4796A.Residuals appear cleaner in both systems compared to their H-band images, potentially due to the fainter disk brightesses at K. Our best-fit model for HD 32297 still has some significant residuals close to the edge of the FPM, and nonlinearity in the KLIP reduction may still be evident to a lesser extent.

Limitations of the Forward-Modeling Approach
From our analysis, not all systems yield low residual model fits.The linear approximation in Equation 4of Pueyo (2016) breaks down in the presence of astrophysical sources that are brighter than or of a similar brightness to local speckles, particularly when using "aggressive" KLIP parameters.Mazoyer et al. (2020) found that this forward modeling nonlinearity can be worse for bright and extended disks.Nonlinearity in the KLIP reduction can hinder the ability of DiskFM to properly treat self-subtraction and PSF convolution and is likely the cause of our poor model fits for HD 32297, the brightest disk in our sample.To mitigate this issue for future characterization of bright debris disks, a different PSF post-processing technique (e.g.NMF; Ren et al. 2018) may be necessary.Our inability to achieve a low residual model for HD 32297 prevents us from making any definitive conclusions on its geometry and dust density distribution properties, and will be excluded from our discussion and comparisons of morphological properties between the other systems in our sample.As one of the primary goals of this study is to compare dust density properties among our sample, we present ensemble results in §7.3 and §7.4.For a detailed discussion of constrained dust density properties to previous studies on a system-by-system basis, we refer the reader to Appendix B.
Although estimations of  and   appear mostly consistent with past studies (see Appendix B) despite a variety of observations from different instruments and near-IR filters, inconsistencies become apparent in estimating dust density distribution properties, evident in the midplane density profiles of previous studies and this work.The uniqueness of approaches and associated limitations in analyzing disks also prevent direct comparisons of structural properties, further complicating ensemble studies of debris disks.
Finally, the manner in which disk morphological and dust density distribution properties have been explored has been relatively limited up until recent studies within the past few years.Many studies create model grids in order to determine the best-fitting parameters and do not explore full parameter spaces.The use of MCMC and multinest sampling analysis have only recently been utilized in disk image modeling (e.g.Duchêne et al. 2020;Esposito et al. 2018;Crotts et al. 2021;Engler et al. 2023;Chen et al. 2020;Olofsson et al. 2022), allowing a more thorough search of parameter spaces.
As we obtain more high-resolution images of debris disk systems, consistent methodologies (e.g.Crotts et al. 2024;Olofsson et al. 2022) are needed to ensure more direct comparisons between systems with unique properties such as host star spectral type and age.While a uniform modeling approach may limit the ability to obtain the lowest residual model fits, the general consistency we find in our analysis to previous studies suggests that loss of fit quality does not severely bias the results of dust morphological modeling and that true astrophysical trends can still be identified.

H-and K1-band Model Comparisons
Excluding the poor model fits for HD 32297, we find that the constrained parameters of HR 4796A between H-and K1-band are consistent to within 3, suggesting that the two wavelengths are probing similar scattering grain populations of dust.This suggests that our methodology is not too sensitive to the choice in dataset as long as the SPF remains largely unchanged between different wavelengths.
Interestingly, use of the measured H-band SPF for HR 4796A can still achieve a best-fit model at K1-band with reasonable residuals.This is supported by findings from Chen et al. (2020), where the extracted HR 4796A scattering phase functions from J-, H-, K1-, and K2-band observations appear consistent with each other within uncertainties.

Vertical and Radial Extents
Overall, we find that vertical aspect ratios are all < 0.14.Although most results are consistent with previous findings, we note that our determination of   appears to suggest extremely narrow profiles for HR 4796A (vertical FWHM ∼ 0.1 au at  C ), unlike what was found in Olofsson et al. (2022) (vertical FWHM ∼ 4 au at  C ).Additionally, the shape of the vertical density distribution may affect the aspect ratio, as Olofsson et al. (2022) kept  vert as a free parameter whereas we fix ours.Our estimated aspect ratios are also notably smaller than what was found in Crotts et al. (2024), but it should be noted that the method used to measure vertical FWHM and aspect ratio in Crotts et al. (2024) did not account for disk projection effects and confusion between radial and vertical extents.
Within this sample, we find no obvious correlations between aspect ratios and spectral type or approximate age.However, this may be an artifact of our small sample size that is limited to A-and F-type stars, with five of our eight systems likely sharing similar ages as members of the Scorpius-Centaurus OB association (de Zeeuw et al. 1999).With a larger sample size of more diverse spectral types, Crotts et al. (2024) does identify a correlation of increasing aspect ratio with increasing host star temperature in the case of radially compact disks, however.Olofsson et al. (2022) investigated the effects of gaseous components of disks on dust density morphologies and found that intermediate levels of gas mass most strongly affect the vertical density distribution at scattered light wavelengths.From their model simulations, increases in the gas mass correlate with a lower aspect ratio or "thinner" disks as gas drag boosts the efficiency of vertical settling.Of the eight targets investigated in our sample, two disks (HD 32297; Greaves et al. 2016, HD 110058;Hales et al. 2022) have measured CO gas detections and one disk (HD 146897;Lieman-Sifry et al. 2016) has a tentative CO gas detection.The largest modeled aspect ratios appear to be around the star HD 110058, where   was found to be greater than 10% in modeling both the entire disk and the SE side independently.This may seem contradictory to the idea that a gasbearing disk would likely be thinner, but the distance and compact nature of HD 110058 could also suggest a lack of resolution along the vertical direction and prevent the formation of a definitive conclusion on the presence of gas and disk vertical extent.Additionally, vertical stirring mechanisms (e.g. from dynamical interactions with substellar companions) may supersede the dampening of inclination from gas drag.
To assess correlations between radial and vertical structure, we calculated the median likelihood relative radial widths from our model posterior distributions, defined as the FWHM of the radial density profile following Equation 2, restricted by  in and  out on either side and divided by the radius of maximum dust density.In Figure 11, we plot   as a function of the relative radial width.In a study of ALMAobserved debris disks, Terrill et al. (2023) identified a tentative trend of increasing relative radial width and aspect ratio, with HD 110058 exhibiting the largest aspect ratio and third largest relative radial width, assuming a 90 • inclination.We do not immediately identify correlations between properties in either comparisons, although the sample size of models is relatively small.In the most narrow systems (HR 4796A and HD 114082), a few mechanisms could induce smaller radial widths, including the sculpting/shepherding of dust along the inner and outer edges of a ring (Boley et al. 2012), truncation from external perturbers (Nesvold et al. 2017), primordial eccentricities from instabilities (Kennedy 2020), and gas/dust interactions (Lyra & Kuchner 2013).As a general case, Chiang et al. (2009) and Rodigas et al. (2014) identified a relation between the mass of a single interior planetary sculptor and the normalized relative width of a ring, with higher mass planets linearly related to a more radially thickened exterior debris ring.
With the exception of HD 110058, most of our constrained disk aspect ratios are consistent with "natural" scale height from radiation pressure (Thébault 2009), although planetary companions can still induce vertical stirring.Dong et al. (2020) investigated model scenarios of massive planets influencing dust grains and found that disk aspect ratios can be enhanced most strongly (∼ 0.05) for multipleplanet systems where the planets lie interior to the disk.Multiple interior planets can also lead to more radially narrow structures depending on orbital and mass configurations.In a case study of HD 106906, a dynamical model analysis in Nesvold et al. (2017) found that increasing the mutual inclination of an external perturbing companion increases the vertical extent of a disk as well.In the case of HD 110058, interior planetary companions could be responsible for the warps seen in scattered-light images, similar to the  Pic system (Heap et al. 2000;Lagrange et al. 2009).The detection of CO gas in HD 110058 (Hales et al. 2022) could also contribute to the radially narrow and warped geometry observed.López et al. (2023) also found that both interior and exterior planetary companions can produce the geometry observed in HD 110058 through secular perturbations, while Stasevic et al. (2023) concluded that only an interior companion could produce the observed warped geometry.Interestingly, the HD 106906 median likelihood whole disk and NW side radial widths are fairly broad (∼50 au), despite its eccentric shape likely being caused by an external perturber (Nesvold et al. 2017).Broad radial widths and heightened aspect ratios both would support dynamical models of planet scattering (Nesvorný 2015), although our aspect ratios for HD 106906 are not large relative to the rest of our sample.For more vertically-thin disks, self-stirring and/or secular interactions may be dominant (Matrà et al. 2019); this behavior may be relevant for systems such as HD 35841 and HD 146897.

Disk Eccentricity
Disk eccentricity in our model can only be parameterized in the form of stellocentric offsets.In general, the most significant offsets we observe are seen in the best-fit models of HD 106906, HD 110058, HD 111520, and HR 4796A, consistent with identified eccentricities and high degrees of asymmetry measured in polarimetric images in Crotts et al. (2024) 2 .Pericenter glow (Wyatt et al. 1999) is postulated as a cause of the asymmetry in HR 4796A, but despite inducing a pericenter glow from an offset, a positive residual remains, suggesting that a localized density enhancement (Schneider et al. 2009;Milli et al. 2019) may be needed to account for this effect.
In HD 106906, studies of dust dynamics suggest that a planetary perturber, either interior or exterior to the debris disk (e.g.HD 106906 b, Bailey et al. 2014) may be the cause of the asymmetry.Additionally, catastrophic collisions from large solid bodies could also explain the needle-like morphology observed in HD 106906 (see Jones et al. 2023 and references therein).Studies of HD 110058 have been limited in scope, but the large aspect ratio found in Crotts et al. (2024) and this study along with the characterized warps in Stasevic et al. (2023) and López et al. (2023) provide tentative evidence of a vertical stirring mechanism such as a companion.Dedicated dynamical sculpting studies of HD 111520 have not been conducted, and the extremely edge-on inclination further limits our ability to assess the true eccentricity of the disk.Giant impact modeling conducted in Jones et al. (2023) could also not reproduce the observed morphology of HD 111520, and suggest that an interior planet is more likely to produce the geometry observed.
In Figure 12, we compare median likelihood eccentricities of all disk modeling analyses to relative radial widths.A tentative positive trend is seen in eccentricity versus relative radial width, and we observe that all asymmetric systems and half mask disk models appear to have higher relative radial widths in general.We also observe a large range of eccentricities over a small range of relative radial widths, similar to the large range in   we observe.We can compare our results to findings from Kennedy (2020), who observed that the highly eccentric debris disks around Fomalhaut and HD 202628 appear narrower than expected from secular perturbation models and assumptions of zero initial eccentricity.The dashed line in Figure 9 of Kennedy (2020) represents a predicted positive-slope relationship of eccentricity and relative radial width in a zero initial eccentricity secular perturbation condition.Although this relationship is positive along with the tentative correlation we observe, the slope of our correlation is shallower.Our most eccentric disk models, associated with the whole and half mask models of HD 106906, appear to disagree with the assessment from Kennedy (2020) that the most eccentric systems have the most radially narrow sizes, although both studies contain small samples and explore vastly different grain size populations.Finally, similar to our comparisons of radial and vertical structure, the range in eccentricities could also suggest multiple stirring mechanisms that may be responsible for eccentric shape.

Implications for Grain Properties and Scattering Phase Functions
Efforts in modeling the SPF for debris disk systems most commonly utilize the HG phase function (Henyey & Greenstein 1941) as an analytical (but not physical) model to all scattering properties or fully parameterize dust grain properties assuming a uniform shape to dust grains such as the Mie model (Mie 1908).Both of these approaches, however, often fail to properly account for all observed scattering behavior and are inconsistent across multiple wavelengths and/or observing modes.In modeling analysis of HST-NICMOS observations of HD 181327, Schneider et al. (2006) identified inconsistencies between a minimum grain size ( min ) determined from observed mean asymmetry factor in the SPF and an  min determined from disk color and thermal SED.A more physically-grounded basis from which to model dust grains may be in the form of complex aggregates of spherical and/or fractal shapes, identified in some studies of solar system dust (e.g.Bentley et al. 2016).Although this treatment appears more realistic, scattered-light modeling of aggregate dust grains is a computationally daunting task (Tazaki et al. 2016;Arnold et al. 2019).The Distribution of Hollow Spheres (DHS; Min et al. 2016) approach attempts to circumvent the issue of computational efficiency of aggregates, but is also not always sufficient in reproducing features seen in both SPFs and polarized intensity surface brightness profiles (Arriaga et al. 2020).Arnold et al. (2022) has also assessed model SPFs for Mie, DHS, and aggregate grains, demonstrating similarities between them at low scattering angles, but deviations at larger scattering angles and across multiple wavelengths.The wavelength-dependent nature observed in Arnold et al. Figure 12.Eccentricities estimated using the posterior distribution of stellocentric offsets and the radii of maximum density as a function of relative radial widths.There appears to be a tentative positive trend between eccentricity and relative radial width.The slope of this trend, however, appears shallower than the slope of the dashed line in Figure 9 of Kennedy (2020), shown in red.This line represents the expected radial width in a zero initial eccentricity secular perturbation scenario.The assessment from Kennedy (2020) suggests that our systems should be roughly circular, even though that is not the case for at least HD 106906.
(2022) seems to contrast with many empirical SPF measurements, which show similar behavior between different wavelength regimes (e.g., HR 4796A, see Chen et al. 2020), and multi-wavelength measurements of SPFs are likely also needed in order to better constrain the grain properties of debris disk systems.
The limited number of debris disk SPF measurements (e.g.Graham et al. 2007;Duchêne et al. 2020) have been shown to follow trends in shape seen in solar system dust environments (Hughes et al. 2018).This generic SPF shape has been reproduced with highly porous grains where the minimum grain size is greater than 1 assuming both Mie theory (Grynko et al. 2004) and clusters of complex aggregates (Bentley et al. 2016) suggesting that debris disk dust shares many similarities to cometary dust and that overall dust composition is not as important of a factor.Our generic SPF uti-lizes measurements from various comets and planetary rings, but other approaches in generating a generic SPF from solar system dust have also been demonstrated (e.g., Schleicher & Bair 2011;Marcus 2007a,b).Additionally, laboratory measurements and simulated dust properties (e.g.Lolachi et al. 2023) could also inform the general shape of a generic SPF, and future studies of debris disk SPFs could incorporate these results in their analysis.
The achievement of low residual model fits for most of the debris disk systems is a result that supports findings and conclusions from Hughes et al. ( 2018) and previous studies that a common SPF could exist between different dust systems with only modest variations.We have demonstrated that we can achieve low residual debris disk models of the HD 35841, HD 106906, HD 110058, HD 111520, and HD 114082 systems, utilizing the same SPF informed from solar system dust.Although we did not attain a low-residual best-fit model of HD 32297, the measurement of the SPF in Duchêne et al. ( 2020) is still consistent with our generic SPF.These commonalities also support predictions that highly porous, aggregate grains may be a major contributing factor to the shape of an SPF.The robustness of a majority of our results demonstrates our methods as a uniform approach to constraining disk geometry without having to parameterize grain properties or utilizing physically limiting assumptions (e.g.ellipse fitting).
Despite these findings, we have identified one system where we could not achieve low residual model fits using the generic SPF.HR 4796A is known to harbor a distinct SPF shape as measured in Milli et al. (2017), with relatively higher scattering at the ansae and back side of the disk than what the generic SPF allows.Modeling of the HR 4796A scattering phase function still suggests that porous aggregates are likely the types of grains inducing scattering, but at a larger minimum grain size (> 5m) than predicted for other dust systems that follow the generic SPF (Milli et al. 2017;Chen et al. 2020).Although our approach with the generic SPF did not achieve low residual models to HR 4796A, the analysis was still informative in supporting that not all debris disk systems have SPFs similar to the generic SPF, and that a low residual model solution cannot be achieved if the scattering properties are treated incorrectly.Additional model tests using the HR 4796A SPF for HD 114082 further supports that SPF choice is an important consideration, as the bestfit model solution preferred model solutions where the front side of the disk was on the Northern side, notably inconsistent with other studies of these systems and our own disk analysis using our generic SPF.This conclusion, however, is founded upon using two very different SPFs, and does not account for more subtle variations, which we leave for a future study.Although we could not achieve a low residual model fit with the generic SPF, use of the measured SPF from Milli et al. (2017) achieved significantly lower residual model fits to HR 4796A, further supporting the utility of this approach in constraining morphological properties with fixed scattering properties.To further assess the sensitivity of our model fits to the choice in SPF, we demonstrated that using the HR 4796A SPF for a system fit well with the generic SPF led to poorer model fits and inconsistent morphological properties.

SUMMARY
With advances in direct imaging instrumentation, more scatteredlight debris disks have been resolved at small separations than ever before.With such a large sample of imaged debris disks, consistent and uniform characterization approaches should be implemented in order to make more direct comparisons between systems.We have performed a uniform forward-modeling analysis of eight bright debris disks imaged with the Gemini Planet Imager in total intensity light.Scattering properties were determined by empirically-informed scattering phase functions and not from parameterizing Henyey-Greenstein functions or specific grain properties.From our results, we were able to identify two families of debris disks: one where the midplane density profile can be described by one power law (disks with extremely sharp inner edges) and the other where the midplane density profile can be defined by two power laws (smoother declines of dust interior and exterior to peak density of the ring).Even though we find consistent results among many prior studies of these systems, inconsistencies shed light on the sensitivity of results to approaches of data reduction and modeling.
In applying the same empirically-informed scattering phase functions used in H-band for K1-band images of the same system, we still achieve consistent modeling results.The range of aspect ratios we find in our modeling analyses are mostly consistent with a "natural" disk scale height associated with radiation pressure, with the exception of HR 4796A and HD 110058.The thin aspect ratio of HR 4796A is inconsistent with prior results and may be related to issues in PSF subtraction and forward modeling.The broader aspect ratio of HD 110058 along with its narrow radial width and structural warps may be indicative of complex stirring mechanisms including perturbing planetary companions.We also identify asymmetric disk systems (HD 106906, HD 110058, HD 111520, HR 4796A) but cannot fully characterize their asymmetry due to the physically-limiting assumptions in our models.The range of eccentricities compared to relative radial width further highlights the diversity of systems and potential stirring mechanisms that cause divisions in relative width and aspect ratio.
We have demonstrated that rigorous morphological modeling can be conducted with a uniform and consistent set of assumptions about grain properties.While it is unlikely that most disks have the same dust populations, the scale-invariant nature of highly porous aggregates may cause SPFs of different dust environments to appear similar, regardless of inherent grain properties such as grain size distribution or composition.The SPFs from aggregate grain models also appear consistent with the generic SPF we determined from Solar System dust measurements.The achievment of well-fitting models in most cases suggest links between porous aggregates, Solar System cometary dust, and extrasolar debris disks.In systems where low residual model fits were not readily found (e.g.HR 4796A), it is likely that the scattering phase function applied was not representative of the system, or another factor (e.g., KLIP nonlinearity for HD 32297) may introduce biases into the modeling approach.
To further improve the robustness of this model, more observations at higher angular resolution are needed.Additionally, other PSF subtraction algorithms should be explored to avoid complications of KLIP nonlinearity and self-subtraction in bright disks.Future observations with next-generation ELTs (e.g.TMT, GMT, E-ELT) may prove critical in resolving finer structures that may be present within these systems at greater SNR.Additionally, further studies in dust grain dynamics may improve the robustness of our treatment of the dust density distribution.Radiative-transfer models in turn will also require more complexity accounting for asymmetric features such as eccentricity and warps.Finally, more empirical measurements of debris disk scattering phase functions are needed in order to support and identify common trends and properties of grains in a variety of different environments.

DATA AVAILABILITY
The raw data files for all observational sequences are available through the Gemini archive (https://archive.gemini.edu/searchform).The reduced data files, modeling data files, and other materials underlying this article will be shared on reasonable request to the corresponding author.The generic SPF used for this model analysis and the updated Python scripts for DiskFM used in this analysis will be available through Hom et al. (2024) and Github 3 .Inquiries regarding the scattering phase function from Milli et al. (2017) should be directed to Julien Milli.GPI H band ADI (KL#: 3) 5,074 iterations (+ 600 burn-in) with 120 walkers: 680,880 models  Esposito et al. (2018).Our analysis prefers a narrow ring with a dust density profile modeled with two power laws, while the analysis from Esposito et al. (2018) appears to prefer a somewhat broader density profile modeled with a single power law.The discrepancies between our results and Esposito et al. (2018) may be due to the differences in modeling approaches applied.Esposito et al. (2018) fit jointly to polarized and total intensity data, where the constraints placed on grain properties from the polarized intensity data can induce distinct brightness distributions between the two images.This in turn influences the morphological parameters of the model fitting that are not present in fitting total intensity data alone.

B2 HD 106906
Table B2 presents a comparison of our constrained parameters for all modeling analyses to previous studies.In our modeling analysis of HD 106906, we find shallow radial density profiles in all modeling scenarios, except for the SE side where  in is unconstrained.Figure B2 compares the midplane density profiles of 100 randomlyselected models from our posterior distributions to midplane density profiles from previous studies.Our whole disk model profile is most consistent with modeling analysis from Crotts et al. (2021), highlighting a broad ring with shallow inner and outer slopes.Overall,   et al. (2022), and this work.This could be due to Lagrange et al. (2016) fixing the inner power law index to a value of 10, forcing a narrower radial profile.Offsets are found to be significant enough to warrant our asymmetric disk modeling approach, and even larger offsets were also found in Crotts et al. (2021).In our analysis of the separate sides of HD 106906, morphological properties appear distinct from each other, with the SE model appearing to favor a larger radial extent and a much narrower aspect ratio, consistent with a needle-like asymmetry.This appears in contrast to observations of the outer halo of HD 106906 seen in HST images of the system (Kalas et al. 2015), although this may be unsurprising if a perturbing companion is influencing outer halo dust grains.
Our best-fit model for HD 106906 does not reproduce the brightness asymmetry observed in images despite containing large stellocentric offsets, likely due to the limited complexity of our model.The strong brightness asymmetry suggests that the disk could have high eccentricity and/or localized density enhancements/dearth.The best-fit models identified from analyzing the separate sides of the disk independently appear to fit the data more cleanly, but are inconsistent in  and  , potentially suggesting the presence of a warp and further emphasizing that the complex geometry of the system is not easily determined from a relatively simple disk model.

B3 HD 110058
HD 110058's resolved debris disk was first presented in Kasper et al. (2015).The most interesting features present in both SPHERE and GPI total intensity images of the system are symmetric warps present at about 0. ′′ 3 from the star oriented in opposite directions, creating an "S"-like shape.Similar to measurements in Crotts et al. (2024), our modeling analyses suggest vertically-broad profiles, which could be indicative of vertical stirring mechanisms such as perturbing planetary companions (see §7.3).Our estimates of modeling the whole disk achieve a similar result for  in as  C in Kasper et al. (2015) and the measurement of   from Crotts et al. (2024).In modeling the two sides of the disk separately, differences in   of 2 • between both analyses provide further evidence of the warps seen in images.The midplane density profiles from our models tend to prefer a narrow ring with a somewhat shallow outer slope.Regardless, our modeling results are tentative at best due to the loose constraints of the posterior distributions in all modeling scenarios of this system.The compact size of the disk and the observed "S"-shape are expected to be the most significant limiting factors in our interpretation of the results, as our model has less information to calculate a likelihood function over and contains no formalism for producing warps.Stasevic et al. (2023) provides quantitative constraints on the strength of the warps and suggest an inner perturbing planetary companion could be responsible, similar to the  Pictoris system.López et al. (2023) determined that the strengths of the warps on each side of the disk are distinct, and that constraints on disk measurements are also difficult to achieve due to the compact nature of the system.

B4 HD 111520
Our best-fit symmetric model solution for HD 111520 suggests a significant stellocentric offset, though it does not fully reproduce the strong brightness asymmetry present in the GPI images (Draper et al. 2016).Thus, a truly eccentric geometry is likely needed to properly describe the disk.In modeling the two sides of the disk separately, an offset along the major axis could not be constrained in the SE model.The NW model constrained a nonzero offset, but with a wide posterior that overlaps with zero.Crotts et al. (2022) did not identify a stellocentric offset from fitting the location of the disk spine, but a surface density asymmetry could produce the strong brightness asymmetry present.No previous studies have attempted to model radial dust density distribution properties, but we do find disk aspect ratios that are similar to other systems in our analysis, and radii roughly consistent with measurements from Crotts et al. (2024).Midplane density profiles from our posterior distributions support a narrow ring structure with shallow outer slopes.In their analysis of the GPI H-spec total intensity and H-pol polarized intensity images, Draper et al. (2016) identified a nearly edge-on geometry from modeling the vertical offset profile but did not incorporate disk models.Crotts et al. (2022) expanded this analysis by also investigating the GPI J-and K1-band total and polarized intensity images, applying an MCMC analysis with a simple inclined ring model.We find that our  estimates for both our whole disk model and asymmetric model approaches are consistent with the estimations from Crotts et al. (2022) to 3, while the   appears only inconsistent in analyzing the SE side of the disk.A warp does appear present along the SE side, also seen in polarized intensity images of the system in Crotts et al. ( 2022), but the reduced SNR on the SE side makes the existence of such a warp tentative.

B5 HD 114082
We achieve well-constrained parameters and relatively low residual model fits for HD 114082.The median likelihood parameters suggest a compact, narrow-ring structure as observed in previous studies (e.g.Wahhaj et al. 2016;Engler et al. 2023), although the posterior distribution functions show two families of distinct profiles, one with a smaller peak radius and truncated slope (single power law-preferred) and one with a slightly higher peak radius and shallower inner slope (two power law-preferred).Figure B3 reveals that our model analy- ses preferred narrow-ring structures similar to the analysis of Engler et al. (2023), albeit with a shallower outer profile.HD 114082 is also one of a handful of debris disks with a measured SPF; Engler et al. (2023) determined the measured SPF of HD 114082 to be similar in shape to other solar system dust environments and other debris disks, including zodiacal light and the Saturn D68 ring which our generic SPF is partially derived from.The aspect ratio appears much broader than what was found in Engler et al. (2023), with a vertical FWHM approximately twice as large.This discrepancy is likely related to the quality of data and differing model approaches.In the SPHERE IFS and IRDIS images of the system, the front and back sides of the ring are resolved, while GPI was only able to resolve the front side of the disk.The compact nature and inclination of the disk may also lead to degeneracies between radial and vertical density properties.

B6 HD 146897
Table B3 presents comparisons of our median likelihood model parameters to previous studies of the system.In comparing the midplane density profiles from our analysis to previous studies (Thalmann et al. 2013;Engler et al. 2017;Goebel et al. 2018) as in Figure B4, we identify a broad range of profiles, suggesting that this system presents large uncertainties in general.Despite this broad range, there does appear to be a steep inner slope around 50 au seen in Goebel et al. (2018) and this work.Thalmann et al. (2013) also has a steep inner slope, although their inner radial density properties were fixed.The analysis from Engler et al. (2017) contrasts from these other profiles, with a broader midplane density profile in general peaking at a larger radius.The compact nature of this system, overall noise levels near the FPM, and inability to resolve the back side of the disk despite the moderate inclination of the system can contribute to confusion among radial and vertical density properties.One explanation for the high residuals could be bright PSF halo features overlapping with the SW side of the disk biasing the likelihood calculation of the MCMC, resulting in a disk model that is unable to match the brightness of the front side of the disk without overestimating the brightness of the back side of the disk.Additionally, HD 146897 could have a distinct SPF similar to HR 4796A, inducing systematic residual structure between the front and back sides of the disk, although the compact nature of the system prevents a definitive conclusion on this behavior.

B7 HR 4796A
Table B4 presents comparisons of our median likelihood model parameters to previous studies of HR 4796A.Our best-fit model using the SPF from Milli et al. (2017) fits the data reasonably well at both H-and K1-band, and we find consistent results among our estimations of radii and offsets to prior studies within 3 (e.g.Milli et al. 2015;Olofsson et al. 2020Olofsson et al. , 2022;;Crotts et al. 2024).The posterior distributions of the radial power law indices suggest steep profiles, which is unsurprising given the sharply defined edges of the disk.Midplane density profiles from our analysis appear consistent to narrow ring width determined in previous studies, as seen in Figure Overall, our midplane density profiles appear consistent to all three studies, suggesting that the system is a narrow ring, although our models prefer a steeper inner cutoff.
B5.The most significant inconsistency is seen in our determination of the aspect ratio, with our posterior distributions suggesting an unexpectedly thin ring compared to analyses by Milli et al. (2017) and Olofsson et al. (2022).Unusually thin vertical structure was also identified in Chen et al. (2020), which analyzed the same datasets we have used.This could suggest systematic differences between the GPI and SPHERE analyses, such as the treatment of the instrumental PSF.
Another consideration is the apparent brightness asymmetry between the NE and SW sides of the front-side of the disk.By introducing a stellocentric offset to the ring, pericenter glow (Wyatt et al. 1999) may be induced.Our best-fit model appears unable to compensate for this brightness asymmetry, even after considering stellocentric offsets along both the major and minor axes of the disk.This result may support findings in Olofsson et al. (2019), Milli et al. (2019), and Chen et al. (2020) that a dust density enhancement at pericenter may be responsible for the asymmetry.

Figure 1 .
Figure1.A diagram demonstrating the forward modeling process using DiskFM.A forward model is generated from a fixed scattering phase function, a set of fixed parameters, and a varying set of free parameters (blue boxes).From the reduced science images, a model PSF and representative noise map are generated.This model PSF is then convolved with all disk models, with the noise map used for likelihood calculation.The reduced science images are post-processed with pyKLIP-ADI, and the KL-basis vectors are saved for projection onto all forward models generated with DiskFM.The KL-basis vectors, model PSF, and representative noise map (gold boxes) are all used for the forward modeling of a disk given a set of inputs.A  2 function is used to inform the MCMC sampler of forward model minimization.The outputs of DiskFM include the maximum likelihood best-fit disk model and posterior probability density distributions of free parameters (green boxes).Purple boxes represent calculation procedures internal to DiskFM.

Figure 2 .
Figure2.Our generic SPF overplotted with measured SPFs from solar system dust and the HR 4796A disk(Milli et al. 2017).A 9th-degree polynomial fit is used to account for the similar shapes observed in many dust environment SPFs.

Figure 3 .
Figure 3. Modeling results from the initial uniform approach process for three targets in the sample, with the white outline representing the shape of the likelihood mask.The KLIP Reduced Data and DiskFM best-fit model are normalized to the maximum signal within the likelihood mask region of the KLIP Reduced Data.The noise-scaled residuals are shown in the right column.The strong residuals in the HD 32297 best-fit model are likely due to the inherent brightness of the disk, creating strong self-subtraction wings in the KLIP-ADI reduction and inducing nonlinearity in the pyKLIP forward modeling process.In the HD 35841 best-fit model, the strongest residuals (darkest red and blue regions in the residual map) are associated with regions closest to the FPM, where noise is expected to be high.In the HD 106906 model, the residual map appears to have structure along the spine of the disk, particularly on the SE side.

Figure 4 .
Figure 4. Same as Figure 3 but for three more targets in the sample.The HD 110058 residual map shows a pattern along the spine of the disk, suggesting mismatches in   for either side.The HD 111520 residual map highlights the strong brightness asymmetry first identified in Draper et al. (2016), where the best-fit model cannot reconcile the brightness between the NW and SE sides.The possible presence of a warp on the SE side also causes a   discrepancy between the data and the best-fit model.The HD 114082 best-fit model appears to have high residual features that do not overlap with the disk region and are likely associated with higher levels of noise in the immediate vicinity of the FPM.

Figure 5 .
Figure 5. Same as Figure3but for two more targets in the sample.The structured and high residuals present in the HD 146897 residual map are likely due to diffuse PSF halo structure biasing the forward modeling process towards a lower inclination solution.Another possible explanation could be the choice in SPF; at a higher inclination, the model would produce too much forward scattering, and the curvature around the ansae would not be resolved.The highly structured residuals in HR 4796A are likely due to an incorrect choice of SPF, which cannot reproduce the distribution of forward and back scattering present in the data.

Figure 6 .Figure 7 .Figure 8 .Figure 9 .Figure 10 .
Figure 6.Best-fit modeling results for separate analyses of the NW and SE sides of HD 106906.The models are distinctly different, most notably seen in the  , , and vertical structure profiles of both best-fit disk models.

Figure 11 .
Figure 11.Aspect ratio as a function of radial disk FWHM divided by the radius of maximum dust density.No trends are observed as a function of   .

Figure B1 .
Figure B1.Normalized midplane dust density profiles from 100 randomlyselected models from our posterior distribution (gray) of HD 35841 models, the maximum likelihood of our models (red), and the maximum likelihood model fromEsposito et al. (2018) (blue).The range along the x-axis is restricted to within our likelihood calculation mask along the disk major axis for consistent comparison.Overall, our analysis prefers a narrow ring, twopower law solution for midplane density, where the analysis from Esposito et al. (2018) prefers a single power-law solution and a relatively broader outer profile.Notably, our models prefer solutions that place dust interior to the inner edge of the profile fromEsposito et al. (2018).

Figure B2 .
Figure B2.Same as Figure B1 but for 100 randomly-selected models from our posterior distribution (gray) of whole disk HD 106906 models, the maximum likelihood of our models (red), and the best fit models from Lagrange et al. (2016) (blue), Crotts et al. (2021) (orange), and Olofsson et al. (2022) (magenta).Overall, our analysis prefers a broad-ring, single power-law, flatter inner profile, and shallow outer slope profile for midplane density, similar toCrotts et al. (2021).

Figure B3 .
Figure B3.Same as Figure B1 but for 100 randomly-selected models from our posterior distribution (gray) of HD 114082 models, the maximum likelihood of our models (red), and the best fit model from Engler et al. (2023) (blue).Both the analysis from Engler et al. (2023) and our own work suggest steep inner truncations at ∼30 au, although Engler et al. (2023) appears to prefer a shallower outer slope for midplane dust density.

Figure B4 .
Figure B4.Same as Figure B1 but for 100 randomly-selected models from our posterior distribution (gray) of HD 146897 models, the maximum likelihood of our models (red), and the best fit models from Thalmann et al. (2013) (blue), Engler et al. (2017) (orange), and Goebel et al. (2018) (magenta).Our midplane density profiles appear most similar to Goebel et al. (2018) albeit with a steeper outer profile, while the slope of the outer density profile is consistent with Engler et al. (2017).

Figure B5 .
Figure B5.Same as Figure B1 but for 100 randomly-selected models from our posterior distribution (gray) of HR 4796A -band models, the maximum likelihood of our models (red), and the best fit models from Milli et al. (2015) (blue), Olofsson et al. (2020) (orange), and Olofsson et al. (2022) (magenta).Overall, our midplane density profiles appear consistent to all three studies, suggesting that the system is a narrow ring, although our models prefer a steeper inner cutoff.

Figure B8 .
Figure B8.Posterior distribution for HD 106906 whole disk H-band models.

Figure B9 .
Figure B9.Posterior distribution for HD 106906 NW H-band models.

Figure B10 .
Figure B10.Posterior distribution for HD 106906 SE H-band models.

Figure B11 .
Figure B11.Posterior distribution for HD 110058 whole disk H-band models.

Figure B12 .
Figure B12.Posterior distribution for HD 110058 NW H-band models.

Figure B13 .
Figure B13.Posterior distribution for HD 110058 SE H-band models.

Figure B14 .
Figure B14.Posterior distribution for HD 111520 whole disk H-band models.

Figure B15 .
Figure B15.Posterior distribution for HD 111520 NW H-band models.

Figure B16 .
Figure B16.Posterior distribution for HD 111520 SE H-band models.

Table 1 .
Target List summary.H and K magnitudes originate from 2MASS photometry

Table 2 .
Summary of observations for the sample.Seeing estimates are not available for HD 114082 and HD 146897 observations, as the DIMM seeing monitor at Gemini-South was inoperable at the time of observation.

Table 3 .
KLIP reduction parameters for all targets, chosen to maximize disk SNR.

Table 5 .
MCMC median likelihood parameters with ±1 error bars and 3 upper/lower limits.
These upper/lower limits span most of the prior range but have low probability tails.

Table B1 .
Comparison of our results to the previous sub-arcsecond resolution study of HD 35841 from Esposito et al. (2018).

Table B2 .
Comparison of our results to previous sub-arcsecond resolution studies of HD 106906.References: 1. Lagrange et al. (2016), 2. Crotts et al. (2021), 3. Olofsson et al. (2022)  These values were not determined as part of the dust density distribution model and were found via other techniques.

Table B3 .
Comparison of our results to previous sub-arcsecond resolution studies of HD 146897.References: 1. Thalmann et al. (2013), 2. Engler et al. (2017), 3. Goebel et al. (2018)  These values were not determined as part of the dust density distribution model and were found via other techniques.