# The Inner Disk Rim of HD 163296: Linking Radiative Hydrostatic Models with Infrared Interferometry

Ondřej Chrenko, Mario Flock, Takahiro Ueda, Antoine Mérand, Myriam

Benisty, Raúl O. Chametla

► To cite this version:

Ondřej Chrenko, Mario Flock, Takahiro Ueda, Antoine Mérand, Myriam Benisty, et al. The Inner Disk Rim of HD 163296: Linking Radiative Hydrostatic Models with Infrared Interferometry. The Astronomical Journal, 2024, 167, 10.3847/1538-3881/ad234d. insu-04836894

HAL Id: insu-04836894

https://insu.hal.science/insu-04836894v1

Submitted on 13 Dec 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The Inner Disk Rim of HD 163296: Linking Radiative Hydrostatic Models with Infrared Interferometry

Ondřej Chrenko¹, Mario Flock² and Takahiro Ueda², Antoine Mérand³, Myriam Benisty⁴,⁵ and Raúl O. Chametla¹

¹Charles Univ, Fac Math & Phys, Astronomical Institute, V Holešovičkách 747/2, 180 00 Prague 8, Czech Republic; chrenko@sirrah.troja.mff.cuni.cz

²Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany

³European Southern Observatory, Karl-Schwarzschild-Str. 2, 85748 Garching, Germany

⁴Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, France

⁵Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France

Received 2023 July 19; revised 2024 January 25; accepted 2024 January 26; published 2024 February 22

## Abstract

Previous studies of the protoplanetary disk HD 163296 revealed that the morphology of its sub-au infrared emission encompasses the terminal sublimation front of dust grains, referred to as the inner rim, but also extends into the (supposedly) dust-free region within it. Here, we present a set of radiative hydrostatic simulations of the inner rim in order to assess how much the rim alone can contribute to the observed interferometric visibilities $V$, half-light radii $R_{hl}$, and fractional disk fluxes $\mathcal{F}$ in the wavelength range 1.5–13 $\mu$m. In our set of models, we regulate the cooling efficiency of the disk via the boundary condition for radiation diffusion and we also modify the shape of the sublimation front. We find that when the cooling efficiency is reduced, the infrared photosphere at the rim becomes hotter, leading to an increase in $R_{hl}$ sufficient to match the observations. However, the near-infrared disk flux is typically too low ($\mathcal{F} \approx 0.25$ at 1.5 $\mu$m), resulting in H-band visibility curves located above the observed data. We show that the match to the H-band observations up to moderate baselines can be improved when a wall-shaped rather than curved sublimation front is considered. Nevertheless, our model visibilities always exhibit a bounce at long baselines, which is not observed, confirming the need for additional emission interior to the rim. In summary, our study illustrates how the temperature structure and geometry of the inner rim need to change in order to boost the rim's infrared emission.

*Unified Astronomy Thesaurus concepts: Protoplanetary disks (1300); Planet formation (1241); Radiative transfer (1335); Interferometry (808); Exoplanets (498)*

## 1. Introduction

Sub-au regions of protoplanetary disks represent the environment that has shaped precursors of terrestrial planets as well as the numerous populations of short-period exoplanets (e.g., Mulders et al. 2018; Petigura et al. 2018) during their early evolution. For instance, the transition between the outer dead zone and the inner zone of active magnetorotational (MRI) turbulence (at $\approx 900$ K) is considered a sweet spot for accumulation of dust grains (e.g., Varnière & Tagger 2006; Dzyurkevich et al. 2010; Ueda et al. 2019; Jankovic et al. 2022) as well as a migration trap for planets (e.g., Masset et al. 2006; Flock et al. 2019), although the local suppression of the dust drift and planet migration seems to strongly depend on disk properties (e.g., Schobert et al. 2019; Jankovic et al. 2021; Chrenko et al. 2022).

One of the possibilities to study sub-au disk regions using observations lies in the emission of the terminal sublimation front of dust grains, hereinafter referred to as the inner disk rim (Dullemond & Monnier 2010). As the grains at the rim equilibrate close to their sublimation temperature, their thermal emission can become an important contributor to the near-infrared (NIR) excess of Herbig Ae and Be stars (Hillenbrand et al. 1992; Lada & Adams 1992; Millan-Gabet et al. 2001; Natta et al. 2001). According to the pioneering models of the inner rim (Dullemond et al. 2001), no dust grains should exist

inwards from the sublimation radius<fnref n="6" />, and the rim should remain exposed to irradiation from the central star and heated, thus becoming wall-shaped and puffed up. The size of the dust-free region within the rim radius then scales with the square root of the stellar luminosity (so-called size–luminosity relation; Monnier & Millan-Gabet 2002).

The spatial morphology of the NIR and mid-infrared (MIR) inner disk emission, accessible through the advent of interferometric techniques (with difficulties related to instrumentation limitations and data sparseness), should in principle trace the inner rim geometry, manifesting itself via a bright ring or a torus. While this is sometimes the case and the torus-like emission is indeed observed (Tuthill et al. 2001; Monnier & Millan-Gabet 2002; Monnier et al. 2005), there are also cases when an additional emission source located somewhere between the magnetospheric cavity and the dust sublimation radius is required (e.g., Eisner et al. 2007; Kraus et al. 2008).

The emission of inner regions of the protoplanetary disk HD 163296, which is the subject of this work, is similarly puzzling. Tannirkulam et al. (2008) and Benisty et al. (2010) found that if only the inner rim torus-like emission is considered, (i) the visibility curve in *H* and *K* bands exhibits a bounce at long baselines inconsistent with observations, (ii) the observed NIR excess can be recovered only partially, (iii) and the closure phases typically become too large. By adding a smooth emission source at radii inside the actual dust rim, they were able to make the visibility curves featureless, increase the

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

<footnote>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.</footnote>

<footnote n="6">⁶ We point out, however, that large dust grains can cool down efficiently and thus can exist inwards from the conventional sublimation radius (Kama et al. 2009; Klarmann 2018).</footnote>
NIR flux, and reduce the closure phase signal. Among possible explanations of the additional emission component is the optically thin emission of either the hot gas (Tannirkulam et al. 2008) or refractory dust grains (Vinković et al. 2006; Benisty et al. 2010), but note that a predictive physical model for neither has ever been put forward. Further evidence for emission interior to the dust rim was obtained by Setterholm et al. (2018) who performed morphological fitting of the CHARA and VLTI interferometric data to constrain the brightness distribution profile of HD 163296 and concluded that the best-fitting model is a Gaussian-like 2D disk centrally peaked at the star location, without any strong indication of a sharp dust sublimation radius. The same conclusion was reached in Kluska et al. (2020) by means of image reconstruction (but one should bear in mind the resolution-related issues of image reconstruction on sub-au scales).

In another line of studies, based primarily on fitting prescribed parametric brightness distributions directly to the visibility data, it was found that the best-fitting model for HD 163296 in the VLTI bands H (Lazareff et al. 2017), K (GRAVITY Collaboration et al. 2019, 2021), L, and N (Varga et al. 2021) is a wide ring with an azimuthal modulation. The azimuthal modulation was found to be evolving with time (Kobus et al. 2020; GRAVITY Collaboration et al. 2021; Varga et al. 2021), possibly pointing to the presence of a vortex, a warp, or a variability in the launching zone of the disk wind (Bans & Königl 2012). The resulting width of the emitting ring was again found to be extending within the sublimation radius of dust grains, reconfirming that the inner rim is not the only contributor to the NIR and MIR excess in HD 163296.

In summary, it is clear that the inner disk emission of HD 163296 can possibly have two components: one arising from the inner rim and one (of unknown origin) from the region inside the rim. However, it remains unclear how the two components compare to one another—is their contribution equally important, or is one of them dominant? The question remains unsettled mostly because the recently used parametric fits (Lazareff et al. 2017; GRAVITY Collaboration et al. 2021; Varga et al. 2021) employ only a handful of parameters to avoid degeneracy and they are also difficult to link directly to physical models. With this in mind, the strategy of our paper is to start from a physical model of the inner rim alone (following the framework developed by Flock et al. 2016, 2019) and see how it compares to the visibility profiles in multiple NIR and MIR bands, to the half-light radii determined in earlier works, and to previously reported fractional disk fluxes. Our objective is to answer what it takes to modify the physical model in order to push some of the synthetic observables closer to the real data. We mostly focus on modifying the cooling efficiency of the disk and the shape of the sublimation front.

The aim of our study is by no means to explain the interferometric observations fully (as we do not model the emission component inside the rim), nor describe the temporal variability of the inner disk asymmetry (as our model is by construction static and symmetric). It is rather to set the groundwork for follow-up studies to help distinguish how much the inner rim can contribute to the interferometric signals. In the future, our models can be readily combined with morphological fitting (e.g., by parameterizing an azimuthal asymmetry on top of one of our base models), or they can help tweak the relative contribution between the rim and the interior

emission when a physical description of the latter becomes available.

The manuscript is structured as follows. We describe the radiative hydrostatic method for deriving the structure of the inner rim in Section 2.1. The list of nominal parameters is given in Section 2.2 where we also summarize our individual models, their boundary conditions, and assumptions for the sublimation temperature of dust grains. Section 2.3 gives an overview of observables and provides a discussion of a theoretical link between interferometric visibilities and half-light radii. Our results are presented in Section 3 and the paper is concluded in Section 4. Appendix A is devoted to demonstrating the importance of boundary conditions and the convergence of our models is discussed in Appendix B.

## 2. Method

### 2.1. Radiative Hydrostatic Disk Structure

We use the radiative hydrostatic approach of Flock et al. (2016, 2019) to calculate the distribution of the disk gas density $\rho$, dust density $\rho_d$, and temperature $T$. Our implementation was done in the FARGO3D code (Benítez-Llambay & Masset 2016), extending the work of Chrenko & Nesvorný (2020). The model relies on a decoupling between the timescales of thermal relaxation (driven mostly by the radiation reprocessing), vertical hydrostatic relaxation (driven by the propagation of sound waves), and disk accretion (driven by the redistribution of the angular momentum).

The decoupling makes it possible to proceed iteratively and one iteration can be summarized as follows:

1. Starting with an initial guess of $T$ and keeping it fixed, find $\rho$ in hydrostatic equilibrium (see Section 2.1.1 for constraints and details).

2. Perform 10 subiterations of:

(a) Radially integrated optical depths to stellar irradiation $\tau$ (which depend on the dust-to-gas ratio $f_{d2g}$ that sets the local optical depth of dust in each grid cell; see Section 2.1.2 and Equation (5)),

(b) Dust-to-gas ratio $f_{d2g}$ (which has to depend on $\tau$ in order to properly resolve irradiation absorption; see Section 2.1.3 and Equation (6))

3. Keeping $\rho$, $\tau$, and $f_{d2g}$ fixed, evolve two-temperature energy Equations (7) and (8) over a reasonably chosen time step $dt$ while accounting for the disk heating due to stellar irradiation and viscous heating (Section 2.1.4).

4. Return to the beginning with the new temperature field.

#### 2.1.1. Density Distribution of Gas

At the beginning of each iteration, we fix the temperature field and solve the equations of the radial-vertical hydrostatic equilibrium in spherical coordinates. In a compact form (e.g., Chrenko & Nesvorný 2020), one can write

$$ \frac{\partial P}{\partial \phi} = r \left( \frac{\partial P}{\partial r} + \rho \frac{GM_*}{r^2} \right) \frac{1}{\tan \phi}, \quad (1) $$

where $P$ is the thermal pressure, $\phi$ is the colatitude, $r$ is the radius, $G$ is the gravitational constant, and $M_*$ is the mass of the central star. The azimuthal dimension $\theta$ is ignored, assuming an axisymmetric solution. The radial spacing of the grid is logarithmic and the vertical spacing is equidistant.
Equation (1) can only be solved along with suitable closure relations, the first one being the ideal gas equation of state

$$P = (\gamma - 1)\epsilon = (\gamma - 1)\rho c_v T, \quad (2)$$

where $\gamma$ is the adiabatic index, $\epsilon$ is the internal energy density of gas, and $c_v$ is the specific heat at constant volume.

Additionally, we assume that the disk is viscously evolving and its mass accretion rate $\dot{M}$ is uniform. Then the equation

$$\dot{M} = 3\pi\nu\Sigma, \quad (3)$$

where $\nu$ is the effective viscosity, provides a constraint on the gas surface density $\Sigma$. The viscosity is parametrized via the Shakura & Sunyaev (1973) prescription $\nu = \alpha c_s^2 / \Omega_K$, where $c_s = \sqrt{\gamma P / \rho}$ is the adiabatic sound speed and $\Omega_K$ is the Keplerian angular frequency. We point out that a density-weighted vertical average of $\nu$ is used when evaluating Equation (3) (see Chrenko & Nesvorný 2020). Furthermore, our $\alpha$-parameterization mimics the ionization transition that separates an inner region where the MRI is active and an outer dead zone (DZ) where the MRI is suppressed. Following Flock et al. (2016), we write

$$\alpha(T) = (\alpha_{\text{MRI}} - \alpha_{\text{DZ}}) \left[ \frac{1 - \tanh\left(\frac{T_{\text{MRI}} - T}{70 \text{ K}}\right)}{2} \right] + \alpha_{\text{DZ}}, \quad (4)$$

where $T_{\text{MRI}}$ is the transition temperature, and $\alpha_{\text{DZ}}$ and $\alpha_{\text{MRI}}$ are $\alpha$-viscosities in the dead and active zones, respectively. Finally, we assume that $\Sigma$ is related to the volume density in the midplane $\rho_{\text{mid}} = \Sigma / (\sqrt{2\pi H})$ as if the disk was vertically isothermal, with $H = c_s / (\sqrt{\gamma} \Omega_K)$ being the pressure scale height ($H$ is evaluated from $c_s$ in the midplane for the purpose of estimating $\rho_{\text{mid}}$).

Equipped with the aforementioned closure relations, it is possible to reconstruct the radial profile of $\rho_{\text{mid}}(r)$ for a fixed temperature field. This serves as a starting point for solving Equation (1) and thus finding $\rho(r, \phi)$ throughout the rest of the disk in each iteration (see Flock et al. 2016; Chrenko & Nesvorný 2020).

### 2.1.2. Opacities

Before the energy (or temperature) is advanced in our iterative scheme, it is necessary to determine the opacities in each cell. As in Flock et al. (2016), we use a simple three-opacity model and we assume that the Planck ($\kappa_{\text{P}}$) and Rosseland ($\kappa_{\text{R}}$) opacities are the same. We define the gas opacity $\kappa_{\text{gas}}$, the dust opacity to its own thermal emission $\kappa_{\text{d}}(T_s)$, and the dust opacity to stellar irradiation $\kappa_{\text{d}}(T_*)$ (where $T_s$ and $T_*$ are the sublimation and stellar temperatures, respectively). The optical depth to stellar irradiation is then calculated along radial rays as

$$\tau(r) = \tau_0 + \int_{r_{\text{in}}}^{r} (\kappa_{\text{gas}}\rho + \kappa_{\text{d}}(T_*)\rho_{\text{d}}) dr', \quad (5)$$

where $\tau_0$ is the optical depth inwards from our computational grid, at $r < r_{\text{in}}$ (see Flock et al. 2016). We point out that $\rho_{\text{d}}$ does not necessarily represent the total dust content but mainly accounts for small grains, which are the dominant opacity contributors.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Top: Opacity of dust grains as a function of the wavelength $\lambda$. We show the absorption opacity $\kappa_{\text{abs}}$ (solid black curve) and the scattering opacity $\kappa_{\text{sca}}$ (dashed gray curve). Dotted blue horizontal lines mark the values of Planck-averaged opacities to thermal emission $\kappa_{\text{d}}(T_s)$ and to stellar irradiation $\kappa_{\text{d}}(T_*)$. Planck-averaged opacities are utilized in our hydrostatic models with radiative diffusion, while the wavelength-dependent opacities are used to ray trace synthetic images using RADMC-3D. Bottom: Stellar spectrum used in RADMC-3D. Individual points represent the frequency sampling.</figcaption></figure>

As for the actual value of $\kappa_{\text{gas}} = 10^{-5} \text{ cm}^2 \text{ g}^{-1}$, we set it very low in order to maintain the innermost dust-free disk regions optically thin (we refer the reader to Appendix B of Flock et al. 2019). To determine $\kappa_{\text{d}}(T_s)$ and $\kappa_{\text{d}}(T_*)$, we first calculated wavelength-dependent dust opacities. We assumed that the dust grains are composed of 62.5% astronomical silicate (Draine 2003) and 37.5% amorphous carbon (Preibisch et al. 1993), having a distribution of physical sizes $f(a) \propto a^{-3.5}$ ranging between $3 \times 10^{-3}$ and $10^2 \mu\text{m}$. Using the OPTOOL code<fnref n="7" /> (Dominik et al. 2021), we obtained the wavelength-dependent opacities shown in the top panel of Figure 1. Subsequently, we calculated $\kappa_{\text{d}}(T_s)$ as the Planck opacity at 1400 K (which is a proxy of the temperature in the dusty disk halo) and $\kappa_{\text{d}}(T_*)$ as the Planck opacity at 9000 K (which is the effective temperature of the irradiating central star).

The composition and size distribution of dust grains at the inner rim are largely unconstrained and for simplicity, we chose them in analogy to some of the previous works (Turner et al. 2014; Flock et al. 2019). To a limited extent, we varied the optical dust properties in Section 3.5.

### 2.1.3. Density Distribution of Dust

We treat the dust grains as passive tracers of the gas and track their volumetric content using the dust-to-gas ratio

<footnote n="7">⁷ Our wavelength-dependent dust opacities can be reproduced with the following command: `optool astrosil 0.625 c-p 0.375 -a 0.003 100.0 3.5 -mie`.</footnote>
<table><caption>**Table 1**
Fiducial Parameters for the Radiative Hydrostatic Model</caption><tr><td>Grid size (radial vertical)</td><td>$N_r \times N_\phi = 4096 \times 256$</td></tr><tr><td>Opening angle of the domain</td><td>$\Delta\phi = \pm 0.24$ rad</td></tr><tr><td>Inner radial boundary</td><td>$r_{\text{in}} = 0.1$ au</td></tr><tr><td>Outer radial boundary</td><td>$r_{\text{out}} = 15$ au</td></tr><tr><td>MRI transition temperature</td><td>$T_{\text{MRI}} = 900$ K</td></tr><tr><td>MRI-active viscosity</td><td>$\alpha_{\text{MRI}} = 10^{-1}$</td></tr><tr><td>Dead-zone viscosity</td><td>$\alpha_{\text{DZ}} = 10^{-3}$</td></tr><tr><td>Dust opacity at $T_s$</td><td>$\kappa_d(T_s) = 751 \text{ cm}^2 \text{ g}^{-1}$</td></tr><tr><td>Dust opacity at $T_*$</td><td>$\kappa_d(T_*) = 1878 \text{ cm}^2 \text{ g}^{-1}$</td></tr><tr><td>Gas opacity</td><td>$\kappa_{\text{gas}} = 10^{-5} \text{ cm}^2 \text{ g}^{-1}$</td></tr><tr><td>Maximum dust-to-gas ratio</td><td>$f_{\text{d2g,max}} = 10^{-3}$</td></tr><tr><td>Minimum dust-to-gas ratio</td><td>$f_{\text{d2g,min}} = 10^{-10}$</td></tr><tr><td>Mean molecular weight</td><td>$\mu = 2.3$</td></tr><tr><td>Adiabatic index</td><td>$\gamma = 1.43$</td></tr><tr><td>Stellar temperature</td><td>$T_* = 9000$ K</td></tr><tr><td>Stellar radius</td><td>$R_* = 1.87 R_\odot$</td></tr><tr><td>Stellar mass</td><td>$M_* = 1.95 M_\odot$</td></tr><tr><td>Stellar luminosity</td><td>$L_* \approx 20.6 L_\odot$</td></tr><tr><td>Mass accretion rate</td><td>$\dot{M} = 3.24 \times 10^{-8} M_\odot \text{ yr}^{-1}$</td></tr></table>

Note. Stellar parameters of HD 163296 are adopted from Wichittanakom et al. (2020) and the mass accretion rate is from Fairlamb et al. (2015).

$f_{\text{d2g}}(r, \phi) = \rho_d / \rho$ calculated as (Flock et al. 2019)

$$f_{\text{d2g}} = \begin{cases} f_{\text{d2g,max}} & T < T_s \land \tau > 3 \\ \left[ \frac{1 - \tanh\left(\left(\frac{T-T_s}{100 \text{ K}}\right)^3\right)}{2} \right] & \text{otherwise} \\ \times \left[ \frac{1 - \tanh\left(\frac{2}{3} - \tau\right)}{2} \right], & \end{cases} \quad (6)$$

where $f_{\text{d2g,max}}$ is the maximum dust-to-gas ratio and $T_s$ is the sublimation temperature of dust grains (see Section 2.2). The value of $f_{\text{d2g,max}}$ is somewhat lower (see Table 1) compared to the canonical value of $10^{-2}$ to reflect the fact that the growth of dust grains depletes the sub-$\mu$m-sized grains (Birnstiel et al. 2012). To prevent numerical problems, we also define a floor value $f_{\text{d2g}} > f_{\text{d2g,min}}$. Additionally, the stability of our method is improved by ramping $f_{\text{d2g,max}}$ from $f_{\text{d2g,min}}$ up to the desired value over the first 25 iterations (similarly to Schobert et al. 2019).

The term $f_{\Delta\tau} = 0.2/(\rho_d\kappa_d\Delta r)$ regulates the maximum increase in the optical depth to stellar irradiation per one grid cell with the radial length $\Delta r$ and allows us to resolve the transition between optically thin and thick regions even with a coarse grid spacing (see also Kama et al. 2009). We also impose an upper limit $f_{\Delta\tau} \le f_{\text{d2g,max}}$ to prevent $f_{\Delta\tau}$ from becoming too large in regions with low $\rho_d$. As $\tau$, $f_{\text{d2g}}$ and opacities are mutually dependent through Equations (5) and (6), we perform 10 subiterations within each iteration to evaluate them.

### 2.1.4. Evolving the Temperature

To finish one iteration, we search for a new temperature field corresponding to the hydrostatic distribution of gas and dust. This is done by integrating the coupled set of energy equations describing the evolution of $\epsilon$ and the energy density of thermal

radiation field $E_R$ (Dobbs-Dixon et al. 2010):

$$\frac{\partial \epsilon}{\partial t} = -\rho \kappa_P [4\sigma T^4 - cE_R] + Q_{\text{irr}} + Q_{\text{visc}}, \quad (7)$$

$$\frac{\partial E_R}{\partial t} + \nabla \cdot \vec{F} = \rho \kappa_P [4\sigma T^4 - cE_R], \quad (8)$$

where $\kappa_P = \kappa_{\text{gas}} + f_{\text{d2g}}\kappa_d(T_s)$, $\sigma$ is the Stefan-Boltzmann constant, $c$ is the speed of light, $Q_{\text{irr}}$ is the irradiation heating rate, $Q_{\text{visc}}$ is the viscous heating rate, and $\vec{F}$ is the radiation flux. The radiative energy is transported using the flux-limited diffusion approximation (Levermore & Pomraning 1981) with the flux limiter of Kley (1989). Equations (7) and (8) are solved in an implicit form (Bitsch et al. 2013; Chrenko & Lambrechts 2019) using a simple successive over-relaxation method. In our iteration scheme, the time step to advance Equations (7) and (8) is $dt_1 = 10^5$ s during the first 100 iterations, followed by 100 iterations with $dt_2 = 10^8$ s. We let $f_{\text{d2g}}$ evolve only during the first 100 iterations; afterward it remains fixed. The number of iterations and time step sizes are chosen empirically: $dt_1$ is short enough to avoid convergence problems when the disk is being gradually filled with dust and $dt_2$ is long enough to bring the most optically thick disk regions to thermal equilibrium by radiation diffusion.

The heating due to the absorption of stellar photons (e.g., Bitsch et al. 2013; Kolb et al. 2013; Chrenko & Nesvorný 2020) is

$$Q_{\text{irr}} = \frac{L_*}{4\pi r^2} (e^{-\tau} - e^{-(\tau+d\tau)}) \frac{S_{\text{cell}}}{V_{\text{cell}}}, \quad (9)$$

where $L_*$ is the stellar luminosity, $d\tau$ is the increment of the optical depth across a grid cell of interest, $S_{\text{cell}}$ is the irradiated cross section of the cell, and $V_{\text{cell}}$ is its volume.

The viscous heating term is (e.g., D'Angelo & Bodenheimer 2013)

$$Q_{\text{visc}} = \frac{1}{2\nu\rho} T_{ij} T^{ij}, \quad (10)$$

where $T_{ij}$ are the components of the viscous stress tensor. During our first experiments with $Q_{\text{visc}}$, we found that a straightforward implementation of this term leads to fluctuating (nonconverging) solutions due to the coupling with Equation (4). The coupling often results in a feedback loop at spurious locations across the inner disk rim—if $Q_{\text{visc}}$ manages to locally increase $T$ so that $\nu(\alpha, T)$ starts to increase, the local surface density starts to drop through $\Sigma = \dot{M}/(3\pi\nu)$, thus changing optical depths and unbalancing the system from thermal equilibrium. Moreover, the temperature fluctuations also directly affect the dust content via Equation (6) (see also Schobert et al. 2019). To circumvent the aforementioned issues, we considered uniform $\alpha = \alpha_{\text{DZ}}$ for the purpose of calculating $Q_{\text{visc}}$ (Schobert et al. 2019). Although this leads to an inner inconsistency in our model, we think it is a reasonable first approximation because $Q_{\text{visc}}$ is calculated correctly in the optically thick regions within the disk interior and the incorrect solution (with too low $\alpha$) applies mostly inwards from the disk rim where we expect $Q_{\text{irr}}$ to dominate anyway (see, for instance, Figure 9 in Flock et al. 2019).
## 2.2. Individual Models

Our main set of models revolves around modifications of the disk’s cooling efficiency and the shape of the dust sublimation front. The former is achieved by modifying the boundary condition for the radiation energy density $E_R$ while the latter is achieved by modifying the prescription for the sublimation temperature of dust grains $T_s$. Starting with $E_R$, we consider two sets of boundary conditions. The first set is referred to as the cold boundary and it sets (e.g., Schobert et al. 2019)

$$E_R = (1 - \exp(-\tau_{bc}))a_R T_{bc}^4, \quad (11)$$

at the inner radial boundary, with $a_R$ being the radiation constant, $\tau_{bc} = 10^{-2}$, and $T_{bc} = T_{\text{thin}}$ representing the temperature of optically thin gas. At the outer radial boundary, we prevent the diffusion of photons. At the boundaries in colatitude, we set $E_R = a_R(5\,\text{K})^4$, assuming a very low ambient temperature. The cold boundary is motivated by the fact that Equation (8) describes the evolution of the diffusing field of photons related to thermal radiation, while the field of irradiating photons is split and treated using an explicit absorption in Equation (7). The cold boundary, therefore, allows the diffusing photons to escape freely.

The second set is referred to as the warm boundary, and it assumes

$$E_R(r, \phi) = (E_R^{\text{thin}} - E_R^{\text{thick}}) \left[ \frac{1 - \tanh\left(\frac{r - r_{\text{fc}}}{0.1 r_{\text{fc}}}\right)}{2} \right] + E_R^{\text{thick}}, \quad (12)$$

where $r_{\text{fc}}$ is the radius of full dust condensation at all heights above the midplane (Equation (17) in Ueda et al. 2017). To evaluate $E_R^{\text{thin}}$, we use Equation (11) where $T_{bc}$ is set to the optically thin temperature of a dusty disk (Equation (1) in Ueda et al. 2017) and $\tau_{bc} = \min(\tau_{\text{mid}}, 1)$, $\tau_{\text{mid}}$ being the optical depth to stellar irradiation in the midplane. Finally, we use $E_R^{\text{thick}} = a_R T_{\text{thick}}^4$, with $T_{\text{thick}}$ corresponding to the surface temperature of an optically thick passively irradiated disk (Equations (11)–(15) in Ueda et al. 2017). The warm boundary sets a shallower gradient of $E_R$ at the grid edge in colatitude, thus reducing the cooling efficiency of the disk.

The purpose and influence of the boundary conditions are further demonstrated and discussed in Appendix A. The cold boundary leads to disks with temperature profiles similar to thermal Monte Carlo simulations. The warm boundary leads to temperature profiles similar to Flock et al. (2016).

Regarding $T_s$, we either consider the dust sublimation temperature of silicate grains (Pollack et al. 1994; Isella & Natta 2005)

$$T_s = 2000\,\text{K} \left( \frac{\rho}{1\,\text{g cm}^{-3}} \right)^{0.0195} \quad (13)$$

or we set it to a uniform and density-independent value of $T_s = 1350\,\text{K}$ (or $1550\,\text{K}$ in Section 3.5). The purpose of Equation (13) is to account for the change of sublimation conditions with the height above the midplane, which then leads to a curved inner rim (Kama et al. 2009), while the

purpose of the uniform sublimation temperature is to produce a wall-shaped rim geometry.

If not specified otherwise, all our models use parameters from Table 1. Differences between individual models are specified in Table 2. Basically, we start from a nominal model M1 with a cold boundary and a curved rim. Then we go to model M2 by switching to the warm boundary. Keeping the warm boundary, we change the rim geometry to wall-like in model M3. Results of models M1–M3 constitute most of Section 3; model M3Fe with a different dust composition and a larger sublimation temperature is discussed in Section 3.5, before concluding the paper.

## 2.3. Diagnostics

### 2.3.1. Synthetic Images

We use the Monte Carlo radiative transfer code RADMC-3D (Dullemond et al. 2012) to post-process the results of our hydrostatic modeling. We use $\rho$, $\rho_d$, and $T$ obtained with the hydrostatic computations as direct inputs for ray tracing synthetic images of the inner disk. Although it is possible to recalculate $T$ in RADMC-3D using the thermal Monte Carlo method, we do not do so as we verified that the resulting temperature would be similar to our models with the cold boundary (see also Appendix A).

We use the same computational grid as for the hydrostatic calculations, thus imposing the axisymmetric approximation and the simplest isotropic scattering mode. We introduce two species in RADMC-3D. The first species with the density $\rho$ represents the gas, for which the absorption opacity at each wavelength is considered uniform, $\kappa_{\lambda} = \kappa_{\text{gas}}$, and the scattering opacity is neglected. The second species with the density $\rho_d = f_{\text{d2g}}\rho$ represents the dust and its $\kappa_{\lambda}$ is shown in Figure 1. RADMC-3D offers the possibility to thermalize all species together, and we apply this option. This approximation is incorrect in the optically thin dusty halo of the inner rim where the dust and gas should be decoupled; however, we apply it for the sake of consistency because our hydrostatic runs are thermalized as well (we use only one temperature to describe the gas and dust).<fnref n="8" />

The spectrum of the irradiating star (bottom panel of Figure 1) is adopted from the BOSZ database of stellar atmospheric models (Mészáros et al. 2012; Bohlin et al. 2017) and corresponds to the stellar parameters shown in Table 1, along with the surface gravity $\log g = 4$ (Wichittanakom et al. 2020), and $\text{Fe/H} = 0.2$ (Tilling et al. 2012). The wavelength sampling of generated photons, represented by the data points in the bottom panel of Figure 1, covers three log-spaced intervals with 100 samples in $0.05-7\,\mu\text{m}$, 100 samples in $7-25\,\mu\text{m}$, and 30 samples in $25-10^4\,\mu\text{m}$. We use $10^9$ photon packages for synthetic image calculations, while the number of scattering photons amounts to $10^8$ (for clarity, we emphasize that scattering is only considered in calculations with RADMC-3D). Synthetic images are produced at $\lambda \in (1.5, 1.75, 2, 2.2, 2.45, 2.8, 3.5, 4.2, 5.5, 8, 10.5, 13)\,\mu\text{m}$, covering NIR and MIR bands of the VLTI instruments, assuming the disk inclination $i = 52^{\circ}$ and position angle $\text{PA} = 143^{\circ}$ (Varga et al. 2021). The

<footnote n="8">⁸ The thermalization in RADMC-3D, together with the incorporation of the gas-representing species, ensures that the temperature profile in the dust-free inner disk and at the edges of the optically thin dusty halo remains consistent between our hydrostatic calculations and thermal Monte Carlo calculations (see Appendix A).</footnote>
<table><caption>**Table 2**
Overview of Individual Models</caption><tr><td>M1</td><td>Parameters from Table 1, cold boundary for E<sub>R</sub> (Section 2.2), T<sub>s</sub> given by Equation (13)</td></tr><tr><td>M2</td><td>Parameters from Table 1, warm boundary for E<sub>R</sub> (Section 2.2), T<sub>s</sub> given by Equation (13)</td></tr><tr><td>M3</td><td>Parameters from Table 1, warm boundary for E<sub>R</sub> (Section 2.2), uniform T<sub>s</sub> = 1350 K</td></tr><tr><td>M3Fe</td><td>As M3, but with a modified dust composition (Section 3.5) and uniform T<sub>s</sub> = 1550 K</td></tr></table>

resolution is $3 \times 10^3$ pixels along the image edge, i.e., 0.1 mas per pixel. Second-order ray tracing of RADMC-3D is utilized.

### 2.3.2. Observables

Using the synthetic images, we calculate the half-light radii $R_{hl}$, fractional disk fluxes $\mathcal{F} = F_{\text{disk}}/F_{\text{tot}} = F_{\text{disk}}/(F_{\text{disk}} + F_{\text{star}})$, and interferometric visibilities $V$ at various $\lambda$ in order to compare them to the real data. The flux from an individual image pixel is calculated simply by multiplying the local emission intensity $I_\nu$ with the pixel surface area and considering the distance of HD 163296 being $d = 101.5$ pc (Wichittanakom et al. 2020). Then, $F_{\text{disk}}$ is an integral over all pixels occupied by our disk model and $F_{\text{tot}}$ is an integral over the entire image.

The half-light radius is defined via (Leinert et al. 2004; Varga et al. 2021)

$$ \frac{F_{\text{disk}}}{2} = \int_{0}^{R_{hl}} 2\pi r' I_{\nu}(r') dr', \quad (14) $$

where $r'$ is the radius from the center of the image plane and the stellar flux is excluded from the calculation.

The synthetic interferometric visibilities are calculated using codes RADMC3DPY<fnref n="9" /> and PMOURED (Mérand 2022). When analyzing the visibilities, we deproject the baselines $B$ using (e.g., Tannirkulam et al. 2008)

$$ B_{\text{eff}} = B \sqrt{\cos \chi^2 + \sin \chi^2 \cos i^2}, \quad (15) $$

where $i$ is the disk's inclination and $\chi = PA_B - PA_{\text{major}}$ is the difference between the position angle of the given baseline configuration and the major axis of the on-sky disk projection.

### 2.3.3. Link between Half-light Radii and Interferometric Visibilities

Although the half-light radius is a secondary interferometric observable, we would like to point out that there is a theoretical argument for a link between $R_{hl}$ and visibilities that has not been fully appreciated in prior works. Assuming the disk is viewed face-on<fnref n="10" /> and only experiences radial variations in intensity, the interferometric visibility amplitude $V$ is the Hankel transform of the disk profile $I_\nu$ (for the spatial frequency $B/\lambda$), combined linearly with the unit visibility of the unresolved central star:

$$ V(B/\lambda) = \frac{\int_{0}^{\infty} 2\pi r I_{\nu}(r) J_{0}(2\pi r^{B}/\lambda) dr + F_{\text{star}}}{F_{\text{disk}} + F_{\text{star}}}. \quad (16) $$

The Bessel function $J_0$ is an oscillating and vanishing function, so for sufficiently large $B/\lambda$, the disk is fully resolved: the integral in Equation (16) vanishes to 0 and the visibility becomes constant as a function of the baseline. On the other

hand, for very small baselines, the disk is unresolved<fnref n="11" /> ($V \sim 1$). So as the baseline increases, the visibility decreases from 1 to the saturation value $F_{\text{star}}/(F_{\text{disk}} + F_{\text{star}})$. For an intermediate spatial frequency $B_2/\lambda$, the visibility reaches a midpoint, which can be measured if a sufficiently wide range of baseline lengths was explored. In that case:

$$ \int_{0}^{\infty} 2\pi r I_{\nu}(r) J_{0}\left(2\pi r \frac{B_{2}}{\lambda}\right) dr = 1/2 \int_{0}^{\infty} 2\pi r I_{\nu}(r) dr \\ = \frac{F_{\text{disk}}}{2}, \quad (17) $$

which is very similar to Equation (14) defining the half-light radius if rewritten as:

$$ \int_{0}^{\infty} 2\pi r I_{\nu}(r) S(r/R_{hl}) dr = \frac{F_{\text{disk}}}{2}, \quad (18) $$

with $S(x) = 1$ for $x \le 1$ and $S(x) = 0$ for $x > 1$. Combining the last 2 Equations (17) and (18)), we get:

$$ \int_{0}^{\infty} r I_{\nu}(r) J_{0}\left(2\pi r \frac{B_{2}}{\lambda}\right) dr = \int_{0}^{\infty} r I_{\nu}(r) S(r/R_{hl}) dr. \quad (19) $$

The exact relation between $B_2$ and $R_{hl}$ in principle depends on the exact profile $I_\nu$. To illustrate whether the dependence is strong or weak, we visualize it in Figure 2 for a variety of intensity profiles that are typically used to analyze interferometric data. Denoting $I_0$ and $R_0$ the unit intensity and radius, respectively, we consider power-law profiles $I(r) = I_0 r^{-k}$ truncated at $R_{\text{max}} = (2, 3, 5)R_0$, a Poisson-like profile $I(r) = I_0 r^{-k} e^{-r}$, flat disk profiles $I(r) = I_0$ extending over $r \in (1, k)R_0$, and Gaussian rings $I(r) = I_0 \exp(-(r-1)^2 k^2)$. Parameter $k \in (1, 4)$ modulates the shape or extent of individual profiles. Figure 2 shows where the intensity profiles reach the half-light radius (left panel), where the corresponding visibility curves reach their midpoint (middle panel), and how $B_2$ relates to $R_{hl}$ over the considered range of $k$ (right panel). We find that the dependence on the exact intensity profile is relatively weak, in the range

$$ R_{hl} [\text{mas}] = 41..51 \frac{\lambda [\mu\text{m}]}{B_2 [\text{m}]} \quad (20) $$

For the sake of clarity, let us emphasize that considerations in this section were based on simple parametric radial intensity profiles, while our physical disk models generally lead to more complex brightness distributions (e.g., Section 3.4) with azimuthal variations due to projection and radiative transfer

<footnote n="9">⁹ Available at https://www.ita.uni-heidelberg.de/~dullemond/software/radmc-3d/manual_rmcpy.</footnote>

<footnote n="10">¹⁰ This discussion is still valid in the case of an inclined disk, but the deprojected baseline should be used instead.</footnote>

<footnote n="11">¹¹ Studies focused on inner disk regions often assume another contribution of a large-scale over-resolved emission component referred to as halo (e.g., Lazareff et al. 2017; Setterholm et al. 2018). Such a component would result in an addition of $F_{\text{halo}}$ in the denominator of Equation (16). It might result in $V \lesssim 1$ at very small baselines, which would slightly modify the subsequent analysis of this section. However, we neglect this halo component and we also caution the reader not to confuse it with the optically thin dusty halo defined later in Section 3.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Comparison of $B_2^1 \times R_{hl}$ for various analytical profiles of $I_ν$. Each radial profile (left) is plotted as a function of $r/R_{hl}$. For each profile, the visibility is computed as a function of $B_2/B_1^2$ (center). Finally, we compare (right) $B_2^1 \times R_{hl}$ as a function of $k$ (which parameterizes the various intensity profiles) and find that it is fairly independent of the function used to represent $I_ν$.</figcaption></figure>

effects. However, one can also look at previously published studies to assess whether the exercise with which we obtained Figure 2 can be generalized. For instance, Varga et al. (2021) showed that the L-band visibility of HD 163296 saturates around 0.1, and $V=0.55$ is reached for $B_2^1[\text{m}]/\lambda[\mu\text{m}] \approx 15$. Their best-fit parametric model, which was a Gaussian ring with an azimuthal modulation (thus a 2D intensity distribution), predicts the half-light radius of 3.28 mas, leading to a constant of 49.2 which falls into our range derived in Equation (20).

## 3. Results

### 3.1. Disk Structure and Temperature Profiles

Figure 3 shows the 2D temperature distribution in disk models M1, M2, and M3. Additionally, it shows several physically distinct surfaces. The gray curve is the inner boundary of the dust halo and it shows where the dust grains start to condense in minor quantities (Flock et al. 2016). The outer edge of the halo is at the black curve where $T=T_s$ and the dust becomes fully condensed. The white curve is the surface where the optical depth given in Equation (5) attains unity and most of the irradiating stellar photons are absorbed. At the inner rim, the irradiation absorption surface nearly overlaps with the front of fully condensed dust. In the outer flaring disk, the irradiation absorption surface delimits the disk atmosphere from the disk interior. The dashed green curve marks the infrared photosphere, i.e., the surface from which most of the detectable thermal emission originates. This surface, however, is in principle dependent on the wavelength and the line of sight—for the purpose of Figure 3, we calculated the optical depth to infrared emission in the direction perpendicular to the midplane (as if the disk was viewed face on) and for the opacity $\kappa_d(T_s)$.

Figure 3 reveals that the inner disk structure of models M1 and M2 is (unsurprisingly) consistent with the general findings of Kama et al. (2009) and Flock et al. (2016), exhibiting a rounded-off irradiated inner rim, an optically thin region inwards from the rim, and a flaring disk (Chiang & Goldreich 1997) outwards from the rim. Model M3 contains the same regions but its sublimation surface has a nearly vertical wall-like shape, similar to the classical rim of Dullemond et al. (2001). Additionally, the dust halo of model M3 is nearly isothermal.

Figure 3 also shows the radial extent of the inner rim (see the blue arrows). The inner edge of the rim, $R_{rim}^{in}$, is defined as the location where the dust fully condenses and the irradiation absorption peaks in the midplane. The outer edge of the rim, $R_{rim}^{out}$, is more difficult to define. Flock et al. (2016) established $R_{rim}^{out}$ as the local maximum of the aspect ratio of the infrared photosphere. Our profile of the infrared photosphere, however, has a monotonically increasing aspect ratio. Therefore, we define $R_{rim}^{out}$ as the radial location where the infrared photosphere has $T=800$ K, which is a typical temperature found in Flock et al. (2016) at the outer edge of the rim. We emphasize that $R_{rim}^{in}$ and $R_{rim}^{out}$ are related to the physical rim size, while the characteristic radius of the infrared emission is defined differently (Section 2.3.2).

The extent of the inner rim differs between models M1 and M2; less so between M2 and M3. We found $R_{rim} = 0.29–0.36$ au for M1, $R_{rim} = 0.3–0.46$ au for M2, and $R_{rim} = 0.3–0.47$ au for M3. Over the extent of the rim, the temperature structure is close to vertically isothermal, suggesting that stellar irradiation dominates. Farther out, the temperature along vertical cuts increases toward the midplane owing to the viscous heating (see also Schobert et al. 2019). When moving from model M1 to M2 and then to M3, we see that the temperature over the rim extent becomes gradually larger, leading to a warmer and warmer infrared photosphere in this region. Similarly, the whole interior of the flaring disk region in models M2 and M3 is puffed up. The main cause for the difference between models M1 and M2 is the boundary condition (Section 2.2, Appendix A), which reduces the cooling efficiency of model M2. The temperature increase found in model M3 is due to the strong frontal irradiation of the wall-like rim and the radial radiation diffusion. It is important to point out that while the classical wall-like rim of Dullemond et al. (2001) has a shadowed region right outside the wall, our model M3 avoids that due to the reduced cooling efficiency combined with viscous heating (see Appendix B where the role of viscous heating is further discussed).

Next, Figure 4 compares several characteristic radial profiles of models M1, M2, and M3. Focusing on the midplane temperature first, we can see that the models differ mostly in the radial range between the dust halo and the outer disk, $r \approx 0.3–0.6$ au. At $r < 0.3$ au, the disk reaches the optically thin temperature while at $r > 0.6$ au, the viscous heating dominates
<figure><img src="image_6.png" /></figure>

**Figure 3.** Temperature profile (represented by filled contours) in the meridional plane of the disk models M1 (top), M2 (middle), and M3 (bottom). Individual curves show where: the dust grains start to condense (solid gray), the dust is fully condensed (solid black), the radially integrated optical depth to stellar irradiation becomes unity (solid white), and the vertically integrated optical depth to infrared emission becomes unity (dashed green). Blue arrows delimit the radial range of the inner disk rim. As we plot $z/r$ on the vertical axis (where z is the height above the midplane and r is the spherical radius), we point out that radially propagating irradiating rays would appear as horizontal lines.

in the midplane, making its temperature independent of the boundary condition. Despite the relative match in the midplane, however, temperature differences do appear in upper disk layers, as already shown in the context of Figure 3 and highlighted in the middle panel of Figure 4, which depicts the temperature profile of the IR photosphere. Clearly, model M3 has the hottest photosphere in the rim region, model M2 is intermediate, and model M1 is the coldest.

## 3.2. Half-light Radii and Infrared Fluxes

Figure 5 compares $R_{hl}$ and $\mathcal{F}$ (see Section 2.3.2 for definitions) derived from our models with observations. It is necessary to point out that the observational data shown in Figure 5 are secondary interferometric quantities in the sense

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

**Figure 4.** Radial profiles of the midplane temperature (top), IR photospheric temperature (middle), and gas surface density $\Sigma$ (bottom) for models M1 (blue curve), M2 (red curve), and M3 (black curve). The top panel distinguishes four characteristic regions of dusty inner rims (see Flock et al. 2016; Ueda et al. 2017, for details).

that they are based on parametric brightness distributions fitted to the interferometric measurements and are therefore model dependent (for an opposing view, however, see Section 2.3.3). The purpose of the comparison here is simply to get a qualitative understanding of how changing various components of our physical model affects the half-light radii and the contribution of the rim to the overall flux. Generally, we see that the increase in $R_{hl}$ and $\mathcal{F}$ follows the increase in the photospheric temperature identified between individual models in the previous Section 3.1.

Starting with $R_{hl}$ (top panel in Figure 5), we can see that all our models are roughly consistent with previously reported values at NIR wavelengths. For model M1, however, $R_{hl}$ increases with $\lambda$ rather weakly, and thus the half-light radius does not grow enough to match the observations at MIR wavelengths. Our models M2 and M3, on the other hand, both exhibit a steepening toward the N band and they seem to match the observations quite well, with model M2 being slightly nearer the data points. It seems that the boost of $R_{hl}$ is mostly driven by the warm boundary because both models M2 and M3 use it and their $R_{hl}$ profiles are quite similar.

Focusing on $\mathcal{F}$ at $\lambda < 5$ µm (bottom panel in Figure 5), we can see that model M1 has the weakest contribution to the flux. By adding the warm boundary (in model M2), the disk flux
<figure><img src="image_8.png" /></figure>

**Figure 5.** Half-light radius (top) and fractional disk flux (bottom) as functions of the wavelength. We show synthetic data corresponding to models M1 (blue), M2 (red), and M3 (black), as well as data corresponding to previous morphological studies of the inner disk emission (black points and error bars; Lazareff et al. 2017; GRAVITY Collaboration et al. 2019; Varga et al. 2021). The wavelength range at the bottom panel is limited to $\lambda < 5$ µm in order to highlight the differences at shorter wavelengths (the dependence at larger wavelengths is asymptotic and therefore less interesting).

increases, but only weakly near $\lambda \sim 1.5$ µm. In model M3, there is yet another flux increase, most prominent at short wavelengths. Therefore, the boost of $F$ at very short wavelengths $\lambda \sim 1.5$ µm can be achieved when the geometry of the sublimation front becomes wall-shaped (because models M1 and M2 have rounded rims and their fractional disk flux for $\lambda \sim 1.5$ µm is nearly the same).

### 3.3. Visibilities

In Figure 6, we plot the visibility amplitude $|V|$ as a function of the deprojected spatial frequency $B_{\text{eff}}/\lambda$ to remove the effect of object inclination. The data points show realistic measurements obtained with VLTI in the bands H, K, and L during 2019, while the colored areas correspond to the visibility profiles of our models, with the boundary curves calculated at the minimum (top boundary of each colored area) and maximum (bottom boundary of each colored area) wavelengths of each band.

First, we notice that the midpoint of the visibility curves (between $|V|=1$ and the first bend; see Section 2.3.3) at a given band does not seem to strongly depend on the model choice, which tells us that the half-light radii up to the L band do not differ very much between different models. This is consistent with what is shown in the inset of Figure 5.

Next, as the disk gets warmer (M1 → M3), the visibility curves start to decay more steeply and they begin to level off at

lower $|V|$. This reflects the fact that $F_{\text{star}}/F_{\text{tot}}$ decreases as the fractional disk flux increases,<fnref n="12" /> as also shown in Figure 5. This change is the most prominent in the H band. The fact that the rim has a torus-like brightness distribution with a sharp edge (Section 3.4) and that it represents a resolved source leads to bounces of the visibility curve, especially at long baselines. While we cannot say much about the presence or absence of bounces in the displayed K- and L-band observations, they are clearly absent in the H-band observations, which confirms earlier works (Benisty et al. 2010; Setterholm et al. 2018) attributing this mismatch to a presence of a smooth emission source filling the region inside the rim in HD 163296 (e.g., the gas continuum or superrefractory dust species inwards from the sublimation radius).

To provide a simple quantitative comparison between the models and the observations, we counted the number of observational data points enclosed<fnref n="13" /> between the model curves of each specific band. We found that model M1 matches 2%, 54%, and 79% of H-band, K-band, and L-band observations, respectively. As for model M2, we found an overlap with 15%, 66%, and 60% of H-band, K-band, and L-band observations, respectively. Finally, model M3 is consistent with 48%, 42%, and 66% of observations in bands H, K, and L, respectively. On average, model M3 leads to visibility curves closest to the observations, although model M1 is better when focusing on the L band alone and model M2 outperforms the others in the K band.

### 3.4. Infrared Emission in Detail

Let us now examine the synthetic images<fnref n="14" /> themselves (Figure 7) and explore the spatial distribution of infrared emission. Overall, all synthetic images exhibit a dominant torus-like emission, with the torus being sharply truncated at the inner edge of the rim (the blue arrows can guide the eye as they correspond to the rim extent shown in Figure 3). However, it is important to point out that the width of the brightest part of the torus is smaller than 1 mas, which is roughly the best possible resolution that the VLTI can reach with its most detailed H-band observations. The images shown here are therefore highly idealized.

Comparing models M1 (top row) and M2 (middle row) first, we can see that their H-band infrared emission (left column at $\lambda = 1.75$ µm) is largely similar, despite the differences in the radial range of the rim. This is because most of the H-band emission comes from the very tip of the rim, which exhibits a similar grazing angle with respect to the incoming irradiation (Figure 3) and also a similar profile of the infrared photosphere at $r \approx 0.3-0.35$ au (Figure 4). This is consistent with the similarity of the H-band fractional disk flux and visibility curves of these two models. On the other hand, the L-band emission (right column at $\lambda = 3.5$ µm) is more radially extended for model M2, for which it covers roughly the entire radial extent of the rim. The additional emission of model M2 (dark-red-colored) compared to model M1 is relatively weak but covers a large enough surface area to increase the L-band flux of model M2 as seen in Figure 5. The cause of this

<footnote n="12">¹² We refer an interested reader to Benisty et al. (2010), Dullemond & Monnier (2010), and Lazareff et al. (2017), where similar interpretations of the rim-induced visibility curves were given.</footnote>

<footnote n="13">¹³ We point out that we take the observational errors of $V^2$ into account.</footnote>

<footnote n="14">¹⁴ Synthetic images of our main models M1, M2, M3, and M3Fe are freely available as FITS files at https://sirrah.troja.mff.cuni.cz/~chrenko/hd163296/.</footnote>
<figure><img src="image_9.png" /></figure>

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

<figure><img src="image_11.png" /><figcaption>**Figure 6.** Visibility $|V|$ as a function of the deprojected spatial frequency. Colored areas represent bundles of synthetic visibility curves in bands H (blue), K (green), and L (red) that are derived from our models M1 (top), M2 (middle), and M3 (bottom). Data points show the VLTI observations from the 2019 epoch taken with the PIONIER (circles), GRAVITY (triangles), and MATISSE (crosses) instruments. The wavelength intervals of synthetic visibilities and observations are the same and are given in the plot legend.</figcaption></figure>

difference can be traced back to Figure 3 where the transition from the tip of the inner rim to the flaring outer disk is more abrupt for model M1 than it is for model M2. Model M2,

instead, has a noticeable transitional region that is less exposed to stellar irradiation than the tip of the rim but more exposed than the outer flaring disk. Similar geometrical differences can be noticed in the infrared photosphere of model M2, as well as a temperature bump between ≈0.35 and 0.6 au in the middle panel of Figure 4.

Synthetic images of model M3 (bottom row of Figure 7) exhibit the largest absolute intensity compared to models M1 and M2. The bright central part matches the frontally irradiated wall-like sublimation front, viewed under the inclination angle of the disk. The overall larger intensity, which also translates to larger fluxes and the previously discussed shifts in the visibility profiles, is yet another manifestation of the hot infrared photosphere.

### 3.5. Maximizing the NIR Flux

Previous findings related to models M1–M3 indicate a close connection between the NIR flux and the rim temperature. In this section, we explore whether it is possible to increase the flux even more by simply considering a larger sublimation temperature of dust grains, thus making the rim hotter (see also Klarmann et al. 2017). To demonstrate this possibility, we computed one additional variation of model M3, designated M3Fe, assuming $T_s = 1550$ K.

Increasing $T_s$ alone while keeping the other model components fixed would allow the dust grains to survive closer to the star and the entire disk rim would shift inwards, making its $R_{hl}$ inconsistent with the observations. To keep $R_{hl}$ comparable to the dependence discussed in Figure 5, it is necessary to modify the dust opacities. We remind the reader that the ratio

$$ \varepsilon = \frac{\kappa_d(T_s)}{\kappa_d(T_*)} \qquad (21) $$

determines the optically thin temperature of isolated dust grains

$$ T_{\text{thin}} = \varepsilon^{-1/4} \sqrt{\frac{R_{\star}}{2r}} T_{\star}, \qquad (22) $$

as well as the radius where dust grains condense in the midplane (Monnier & Millan-Gabet 2002; Ueda et al. 2017):

$$ R_{\text{rim}}^{\text{in}} = \frac{1}{2} R_{\star} \varepsilon^{-1/2} \left( \frac{T_{\star}}{T_{\text{halo}}} \right)^2, \qquad (23) $$

where $T_{\text{halo}}$ is the temperature in the halo of the rim.

After testing several dust compositions, we found that $R_{hl}$ can be preserved when dust grains composed of solid metallic iron (Henning & Stognienko 1996; Woitke et al. 2018) are considered.<fnref n="15" /> with a size distribution $f(a) \propto a^{-3.5}$ ranging from 0.05 to 1 µm. Then, $\kappa_d(T_s) = 873$ and $\kappa_d(T_*) = 4187 \text{ cm}^2 \text{ g}^{-1}$, yielding $\varepsilon \approx 1/5$ and $T_{\text{halo}} \approx 1700$ K (the latter was found directly from our simulation). The sublimation temperature of metallic iron sensitively depends on the local chemical conditions but can in principle reach the assumed value $T_s = 1550$ K (e.g., Brož et al. 2021). We also point out that recent laboratory experiments (Bogdan et al. 2023) show that metallic iron efficiently and “automatically” forms from

<footnote n="15">¹⁵ Our wavelength-dependent opacities of metallic iron can be reproduced with the OPTOOL code as follows: optool iron 1.0 -a 0.05 1.0 3.5 -mie.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 7.** Synthetic images (emission intensity maps) derived from models M1 (top), M2 (middle), and M3 (bottom) at wavelengths $\lambda = 1.75$ and 3.5 µm (left and right column, respectively). We point out that the intensity is shown in the linear scale and the range of color bars slightly differs between the left and right column. The emission of the central star is masked. The blue arrow marks the approximate radial range of the disk rim defined in Figure 3. The white scale bar corresponds to 0.5 au.</figcaption></figure>

silicates at $T > 1200$ K, and thus it can indeed be present at the inner disk rim at large abundances.

The results for model M3Fe are presented in Figure 8. The temperature distribution in the meridional plane (top panel) is similar to model M3 but the halo and the rim itself are hotter, while the flared outer disk is puffed up even more. Looking at

the fractional disk flux (middle panel), it is clear that the model now matches observations even at the shortest H- and K-band wavelengths. As for the visibility trend (bottom panel), we can see that model M3Fe is a logical continuation of the sequence of panels shown in Figure 6: the visibility curves reflect the increase in the disk flux and so they decrease more steeply. The
<figure><img src="image_13.png" /><figcaption>**Figure 8.** Summary of results for model M3Fe. Top: temperature distribution in the meridional plane of the disk (compare with Figure 3). Middle: fractional disk flux as a function of the wavelength (compare with Figure 5). Bottom: synthetic visibility curves compared to observations (compare with Figure 6).</figcaption></figure>

H-band synthetic curves are now positioned partially below the set of observations. The K-band and L-band synthetic curves continue to depart from the real data, overlapping 27% and 40% of observations, respectively (the fewest of all models). The bounce at long wavelengths is the smallest when compared to models M1–M3, yet it is still present.

## 4. Conclusions

Infrared emission of the protoplanetary disk HD 163296 can possibly arise from the sublimation front of dust grains, known as the inner rim, as well as from the dust-free region interior to the rim (e.g., Tannirkulam et al. 2008; Benisty et al. 2010; Setterholm et al. 2018; GRAVITY Collaboration et al. 2021;

Varga et al. 2021). In this paper, our strategy was to calculate various physical models of the inner rim in order to assess how they compare to interferometric observables. We used radiative hydrostatic modeling (following Flock et al. 2016) to derive the structure of the inner rim, we calculated synthetic images of the NIR and MIR emission, and we compared the half-light radii $R_{hl}$, fractional disk fluxes $\mathcal{F} = F_{\text{disk}}/F_{\text{tot}}$, and interferometric visibilities $V$ with the VLTI multiband data (e.g., Lazareff et al. 2017; GRAVITY Collaboration et al. 2019; Varga et al. 2021). Of the three quantities, $V$ is the most and $\mathcal{F}$ is the least robust. Interestingly, we found theoretical arguments for $R_{hl}$ to be more robust than previously thought, as discussed in Section 2.3.3.

In our set of models, we started from a nominal model (M1) and we gradually increased the temperature of the infrared photosphere near the inner rim by reducing the cooling efficiency of the disk (model M2), changing the rim geometry from rounded to wall-like (model M3), and allowing the dust grains to survive at larger temperatures (model M3Fe). We concluded that model M3 is the one closest to observations because it provides the best match to the visibility curves and it reproduces previously reported $R_{hl}$. It can also match $\mathcal{F}$ fairly well, with the exception of bands H and K for which earlier morphological fits (Lazareff et al. 2017; GRAVITY Collaboration et al. 2019) predict fractional disk contributions larger by a factor of 1.4 and 1.1, respectively, compared to our model.

However, we pointed out that our model visibility curves always exhibit a bounce at long baselines, most notably in the H band, due to the fact that the rim emits as a narrow torus. Such bounce is not observed, confirming the need for an additional emission component in the disk (e.g., Benisty et al. 2010; Setterholm et al. 2018). Additionally, matching the visibility curves across multiple bands with a rim model alone is clearly challenging because even though model M3 provides the best match on average, other models outperform it when focusing on single bands (e.g., model M1 is better in the L band; model M2 is better in the K band). In other words, when one tries to modify the physical model to improve the match in a single band, the match in other bands might actually become worse.

In model M3, the reduced cooling efficiency was achieved by setting a warmer boundary condition for the escape of photons by radiative diffusion, and the wall-like shape was obtained owing to the uniform sublimation temperature of dust grains. The realism of both model ingredients is debatable. The reduced cooling efficiency was used in similar forms in the majority of recent works oriented on the inner rim (Flock et al. 2016; Schobert et al. 2019) but we pointed out (Appendix A) that it leads to temperature profiles warmer than, and therefore inconsistent with, those resulting from frequency-dependent Monte Carlo calculations. The uniform sublimation temperature, on the other hand, used to be assumed in the classical inner rim models (e.g., Dullemond et al. 2001) but was later abandoned as the local conditions in terms of vapor densities at saturation pressures are expected to change with height above the disk midplane (Isella & Natta 2005; Kama et al. 2009). Nevertheless, even if the reduced cooling efficiency or the uniform sublimation temperature turns out to be inadequate, they might still point to the correct disk structure, which future models could strive to reproduce by considering additional physical processes, for instance, thermo-chemical effects or high-energy nonequilibrium heating of the gas atmosphere.
To summarize the physical features of model M3, the rim extends between 0.3 and 0.47 au, with $R_{hl} \simeq (0.28, 0.29, 0.31, 0.97)$au and $\mathcal{F} \simeq (0.32, 0.75, 0.94, 1.0)$ at $\lambda = (1.5, 2.2, 3.5, 10.5)$ µm, respectively.

The impact and future applicability of our study can be threefold. First, while parametric morphological fitting is by far the most common approach to interpret sparse interferometric observations of sub-au disk regions, it is rarely done in a multiband manner and if so, it often lacks a link to physical models. Therefore, synthetic data from physical models could be used to calibrate morphological fits (by checking if the fit can retrieve the important features of the physical model) before applying those to real observations. Second, the emission component interior to the rim in HD 163296 has not yet been described with a physical model. When such a description becomes available, our study can help tweak the relative contribution between the rim and the additional interior source to the overall signal. Third, radiative hydrostatic models of the inner rim are notoriously known for not producing enough NIR flux (Vinković et al. 2006; Turner et al. 2014; Flock et al. 2016) and our findings provide ways how to boost the flux when needed (note that model M3Fe matches the observed $\mathcal{F}$ very well).

## Acknowledgments

Based on observations collected at the European Southern Observatory under ESO program(s) 0101.C-0896(D), 0103.C-0347(A,C), 0103.C-0915(A,B,C), 0103.D-0153(C,D), and 0103.D-0294(A). This work was supported by the Czech Science Foundation (grant 21-23067M), the Charles University Research Centre program (No. UNCE/24/SCI/005), and the Ministry of Education, Youth and Sports of the Czech Republic through the e-INFRA CZ (ID:90254). M.F. acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 757957). T.U. acknowledges the support of the DFG-Grant "INSIDE: The INner regions of protoplanetary disks: Simulations and observations" (project number 465962023). M.B. has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 101002188). We sincerely thank Kees Dullemond and Roy van Boekel for motivating discussions, and Jozsef Varga for sharing the data from Varga et al. (2021) with us. We wish to thank the two anonymous referees whose valuable comments allowed us to greatly improve this paper.

**Facility:** VLTI (PIONIER, GRAVITY, MATISSE)

**Software:** ASTROPY (Astropy Collaboration et al. 2013, 2018), FARGO3D (Benítez-Llambay & Masset 2016), MATPLOTLIB (Hunter 2007), NUMPY (Harris et al. 2020), OPTOOL (Dominik et al. 2021), PMOIRED (Mérand 2022), RADIATION_CODE (Schobert et al. 2019), RADMC-3D (Dullemond et al. 2012), RADMC3DPY (https://www.ita.uni-heidelberg.de/dullemond/software/radmc-3d/manual_rmcpy).

## Appendix A

### Testing Boundary Conditions for the Radiation Energy Density

The aim of this appendix is to exemplify the influence of boundary conditions for the radiation energy density (Section 2.2) on the temperature profile of the inner disk, as

<figure><img src="image_14.png" /><figcaption>**Figure 9.** Radial profiles of the midplane temperature based on the model S100 of Flock et al. (2016). The individual panels demonstrate the influence of the boundary conditions for the radiation energy density (Section 2.2)—the warm boundary is used in the top panel while the cold boundary is used in the bottom panel. Individual curves are labeled in the plot legend and described in detail in Appendix A. To guide the eye, we advise the reader to mainly follow the changes of the blue and green curves (black and gray curves remain fixed, while the red curve changes only slightly).</figcaption></figure>

well as to provide a link to previous studies. To do so, we followed Flock et al. (2016) and tried to reproduce their model designated S100, in which the surface density of the gas is nonevolving and uniform, $\Sigma(r) = 100 \text{ g cm}^{-2}$.

The results are summarized in Figure 9. Let us first explain the meaning of individual curves. The gray dashed curves show the optically thin temperature of gas for reference, $T_{\text{thin}} = \sqrt{R_{\star}/(2r)} T_{\star}$. The black curves show the original result of Flock et al. (2016). The blue curves are the results of calculations with the code used in this study; the top panel corresponds to the warm boundary, while the bottom panel shows the result for the cold boundary. In a similar fashion, the green curves show the results of calculations that we performed using the RADIATION_CODE of Schobert et al. (2019), the aim being to provide an independent sanity check. Finally, the red curves show temperature profiles obtained with thermal frequency-dependent Monte Carlo simulations using RADMC-3D, in which the input gas and dust densities were taken from our hydrostatic calculations (those represented by blue curves).

Ideally, all solid curves should overlap. However, we see that this is only true in the innermost dust-free disk and in the adjacent dusty halo. The midplane temperature of the optically thick regions ($r \gtrsim 0.45$ au) exhibits differences. When the boundary condition is warm (or see Schobert et al. 2019, for their default boundary condition), the temperature profiles based on our hydrostatic calculations match that of Flock et al. (2016). However, the thermal Monte Carlo calculation with RADMC-3D leads to substantially lower temperatures across the
disk rim and the outer disk, even though the gas and dust density is directly adopted from the hydrostatic calculations.

If, on the other hand, the cold boundary condition is used, there is an agreement between our hydrostatic calculations and the thermal Monte Carlo run, but all these temperature profiles depart from that of Flock et al. (2016). For completeness, we point out that the waves at $r \gtrsim 1$ au in the bottom panel of Figure 9 are manifestations of the irradiation instability (e.g., Watanabe & Lin 2008; Wu & Lithwick 2021; Melon Fuksman & Klahr 2022).

It is difficult to assess which of these boundary conditions is more realistic. In general, it seems that the cold boundary is more common in models with radiative diffusion (e.g., Bitsch et al. 2013; Kolb et al. 2013) and it also leads to a better match with the Monte Carlo multifrequency approach, which is physically superior to the simple radiative diffusion. However, it is still instructive not to disregard the warm boundary because, as shown in the main text, it can sometimes lead to a better match with observations, thus laying valuable ground-work for future studies.

# Appendix B

## On the Role of Viscous Heating and Convergence of the Hydrostatic Method

Models presented in our work include a simplified treatment for viscous heating, as explained in Section 2.1.4. In terms of optically thick regions of the disk, viscous heating is implemented correctly in the dead zone of the disk but its influence is underestimated at the very tip of the rim (because we use $\alpha_{DZ}$ to evaluate the viscous heating term). Nevertheless, it plays an important role in the stability and convergence of our models, especially those employing the cold boundary condition for $E_R$ (see Section 2.2). To illustrate that, the left panel of Figure 10 shows a variation of model M1 computed without viscous heating, meaning that the disk is only passively irradiated. Without viscous heating, the disk region outward from the rim falls into a shadow and progressively becomes colder and colder. The reason is that at one point during the iterative sequence, the surface density outwards from the rim becomes so large that vertical cooling starts to act more efficiently than heating by radial radiative diffusion. As no other heat source is operating in this region (due to the shadowing by the rim), the temperature slightly decreases, which means that the local viscosity $\nu$ (see Section 2.1.1)

decreases as well. Because we impose constant $\dot{M}$ through the disk, Equation (3) dictates that the local surface density $\Sigma$ has to increase to compensate for lower $\nu$, leading to a feedback loop that creates a cold spot visible in Figure 10. Models that “fall” into this loop cannot be reliably converged and it is questionable whether they are physically realistic (we expect that the local density peak $\Sigma \sim 3 \times 10^4$ g cm$^{-2}$ overlapping with the cold spot would become Rossby-unstable in a hydrodynamic run). Viscous heating, even in our simplified form, helps circumvent such convergence issues.

Our models discussed in the main body of the paper are well converged, reaching the relative change in the temperature of the order of $10^{-5}$ during the last iterations. The right panel of Figure 10 shows the evolution of several characteristic surfaces in model M2 during selected iterations (iteration number 200 is the last one). We can see that the front of full dust evaporation, which is the inner boundary of the dusty halo, has converged already after ~100 iterations as the solid blue curve is hidden underneath the solid black curve. The IR photosphere is converged after ~150 iterations because the dashed purple curve is indistinguishable from the dashed black one (the relative difference between the two curves is below $10^{-3}$).

Finally, let us verify our assumptions concerning the characteristic timescales in the disk (Section 2.1). We focus on model M2 and the characteristic radial distance $r_c = 0.4$ au, roughly corresponding to the middle of the rim extent. The dynamical timescale is $t_{dyn} = 1/\Omega_K(r_c) \simeq 10^6$ s and the viscous timescale is $t_{vis} = 1/(\alpha_{MR}h^2\Omega_K(r_c)) \simeq 10^{10}$ s (with the local aspect ratio $h = 0.027$). The timescale of thermal relaxation is determined by the radiative diffusion in the vertical direction from the midplane toward the infrared photosphere. The thermal diffusivity due to radiation in the optically thick limit is (e.g., Lin & Youdin 2015; Jiménez & Masset 2017)

$$ \chi_t = \frac{16\sigma T^3}{3\rho^2 c_v \kappa_d(T_s) f_{d2g,max}}, \qquad (B1) $$

leading to $\chi_t \simeq 1.7 \times 10^{16}$ cm$^2$ s$^{-1}$ for $T \simeq 860$ K and $\rho \simeq 10^{-9}$ g cm$^{-3}$. Taking the local height of the infrared photosphere $H_{IR}(r_c) \simeq 2.4H \simeq 0.026$ au, we obtain $t_{rad} = H_{IR}^2/\chi_t \simeq 10^7$ s. Our estimates yield the inequality $t_{dyn} < t_{rad} \ll t_{vis}$, which is consistent with Flock et al. (2016).

<figure><img src="image_15.png" /><figcaption>Figure 10. Left: Temperature map of model M1 with neglected viscous heating. There is a shadowed region outwards from the rim and an apparent cold spot. Right: Evolution of surfaces of full dust evaporation (solid curve) and IR emission (dashed curve) during the full iteration cycle of model M2.</figcaption></figure>
ORCID iDs

Ondřej Chrenko https://orcid.org/0000-0001-7215-5026

Mario Flock https://orcid.org/0000-0002-9298-3029

Takahiro Ueda https://orcid.org/0000-0003-4902-222X

Antoine Mérand https://orcid.org/0000-0003-2125-0183

Myriam Benisty https://orcid.org/0000-0002-7695-7605

Raúl O. Chametla https://orcid.org/0000-0002-4327-7857

References

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

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

Bans, A., & Königl, A. 2012, ApJ, 758, 100

Benisty, M., Natta, A., Isella, A., et al. 2010, A&A, 511, A74

Benítez-Llambay, P., & Masset, F. S. 2016, ApJS, 223, 11

Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148

Bitsch, B., Crida, A., Morbidelli, A., Kley, W., & Dobbs-Dixon, I. 2013, A&A, 549, A124

Bogdan, T., Pillich, C., Landers, J., Wende, H., & Wurm, G. 2023, A&A, 670, A6

Bohlin, R. C., Mészáros, S., Fleming, S. W., et al. 2017, AJ, 153, 234

Brož, M., Chrenko, O., Nesvorný, D., & Dauphas, N. 2021, NatAs, 5, 898

Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368

Chrenko, O., Chametla, R. O., Nesvorný, D., & Flock, M. 2022, A&A, 666, A63

Chrenko, O., & Lambrechts, M. 2019, A&A, 626, A109

Chrenko, O., & Nesvorný, D. 2020, A&A, 642, A219

D'Angelo, G., & Bodenheimer, P. 2013, ApJ, 778, 77

Dobbs-Dixon, I., Cumming, A., & Lin, D. N. C. 2010, ApJ, 710, 1395

Dominik, C., Min, M., & Tazaki, R., 2021 OpTool: Command-line driven tool for creating complex dust opacities, Astrophysics Source Code Library, ascl:2104.010

Draine, B. T. 2003, ApJ, 598, 1017

Dullemond, C. P., Dominik, C., & Natta, A. 2001, ApJ, 560, 957

Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, RADMC-3D: A multi-purpose Radiative Transfer Tool, Astrophysics Source Code Library, ascl:1202.015

Dullemond, C. P., & Monnier, J. D. 2010, ARA&A, 48, 205

Dzyurkevich, N., Flock, M., Turner, N. J., Klahr, H., & Henning, T. 2010, A&A, 515, A70

Eisner, J. A., Chiang, E. I., Lane, B. F., & Akeson, R. L. 2007, ApJ, 657, 347

Fairlamb, J. R., Oudmaijer, R. D., Mendigutía, I., Ilee, J. D., & van den Ancker, M. E. 2015, MNRAS, 453, 976

Flock, M., Fromang, S., Turner, N. J., & Benisty, M. 2016, ApJ, 827, 144

Flock, M., Turner, N. J., Mulders, G. D., et al. 2019, A&A, 630, A147

GRAVITY Collaboration, Perraut, K., Labadie, L., et al. 2019, A&A, 632, A53

GRAVITY Collaboration, Sanchez-Bermudez, J., Caratti O Garatti, A., et al. 2021, A&A, 654, A97

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Henning, T., & Stognienko, R. 1996, A&A, 311, 291

Hillenbrand, L. A., Strom, S. E., Vrba, F. J., & Keene, J. 1992, ApJ, 397, 613

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

Isella, A., & Natta, A. 2005, A&A, 438, 899

Jankovic, M. R., Mohanty, S., Owen, J. E., & Tan, J. C. 2022, MNRAS, 509, 5974

Jankovic, M. R., Owen, J. E., Mohanty, S., & Tan, J. C. 2021, MNRAS, 504, 280

Jiménez, M. A., & Masset, F. S. 2017, MNRAS, 471, 4917

Kama, M., Min, M., & Dominik, C. 2009, A&A, 506, 1199

Klarmann, L. 2018, in Take a Closer Look (München: ESO), 84

Klarmann, L., Benisty, M., Min, M., et al. 2017, A&A, 599, A80

Kley, W. 1989, A&A, 208, 98

Kluska, J., Berger, J. P., Malbet, F., et al. 2020, A&A, 636, A116

Kobus, J., Wolf, S., Ratzka, T., & Brunngräber, R. 2020, A&A, 642, A104

Kolb, S. M., Stute, M., Kley, W., & Mignone, A. 2013, A&A, 559, A80

Kraus, S., Hofmann, K. H., Benisty, M., et al. 2008, A&A, 489, 1157

Lada, C. J., & Adams, F. C. 1992, ApJ, 393, 278

Lazareff, B., Berger, J. P., Kluska, J., et al. 2017, A&A, 599, A85

Leinert, C., van Boekel, R., Waters, L. B. F. M., et al. 2004, A&A, 423, 537

Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321

Lin, M.-K., & Youdin, A. N. 2015, ApJ, 811, 17

Masset, F. S., Morbidelli, A., Crida, A., & Ferreira, J. 2006, ApJ, 642, 478

Melon Fuksman, J. D., & Klahr, H. 2022, ApJ, 936, 16

Mérand, A. 2022, Proc. SPIE, 12183, 121831N

Mészáros, S., Allende Prieto, C., Edvardsson, B., et al. 2012, AJ, 144, 120

Millan-Gabet, R., Schloerb, F. P., & Traub, W. A. 2001, ApJ, 546, 358

Monnier, J. D., & Millan-Gabet, R. 2002, ApJ, 579, 694

Monnier, J. D., Millan-Gabet, R., Billmeier, R., et al. 2005, ApJ, 624, 832

Mulders, G. D., Pascucci, I., Apai, D., & Ciesla, F. J. 2018, AJ, 156, 24

Natta, A., Prusti, T., Neri, R., et al. 2001, A&A, 371, 186

Petigura, E. A., Marcy, G. W., Winn, J. N., et al. 2018, AJ, 155, 89

Pollack, J. B., Hollenbach, D., Beckwith, S., et al. 1994, ApJ, 421, 615

Preibisch, T., Ossenkopf, V., Yorke, H. W., & Henning, T. 1993, A&A, 279, 577

Schobert, B. N., Peeters, A. G., & Rath, F. 2019, ApJ, 881, 56

Setterholm, B. R., Monnier, J. D., Davies, C. L., et al. 2018, ApJ, 869, 164

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 500, 33

Tannirkulam, A., Monnier, J. D., Millan-Gabet, R., et al. 2008, ApJL, 677, L51

Tilling, I., Woitke, P., Meeus, G., et al. 2012, A&A, 538, A20

Turner, N. J., Benisty, M., Dullemond, C. P., & Hirose, S. 2014, ApJ, 780, 42

Tuthill, P. G., Monnier, J. D., & Danchi, W. C. 2001, Natur, 409, 1012

Ueda, T., Flock, M., & Okuzumi, S. 2019, ApJ, 871, 10

Ueda, T., Okuzumi, S., & Flock, M. 2017, ApJ, 843, 49

Varga, J., Hogerheijde, M., van Boekel, R., et al. 2021, A&A, 647, A56

Varnière, P., & Tagger, M. 2006, A&A, 446, L13

Vinković, D., Ivezić, Ž., Jurkić, T., & Elitzur, M. 2006, ApJ, 636, 348

Watanabe, S.-i., & Lin, D. N. C. 2008, ApJ, 672, 1183

Wichittanakom, C., Oudmaijer, R. D., Fairlamb, J. R., et al. 2020, MNRAS, 493, 234

Woitke, P., Helling, C., Hunter, G. H., et al. 2018, A&A, 614, A1

Wu, Y., & Lithwick, Y. 2021, ApJ, 923, 123