# Molecules with ALMA at Planet-forming Scales (MAPS). XVII. Determining the 2D Thermal Structure of the HD 163296 Disk

Jenny K. Calahan, Edwin A. Bergin, Ke Zhang, Kamber R. Schwarz, Karin I.

Öberg, Viviana V. Guzmán, Catherine Walsh, Yuri Aikawa, Felipe Alarcón,

Sean M. Andrews, et al.

► To cite this version:

Jenny K. Calahan, Edwin A. Bergin, Ke Zhang, Kamber R. Schwarz, Karin I. Öberg, et al.. Molecules with ALMA at Planet-forming Scales (MAPS). XVII. Determining the 2D Thermal Structure of the HD 163296 Disk. The Astrophysical Journal Supplement Series, 2021, 257, 10.3847/1538-4365/ac143f. insu-03672364

HAL Id: insu-03672364

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

Submitted on 19 Sep 2025

**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>

# Molecules with ALMA at Planet-forming Scales (MAPS). XVII. Determining the 2D Thermal Structure of the HD 163296 Disk

Jenny K. Calahan<fnref n="1" />, Edwin A. Bergin<fnref n="1" />, Ke Zhang<fnref n="1" />,<fnref n="2" />,<fnref n="23" />, Kamber R. Schwarz<fnref n="3" />,<fnref n="24" />, Karin I. Öberg<fnref n="4" />, Viviana V. Guzmán<fnref n="5" />, Catherine Walsh<fnref n="6" />, Yuri Aikawa<fnref n="7" />, Felipe Alarcón<fnref n="8" />, Sean M. Andrews<fnref n="4" />, Jaehan Bae<fnref n="9" />,<fnref n="10" />,<fnref n="24" />, Jennifer B. Bergner<fnref n="11" />,<fnref n="24" />, Alice S. Booth<fnref n="6" />,<fnref n="12" />, Arthur D. Bosman<fnref n="1" />, Gianni Cataldi<fnref n="7" />,<fnref n="13" />, Ian Czekala<fnref n="14" />,<fnref n="15" />,<fnref n="16" />,<fnref n="17" />,<fnref n="18" />,<fnref n="24" />, Jane Huang<fnref n="1" />,<fnref n="4" />,<fnref n="24" />, John D. Ilee<fnref n="6" />, Charles J. Law<fnref n="4" />, Romane Le Gal<fnref n="4" />,<fnref n="19" />,<fnref n="20" />,<fnref n="21" />, Feng Long<fnref n="4" />, Ryan A. Loomis<fnref n="22" />, François Ménard<fnref n="19" />, Hideko Nomura<fnref n="13" />, Chunhua Qi<fnref n="4" />, Richard Teague<fnref n="4" />, Merel L. R. van't Hoff<fnref n="1" />, David J. Wilner<fnref n="4" />, and Yoshihide Yamato<fnref n="7" />

Received 2021 February 17; revised 2021 July 2; accepted 2021 July 12; published 2021 November 3

## Abstract

Understanding the temperature structure of protoplanetary disks is key to interpreting observations, predicting the physical and chemical evolution of the disk, and modeling planet formation processes. In this study, we constrain the two-dimensional thermal structure of the disk around the Herbig Ae star HD 163296. Using the thermochemical code RAC2D, we derive a thermal structure that reproduces spatially resolved Atacama Large Millimeter/submillimeter Array observations (~0.′′12 (13 au)–0.′′25 (26 au)) of $^{12}$CO $J=2-1$, $^{13}$CO $J=1-0, 2-1$, $^{C^{18}O}J=1-0, 2-1$, and $^{C^{17}O}J=1-0$, the HD $J=1-0$ flux upper limit, the spectral energy distribution (SED), and continuum morphology. The final model incorporates both a radial depletion of CO motivated by a timescale shorter than typical CO gas-phase chemistry (0.01 Myr) and an enhanced temperature near the surface layer of the inner disk ($z/r \ge 0.21$). This model agrees with the majority of the empirically derived temperatures and observed emitting surfaces derived from the $J=2-1$ CO observations. We find an upper limit for the disk mass of 0.35 $M_\odot$, using the upper limit of the HD $J=1-0$ and $J=2-1$ flux. With our final thermal structure, we explore the impact that gaps have on the temperature structure constrained by observations of the resolved gaps. Adding a large gap in the gas and small dust additionally increases gas temperature in the gap by only 5%–10%. This paper is part of the MAPS special issue of the Astrophysical Journal Supplement.

Unified Astronomy Thesaurus concepts: Protoplanetary disks (1300); Astrochemistry (75)

Supporting material: figure set

## 1. Introduction

The two-dimensional (radial + vertical) thermal structure of protoplanetary disks has been a long sought-after property due to its importance in interpreting of observations and its importance for disk evolution. Disk temperature sets the physical and chemical evolution of disk material, with subsequent implications for planet formation. The translation of observations into

fundamental disk properties, namely gas mass, strongly depends on the assumed underlying disk temperature at which the observed molecule emits. The mass tracer HD (Bergin et al. 2013; McClure et al. 2016; Bergin & Williams 2017; Trapman et al. 2017; Kama et al. 2020) especially requires a well-defined disk thermal structure as its $J=1-0$ transition has an $E_{up}$ of 128.5 K. This approaches typical temperatures within the warm molecular layer of protoplanetary disks. The temperature throughout the disk regulates physical evolution by setting local sound speeds (Bergin & Williams 2017), turbulent viscosity

<footnote n="1">¹ University of Michigan, 323 West Hall, 1085 S. University Ave., Ann Arbor, MI 48109, USA; jcalahan@umich.edu</footnote>

<footnote n="2">²Department of Astronomy, University of Wisconsin–Madison, 475 N. Charter St., Madison, WI 53706, USA</footnote>

<footnote n="3">³Lunar and Planetary Laboratory, University of Arizona, 1629 E. University Blvd, Tucson, AZ 85721, USA</footnote>

<footnote n="4">⁴Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA</footnote>

<footnote n="5">⁵Instituto de Astrofísica, Pontificia Universidad Católica de Chile, Av. Vicua Mackenna 4860, 7820436 Macul, Santiago, Chile</footnote>

<footnote n="6">⁶School of Physics and Astronomy, University of Leeds, Leeds, LS2 9JT, UK</footnote>

<footnote n="7">⁷Department of Astronomy, Graduate School of Science, The University of Tokyo, Tokyo 113-0033, Japan</footnote>

<footnote n="8">⁸Department of Astronomy, University of Michigan, 323 West Hall, 1085 S. University Ave., Ann Arbor, MI 48109, USA</footnote>

<footnote n="9">⁹Earth and Planets Laboratory, Carnegie Institution for Science, 5241 Broad Branch Road NW, Washington, DC 20015, USA</footnote>

<footnote n="10">¹⁰Department of Astronomy, University of Florida, Gainesville, FL 32611, USA</footnote>

<footnote n="11">¹¹University of Chicago, Department of the Geophysical Sciences, Chicago, IL 60637, USA</footnote>

<footnote n="12">¹²Leiden Observatory, Leiden University, 2300 RA Leiden, The Netherlands</footnote>

<footnote n="13">¹³National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan</footnote>

<footnote n="14">¹⁴Department of Astronomy and Astrophysics, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA</footnote>

<footnote n="15">¹⁵Center for Exoplanets and Habitable Worlds, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA</footnote>

<footnote n="16">¹⁶Center for Astrostatistics, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA</footnote>

<footnote n="17">¹⁷Institute for Computational & Data Sciences, The Pennsylvania State University, University Park, PA 16802, USA</footnote>

<footnote n="18">¹⁸Department of Astronomy, 501 Campbell Hall, University of California, Berkeley, CA 94720-3411, USA</footnote>

<footnote n="19">¹⁹Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France</footnote>

<footnote n="20">²⁰IRAP, Université de Toulouse, CNRS, CNES, UT3, 31400 Toulouse, France</footnote>

<footnote n="21">²¹IRAM, 300 rue de la piscine, F-38406 Saint-Martin d'Hères, France</footnote>

<footnote n="22">²²National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903, USA</footnote>

<footnote n="23">²³NASA Hubble Fellow.</footnote>

<footnote n="24">²⁴NASA Hubble Fellowship Program Sagan Fellow.</footnote>
(Shakura & Sunyaev 1973; Penna et al. 2013), and vertical
flaring of the disk (Kenyon & Hartmann 1987), which can in
turn set the angle of incidence of stellar irradiation, further
changing the thermal structure. From a chemical perspective, gas
temperature controls the rate of gas-phase exothermic reactions.
Temperature also dictates the rates of gas-phase deposition and
thermal sublimation, effectively prescribing the relative radial
chemical composition of ices in the midplane. These sublimation
fronts, or snowlines, have been theorized to be favorable sites for
planet formation (Hayashi 1981; Stevenson & Lunine 1988;
Zhang et al. 2015; Schoonenberg & Ormel 2017), and may
affect the chemical composition of a planetary embryo or
accreting atmosphere (e.g., C/O; Öberg et al. 2011; Öberg &
Bergin 2016; Cridland et al. 2020).

The most robust attempts to uncover the temperature profile of protoplanetary disks involve a convergence of observations of multiple gas emission lines and thermochemical modeling (e.g Kama et al. 2016; Woitke et al. 2019; Rab et al. 2020; Calahan et al. 2021). With the advent of the Atacama Large Millimeter/submillimeter Array (ALMA), spatially resolved observations of disks became feasible, providing the first empirical measurements of the radial distribution of dust and gas at high angular resolution. This gave much-needed constraints for protoplanetary disk models. The commonly observed ¹²CO J = 2 – 1 is optically thick in most gas-rich disks due to it being the second most abundant molecule in the gas phase after H₂. Lower-energy transitions and less-abundant isotopologues emit from lower vertical layers, resulting in intensity profiles that are affected by temperature and CO surface density/abundance where that species is emitting. Thus, joint modeling of multiple spatially resolved CO isotopologues provides radial and vertical constraints on the temperature and CO chemistry.

This study determines the 2D temperature structure of the disk
around HD 163296, as part of the ALMA-MAPS project, which
observed five protoplanetary disks at high resolution. Each disk is
detailed in Öberg et al. (2021). HD 163296 is a Herbig Ae star
with a stellar mass of 2.0 M<sub>☉</sub> (Andrews et al. 2018; Wichitana-
kom et al. 2020) surrounded by a massive disk (8 × 10<sup>-3</sup>
−5.8 × 10<sup>-1</sup> M<sub>☉</sub>; Kama et al. 2020, and references therein) 101 pc
away (Bailer-Jones et al. 2018; Gaia Collaboration et al. 2018).
Both millimeter continuum and scattered-light observations show
a significant substructure (Grady et al. 2000; Wisniewski et al.
2008; Benisty et al. 2010; Garufi et al. 2014; Monnier et al. 2017;
Muro-Arena et al. 2018; Isella et al. 2018; Dent et al. 2019; Rich
et al. 2020) including three gaps in millimeter continuum
observations at 10, 48, and 86 au. These gaps have been prime
targets to test planet formation theories; planets ranging in mass
from 0.07–4.45 M<sub>J</sub> have been theorized to exist within the gaps in
HD 163296 (Zhang et al. 2018). The temperature of both the gas
and dust populations in the gaps can be used to directly inform
planet formation models. Kinematic studies of HD 163296
suggest the existence of a ~2 M<sub>J</sub> planet at 260 au (Pinte et al.
2018, 2020). Further analysis of the velocity structure by Teague
et al. (2019) mapped out the 3D gradient of velocity using ¹²CO
J = 2 − 1 from the DSHARP project (Andrews et al. 2018), and
in Teague et al. (2021). Meridional flows are found at the location
of the continuum gaps, indicative of ongoing planet formation.
The protoplanetary disk around HD 163296 presents an excellent
case for determining the 2D thermal profile using some of the
highest resolution data available for a disk in which planets are
believed to be actively forming.

Two-dimensional temperature structures for HD 163296 have been modeled previously using a combination of spatially unresolved observations of the CO rotational ladder, and other higher-energy molecular and atomic transitions (Qi et al. 2011; de Gregorio-Monsalvo et al. 2013; Rosenfeld et al. 2013) with a few studies utilizing thermochemical modeling to match numerous observations (Tilling et al. 2012; Fedele et al. 2016; Woitke et al. 2019). These studies found that HD 163296 is a relatively cooler disk compared to other disks surrounding Herbig stars. Resolved observations of ¹²CO J = 2 − 1 have also revealed insight into the temperature within the disk gaps (Rab et al. 2020).

The observations from the ALMA-MAPS program consist of highly resolved (~0.′′12–0.′′25), and high signal-to-noise ratio observations of ¹²CO J = 2 − 1, ¹³CO J = 2 − 1, 1 − 0, C¹⁸O J = 2 − 1, 1 − 0, and C¹⁷O J = 1 − 0. These observations show structure in the radial intensity profiles (Law et al. 2021a) and CO column density (Zhang et al. 2021; hereafter Z21) and allow for a new derivation of the 2D thermal structure.

This paper is organized as follows: We detail our modeling
procedure and describe the observations in Section 2. In Section 3,
we present our best-fit thermal model and the necessary alterations
needed to fit each observed line. In Section 4, comparisons of the
final model to the empirically derived temperature structure and
emitting surfaces are discussed along with our derivation of an
upper mass limit for the HD 163296 disk and we explore CO/H₂
degeneracy and temperature effects within the two largest gaps of
the HD 163296 disk. We summarize our findings in Section 5.

## 2. Methods and Observations

2.1. *Observations*

We use observations of ¹²CO J = 2 − 1, C¹⁸O J = 2 − 1, 1 − 0, ¹³CO J = 2 − 1, 1 − 0, and C¹⁷O J = 1 − 0 toward the HD 163296 disk. For this study, we use images which have an average spatial resolution of 26 au for the 1 − 0 transitions (0.′′25) and 13 au for the 2 − 1 transitions (Öberg et al. 2021, and Table 1). The final image cubes used for this study had a spectral resolution of 200 m s⁻¹, and used a robust = 0.5 weighting. The robust = 0.5 images were used as they provided the highest absolute resolution. These images were “JvM corrected” (Jorsater & van Moorsel 1995), which refers to the method of scaling the residuals in the image cube to have identical units to the CLEAN model. This ensures that the final image used for moment zero maps, and subsequent radial profiles, have the correct units (Jy CLEAN beam⁻¹); see Czekala et al. (2021) for a detailed explanation. The MAPS program summary, along with data reduction and calibration specifics can be found in Öberg et al. (2021), and the methods of the imaging process are thoroughly described in Czekala et al. (2021). Radial intensity profiles are created for each CO line for comparison with our model results. Radial profile derivations are presented in Law et al. (2021a). The package BETTERMOMENTS<fnref n="25" /> (Teague & Foreman-Mackey 2018) is used to extract the moment zero map which is then used to calculate the radial intensity profile using GOFISH<fnref n="26" /> (Teague 2019). We follow the same methods in creating our simulated radial intensity profiles.

In addition to ALMA observations of CO, we also compared
our model to the Herschel Space Observatory observation of HD

<footnote n="25">²⁵ https://github.com/richteague/bettermoments</footnote>

<footnote n="26">²⁶ https://github.com/richteague/gofish</footnote>
<table><caption>**Table 1**
ALMA Observations Summary</caption><thead><tr><th>Transition</th><th>Frequency (GHz)</th><th>E<sub>up</sub> (K)</th><th>Beam (" × ")</th><th>PA (deg)</th><th>rms (mJy beam<sup>-1</sup>)</th></tr></thead><tbody><tr><td><sup>12</sup>CO J = 2 − 1</td><td>230.538</td><td>16.6</td><td>0.14 × 0.10</td><td>−76.3</td><td>0.561</td></tr><tr><td><sup>13</sup>CO J = 1 − 0</td><td>110.201</td><td>5.3</td><td>0.28 × 0.22</td><td>−89.1</td><td>0.434</td></tr><tr><td><sup>13</sup>CO J = 2 − 1</td><td>220.399</td><td>15.8</td><td>0.14 × 0.11</td><td>−76.6</td><td>0.541</td></tr><tr><td>C<sup>18</sup>O J = 1 − 0</td><td>109.782</td><td>5.3</td><td>0.28 × 0.22</td><td>−88.8</td><td>0.449</td></tr><tr><td>C<sup>18</sup>O J = 2 − 1</td><td>219.560</td><td>15.8</td><td>0.14 × 0.11</td><td>−76.5</td><td>0.385</td></tr><tr><td>C<sup>17</sup>O J = 1 − 0</td><td>112.359</td><td>5.4</td><td>0.28 × 0.21</td><td>−89.5</td><td>0.528</td></tr></tbody></table>

Note. These data were taken from values in Öberg et al. (2021) and Czekala et al. (2021), where details regarding the data reduction can also be found.

$J=1-0$ and $2-1$ toward HD 163296. The HD observations were made using the Photoconductor Array Camera and Spectrometer (PACS) instrument and were spatially and spectral unresolved. Kama et al. (2020) used archival data from the DIGIT program (PI N. J. Evans) to determine the upper flux limits for 15 Herbig Ae/Be disks, including HD 163296. They found an upper limit for the flux of the $J=1-0$ line of $\le 0.9 \times 10^{-17} \text{ W m}^{-2}$ and for the $J=2-1$ they determined a flux of $\le 2.4 \times 10^{-17} \text{ W m}^{-2}$.

## 2.2. RAC2D Physical Structure

In this study, we use the time-dependent thermochemical code RAC2D<sup>27</sup> (Du & Bergin 2014) to model the thermal structure of the HD 163296 disk, and create simulated observations to compare with the ALMA data. A brief description of the physical code of RAC2D is given below; a detailed description of the code can be found in the aforementioned paper.

RAC2D takes into account both the disk gas and dust structure while simultaneously computing the temperature and chemical structure over time. Our model consists of three mass components: gas, small-dust ($\le$ micrometer), and large-dust grains ($\le$ millimeter). The distribution of each component is described by a global surface density distribution (Lynden-Bell & Pringle 1974), which is widely used in the modeling of protoplanetary disks and corresponds to the self-similar solution of a viscously evolved disk.

$$ \Sigma(r) = \Sigma_c \left(\frac{r}{r_c}\right)^{-\gamma} \exp\left[-\left(\frac{r}{r_c}\right)^{2-\gamma}\right], \quad (1) $$

$r_c$ is the characteristic radius at which the surface density is $\Sigma_c/e$ and $\gamma$ is the power-law index that describes the radial behavior of the surface density.

A density profile for the gas and dust components can be derived from the surface density profile and a scale height:

$$ \rho(r, z) = \frac{\Sigma(r)}{\sqrt{2\pi}h(r)} \exp\left[-\frac{1}{2}\left(\frac{z}{h(r)}\right)^2\right], \quad (2) $$

$$ h = h_c \left(\frac{r}{r_c}\right)^{\Psi}, \quad (3) $$

where $h_c$ is the scale height at the characteristic radius, and $\Psi$ is a power index that characterizes the flaring of the disk structure.

Both dust populations follow a Mathis–Rumpl–Nordsieck grain distribution $n(a) \propto a^{-3.5}$ (Mathis et al. 1977). The small-

dust grains have radii between $5 \times 10^{-3}-1$ µm, and the large grains have radii between $5 \times 10^{-3}-10^3$ µm. The following description of the dust-grain population used in this study is adopted from Z21, which used both the HD 163296 spectral energy distribution (SED) and millimeter continuum observations to constrain the dust population. The gas and small-dust grains are spatially coupled and extend to 600 au. The large-dust population is settled in the midplane with a smaller vertical extent ($h=1.69$ au) and radial extent (240 au). This settled large-grain population is the result of dust evolution, namely vertical settling to the midplane and radial drift. The large-grain population has a unique, nonsmooth, surface density profile that reproduces the millimeter continuum observations of the HD 163296 disk (Z21). The dust composition adopted is consistent across the MAPS Collaboration, and opacity values are calculated based on Birnstiel et al. (2018). Large dust grains consist of water ice (Warren & Brandt 2008), silicates (Draine 2003), troilites, and refractory organics (Henning & Stognienko 1996). Small-dust grains consist of 50% silicates and 50% refractory organics.

## 2.3. RAC2D Dust and Gas Temperature

The code computes a dust and gas temperature after initializing RAC2D with a model density structure for each population, a stellar spectrum, and external UV radiation field ($G_0=1$). The stellar spectrum we use is a composite of multiple UV and X-ray observations, and is further detailed in Z21. The determination of dust and gas temperature is an iterative process that is allowed to change over time due to the evolving chemical composition. The dust thermal structure is calculated first using a Monte Carlo method to simulate photo absorption and scattering events. This also results in a radiation field spanning from centimeter to X-ray wavelengths.

We adopt the reaction rates from the UMIST 2006 database (Woodall et al. 2007) for the gas-phase chemistry with additional rates considering the self-shielding of CO, H$_2$, H$_2$O, and OH, dust-grain surface chemistry driven by temperature, UV, cosmic rays, and two-body chemical reactions on dust-grain surfaces (see references given by Du & Bergin 2014). Chemical processes also provide heating or cooling to the surrounding gas. These mechanisms, along with stellar and interstellar radiation, drive the thermal gas structure. Our study explores models that account for 0.01 Myr of chemical evolution. This is a relatively short period of time compared to disk lifetimes and is motivated by an average CO-processing timescale, which is found to be on the order of 1 Myr (Bergin et al. 2014; Bosman et al. 2018; Schwarz et al. 2018). The calculated depletion profile of CO is motivated by observations, and thus accounts for any long-timescale chemical effects that have taken place in the disk’s history; we

<footnote n="27"><sup>27</sup> https://github.com/fjdu/rac-2d</footnote>
did not want to further process CO by running the disk model for longer than 0.01 Myr. It is also worth noting that the exact timescale will not affect our final temperature structure nor CO column density distribution; it will only affect the CO depletion profile. CO will be created or destroyed in our chemical network over time, and the CO depletion profile acts to counter any over- or underabundance of CO being produced. While localized variations in the CO abundance may affect the gas temperature, the global temperature structure is not strongly affected by the CO abundance. Finally, at the end of a given run, we extracted the dust and gas thermal profiles.

Simulated line images for CO, ¹³CO, C¹⁸O, C¹⁷O, and HD are necessary for our comparison to observations. We did not model isotopologue fractionation in this chemical network, as fractionation of CO is not significant in a massive disk like HD 163296 (Miotello et al. 2014). Thus, we computed ¹³CO and C¹⁸O abundances based on interstellar medium ratios of ¹²CO/¹³CO = 69, ¹²CO/C¹⁸O = 557, C¹⁸O/C¹⁷O = 3.6 (Wilson 1999). Given these abundances and the local gas temperature, RAC2D computes synthetic channel maps using a ray-tracing technique. We then convolved these simulated observations with the corresponding ALMA CLEAN beam to make direct comparisons to data. To directly compare to observations, we created simulated integrated radial intensity profiles of each of the CO lines, following the methods from Law et al. (2021a). Our simulated observations have the same beam sizes, spectral resolution, and pixel size as what is reported in Öberg et al. (2021). To recreate the unresolved integrated flux measurements for HD, we convolved our model HD lines with a Gaussian profile corresponding to the velocity resolution of the Herschel PACS instrument: ~300 km s⁻¹ at ~112 µm and ~100 km s⁻¹ at ~51 µm (Poglitsch et al. 2010).

We began with a model of the HD 163296 disk from Z21, which starts with an initial set of disk parameters taken from the literature (see references within Zhang et al. 2021) and applied the observed UV, X-ray, and photosphere stellar spectra. These authors then determined a gas and dust density and dust temperature structure by matching the SED and ALMA continuum image using RADMC3D (Dullemond et al. 2012). Given this initial density and dust temperature distribution, RAC2D is then used to compute the gas temperature and disk chemistry. Z21 found that in order for the simulated radial profiles to match what is observed, they must modify the CO abundance relative to H₂ as a function of radius. They use the difference between C¹⁸O *J* = 2 – 1 as observed and as simulated to predict a CO depletion profile. That depletion profile is shown in Figure 1. The initial chemical abundances used in Z21 are shown in Table 2 and final model parameters are summarized in Table 3. In the study described here, we use the Z21 model as a starting point to derive a 2D temperature map incorporating additional constraints from CO isotopologues, multiple higher-level CO transitions, and HD flux. We note that while any revisions made were important to constrain the disk temperature structure, they do not affect the results reported in Z21, and the final CO column densities from this model agree well with what is found in Z21.

## 3. Results

### 3.1. CO Depletion

In Z21, depletion of CO was calculated based on a model of the disk based on previous determinations of disk parameters

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Multiplicative depletion factor for the initial CO abundance (CO/H = 1.4 × 10⁻⁴) as shown in Table 2. Light gray lines indicate the location of the dust rings (solid) and gaps (dashed). The line in blue shows the CO depletion as derived by Z21 and is determined after an initial thermochemical model of the HD 163296 disk and is based on the C¹⁸O *J* = 2 – 1 observation, and acts as a starting point for our CO depletion determination. The green line represents the final CO depletion, also motivated by the C¹⁸O *J* = 2 – 1 flux, but it is taken into account at the beginning of the chemical and thermal calculations.</figcaption></figure>

<table><caption>**Table 2.**
Initial Chemical Abundances for Final Model</caption><thead><tr><th></th><th>Abundance Relative to Total<br>Hydrogen Nuclei</th></tr></thead><tbody><tr><td>H<sub>2</sub></td><td>5 × 10<sup>-1</sup></td></tr><tr><td>HD</td><td>2 × 10<sup>-5</sup></td></tr><tr><td>He</td><td>0.09</td></tr><tr><td>CO<sup>a</sup></td><td>See Figure 1</td></tr><tr><td>N</td><td>7.5 × 10<sup>-5</sup></td></tr><tr><td>H<sub>2</sub>O (ice)</td><td>1.8 × 10<sup>-4</sup></td></tr><tr><td>S</td><td>8 × 10<sup>-8</sup></td></tr><tr><td>Si<sup>+</sup></td><td>8 × 10<sup>-9</sup></td></tr><tr><td>Na<sup>+</sup></td><td>2 × 10<sup>-8</sup></td></tr><tr><td>Mg<sup>+</sup></td><td>7 × 10<sup>-9</sup></td></tr><tr><td>Fe<sup>+</sup></td><td>3 × 10<sup>-9</sup></td></tr><tr><td>P</td><td>3 × 10<sup>-9</sup></td></tr><tr><td>F</td><td>2 × 10<sup>-8</sup></td></tr><tr><td>Cl</td><td>4 × 10<sup>-9</sup></td></tr></tbody></table>

**Note.**

and the C¹⁸O *J* = 2 – 1 radial intensity profile. The CO was depleted throughout the disk at the start of the ray-tracing method that creates the simulated observations, i.e., after the temperature and chemical evolution of the disk. This depletion means the CO abundance is reduced from its expected value (~10⁻⁴ relative to H₂) in layers where the dust temperature is above the CO sublimation point and the gas is self-shielded from radiation. Since CO is a significant coolant in the disk-surface layers where the dust and gas are thermally decoupled, any depletion of CO may affect the thermal structure. Thus, for this study, we recalculated the CO depletion factors. We started by applying the Z21 CO depletion to the initial CO abundance. After running the model for 0.01 Myr, there were small inconsistencies between the simulated and observed C¹⁸O *J* = 2 – 1 radial emission profiles (see Figure A1). We then calculated a new CO depletion profile by determining what

<footnote>a CO starts with an abundance of 1.4 × 10⁻⁴ with an imposed radial depletion profile as shown in Figure 1.</footnote>
<table><caption>**Table 3**
Gas and Dust Population Parameters: Best-fit Model Values</caption><thead><tr><th></th><th>Gas</th><th>Small Dust</th><th>Large Dust</th></tr></thead><tbody><tr><td>Mass (M<sub>c</sub>)</td><td>0.14</td><td>2.6 × 10<sup>-4</sup></td><td>2.4 × 10<sup>-3</sup></td></tr><tr><td>Ψ</td><td>1.08</td><td>1.08</td><td>1.08</td></tr><tr><td>γ</td><td>0.8</td><td>0.8</td><td>0.1</td></tr><tr><td>h<sub>c</sub> (au)</td><td>14.5</td><td>14.5</td><td>...</td></tr><tr><td>r<sub>c</sub> (au)</td><td>165</td><td>165</td><td>...</td></tr><tr><td>r<sub>in</sub> (au)</td><td>0.45</td><td>0.45</td><td>0.45</td></tr><tr><td>r<sub>out</sub> (au)</td><td>600</td><td>600</td><td>240</td></tr></tbody></table>

Note. Final values of the HD 163296 model that reproduces the CO, HD, and SED observations. Small dust-grain sizes range from $5 \times 10^{-3}-1$ µm, large dust-grain sizes range from $5 \times 10^{-3}-1 \times 10^{3}$ µm. The large-dust population does not have associated $h_c$ nor $r_c$ because the surface density is set by continuum observations, thus it is nonsmooth and not dictated by Equations (1)–(3).

factor of increase or decrease in CO flux would be needed to reproduce the observation at each radii (using the same radial resolution as Z21). That factor was then applied to the original CO depletion profile at the corresponding radius. The model was then run again with the new CO depletion profile. We needed to iterate this process three times, and stopped iterating once the $\chi^2$ value (using stat.chisquare function from the scipy package and comparing the simulated and observed C<sup>18</sup>O $J=2-1$ radial profile intensity values) had dropped well below that of the first attempts. The first three attempts had a $\chi^2$ value of 8.06–9.72, and the final model $\chi^2 = 2.47$. The final depletion profile is shown in Figure 1, and on average the newly calculated CO depletion is 33% more depleted than that derived in Z21. Most of the difference is at large radii where the profiles vary by only 16% within 200 au.

There is a strong decrease in CO abundance beyond ~250 au, which is not seen in the CO depletion profile calculated by Z21. The CO depletion profile from Z21 accounts for an initial CO depletion in HD 163296 plus 1 Myr of CO processing accounted for in the RAC2D chemical evolution. In our model, we include only 0.01 Myr of additional chemical processing, thus we do not attempt to model the full evolution of CO. In the Z21 model, their derived CO depletion profile takes into account chemical processing that occurs over the course of 1 Myr. The difference in CO abundance for our HD 163296 disk model as compared to the Z21 model is then attributed to chemical CO depletion mechanisms that occur over 1 Myr. In the Z21 model, fully self-shielded CO becomes frozen out and interacts with OH that is frozen out onto grains and creates CO<sub>2</sub>. In our model, we treat past CO chemistry as an initial condition. We find it necessary to additionally deplete CO beyond ~240 au in the model that runs for only 0.01 Myr, since the pathway to convert CO and CO<sub>2</sub> is unavailable within that time.

## 3.2. Additional Heating

While our initial model with the updated CO depletion profile reproduced most of the observed CO lines, <sup>12</sup>CO $J=2-1$ was underpredicted on average by a factor of 1.6 within 100 au, and a factor of 2.2 within 35 au. We completed a thorough exploration of the parameter space including gas mass, small-dust mass, surface density power-law index ($\gamma$), flaring index ($\Psi$), scale height ($h_c$), critical radius ($r_c$), and outer radius (see Appendix A). However, we found no combination of parameters that could

improve the <sup>12</sup>CO $J=2-1$ flux agreement while leaving the other lines consistent with observations. This result, together with the high optical depth of the <sup>12</sup>CO $J=2-1$ transition, suggests that the discrepancy is due to a higher temperature in the CO-emitting layer. This has also been seen in de Gregorio-Monsalvo et al. (2013); <sup>12</sup>CO $J=3-2$ was observed to be brighter than their model, which matches the outer disk. To solve this issue they increased the gas temperature beyond the dust temperature within 80 au in their HD 163296 disk model.

This higher gas temperature requires additional heating, beyond the level generated by the UV and X-ray field in our model. We began by isolating the layers in which the gas and dust thermal structure are decoupled. Within those regions, with the goal of representing excess heating due to radiation from the star, we artificially increased the gas temperature after the thermochemical calculation, and simulated new CO observations. We increased the temperature in this region following a power law dependent on radius, and by a constant amount vertically (see the equation below). There was no realistic amount of heating within the decoupled layer that would increase the intensity of the <sup>12</sup>CO $J=2-1$ line to match what is observed; any excess heating only affected emission from the inner few astronomical units. This then suggests that the <sup>12</sup>CO $J=2-1$ primarily emits from the region lower in the disk, in which the dust and gas temperatures are coupled, and the necessary excess heating would decouple the gas and dust temperatures.

We increase the gas temperature radially following the dust thermal profile. Increasing the gas temperature as function of radius following a power law,

$$T_{\text{new}} = T_{\text{original}} \times \left( \frac{r}{100 \text{ au}} \right)^{-0.4}$$

above a constant ratio of height over radius ($z/r$) = 0.21 reproduces the observed <sup>12</sup>CO $J=2-1$ radial profile. Outside of this region, the gas and dust temperatures remain coupled. Various $z/r$ limits were explored, and $z/r$ = 0.21 produced a result that appeared to significantly affect the <sup>12</sup>CO $J=2-1$ line while not altering any other lines (see Figure 2 for a sample of the parameter exploration results). It can then be presumed that the emitting surface of <sup>12</sup>CO $J=2-1$ exists at or above $z/r$ = 0.21 within 100 au, and all other lines emit from below these heights (this is supported by further analysis probing the emitting surfaces of each $J=2-1$ transition, see Section 4.1 and Law et al. (2021b)). Higher transitions of CO, up to <sup>12</sup>CO $J=23-22$ have been observed using the SPIRE instrument on Herschel and are compiled in Fedele et al. (2016). These transitions will primarily emit from the very inner regions of the disk, and would be affected by the increase in temperature. Our final model reproduces all of the observed fluxes of the 19 higher-level transitions of CO, most within $1\sigma$ uncertainty, and five of the transitions fit within $2\sigma$ or $3\sigma$ uncertainty (see Table B1). While the flux predictions from the model all fall within $3\sigma$, all are on the fainter end of the observed flux. Additionally, the HD $J=2-1$ flux does not change significantly, increasing only by 7%, remaining well below the observed upper limit. A model without the additional heat produces a CO flux from the $J=11-10$ transition that is more than twice as small as our final model. On average, for the transitions between $J=11-10$ and $J=23-22$, a model
<figure><img src="image_4.png" /><figcaption>**Figure 2.** These figures represent the effect of excess heating on the ¹²CO J = 2–1 line profile given different height over radius (z/r) limits (0.17, 0.21, 0.25 columns) and R values (100 and 200 au, dashed line and dotted line) using (r/R₀)⁻⁰.⁴ to calculate the excess factor of heat in the regions above a certain z/r. We sought to find a z/r and R₀ value that would reproduce the observed ¹²CO J = 2–1 radial profile (solid red lines). The above is just a sample of the z/r and R₀ values explored, with our final model using a z/r limit of 0.21 and R₀ = 100 au. The bottom left plot shows the observed emitting surface of ¹²CO J = 2–1 (blue) and the z/r values explored. The bottom right plot shows the dust temperature over radius at different heights in the disk, which follow an inverse power-law function, and each dashed black line is a power law following r⁻⁰.⁴, which is what we use to determine the excess heat factor.</figcaption></figure>

without the additional heat produces a CO flux over four times lower than our current model.

The majority of the excess heat is added within 25 au where gas temperatures are over twice the original output from our thermochemical model. A likely source of this excess heat is mechanical heating from processes such as stellar winds or the dissipation of gravitational energy from accretion through the disk onto the star (e.g., Glassgold et al. 2004; Najita & Ádámko-vics 2017). Mechanical heating from such phenomena is expected to be prominent in the inner disk (<10 au), and is not accounted for in RAC2D. Another possible heating source originates from polycyclic aromatic hydrocarbons (PAHs), as photoelectric heating of small grains (Draine & Li 2001; Li & Draine 2001) is one of the main heating mechanisms in this region after direct UV heating from the star. The PAH abundance might be enhanced in the inner disk if dust sintering is taken into account (Okuzumi et al. 2016). As dust grains pass the sublimation front of their constituent material, PAHs can be released, enhancing the effect of photoelectric heating near which ¹²CO 2–1 emits. PAH emission toward HD 163296 has been observed as a part of the Infrared Space Observatory Short Wavelength Spectrometer atlas (Sloan et al. 2003) and modeled by Seok & Li (2017). Their best-fit model uses a PAH abundance of 8 × 10⁻⁷ Mₐ, which is 1.5× the abundance used in our RAC2D model. In our model and assumed physical setup, the excess heat is necessary in a region where the gas and dust are thermally coupled. Any increase in PAH abundance in our model will not change the temperature at which the majority of ¹²CO J = 2–1 emits. It only affects the temperature in thermally decoupled layers where dust densities are low and gas and dust collisions do not occur often. However, it remains a possibility given an alternate setup of HD 163296 in which ¹²CO J = 2–1 emits in the thermally decoupled layers.

### 3.3. Final Model

The model that best reproduces the CO radial profiles observed toward HD 163296 was achieved by altering the CO depletion profile derived by Z21, and increasing the gas temperature above a z/r = 0.21 within the inner 100 au of the disk. The CO radial profiles derived from the final model are shown in Figure 3 together with the observed profiles. The final gas and dust thermal structures are shown in Figure 4. Comparing the gas and dust temperatures, we find that below a z/r ≃ 0.25 the dust and gas are thermally coupled (with the exception of the increased gas temperature added artificially). Right above z/r ≃ 0.25 the dust is hotter than the gas by factors of 25%–50%. This region has been referred to as the “undershoot” region (Kamp & Dullemond 2004) and occurs in disks where the gas density increases sufficiently for line cooling to become very efficient, and occurs around the τ = 1 surface. In our model, atomic oxygen dominates cooling in this region, and cooling via gas-grain collisions is particularly inefficient (see Figure 5). At the very surface of the disk, the gas temperature then increases drastically, overtaking the dust temperature significantly, which tends to plateau above the undershoot region. Directly above the undershoot region, photoelectric heating and the vibrational transitions of H₂ are the most significant heating processes, with direct heating of the gas via X-ray radiation overtaking them beyond z/r ≃ 0.45–0.5. This gas/dust temperature decoupling has been predicted (Kamp & Dullemond 2004; D’Alessio et al. 2005) and seen in models before (Tilling et al. 2012; Woitke et al. 2019; Rab et al. 2020).

It is worth noting that the C¹⁷O and C¹⁸O 1–0 lines are the least well-fit to the observations. We also find that these transitions originate from deep within the disk, only slightly above the midplane, and CO snow surface. In Z21 it was noted that the C¹⁷O column densities appear to be higher than those
Integrated Intensity Profiles of CO

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Integrated radial intensity profiles of $^{12}$CO and its isotopologues $^{13}$CO, $^{18}$O, and $^{17}$O as observed (solid line) and as simulated by the thermochemical code RAC2D (dashed line). The gray shaded regions correspond to the FWHM of the corresponding beam for each observed line. This presents our best-fit model to the observations of the HD 163296 disk and represents the thermal structure shown in Figure 4.</figcaption></figure>

estimated from $^{18}$O (correcting for optical depth). This is also the case for an earlier $^{13}$C$^{18}$O detection, which appears to require a higher CO column than inferred from $^{18}$O $J=2-1$ (Z21). A similar analysis is performed for the MAPS data in Z21. It is suggested that the CO abundance might be higher in the midplane as the icy pebbles drift inwards and sublimate CO inside the snowline. Booth et al. (2019) estimate a disk mass of 0.21 $M_{\odot}$ using $^{13}$C$^{17}$O $J=3-2$, a value 50% more massive than our final model. However, if CO is enhanced near the midplane, a larger mass may not be necessary. Such an effect is not accounted for in our analysis, but, at face value, would be consistent with our results in that a locally higher abundance of CO in the midplane, perhaps inside the CO snowline, would increase the emission of optically thin tracers (such as $^{17}$O) while having a reduced effect on the emission of the more abundant isotopologues. This process has a negligible effect on the thermal structure as the dust and gas are strongly coupled in these layers.

### 3.3.1. CO and CO$_2$ Snowline

To determine the location of the CO snowline in our model, we determine the radial location at which there are equal parts CO frozen onto the dust and in the gas phase. RAC2D includes

adsorption of CO (and other species) onto dust grains as well as desorption from the dust surface either thermally, or via UV photons or cosmic rays (see Du & Bergin 2014). Using this metric, our CO snowline is located at 60 au at a temperature of 18 K. This falls between the largest gaps in the continuum, and is consistent with Z21. Qi et al. (2015) used observations of $N_2H^+$ toward HD 163296 to determine the CO snowline as $N_2H^+$ formation is inhibited in the presence of gas-phase CO. When CO is frozen out $N_2H^+$ can exist, emitting as a ring with the innermost edge corresponding to the CO snowline. Using this method, and a model of the HD 163296 disk, they found a CO snowline at $74_{-5}^{+7}$ au (corrected for the pre-Gaia distance). Our midplane snowline location is in good agreement with the $N_2H^+$ derived snowline location, especially when considering findings from van't Hoff et al. (2017) that show the $N_2H^+$ column density peaks at 5 au or more beyond the midplane CO snowline. Additionally, we can predict the CO$_2$ snowline to be at 4 au, at a temperature of 65 K. These freeze-out temperatures depend on the given desorption energy. RAC2D uses values from Garrod et al. (2008), which assumed an amorphous ice surface. Assuming a different desorption energy will alter the derived snowline locations. It is important to also note that these snowlines are based on a thermochemical model
<figure><img src="image_6.png" /><figcaption>**Figure 4.** The two-dimensional profile of gas temperature, dust temperature, percentage difference between the two, and CO abundance. The dashed black or white lines in each plot show where $z/r = 0.21$, which corresponds with the $^{12}$CO $J = 2 - 1$ surface within 100 au. This represents the model that best represents the observed radial profile, SED, and an HD flux that agrees with the currently derived upper limit. In the percentage difference plot, one can see where the dust and gas temperatures are coupled (gray), the “undershoot” region (Kamp & Dullemond 2004), where dust is warmer than the gas (blue), and where gas is warmer high up in the atmosphere (red).</figcaption></figure>

influenced by observations that do not directly probe the midplane.

# 4. Discussion and Analysis

## 4.1. Emitting Surface

An additional observational constraint on our model is the resolved emission heights of $^{12}$CO, $^{13}$CO, and C$^{18}$O $J = 2 - 1$ as measured in Law et al. (2021b). We calculate the emitting surfaces of $^{12}$CO, $^{13}$CO, and C$^{18}$O using simulated image cubes and the same methods as applied to the observations by Law et al. (2021b). This method depends on the ability to resolve the front and back side of a given disk over multiple channels of an image cube, and assumes azimuthal symmetry and gas rotating in circular orbits (Pinte et al. 2018). The emitting height extraction tools are found in the Python package `DISKPROJEC-TION`.<fnref n="28" /> A comparison of the emitting layers of our model and observations is shown in Figure 6.

The calculated emitting layers for both simulated and observed data from Law et al. (2021b) exist at similar heights, and agree within $1\sigma$ at most radii. No additional changes to the model were necessary to arrive at this agreement between the model and observation. Our modeled emission height for $^{12}$CO $J = 2 - 1$ reproduces what was derived using the ALMA observations up

until ~125 au; beyond this radius and up to 250 au, the model's $^{12}$CO emits on average 11 au lower in the disk. When looking at the $^{13}$CO-emitting height, the model reproduces the observed heights beyond ~125 au. This suggests that our model does not completely represent the CO vertical distribution, as it systematically produces a lower $^{12}$CO $J = 2 - 1$ emission height beyond ~100 au. However, this fit is quite good considering the detailed physics and chemistry included in the model. The C$^{18}$O emitting heights agree well with what is observed, especially within ~100 au, where both our model and the observations pick up on some substructure, a bump at 70 au, and subsequent dip at 85 au. Our model shows another increase in the C$^{18}$O emission height at 110 au, not present in the ALMA observations, a feature created due to relative CO abundance based on our depletion profile (see Figure 1).

## 4.2. Comparison to Empirical Temperature Structure

In Law et al. (2021b), radial brightness temperatures ($T_B$) are calculated for each of the CO $J = 2 - 1$ isotopologue lines, with a resolution of a quarter of the beam size. Here, we derive the brightness temperature of our model using an identical procedure to enable a consistent comparison. We compare the empirically derived temperatures for each of the $J = 2 - 1$ lines in Figure 7. At most radii, the model's derived brightness temperature agrees within 10% of the observed brightness

<footnote n="28">²⁸ https://github.com/richteague/diskprojection</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** The seven most significant heating (right) and cooling (left) mechanisms in our final model at a disk radius of 150 au, as a function of disk height, normalized according to the highest heating/cooling value at this radius. The black dashed line corresponds to the height (40 au) at which the dust and gas temperatures become decoupled.</figcaption></figure>

temperature. Regions where the brightness model temperature is less than ~20% of the observed $T_B$ corresponds to regions where the emitting heights diverge. For example, in the $^{12}\text{CO}$ comparison, the brightness temperature derived from the model is 20% less than the observed value starting at ~125 au, precisely where the $^{12}\text{CO}$-emitting heights diverge and the model $^{12}\text{CO}$-emitting height is 11 au deeper in the disk, where temperatures are cooler. Brightness temperatures measured from $^{13}\text{CO}$ are in good agreement, up until ~250 au. Beyond ~250 au the model's emitting height sharply decreases while the observed emitting heights stay relatively constant (although with a large uncertainty), which explains the large disconnect between the model and observed brightness temperatures. The average brightness temperature derived from ALMA observations of $\text{C}^{18}\text{O}$ is 24 K, thus at radii where the model differs from the observed brightness temperature the most (~20% less than observed), it is by just 4–5 K. Based on these comparisons, we show that this temperature structure based on a model derived by matching radial intensity profiles of CO, the disk SED, and unresolved fluxes of HD and upper transitions of CO, reproduces the observed brightness temperatures relatively well.

### 4.3. Implications for Disk Mass

This model uses a disk mass of 0.14 $M_\odot$. However, there is a degeneracy between the CO and $H_2$ surface densities. Increasing one parameter while decreasing the other by the same factor produces a model that reproduces radial profiles in CO nearly identical to the original model (Calahan et al. 2021). Any change in heating based on the $H_2$ or CO mass (heating via $H_2$ formation, $H_2$ and CO self-shielding) will not significantly affect the regions where the dust and gas temperatures are coupled. This applies to the vast majority of the disk, and where the bulk of the HD flux originates. Observing the flux of the HD molecule is a way to break this degeneracy and directly probe mass, as its ratio to $H_2$ is well understood (Linsky 1998). Unfortunately, there was a non-detection of HD toward HD 163296 when observed with the PACS instrument on Herschel (Pilbratt et al. 2010; Poglitsch et al. 2010). Kama et al. (2020) provides $3\sigma$ line flux upper limits for the HD 1–0 and 2–1 lines, $6.0 \times 10^{-18} \text{ W m}^{-2}$ and $3.0 \times 10^{-18} \text{ W m}^{-2}$, respectively. Our model predicts an HD 1–0 flux of $3.3 \times 10^{-18} \text{ W m}^{-2}$ and a 2–1 flux of

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Derived emitting surfaces and calculated uncertainties of the 2–1 transitions of $^{12}\text{CO}$, $^{13}\text{CO}$, and $\text{C}^{18}\text{O}$. Observed emitting surfaces with calculated errors are shown in blue, red, and green from Law et al. (2021b) while the model emitting surfaces are shown by the brighter blue, red, and green lines without accompanying errors.</figcaption></figure>

$8.26 \times 10^{-19} \text{ W m}^{-2}$. While it remains below the flux upper limit, this results in considerable uncertainty in the mass.

To determine an upper mass limit, we use our final model and increase the disk gas mass while decreasing the CO abundance by the same factor. Due to the degeneracy between $H_2$ and CO, we can continue to reproduce the observed CO radial profiles while increasing the mass of the disk. An HD 163296 disk model with a mass of $0.35 M_\odot$ produces a predicted HD 1–0 flux of $5.9 \times 10^{-18} \text{ W m}^{-2}$ and a 2–1
<figure><img src="image_9.png" /><figcaption>**Figure 7.** A comparison of the derived brightness temperature from the HD 163296 disk model and ALMA observations (Law et al. 2021b). The y-axis shows the percentage difference between derived brightness temperature from the model and the ALMA results ($\frac{T_{\text{fuc2d}} - T_{\text{obs}}}{T_{\text{obs}}} \times 100$). The error is plotted as the shaded region and is motivated by the $1\sigma$ uncertainty in the observed brightness temperature ($\frac{T_{\text{fuc2d}} - T_{\text{obs}} \pm \Delta T_{\text{obs}}}{T_{\text{obs}}} \times 100$).</figcaption></figure>

flux of $2.0 \times 10^{-18}$ W m$^{-2}$. This is our upper mass limit, as constrained by the HD flux. This upper mass limit is higher than the majority of mass estimates for HD 163296, and lower than the largest estimate for HD 163296 which is $0.58 M_{\odot}$ from Woitke et al. (2019). A more recent estimation for HD 163296 is presented in Booth et al. (2019) using the optically thin CO isotopologue, $^{13}$C$^{17}$O, and predicts a mass of $0.21 M_{\odot}$. This is a higher mass than we use in our model, but is also consistent with the HD flux limits and our thermochemical model. Using their derived HD flux limits, Kama et al. (2020) determined an mass limit for the disk of $\le 0.067 M_{\odot}$. They created a disk model that reproduced dust observations and a number of unresolved rotational transitions of $^{12}$CO as well as a resolved observation of $^{12}$CO J=3–2. There are two significant differences between our physical models that could explain the disparity between our calculated upper mass limits. Their HD 163296 model uses a star with a luminosity of 31 $L_{\odot}$ (Folsom et al. 2012) while this model uses a much lower luminosity of 17 $L_{\odot}$ (Fairlamb et al. 2015). Additionally, they do not deplete CO which acts as a gas coolant, namely in the decoupled thermal regions, leading to a slightly cooler disk than ours. The contrast between our two models using the same HD flux observation to derive different disk-mass estimates highlights the importance of a well-defined temperature structure.

### 4.4. Implications for Gap Thermochemistry

The temperature of gas in potentially planet-carved gaps can provide an insight into planet formation models directly linking protoplanetary disk environments to planet formation theories. The effects of disk geometry on temperature have been studied previously and it has been shown that the properties of local perturbations in the disk, including gap size and depth, radial location from the star, and disk inclination, will have an effect on the temperature structure (Jang-Condell 2008). Gaps have been found to be either cooler or warmer than the surrounding medium, dependent on disk geometry and other model assumptions. A decrease in the gas and small dust surface densities exposes material in and near the midplane to more UV flux, allowing for an increase in the temperature of both the gas and dust (van der Marel et al. 2018; Alarcón et al. 2020), however puffed-up walls can produce shadows cooling the disk midplane (Nealon et al. 2019) or the gas-to-dust ratio can be low enough to decouple gas and dust temperatures heating only the dust (Facchini et al. 2018). Using our final model, we

explored the effect that corresponding gaps in the gas and small dust in HD 163296 may have on the presumed CO depletion, and subsequent gap-temperature effects.

The gaps observed in the continuum emission in the HD 163296 disk are accounted for in the surface density of the large grains in our model. Meanwhile, the gas and small dust surface densities are smooth, and do not account for these gaps. This is a deliberate decision, as a primary goal of this study is to constrain the global thermal properties of the gas disk in general. We find that even with a smooth H$_2$ gas distribution, the location of the gaps in the large-dust population are warmer than outside of the gap by an order of a few kelvin. This is supported by previous studies of the HD 163296 disk (van der Marel et al. 2018; Rab et al. 2020) which also find an increased temperature at gap locations using thermochemical models and observations of the gas and dust.

Gaps are often attributed to ongoing planet formation. To determine the gas and small-dust depletion level within the gap, we created four models representing two of the gaps and two possible planet masses. The most prevalent gaps in the continuum (the widest and highest contrast gaps) are located at 48 au and 86 au. The location, depth, and width of the gas and small-dust gap used in our model is motivated by the measured gap widths in continuum from Isella et al. (2018) and the predicted planet masses from Zhang et al. (2018) that are dependent on dust size distribution and viscosity. We set gap depths, parameterized by a depletion of H$_2$ gas within the gap, using two assumed planet masses, $1 M_J$, which represents a typical planet mass estimate for the HD 163296 disk (Pinte et al. 2018; Teague et al. 2021), and $4.45 M_J$, which represents the highest-mass planet predicted within HD 163296 (Zhang et al. 2018). We used Equation (5) in Kanagawa et al. (2015) to determine gap depth using our model’s density distribution, viscosity ($\alpha = 10^{-2}$), and planet mass. A $1 M_J$ planet at 48 au leads to a gas depletion depth of 23%, and at 86 au, 19%. A $4.45 M_J$ planet at 48 au leads to a gas depletion depth of 85%, and at 86 au, 82%.

We use our final HD 163296 disk model and rerun it with the new gas surface density dependent on gap location and planet mass. The two models with a gap from a $1 M_J$ planet produce a negligible change in the observed radial intensity profile, while a model with gap depths corresponding to a $4.45 M_J$ planet shows a significant decrease in flux, see Figure 8. We then calculate a new CO depletion profile in order for the models with a $4.45 M_J$ planet-induced gap to match the observed radial
<figure><img src="image_10.png" /><figcaption>**Figure 8.** C<sup>18</sup>O *J* = 2 – 1 radial profiles for a smooth model (black), a model with a gap in the small-dust and gas population that corresponds to a 1 *M*<sub>J</sub> planet (blue dashed), and a 4.45 *M*<sub>J</sub> planet (teal dashed). The left panel shows the results for a gap centered at 48 au, and the right panel is for a gap at 86 au.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** The calculated CO depletion for a smooth surface density model (green), and for a model with a gap at 48 au (top, pink) and 86 au (bottom, dark pink). The dashed gray lines correspond to the center of a gap, and the solid gray line corresponds to a ring. When imposing a gap in the total surface density, the CO depletion is flat or enhanced as opposed to being a local minimum at the location of the gap.</figcaption></figure>

emission profile. Those profiles are shown in Figure 9. Previously, with a smooth gas distribution, the CO depletion profile had local minima at the location of the gaps. In the case of gas and small dust surface densities with deep gaps, the CO abundance increases at the location of the gaps, bringing it to a level on par with the depletion factors just outside of the gaps. Our method of CO depletion makes it difficult to disentangle chemical/gap effects on the CO abundance, so it is impossible to discern if this relative increase in CO in the gaps is an enhancement or a leveling off to a more constant CO depletion across the radius. Studies such as Alarcón et al. (2020) and van der Marel et al. (2016) have predicted that CO enhancement in the gaps is needed to match observations. One possible mechanism to enhance the CO abundance is the presence of

<figure><img src="image_12.png" /><figcaption>**Figure 10.** A comparison of gas temperatures between a model with a smooth gas surface density, and a model that has a gap in the gas at 48 au or 86 au. Both models have gaps in the large-dust population and both models match observed radial profiles of CO, following CO depletion profiles seen in Figure 9. Contour lines show gas temperatures of the smooth gas surface density model. The background corresponds to percentage difference from the smooth surface density model. Positive values imply a hotter temperature in the gas-gap models, negative values indicate a hotter temperature in the smooth model. The dashed black lines correspond to the gap width.</figcaption></figure>

meridional flows found at the gaps in HD 163296 by Teague et al. (2019), which can transport gas and small grains from the upper atmosphere, bringing CO sublimating from grains and a rich chemistry to an otherwise chemically inert midplane. Whether or not CO is enhanced locally, is dependent on the depth of the gas gap. At present, the best method to determine gas depletion in gaps is to use kinematics to constrain the H<sub>2</sub> pressure gradients (Teague & Foreman-Mackey 2018) and thus the amount of CO chemical processing that might be happening (Alarcón et al. 2021).

Using the new CO depletion profiles for our deep gas-gap models, we find that the gap temperature increases by upwards of 10% compared to a model with a smooth gas surface density. Figure 10 shows the difference in temperature between our smooth gas model and our gapped models, both of which match the observed CO radial profiles. The increase in temperature in the gap arises due to the decreased UV opacity from the depletion of small grains, allowing more flux to enter the region. There is a region in the atmosphere above the heated layer that is cooler than a smooth gas surface density model by over 10%. This occurs at the undershoot region where the gas temperature is lower than the dust temperature due to atomic lines becoming major coolants. In the gas-gap models, this cooler, optically thick layer is at a slightly lower height.
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Emitting surface of ¹²CO, ¹³CO, and C¹⁸O J = 2 – 1 in a model that uses a depletion of CO (solid line), and a model that has a gap in the gas surface density at 86 au (dashed black line) as compared to observations (thick line).</figcaption></figure>

We also compare the observed emitting surface of the CO J = 2 – 1 line to the model with a deep gap at 86 au in Figure 11. There is limited emitting-surface information for C¹⁸O J = 2 – 1 at the location of the 48 au gap, thus we focus on the gap farthest out. The models with 1 *M*<sub>J</sub> planet gaps showed very little variation compared to the smooth model. The 4.45 *M*<sub>J</sub> model provides the best insight into how a significant deviation in the gas surface density could affect the emitting surface and the degeneracy between the H<sub>2</sub> surface density and CO abundance. At the location of the gap, each CO isotopologue 2 – 1 transition emits from a lower layer, with ¹³CO and C¹⁸O showing the strongest change in height (<6 au). This suggests that in highly gas-depleted gaps the degeneracy between the H<sub>2</sub> surface density and CO abundance can be broken. An even more significant difference between the two models is highlighted in the C¹⁸O J = 2 – 1 emitting surface. Just beyond the gap at 86 au, the model with a deep gap in the surface density follows the observed emitting surface more so than the model with a smooth surface density and structured CO depletion profile. Exploring the limits of the degeneracy between CO and H<sub>2</sub> is beyond the scope of this study, however extracting and comparing observed and simulated emitting surfaces may be an interesting tool to use in the future.

**5. Conclusions**

Using high spatial resolution observations of ¹²CO J = 2 – 1, ¹³CO J = 2 – 1, 1 – 0, C¹⁸O J = 2 – 1, 1 – 0, and C¹⁷O J = 0 from the ALMA-MAPS large program in concert with the thermochemical code RAC2D, we derive a 2D thermal structure for the disk around the Herbig Ae star HD 163296. Our conclusions are the following:

1. We derived a 2D thermal structure for the disk around the Herbig Ae star HD 163296 that reproduces six spatially resolved rotational transition lines of CO and its isotopologues, in addition to the observed SED and structures

observed in the continuum, the predicted HD flux remains below the observed upper limit, and we reproduce observed fluxes of 19 higher J-level transitions of ¹²CO.

2. The derived temperature agrees well with empirically derived temperatures and calculated emitting heights. This temperature structure represents the full 2D thermal structure, filling in temperature information at spatial locations, which observations do not directly probe.

3. We calculate a CO depletion profile which shows a relative enhancement of CO within the CO snowline, a slightly lower depletion value between the CO snowline and 250 au, and then a significant drop off that can be explained via CO chemical processing.

4. Using the derived thermal structure, we predict a midplane snowline location for CO of 60 au which corresponds to a freeze-out temperature of 18 K. We also find an upper mass limit of 0.35 *M*<sub>☉</sub>.

5. The temperature is locally enhanced in the millimeter continuum gaps at 48 and 86 au. This is true both for gaps consisting of only large-grain depletion, and those with large-grain, small-grain and gas depletion. The temperature within the gap increases slightly when gas and small dust are depleted in addition to large grains, by at most 10% in the case of significant depletion.

This paper makes use of the following ALMA data: ADS/JAO.ALMA#2018.1.01055.L. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

J.K.C. acknowledges support from the National Science Foundation Graduate Research Fellowship under grant No. DGE 1256260 and the National Aeronautics and Space Administration FINESST grant, under grant No. 80NSSC19K1534. E.A.B. acknowledges support from NSF AAG grant #1907653. K.Z. acknowledges the support of the Office of the Vice Chancellor for Research and Graduate Education at the University of Wisconsin–Madison with funding from the Wisconsin Alumni Research Foundation, and the support of NASA through Hubble Fellowship grant HST-HF2-51401.001 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. K.R.S. acknowledges the support of NASA through Hubble Fellowship Program grant HST-HF2-51419.001, awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. K.I.Ö. acknowledges support from the Simons Foundation (SCOL #321183) and an NSF AAG grant (#1907653). R.T. and F.L. acknowledge support from the Smithsonian Institution as a Submillimeter Array (SMA) Fellow. J.D.I. acknowledges support from the Science and Technology Facilities Council of the United Kingdom (STFC) under ST/T000287/1. C.J.L. acknowledges funding from the National Science Foundation Graduate Research Fellowship under grant DGE1745303. R.L.G. acknowledges support from a CNES fellowship grant. M.L.R.H. acknowledges support from the Michigan Society of Fellows. Y.A. acknowledges support by NAOJ ALMA Scientific Research grant code 2019-13B
and Grant-in-Aid for Scientific Research (S) 18H05222. S.M.A. and J.H. acknowledge funding support from the National Aeronautics and Space Administration under grant No. 17-XRP17 2-0012 issued through the Exoplanets Research Program. Support for this work was provided by NASA through the NASA Hubble Fellowship grant \#HST-HF2-51460.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. J.B. acknowledges support by NASA through the NASA Hubble Fellowship grant \#HST-HF2-51427.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. C.W. acknowledges financial support from the University of Leeds, STFC, and UKRI (grant Nos. ST/R000549/1, ST/T000287/1, MR/T040726/1). I.C. was supported by NASA through the NASA Hubble Fellowship grant HST-HF2-51405.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. F.M. acknowledges support from ANR of France under contract ANR-16-CE31-0013 (Planet-Forming-Disks) and ANR-15-IDEX-02 (through CDP "Origins of Life"). G.C. is supported by the NAOJ ALMA Scientific Research grant code 2019-13B. V. V.G. acknowledges support from FONDECYT Iniciación 11180904 and ANID project Basal AFB-170002. F.A. acknowledges support from NSF AAG grant No. 1907653. J.B.B.

acknowledges support from NASA through the NASA Hubble Fellowship grant No. HST-HF2-51429.001-A, awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. A.S.B acknowledges the studentship funded by the Science and Technology Facilities Council of the United Kingdom (STFC). A.D.B. acknowledges support from NSF AAG grant No. 1907653. H.N. acknowledges support by NAOJ ALMA Scientific Research grant code 2018-10B and Grant-in-Aid for Scientific Research 18H05441. Y.Y. is supported by IGPEES, WINGS Program, the University of Tokyo.

We also thank the referee for a careful read of the text and insightful comments.

**Facility:** ALMA

**Software:** RAC2D (Du & Bergin 2014), bettermoments (Teague & Foreman-Mackey 2018), GoFish (Teague et al. 2019), diskprojection (https://github.com/richteague/diskprojection), CASA (McMullin et al. 2007), RADMC3D (Dullemond et al. 2012), scipy (Virtanen et al. 2020), Matplotlib (Hunter 2007), Astropy (Astropy Collaboration et al. 2013, 2018)

## Appendix A
### Parameter Exploration Analysis

Our model of the HD 163296 disk based on the density structure from Z21 with an updated CO depletion profile reproduces radial profiles from all lines except $^{12}$CO $J=2-1$ relatively well (see Figure A1). Thus, this initial model appears

<figure><img src="image_14.png" /><figcaption>**Figure A1.** Radial emission profiles as predicted by the initial HD 163296 model based on Z21 compared to the observed profiles (solid line). Parameters are listed in Table 3. This model uses the CO depletion calculated in Z21, applied before the temperature calculation. The chemistry runs for 0.01 Myr and this model lacks the excess heating found necessary to reproduce the $^{12}$CO $J=2-1$ observation.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure A2.** Modeled radial emission profiles of the HD 163296 model based on Z21 compared to the observed profiles (solid line). These models exhibit a varying gas mass of 0.2 $M_\odot$, which is a mass prediction for HD 163296 by Booth et al. (2019), and 0.008 $M_\odot$, which is the smallest predicted mass in the literature. The radial profiles of models with the highest and lowest values in our parameter exploration are compared to observations.</figcaption></figure>

(The complete figure set (7 images) is available.)

to roughly represent the 2D thermal structure, but cannot reproduce the disk layers at which $^{12}$CO J = 2 − 1 emission originates, as this line is underpredicted by a factor of about 2 within the brightest region (within 50 au). However, improvements can be made to all lines except for $^{18}$O J = 2 − 1 (and arguably $^{13}$CO 1 − 0) because they are slightly underpredicted in the inner 75 au, while slightly overpredicted in the outer disk. In order to achieve a more complete and accurate thermal profile, it is worth exploring the sensitivity of the radial emission profiles to the disk physical parameters.

Throughout the parameter exploration, we did not alter the CO depletion profile. We find that the CO depletion is only degenerate with disk physical parameters that affect the total gas mass or mass distribution (γ). When the scale height, characteristic radius, or flaring parameter are changed, we find that there is no single CO depletion profile that brings all lines of CO toward what is observed.

We ran a set of models that explored the parameter space of gas mass, small-dust mass, flaring parameter (Ψ), surface density power index (γ), characteristic radius ($r_c$), critical height ($h_c$), and radial extent and allowed for the temperature to evolve in a new physical environment. Because HD 163296 has been widely studied, we determined the range exploration based on literature values. There has been a wide array of gas-mass estimates, from as low as $8 \times 10^{-3}$ up to 0.58 $M_\odot$ (Isella et al. 2007; Williams & Best 2014; Boneberg et al. 2016; Miotello et al. 2016; Williams & McPartland 2016; Booth et al. 2019; Powell et al. 2019; Woitke et al. 2019; Kama et al. 2020). Our initial gas mass is 0.14 $M_\odot$, which is on the higher end, thus we created three models that are lower in mass, and one higher in mass: $8 \times 10^{-3}$, $1 \times 10^{-2}$, $6.7 \times 10^{-2}$ (predicted gas upper limit based on the HD 1 − 0 flux

(Kama et al. 2020)), and 0.21 $M_\odot$ (the HD 163296 disk-mass estimate from Booth et al. 2019). In terms of small-dust mass, the SED constrains the range of exploration. We start with a mass of $1 \times 10^{-4} M_\odot$, and explore values above and below this with the extremes being clear under- and overpredictions of the SED flux beyond $10 \ \mu\text{m}$: $1 \times 10^{-6} M_\odot$, $1 \times 10^{-5} M_\odot$, $1 \times 10^{-3} M_\odot$, and $1 \times 10^{-2} M_\odot$. Flaring for HD 163296 has also had a wide range of estimates throughout the literature, and different studies assume both flat and flared disks. We explore models that use $\Psi$ values above and below what we have used (1.08). This includes 0.05, 1.0, 1.6, and 2.2 (Tilling et al. 2012; de Gregorio-Monsalvo et al. 2013; Woitke et al. 2019; Kama et al. 2020). The surface density index, γ, has a natural limit in our description of surface density (see Equation (1)). We explore values above and below our initial value of 0.8: 0.2, 0.6, 1.2, and 1.8. There is $^{12}$CO J = 2 − 1 detected out to 600 au, thus we only explored values above that for $r_{\text{out}}$: 700, 1000, and 1200 au. For $r_c$ and $h_c$ we explore four values, the lowest of which corresponds to the millimeter-dust distribution, and the largest being double our initially inferred value.

Changes in gas mass, small-dust mass, and scale height have similar effects on the CO radial profiles, increasing or decreasing the CO intensity with increasing/decreasing gas mass and scale height and decreasing/increasing small-dust mass. Changing the scale height will increase or decrease the CO flux along all radii for all transitions (see Figure A2), while small-dust mass and gas mass affect some lines more than others (see Figure A2). For example, changes in the mass of the gas or small-dust populations do not have as strong an effect on $^{12}$CO J 2 − 1 as they do on nearly every other line. While mass changes in these two populations appear to have similar effects, any degeneracy between the two can be broken as the small-dust population is constrained by the SED.
As the flaring parameter, $\Psi$, increases, emission is enhanced in the outer disk (the divide between “outer” and “inner” disk depends on the characteristic radius). As $\Psi$ decreases, there is a significant increase in the inner disk (see Figure A2). We found for this model that even small changes in $\Psi$ strongly affect the final radial profile across all lines. The surface density power index, $\gamma$, tends to leave the flux in the inner few astronomical units unaffected, but with a smaller $\gamma$ more emission can be found farther out in the disk (see Figure A2). This is due to the fact that a smaller $\gamma$ produces a population that is more evenly distributed. The characteristic radius affects both the height and surface density of a given population, with lower $r_c$ values producing a rapid increase in height and a turn over the surface density at a shorter radius (see Figure A2). The combination of these two effects produces radial profiles that are brighter for smaller critical radii. The outer radius values we explored did not produce significantly different CO radial profiles, due to the fact that the majority of the emission from all lines exists within 400 au, thus an outer radius cutoff well beyond this limit does not affect the observed emission (see Figure A2).

After this exploration of parameter space, and subsequently creating a number of models that altered multiple parameters simultaneously, we determined that the best set of parameters were the ones that were used by Z21. Keeping with these values, the derived thermal structure matches five out of six of the radial intensity profiles relatively well and is the best $\chi^2$ fit for the SED (as per Z21).

## Appendix B
Gas-temperature Structures in HD 163296 Models from the Literature

There has been one other recent attempt to characterize the 2D thermal structure of the HD 163296 disk specifically, using a thermochemical code, matching multiple line observations

and the SED. Those results are presented in Woitke et al. (2019). They use the thermochemical code ProDiMo (PROto-planetary DISk MOdel; Woitke et al. 2009; Kamp et al. 2010; Woitke et al. 2016) and derive a disk model that reproduces observed line fluxes from infrared to millimeter wavelengths within a factor of about 2, along with the observed SED. The model outputs from Woitke et al. (2019) are available publicly, and we compare our final thermal structure to their results in Figure B1. The dust temperatures in both models are very similar, with the disk in our model being slightly more flared. The gas temperatures are very similar within ~200 au, while past 300 au the ProDiMo model has a much hotter disk than what this study predicts. That relatively hot temperature most likely would affect the spatially resolved radial profiles, namely $^{12}$CO $J=2-1$ which emits beyond 300 au.

There are a few possible reasons as to why that study did not find it necessary to invoke additional heating in the layers at which $^{12}$CO $J=2-1$ primarily emits. A key difference between the two models is the underlying gas mass: 0.58 $M_{\odot}$ from Woitke et al. (2019) versus 0.14 $M_{\odot}$ in this study (although the ProDiMo model used a pre-Gaia distance of 119 pc, as opposed to the Gaia-determined 101 pc). Additionally, the ProDiMo model utilizes an enhanced gas/dust ratio of >100 throughout the whole disk, and an even higher ratio within the inner few astronomical units. In our case, depleting small dust within the inner disk (tens of astronomical units) in our model did not appear to rectify the difference between predicted and observed $^{12}$CO $J=2-1$. The gas temperature in the inner 4 au of the ProDiMo model is on the order of 500 K or more, which we do not see in our model, and none of our observations constrain gas temperatures at such small scales. ProDiMo specifically models the inner disk separately from the outer disk, and the temperature within this region is significantly warmer than the outer disk region.

### RAC2D: This Work

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

### ProDiMo: Woitke et al. 2019

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

**Figure B1.** A comparison between this work's HD 163296 thermal structure that best reproduces the CO radial profiles, and the HD 163296 specific model from Woitke et al. (2019). The two contours in black follow 19 and 25 K.
<table><caption>Table B1
CO Upper Level Transitions</caption><thead><tr><th>Transition</th><th>Model Flux<br>Integrated Intensity<br>(1 × 10<sup>-17</sup> W<br>m<sup>-2</sup>)</th><th>Observed Flux<br>Model in<br>Observed Range?</th></tr></thead><tbody><tr><td>5 – 4</td><td>0.504</td><td>1.04 ± 0.4</td><td>2σ</td></tr><tr><td>6 – 5</td><td>0.677</td><td>0.74 ± 0.29</td><td>1σ</td></tr><tr><td>7 – 6</td><td>0.807</td><td>0.9 ± 0.3</td><td>1σ</td></tr><tr><td>8 – 7</td><td>0.847</td><td>1.24 ± 0.55</td><td>1σ</td></tr><tr><td>9 – 8</td><td>0.766</td><td>0.91 ± 0.4</td><td>1σ</td></tr><tr><td>10 – 9</td><td>0.628</td><td>1.17 ± 0.35</td><td>2σ</td></tr><tr><td>11 – 10</td><td>0.566</td><td>1.13 ± 0.35</td><td>2σ</td></tr><tr><td>12 – 11</td><td>0.586</td><td>1.17 ± 0.35</td><td>2σ</td></tr><tr><td>13 – 12</td><td>0.606</td><td>1.52 ± 0.4</td><td>3σ</td></tr><tr><td>14 – 13</td><td>0.626</td><td>&lt;1.6</td><td>True</td></tr><tr><td>15 – 14</td><td>0.648</td><td>1.03 ± 0.5</td><td>1σ</td></tr><tr><td>16 – 15</td><td>0.668</td><td>&lt;1.3</td><td>True</td></tr><tr><td>17 – 16</td><td>0.680</td><td>0.75 ± 0.25</td><td>1σ</td></tr><tr><td>18 – 17</td><td>0.682</td><td>&lt;0.9</td><td>True</td></tr><tr><td>19 – 18</td><td>0.671</td><td>&lt;0.9</td><td>True</td></tr><tr><td>20 – 19</td><td>0.65</td><td>&lt;0.9</td><td>True</td></tr><tr><td>21 – 20</td><td>0.621</td><td>&lt;0.9</td><td>True</td></tr><tr><td>22 – 21</td><td>0.587</td><td>&lt;0.9</td><td>True</td></tr><tr><td>23 – 22</td><td>0.549</td><td>&lt;0.9</td><td>True</td></tr></tbody></table>

Note. Observed flux values are from Fedele et al. (2016) and references therein. Quoted errors are the 1σ noise measured in the continuum nearby to each line. The rightmost column shows whether the modeled flux agrees within 1σ, 2σ, or 3σ of the observation, or “True” if the modeled flux is below an upper limit.

## ORCID iDs

Jenny K. Calahan https://orcid.org/0000-0002-0150-0125

Edwin A. Bergin https://orcid.org/0000-0003-4179-6394

Ke Zhang https://orcid.org/0000-0002-0661-7517

Kamber R. Schwarz https://orcid.org/0000-0002-6429-9457

Karin I. Öberg https://orcid.org/0000-0001-8798-1347

Viviana V. Guzmán https://orcid.org/0000-0003-4784-3040

Catherine Walsh https://orcid.org/0000-0001-6078-786X

Yuri Aikawa https://orcid.org/0000-0003-3283-6884

Felipe Alarcón https://orcid.org/0000-0002-2692-7862

Sean M. Andrews https://orcid.org/0000-0003-2253-2270

Jaehan Bae https://orcid.org/0000-0001-7258-770X

Jennifer B. Bergner https://orcid.org/0000-0002-8716-0482

Alice S. Booth https://orcid.org/0000-0003-2014-2121

Arthur D. Bosman https://orcid.org/0000-0003-4001-3589

Gianni Cataldi https://orcid.org/0000-0002-2700-9676

Ian Czekala https://orcid.org/0000-0002-1483-8811

Jane Huang https://orcid.org/0000-0001-6947-6072

John D. Ilee https://orcid.org/0000-0003-1008-1142

Charles J. Law https://orcid.org/0000-0003-1413-1776

Romane Le Gal https://orcid.org/0000-0003-1837-3772

Feng Long https://orcid.org/0000-0002-7607-719X

Ryan A. Loomis https://orcid.org/0000-0002-8932-1219

François Ménard https://orcid.org/0000-0002-1637-7393

Hideko Nomura https://orcid.org/0000-0002-7058-7682

Chunhua Qi https://orcid.org/0000-0001-8642-1786

Richard Teague https://orcid.org/0000-0003-1534-5186

Merel L. R. van't Hoff https://orcid.org/0000-0002-2555-9869

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

Yoshihide Yamato https://orcid.org/0000-0003-4099-6941

## References

Alarcón, F., Bosman, A. D., Bergin, E. A., et al. 2021, ApJS, 257, 8

Alarcón, F., Teague, R., Zhang, K., Bergin, E. A., & Barraza-Alfaro, M. 2020, ApJ, 905, 68

Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJL, 869, L41

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

Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Mantelet, G., & Andrae, R. 2018, AJ, 156, 58

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

Bergin, E. A., Cleeves, L. I., Crockett, N., & Blake, G. A. 2014, FaDi, 168, 61

Bergin, E. A., Cleeves, L. I., Gorti, U., et al. 2013, Natur, 493, 644

Bergin, E. A., & Williams, J. P. 2017, in The Determination of Protoplanetary Disk Masses, ed. M. Pessah & O. Gressel (Cham: Springer), 1

Birnstiel, T., Dullemond, C. P., Zhu, Z., et al. 2018, ApJL, 869, L45

Boneberg, D. M., Panić, O., Haworth, T. J., Clarke, C. J., & Min, M. 2016, MNRAS, 461, 385

Booth, A. S., Walsh, C., Ilee, J. D., et al. 2019, ApJL, 882, L31

Bosman, A. D., Walsh, C., & van Dishoeck, E. F. 2018, A&A, 618, A182

Calahan, J., Bergin, E., Zhang, K., et al. 2021, ApJ, 908, 8

Cridland, A. J., Bosman, A. D., & van Dishoeck, E. F. 2020, A&A, 635, A68

Czekala, I., Loomis, R. A., Teague, R., et al. 2021, ApJS, 257, 2

D'Alessio, P., Calvet, N., & Woolum, D. S. 2005, in ASP Conf. Ser. 341, Chondrites and the Protoplanetary Disk, ed. A. N. Krot, E. R. D. Scott, & B. Reipurth (San Francisco, CA: ASP), 353

de Gregorio-Monsalvo, I., Ménard, F., Dent, W., et al. 2013, A&A, 557, A133

Dent, W. R. F., Pinte, C., Cortes, P. C., et al. 2019, MNRAS, 482, L29

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

Draine, B. T., & Li, A. 2001, ApJ, 551, 807

Du, F., & Bergin, E. A. 2014, ApJ, 792, 2

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

Facchini, S., Pinilla, P., van Dishoeck, E. F., & de Juan Ovelar, M. 2018, A&A, 612, A104

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

Fedele, D., van Dishoeck, E. F., Kama, M., Bruderer, S., & Hogerheijde, M. R. 2016, A&A, 591, A95

Folsom, C. P., Bagnulo, S., Wade, G. A., et al. 2012, MNRAS, 422, 2072

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1

Garrod, R. T., Widicus Weaver, S. L., & Herbst, E. 2008, ApJ, 682, 283

Garufi, A., Quanz, S. P., Schmid, H. M., et al. 2014, A&A, 568, A40

Glassgold, A. E., Najita, J., & Igea, J. 2004, ApJ, 615, 972

Grady, C. A., Devine, D., Woodgate, B., et al. 2000, ApJ, 544, 895

Hayashi, C. 1981, PThPS, 70, 35

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

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

Isella, A., Huang, J., Andrews, S. M., et al. 2018, ApJL, 869, L49

Isella, A., Testi, L., Natta, A., et al. 2007, A&A, 469, 213

Jang-Condell, H. 2008, ApJ, 679, 797

Jorsater, S., & van Moorsel, G. A. 1995, AJ, 110, 2037

Kama, M., Bruderer, S., van Dishoeck, E. F., et al. 2016, A&A, 592, A83

Kama, M., Trapman, L., Fedele, D., et al. 2020, A&A, 634, A88

Kamp, I., & Dullemond, C. P. 2004, ApJ, 615, 991

Kamp, I., Tilling, I., Woitke, P., Thi, W. F., & Hogerheijde, M. 2010, A&A, 510, A18

Kanagawa, K. D., Muto, T., Tanaka, H., et al. 2015, ApJL, 806, L15

Kenyon, S. J., & Hartmann, L. 1987, ApJ, 323, 714

Law, C. J., Loomis, R. A., & Teague, R. 2021a, ApJS, 257, 3

Law, C. J., Teague, R., & Loomis, R. A. 2021b, ApJS, 257, 4

Li, A., & Draine, B. T. 2001, ApJ, 554, 778

Linsky, J. L. 1998, SSRv, 84, 285

Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603

Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425

McClure, M. K., Bergin, E. A., Cleeves, L. I., et al. 2016, ApJ, 831, 167

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in ASP Conf. Ser. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell (San Francisco, CA: ASP), 127

Miotello, A., Bruderer, S., & van Dishoeck, E. F. 2014, A&A, 572, A96

Miotello, A., van Dishoeck, E. F., Kama, M., & Bruderer, S. 2016, A&A, 594, A85

Monnier, J. D., Harries, T. J., Aarnio, A., et al. 2017, ApJ, 838, 20

Muro-Arena, G. A., Dominik, C., Waters, L. B. F. M., et al. 2018, A&A, 614, A24
Najita, J. R., & Ádámkovics, M. 2017, ApJ, 847, 6

Nealon, R., Pinte, C., Alexander, R., Mentiplay, D., & Dipierro, G. 2019, MNRAS, 484, 4951

Öberg, K., Guzmán, V. V., Walsh, C., et al. 2021, ApJS, 257, 1

Öberg, K. I., & Bergin, E. A. 2016, ApJL, 831, L19

Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011, ApJL, 743, L16

Okuzumi, S., Momose, M., Sirono, S.-i., Kobayashi, H., & Tanaka, H. 2016, ApJ, 821, 82

Penna, R. F., Sądowski, A., Kulkarni, A. K., & Narayan, R. 2013, MNRAS, 428, 2255

Pilbratt, G. L., Riedinger J. R., Passvogel, T., et al. 2010, A&A, 518, L1

Pinte, C., Price, D. J., Ménard, F., et al. 2018a, ApJL, 860, L13

Pinte, C., Price, D. J., Ménard, F., et al. 2020, ApJL, 890, L9

Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A, 518, L2

Powell, D., Murray-Clay, R., Pérez, L. M., Schlichting, H. E., & Rosenthal, M. 2019, ApJ, 878, 116

Qi, C., D'Alessio, P., Öberg, K. I., et al. 2011, ApJ, 740, 84

Qi, C., Öberg, K. I., Andrews, S. M., et al. 2015, ApJ, 813, 128

Rab, C., Kamp, I., Dominik, C., et al. 2020, A&A, 642, A165

Rich, E. A., Wisniewski, J. P., Sitko, M. L., et al. 2020, ApJ, 902, 4

Rosenfeld, K. A., Andrews, S. M., Hughes, A. M., Wilner, D. J., & Qi, C. 2013, ApJ, 774, 16

Schoonenberg, D., & Ormel, C. W. 2017, A&A, 602, A21

Schwarz, K. R., Bergin, E. A., Cleeves, L. I., et al. 2018, ApJ, 856, 85

Seok, J. Y., & Li, A. 2017, ApJ, 835, 291

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

Sloan, G. C., Kraemer, K. E., Price, S. D., & Shipman, R. F. 2003, ApJS, 147, 379

Stevenson, D. J., & Lunine, J. I. 1988, Icar, 75, 146

Teague, R. 2019, JOSS, 4, 1632

Teague, R., Bae, J., & Aikawa, Y. 2021, ApJS, 257, 18

Teague, R., Bae, J., & Bergin, E. A. 2019, Natur, 574, 378

Teague, R., & Foreman-Mackey, D. 2018, RNAAS, 2, 173

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

Trapman, L., Miotello, A., Kama, M., van Dishoeck, E. F., & Bruderer, S. 2017, A&A, 605, A69

van der Marel, N., van Dishoeck, E. F., Bruderer, S., et al. 2016, A&A, 585, A58

van der Marel, N., Williams, J. P., & Bruderer, S. 2018, ApJL, 867, L14

van't Hoff, M. L. R., Walsh, C., Kama, M., Facchini, S., & van Dishoeck, E. F. 2017, A&A, 599, A101

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Warren, S. G., & Brandt, R. E. 2008, JGRD, 113, D14220

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

Williams, J. P., & Best, W. M. J. 2014, ApJ, 788, 59

Williams, J. P., & McPartland, C. 2016, ApJ, 830, 32

Wilson, T. L. 1999, RPPh, 62, 143

Wisniewski, J. P., Clampin, M., Grady, C. A., et al. 2008, ApJ, 682, 548

Woitke, P., Kamp, I., Antonellini, S., et al. 2019, PASP, 131, 064301

Woitke, P., Kamp, I., & Thi, W. F. 2009, A&A, 501, 383

Woitke, P., Min, M., Pinte, C., et al. 2016, A&A, 586, A103

Woodall, J., Agúndez, M., Markwick-Kemper, A. J., & Millar, T. J. 2007, A&A, 466, 1197

Zhang, K., Booth, A. S., & Law, C. J. 2021, ApJ, 257, 5

Zhang, K., Blake, G. A., & Bergin, E. A. 2015, ApJL, 806, L7

Zhang, S., Zhu, Z., Huang, J., et al. 2018, ApJL, 869, L47