# Gemini Planet Imager Observations of a Resolved Low-inclination Debris Disk around HD 156623

Briley Lewis, Michael Fitzgerald, Thomas Esposito, Pauline Arriaga, Ronald López, Katie Crotts, Gaspard Duchêne, Katherine Follette, Justin Hom, Paul Kalas, et al.

<figure><img src="image_1.png" /></figure>

To cite this version:

Briley Lewis, Michael Fitzgerald, Thomas Esposito, Pauline Arriaga, Ronald López, et al. Gemini Planet Imager Observations of a Resolved Low-inclination Debris Disk around HD 156623. The Astronomical Journal, 2024, 168 (4), pp.142. 10.3847/1538-3881/ad67e1. hal-04749968

HAL Id: hal-04749968

https://hal.science/hal-04749968v1

Submitted on 24 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
OPEN ACCESS

<figure><img src="image_3.png" /></figure>

CrossMark

# Gemini Planet Imager Observations of a Resolved Low-inclination Debris Disk around HD 156623

Briley L. Lewis¹, Michael P. Fitzgerald¹, Thomas M. Esposito²,³, Pauline Arriaga¹,⁴, Ronald López¹,⁵, Katie A. Crotts⁶, Gaspard Duchêne²,⁷, Katherine B. Follette⁸, Justin Hom⁹, Paul Kalas³,¹⁰,¹¹

Brenda C. Matthews⁶,¹², Maxwell Millar-Blanchaer¹³, David J. Wilner¹⁴, Johan Mazoyer¹⁵, and Bruce Macintosh¹⁶,¹⁷

¹Department of Physics and Astronomy, UCLA, Los Angeles, CA 90024, USA; blewis@astro.ucla.edu

²Department of Astronomy, University of California Berkeley, Berkeley, CA 94720, USA

³SETI Institute, Carl Sagan Center, 339 Bernardo Ave., Suite 200, Mountain View, CA 94043, USA

⁴Arete, 9301 Corbin Ave. STE 2000 Northridge, CA 91324, USA

⁵Department of Physics, UCSB, Santa Barbara, CA 93106, USA

⁶Physics & Astronomy Department, University of Victoria, 3800 Finnerty Rd. Victoria, BC, V8P 5C2, Canada

⁷Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France

⁸Department of Physics and Astronomy, Amherst College, Amherst, MA 01002, USA

⁹Steward Observatory and Department of Astronomy, University of Arizona, Tucson, AZ 85721, USA

¹⁰Astronomy Department, University of California, Berkeley, CA 94720, USA

¹¹Institute of Astrophysics, FORTH, GR-71110 Heraklion, Greece

¹²Herzberg Astronomy and Astrophysics, National Research Council of Canada, 5071 West Saanich Rd., Victoria, BC V9E 2E7, Canada

¹³Department of Physics, University of California Santa Barbara, Santa Barbara, CA 93106, USA

¹⁴Center for Astrophysics, Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138-1516, USA

¹⁵LESIA, Observatoire de Paris, Université PSL, Sorbonne Université, Université Paris Cité, CNRS, 5 place Jules Janssen, 92195 Meudon, France

¹⁶University of California Santa Cruz, 1156 High St., Santa Cruz, CA 95064, USA

¹⁷Kavli Institute for Particle Astrophysics and Cosmology, Department of Physics, Stanford University, Stanford, CA 94305, USA

Received 2024 April 18; revised 2024 July 10; accepted 2024 July 22; published 2024 September 4

## Abstract

The 16 Myr old AOV star HD 156623 in the Scorpius-Centaurus association hosts a high-fractional-luminosity debris disk, recently resolved in scattered light for the first time by the Gemini Planet Imager (GPI) in polarized intensity only. We present a new analysis of the GPI H-band polarimetric detection of the HD 156623 debris disk, with particular interest in its unique morphology. This debris disk lacks a visible inner clearing, unlike the majority of low-inclination disks in the GPI sample and in Sco-Cen, and it is known to contain CO gas, positioning it as a candidate “hybrid” or “shielded” disk. We use radiative transfer models to constrain the geometric parameters of the disk based on scattered light data and thermal models to constrain the unresolved inner radius based on the system’s spectral energy distribution (SED). We also compute a measurement of the polarized scattering phase function, adding to the existing sample of empirical phase function measurements. We find that HD 156623’s debris disk inner radius is constrained to less than 26.6 au from scattered light imagery and less than 13.4 au from SED modeling at a 99.7% confidence interval and suggest that gas drag may play a role in retaining sub-blowout-size dust grains so close to the star.

*Unified Astronomy Thesaurus concepts: Exoplanet formation (492); Planetary system formation (1257); Planet formation (1241); Circumstellar disks (235); Debris disks (363); High contrast techniques (2369); Direct imaging (387); Coronagraphic imaging (313); Polarimetry (1278)*

## 1. Introduction

High-contrast imaging is extremely useful for spatially resolving debris disks, the extrasolar analogs of our solar system’s interplanetary dust/Kuiper Belt (Mann et al. 2006). Debris disks represent a later stage of a system’s formation and evolution, wherein the circumstellar disk has evolved and lost most of its gas content. They also have low optical depth, and their dust is thought to be replenished by collisions and/or sublimation of planetesimals (Wyatt 2008). Imaging of debris disks provides significant information about the system beyond what an infrared excess measurement alone can provide, such as grain properties and disk morphology (Hughes et al. 2018). Longer-wavelength observations (e.g., those from ALMA, Atacama Large Millimeter/submillimeter Array) are able to

resolve the thermal emission of disks, providing information on the distribution of large grains and the presence of gas in a system. In the near-infrared and optical, however, we are able to trace smaller grains, and polarimetric scattered-light observations can reveal additional information about a disk’s composition and structure. Morphological information can even reveal signs of a planetary companion, as demonstrated by the disk warp that led to the discovery of β Pictoris b (Roques et al. 1994; Lagrange et al. 2010).

The Gemini Planet Imager (GPI; Macintosh et al. 2008, 2014), a ground-based adaptive-optics instrument dedicated to high-contrast discovery and characterization of exoplanets (formerly at Gemini South), has also imaged multiple debris disks with its integral-field spectrometer and dual-channel polarimeter (Esposito et al. 2020). In addition to the planet search component of the GPI Exoplanet Survey, whose results are reported in Nielsen et al. (2019), Esposito et al. (2020) report on a survey of debris disks with GPI. This survey resolved 26 debris disks and 3 protoplanetary/transitionary

<figure><img src="image_4.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
disks, several of which were resolved in scattered light for the first time.

One of these newly resolved disks was HD 156623, a debris disk around a 16 Myr old A0V star in Sco-Cen, at a distance of ~112 pc (Houk & Fuentes-Williams 1982; Gaia Collaboration 2018). HD 156623 is a low-inclination disk (≈35°), detected by GPI only in polarized light (Esposito et al. 2020). This system lacks a visible inner clearing unlike many other imaged debris disks, including others in Sco-Cen (Bonnefoy et al. 2021); close to the star, dust is usually cleared by the Poynting-Robertson drag, accretion, stellar radiation pressure, and/or blowout from stellar winds on short timescales (Wyatt 2008; Hughes et al. 2018), and a planet, if present, can sculpt a sharp inner edge of a debris disk (e.g., Engler et al. 2020; Pearce et al. 2024). Additionally, prior observations in thermal emission via ALMA imaging from Lieman-Sifry et al. (2016) reported evidence for an inner clearing and showed significant CO gas mass, leading to the “hybrid” disk label given to this system. There is still significant debate on the mechanism that leads to the observed gas mass; either primordial gas has been retained (i.e., due to shielding by carbon atoms from photodissociation of CO; Kral et al. 2019; Moór et al. 2019; Marino et al. 2020) or secondary gas has been regenerated via collisions, sublimation, or another process (Beust et al. 1990; Moór et al. 2017; Péricaud et al. 2017; Hales et al. 2019).

In this work, we aim to perform more detailed modeling of the disk around HD 156623 as seen in polarized infrared scattered light by GPI, with particular interest in constraining the presence and location of the hidden inner edge of the debris disk. In Sections 2 and 3, respectively, we include information on the host star HD 156623 and observations done by the GPI with relevant information on the data reduction. In Section 4, we present scattered-light imagery models from MCFOST and the resulting best-fit parameters for the disk geometry (Pinte et al. 2006, 2009; Pinte 2022). In Section 5, we present empirical measurements of radial brightness profiles, scattering phase function (SPF), and polarization fraction. We then present analysis of the system’s spectral energy distribution (SED) in Section 6. These measurements and insights into the inner edge location and minimum grain size in the context of this gas-rich disk are then discussed in Section 7.

## 2. Previous Observations of HD 156623

HD 156623 (HIP 84881) is an A0V star with a parallax of 9.2 ± 0.03 mas (Houk & Fuentes-Williams 1982; Gaia Collaboration 2018; Brown et al. 2021). Stellar parameters are summarized in Table 1. It exhibits pulsations consistent with zero-age main-sequence δ Scutis stars, and it is a member of the Upper Centaurus-Lupus subgroup of the Sco-Cen association, which contains stars of a wide range of ages from ~3 to 19 Myr (Pecaut & Mamajek 2016; Mellon et al. 2019; Ratzenböck et al. 2023). This star system appears to be relatively young—16 ± 7 Myr according to Mellon et al. (2019).

HD 156623 has been covered by many ground-based photometric surveys and space-based missions—including Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010; Cutri et al. 2012), 2MASS (Skrutskie et al. 2006), IRAS (Neugebauer et al. 1984; Beichman et al. 1988), and Tycho (Høg et al. 2000; Pickles & Depagne 2010)—which is useful for constructing an SED as described in Section 6. Data

<table><caption>**Table 1**
Summary of Parameters for HD 156623</caption><thead><tr><th>Parameter</th><th>Value</th><th>Reference</th></tr></thead><tbody><tr><td><i>D</i> (pc)</td><td>111.75 ± 0.96</td><td>Gaia et al. (2018)</td></tr><tr><td>Age (Myr)</td><td>16 ± 7</td><td>Mellon et al. (2019)</td></tr><tr><td><i>T</i><sub>eff</sub> (K)</td><td>9040<sup>+240</sup><sub>-160</sub></td><td>Mellon et al. (2019)</td></tr><tr><td><i>M</i><sub>*</sub> (<i>M</i><sub>☉</sub>)</td><td>1.90<sup>+0.05</sup><sub>-0.04</sub></td><td>Esposito et al. (2020)</td></tr><tr><td><i>L</i><sub>*</sub> (<i>L</i><sub>☉</sub>)</td><td>13.06 ± 1.80</td><td>Esposito et al. (2020)</td></tr><tr><td><i>R</i><sub>*</sub> (<i>R</i><sub>☉</sub>)</td><td>1.51 ± 0.09</td><td>Mellon et al. (2019)</td></tr><tr><td><i>I</i> (mag)</td><td>7.1</td><td>Esposito et al. (2020)</td></tr><tr><td><i>H</i> (mag)</td><td>7.0</td><td>Esposito et al. (2020)</td></tr><tr><td><i>L</i><sub>IR</sub>/<i>L</i><sub>*</sub></td><td>4.3 · 10<sup>-3</sup></td><td>Rizzuto et al. (2012)</td></tr></tbody></table>

**Note.** Magnitudes are synthetic photometry in Cousins *I* (Bessell et al. 1998) and 2MASS *H* (Skrutskie et al. 2006) as calculated by Esposito et al. (2020).

from WISE indicate an infrared excess for HD 156623 of $L/L_* = 4.33 \times 10^{-3}$ (Rizzuto et al. 2012).

The disk around HD 156623 was also previously observed with ALMA, which marginally resolved the disk and found a strong $^{12}$CO(2–1) signature, indicating the presence of a mass of CO gas $\ge 3.9 \cdot 10^{-4} M_{\oplus}$ (Lieman-Sifry et al. 2016). Lieman-Sifry et al. (2016) also report a 1240 µm continuum flux for HD 156623 of $720 \pm 110$ µJy. Further detections of CO gas with APEX and IRAM were also completed in Péricaud et al. (2017).

## 3. Observations and Data Reduction

### 3.1. GPI Observations

GPI observations resolved the disk in scattered light using GPI’s polarimetric mode, which uses a half-wave plate for modulation and a Wollaston prism beam splitter for analysis (Perrin et al. 2010, 2014). The polarimetric data presented herein were first published as part of a larger survey in Esposito et al. (2020). In total intensity, GPI uses angular differential imaging (ADI) to create angular diversity that can be exploited to improve the images’ final contrast (Marois et al. 2006). Observations in polarimetric mode result in measurements of the Stokes *I*, *Q*, and *U* vectors. Most data products for debris disks—including those used in this work—use azimuthal Stokes vectors, $Q_{\phi}$ and $U_{\phi}$ (described below in Equations (1) and (2), where $\phi$ is the azimuthal angle around the star measured counterclockwise starting from the north), which is useful for understanding polarization relative to the central star (Schmid et al. 2006; Millar-Blanchaer et al. 2015):

$$Q_{\phi} = Q \cos 2\phi + U \sin 2\phi \quad (1)$$

$$U_{\phi} = -Q \sin 2\phi + U \cos 2\phi. \quad (2)$$

$Q_{\phi}$ then contains all polarized light from the debris disk that is aligned perpendicular or parallel to a vector from the host star to a given pixel. $U_{\phi}$ then contains polarized light oriented 45° to this vector, which we expect to be devoid of signal for an optically thin debris disk (discussed further in Esposito et al. 2020). GPI images also include four “satellite” spots, reference images of the occulted star used for astrometric and photometric calibrations (Sivaramakrishnan & Oppenheimer 2006; Esposito et al. 2020).
<table><caption>**Table 2**
GPI Observations of HD 156623</caption><thead><tr><th>Date (UT)</th><th>Mode</th><th>t<sub>exp</sub> (s)</th><th>t<sub>int</sub> (s)</th><th>ΔPA (deg)</th></tr></thead><tbody><tr><td>2019 Apr 27</td><td>Pol</td><td>88.74</td><td>2129.81</td><td>28.2</td></tr><tr><td>2017 Apr 21</td><td>Pol</td><td>59.65</td><td>954.34</td><td>11.2</td></tr></tbody></table>

**Note.** The debris disk was detected in both polarimetric observations, and for the purposes of this investigation, we focus on the higher-quality, longer-exposure data from 2019. The following observation parameters are provided: exposure time ($t_{exp}$), total integration time ($t_{int}$), and field rotation (ΔPA, degrees). All observations were taken in the H band.

The effective outer radius of the full circular field of view in reduced images is 1"4 from the star, with radial separations up to 1"8 visible over limited ranges of position angle (PA). The field-of-view size of each frame is 2"6 x 2"6, and we use the GPI pixel scale value of $14.166 \pm 0.007$ mas lenslet$^{-1}$ (Konopacky et al. 2014; De Rosa et al. 2015, 2020b). The minimum projected separation viewable in the data is 0"123, set by the edge of the focal-plane mask (FPM; De Rosa et al. 2015).

We observed HD 156623 on 2017 April 21 and 2019 April 27. A summary of the observations are listed in Table 2. The debris disk around HD 156623 was detected by GPI in H-band ($\lambda_c = 1.65$ µm) polarized light ($Q_\phi$), in both 2017 and 2019. In 2017, sixteen 60 s exposures were taken in H-band polarimetry mode with the wave plate rotating between exposures (typically 22°.5 per move as described by Perrin et al. 2015), with 11°.2 of field rotation and approximately 1"15 seeing. In 2019, twenty-four 88 s exposures were taken in H-band polarimetry mode with 28°.2 of field rotation and ~0"39 seeing. The 2019 data was higher quality and is what will be used in our analysis.

## 3.2. Data Reduction and Polarized Differential Imaging

Data for HD 156623 were reduced using the standard GPI data reduction pipeline (DRP), documented in Perrin et al. (2014, 2016) and Wang et al. (2018). The data reduction procedure for polarimetry data in particular, as well as sources of noise in GPI polarimetric data, are further described in Perrin et al. (2015) and Section 4.1 of Esposito et al. (2020); this reduction, however, is an improved reduction from that presented in Esposito et al. (2020) due to the quadropole correction described later in this section. The GPI DRP completes dark subtraction, flexure correction, field distortion correction, destriping, bad pixel corrections, and other calibrations specific to the GPI instrument, such as a cleaning procedure to account for biasing between orthogonal polarization channels (Konopacky et al. 2014; Perrin et al. 2014, 2015; Wiktorowicz et al. 2014). After the standard DRP, the series of reduced polarimetric observations is converted into a single data cube containing the Stokes vectors $I$, $Q$, $U$, and $V$, before then converting $Q$ and $U$ to radial Stokes vectors ($Q_\phi$ and $U_\phi$) as previously described (Perrin et al. 2015). This DRP produces images in units of analog-to-digital units per coadd (ADU coadd$^{-1}$); using the measured flux of the satellite spots, which are themselves images of the host star and the star's known H-band magnitude (see Table 1), we converted those detector units into surface brightness units of millijanskys per square arcsecond for the subsequent analysis, as in Wang et al. (2014), Hung et al. (2016), and Esposito et al. (2020). The calibration

factor for this data reduction was $6.41 \times 10^{-9} \pm 0.27 \times 10^{-9}$ Jy/(ADU/coadd). Additionally, an extra instrumental polarization subtraction—a quadrupole modulation of the total intensity pattern scaled to data and subtracted for each wave-plate modulation—was applied to the data to remove an instrumental polarization effect seen in GPI images (De Rosa et al. 2020a).

Polarized intensity frames from the standard data reduction readily show the disk due to the natural separation of unpolarized stellar light and polarized light scattered from the disk. For total intensity, we must complete stellar point-spread function (PSF) subtraction, for which we use the `pyKLIP` package (Wang et al. 2015). `pyKLIP` uses Karhunen–Loève Image Processing (KLIP), a form of principal components analysis, to create a model of the stellar PSF using its dominant modes of variation for PSF subtraction (Soummer et al. 2012). For polarimetric GPI data, the total intensity is the sum of the two orthogonal polarization states, plus any unpolarized light. We employ ADI in our reduction for total intensity data taken with GPI polarimetry mode. Parameters for the `pyKLIP` reduction generally must be chosen with care to maximize the signal from the disk, avoiding self-subtraction (Milli et al. 2012) while also reducing contamination from the stellar PSF. For this particular disk, which is faint relative to the stellar halo in total intensity, we use a particularly unaggressive implementation of KLIP, using one annulus and one subsection in the subtraction as recommended by `pyKLIP` documentation. We additionally test multiple Karhunen–Loève modes (1, 5, 10, 20, 50) and use `pyKLIP`'s minimum rotation ($\minrot=3$) criterion for basis vector selection, as is recommended to avoid self-subtraction. Even after applying PSF subtraction, however, we were unable to recover the disk in total intensity light, as shown in Figure 1. Although other recent algorithms have focused on improving this issue of self-subtraction (e.g., nonnegative matrix factorization; Ren et al. 2018), we do not pursue additional reductions in this work; although such algorithms may be able to recover the disk's general morphology, we choose to focus on the polarized intensity detection that already provides such morphological information. Additionally, low-inclination and radially extended disks are particularly challenging to retrieve in total intensity with an ADI approach. The point sources in total intensity were already determined by Esposito et al. (2020) to be unbound background objects and are thus not considered in our analysis of this system.

## 4. Disk Morphology

In this section, we reproduce the modeling originally presented in Esposito et al. (2020) to show the scattered-light model disk images and quality of fit to the data, which were not provided in that work, with the goal of constraining the inner edge from our imagery. This modeling is nearly identical to that of Esposito et al. (2020) but run with a 40% longer chain. This leads to a small improvement in the quality of the results due to more total samples of the posterior. We use these morphological results (namely, PA, inclination, and radius constraints) to guide our measurements of the disk's brightness profiles and phase function, calculated in Section 5.

As in Esposito et al. (2020), morphological parameters of the disk were estimated by comparison to models from the Monte Carlo radiative transfer code `MCFOST` (Pinte et al. 2006, 2009). `MCFOST` is a radiative transfer code that uses Monte Carlo and
<figure><img src="image_5.png" /><figcaption>**Figure 1.** PSF-subtracted total-intensity GPI H-band image of HD 156623 using five KL modes (top) and polarized intensity GPI H-band image of HD 156623 smoothed with a 2 pixel wide Gaussian kernel for visualization (bottom). The debris disk is only detected in polarized light, with a signal-to-noise ratio of ~2–3 per resolution element. The point sources visible in the total intensity image were determined to be unbound background stars in Esposito et al. (2020).</figcaption></figure>

ray-tracing methods to produce models of observables (scattered-light imagery, millimeter and infrared visibilities, line maps, SEDs, and dust temperature distributions) for circumstellar environments such as protoplanetary and debris disks. Models were fit to Stokes $Q_\phi$, $Q$, and $U$ data simultaneously, as $Q_\phi$ is the only available data with the quadrupole noise signal removed. Following the procedure in Esposito et al. (2018), our models assume Mie scattering from spherical grains in an azimuthally symmetric, optically thin disk centered on the star, observed at 1.65 µm (the H-band peak wavelength). We additionally use the same disk parameters and procedure as described in Esposito et al. (2020), as summarized below. The disk has a vertical structure where the scale height is a constant fraction of the radius, set to 0.055, a value consistent with measured values for similar disks in the literature and consistent with modeled vertical structure assuming the effects of radiation pressure and collisions (Krist et al. 2005; Thébault 2009; Millar-Blanchaer et al. 2015;

Esposito et al. 2018). The surface density is described by a smoothed broken/two-component power law, with two different fixed power-law indices: one for the disk interior to a critical radius ($R_c$), set to $\alpha_{in} = -3.5$, and one for the disk exterior to $R_c$, set to $\alpha_{out} = 1.5$ as in Esposito et al. (2020); these values were manually tuned to a reasonable by-eye match in that work, and we have adopted the same parameter choice for consistency. The following model parameters were left free to vary in the fit: the minimum grain size ($a_{min}$), grain porosity, dust mass (within the grain size distribution used) $M_d$, disk PA, inclination ($i$), dust inner radius $R_{in}$, dust outer radius $R_{out}$, and the critical radius $R_c$ at which the inner and outer surface density profiles cross. Although the $\alpha$ power-law exponents can theoretically be used to define a disk edge where the signal falls below the noise, $R_{in}$ and $R_{out}$ are included as free parameters to be able to define an edge that is more abrupt than that from a smoothly declining power law, e.g., from planet sculpting. Additionally, we do not expect the dust mass $M_d$ to be physically meaningful; instead, it essentially functions as a scaling factor for the global surface brightness of the disk. PA is defined as the angle counterclockwise from north to the projected major axis of the disk.

We used a parallel-tempered Markov Chain Monte Carlo (MCMC) with three temperatures, each with 120 walkers. Walkers were initialized randomly from a uniform distribution and then simulated for 3500 iterations. We discarded the first 2900 iterations as “burn-in,” where the ensemble of walkers had not yet converged to their peak in posterior space, and used the final 600 iterations of the zeroth temperature walkers (equivalent to $2.16 \times 10^5$ models) for the disk parameter value estimates, which are summarized in Table 3. “Max Lk Value” contains the values for the maximum-likelihood model, which is shown in Figure 2, and the adjacent column shows the median values with a $1\sigma$ confidence interval. A corner plot of posteriors from the MCMC is shown in Figure 3. Multiple parameters have slightly multimodal or non-Gaussian posteriors, and as such, the $1\sigma$ values should be interpreted with caution for those parameters; parameters with strongly non-Gaussian posteriors are reported as either maximum-likelihood values or 99.7% confidence upper/lower limits. Similar to Esposito et al. (2020), we aim for a model that constrains the necessary morphological parameters with sufficient accuracy for our purposes, while acknowledging that a more complex/physically motivated treatment of grain properties is necessary to truly reproduce the observed scattered-light distribution.

# 5. Brightness and Scattering

## 5.1. Radial Brightness Profiles

Radial brightness profiles allow for further investigation of both the inner clearing and the outer radius of the disk, e.g., determining when the observed data is consistent with zero disk signal, in addition to the above modeling with MCFOST.

To that end, we measure radial brightness profiles (Figure 5) for various PAs. To compute these profiles, we first deproject the disk based on our determined value of inclination, such that each radial cut traverses the same distance from the center of the disk, using the Python package `diskmap` (Stolker et al. 2016). Zero PA is defined counterclockwise from the north, as described earlier. Radial brightness profiles are the average of a wedge at the given angle range that is 5 pixels long in the radial direction and 45° (~5–80 pixels depending on the radial distance from the
<table><caption>**Table 3**
Summary of Retrieved Disk Parameters from MCFOST from the Final 600 Iterations (i.e., Excluding Burn-in)</caption><thead><tr><th colspan="3">Best-fit Morphological Parameters</th></tr><tr><th>Parameter</th><th>Max. Lk. Value</th><th>Value from Esposito et al. (2020)</th></tr></thead><tbody><tr><td>PA (deg)</td><td>102.081</td><td>100.9<sup>+1.6</sup><sub>-2.2</sub></td></tr><tr><td>Inclination (deg)</td><td>38.978</td><td>36.6<sup>+2.1</sup><sub>-9.9</sub></td></tr><tr><td>R<sub>in</sub> (au)</td><td>12.707</td><td>&lt; <strong>26.6</strong></td></tr><tr><td>R<sub>out</sub>(au)</td><td>184.482</td><td>&gt; <strong>139.2</strong></td></tr><tr><td>R<sub>c</sub>(au)</td><td>67.113</td><td>65.2<sup>+1.3</sup><sub>-2.3</sub></td></tr></tbody></table>

<table><thead><tr><th colspan="4">Other Parameters Used</th></tr><tr><th>Parameter</th><th>Free/ Fixed</th><th>Max Lk. or Fixed Value</th><th>Value from Esposito et al. (2020)</th></tr></thead><tbody><tr><td>a<sub>min</sub> (μm)</td><td>Free</td><td>1.13</td><td>1.19</td></tr><tr><td>M<sub>d</sub> (logM<sub>☉</sub>)</td><td>Free</td><td>-6.31</td><td>-6.31<sup>+0.21</sup><sub>-0.05</sub></td></tr><tr><td>Porosity (%)</td><td>Free</td><td>0.82</td><td>0.84</td></tr><tr><td>α<sub>in</sub></td><td>Fixed</td><td>1.5</td><td>1.5</td></tr><tr><td>α<sub>out</sub></td><td>Fixed</td><td>-3.5</td><td>-3.5</td></tr><tr><td>β</td><td>Fixed</td><td>1</td><td>1</td></tr><tr><td>H<sub>0</sub>/r<sub>H</sub></td><td>Fixed</td><td>0.055</td><td>0.055</td></tr></tbody></table>

**Note.** The top section outlines the best-fit morphological parameters adopted for our analysis, alongside the values from Esposito et al. (2020), which are generally in agreement. The bottom section describes parameters that are not necessarily physically meaningful in this application and should be interpreted with caution but are included here for reproducibility of our models.

center) wide, with its error as the variance of that slice of the wedge. Pixels within 15 au of the center are masked, as this region is within the influence of the FPM and therefore contains large errors. To determine a rough estimate for the visual outer radius of the disk—i.e., where the disk signal is dominant over the noise—in our data, we rebin the brightness profiles to 5 pixel wide annuli and determine where the signal becomes consistent with zero by $2\sigma$ in our images at each angle. We average those measurements across azimuthal angles to find a visual outer radius of ~78 au. This visual outer radius is not an inherent property of the disk but instead a quantitatively derived limit for where the image is dominated by the signal instead of the noise for use in further calculations, e.g., the SPFs in Section 5.2. We additionally compute radial surface density profiles by correcting the brightness profiles by a factor of $1/r^2$, as shown in Figure 4. As an independent check on the disk morphology, we fit a smoothly broken power law (as defined via MCFOST and below in Equation (3)) to the average measured surface density profile using least-squares, shown in Figure 4:

$$ \rho(r, z) \propto \frac{\exp[-|z|/H(r)^\gamma]}{[(r/r_c)^{-2\alpha_{in}} + (r/r_c)^{-2\alpha_{out}}]^{1/2}} \quad (3) $$

We set the z coordinate to a constant for simplicity; use the same formalism for $H(r) = H_0(r/r_H)^\beta$ as Esposito et al. (2020), where $\beta=0$ and $H_0/r_H = 0.055$; set the exponent $\gamma=1$; and leave $r_c$, $\alpha_{in}$, and $\alpha_{out}$ as free parameters, along with

an arbitrary scaling/normalization factor $N$. The best-fit values for the average radial surface density profile are $N=1.73 \pm 0.13$, $r_c = 78.8 \pm 6.0$ au, $\alpha_{in} = 0.58 \pm 0.15$, and $\alpha_{out} = -1.31 \pm 0.15$. These differ significantly from the scattered-light modeling values for $r_c$, $\alpha_{in}$, and $\alpha_{out}$. These differences may be explained simply by the fact that $\alpha_{in}$ and $\alpha_{out}$ were fixed in the MCFOST models, as $r_c$ is expected to vary with any change in the power-law exponents. Additionally, this average profile may be biased due to a lower SNR detection of the backscattering half of the disk, where the outer edge is much less well defined, or a different SPF than that assumed by MCFOST, which may change the ratio of surface brightness to surface density.

## 5.2. Observed Polarized SPFs

SPFs describe the fraction of light scattered into a solid angle by a single dust grain, and their shapes are determined by the properties of the dust itself; therefore, they are often used as a probe of the grain size, porosity, composition, and shape of a disk’s dust grains.

We make several assumptions in this calculation, similar to those used for geometric modeling; that is, we take the disk to be azimuthally symmetric (but radially varying, related to the observed surface brightness) and optically thin, with the previously described morphological parameters inferred via MCFOST. We also assume that scattering is independent of azimuthal angle (e.g., only the scattering angle $\theta$ is a variable in our equations), which is true for a spherical grain by symmetry and a good approximation for randomly oriented nonspherical grains. The expression for a normalized SPF is then as follows:

$$ f(\theta) \sim \frac{B_{\text{obs}}(r, \theta)}{\rho(r)} \cdot r^2, \quad (4) $$

where $f(\theta)$ is the SPF, solely dependent on scattering angle and radially varying properties; $B_{\text{obs}}(r, \theta)$ is the observed brightness at a given radius $r$ and PA $\theta$; $\rho(r)$ is the surface density of the disk, derived from the observed radial brightness profile; and $r$ is the radius from the central star.

Given that this disk does not have a clearly defined single ring at which to compute an SPF, SPFs were instead determined by taking measurements of $B_{\text{obs}}(r, \theta)$ at multiple radii around the disk, then averaging over the varying radii. For this specific calculation, we limit the range over which we compute the SPF to the annuli beyond the FPM and within the visual outer radius as determined by the radial brightness curves; this range is 15 au (at the FPM limit) to 78 au. Since we do assume the disk’s surface density can vary radially, we retain $\rho(r)$ and $r$ in our calculations even though they are constants for a given radius. The disk surface density, $\rho(r)$, was approximated using the average best-fit radial brightness profile where the peak surface brightness was normalized to 1.

Scattering angle ($\theta$) for each point in the disk is determined analytically, as shown in Figure 6, following the above geometric assumptions. This disk is most forward scattering around ~200° PA, with a scattering angle of 54°, and most backward scattering near ~10° PA with a scattering angle of 126°.

HD 156623’s debris disk appears slightly asymmetric in visual inspection, particularly in the northwest quadrant, so SPFs were measured separately for the east and west halves of the disk (divided approximately along the minor axis), as shown in Figure 7. Errors were derived by the variance of
<figure><img src="image_6.png" /><figcaption>**Figure 2.** MCFOST maximum-likelihood scattered-light model of HD 156623 with data and residuals for Stokes $Q_\phi$, $U_\phi$, $Q$, and $U$, drawn from only the last 600 iterations (i.e., excluding burn-in).</figcaption></figure>

pixels in each annulus, similar to the procedure for the radial brightness profiles, and propagated appropriately through the calculation. There is a slight asymmetry that is most pronounced near the eastern ansae (~80°–95°), consistent with the brightness asymmetry described in Crotts et al. (2024) that suggests the disk may be eccentric or experiencing more

frequent collisions on the eastern half. However, it is also worth noting that in our calculation, there is an implicit assumption that the density is azimuthally constant and radially varying. As a result, even if there is a true increase in brightness at those angles, there may not be a corresponding overbrightness in the SPF if we do not assume this symmetry.
<figure><img src="image_7.png" /><figcaption>Figure 3. Corner plot of posteriors from MCFOST modeling of scattered-light imaging.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 4. The average radial surface density profile (gray points with error bars) fit with a smoothly broken two-component power law (the dashed blue line). The hatched region on the left marks the limits from the FPM as in other figures, and the dotted-dashed light purple curve represents the radial profile determined by scattered-light modeling. The values derived from this average surface density profile for r_c, α_in, and α_out are substantially different from those determined via scattered-light modeling, perhaps due to a difference in the SPF or bias due to the low SNR detection of the backscattering portion of the disk; alternatively, these differences may simply be an artifact of the different modeling approaches.</figcaption></figure>

We additionally fit a Henyey-Greenstein (HG) SPF (Henyey & Greenstein 1941) to each of our observed scattering SPFs divided by a Rayleigh polarization curve (approximating the total intensity SPF), also shown in Figure 7. The HG function p (θ) is an empirical description, not based in Mie theory or another physical scenario, and it is calculated as

$$p(\theta) = \frac{1}{4\pi} \frac{1 - g^2}{[1 + g^2 - 2g \cos \theta]^{3/2}}, \quad (5)$$

where g is the HG asymmetry parameter (g = cos θ, ranging from -1 to 1) and θ is the scattering angle. For debris disks, which are typically forward scattering, g is greater than 0. For the eastern half of the disk, the HG function is a poor fit, even with the best-fit value g = 0.33; however, the HG function with g = 0.45 is a reasonable approximation for the western half. These g values are broadly consistent with similar measurements for other low-inclination debris disks as summarized by Hughes et al. (2018), supporting the ideas that disk dust may not obey one HG SPF at all scattering angles and that best-fit g values may correlate with the range of scattering angles probed. Although the HG function is not inherently physically informative, we do not pursue physical dust model fits to our
<figure><img src="image_9.png" /><figcaption>**Figure 5.** (Left) Measured brightness curves along the radial direction of HD 156623's deprojected disk for various angles. The region covered by the FPM is blocked out (hatched in light purple, labeled FPM), and it is clear that in most cases, the inner hole is unresolved, with the exception of possibly PA 45°–90° and PA 90°–235°. Each profile is arbitrarily offset for clarity, and the dashed lines correspond to zero for each profile depicted. (Right) Surface density (brightness corrected by a factor of 1/r²) profiles along the radial direction of HD 156623's deprojected disk for various angles. The region covered by the FPM is blocked out (hatched in light purple, labeled FPM). Each profile is arbitrarily offset for clarity, and the dashed lines correspond to zero for each profile depicted.</figcaption></figure>

<figure><img src="image_10.png" /></figure>

phase function given the limited range of scattering angles for the available data. Additionally, there is still an asymmetry in the derived total intensity SPF—albeit less pronounced—again mirroring findings from Crotts et al. (2024).

## 6. Thermal Emission

The main goal of our investigation of HD 156623's SED is to constrain the location of the disk's unseen inner edge. Photometric measurements of HD 156623 were acquired from a variety of sources, ranging from visible to mid-infrared to millimeter ALMA observations. Most measurements were obtained in magnitudes from the IRSA or Vizier databases and then converted to Janskys using the zero-points provided by each mission/project's documentation. These observations and their references are summarized in Table 4; all data listed here are before color corrections as we instead treat points with large bandpasses (namely W3/W4 and AKARI S9W) using synthetic photometry as described below. Although photometric data from IRAS were available for this object, they were not used due to poor-data-quality flags present in the catalog. Given the wide bandpasses of the WISE W3/W4 (Wright et al. 2010) and AKARI S9W (Fujiwara et al. 2010) filters, we compute synthetic photometry by integrating the model flux multiplied by the transmission profile for each filter (from Rodrigo & Solano 2020) over its bandpass via the following

<figure><img src="image_11.png" /><figcaption>**Figure 6.** Diagram with the scattering angle across the disk of HD 156623, assuming the geometric parameters defined in Section 4. The most forward scattering occurs in the southwest, and the most backscattering occurs in the northeast.</figcaption></figure>

equation as in Tokunaga & Vacca (2005):

$$ \langle F_{\nu} \rangle = \frac{\int F_{\nu}(\nu) S(\nu) / \nu \, d\nu}{\int S(\nu) / \nu \, d\nu}, \qquad (6) $$

where $\langle F_{\nu} \rangle$ is the monochromatic flux density (ergs per second per square centimeter per Hertz) for our synthetic photometry,
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 7.** (Left) Measured polarized light SPF for the disk, split between the eastern and western halves of the disk as defined by the approximate minor axis. We observe an asymmetry between the eastern and western halves, consistent with recent findings in Crotts et al. (2024) suggesting this disk may be slightly eccentric. (Right) The measured east/west SPFs divided by a Rayleigh polarization curve, along with HG functions fit to each SPF. There is a lesser but still distinct asymmetry between the two halves of the disk, and the HG function is a poor fit for the eastern half.</figcaption></figure>

<table><caption>**Table 4**
Summary of Photometry Used in Creating the Spectral Energy Distribution for HD 156623</caption><thead><tr><th>Band</th><th>Source</th><th>λ<sub>eff</sub><br>(μm)</th><th>Flux<br>(Jy)</th><th>References</th></tr></thead><tbody><tr><td><i>B</i></td><td>Tycho 2</td><td>0.44</td><td>4.75 ± 0.10</td><td>Hog et al. (2000);<br>Cox (2015)</td></tr><tr><td><i>V</i></td><td>Tycho 2</td><td>0.55</td><td>4.71 ± 0.1</td><td>Hog et al. (2000);<br>Cox (2015)</td></tr><tr><td><i>J</i></td><td>2MASS</td><td>1.24</td><td>2.33 ± 0.05</td><td>Cohen et al. (2003);<br>Skrutskie et al.<br>(2006)</td></tr><tr><td><i>H</i></td><td>2MASS</td><td>1.66</td><td>1.47 ± 0.03</td><td>Cohen et al. (2003);<br>Skrutskie et al.<br>(2006)</td></tr><tr><td><i>K</i><sub>s</sub></td><td>2MASS</td><td>2.16</td><td>1.05 ± 0.02</td><td>Cohen et al. (2003);<br>Skrutskie et al.<br>(2006)</td></tr><tr><td><i>W1</i></td><td>ALLWISE</td><td>3.35</td><td>0.502 ± 0.027</td><td>Wright et al. (2010);<br>Mainzer et al.<br>(2011)</td></tr><tr><td><i>W2</i></td><td>ALLWISE</td><td>4.6</td><td>0.30 ± 0.010</td><td>Wright et al. (2010);<br>Mainzer et al.<br>(2011)</td></tr><tr><td><i>S9W</i></td><td>AKARI</td><td>9</td><td>0.215 ± 0.015</td><td>Murakami et al.<br>(2007); Ishihara<br>et al. (2010)</td></tr><tr><td><i>W3</i></td><td>ALLWISE</td><td>11.6</td><td>0.155 ± 0.007</td><td>Wright et al. (2010);<br>Mainzer et al.<br>(2011)</td></tr><tr><td><i>W4</i></td><td>ALLWISE</td><td>22.1</td><td>0.626 ± 0.036</td><td>Wright et al. (2010);<br>Mainzer et al.<br>(2011)</td></tr><tr><td>Band 6</td><td>ALMA</td><td>1240</td><td>0.00072 ± 0.00011</td><td>Lieman-Sifry et al.<br>(2016)</td></tr></tbody></table>

$F_ν$ is the monochromatic flux density (ergs per second per square centimeter per Hertz) from the object, $S(ν)$ is the filter transmission function, and $ν$ represents the frequencies in the

bandpass. Additional absolute calibration uncertainties were incorporated for the WISE data points (Cutri et al. 2012) and AKARI data points (Tanabé et al. 2008).

We initially fit a simple model consisting of two single-temperature blackbody curves for the stellar emission plus the disk's thermal emission to these monochromatic flux densities using a Levenberg-Marquardt least-squares fit implemented in `scipy`. The free parameters included are disk temperature ($T_d$) and a parameter proportional to the disk-emitting area ($A$). The stellar radius ($R_*$) and temperature ($T_*$) are set using the values in Table 1 from Mellon et al. (2019). The best-fit dust temperature is $184 \pm 23$ K with a reduced chi-squared $\chi^2$ value of $1.088$, corresponding to a radius of $8.49^{+2.60}_{-1.78}$ au if treated as an equilibrium temperature for a grain radiating as a blackbody. It is additionally worth noting that radii derived from SEDs are known to be underestimates, i.e., they are often smaller than the inner radii derived from resolved observations (Wyatt 2008; Booth et al. 2013; Pawellek & Krivov 2015; Hom et al. 2020; Pawellek et al. 2021).

Given the poor fit at long wavelengths while using a single-temperature blackbody to represent thermal emission from dust, we subsequently tested a model for the SED using a single-temperature blackbody plus an additional modified wavelength-dependent emissivity, following the formalism in Backman et al. (1992). Again, this model does not assume a spatial distribution or grain size distribution. In this model, grains of a characteristic size ~$\lambda_0/2\pi$ radiate with a modified emissivity $\epsilon_ν$ as follows:

$$ \epsilon_{\nu} = \begin{cases} \frac{\lambda_0}{\lambda} & \lambda > \lambda_0 \\ 1 & \text{otherwise} \end{cases}. \qquad (7) $$

This model was again fit using the `scipy` implementation of a Levenberg-Marquardt least-squares fit, now for both temperature and $\lambda_0$. The best-fit disk temperature for this model is $140 \pm 23$ K with a $\chi^2$ value of 1.21, corresponding to a radius of $14.7^{+6.34}_{-3.85}$ au, with $\lambda_0$ unconstrained (best-fit value of 219 µm, corresponding to a characteristic grain size of ~35 µm). These models are visualized in Figure 8. Neither model is able to capture the 9 µm AKARI point, and a modified emissivity is necessary to reproduce the flux observed by ALMA.
<figure><img src="image_14.png" /><figcaption>**Figure 8.** Single-component blackbody models using best-fit values to available photometry for HD 156623, with fractional residuals shown below. The dust's blackbody contribution to the SED is plotted in red. (Left) For a simple single-component blackbody SED, the photometry is best fit with temperature $T = 184 \pm 23$ K. (Right) For a modified-emissivity single-component blackbody SED, the photometry is best fit with $T = 140 \pm 23$ K and $\lambda_0 = 219$ µm. Neither model is able to capture the 9 µm AKARI point, and a modified emissivity is necessary to reproduce the flux observed by ALMA.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 9.** Maximum-likelihood SED model from MCFOST (as a dashed black line, with synthetic photometric points marked) compared to available photometry (colored according to its source). The best-fit model requires sub-blowout size grains and a close-in inner radius around ~3 au. Five hundred additional SED models randomly drawn from the post-burn-in MCMC are plotted in gray.</figcaption></figure>

However, these simple models do not explicitly fit for the inner radius we are interested in as they assume a single temperature across the entire disk. To constrain the inner radius, we then compute the SED for a spatially resolved, azimuthally symmetric disk using MCFOST. We set up this resolved model using a Kurucz 9000 K model atmosphere (Castelli & Kurucz 2004) and the same parameters and assumptions (e.g., Mie scattering, dust grain size and composition) as the scattered-light model in Section 4, leaving only $r_{in}$, $a_{min}$, and $M_d$ as free parameters. We explore the posterior for those parameters again with a Goodman & Weare (2010) affine invariant MCMC ensemble sampler using the emcee package with 16 walkers and uniform priors, with bounds on the priors for $r_{in}$ ($0 < r_{in} < 15$) and $\log a_{min}$ ($-2 < \log a_{min} < 1$). Walkers were initialized randomly from a uniform distribution and then simulated for 1500 iterations. We discarded the first 100 iterations as “burn-in” where the ensemble of walkers had not yet converged to their peak in

<figure><img src="image_16.png" /><figcaption>**Figure 10.** Corner plot for MCMC fit with photometry and MCFOST SED models. The reported value for $r_{in}$ is an upper bound ($<13.43$ au at 99.7% confidence), as seen in the posterior distribution. Additionally, the value for $\log a_{min}$ is also only an upper limit of $-0.685$ (0.21 µm); this suggests that different assumptions for the grain composition and/or size distribution may be needed to fully understand the minimum grain size of this system.</figcaption></figure>

posterior space and used the remaining 1400 iterations (equivalent to $2.2 \times 10^4$ models) for the disk parameter value estimates. The maximum-likelihood model SED is shown compared to available photometry in Figure 9, and a corner plot of the approximate posteriors is shown in Figure 10.

Scattered-light image modeling is limited by the focal-plane mask in its determination of $r_{in}$, resulting in a maximum-likelihood value that corresponds to the FPM edge (~12–15 au) and a $3\sigma$ upper limit of $<26.6$ au; SED modeling, on the other hand, returns a maximum-likelihood value of 2.37 au and is limited by the prior bound at zero, resulting in an
upper limit of <13.4 at 99.7% confidence (<8.85 au at 95% confidence). The value for $log(M_d)$ should be considered as a scaling factor more than a physically relevant value as we assume a sparse grain density and do not consider dynamical effects like collisions; for the sake of reproducing our modeling, the value for $log(M_d)$ for the maximum-likelihood model shown in Figure 9 is -7.02. The minimum grain size $a_{min}$ is much smaller as suggested by SED modeling, <0.21 µm as compared to ~1 µm from scattered-light image modeling. For comparison, the blowout size for the system is around 3.5 µm, as further discussed in Section 7. This $a_{min}$ value from SED modeling reached a bound on the prior at 0.01 µm, suggesting that the Dohnanyi size distribution or other grain characteristics we have assumed may not fully describe this system; however, given the limited number of photometric data points currently available (four data points—AKARI 9 µm, W3, and W4 and ALMA—and currently three free parameters), we choose not to add further free parameters to our model in this work. This best-fit MCFOST SED model has $\chi^2 = 1.11$, similar to the simpler models above. Again, the AKARI and ALMA points are not fit well by this model; we discuss this discrepancy—as well as the discrepancy in $a_{min}$—further in Section 7. The maximum temperatures calculated by MCFOST at $r_{in}$ and $R_c$ at the disk midplane in the best-fit model are ~500 K and 355 K respectively.

## 7. Discussion

All SED fits are unable to capture the AKARI data point (the red point in Figure 8), and no single model we have tested is able to satisfactorily explain the entire set of photometric data points available. MCFOST modeling similarly has trouble reproducing the AKARI point, as well as the long wavelength flux measured by ALMA, indicating that different assumptions for grain composition and/or size distribution may be needed; given the limited number of photometric points available, exploring those as free parameters in a model is beyond the scope of this work. However, based on our modeling, it appears that—under the assumptions described in our work—the thermal signature of the disk is best explained by a population of close-in warm dust with sub-blowout size grains and a modified emissivity trend at longer wavelengths. One could additionally imagine that the debris disk has two components: a cold outer disk imaged by GPI and ALMA and a separate warm inner belt that is behind the GPI FPM and unresolved by ALMA but also boosts the 9 µm emission as seen by AKARI. Future work may include exploring two-component models, and further observations (e.g., JWST/MIRI spectroscopy) are necessary to better define the shape of the spectrum through the mid-infrared.

The radiative transfer modeling of the scattered-light data is limited by the focal-plane mask, initially constraining $r_{in}$ to be less than 26.6 au. Single-component blackbody fits to the SED suggest an inner edge of ~8–15 au (not accounting for nonblackbody behavior of real dust grains or the possibility of grain overheating; e.g., Morales et al. 2011; Pawellek & Krivov 2015) and disk models with MCFOST fit to the SED constrain the edge to <13.4 au. However, it is again worth noting that even the best-fit SED model is unable to reproduce both the 9 µm excess and the long wavelength flux observed by ALMA, indicating that different assumptions on the grain properties or surface density distribution (e.g., inner and outer power laws) are necessary to fully explain this system; given

the lack of photometry at this time, we leave this endeavor for future work. Lieman-Sifry et al. (2016) similarly constrain the inner radius of this disk using ALMA data and SED modeling, with a best-fit value of $r_{in} = 10$ au from a visibility-only fit and $r_{in} = 2.99$ au from a simultaneous visibility and SED fit. Our result is in agreement with these values within $1\sigma$. This <13.4 au inner edge is certainly unique; HD 156623 is one of few known low-inclination debris disks in Sco-Cen imaged in scattered light (which is biased toward edge-on disks) where the inner clearing is not visible due to its location behind the FPM (Bonnefoy et al. 2021). Additionally, this modeling assumes a certain power law for surface density; different assumptions in that value may lead to significant changes in the value of $r_{in}$, possibly explaining the slight differences between our estimate and that of Lieman-Sifry et al. (2016).

The blowout grain size of this system is estimated as ~3.5 µm, according to Arnold et al. (2019), which uses the equation below (Burns et al. 1979) to compute the blowout size for HR 4796, a similar temperature A-type star:

$$a_{\text{BO}} = 1.152 \langle Q_{\text{pr}} \rangle \left( \frac{L_{*}}{L_{\odot}} \right) \left( \frac{M_{*}}{M_{\odot}} \right)^{-1} (1 - \mathcal{P}) \left( \frac{\rho}{g \text{ cm}^{-3}} \right)^{-1}, \quad (8)$$

where $L_*$ is the stellar luminosity in solar luminosities, $M_*$ is the stellar mass in solar masses, $\mathcal{P}$ is the grain porosity (assumed to be zero), $\rho$ is the grain density (assumed to be $3.3 \text{ g cm}^{-3}$ from Draine & Lee 1984; Draine 2003), and $\langle Q_{\text{pr}} \rangle$ is the average radiation pressure efficiency, calculated from grain-scattering properties through

$$\langle Q_{\text{pr}} \rangle = \frac{1}{n} \sum_{i=\nu_{\text{min}}}^{\nu_{\text{max}}} [F_{\nu_{*},*} Q_{\text{abs},i} + (1-g) Q_{\text{sc},i}], \quad (9)$$

where $\nu_{min}$ is the lowest frequency used in the calculation, $\nu_{max}$ is the highest frequency used in the calculation, $F_{\nu_{*},*}$ is the stellar spectrum at a given frequency, $Q_{\text{abs},i}$ is the absorption coefficient for a given frequency, $g$ is the asymmetry parameter, $Q_{\text{sc},i}$ is the scattering coefficient for a given frequency, and $n$ is the number of frequency samples used. We will adopt the blowout size of 3.5 µm for the remainder of this discussion, which assumes astrosilicate grains as we have in our modeling.

Values for the minimum grain size from both scattered-light image modeling and SED modeling are in significant disagreement with the calculated blowout size; additionally, these values constrain the minimum grain size from different data and different assumptions on the physics of grain interactions with light as the image modeling value is derived from the polarized SPF and the SED modeling value is derived from thermal equilibrium. MCFOST scattered-light image modeling reports the maximum-likelihood value for $a_{min}$ as 1.13 µm, with smaller peaks in the probability distribution around 1.35–1.82 µm as shown in Figure 3, and SED modeling suggests <1 µm); however, these values are much smaller than the minimum dust grain size predicted by the blowout size of 3.5 µm. Lieman-Sifry et al. (2016) provide additional evidence for sub-blowout size grains via their simultaneous model fits to their observed ALMA visibilities and the SED, finding a minimum grain size around ~1 µm, in agreement with our scattered-light modeling. This discrepancy from the blowout size could indicate a rapid resupply of small grains, possibly from ongoing collisions; however, this blowout size could
alternatively indicate that we have underestimated the mini-
mum grain size from radiative transfer and grain modeling due
to our assumptions on grain properties (Grigorieva et al. 2007)
and/or the use of Mie theory, which is known to be a poor
representation of true grain properties. However, if the grains
are, in fact, porous as indicated by our radiative transfer
modeling, the blowout size could be significantly larger,
making this discrepancy even larger. The blowout timescale is
approximately the free-fall time, which is short, on the order of
orbital timescales. As a result, the blowout timescale should
be the dominant timescale for this system. The collision
timescale, which dominates the dust replenishment, is
~$10^2$ yr—substantially longer than the blowout timescale.

As HD 156623 is known to host a gas-rich disk, gas drag
could change the blowout size and associated timescales—that
is, small dust grains may be coupled to the gas, enabling their
retention and explaining the observed inner edge and sub-
blowout size grains. In Lieman-Sifry et al. (2016), HD 156623
has one of the strongest CO $J = 2-1$ detections out of their
sample in Sco-Cen; perhaps this system has a particularly high
gas mass and therefore greater influence of gas drag on the dust
dynamics. The sample of Sco-Cen disks from Lieman-Sifry
et al. (2016) also suggests a correlation between the presence of
gas and sub-blowout grains, as seen in this system. Further
modeling is needed to fully understand the impact of the
presence of gas in this disk.

8. Conclusions

We have (1) presented analysis of the morphology of the
possible “hybrid” debris disk around HD 156623 in H-band
polarized scattered light, (2) presented an empirical measure-
ment of the brightness/surface density profile and SPF for
HD 156623’s debris disk, and (3) determined that the inner
radius of the disk is likely close to the star to <13.4 au via SED
modeling and imaging. Our work suggests that HD 156623
retains sub-blowout size grains, and gas drag may be a key
influence on the dust population of this debris disk, in line with
findings from Lieman-Sifry et al. (2016). However, we
acknowledge the numerous degeneracies present in this type
of debris disk data as the observed scattered-light distribution is
a complicated and difficult-to-disentangle combination of disk
geometry, grain size distribution, grain properties (e.g.,
porosity, scattering efficiency), and more (Stark et al. 2014).
Future work should consider obtaining additional mid-infrared
photometry (e.g., from JWST) to further constrain the debris
disk SED and testing more complex models, such as a two-
component dust population or simultaneously modeling the gas
and dust components of the disk to explore how gas drag
impacts small grain retention in this system.

Acknowledgments

This material is based upon work supported by the National
Science Foundation Graduate Research Fellowship under grant
No. 2016-21 DGE-1650604 and 2021-25 DGE-2034835. Any
opinions, findings, and conclusions or recommendations
expressed in this material are those of the authors(s) and do
not necessarily reflect the views of the National Science
Foundation. K.A.C. and B.C.M. acknowledge a Discovery
Grant from the Natural Science and Engineering Research
Council of Canada.

This work is based on data obtained by the international
Gemini Observatory, a program of NSF’s NOIRLab, which is
managed by the Association of Universities for Research in
Astronomy (AURA) under a cooperative agreement with the
National Science Foundation on behalf of the Gemini
Observatory partnership: the National Science Foundation
(United States), National Research Council (Canada), Agencia
Nacional de Investigación y Desarrollo (Chile), Ministerio de
Ciencia, Tecnología e Innovación (Argentina), Ministério da
Ciência, Tecnologia, Inovações e Comunicações (Brazil), and
Korea Astronomy and Space Science Institute (Republic of
Korea).

This publication makes use of data from the following: (1) the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, and NEOWISE, which is a project of the Jet Propulsion Laboratory/California Institute of Technology. WISE and NEOWISE are funded by the National Aeronautics and Space Administration; (2) the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation; (3) 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 is based on observations with AKARI, a JAXA project with the participation of ESA. The Infrared Astronomical Satellite (IRAS) was a joint project of the US, UK and the Netherlands; and (4) ADS/JAO.ALMA#2012.1.00688.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

This research has made use of NASA’s Astrophysics Data
System Bibliographic Services, the VizieR catalog access tool,
CDS, Strasbourg, France, and, for the following data sets, the
NASA/IPAC Infrared Science Archive, which is funded by
the National Aeronautics and Space Administration and
operated by the California Institute of Technology: Cutri
et al. (2003), Neugebauer et al. (2019), Wright et al. (2019),
AKARI (2020). This research has made use of the Spanish
Virtual Observatory (https://svo.cab.inta-csic.es) project
funded by MCIN/AEI/10.13039/501100011033/ through
grant PID2020-112949GB-I00.

This work also made use of Astropy: a community-
developed core Python package and an ecosystem of tools
and resources for astronomy (Astropy Collaboration et al.
2013, 2018, 2022) and the imgmasks package, developed by
B. Lewis, R. Lopez, E. Teng, and J. Williams as part of the
2021 Code/Astro workshop. Thanks to Elizabeth and Jason for
their work on this package, and to the Code/Astro leads
(namely Sarah Blunt and Jason Wang) for a great workshop.
The authors would also like to thank Dr. A. Meredith Hughes and Dr. John Carpenter for their helpful correspondence regarding the data used in Lieman-Sifry et al. (2016), Marcos M. Flores for showing B. Lewis that *np.arctan2* exists, Quinn Pratt and Joseph Marcinik for helpful discussions on interpolation and numerical integration, and Minghan Chen for helpful discussions on fitting for disk geometry. Thanks to Carlos Quiroz and Fredrik Rantakyro for assistance in obtaining seeing monitor records from Gemini Observatory for our data from 2019. B.L. would also like to acknowledge the invaluable labor of the maintenance and clerical staff at her institution and the observatories, whose contributions make scientific discovery a reality. This research was partially conducted on the traditional lands of the Gabrielino-Tongva people.

*Software:* NumPy (van der Walt et al. 2011), IPython (Perez & Granger 2007), Jupyter Notebooks (Kluyver et al. 2016), Matplotlib (Hunter 2007), Astropy (Astropy Collaboration et al. 2013, 2018, 2022), SciPy (Jones et al. 2001), pyKLIP (Wang et al. 2015), MCFOST (Pinte et al. 2006, 2009), pycmfost (Pinte 2022), diskmc (Esposito et al. 2018), diskmap (Stolker et al. 2016), imgmasks (Lewis et al. 2021), emcee (Foreman-Mackey et al. 2013, 2019), mcfost-python (S. Wolff et al. 2017, S. G. Wolff et al. 2017).

## ORCID iDs

Briley L. Lewis https://orcid.org/0000-0002-8984-4319

Michael P. Fitzgerald https://orcid.org/0000-0002-0176-8973

Thomas M. Esposito https://orcid.org/0000-0002-0792-3719

Pauline Arriaga https://orcid.org/0000-0001-6364-2834

Ronald López https://orcid.org/0000-0001-7646-4184

Katie A. Crotts https://orcid.org/0000-0003-4909-256X

Gaspard Duchêne https://orcid.org/0000-0002-5092-6464

Katherine B. Follette https://orcid.org/0000-0002-7821-0695

Justin Hom https://orcid.org/0000-0001-9994-2142

Paul Kalas https://orcid.org/0000-0002-6221-5360

Brenda C. Matthews https://orcid.org/0000-0003-3017-9577

Maxwell Millar-Blanchaer https://orcid.org/0000-0001-6205-9233

David J. Wilner https://orcid.org/0000-0003-1526-7587

Johan Mazoyer https://orcid.org/0000-0002-9133-3091

Bruce Macintosh https://orcid.org/0000-0003-1212-7538

## References

* AKARI, T. 2020, IPAC, 10, 26131

* Arnold, J. A., Weinberger, A. J., Videen, G., & Zubko, E. S. 2019, AJ, 157, 157

* Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167

* Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

* Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

* Backman, D., Gillett, F., & Witteborn, F. 1992, ApJ, 385, 670

* Beichman, C., Neugebauer, G., Habing, H., Clegg, P., & Chester, T. J. 1988, Infrared Astronomical Satellite (IRAS) Catalogs and Atlases Explanatory Supplement, Caltech

* Bessell, M. S., Castelli, F., & Plez, B. 1998, A&A, 333, 231

* Beust, H., Lagrange-Henri, A., Vidal-Madjar, A., & Ferlet, R. 1990, A&A, 236, 202

* Bonnefoy, M., Milli, J., Menard, F., et al. 2021, A&A, 655, A62

* Booth, M., Kennedy, G., Sibthorpe, B., et al. 2013, MNRAS, 428, 1263

* Brown, A. G., Vallenari, A., Prusti, T., et al. 2021, A&A, 649, A1

* Burns, J. A., Lamy, P. L., & Soter, S. 1979, Icar, 40, 1

* Castelli, F., & Kurucz, R. L. 2004, IAU Symp. 210, Modelling of Stellar Atmospheres, ed. N. Piskunov et al. San Francisco, CA (ASP) A20

* Cohen, M., Wheaton, W. A., & Megeath, S. 2003, AJ, 126, 1090

* Cox, A. N. 2015, Allen's Astrophysical Quantities (Berlin: Springer)

* Crotts, K. A., Matthews, B. C., Duchêne, G., et al. 2024, ApJ, 961, 245

* Cutri, R., Skrutskie, M., Van Dyk, S., et al. 2003, The IRSA 2MASS All-Sky Point Source Catalog (Washington, DC: NASA)

* Cutri, R., Wright, E. L., Conrow, T., et al. 2012, yCat, II/311

* De Rosa, R. J., Esposito, T. M., Gibbs, A., et al. 2020a, Proc. SPIE, 11447, 114475A

* De Rosa, R. J., Nguyen, M. M., Chilcote, J., et al. 2020b, JATIS, 6, 015006

* De Rosa, R. J., Nielsen, E. L., Blunt, S. C., et al. 2015, ApJL, 814, L3

* Draine, B. T. 2003, ARA&A, 41, 241

* Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89

* Engler, N., Lazzoni, C., Gratton, R., et al. 2020, A&A, 635, A19

* Esposito, T. M., Duchêne, G., Kalas, P., et al. 2018, AJ, 156, 47

* Esposito, T. M., Kalas, P., Fitzgerald, M. P., et al. 2020, AJ, 160, 24

* Foreman-Mackey, D., Farr, W. M., Sinha, M., et al. 2019, JOSS, 4, 1864

* Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

* Fujiwara, H., Ishihara, D., Oyabu, S., et al. 2010, AKARI/IRC Catalogue Release Note

* Gaia, C., Brown, A., Vallenari, A., et al. 2018, A&A, 616, 22

* Gaia Collaboration 2018, yCat, I/345

* Goodman, J., & Weare, J. 2010, Camcos, 5, 65

* Grigorieva, A., Artymowicz, P., & Thébault, P. 2007, A&A, 461, 537

* Hales, A., Gorti, U., Carpenter, J. M., Hughes, M., & Flaherty, K. 2019, ApJ, 878, 113

* Henyey, L. G., & Greenstein, J. L. 1941, ApJ, 93, 70

* Hog, E., Fabricius, C., Makarov, V. V., et al. 2000, A&A, 355, L27

* Hög, E., Fabricius, C., Makarov, V. V., et al. 2000, A&A, 355, L27

* Hom, J., Patience, J., Esposito, T. M., et al. 2020, AJ, 159, 31

* Houk, N., & Fuentes-Williams, T. 1982, BAAS, 14, 615

* Hughes, A. M., Duchêne, G., & Matthews, B. C. 2018, ARA&A, 56, 541

* Hung, L.-W., Bruzzone, S., Millar-Blanchaer, M. A., et al. 2016, Proc. SPIE, 9908, 1044

* Hunter, J. D. 2007, CSE, 9, 90

* Ishihara, D., Onaka, T., Kataza, H., et al. 2010, A&A, 514, A1

* Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open Source Scientific Tools for Python. http://www.scipy.org/

* Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in Positioning and Power in Academic Publishing: Players, Agents and Agendas, ed. F. Loizides & B. Schmidt (Amsterdam: IOS Press), 87

* Konopacky, Q. M., Thomas, S. J., Macintosh, B. A., et al. 2014, Proc. SPIE, 9147, 2585

* Kral, Q., Marino, S., Wyatt, M. C., Kama, M., & Matra, L. 2019, MNRAS, 489, 3670

* Krist, J. E., Ardila, D., Golimowski, D., et al. 2005, AJ, 129, 1008

* Lagrange, A. M., Bonnefoy, M., Chauvin, G., et al. 2010, Sci, 329, 57

* Lewis, B. L., Teng, E., & rlopezcla 2021, briley-lewis/elliptical-mask, v0.2-beta, Zenodo, doi:10.5281/zenodo.502948

* Lieman-Sifry, J., Hughes, A. M., Carpenter, J. M., et al. 2016, ApJ, 828, 25

* Macintosh, B., Graham, J. R., Ingraham, P., et al. 2014, PNAS, 111, 12661

* Macintosh, B. A., Graham, J. R., Palmer, D. W., et al. 2008, Proc. SPIE, 7015, 701518

* Mainzer, A., Bauer, J., Grav, T., et al. 2011, ApJ, 731, 53

* Mann, I., Köhler, M., Kimura, H., Cechowski, A., & Minato, T. 2006, A&ARV, 13, 159

* Marino, S., Flock, M., Henning, T., et al. 2020, MNRAS, 492, 4409

* Marois, C., Lafreniere, D., Doyon, R., Macintosh, B., & Nadeau, D. 2006, ApJ, 641, 556

* Mellon, S. N., Mamajek, E. E., Zwintz, K., et al. 2019, ApJ, 870, 36

* Millar-Blanchaer, M. A., Graham, J. R., Pueyo, L., et al. 2015, ApJ, 811, 18

* Milli, J., Mouillet, D., Lagrange, A.-M., et al. 2012, A&A, 545, A111

* Moór, A., Curé, M., Kóspál, Á., et al. 2017, ApJ, 849, 123

* Moór, A., Kral, Q., Ábrahám, P., et al. 2019, ApJ, 884, 108

* Morales, F. Y., Rieke, G., Werner, M., et al. 2011, ApJL, 730, L29

* Murakami, H., Baba, H., Barthel, P., et al. 2007, PASJ, 59, S369

* Neugebauer, G., Habing, H., van Duinen, R., et al. 2019, IRAS Point Source Catalog, v2.1, IPAC, doi:10.26131/IRSA4

* Neugebauer, G., Habing, H. J., van Duinen, R., et al. 1984, ApJL, 278, L1

* Nielsen, E. L., De Rosa, R. J., Macintosh, B., et al. 2019, AJ, 158, 13
Pawellek, N., & Krivov, A. V. 2015, MNRAS, 454, 3207

Pawellek, N., Wyatt, M., Matra, L., Kennedy, G., & Yelverton, B. 2021, MNRAS, 502, 5390

Pearce, T. D., Krivov, A. V., Sefilian, A. A., et al. 2024, MNRAS, 527, 3876

Pecaut, M. J., & Mamajek, E. E. 2016, MNRAS, 461, 794

Perez, F., & Granger, B. E. 2007, CSE, 9, 21

Péricaud, J., Di Folco, E., Dutrey, A., Guilloteau, S., & Piétu, V. 2017, A&A, 600, A62

Perrin, M. D., Duchene, G., Millar-Blanchaer, M., et al. 2015, ApJ, 799, 182

Perrin, M. D., Graham, J. R., Larkin, J. E., et al. 2010, Proc. SPIE, 7736, 1964

Perrin, M. D., Ingraham, P., Follette, K. B., et al. 2016, Proc. SPIE, 9908, 1010

Perrin, M. D., Maire, J., Ingraham, P., et al. 2014, Proc. SPIE, 9147, 1168

Pickles, A., & Depagne, É. 2010, PASP, 122, 1437

Pinte, C. 2022, pymcfost, GitHub, https://github.com/cpinte/pymcfost

Pinte, C., Harries, T., Min, M., et al. 2009, A&A, 498, 967

Pinte, C., Ménard, F., Duchêne, G., & Bastien, P. 2006, A&A, 459, 797

Ratzenböck, S., Großschedl, J. E., Alves, J., et al. 2023, A&A, 678, A71

Ren, B., Pueyo, L., Zhu, G. B., Debes, J., & Duchêne, G. 2018, ApJ, 852, 104

Rizzuto, A. C., Ireland, M. J., & Zucker, D. B. 2012, MNRAS: Letters, 421, L97

Rodrigo, C., & Solano, E. 2020, in XIV. 0 Scientific Meeting (virtual) of the Spanish Astronomical Society (Madrid: SEA), 182

Roques, F., Scholl, H., Sicardy, B., & Smith, B. A. 1994, Icar, 108, 37

Schmid, H., Joos, F., & Tschan, D. 2006, A&A, 452, 657

Sivaramakrishnan, A., & Oppenheimer, B. R. 2006, ApJ, 647, 620

Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163

Soummer, R., Pueyo, L., & Larkin, J. 2012, ApJL, 755, L28

Stark, C. C., Schneider, G., Weinberger, A. J., et al. 2014, ApJ, 789, 58

Stolker, T., Dominik, C., Min, M., et al. 2016, A&A, 596, A70

TanabÉ, T., Sakon, I., Cohen, M., et al. 2008, PASJ, 60, S375

Thébault, P. 2009, A&A, 505, 1269

Tokunaga, A., & Vacca, W. 2005, PASP, 117, 421

van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, CSE, 13, 22

Wang, J. J., Perrin, M. D., Savransky, D., et al. 2018, JATIS, 4, 018002

Wang, J. J., Rajan, A., Graham, J. R., et al. 2014, Proc. SPIE, 9147, 1678

Wang, J. J., Ruffio, J.-B., De Rosa, R. J., et al., 2015 pyKLIP: PSF Subtraction for Exoplanets and Disks, Astrophysics Source Code Library, ascl:1506.001

Wiktorowicz, S. J., Millar-Blanchaer, M., Perrin, M. D., et al. 2014, Proc. SPIE, 9147, 2574

Wolff, S., Perrin, M., Ren, B., & Pinte, C. 2017, swolff9/mcfost-python: mcfost-python-v1.0, Zenodo, doi:10.5281/ZENODO.839863

Wolff, S. G., Perrin, M. D., Stapelfeldt, K., et al. 2017, ApJ, 851, 56

Wright, E. L., Eisenhardt, P. R., Mainzer, A. K., Ressler, M. E., & Cutri, R. M. 2019, AllWISE Source Catalog, IPAC, doi:10.26131/IRSA1

Wright, E. L., Eisenhardt, P. R., Mainzer, A. K., et al. 2010, AJ, 140, 1868

Wyatt, M. C. 2008, ARA&A, 46, 339