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

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 F in the wavelength range 1.5–13 μ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 ( F≃0.25 at 1.5 μ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.


INTRODUCTION
Sub-au regions of protoplanetary disks represent the environment that has shaped precursors of terrestrial planets as well as the numerous population of shortperiod 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 ≈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 radius1 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 NIR flux, and reduce the closure phase signal.Among possible explanations of the additional emission component are 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;Varga et al. 2021;GRAVITY Collaboration et al. 2021), possibly pointing to a 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, re-confirming 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. 2016Flock et al. , 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 (since we do not model the emission component inside the rim), nor describe the temporal variability of the inner disk asymmetry (since our model is by construction static and symmetric).It is rather to set groundwork for followup studies to help to distinguish how much the inner rim can contribute to the interferometric signals.In future, our models can be readily combined with morphological fitting (e.g. by parametrizing an azimuthal asymmetry on top of one of our base models), or they can help tweaking 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 the 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.

Radiative hydrostatic disk structure
We use the radiative hydrostatic approach of Flock et al. (2016Flock et al. ( , 2019) ) to calculate the distribution of the disk gas density ρ, dust density ρ 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 ρ in hydrostatic equilibrium (see Section 2.1.1 for constraints and details).
2. Perform ten sub-iterations of: (a) Radially integrated optical depths to stellar irradiation τ (which depend on the dust-togas ratio f d2g that sets the local optical depth of dust in each grid cell; see Section 2.1.2and Equation 5), (b) Dust-to-gas ratio f d2g (which has to depend on τ in order to properly resolve irradiation absorption; see Section 2.1.3and Equation 6) 3. Keeping ρ, τ , and f d2g fixed, evolve twotemperature 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).Note-Stellar parameters of HD 163296 are adopted from Wichittanakom et al. (2020) and the mass accretion rate is from Fairlamb et al. (2015).
4. Return to the beginning with the new temperature field.

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 where P is the thermal pressure, ϕ is the colatitude, r is the radius, G is the gravitational constant, and M ⋆ is the mass of the central star.The azimuthal dimension θ 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 where γ is the adiabatic index, ϵ 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 Ṁ is uniform.Then the equation where ν is the effective viscosity, provides a constraint on the gas surface density Σ.The viscosity is parametrized via the Shakura & Sunyaev (1973) prescription ν = αc 2 s /Ω K , where c s = γP/ρ is the adiabatic sound speed and Ω K is the Keplerian angular frequency.We point out that a density-weighted vertical average of ν is used when evaluating Equation (3) (see Chrenko & Nesvorný 2020).Furthermore, our α-parametrization 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 ) where T MRI is the transition temperature, α DZ and α MRI are α-viscosities in the dead and active zones, respectively.Finally, we assume that Σ is related to the volume density in the midplane ρ mid = Σ/( √ 2πH) as if the disk was vertically isothermal, with H = c s /( √ γΩ K ) being the pressure scale height (H is evaluated from c s in the midplane for the purpose of estimating ρ mid ).
Equipped with the aforementioned closure relations, it is possible to reconstruct the radial profile of ρ mid (r) for a fixed temperature field.This serves as a starting point for solving Equation (1) and thus finding ρ(r, ϕ) throughout the rest of the disk in each iteration (see Flock et al. 2016;Chrenko & Nesvorný 2020).

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 (κ P ) and Rosseland (κ R ) opacities are the same.We define the gas opacity κ gas , the dust opacity to its own thermal emission κ d (T s ), and the dust opacity to stellar irradiation κ 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 where τ 0 is the optical depth inwards from our computational grid, at r < r in (see Flock et al. 2016).We point out that ρ d does not necessarily represent the total dust content but mainly accounts for small grains, which are the dominant opacity contributors.
As for the actual value of κ gas = 10 −5 cm 2 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 κ d (T s ) and κ d (T ⋆ ), we first calculated wavelength-dependendent 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) ∝ a −3.5 ranging between 3 × 10 −3 and 10 2 µm.Using the optool code 2 (Dominik et al. 2021), we obtained the wavelength-dependent opacities shown in the top panel of Figure 1.Subsequently, we calculated κ d (T s ) as the Planck opacity at 1400 K (which is a proxy of the temperature in the dusty disk halo) and κ 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.

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 f d2g (r, ϕ) = ρ d /ρ calculated as (Flock et al. 2019) otherwise.
(6) where f 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 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-µm-sized grains (Birnstiel et al. 2012).
To prevent numerical problems, we also define a floor value f d2g > f d2g,min .Additionally, the stability of our method is improved by ramping f d2g,max from f d2g,min up to the desired value over the first 25 iterations (similarly to Schobert et al. 2019).
The term f ∆τ = 0.2/(ρ d κ d ∆r) regulates the maximum increase of the optical depth to stellar irradiation per one grid cell with the radial length ∆r and allows 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 2 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.
f ∆τ ≤ f d2g,max to prevent f ∆τ from becoming too large in regions with low ρ d .Since τ , f d2g and opacities are mutually dependent through Equations ( 5) and ( 6), we perform 10 sub-iterations within each iteration to evaluate them.

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 ϵ and the energy density of thermal radiation field E R (Dobbs-Dixon et al. 2010): where κ P = κ gas +f d2g κ d (T s ), σ is the Stefan-Boltzmann constant, c is the speed of light, Q irr is the irradiation heating rate, Q visc is the viscous heating rate, and ⃗ 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 d2g evolve only during the first 100 iterations; afterwards 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 where L ⋆ is the stellar luminosity, dτ is the increment of the optical depth across a grid cell of interest, S cell is the irradiated cross-section of the cell, and V cell is its volume.
The viscous heating term is (e.g.D'Angelo & Bodenheimer 2013) where T ij are the components of the viscous stress tensor.During our first experiments with Q visc , we found that a straightforward implementation of this term leads to fluctuating (non-converging) solutions due to the coupling with Equation ( 4).The coupling often results in a feedback loop at spurious locations accross the inner disk rim-if Q visc manages to locally increase T so that ν(α, T ) starts to increase, the local surface density starts to drop through Σ = Ṁ /(3πν), 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 α = α DZ for the purpose of calculating Q 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 visc is calculated correctly in the optically thick regions within the disk interior and the incorrect solution (with too low α) applies mostly inwards from the disk rim where we expect Q irr to dominate anyway (see, for instance, figure 9 in Flock et al. 2019).

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) at the inner radial boundary, with a R being the radiation constant, τ bc = 10 −2 , and T bc = T 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 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 where r fc is the radius of full dust condensation at all heights above the midplane (equation 17 in Ueda et al. 2017).To evaluate E thin R , 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 τ bc = min(τ mid , 1), τ mid being the optical depth to stellar irradiation in the midplane.Finally, we use thick , with T 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 is 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) or we set it to a uniform and density-independent value of T s = 1350 K (or 1550 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.

Synthetic images
We use the Monte Carlo radiative transfer code Radmc-3D (Dullemond et al. 2012) to post-process the results of our hydrostatic modelling.We use ρ, ρ 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 since 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 ρ represents the gas, for which the absorption opacity at each wavelength is considered uniform, κ λ = κ gas , and the scattering opacity is neglected.The second species with the density ρ d = f d2g ρ represents the dust and its κ λ is shown in Figure 1.Radmc-3D offers a 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)3 .
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 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 µm, 100 samples in 7-25 µm, and 30 samples in 25-104 µ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 λ ∈ (1. 5, 1.75, 2, 2.2, 2.45, 2.8, 3.5, 4.2, 5.5, 8, 10.5, 13) µm, covering NIR and MIR bands of the VLTI instruments, assuming the disk inclination i = 52 • and position angle PA = 143 • (Varga et al. 2021).The resolution is 3 × 10 3 pixels along the image edge, i.e. 0.1 mas per pixel.Second-order ray tracing of Radmc-3D is utilized.

Observables
Using the synthetic images, we calculate the halflight radii R hl , fractional disk fluxes F = F disk /F tot = F disk /(F disk +F star ), and interferometric visibilities V at various λ 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 ν with the pixel surface area and considering the distance of HD 163296 being d = 101.5 pc (Wichittanakom et al. 2020).Then, F disk is an integral over all pixels occupied by our disk model and F tot is an integral over the entire image.
The half-light radius is defined via (Leinert et al. 2004;Varga et al. 2021) where r ′ is the radius from the centre of the image plane and the stellar flux is excluded from the calculation.
The synthetic interferometric visibilities are calculated using codes radmc3dPy 4 and pmoired (Mérand 2022).When analyzing the visibilities, we deproject the baselines B using (e.g.Tannirkulam et al. 2008) where i is the disk's inclination and χ = PA B − PA major is the difference between the position angle of the given baseline configuration and the major axis of the on-sky disk projection.

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-on5 and only experiences radial variations in intensity, the interferometric visibility amplitude V is the Hankel transform of the disk profile I ν (for the spatial frequency B/λ), combined linearly with the unit visibility of the unresolved central star: The Bessel function J 0 is an oscillating and vanishing function, so for sufficiently large B/λ, the disk is fully resolved: the integral in Equation ( 16) vanishes to 0 and the visibility becomes constant as function of the baseline.On the other hand, for very small baselines, the  disk is unresolved 6 (V ∼ 1).So as the baseline increases, the visibility decreases from 1 to the saturation value F star /(F disk + F star ).For an intermediate spatial frequency B 1 2 /λ, the visibility reaches a mid-point, which can be measured if a sufficiently wide range of baselines length was explored.In that case: (17) which is very similar to Equation ( 14) defining the halflight radius, if re-written as: with S(x) = 1 for x ≤ 1 and S(x) = 0 for x > 1.
Combining the last 2 equations (17 and 18), we get: (19) The exact relation between B 1 2 and R hl in principle depends on the exact profile I ν .To illustrate whether the 6 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 halo in the denominator of Equation 16.It might result in V ≲ 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.
dependence is strong or weak, we visualize it in Figure 2 for a variety of intensity profiles that are typically used to analyse interferometric data.Denoting I 0 and R 0 the unit intensity and radius, respectively, we consider power-law profiles Figure 2 shows where the intensity profiles reach the half-light radius (left panel), where the corresponding visibility curves reach their mid-point (middle panel), and how B 1 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 For the sake of clarity, let us emphasize that considerations in this sections 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 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 1 2 [m]/λ[µm] ≈ 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).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.Since 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.

Disk structure and temperature profiles
Figure 3 shows the two-dimensional 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 quanti-  ties (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 κ d (T s ).
Figure 3 reveals that the inner disk structure of models M1 and M2 is (unsurprisingly) consistent with 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 in rim , 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 out rim , is more difficult to define.Flock et al. (2016) established R out rim 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 out rim 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 in rim and R out rim 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.36au for M1, R rim = 0.3-0.46au 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 towards 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).2019;Varga et al. 2021).The wavelength range at the bottom panel is limited to λ < 5 µm in order to highlight the differences at shorter wavelengths (the dependence at larger wavelengths is asymptotic and therefore less interesting).
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 ≃ 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 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.

Half-light radii and infrared fluxes
Figure 5 compares R hl and 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 a sense that they are based on parametric brightness distributions fitted to the interferometric measurements and are therefore model-dependent (cf.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 of R hl and F follows the increase of 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 λ 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 towards 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 F at λ < 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 increases, but only weakly near λ ∼ 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 λ ∼ 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 λ ∼ 1.5 µm is nearly the same).

Visibilities
In Figure 6, we plot the visibility amplitude |V | as a function of the deprojected spatial frequency B eff /λ 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 mid-point of the visibility curves (between |V | = 1 and the first bend; see Sec-  tion 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 star /F tot decreases as the fractional disk flux increases7 , 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 of absence of bounces in the displayed K-and L-band observations, they are clearly absent in the H-band observation, 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 super-refractory 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 enclosed8 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.

Infrared emission in detail
Let us now examine the synthetic images9 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 corresponds 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 λ = 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 ≃ 0.3-0.35au (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 λ = 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 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-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.

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  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).
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 determines the optically thin temperature of isolated dust grains as well as the radius where dust grains condense in the midplane (Monnier & Millan-Gabet 2002;Ueda et al. 2017): where T 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 considered10 , with a size distribution f (a) ∝ a −3.5 ranging from 0.05 to 1 µm.Then, κ d (T s ) = 873 and κ d (T ⋆ ) = 4187 cm 2 g −1 , yielding ε ≃ 1/5 and T halo ≃ 1700 K (the latter was found directly from our simulation).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 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 of the disk flux and so they decrease more steeply.The 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 (fewest of all models).The bounce at long wavelengths is the smallest when compared to models M1-M3, yet it is still present.

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;Varga et al. 2021;GRAVITY Collaboration 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 modelling (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 F = F disk /F tot , and interferometric visibilities V with the VLTI multi-band 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 F is the least robust.Interestingly, we found theoretical arguments for R hl being 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 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 frequencydependent 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 since 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 turn 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 non-equilibrium heating of the gas atmosphere.
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 multi-band 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 tweaking 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 F very well).(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 labelled 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 grey curves remain fixed, while the red curve changes only slightly).
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 ≳ 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 ≳ 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 multi-frequency 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 groundwork for future studies.

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 α DZ to evaluate the viscous heating term).Nevertheless, it plays an important role for 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 outwards 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 the vertical cooling starts to act more efficiently than heating by radial radiative diffusion.Since 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 ν (see Section 2.1.1)decreases as well.Because we impose constant Ṁ through the disk, Equation (3) dictates that the local surface density Σ has to increase to compensate for lower ν, 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 Σ ∼ 3 × 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 circumventing such convergence issues.
Our models discussed in the main body of the paper are well converged, reaching the relative change in the temperature during last iterations of the order of 10 −5 .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, is converged already after ∼100 iterations since 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/Ω K (r c ) ≃ 10 6 s and the viscous timescale is t vis = 1/(α MRI h 2 Ω K (r c )) ≃ 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 towards the infrared photosphere.The thermal diffusivity due to radiation in the optically thick limit is (e.g.Lin & Youdin 2015;Jiménez & Masset 2017) leading to χ t ≃ 1.7 × 10 16 cm 2 s −1 for T ≃ 860 K and ρ ≃ 10 −9 g cm −3 .Taking the local height of the infrared photosphere H IR (r c ) ≃ 2.4H ≃ 0.026 au, we obtain t rad = H 2 IR /χ t ≃ 10 7 s.Our estimates yield the inequality t dyn < t rad ≪ t vis , which is consistent with Flock et al. (2016).

Figure 2 .
Figure 2. Comparison of B 1 2 × 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/B 1 2 (center).Finally, we compare (right) B 1 2× 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ν .

Figure 3 .
Figure3.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.Since 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.

Figure 5 .
Figure5.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 λ < 5 µm in order to highlight the differences at shorter wavelengths (the dependence at larger wavelengths is asymptotic and therefore less interesting).

Figure 6 .
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), GRAV-ITY (triangles), and MATISSE (crosses) instruments.The wavelength intervals of synthetic visibilities and observations are the same and are given in the plot legend.

Figure 7 .Figure 8 .
Figure 7. Synthetic images (emission intensity maps) derived from models M1 (top), M2 (middle), and M3 (bottom) at wavelengths λ = 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.

Figure 9 .
Figure9.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 labelled 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 grey curves remain fixed, while the red curve changes only slightly).

Figure 10 .
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.

Table 1 .
Fiducial parameters for the radiative hydrostatic model.

Table 2 .
Overview of individual models.