A uniform analysis of debris discs with the Gemini Planet Imager II: constraints on dust density distribution using empirically informed scattering phase functions

ABSTRACT

Spatially resolved images of debris discs are necessary to determine disc morphological properties and the scattering phase function (SPF) thatantifies the brightness of scattered light as a function of phase angle. Current high-contrast imaging instruments have successfully resolved several dozens of debris discs around other stars, but few studies have investigated trends in the scattered-light, resolved population of debris discs in a uniform and consistent manner. We have combined Karhunen-Loeve Image Projection (KLIP) with radiative-transfer disc forward modelling in order to obtain the highest-quality image reductions and constrain disc morphological properties of eight debris discs 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 discs that may have resulted from substellar companion perturbations, along with a tentative positive trend in disc eccentricity with relative disc width. We also find that using a common SPF can achieve reasonable model fits for discs that are axisymmetric and asymmetric when fitting models to each side of the disc independently, suggesting that scattering behaviour from debris discs may be similar to Solar system dust.

scattering, techniques: high angular resolution, stars: circumstellar matter, infrared: planetary systems

1 INTRODUCTION

Since the first indirect detection of a circumstellar debris disc 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 discs have been observed in sizes ranging from tens to hundreds of au in diameter (e.g. Schneider et al. 1999). Due to the short time-scale processes of Poynting–Robertson drag and radiation pressure, debris discs 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 disc system. While spectral energy distributions (SEDs) can place some constraints on the radial extent, dust mass, and composition of system dust, debris disc images can more directly constrain overall properties of disc 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 disc dynamics and evolution (e.g. Lee & Chiang 2016), resolved debris disc images can also constrain the gravitational interactions between planets and discs (e.g. Liou & Zook 1999; Kuchner & Holman 2003; Quillen & Faber 2006; Wyatt 2006).

Direct imaging studies at the scales of 0|${_{.}^{\prime\prime}}$|01–1 arcsec have been explored at optical and near-infrared (IR) wavelengths with instruments such as the Space Telscope Imaging Spectrograph (e.g. Schneider et al. 2014, 2016) and ground-based adaptive optics (AO; 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 disc brightness from dust scattering is typically ∼10⁶ times fainter than the host star brightness. Current generation direct imaging instruments such as the Gemini Planet Imager (GPI; Macintosh et al. 2014) and the Spectro-Polarimetric High-contrast Exoplanet REsearch (SPHERE; Beuzit et al. 2019) utilize 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 discs for the first time (e.g. Currie et al. 2015; Hung et al. 2015; Kalas et al. 2015; Millar-Blanchaer et al. 2016; Esposito et al. 2018; Hom et al. 2020) at inner working angles (IWA) and resolutions not accessible by previous generation instruments. Combined with an additional sample of newly resolved scattered-light debris discs from SPHERE and other instruments (e.g. Thalmann et al. 2013; Wahhaj et al. 2016; Engler et al. 2018, 2020; Bonnefoy et al. 2021) and previously resolved systems (e.g. Schneider et al. 1999; Liu 2004; Schneider, Silverstone & Hines 2005; Hines et al. 2007; Kalas, Fitzgerald & Graham 2007; Soummer et al. 2014; Choquet et al. 2016; Padgett & Stapelfeldt 2016), the ensemble of scattered-light imaged discs 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 discs, few group/population studies have been conducted. Ren et al. (2023) investigated debris disc colour through HST observations, identifying a pre-dominantly blue colour suggesting higher scattering efficiency at shorter wavelengths. Esposito et al. (2020) first reported on the statistics and properties of the 29 circumstellar discs imaged as a part of the Gemini Planet Imager Exoplanet Survey (GPIES; GS-2015A-Q-500, PI B. Macintosh). While the study reported and characterized the geometries of a few newly resolved discs, 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 disc 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 disc system (Milli et al. 2012), biasing morphological characterization efforts. Different approaches to disc modelling 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 disc is an infinitely narrow ring (e.g. Crotts et al. 2024). The treatment of disc 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 parametrization 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 disc 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 modelling approach to facilitate more direct comparisons between morphological properties of a sample of debris disc 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 disc targets. Rather than modelling 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 modelling results of eight debris disc targets imaged by GPI. In Section 2, we describe the samples from which our observations are derived and the additional criteria employed to construct our final sample. In Section 3, we describe the properties of the observations of our target sample. In Section 4, we describe the data reduction approach for all data sets in our sample. In Section 5, we describe our forward modelling setup. In Section 6, we present the images and constrained parameters of our model analysis. In Section 7, we describe the interpretations of our results in modelling and understanding debris disc systems. In Section 8, we summarize our findings and discuss future implications for the results of our analysis.

2 TARGET SAMPLE

Selected targets in this analysis originate from three distinct observational programs. The majority of the data originate from the GPIES campaign (GS-2015A-Q-500, PI B. Macintosh), a 600-star direct imaging survey utilizing the 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 discs as well as stars with notable IR-excesses exceeding 10⁻⁵ 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 Disc 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 discs at higher contrast (∼10⁻⁶) and smaller angular separations (tens of mas) than were possible with previous generation instruments.

From the extensive GPIES Disc and GPI-LLP Disc samples, eight targets were selected. These targets satisfied three main criteria: (1) data from a GPI spectral mode data set, (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 Section 5.2, a simulated PSF core unique to each data set is necessary for our analysis, as PSF structure and intensity can vary with observing conditions and target brightness. In GPI data sets, this PSF core is modelled 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 data sets, these spots have well-defined shapes and flux ratios. Although GPI has resolved an extensive sample of discs in polarized intensity (Esposito et al. 2020; Crotts et al. 2024), they are not considered in this analysis. In polarimetric data sets, the satellite spots necessary for our PSF core modelling routine are elongated and cannot be used to accurately generate a model PSF core for disc forward modelling. Finally, although polarized intensity phase functions contain notable features among scattered-light resolved discs, 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 scattered-light disc can vary distinctly between systems and would introduce an additional parameter space of investigation that is not relevant to studying morphological properties of discs (e.g. Hadamcik & Levasseur-Regourd 2003; Hadamcik et al. 2007). 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 Disc Sample, only K1-band observations of HD 32 297 and HR 4796A were considered due to the SNR threshold we applied. Summary information regarding our sample targets is shown in Table 1.

Table 1.

Target List summary. H and K magnitudes originate from 2MASS photometry (Cutri et al. 2003). Distances are retrieved from Gaia DR2 (Gaia Collaboration et al. 2016, 2018). Other references: (1) Nielsen et al. (2019), (2) Torres et al. (2006), (3) Schneider, Silverstone & Hines (2005), (4) Bell, Mamajek & Naylor (2015), (5) Houk & Smith-Moore (1988), (6) Torres et al. (2008), (7) Soummer et al. (2014), (8) Pecaut & Mamajek (2016), (9) Houk & Cowley (1975), (10) de Zeeuw et al. (1999), (11) Kalas et al. (2015), (12) Houk (1978), (13) Kasper et al. (2015), (14) Padgett & Stapelfeldt (2016), (15) Wahhaj et al. (2016), (16) Thalmann et al. (2013), (17) Houk (1982), (18) Schneider et al. (1999).

Name	H	K	d	Age	Spectral	Moving	First resolved
	(mag)	(mag)	(pc)	(Myr)	Type	Group	Detection
HD 32 297	7.6	7.6	132.79 ± 1.06	15–45 (1)	A0V (2)	None	(3)
HD 35 841	7.8	7.8	103.68 ± 0.30	38–48 (4)	F3V (5)	Columba (6)	(7)
HD 106 906	6.8	6.7	103.33 ± 0.46	12–18 (8)	F5V (9)	LCC (10)	(11)
HD 110 058	7.5	7.6	129.98 ± 1.33	12–18 (8)	A0V (12)	LCC (10)	(13)
HD 111 520	7.7	7.7	108.94 ± 0.65	12–18 (8)	F5/6V (12)	LCC (10)	(14)
HD 114 082	7.2	7.2	95.65 ± 0.45	12–18 (8)	F3V (9)	LCC (10)	(15)
HD 146 897	7.8	7.8	131.50 ± 0.93	7–13 (8)	F3V (5)	US (10)	(16)
HR 4796A	5.8	5.8	72.78 ± 1.75	7–13 (4)	F3V (17)	TWA (4)	(18)

3 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 Disc 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 disc 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.

Table 2.

Summary of observations for the sample. Seeing estimates are not available for HD 114 082 and HD 146 897 observations, as the DIMM seeing monitor at Gemini-South was inoperable at the time of observation.

Name	Filter	texp	N	ΔPA	Air Mass	Seeing	Date	Program
		(s)		(°)		(arcsec)
HD 32 297	H	59.65	38	16.7	1.26–1.29	0.63	2016 Dec 20	GS-2015B-Q-500
HD 32 297	K1	88.74	51	33.1	1.26–1.30	0.72	2015 Nov 30	GS-2015B-LP-6
HD 35 841	H	59.65	50	46.9	1.01–1.06	1.01	2016 Feb 28	GS-2015B-Q-500
HD 106 906	H	59.65	42	25.3	1.11–1.12	0.86	2015 May 4	GS-2015A-Q-500
HD 110 058	H	59.65	38	29.6	1.06	0.71	2016 Mar 19	GS-2015B-Q-500
HD 111 520	H	59.65	42	34.8	1.06–1.07	1.01	2015 Jul 02	GS-2015A-Q-500
HD 114 082	H	59.65	47	25.8	1.16	–	2018 Jan 29	GS-2017B-Q-500
HD 146 897	H	59.65	38	59.5	1.01–1.02	–	2018 Aug 15	GS-2017B-Q-500
HR 4796A	H	59.65	37	53.0	1.01–1.02	0.73	2016 Mar 18	GS-2015B-Q-500
HR 4796A	K1	88.74	46	78.5	1.01–1.09	0.55	2015 Apr 3	GS-2015A-Q-27

For all observations, the primary star of each system was centred behind a focal plane mask (FPM), allowing for higher contrast to be reached in the immediate area surrounding the star. The IWA of the GPI FPM is ∼0|${_{.}^{\prime\prime}}$|1, and the field of view (FOV) of GPI is 2.8 × 2.8 arcsec². GPI has a spatial sampling of ∼14 mas pixel⁻¹. For K1-band observations, five sky exposures every hour for each target were acquired for thermal/sky background subtraction.

Although the targets were observed in GPI’s spectral mode that provides spatially resolved low-resolution spectra, spectral properties of targets are not investigated in this work and the full set of wavelength channels per science image was collapsed into single broadband images to increase disc SNR. Faint disc surface brightness in individual wavelength channels can be challenging to measure, and Esposito et al. (2020) notes that the measured H-band spectra of GPI-imaged debris discs tend to be relatively featureless. Further details of H-band observations are given in Esposito et al. (2020).

4 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 Section 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 data cube 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 centre located at the centre 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, Pueyo & Larkin 2012) in combination with ADI. ADI takes advantage of the fact that instrumental PSF artefacts do not rotate with the FOV in pupil-stabilized observations. Any astrophysical object within the FOV will rotate with respect to the centre 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 discs) 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 (N_δ 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 disc SNR, and minimizing overlap of the disc 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 disc forward models for later analysis, as described in Section 5.

Table 3.

KLIP reduction parameters for all targets, chosen to maximize disc SNR.

Name	Filter	KL modes	Nδ (°)
HD 32 297	H	3	4
HD 32 297	K1	5	4
HD 35 841	H	7	6
HD 106 906	H	5	6
HD 110 058	H	5	6
HD 111 520	H	5	6
HD 114 082	H	5	6
HD 146 897	H	5	6
HR 4796A	H	7	6
HR 4796A	K1	5	6

5 MODEL METHODOLOGY

As high-contrast images of debris discs reduced with ADI often suffer from severe self-subtraction (Milli et al. 2012), forward-modelling is often necessary to constrain the properties of a system. The key steps in the forward modelling methodology employed in this study are shown in Fig. 1. To assess morphological properties, we adopt a forward modelling approach using DiskFM (Mazoyer et al. 2020) for comparison with the science data. Disc 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 disc models per target.

Figure 1.

A diagram demonstrating the forward modelling 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 disc 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 modelling of a disc given a set of inputs. A χ² function is used to inform the MCMC sampler of forward model minimization. The outputs of DiskFM include the maximum likelihood best-fitting disc model and posterior probability density distributions of free parameters (green boxes). Purple boxes represent calculation procedures internal to DiskFM.

5.1 Disc model

The disc 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 η(r, z) follows a smoothly connected two power-law structure described in more detail in Augereau et al. (1999) and given in equation (1):

$$\begin{eqnarray} \eta (r, z) \propto R(r) Z(r, z) \end{eqnarray}$$

(1)

where R(r) is given as:

$$\begin{eqnarray} R(r) \propto \left\lbrace \left(\frac{r}{R_{\rm C}}\right)^{-2\alpha _{\rm in}} + \left(\frac{r}{R_{\rm C}}\right)^{-2\alpha _{\rm out}}\right\rbrace ^{-\frac{1}{2}} \end{eqnarray}$$

(2)

where r is the radial distance from the host star and R_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 disc, respectively. Z(r, z) is given as

$$\begin{eqnarray} Z(r,z) \propto \exp \left(-\left(\frac{|z|}{h(r)}\right)^{\gamma _{\rm vert}}\right) \end{eqnarray}$$

(3)

where z is the distance from the disc mid-plane, γ_vert dictates the shape of the vertical density distribution, and h(r) is the height above the disc mid-plane as a function of r. 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 h(r) is given by a_r × r, where a_r is the constant aspect ratio h/r, assuming a ‘bow-tie’ shape for every disc. Both inner (R_in) and outer (R_out) cutoff radii define where the dust density η(r, z) is treated as zero; only R_in is treated as a free parameter. All discs except for HD 110 058 have evidence of extensive halo emission observed at larger FOV than GPI (e.g. Schneider, Silverstone & Hines 2005; Soummer et al. 2014), therefore R_out is set to be the outer working angle of the GPI FOV for these systems. For HD 110058, R_out is loosely defined to exist outside of the MCMC prior range for R_C. The proportionality constant is not solved for in our 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 third-order spline function.

One of the premises of this study is the observation by Hughes, Duchêne & Matthews (2018) that scattering in most debris discs 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, Millis & Birch 1998; Moreno et al. 2012; Hui 2013). To avoid one particular data set 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 114 082 from Engler et al. (2023). Finally, we fitted a ninth-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 Fig. 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).

Figure 2.

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

We note that most of the SPFs used in this process were observed in the optical (the main exception being the J 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 discs.

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, Duchêne & Matthews (2018) and was used as a test case for the HD 114 082 and HR 4796A data sets to demonstrate the differences in best-fitting model appearances and noise-scaled residuals (hereafter called residuals) when using different SPFs. Although SPFs have been measured for HD 32 297 (Duchêne et al. 2020), HD 35 841 (Esposito et al. 2018), and HD 114 082 (Engler et al. 2023), they are not considered for this analysis as their shapes were found to be largely similar to the trends observed in Hughes, Duchêne & Matthews (2018). SPFs for the other targets in this sample have not been measured.

Physical disc models are axisymmetric and generated in a three-dimensional 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 disc model is properly rotated, the three-dimensional disc model cube is collapsed along the line-of-sight direction, creating a two-dimensional image of the disc.

5.2 Forward modelling process

DiskFM (Mazoyer et al. 2020) was used to perform forward modelling on all disc 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 disc 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 disc 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 data set 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 artefacts than the reduced science data set.

To generate data set-specific noise maps for calculating forward modelling 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 disc signal overlaps between individual science frames. PyKLIP-ADI is then utilized on the randomly rotated science data set to eliminate traces of the disc signal entirely and create a noise map that has been reduced with the same KLIP parameters as the non-rotated science data set. 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 disc parameters, although they only tested this approach for a small sample of model systems. As shown in Section 6, 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 |$\chi ^2_{\rm red}$|⁠. This noise scaling factor is likely not consistent between different data sets, but we retain the factor of 3 to be consistent to previous high-contrast imaging studies and so that all disc analyses are uniform. A synopsis of the steps and components of DiskFM is shown in the gold and purple boxes of Fig. 1.

5.3 Likelihood calculation

The forward model is compared to the data by measuring

$$\begin{eqnarray} \chi ^2 = \sum _S \frac{(Data-Forward\,\,Model)^2}{Uncertainty^2} \end{eqnarray}$$

(4)

where S 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 |$e^{-\chi ^2/2}$| in a masked region unique to each disc in the sample for both the disc forward model and science data. This calculation is performed until chain convergence. For all disc target analyses, 120 walkers were used, calculated for at least 5000 iterations. More iterations were added as necessary to confirm converged behaviour within the MCMC chains. After at least 300 iterations of converged behaviour, 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 disc flux is most apparent and some background with no apparent disc signal is present. For systems where we tentatively resolve the back sides of discs – HD 35841, HD 106906, and HR 4796A – concentric rings enclosing the spine of the disc are created, with some space given between the apparent edge of the disc and the edge of the mask for the inclusion of the noise background in the likelihood calculation. If the front side of the disc 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 Section 5.2. Although spatially and spectrally correlated noise among pixels is a concern for direct imaging IFS data sets (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.

Disc models have either eight or nine physical free parameters depending on disc inclination and apparent disc thickness. The eight free parameters common to all discs are the inner cutoff radius R_in, critical radius R_C, surface dust density inner power law index α_in, surface dust density outer power law index α_out, aspect ratio a_r, inclination i, position angle PA, stellocentric offset in the disc plane along the projected major axis dy (where positive dy corresponds to a stellocentric offset in au in the disc 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 discs where we resolve the front and back sides of the disc, an additional free parameter dx is parametrized (where positive dx corresponds to a stellocentric offset in au in the disc plane pointing toward the observer along the minor axis), describing the stellocentric offset in the disc plane along the projected minor axis. Although we do not resolve the back side of HD 106906, we also parametrize dx 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 dx.

Table 4.

MCMC initial parameters and unique prior ranges for each target. For all values of α_in, the prior range extends from [0...10]. For α_out, the prior range is from [−10...0]. For all values of a_r, the prior range is [0.001...0.3]. For all values of dx and dy, the prior range extends from [−5...5] au. References for initial parameter selection: (1) Duchêne et al. (2020), (2) Esposito et al. (2018), (3) Lagrange et al. (2016), (4) Kalas et al. (2015), (5) Kasper et al. (2015), (6) Esposito et al. (2020), (7) Draper et al. (2016), (8) Wahhaj et al. (2016), (9) Engler et al. (2017), (10) Thalmann et al. (2013), and (11) Perrin et al. (2015).

Name	Rin	RC	αin	αout	ar	i	PA	dx	dy
	(au)	(au)				(°)	(°)	(au)	(au)	Pos. dy
HD 32 297	50 [25...100] (1)	98.4 [75...150] (1)	5	−5	0.001	88.4 [80...90] (1)	47.9 [45...55] (1)	–	0.1	SW
HD 35 841	60.3 [40...70] (2)	60.3 [50...90] (2)	5	−5	0.01	84.9 [80...90] (2)	165.8 [160...170] (2)	0.1	0.1	SE
HD 106 906	66.6 [40...80] (3)	72.3 [50...90] (3)	5	−5	0.001	84.6 [80...90] (4)	284.2 [280...290] (4)	0.1	0.1	NW
HD 110 058	39.0 [20...60] (5)	39.0 [20...80] (5)	5	−5	0.1	84.0 [80...90] (6)	155.0 [150...160]a (5)	–	0.1	SE
HD 111 520	71.0 [50...90] (7)	81.0 [60...100] (7)	5	−5	0.001	88.0 [80...90] (7)	165.0 [160...170] (7)	–	0.1	SE
HD 114 082	28.7 [10...50] (8)	30.7 [15...70] (8)	5	−5	0.01	83.3 [80...90] (8)	105.7 [100...110] (8)	–	0.1	SE
HD 146 897	50.0 [25...100] (9)	65.0 [40...115] (9)	5	−5	0.1	84.0 [80...90] (10)	113.9 [110...120] (6)	–	0.1	SE
HR 4796A	74.4 [60...90] (11)	78.5 [70...100] (11)	5	−5	0.01	76.5 [70...80] (11)	26.1 [20...30] (11)	0.1	0.1	SW

^aThis prior was expanded from [150...160] to [150...165] for models fitted to the NW side of the HD 110 058 disc only, after initial analysis suggested a PA greater than the prior limit.

MCMC initial parameters for R_in, R_C, i, and PA are set as the constrained parameters of these disc targets as described in previous studies of scattered-light resolved imaging of discs, 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 a_r parameters for HD 35841, HD 114082, HD 146897, and HR 4796A are derived to be in similar orders of magnitudes of previous studies (Milli et al. 2015; Esposito et al. 2018; Goebel et al. 2018; Olofsson et al. 2022; Engler et al. 2023). For HD 32297, HD 106906, HD 110058, and HD 111 520, we opt to choose initial a_r 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 disc flaring index, and the outer radial cutoff R_out are fixed. A summary of initial MCMC parameters and prior ranges are shown for each disc in Table 4. A summary of all DiskFM inputs (reduced science data set, free parameters, fixed parameters, and the SPF) are shown in the blue boxes of Fig. 1. The outputs of DiskFM, shown in green boxes, include the maximum likelihood best-fitting model and posterior probability density distributions of free parameters.

Three disc systems – HD 106906, HD 110058, and HD 111520 – have 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 disc model utilized in this analysis does not have any parametrization that accounts for the strong disc 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 Section 6.1, strong residuals are present in the best model fits. As a result, we choose to conduct independent disc forward-modelling analyses on each side of these asymmetric discs (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 Section 5.1 to debris disc systems in general.

6 RESULTS

6.1 Uniform disc model results

The best-fitting models for each disc based on the uniform approach with the generic SPF, a continuous likelihood mask, and their residuals are shown in Figs 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.¹ Posterior distributions for all model analyses are reported in Appendix  A. Fig. 3 shows best-fitting models of HD 32297, HD 35841, and HD 106906. Fig. 4 shows the best-fitting models of HD 110058, HD 111520, and HD 114082. Fig. 5 shows the best-fitting models of HD 146 897 and HR 4796A.

Figure 3.

Modelling 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-fitting 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 32 297 best-fitting model are likely due to the inherent brightness of the disc, creating strong self-subtraction wings in the KLIP-ADI reduction and inducing non-linearity in the pyKLIP forward modelling process. In the HD 35 841 best-fitting 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 106 906 model, the residual map appears to have structure along the spine of the disc, particularly on the SE side.

Figure 4.

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

Figure 5.

Same as Fig. 3 but for two more targets in the sample. The structured and high residuals present in the HD 146 897 residual map are likely due to diffuse PSF halo structure biasing the forward modelling 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.

Table 5.

MCMC median likelihood parameters with ±1σ error bars and 3σ upper/lower limits.

Name	Rin	RC	αin	αout	ar	i	PA	dx	dy	|$\chi ^{2}_{\rm red}$|
	(au)	(au)				(°)	(°)	(au)	(au)
Generic scattering phase function
HD 32 297	|$45.24^{+0.64}_{-0.94}$|	>142.54	<0.26	|$-7.40^{+0.43}_{-0.41}$|	0.01 ± 0.001	87.63 ± 0.03	47.29 ± 0.01	–	−0.46 ± 0.15	2.20
HD 32 297 (K1)	|$40.34^{+1.34}_{-1.51}$|	>142.07	<0.27	−5.07 ± 0.29	<0.01	|$87.58^{+0.05 }_{ -0.04}$|	47.13 ± 0.02	–	0.89 ± 0.28	0.86
HD 35 841	<55.97	|$62.15^{+1.32 }_{ -1.61}$|	|$2.28^{+1.18 }_{ -0.79}$|	<−6.09	|$0.04^{+0.004 }_{ -0.005}$|	|$82.69^{+0.30 }_{ -0.28}$|	165.49 ± 0.15	|$1.94^{+1.11 }_{ -1.09}$|	|$0.04^{+0.49 }_{ -0.48}$|	0.54
HD 106 906	<49.02	<75.88	|$3.28^{+2.23}_{-1.66}$|	|$-1.97^{+0.22 }_{ -0.29}$|	|$0.03^{+0.004 }_{ -0.005}$|	|$83.94^{+0.21 }_{ -0.19}$|	|$285.06^{+0.08 }_{ -0.09}$|	>4.16	|$-2.38^{+0.49}_{-0.52}$|	0.43
HD 110 058	<42.18	<67.10	Unconst.	|$-3.07^{+0.47 }_{ -0.68}$|	0.12 ± 0.01	>84.62	|$158.67^{+0.27}_{-0.26}$|	–	|$-2.04^{+0.59}_{-0.62}$|	0.18
HD 111 520	<54.13	|$72.12^{+7.48 }_{ -7.12}$|	<6.77	|$-4.19^{+0.63 }_{ -0.96}$|	0.07 ± 0.004	>89.82	166.40 ± 0.09	–	<−4.64	0.33
HD 114 082	|$31.11^{+0.26 }_{ -15.10}$|	<33.51	>0.04a	|$-2.99^{+0.12 }_{ -0.13}$|	0.06 ± 0.004	|$83.07^{+0.19 }_{ -0.71}$|	106.16 ± 0.07	–	0.74 ± 0.11	0.90
HD 146 897	<49.25	|$50.10^{+1.03 }_{ -6.87}$|	>0.75a	|$-2.45^{+0.12 }_{ -0.13}$|	<0.04	|$82.51^{+0.27 }_{ -0.26}$|	115.00 ± 0.17	–	|$0.69^{+0.35 }_{ -0.34}$|	1.09
Half mask model analysis
HD 106 906 (NW)	<43.24	|$67.84^{+7.95 }_{ -8.59}$|	<4.36	|$-2.73^{+0.55 }_{ -0.67}$|	0.06 ± 0.005	83.07 ± 0.22	287.28 ± 0.14	<−1.92	>2.23	0.37
HD 106 906 (SE)	|$55.21^{+1.48 }_{ -5.04}$|	<74.60	Unconst.	|$-2.56^{+0.27 }_{ -0.32}$|	<0.01	85.21 ± 0.11	283.32 ± 0.09	<−1.40	<−3.75	0.20
HD 110 058 (NW)	<40.57	<68.37	<9.98a	|$-2.35^{+0.59 }_{ -1.00}$|	<0.14	|$85.05^{+0.64 }_{ -0.77}$|	|$159.97^{+0.52}_{-0.56}$|	–	<4.60	0.29
HD 110 058 (SE)	|$34.83^{+3.64 }_{ -5.44}$|	<57.59	Unconst.	|$-3.77^{+0.60 }_{ -0.76}$|	|$0.12^{+0.02 }_{ -0.01}$|	>84.17	|$157.96^{+0.34 }_{ -0.35}$|	–	|$0.08^{+2.98 }_{ -2.91}$|	0.14
HD 111 520 (NW)	<56.64	<90.35	<9.89a	|$-3.39^{+0.31 }_{ -0.49}$|	|$0.06^{+0.004 }_{ -0.003}$|	>89.62	167.08 ± 0.08	–	|$-1.11^{+2.24 }_{ -2.21}$|	0.17
HD 111 520 (SE)	<74.06	<92.60	Unconst.	|$-4.33^{+0.80 }_{ -1.02}$|	0.07 ± 0.01	>88.11	164.00 ± 0.23	–	Unconst.	0.14
Milli et al. (2017) scattering phase function
HR 4796A	|$73.86^{+0.32 }_{ -0.09}$|	|$73.79^{+0.54 }_{ -2.11}$|	>2.41	<−9.81	<0.004	76.89 ± 0.03	26.75 ± 0.02	−2.20 ± 0.09	0.84 ± 0.05	0.41
HR 4796A (K1)	74.43 ± 0.10	<73.05	>1.96	<−9.78	<0.007	76.95 ± 0.05	26.55 ± 0.04	−2.41 ± 0.16	1.79 ± 0.08	0.21

^aThese upper/lower limits span most of the prior range but have low-probability tails.

Tight constraints (1σ < 0.3°) on i 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 i solution. PA was also tightly constrained (1σ < 0.3°) for all systems. In general, i and PA are typically the easiest parameters to constrain in a debris disc model, mostly due to clearly resolved positions.

Our constraints on the radial dust density distribution suggest two families of results: discs where the radial dust density distribution can be described by one power law (R_C < R_in or R_C > R_out) and discs where the radial dust density distribution can be described by two power laws (R_in < R_C, R_C and α_in and/or α_out well-defined). Inner disc properties (R_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 disc 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 discs around HD 35 841 and HD 111 520 can be modelled with two power laws, while the dust density distribution of discs around HD 106906, HD 110058, HD 114082, HD 146897, and HR 4796A can be modelled with a single power law.

1σ constraints were also found for the aspect ratios of all systems except for HD 146 897 and HR 4796A, for which we achieve upper limits. The less-inclined and compact nature of HD 146 897 as well as the bright PSF halo feature overlapping the front side of the disc can bias estimations of inclination and vertical extent. As HR 4796A is also relatively bright, non-linearity in the KLIP reduction may also similarly bias estimates of some structural properties of the disc such as the vertical density distribution.

Narrow constraints (1σ < 0.3 au) on dy were identified for all systems except for HD 111 520 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 dx 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 (Lagrange et al. 2009; Kalas et al. 2015; 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 (dy) and minor (dx) axes of the disc. However, residuals are still present, suggesting that stellocentric offsets alone may not be enough to account for such asymmetries.

We also calculate |$\chi ^{2}_{\rm 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 disc signal in either model or data, naturally reducing |$\chi ^{2}_{\rm 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 data sets. 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 35 841 and HD 114 082 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 disc signal.

For the asymmetric systems HD 106906, HD 110058, and HD 111520, structured, positive residuals (≳ 1σ) appear in regions that overlap with disc signal. In the case of HD 106906, DiskFM prefers solutions favours a best-fitting model that places the centre of the ring SE of the FPM. For HD 111520, the best-fitting 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 best-fitting model itself appears to have an offset in PA compared to the perceived location of the disc. 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|${_{.}^{\prime\prime}}$|2 of the inner edge of the likelihood calculation mask, centred at the mid-plane of the disc. This region overlaps with the highest SNR values in the reduced data, and are likely inducing non-linearity in the KLIP-ADI reduction (Pueyo 2016), breaking down the linear expansion assumed in the forward modelling process. Regions where the absolute value of the residuals exceed 1 are also present far from the star and throughout most of the mask.Figure 7

For HD 146897, significant positive and negative residuals are structured on the front and back sides of the disc 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 disc signal. Because of the halo structure biasing the brightness of the front side of the disc, 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 azimuthally structured 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 disc and high negative residuals present on the front side of the disc. This behaviour 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 disc and little back-scattering comparatively. Strong positive residuals (Res. ≳ 2σ) are shown on the back side and ansae of the HR 4796A disc.

6.2 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 (dx and dy) 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 discs (HD 106906, HD 110058, and HD 111520), treating each each side of the discs separately, as seen in Figs 6, 7, 8. These test cases assume that the disc 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 disc is not being considered. For example, disc spine curvature is typically well-defined even when only considering half of the disc. For these analyses, the initial MCMC parameters and priors are consistent when modelling both sides of each disc, with the exception of the NW side of HD 110058, where initial investigations of the system preferred model solutions in the PA 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 PA 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 disc that may appear to have different morphological properties compared to the other side. We chose not to conduct a half-mask modelling analysis for HR 4796A despite its asymmetric nature, as the median likelihood value of dy corresponds to an offset less than one pixel in length.

Figure 6.

Best-fitting modelling results for separate analyses of the NW and SE sides of HD 106906. The models are distinctly different, most notably seen in the PA, i, and vertical structure profiles of both best-fitting disc models.

Figure 7.

Best-fitting modelling results for separate analyses of the NW and SE sides of HD 110058. Differences in the models are subtle but most notably seen in the PA and stellocentic offsets of both best-fitting disc models.

Figure 8.

Best-fitting modelling results for separate analyses of the NW and SE sides of HD 111520. The models are drastically different, especially in terms of overall brightness. Other notable differences include the PA and overall radial extent of the models.

When determining best-fitting models for each side separately, we are able to achieve lower residuals without systematic patterns of disc 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 110 058 and HD 111 520 where the entire likelihood mask is used. These two systems do not appear to have morphological asymmetries as significant as HD 106906, with HD 110 058 containing diffuse asymmetrical features around the ansae and HD 111 520 containing a strong brightness asymmetry and tentative warp on the SE side (Draper et al. 2016).

In modelling the two sides of HD 106 906 separately, we are able to constrain positional and outer disc properties and achieve upper/lower limits for stellocentric offsets. The most significant differences are evident in the median likelihood values of i, PA, a_r, and dy. While i may be somewhat degenerate with a_r, 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-fitting models for the disc even when analysing both sides separately. To more robustly characterize the morphology of this system, a disc model that can create truly eccentric discs 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 modelling 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 PA, likely due to the asymmetric warps on both sides of the disc (López et al. ; Stasevic et al. 2023). Appendix B3 contains a more detailed discussion of the properties of the HD 110 058 disc and comparisons to previous literature.

In the two-sided model analysis of the HD 111 520 disc, residuals no longer contain structure and are less (|Res.| ≲ 1σ) compared to modelling the whole disc. 1σ constraints in α_out, a_r, and PA were found for both sides. Constrained parameters for both analyses are consistent with each other except for PA, likely due to the tentative warp observed on the SE side.

While constraints on some density properties of modelling either side of a disc are consistent with each other, namely α_out and a_r, the greatest differences are highlighted in the constraints of PA and i (and a_r 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.

6.3 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-fitting model is shown in Fig. 9.

Figure 9.

Best-fitting modelling results for analysis of HR 4796A and HD 114 082 with the measured SPF from Milli et al. (2017). While we achieve a lower-residual fit for HR 4796A, positive residuals are still present particularly along the front and back sides of the disc on the Northern sides. While a solution with a stellocentric offset is preferred in the MCMC analysis, the presence of relatively strong residuals may be explained by a dust density enhancement as observed in Olofsson et al. (2019), Milli et al. (2019), and Chen et al. (2020). To achieve a reasonable morphological model for HD 114082, the PA of the disc had to be increased by 180°, as DiskFM prefers to align the brighter back side of the disc model with the apparent front side of the disc from the data. The forward-scattering lobe of the SPF located on the front side of the disc is also behind the FPM in this data set, further encouraging this behaviour.

The improvement in the best-fitting model for HR 4796A is evident from the significantly lower residuals of the best-fitting model, particularly at the back side and ansae of the disc (|Res.| ≲ 2σ), in addition to the relatively lower |$\chi ^2_{\rm red}$| (0.41 for the Milli et al. 2017 SPF compared 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 discs can still be determined with a fixed SPF. Additionally, all median likelihood parameters and limits except for the aspect ratio were consistent between both modelling 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 disc system. Our findings suggest much steeper radial profiles than our prior limits, which is consistent with the sharply defined edges observed in the disc. The vertical structure profile of the disc 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 114 082 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 PA, as discussed in this section). Given the high degree of back-scattering in the HR 4796A SPF compared to our generic SPF, the faint and marginally resolved back side of the disc, and the forward scattering peak of the disc being located behind the FPM, a reasonable model for this system would require ‘flipping’ the orientation of the disc such that the apparent front side of the disc 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 PA, simulating placing the front side of the disc on the North side. While we achieve a reasonable maximum likelihood model (Fig. 9), the fit is still poor compared to our initial analysis with the generic SPF and the front side of the disc on the South side, with a |$\chi ^2_{\rm 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 disc on the Northern side of the image, inconsistent with previous studies of the system (Engler et al. 2023; Crotts et al. 2024). Additionally, we find that most constrained parameters are consistent between the generic and HR 4796A SPF analyses, aside from i, α_in, and a_r.

Both of the HD 114 082 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 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.

6.4 K1-band modelling results

For the two systems with K1-band data sets meeting our SNR threshold, best-fitting model results are shown in Fig. 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).

Figure 10.

Best-fitting modelling results for analyses of K1-band images of HD 32 297 and HR 4796A. The generic SPF was used for HD 32 297 and the Milli et al. (2017) SPF was used for HR 4796A. Residuals appear cleaner in both systems compared to their H-band images, potentially due to the fainter disc brightness at K. Our best-fitting model for HD 32 297 still has some significant residuals close to the edge of the FPM, and non-linearity in the KLIP reduction may still be evident to a lesser extent.

Compared to the H-band reduced image of HD 32 297 and best-fitting model, the K1-band reduced image of HD 32 297 and best-fitting model appear less affected by the effects of self-subtraction from PSF subtraction and non-linearity in forward modelling, although high residuals (|Res.| ≳ 2σ) are still present. As the emission from the disc at K1-band is relatively less than emission at H-band, self-subtraction and non-linearity in the forward modelling process is likely to be less severe. Additionally, similarities in relative errors among parameters with a lower disc SNR could suggest lower residual amplitudes. The broader PSF at K1-band may also contribute to the smoothing of finer structural features in the disc, allowing the model to more easily match the data. Regardless, our K1-band model fit still has significant residuals within 0|${_{.}^{\prime\prime}}$|1 of the inner edge of our likelihood mask.

The HR 4796A best-fitting 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 disc, particularly around the ansae and the front side of the disc. 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 disc 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 R_C > R_in, and the lack of distinction between R_in and R_C suggests that the inner radial profile of the disc is sharply defined. Similar to the H-band result, our posteriors suggest steep radial density profiles.

7 DISCUSSION

7.1 Limitations of the forward-modelling approach

From our analysis, not all systems yield low residual model fits. The linear approximation in equation 4 of 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 modelling non-linearity can be worse for bright and extended discs. Non-linearity 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 disc in our sample. To mitigate this issue for future characterization of bright debris discs, 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 32 297 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 Sections 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 i and PA 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 mid-plane density profiles of previous studies and this work. The uniqueness of approaches and associated limitations in analysing discs also prevent direct comparisons of structural properties, further complicating ensemble studies of debris discs.

Finally, the manner in which disc 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 disc image modelling (e.g. Esposito et al. 2018; Chen et al. 2020; Duchêne et al. 2020; Crotts et al. 2021; Olofsson et al. 2022; Engler et al. 2023), allowing a more thorough search of parameter spaces.

As we obtain more high-resolution images of debris disc systems, consistent methodologies (e.g. Olofsson et al. 2022; Crotts et al. 2024) are needed to ensure more direct comparisons between systems with unique properties such as host star spectral type and age. While a uniform modelling 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 modelling and that true astrophysical trends can still be identified.

7.2 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 data set 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-fitTING 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.

7.3 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 a_r appears to suggest extremely narrow profiles for HR 4796A [vertical full width at half-maximum (FWHM) ∼ 0.1 au at R_C], unlike what was found in Olofsson et al. (2022; vertical FWHM ∼ 4 au at R_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 while 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 disc 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 artefact 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 discs, however.

Olofsson et al. (2022) investigated the effects of gaseous components of discs 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’ discs as gas drag boosts the efficiency of vertical settling. Of the eight targets investigated in our sample, two discs (HD 32297; Greaves et al. 2016, HD 110058; Hales et al. 2022) have measured CO gas detections and one disc (HD 146897; Lieman-Sifry et al. 2016) has a tentative CO gas detection. The largest modelled aspect ratios appear to be around the star HD 110058, where a_r was found to be greater than 10  per cent in modelling both the entire disc and the SE side independently. This may seem contradictory to the idea that a gas-bearing disc would likely be thinner, but the distance and compact nature of HD 110 058 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 disc 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 R_in and R_out on either side and divided by the radius of maximum dust density. In Fig. 11, we plot a_r as a function of the relative radial width. In a study of ALMA-observed debris discs, Terrill et al. (2023) identified a tentative trend of increasing relative radial width and aspect ratio, with HD 110 058 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, Naoz & Fitzgerald 2017), primordial eccentricities from instabilities (Kennedy 2020), and gas/dust interactions (Lyra & Kuchner 2013). As a general case, Chiang et al. (2009) and Rodigas, Malhotra & Hinz (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.

Figure 11.

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

With the exception of HD 110058, most of our constrained disc 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 disc aspect ratios can be enhanced most strongly (∼0.05) for multiple-planet systems where the planets lie interior to the disc. 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, Naoz & Fitzgerald (2017) found that increasing the mutual inclination of an external perturbing companion increases the vertical extent of a disc 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 110 058 (Hales et al. 2022) could also contribute to the radially narrow and warped geometry observed. López et al. () also found that both interior and exterior planetary companions can produce the geometry observed in HD 110 058 through secular perturbations, while Stasevic et al. (2023) concluded that only an interior companion could produce the observed warped geometry. Interestingly, the HD 106 906 median likelihood whole disc and NW side radial widths are fairly broad (∼50 au), despite its eccentric shape likely being caused by an external perturber (Nesvold, Naoz & Fitzgerald 2017). Broad radial widths and heightened aspect ratios both would support dynamical models of planet scattering (Nesvorný 2015), although our aspect ratios for HD 106 906 are not large relative to the rest of our sample. For more vertically thin discs, self-stirring and/or secular interactions may be dominant (Matrà et al. 2019); this behaviour may be relevant for systems such as HD 35 841 and HD 146897.

7.4 Disc eccentricity

Disc eccentricity in our model can only be parametrized in the form of stellocentric offsets. In general, the most significant offsets we observe are seen in the best-fitting 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).² Pericentre glow (Wyatt et al. 1999) is postulated as a cause of the asymmetry in HR 4796A, but despite inducing a pericentre 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 disc (e.g. HD 106 906 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 106 906 (see Jones et al. 2023 and references therein). Studies of HD 110 058 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. (in preparation) provide tentative evidence of a vertical stirring mechanism such as a companion. Dedicated dynamical sculpting studies of HD 111 520 have not been conducted, and the extremely edge-on inclination further limits our ability to assess the true eccentricity of the disc. Giant impact modelling 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 Fig. 12, we compare median likelihood eccentricities of all disc modelling 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 disc 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 a_r we observe. We can compare our results to findings from Kennedy (2020), who observed that the highly eccentric debris discs around Fomalhaut and HD 202 628 appear narrower than expected from secular perturbation models and assumptions of zero initial eccentricity. The dashed line in fig. 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 disc 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.

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 fig. 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.

7.5 Implications for grain properties and scattering phase functions

Efforts in modelling the SPF for debris disc systems most commonly utilize the HG phase function (Henyey & Greenstein 1941) as an analytical (but not physical) model to all scattering properties or fully parametrize 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 behaviour and are inconsistent across multiple wavelengths and/or observing modes. In modelling analysis of HST-NICMOS observations of HD 181327, Schneider et al. (2006) identified inconsistencies between a minimum grain size (a_min) determined from observed mean asymmetry factor in the SPF and an a_min determined from disc colour 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 modelling 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. (2022) seems to contrast with many empirical SPF measurements, which show similar behaviour between different wavelength regimes (e.g. HR 4796A, see Chen et al. 2020), and multiwavelength measurements of SPFs are likely also needed in order to better constrain the grain properties of debris disc systems.

The limited number of debris disc SPF measurements (e.g. Graham, Kalas & Matthews 2007; Duchêne et al. 2020) have been shown to follow trends in shape seen in solar system dust environments (Hughes, Duchêne & Matthews 2018). This generic SPF shape has been reproduced with highly porous grains where the minimum grain size is greater than 1 μm assuming both Mie theory (Grynko, Jockers & Schwenn 2004) and clusters of complex aggregates (Bentley et al. 2016) suggesting that debris disc dust shares many similarities to cometary dust and that overall dust composition is not as important of a factor. Our generic SPF utilizes 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. Marcus 2007a, b; Schleicher & Bair 2011). 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 disc SPFs could incorporate these results in their analysis.

The achievement of low residual model fits for most of the debris disc systems is a result that supports findings and conclusions from Hughes, Duchêne & Matthews (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 disc models of the HD 35841, HD 106906, HD 110058, HD 111520, and HD 114 082 systems, utilizing the same SPF informed from solar system dust. Although we did not attain a low-residual best-fitting 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 disc geometry without having to parametrize 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 disc than what the generic SPF allows. Modelling 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 (>5 μm) 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 disc 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 114 082 further supports that SPF choice is an important consideration, as the best-fitting model solution preferred model solutions where the front side of the disc was on the Northern side, notably inconsistent with other studies of these systems and our own disc 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.

8 SUMMARY

With advances in direct imaging instrumentation, more scattered-light debris discs have been resolved at small separations than ever before. With such a large sample of imaged debris discs, consistent and uniform characterization approaches should be implemented in order to make more direct comparisons between systems. We have performed a uniform forward-modelling analysis of eight bright debris discs imaged with the GPI in total intensity light. Scattering properties were determined by empirically informed scattering phase functions and not from parametrizing HG functions or specific grain properties. From our results, we were able to identify two families of debris discs: one where the mid-plane density profile can be described by one power law (discs with extremely sharp inner edges) and the other where the mid-plane 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 modelling.

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 modelling results. The range of aspect ratios we find in our modelling analyses are mostly consistent with a ‘natural’ disc 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 modelling. The broader aspect ratio of HD 110 058 along with its narrow radial width and structural warps may be indicative of complex stirring mechanisms including perturbing planetary companions. We also identify asymmetric disc 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 modelling can be conducted with a uniform and consistent set of assumptions about grain properties. While it is unlikely that most discs 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 achievement of well-fitting models in most cases suggest links between porous aggregates, Solar System cometary dust, and extrasolar debris discs. 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 modelling 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 non-linearity and self-subtraction in bright discs. 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 disc scattering phase functions are needed in order to support and identify common trends and properties of grains in a variety of different environments.

ACKNOWLEDGEMENTS

The authors acknowledge Research Computing at Arizona State University for providing the computational HPC resources from the Sol supercomputer (Jennewein et al. 2023) that have contributed to the research results reported within this paper. We would like to acknowledge Jason Yalim from the Research Computing Core Facilities at Arizona State University for their assistance in optimizing HPC usage. The authors would also like to acknowledge the ASU Research Computing team as a whole for their assistance in setting up and assisting in all other matters related to high power computing, along with Dr Michael Line, Dr Aishwarya Iyer, and Peter Smith for their perspectives and conversations related to MCMC sampling and forward modelling practices. The authors also acknowledge Julien Milli for the measured HR 4796A SPF utilized in this analysis. Finally, the authors would like to thank and acknowledge the anonymous referee for their feedback and helpful suggestions.

This work is based on observations obtained at the Gemini Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), the National Research Council (Canada), CONICYT (Chile), Ministerio de Ciencia, Tecnología e Innovación Productiva (Argentina), and Ministério da Ciência, Tecnologia e Inovação (Brazil). This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research has made use of the SIMBAD and VizieR databases, operated at CDS, Strasbourg, France.

JH’s work was supported by the Space Telescope Science Institute’s Director’s Discretionary Fund and the Arizona State University Graduate Completion Fellowship. AMH’s work is supported by the National Science Foundation under Grant No. AST-2307920. ELN and FM are funded through the NASA 80NSSC17K0535.

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, modelling 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, Mazoyer & Millar-Blanchaer (2024) and Github.³ Inquiries regarding the scattering phase function from Milli et al. (2017) should be directed to Julien Milli.

Footnotes

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 Section 6.3) does not provide meaningful constraints on disc morphology: |$R_{\rm in}=74.17^{+0.08}_{-0.07}$| au, |$R_{\rm C}=70.14^{+0.22}_{-0.10}$| au, α_in > 1.34, α_out < −2.35, a_r = 0.03 ± 0.001, i = 77°.01 ± 0°.04, PA = 26°.75 ± 0°.03, dx = −2.13 ± 0.13 au, dy = 0.82 ± 0.06 au, |$\chi ^{2}_{\rm red}=1.90$|⁠.

2

Crotts et al. (2024) also identifies as HD 146 897 as an eccentric system, although we do not find strong evidence of this in our analysis.

3

https://github.com/jrhom1/DiskFMforGPIES

APPENDIX A: POSTERIOR DISTRIBUTIONS OF MODEL RUNS

Here, we present the posterior distribution functions of our model run of HD 32297 (H-band), generated with corner.py (Foreman-Mackey 2017) and shown in Fig. A1. All dashed lines in the diagonal histograms show the 16th, 50th, and 84th percentiles. Off-diagonal plots display joint probability distributions with contour levels at the same percentiles. All other posterior distribution functions are available online as supplemental material.

Figure A1.

Posterior distribution for HD 32297 H-band models.

APPENDIX B: COMPARISONS TO PREVIOUS SCATTERED-LIGHT RESOLVED MEASUREMENTS AND MODEL ESTIMATIONS

Here, we compare our constrained dust density distribution parameters with properties found in previous studies. As the best model fits in both H- and K1-bands do not model the data sufficiently, we do not discuss constrained parameters for HD 32297. Additionally, constrained parameters for the HR 4796A disc analysis using our generic SPF and the HD 114082 disc analysis using the Milli et al. (2017) SPF are not discussed as they provided poor best model fits. For many comparisons, caution must be given in comparing some dust density distribution properties such as radii from different approaches, as they are not always directly comparable. To best assess the comparisons to previous studies, we opt to compare mid-plane dust density profiles when available.

B1 HD 35841

Overall, our best-fitting model matches the data fairly well, with residual structure appearing mostly from background noise sources. Although we only find an upper limit to R_in and α_out, the remaining parameters are well constrained. From our analysis, the system appears to be a symmetric ring, well explained with a two-power law radial dust density distribution. Esposito et al. (2018) was the first study to resolve and analyse this system, characterizing the disc from GPI H-band total and polarized intensity data. Table B1 presents comparisons of parameters determined from our analysis to Esposito et al. (2018). Overall, we find good agreement in estimations of a_r within 1σ, and similar estimations for the radius of the ring. Fig. B1 compares the radial dust mid-plane density profiles of 100 randomly selected models from our posterior distribution to the median likelihood model in Esposito et al. (2018). Our analysis prefers a narrow ring with a dust density profile modelled with two power laws, while the analysis from Esposito et al. (2018) appears to prefer a somewhat broader density profile modelled with a single power law. The discrepancies between our results and Esposito et al. (2018) may be due to the differences in modelling 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.

Figure B1.

Normalized mid-plane dust density profiles from 100 randomly selected models from our posterior distribution (grey) of HD 35841 models, the maximum likelihood of our models (red), and the maximum likelihood model from Esposito et al. (2018; blue). The range along the x-axis is restricted to within our likelihood calculation mask along the disc major axis for consistent comparison. Overall, our analysis prefers a narrow ring, two-power law solution for mid-plane 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 from Esposito et al. (2018).

Table B1.

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

	This work	Esposito et al. (2018)
Instrument	GPI	GPI
Filter	H	H
Mode	Tot.	Tot & Pol
Sampler	MCMC	MCMC
Scattering	Generic	Mie
Rin	<55.97	|$59.8^{+1.1 }_{ -2.1}$|
RC	|$62.15^{+1.32 }_{ -1.61}$|	<57
αin	|$2.28^{+1.18 }_{ -0.79}$|	>−1.6
αout	<−6.09	|$-3.0^{+0.2 }_{ -0.2}$|
ar	|$0.04^{+0.005 }_{ -0.004}$|	|$0.045^{+0.023 }_{ -0.005}$|
i	|$82.69^{+0.30 }_{ -0.28}$|	84.9 ± 0.2
PA	165.49 ± 0.15	|$165.8^{+0.1}_{-0.2}$|
dx	|$1.94^{+1.11 }_{ -1.09}$|	N/A
dy	|$0.05^{+0.49 }_{ -0.48}$|	N/A

B2 HD 106906

Table B2 presents a comparison of our constrained parameters for all modelling analyses to previous studies. In our modelling analysis of HD 106906, we find shallow radial density profiles in all modelling scenarios, except for the SE side where α_in is unconstrained. Fig. B2 compares the mid-plane density profiles of 100 randomly selected models from our posterior distributions to mid-plane density profiles from previous studies. Our whole disc model profile is most consistent with modelling analysis from Crotts et al. (2021), highlighting a broad ring with shallow inner and outer slopes. Overall, the analysis from Lagrange et al. (2016) suggests a narrow mid-plane density unlike what was determined in Crotts et al. (2021), Olofsson 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 disc modelling 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 favour 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.

Figure B2.

Same as Fig. B1 but for 100 randomly selected models from our posterior distribution (grey) of whole disc HD 106906 models, the maximum likelihood of our models (red), and the best-fitting 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 mid-plane density, similar to Crotts et al. (2021).

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).

	This work	This work (NW)	This work (SE)	(1)	(2)	(3)
Instrument	GPI	GPI	GPI	SPHERE	GPI	SPHERE
Filter	H	H	H	IRDIS-H2	H	IRDIS-H
Mode	Tot. I.	Tot. I.	Tot. I.	Tot. I.	Pol. I.	Pol. I.
Sampler	MCMC	MCMC	MCMC	Grid	MCMC	Multinest
Scattering	Generic	Generic	Generic	1-HG	Mie	see [4]
Rin	<49.02	<43.24	|$55.21^{+1.48 }_{ -5.04}$|	N/A	N/A	N/A
RC	<75.88	|$67.84^{+7.95 }_{ -8.59}$|	<74.60	66.0 ± 1.8	|$72.21^{+3.10}_{-3.29}$|	87.02 ± 3.07
αin	|$3.28^{+2.23}_{-1.66}$|	<4.36	Unconst.	10 (fixed)	|$1.03^{+0.29}_{-0.24}$|	1.7 ± 0.2
αout	|$-1.97^{+0.22 }_{ -0.29}$|	|$-2.73^{+0.55 }_{ -0.67}$|	|$-2.56^{+0.27 }_{ -0.32}$|	−4.5 ± 0.3	|$-2.26^{+0.10}_{-0.08}$|	−3.7 ± 0.2
ar	|$0.03^{+0.004 }_{ -0.005}$|	0.06 ± 0.005	<0.01	0.008 (fixed)	|$0.04^{+0.0004}_{-0.0034}$|	0.047 ± 0.006
i	|$83.94^{+0.21 }_{ -0.19}$|	83.07 ± 0.22	85.21 ± 0.11	85.4 ± 0.1	|${85.57^{+0{_{.}^{\prime\prime}} 14}_{-0{_{.}^{\prime\prime}} 15}}^{a}$|	84.3 ± 0.4
PA	|$285.06^{+0.08 }_{ -0.09}$|	287.28 ± 0.14	283.32 ± 0.09	284.4 ± 0.3	|${284.31^{+0{_{.}^{\prime\prime}} 16}_{-0{_{.}^{\prime\prime}} 30}}^{a}$|	285.3 ± 0.6
dx	>4.16	<−1.92	<−1.40	N/A	|$-2.99^{+0.28}_{-0.29}$|	N/A
dy	|$-2.38^{+0.49}_{-0.52}$|	>2.23	<−3.75	N/A	|$16.48^{+2.04}_{-1.76}$|	N/A

^aThese values were not determined as part of the dust density distribution model and were found via other techniques.

Our best-fitting 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 disc could have high eccentricity and/or localized density enhancements/dearth. The best-fitting models identified from analysing the separate sides of the disc independently appear to fit the data more cleanly, but are inconsistent in i and PA, 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 disc model.

B3 HD 110058

HD 110058’s resolved debris disc 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|${_{.}^{\prime\prime}}$|3 from the star oriented in opposite directions, creating an ‘S’-like shape. Similar to measurements in Crotts et al. (2024), our modelling analyses suggest vertically broad profiles, which could be indicative of vertical stirring mechanisms such as perturbing planetary companions (see Section 7.3). Our estimates of modelling the whole disc achieve a similar result for R_in as R_C in Kasper et al. (2015) and the measurement of R_d from Crotts et al. (2024). In modelling the two sides of the disc separately, differences in PA of 2° between both analyses provide further evidence of the warps seen in images. The mid-plane density profiles from our models tend to prefer a narrow ring with a somewhat shallow outer slope. Regardless, our modelling results are tentative at best due to the loose constraints of the posterior distributions in all modelling scenarios of this system. The compact size of the disc 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. () determined that the strengths of the warps on each side of the disc are distinct, and that constraints on disc measurements are also difficult to achieve due to the compact nature of the system.

B4 HD 111520

Our best-fitting 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 disc. In modelling the two sides of the disc separately, an offset along the major axis could not be constrained in the SE model. The NW model constrained a non-zero 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 disc 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 disc aspect ratios that are similar to other systems in our analysis, and radii roughly consistent with measurements from Crotts et al. (2024). Mid-plane 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 modelling the vertical offset profile but did not incorporate disc 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 i estimates for both our whole disc model and asymmetric model approaches are consistent with the estimations from Crotts et al. (2022) to 3σ, while the PA appears only inconsistent in analysing the SE side of the disc. 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). Fig. B3 reveals that our model analyses 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 discs 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 discs, 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 disc. The compact nature and inclination of the disc may also lead to degeneracies between radial and vertical density properties.

Figure B3.

Same as Fig. B1 but for 100 randomly selected models from our posterior distribution (grey) of HD 114082 models, the maximum likelihood of our models (red), and the best-fitting 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 mid-plane dust density.

B6 HD 146897

Table B3 presents comparisons of our median likelihood model parameters to previous studies of the system. In comparing the mid-plane density profiles from our analysis to previous studies (Thalmann et al. 2013; Engler et al. 2017; Goebel et al. 2018) as in Fig. 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 mid-plane 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 disc 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 disc biasing the likelihood calculation of the MCMC, resulting in a disc model that is unable to match the brightness of the front side of the disc without overestimating the brightness of the back side of the disc. Additionally, HD 146897 could have a distinct SPF similar to HR 4796A, inducing systematic residual structure between the front and back sides of the disc, although the compact nature of the system prevents a definitive conclusion on this behaviour.

Figure B4.

Same as Fig. B1 but for 100 randomly selected models from our posterior distribution (grey) of HD 146897 models, the maximum likelihood of our models (red), and the best-fitting models from Thalmann et al. (2013; blue), Engler et al. (2017; orange), and Goebel et al. (2018; magenta). Our mid-plane 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).

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).

	This work	(1)	(2)	(3)
Instrument	GPI	HiCIAO	SPHERE	CHARIS
Filter	H	H	ZIMPOL-I	1.13 − 2.39μm
Mode	Tot. I.	Tot. I.	Pol. I.	Tot. I.
Sampler	MCMC	Grid	Grid	Grid
Scattering	Generic	1-HG	1-HG	1-HG
Rin	<49.25	N/A	N/A	N/A
RC	|$50.10^{+1.03 }_{ -6.87}$|	40 (fixed)	73 ± 8	53
αin	>0.75a	20 (fixed)	5.0 ± 1.4	6
αout	|$-2.45^{+0.12 }_{ -0.13}$|	−1.7	−2.5 ± 0.7	−1.5
ar	<0.04	N/A	0.03 ± 0.01	0.06
i	|$82.51^{+0.27 }_{ -0.26}$|	84	84.6 ± 0.85	84.6 (fixed)
PA	115.00 ± 0.17	114	114.5 ± 0.6a	114.59 ± 0.4a
dy	|$0.69^{+0.35 }_{ -0.34}$|	3	N/A	N/A

^aThese values were not determined as part of the dust density distribution model and were found via other techniques.

B7 HR 4796A

Table B4 presents comparisons of our median likelihood model parameters to previous studies of HR 4796A. Our best-fitting 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. 2020, 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 disc. Mid-plane density profiles from our analysis appear consistent to narrow ring width determined in previous studies, as seen in Fig. 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 analysed the same data sets we have used. This could suggest systematic differences between the GPI and SPHERE analyses, such as the treatment of the instrumental PSF.

Figure B5.

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

Table B4.

Comparison of our results to previous sub-arcsecond resolution studies of HR 4796A. References: (1) Milli et al. (2015), (2) Olofsson et al. (2020), (3) Olofsson et al. (2022).

	This work	This work	(1)	(2)	(3)
Instrument	GPI	GPI	NACO	SPHERE	SPHERE
Filter	H	K1	Ks	ZIMPOL-I	ZIMPOL-I
Mode	Tot. I.	Tot. I.	Pol. I.	Pol. I.	Pol. I.
Sampler	MCMC	MCMC	Grid	Multinest	Multinest
Scattering	Milli 2017	Milli 2017	1-HG	see [2]	see [2]
Rin	|$73.86^{+0.32 }_{ -0.09}$|	74.43 ± 0.10	N/A	N/A	N/A
RC	|$73.79^{+0.54 }_{ -2.11}$|	<73.07	|$75.3^{+2.0 }_{ -2.2}$|	75.44 ± 0.07	73.6 ± 1.42
αin	>2.41	>1.98	35 (fixed)	25 (fixed)	35.0 ± 0.1
αout	<−9.81	<−9.78	−10 (fixed)	−11.78 ± 0.20	−9.5 ± 0.3
ar	<0.004	<0.007	0.013 (fixed)	0.035 ± 0.001	0.041 ± 0.002
i	76.89 ± 0.03	76.95 ± 0.05	|$75.5^{+1.3 }_{ -1.7}$|	77.60 ± 0.06	77.7 ± 0.2
PA	26.75 ± 0.02	26.55 ± 0.04	26.7 ± 1.6	27.9 (fixed)	28.4 ± 0.3
dx	−2.20 ± 0.09	−2.41 ± 0.16	−5.8 ± 8.3	N/A	N/A
dy	0.84 ± 0.05	1.79 ± 0.08	|$0.8^{+3.5}_{-3.0}$|	N/A	N/A

Another consideration is the apparent brightness asymmetry between the NE and SW sides of the front-side of the disc. By introducing a stellocentric offset to the ring, pericentre glow (Wyatt et al. 1999) may be induced. Our best-fitting model appears unable to compensate for this brightness asymmetry, even after considering stellocentric offsets along both the major and minor axes of the disc. 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 pericentre may be responsible for the asymmetry.

Finally, the consistency of results and similarities in residual maps between H- and K1-band using the same measured SPF suggests that there is no difference in the SPF as a function of wavelength, despite evidence from other dust systems that more forward scattering occurs at redder wavelengths (Schröder et al. 2014).