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

Jenny K. Calahan University of Michigan, 323 West Hall, 1085 South University
Avenue, Ann Arbor, MI 48109, USA  Edwin A. Bergin University of Michigan, 323
West Hall, 1085 South University Avenue, Ann Arbor, MI 48109, USA  Ke Zhang
NASA Hubble Fellow  Department of Astronomy, University of Wisconsin-Madison,
475 N Charter St, Madison, WI 53706  University of Michigan, 323 West Hall,
1085 South University Avenue, Ann Arbor, MI 48109, USA  Kamber R. Schwarz NASA
Hubble Fellowship Program Sagan Fellow  Lunar and Planetary Laboratory,
University of Arizona, 1629 E. University Blvd, Tucson, AZ 85721, USA  Karin
I. Öberg Center for Astrophysics |Harvard & Smithsonian, 60 Garden St.,
Cambridge, MA 02138, USA  Viviana V. Guzmán Instituto de Astrofísica, Ponticia
Universidad Católica de Chile, Av. Vicuña Mackenna 4860, 7820436 Macul,
Santiago, Chile  Catherine Walsh School of Physics and Astronomy, University
of Leeds, Leeds, UK, LS2 9JT  Yuri Aikawa Department of Astronomy, Graduate
School of Science, The University of Tokyo, Tokyo 113-0033, Japan  Felipe
Alarcón  Department of Astronomy, University of Michigan, 323 West Hall, 1085
S. University Avenue, Ann Arbor, MI 48109, USA  Sean M. Andrews Center for
Astrophysics |Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA
Jaehan Bae NASA Hubble Fellowship Program Sagan Fellow  Earth and Planets
Laboratory, Carnegie Institution for Science, 5241 Broad Branch Road NW,
Washington, DC 20015, USA  Department of Astronomy, University of Florida,
Gainesville, FL 32611, USA  Jennifer B. Bergner University of Chicago
Department of the Geophysical Sciences, Chicago, IL 60637, USA  NASA Hubble
Fellowship Program Sagan Fellow  Alice S. Booth Leiden Observatory, Leiden
University, 2300 RA Leiden, the Netherlands  School of Physics and Astronomy,
University of Leeds, Leeds, UK, LS2 9JT  Arthur D. Bosman University of
Michigan, 323 West Hall, 1085 South University Avenue, Ann Arbor, MI 48109,
USA  Gianni Cataldi National Astronomical Observatory of Japan, Osawa 2-21-1,
Mitaka, Tokyo 181-8588, Japan  Department of Astronomy, Graduate School of
Science, The University of Tokyo, Tokyo 113-0033, Japan  Ian Czekala NASA
Hubble Fellowship Program Sagan Fellow  Department of Astronomy and
Astrophysics, 525 Davey Laboratory, The Pennsylvania State University,
University Park, PA 16802, USA  Center for Exoplanets and Habitable Worlds,
525 Davey Laboratory, The Pennsylvania State University, University Park, PA
16802, USA  Center for Astrostatistics, 525 Davey Laboratory, The Pennsylvania
State University, University Park, PA 16802, USA  Institute for Computational
& Data Sciences, The Pennsylvania State University, University Park, PA 16802,
USA  Department of Astronomy, 501 Campbell Hall, University of California,
Berkeley, CA 94720-3411, USA  Jane Huang NASA Hubble Fellowship Program Sagan
Fellow  University of Michigan, 323 West Hall, 1085 South University Avenue,
Ann Arbor, MI 48109, USA  Center for Astrophysics |Harvard & Smithsonian, 60
Garden St., Cambridge, MA 02138, USA  John D. Ilee School of Physics &
Astronomy, University of Leeds, Leeds LS2 9JT, UK  Charles J. Law Center for
Astrophysics |Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA
Romane Le Gal Center for Astrophysics |Harvard & Smithsonian, 60 Garden St.,
Cambridge, MA 02138, USA  IRAP, Université de Toulouse, CNRS, CNES, UT3, 31400
Toulouse, France  IPAG, Université Grenoble Alpes, CNRS, IPAG, F-38000
Grenoble, France  IRAM, 300 rue de la piscine, F-38406 Saint-Martin d’Hères,
France  Feng Long Center for Astrophysics |Harvard & Smithsonian, 60 Garden
St., Cambridge, MA 02138, USA  Ryan A. Loomis National Radio Astronomy
Observatory, 520 Edgemont Rd., Charlottesville, VA 22903, USA  François Ménard
Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France  Hideko Nomura
Division of Science, National Astronomical Observatory of Japan, Osawa 2-21-1,
Mitaka, Tokyo 181-8588, Japan  Chunhua Qi Center for Astrophysics |Harvard &
Smithsonian, 60 Garden St., Cambridge, MA 02138, USA  Richard Teague Center
for Astrophysics |Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138,
USA  Merel L.R. van ’t Hoff University of Michigan, 323 West Hall, 1085 South
University Avenue, Ann Arbor, MI 48109, USA  David J. Wilner Center for
Astrophysics |Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA
Yoshihide Yamato Department of Astronomy, Graduate School of Science, The
University of Tokyo, Tokyo 113-0033, Japan

###### 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 Herbig Ae star HD
163296. Using the thermo-chemical code RAC2D, we derive a thermal structure
that reproduces spatially resolved ALMA observations ($\sim$0.12 $\arcsec$ (13
au) - 0.25$\arcsec$ (26 au)) of ${}^{12}\textrm{CO}$ $J$ = 2-1,
${}^{13}\textrm{CO}$ $J$ = 1-0, 2-1, $\textrm{C}^{18}\textrm{O}$ $J$ = 1-0,
2-1, and $\textrm{C}^{17}\textrm{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 time
scale shorter than typical CO gas-phase chemistry (0.01 Myr) and an enhanced
temperature near the surface layer of the the inner disk (z/r $\geq$ 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 $\textrm{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.

protoplanetary disk, astrochemistry

††journal: AJS††facilities: ALMA††software: RAC2D (Du & Bergin, 2014),
bettermoments (Teague & Foreman-Mackey, 2018), GoFish (Teague et al., 2019),
disksurf (https://github.com/richteague/disksurf) , CASA (McMullin et al.,
2007), RADMC3D (Dullemond et al., 2012)

##  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 on 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 Eup 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 (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
thermo-chemical 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
(sub-)Millimeter 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
${}^{12}\textrm{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
$\textrm{H}_{2}$. 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 at which that species
is emitting. Thus, joint modeling of multiple spatially resolved CO
isotopologues provides radial and vertical constraints on the temperature and
CO chemistry.

Table 1: ALMA Observations Summary

Transition | Frequency | Eup | Beam | PA | rms
---|---|---|---|---|---
| (GHz) | (K) | ($\arcsec\times\arcsec$) | (∘) | (mJy beam-1)
${}^{12}\textrm{CO}$ $J$=2-1 | 230.538 | 16.6 | $0.14\times 0.10$ |  $-$76.3 | 0.561
${}^{13}\textrm{CO}$ $J$=1-0 | 110.201 | 5.3 | $0.28\times 0.22$ |  $-$89.1 | 0.434
${}^{13}\textrm{CO}$ $J$=2-1 | 220.399 | 15.8 | $0.14\times 0.11$ |  $-$76.6 | 0.541
$\textrm{C}^{18}\textrm{O}$ $J$=1-0 | 109.782 | 5.3 | $0.28\times 0.22$ |  $-$88.8 | 0.449
$\textrm{C}^{18}\textrm{O}$ $J$=2-1 | 219.560 | 15.8 | $0.14\times 0.11$ |  $-$76.5 | 0.385
$\textrm{C}^{17}\textrm{O}$ $J$=1-0 | 112.359 | 5.4 | $0.28\times 0.21$ |  $-$89.5 | 0.528

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

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 Oberg et al. (2021). HD
163296 is a Herbig Ae star with a stellar mass of 2.0
$\textrm{M}_{\odot}$(Andrews et al., 2018; Wichittanakom et al., 2020)
surrounded by a massive disk ($8\times 10^{-3}$ \- $5.8\times 10^{-1}$
$\textrm{M}_{\odot}$) (Kama et al., 2020, and references therein) 101 pc away
(Gaia Collaboration et al., 2018; Bailer-Jones et al., 2018). Both millimeter
continuum and scattered light observations show 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 0.07-4.45
MJ have been theorized to exist within the gaps in HD 163296 (Zhang et al.,
2018). The temperature, both of 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 $\sim$ 2 MJ planet at 260 au (Pinte et al.,
2018a, 2020). Further analysis of the velocity structure by Teague et al.
(2019) mapped out the 3D gradient of velocity using ${}^{12}\textrm{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; Rosenfeld et al., 2013; de Gregorio-Monsalvo et al., 2013)
with a few studies utilizing thermo-chemical 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
${}^{12}\textrm{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 consists of highly resolved
($\sim$ 0.12$\arcsec$ \- 0.25$\arcsec$), and high S/N observations of
${}^{12}\textrm{CO}$ $J$ = 2-1, ${}^{13}\textrm{CO}$ $J$ = 2-1, 1-0,
$\textrm{C}^{18}\textrm{O}$ $J$ = 2-1, 1-0, and $\textrm{C}^{17}\textrm{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 our derivation
of an upper mass limit for the HD 163296 disk and explore CO/$\textrm{H}_{2}$
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 ${}^{12}\textrm{CO}$ $J$=2-1,
$\textrm{C}^{18}\textrm{O}$ $J$=2-1,
1-0, ${}^{13}\textrm{CO}$ $J$=2-1, 1-0, and $\textrm{C}^{17}\textrm{O}$
$J$=1-0 towards 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$\arcsec$)
and 13 au for the 2-1 (0.12$\arcsec$) transitions (Oberg et al. (2021), and
Table 1). The final image cubes used for this study had a spectral resolution
of 200 m s-1, 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 as the CLEAN
model. This ensures that the final image used for moment zero maps, and
subsequent radial profiles, have the correct units (Jy
$\mathrm{CLEAN\,beam}^{-1}$); see Czekala et al. (2021) for a detailed
explanation. The MAPS program summary, along with data reduction and
calibration specifics can be found in Oberg et al. (2021), and the methods of
the imaging process are throughly described in Czekala et al. (2021). Radial
intensity profiles are created for each CO line for comparison to our model
results. Radial profile derivations are presented in Law et al. (2021a). The
package bettermoments111https://github.com/richteague/bettermoments (Teague &
Foreman-Mackey, 2018) is used to extract the moment zero map which is then
used to calculate the radial intensity profile using
GoFish222https://github.com/richteague/gofish (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 $J$= 1-0 and 2-1 towards HD
163296. The HD observations were made using the 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 find an upper limit for the flux
of the $J$= 1-0 line of $\leq$ 0.9 $\times$ 10${}^{-17}~{}\rm{W}/\rm{m^{2}}$
and for the $J$= 2-1 they determine a flux of $\leq$ 2.4 $\times$
10${}^{-17}~{}\rm{W}/\rm{m^{2}}$.

###  2.2 RAC2D Physical Structure

In this study, we use the time-dependent thermo-chemical code
RAC2D333https://github.com/fjdu/rac-2d (Du & Bergin, 2014) to model the
thermal structure of the HD 163296 disk, and create simulated observations to
compare to 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 ($\leq$ µm), and
large dust grains ($\leq$ mm). The distribution of each component is described
by a global surface density distribution (Lynden-Bell & Pringle, 1974), which
is widely used in 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]},$
(1)

$r_{c}$ is the characteristic radius at which the surface density is
$\Sigma_{c}/\rm{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]},$
(2)

$h=h_{c}\left(\frac{r}{r_{c}}\right)^{\Psi},$ (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 an MRN grain distribution $n(a)\propto a^{-3.5}$
(Mathis et al., 1977). The small dust grains have radii between $5\times
10^{-3}$ \- $1\mu$m, and the large grains have radii between $5\times 10^{-3}$
\- $10^{3}\mu$m. The following description of the dust grain population used
in this study is adopted from Z21, which used both the HD 163296 SED and mm
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, non-smooth, 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 (G0=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, 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 cm 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, $\textrm{H}_{2}$, $\textrm{H}_{2}\textrm{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 time-scale chemical effects that have taken place in the disk’s
history; we 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 time scale 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 under-abundance 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, ${}^{13}\textrm{CO}$,
$\textrm{C}^{18}\textrm{O}$, $\textrm{C}^{17}\textrm{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 ${}^{13}\textrm{CO}$ and $\textrm{C}^{18}\textrm{O}$ abundances based
on ISM ratios of ${}^{12}\textrm{CO}$/${}^{13}\textrm{CO}$ = 69,
${}^{12}\textrm{CO}$/$\textrm{C}^{18}\textrm{O}$ = 557,
$\textrm{C}^{18}\textrm{O}$/$\textrm{C}^{17}\textrm{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 Oberg 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: $\sim$300 ${\rm
km\,s^{-1}}$ at $\sim$112 $\mu$m and $\sim$100 ${\rm km\,s^{-1}}$ at $\sim$51
$\mu$m (Poglitsch et al., 2010)

Table 2: Initial Chemical Abundances for Final Model

| Abundance Relative to Total
---|---
| Hydrogen Nuclei
$\textrm{H}_{2}$ | $5\times 10^{-1}$
HD | $2\times 10^{-5}$
He | 0.09
COa | see Figure 1
N | $7.5\times 10^{-5}$
$\textrm{H}_{2}\textrm{O (ice)}$ | $1.8\times 10^{-4}$
S | $8\times 10^{-8}$
Si+ | $8\times 10^{-9}$
Na+ | $2\times 10^{-8}$
Mg+ | $7\times 10^{-9}$
Fe+ | $3\times 10^{-9}$
P | $3\times 10^{-9}$
F | $2\times 10^{-8}$
Cl | $4\times 10^{-9}$

imposed radial depletion profile as shown in Figure 1

Note. — aCO starts with an abundance of $1.4\times 10^{-4}$ with an

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
$\textrm{H}_{2}$ as a function of radius. They use the difference between
$\textrm{C}^{18}\textrm{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.

Table 3: Gas and Dust Population Parameters: Best-fit Model Values

| Gas | Small Dust | Large Dust
---|---|---|---
Mass ($\textrm{M}_{\odot}$) | 0.14 | $2.6\times 10^{-4}$ | $2.4\times 10^{-3}$
$\Psi$ | 1.08 | 1.08 | 1.08
$\gamma$ | 0.8 | 0.8 | 0.1
${\rm h}_{\rm c}$ [au] | 14.5 | 14.5 | $-$
${\rm r}_{\rm c}$ [au] | 165 | 165 | $-$
$r_{\rm in}$ [au] | 0.45 | 0.45 | 0.45
$r_{\rm out}$ [au] | 600 | 600 | 240

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$\micron$, large dust grain sizes range from $5\times 10^{-3}$ \- $1\times
10^{3}$$\micron$. The large dust population does not have associated ${\rm
h}_{\rm c}$ nor ${\rm r}_{\rm c}$ because the surface density is set by
continuum observations, thus is non-smooth and not dictated by Equation 1-3

Figure 1: Multiplicative depletion factor for the initial CO abundance
(CO/H=$1.4\times 10^{-4}$) 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 thermo-
chemical model of the HD 163296 disk and is based on the
$\textrm{C}^{18}\textrm{O}$ $J$=2-1 observation, and acts as a starting point
for our CO depletion determination. The green represents the final CO
depletion, also motivated by the $\textrm{C}^{18}\textrm{O}$ $J$=2-1 flux but
is taken into account at the beginning of the chemical and thermal
calculations.

##  3 Results

Figure 2: These figures represent the effect of excess heating on the
${}^{12}\textrm{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{{}_{0}})^{-0.4}$ to calculate the excess
factor of heat in the regions above a certain z/r. We sought to find a z/r and
R0 value that would reproduce the observed ${}^{12}\textrm{CO}$ $J$ = 2-1
radial profile (solid red lines). The above is just a sample of the z/r and R0
values explored, with our final model using z/r limit of 0.21 and R0 = 100 au.
The bottom left plot shows the observed emitting surface of
${}^{12}\textrm{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-0.4, which is what we use to determine the excess
heat factor.

###  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 and the $\textrm{C}^{18}\textrm{O}$
$J$ = 2-1 radial intensity profile. The CO was depleted throughout the disk at
the start of the raytracing method which 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 ($\sim 10^{-4}$
relative to $\textrm{H}_{2}$) 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 $\textrm{C}^{18}\textrm{O}$ $J$=2-1 radial emission profiles (see
Figure 12). We then calculated a new CO depletion profile by determining what
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 the 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) comparing the simulated and observed $\textrm{C}^{18}\textrm{O}$
$J$=2-1 radial profile intensity values had dropped well below that of the
first attempts. The first three attempts had $\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 $\sim$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 Zhang model is then attributed to chemical CO depletion
mechanisms that occur over 1 Myr. In the Zhang model, fully self-shielded CO
becomes frozen out and interacts with OH frozen out onto grains and creates
CO2. In our model, we treat past CO chemistry as an initial condition. We find
it necessary to additionally deplete CO beyond $\sim$240 au in the model that
runs for only 0.01 Myr, since the pathway to convert CO and CO2 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, ${}^{12}\textrm{CO}$ $J$=2-1 was under-predicted 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 (hc), critical radius (rc), and outer radius (see
Appendix Section A). However, we found no combination of parameters that could
improve the ${}^{12}\textrm{CO}$ $J$=2-1 flux agreement while leaving the
other lines consistent with observations. This result together with the high
optical depth of the ${}^{12}\textrm{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); ${}^{12}\textrm{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 thermo-
chemical 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 equation below). There was no realistic amount
of heating within the decoupled layer that would increase the intensity of the
${}^{12}\textrm{CO}$ $J$=2-1 line to match what is observed, any excess
heating only affected emission from the inner few au. This then suggests that
the ${}^{12}\textrm{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_{\rm{new}}=T_{\rm{original}}\times(\frac{\rm{r}}{100\rm{au}})^{-0.4}$

above a constant ratio of height over radius (z/r) = 0.21 reproduces the
observed ${}^{12}\textrm{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 ${}^{12}\textrm{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 ${}^{12}\textrm{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 2-1 transition, see Section 4.1, and Law et al. (2021b)).
Higher transitions of CO, up to ${}^{12}\textrm{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$ uncertainly, and five of the transitions fit
within 2 or 3$\sigma$ uncertainty (see Table 4). 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 over two times less than our final model. On average, for the
transitions between $J$=11-10 and $J$=23-22, a model 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 beyond twice the original output from our thermo-chemical model. A likely
source of this excess heat is mechanical heating from processes such as
stellar winds or dissipation of gravitational energy from accretion through
the disk onto the star (e.g. Glassgold et al., 2004; Najita & Ádámkovics,
2017). Mechanical heating from such phenomena is expected to be prominent in
the inner disk ($<$ 10 au), and are not accounted for in RAC2D. Another
possible heating source originates from 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 12CO 2-1 emits. PAH emission towards HD
163296 has been observed as a part of the ISO/SWS atlas (Sloan et al., 2003)
and modeled by Seok & Li (2017). Their best fit model uses a PAH abundance of
8$\times$10-7 M⊕, which is 1.5$\times$ the abundance used in our RAC2D model.
In our model and assumed physical set-up, 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 ${}^{12}\textrm{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 ${}^{12}\textrm{CO}$ $J$=2-1 emits in
the thermally decoupled layers.

Figure 3: Integrated radial intensity profiles of ${}^{12}\textrm{CO}$ and its
isotopologues ${}^{13}\textrm{CO}$, $\textrm{C}^{18}\textrm{O}$, and
$\textrm{C}^{17}\textrm{O}$ as observed (solid line) and as simulated by the
thermo-chemcial 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.  Figure 4: The two-dimensional
profile of gas temperature, dust temperature, percentage difference between
the two, and CO abundance. The dashed black or white line in each plots shows
where z/r=0.21 which corresponds with the ${}^{12}\textrm{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 (grey), 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). 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.

###  3.3 Final Model

The model that best reproduces the CO radial profiles observed towards 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 $\simeq$ 0.25 the dust and gas are thermally coupled
(with the exception of the increased gas temperature added artificially).
Right above z/r $\simeq$ 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
$\tau=1$ surface. In our model, atomic oxygen dominates cooling in this
region, and cooling via gas-grain collisions are 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 $\textrm{H}_{2}$ are
the most significant heating processes, with direct heating of the gas via
X-ray radiation overtaking beyond z/r $\simeq$ 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 $\textrm{C}^{17}\textrm{O}$ and
$\textrm{C}^{18}\textrm{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 snowsurface. In Z21 it was
noted that the $\textrm{C}^{17}\textrm{O}$ column densities appear to be
higher than those estimated from $\textrm{C}^{18}\textrm{O}$ (correcting for
optical depth). This is also the case for an earlier 13C18O detection, which
appears to require higher CO column than inferred from
$\textrm{C}^{18}\textrm{O}$ $J$ = 2-1 (Z21). A similar analysis is preformed
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) estimates a disk mass of 0.21
$\textrm{M}_{\odot}$ using 13C17O $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
$\textrm{C}^{17}\textrm{O}$) while having a reduced effect on the emission of
the more abundant isotopologues. This process has negligible effect on the
thermal structure as the dust and gas are strongly coupled in these layers.

####  3.3.1 CO and CO2 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 N2H+ towards HD 163296 to determine the CO snowline as
N2H+ formation is inhibited in the presence of gas-phase CO. When CO is
frozen-out N2H+ can exist, emitting as a ring with the inner most edge
corresponding to the CO snowline. Using this method, and a model of the HD
163296 disk, they found a CO snowline at 74${}^{+7}_{-5}$ au (corrected for
pre-Gaia distance). Our midplane snowline location is in good agreement with
the N2H+ derived snowline location, especially when considering findings from
van ’t Hoff et al. (2017) that show the N2H+ column density peaks 5 au or more
beyond the midplane CO snowline. Additionally we can predict the CO2 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 note also, that
these snowlines are based on a thermo-chemical model influenced by
observations that do not directly probe the midplane.

##  4 Discussion & Analysis

###  4.1 Emitting Surface

An additional observational constraint on our model is the resolved emission
heights of ${}^{12}\textrm{CO}$, ${}^{13}\textrm{CO}$, and
$\textrm{C}^{18}\textrm{O}$ $J$ = 2-1 as measured in Law et al. (2021b). We
calculate the emitting surfaces of ${}^{12}\textrm{CO}$, ${}^{13}\textrm{CO}$,
and $\textrm{C}^{18}\textrm{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., 2018a). The emitting height
extraction tools are found in the Python package
disksurf444https://github.com/richteague/disksurf. A comparison of the
emitting layers of our model and observations are shown in Figure 6.

Figure 6: Derived emitting surfaces and calculated uncertainties of the 2-1
transitions of ${}^{12}\textrm{CO}$, ${}^{13}\textrm{CO}$ , and
$\textrm{C}^{18}\textrm{O}$. Observed emitting surfaces with calculated error
is shown in blue, red, and green from Law et al. (2021b) while the model
emitting surfaces are in the brighter blue, red, and green without accompanied
error. 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_{\rm{rac2d}}-T_{\rm{obs}}}{T_{\rm{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_{\rm{rac2d}}-T_{\rm{obs}}\pm\Delta T_{obs}}{T_{\rm{obs}}}\times
100$).

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 model and observation. Our modeled emission height for
${}^{12}\textrm{CO}$ $J$=2-1 reproduces what was derived using the ALMA
observations up until $\sim$125 au; beyond this radius and up to 250 au, the
model’s ${}^{12}\textrm{CO}$ emits on average 11 au lower in the disk. Looking
at ${}^{13}\textrm{CO}$ emitting height, beyond $\sim$125 au the model
reproduces the observed heights. This suggests that our model does not
completely represent the CO vertical distribution, as it systematically
produces a lower ${}^{12}\textrm{CO}$ $J$=2-1 emission height beyond
$\sim$100au. However, this fit is quite good considering the detailed physics
and chemistry included in the model. The $\textrm{C}^{18}\textrm{O}$ emitting
heights agree well with what is observed, especially within $\sim$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
$\textrm{C}^{18}\textrm{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 (TB) 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 temperature. Regions where the brightness model
temperature is less than $\sim$20% of the observed TB corresponds to regions
where the emitting heights diverge. For example, in the ${}^{12}\textrm{CO}$
comparison, the brightness temperature derived from the model is 20% less than
the observed value starting at $\sim$ 125 au, precisely where the
${}^{12}\textrm{CO}$ emitting heights diverge and the model
${}^{12}\textrm{CO}$ emitting height is 11 au deeper in the disk, where
temperatures are cooler. Brightness temperatures measured from
${}^{13}\textrm{CO}$ are in good agreement, up until $\sim$250 au. Beyond
$\sim$250 au the model’s emitting height sharply decreases while the observed
emitting heights stay relatively constant (although with large uncertainty),
which explains the large disconnect between the model and observed brightness
temperatures. The average brightness temperature derived from ALMA
observations of $\textrm{C}^{18}\textrm{O}$ is 24 K, thus at radii where the
model differs from the observed brightness temperature the most ($\sim$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 $\textrm{M}_{\odot}$. However, there is a
degeneracy between the CO and $\textrm{H}_{2}$ surface densities. Increasing
one parameter while decreasing the other by the same factor produces a model
that reproduces nearly identical radial profiles in CO as the original model
(Calahan et al., 2021). Any change in heating based on the $\textrm{H}_{2}$ or
CO mass (heating via $\textrm{H}_{2}$ formation, $\textrm{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 molecule
HD is a way to break this degeneracy and directly probe mass, as its ratio to
$\textrm{H}_{2}$ is well understood (Linsky, 1998). Unfortunately, there was a
non-detection of HD towards HD 163296 when observed with the PACS instrument
on Herschel (Poglitsch et al., 2010; Pilbratt 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}$ W/m2 and $3.0\times 10^{-18}$ W/m2 respectively. Our
model predicts an HD 1-0 flux of $3.3\times 10^{-18}$ W/m2 and a 2-1 flux of
8.26 $\times 10^{-19}$ W/m2. 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 $\textrm{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 $\textrm{M}_{\odot}$ produces a
predicted HD 1-0 flux of 5.9 $\times 10^{-18}$ W/m2 and 2-1 flux of 2.0
$\times 10^{-18}$ W/m2. 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
$\textrm{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}\rm{C}^{17}\rm{O}$, and predict a mass of 0.21
$\textrm{M}_{\odot}$. This is a higher mass than we use in our model, but is
also consistent with the HD flux limits and our thermo-chemical model. Using
their derived HD flux limits, Kama et al. (2020) determined an mass limit for
the disk of $\leq$ 0.067 $\textrm{M}_{\odot}$. They created a disk model that
reproduced dust observations and a number of unresolved rotational transitions
of ${}^{12}\textrm{CO}$ as well as a resolved observation of
${}^{12}\textrm{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⊙ (Folsom et al., 2012) while this model uses a much lower luminosity of 17
L⊙ (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.

Figure 8: $\textrm{C}^{18}\textrm{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 MJ planet (blue dashed), and 4.45 MJ planet (teal dashed).
The left plot shows the results for a gap centered at 48 au, and the right a
gap at 86 au.  Figure 9: The calculated CO depletion for a smooth surface
density model (green), and for a models with a gap at 48 au (top, pink) and 86
au (bottom, dark pink). The dashed grey lines correspond to the center of a
gap, and the solid grey link 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. 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, and follow a 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 hotter temperature in the smooth model. The dashed black lines
correspond to the gap width.

###  4.4 Implications for Gap Thermo-Chemistry

The temperature of gas in potentially planet-carved gaps can provide 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 density 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 $\textrm{H}_{2}$ gas distribution,
the location of the gaps in the large dust population are warmer than outside
of the gap by on 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 thermo-chemical
models and observations of the gas and dust.

Gaps are often attributed to on-going 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) dependent on dust size distribution and viscosity. We set gap
depths, parameterized by a depletion of $\textrm{H}_{2}$ gas within the gap,
using two assumed planet masses, 1 MJ, which represents a typical planet mass
estimate for the HD 163296 disk (Pinte et al., 2018b; Teague et al., 2021),
and 4.45 MJ, 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 MJ planet at 48 au leads to a gas
depletion depth of 23%, and at 86 au, 19%. A 4.45 MJ 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 re-run it with the new gas surface
density dependent on gap location and planet mass. The two models with a gap
from a 1 MJ planet produces a negligible change in the observed radial
intensity profile while a model with gap depths corresponding to 4.45 MJ
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 MJ planet induced
gap to match the observed radial emission profile. Those profiles are show 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 a gas and small
dust surface density 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
decern if this relative increase in CO in the gaps is an enhancement or a
leveling off to a more constant CO depletion across 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 meridional flows
found at the gaps in HD 163296 by Teague et al. (2019), which can transport
gas and small grains from 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 $\textrm{H}_{2}$ 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.

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 $\textrm{C}^{18}\textrm{O}$ $J$ = 2-1 at the location of the
48 au gap, thus we focus on the gap farthest out. The models with 1 MJ planet
gaps showed very little variation compared to the smooth model. The 4.45 MJ
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 $\textrm{H}_{2}$ surface density and CO abundance. At the location of the
gap, each CO isotopologue 2-1 transition emits from a lower layer, with
${}^{13}\textrm{CO}$ and $\textrm{C}^{18}\textrm{O}$ showing the strongest
change in height ($<$ 6 au). This suggests that in highly gas-depleted gaps
the degeneracy between the $\textrm{H}_{2}$ surface density and CO abundance
can be broken. An even more significant difference between the two models is
highlighted in the $\textrm{C}^{18}\textrm{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 $\textrm{H}_{2}$ 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.

Figure 11: Emitting surface of ${}^{12}\textrm{CO}$, ${}^{13}\textrm{CO}$, and
$\textrm{C}^{18}\textrm{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).

##  5 Conclusion

Using high spatial resolution observations of ${}^{12}\textrm{CO}$ $J$=2-1,
${}^{13}\textrm{CO}$ $J$=2-1, 1-0, $\textrm{C}^{18}\textrm{O}$ $J$=2-1, 1-0,
and $\textrm{C}^{17}\textrm{O}$ 1-0 from the ALMA MAPS large program in
concert with the thermo-chemical code RAC2D, we derive a 2D thermal structure
for the disk around Herbig Ae star HD 163296. Our conclusions are the
following:

  1. 1.

We derived a 2D thermal structure for the disk around Herbig Ae star HD 163296
that reproduces six spatially resolved rotational transition lines of CO and
its isotopologues, in addition to the observed SED, 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
${}^{12}\textrm{CO}$.

  2. 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. 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 of that can be explained via
CO chemical processing.

  4. 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 $\textrm{M}_{\odot}$.

  5. 5.

The temperature is locally enhanced in the mm continuum gaps at 48 and 86 au.
This is true both for gaps consisting of only large grain depletion, and those
with large, 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.

##  6 Acknowledgments

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 support of 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 grantHST-HF2-51419.001, awarded by the Space Telescope Science
Institute,which is operated by the Association of Universities forResearch 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). V.V.G. acknowledges support from FONDECYT Iniciación 11180904 and
ANID project Basal AFB-170002. F.A. acknowledges support from NSF AAG Grant
#1907653. R.T. and F.L. acknowledges support from the Smithsonian Institution
as a Submillimeter Array (SMA) Fellow. J.B.B. acknowledges support from NASA
through the NASA Hubble Fellowship grant #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). 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 Numbers 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. A.D.B. acknowledges
support from NSF AAG Grant #1907653. 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
and from the Science and Technology Facilities Council (grant numbers
ST/R000549/1 and ST/T000287/1). H.N. acknowledges support by NAOJ ALMA
Scientific Research Grant Code 2018-10B and Grant-in-Aid for Scientific
Research 18H05441. 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 Number 2019-13B. 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.

, 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}\textrm{CO}$ $J$=2-1 relatively well (see Figure 12). Thus,
this initial model appears to roughly represent the 2D thermal structure, but
cannot reproduce the disk layers at which ${}^{12}\textrm{CO}$ $J$=2-1
emission originates, as this line is underpredicted by a factor of about two
within the brightest region (within 50 au). However, improvements can be made
to all lines except for $\textrm{C}^{18}\textrm{O}$ $J$=2-1 (and arguably
${}^{13}\textrm{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 ($\gamma$).
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 towards what is observed.

We ran a set of models that explored the parameter space of gas mass, small
dust mass, flaring parameter ($\Psi$), surface density power-index ($\gamma$),
characteristic radius (rc), critical height (hc), 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 $\textrm{M}_{\odot}$ (Isella et al.,
2007; Williams & Best, 2014; Boneberg et al., 2016; Miotello et al., 2016;
Williams & McPartland, 2016; Woitke et al., 2019; Booth et al., 2019; Powell
et al., 2019; Kama et al., 2020). Our initial gas mass is 0.14
$\textrm{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 $\textrm{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}$
$\textrm{M}_{\odot}$, and explore values above and below this with the
extremes being clear under and over predictions of the SED flux beyond
10$\mu$m $\micron$: $1\times 10^{-6}$ $\textrm{M}_{\odot}$, $1\times 10^{-5}$
$\textrm{M}_{\odot}$, $1\times 10^{-3}$ $\textrm{M}_{\odot}$, $1\times
10^{-2}$ $\textrm{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, 2.2 (Tilling et al.,
2012; de Gregorio-Monsalvo et al., 2013; Woitke et al., 2019; Kama et al.,
2020). The surface density index, $\gamma$, 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}\textrm{CO}$ $J$=2-1 detected out to 600 au, thus we only explored
values above that for rout: 700, 1000, and 1200 au. For rc, hc we explore four
values, the lowest of which corresponds to the mm-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 17), while small dust mass and
gas mass effect some lines more than others (see Figures 13 and Figures 14).
For example, changes in the mass of the gas or small dust populations do not
have as strong an effect on ${}^{12}\textrm{CO}$ $J$ 2-1 as it does 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 16). 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 au unaffected, but with a smaller
$\gamma$ more emission can be found farther out in the disk (see Figure 15).
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 rc values producing a rapid
increase in height and a turn over the surface density at a shorter radius
(see Figure 18). 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 cut off well beyond this limit does not affect the
observed emission (see Figure 19).

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 (per Z21).

Figure 12: Radial emission profiles as predicted by the initial HD 163296
model based on Z21 (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}\textrm{CO}$ $J$=2-1 observation Figure 13: Modeled radial
emission profiles of HD 163296 model based on Z21 compared to the observed
profiles (solid line). These models exhibit varying gas mass: 0.2
$\textrm{M}_{\odot}$ which is a mass prediction for HD 163296 by Booth et al.
(2019) and 0.008 $\textrm{M}_{\odot}$ which is the smallest predicted mass in
literature. Figure 14: Modeled radial emission profiles of HD 163296 model
based on Z21 compared to the observed profiles (solid line). These models vary
in small dust mass: $2\times 10^{-6}$ $\textrm{M}_{\odot}$ and $2\times
10^{-2}$ $\textrm{M}_{\odot}$. The small dust mass is constrained by the SED,
and these limits represent dust mass values that are just under and over
predicting the SED. Figure 15: Modeled radial emission profiles of HD 163296
model based on Z21 compared to the observed profiles (solid line). The models
vary in surface density index, $\gamma=$0.2 and 1.8. There is a limit for
$\gamma$ of 0 - 2 based on our definition of surface density, thus we explored
values at the upper and lower ends of that natural range. Figure 16: Modeled
radial emission profiles of HD 163296 model based on Z21 compared to the
observed profiles (solid line). These models exhibit varying gas mass: 0.2
$\textrm{M}_{\odot}$ which is a mass prediction for HD 163296 by Booth et al.
(2019) and 0.008 $\textrm{M}_{\odot}$ which is the smallest predicted mass in
literature. Figure 17: Modeled radial emission profiles of HD 163296 model
based on Z21 compared to the observed profiles (solid line). These models vary
in scale height. We explore scale height values as low as what has been used
to describe the large grain population, to as high as twice as what we
originally predicted. Figure 18: Modeled radial emission profiles of HD 163296
model based on Z21 compared to the observed profiles (solid line). These
models vary in characteristic radius. We explore characteristic radius values
as low as what has been used to describe the large grain population, to as
high as twice as what we originally predicted. Figure 19: Modeled radial
emission profiles of HD 163296 model based on Z21 compared to the observed
profiles (solid line). These models vary in outer radius. We explore radius
values as low as 600 au, which is the largest radial extent of the observed
${}^{12}\textrm{CO}$, to twice that size.

##  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 thermo-chemical code,
matching multiple line observations and the SED. Those results are presented
in Woitke et al. (2019). They use the thermo-chemical code ProDiMo
(PROtoplanetary 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 two, 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 20. 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 $\sim$ 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}\textrm{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}\textrm{CO}$ $J$=2-1 primarily emits. A key difference between the two
models is the underlying gas mass; 0.58 $\textrm{M}_{\odot}$ from Woitke et
al. (2019) vs. 0.14 $\textrm{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 au . In our case, depleting small dust within the inner tens of 10 au in
our model did not appear to rectify the difference between predicted and
observed ${}^{12}\textrm{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.

Figure 20: 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 19K and 25K. Table 4:
CO upper level transitions

Transition | Model Flux | Observed Flux | Model in
---|---|---|---
| Integrated Intensity [1 $\times$ 10-17 W m-2] | Observed Range?
5-4 | 0.504 | 1.04 $\pm$ 0.4 | 2$\sigma$
6-5 | 0.677 | 0.74 $\pm$ 0.29 | 1$\sigma$
7-6 | 0.807 | 0.9 $\pm$ 0.3 | 1$\sigma$
8-7 | 0.847 | 1.24 $\pm$ 0.55 | 1$\sigma$
9-8 | 0.766 | 0.91 $\pm$ 0.4 | 1$\sigma$
10-9 | 0.628 | 1.17 $\pm$ 0.35 | 2$\sigma$
11-10 | 0.566 | 1.13 $\pm$ 0.35 | 2$\sigma$
12-11 | 0.586 | 1.17 $\pm$ 0.35 | 2$\sigma$
13-12 | 0.606 | 1.52 $\pm$ 0.4 | 3$\sigma$
14-13 | 0.626 | < 1.6 | True
15-14 | 0.648 | 1.03 $\pm$ 0.5 | 1$\sigma$
16-15 | 0.668 | < 1.3 | True
17-16 | 0.680 | 0.75 $\pm$ 0.25 | 1$\sigma$
18-17 | 0.682 | < 0.9 | True
19-18 | 0.671 | < 0.9 | True
20-19 | 0.65 | < 0.9 | True
21-20 | 0.621 | < 0.9 | True
22-21 | 0.587 | < 0.9 | True
23-22 | 0.549 | < 0.9 | True

Note. — Observed flux values are from Fedele et al. (2016) and references
therein. Quoted errors are the 1$\sigma$ noise measured in the continuum
nearby each line. The right-most column shows whether the modeled flux agrees
within 1, 2, or 3$\sigma$ of the observation, or ‘True’ if the modeled flux is
below an upper limit.

## References

  * Alarcón et al. (2020) Alarcón, F., Teague, R., Zhang, K., Bergin, E. A., & Barraza-Alfaro, M. 2020, ApJ, 905, 68, doi: 10.3847/1538-4357/abc1d6
  * Alarcón et al. (2021) Alarcón, F., Bosman, A., Bergin, E., et al. 2021, arXiv e-prints, arXiv:2109.06263.  https://arxiv.org/abs/2109.06263
  * Andrews et al. (2018) Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41, doi: 10.3847/2041-8213/aaf741
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Bailer-Jones et al. (2018) Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Mantelet, G., & Andrae, R. 2018, AJ, 156, 58, doi: 10.3847/1538-3881/aacb21
  * Benisty et al. (2010) Benisty, M., Natta, A., Isella, A., et al. 2010, A&A, 511, A74, doi: 10.1051/0004-6361/200912898
  * Bergin et al. (2014) Bergin, E. A., Cleeves, L. I., Crockett, N., & Blake, G. A. 2014, Faraday Discussions, 168, 61, doi: 10.1039/C4FD00003J
  * Bergin & Williams (2017) Bergin, E. A., & Williams, J. P. 2017, Astrophysics and Space Science Library, Vol. 445, The Determination of Protoplanetary Disk Masses, ed. M. Pessah & O. Gressel, 1, doi: 10.1007/978-3-319-60609-5_1
  * Bergin et al. (2013) Bergin, E. A., Cleeves, L. I., Gorti, U., et al. 2013, Nature, 493, 644, doi: 10.1038/nature11805
  * Birnstiel et al. (2018) Birnstiel, T., Dullemond, C. P., Zhu, Z., et al. 2018, ApJ, 869, L45, doi: 10.3847/2041-8213/aaf743
  * Boneberg et al. (2016) Boneberg, D. M., Panić, O., Haworth, T. J., Clarke, C. J., & Min, M. 2016, MNRAS, 461, 385, doi: 10.1093/mnras/stw1325
  * Booth et al. (2019) Booth, A. S., Walsh, C., Ilee, J. D., et al. 2019, ApJ, 882, L31, doi: 10.3847/2041-8213/ab3645
  * Bosman et al. (2018) Bosman, A. D., Walsh, C., & van Dishoeck, E. F. 2018, A&A, 618, A182, doi: 10.1051/0004-6361/201833497
  * Calahan et al. (2021) Calahan, J. K., Bergin, E., Zhang, K., et al. 2021, ApJ, 908, 8, doi: 10.3847/1538-4357/abd255
  * Cridland et al. (2020) Cridland, A. J., Bosman, A. D., & van Dishoeck, E. F. 2020, A&A, 635, A68, doi: 10.1051/0004-6361/201936858
  * Czekala et al. (2021) Czekala, I., Loomis, R. A., Teague, R., et al. 2021, arXiv e-prints, arXiv:2109.06188.  https://arxiv.org/abs/2109.06188
  * D’Alessio et al. (2005) D’Alessio, P., Calvet, N., & Woolum, D. S. 2005, in Astronomical Society of the Pacific Conference Series, Vol. 341, Chondrites and the Protoplanetary Disk, ed. A. N. Krot, E. R. D. Scott, & B. Reipurth, 353
  * de Gregorio-Monsalvo et al. (2013) de Gregorio-Monsalvo, I., Ménard, F., Dent, W., et al. 2013, A&A, 557, A133, doi: 10.1051/0004-6361/201321603
  * Dent et al. (2019) Dent, W. R. F., Pinte, C., Cortes, P. C., et al. 2019, MNRAS, 482, L29, doi: 10.1093/mnrasl/sly181
  * Draine (2003) Draine, B. T. 2003, ARA&A, 41, 241, doi: 10.1146/annurev.astro.41.011802.094840
  * Draine & Li (2001) Draine, B. T., & Li, A. 2001, ApJ, 551, 807, doi: 10.1086/320227
  * Du & Bergin (2014) Du, F., & Bergin, E. A. 2014, ApJ, 792, 2, doi: 10.1088/0004-637X/792/1/2
  * Dullemond et al. (2012) Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, RADMC-3D: A multi-purpose radiative transfer tool, Astrophysics Source Code Library.  http://ascl.net/1202.015
  * Facchini et al. (2018) Facchini, S., Pinilla, P., van Dishoeck, E. F., & de Juan Ovelar, M. 2018, A&A, 612, A104, doi: 10.1051/0004-6361/201731390
  * Fairlamb et al. (2015) Fairlamb, J. R., Oudmaijer, R. D., Mendigutía, I., Ilee, J. D., & van den Ancker, M. E. 2015, MNRAS, 453, 976, doi: 10.1093/mnras/stv1576
  * Fedele et al. (2016) Fedele, D., van Dishoeck, E. F., Kama, M., Bruderer, S., & Hogerheijde, M. R. 2016, A&A, 591, A95, doi: 10.1051/0004-6361/201526948
  * Folsom et al. (2012) Folsom, C. P., Bagnulo, S., Wade, G. A., et al. 2012, MNRAS, 422, 2072, doi: 10.1111/j.1365-2966.2012.20718.x
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1, doi: 10.1051/0004-6361/201833051
  * Garrod et al. (2008) Garrod, R. T., Widicus Weaver, S. L., & Herbst, E. 2008, ApJ, 682, 283, doi: 10.1086/588035
  * Garufi et al. (2014) Garufi, A., Quanz, S. P., Schmid, H. M., et al. 2014, A&A, 568, A40, doi: 10.1051/0004-6361/201424262
  * Glassgold et al. (2004) Glassgold, A. E., Najita, J., & Igea, J. 2004, ApJ, 615, 972, doi: 10.1086/424509
  * Grady et al. (2000) Grady, C. A., Devine, D., Woodgate, B., et al. 2000, ApJ, 544, 895, doi: 10.1086/317222
  * Hayashi (1981) Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35, doi: 10.1143/PTPS.70.35
  * Henning & Stognienko (1996) Henning, T., & Stognienko, R. 1996, A&A, 311, 291
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Isella et al. (2007) Isella, A., Testi, L., Natta, A., et al. 2007, A&A, 469, 213, doi: 10.1051/0004-6361:20077385
  * Isella et al. (2018) Isella, A., Huang, J., Andrews, S. M., et al. 2018, ApJ, 869, L49, doi: 10.3847/2041-8213/aaf747
  * Jang-Condell (2008) Jang-Condell, H. 2008, ApJ, 679, 797, doi: 10.1086/533583
  * Jorsater & van Moorsel (1995) Jorsater, S., & van Moorsel, G. A. 1995, AJ, 110, 2037, doi: 10.1086/117668
  * Kama et al. (2016) Kama, M., Bruderer, S., van Dishoeck, E. F., et al. 2016, A&A, 592, A83, doi: 10.1051/0004-6361/201526991
  * Kama et al. (2020) Kama, M., Trapman, L., Fedele, D., et al. 2020, A&A, 634, A88, doi: 10.1051/0004-6361/201937124
  * Kamp & Dullemond (2004) Kamp, I., & Dullemond, C. P. 2004, ApJ, 615, 991, doi: 10.1086/424703
  * Kamp et al. (2010) Kamp, I., Tilling, I., Woitke, P., Thi, W. F., & Hogerheijde, M. 2010, A&A, 510, A18, doi: 10.1051/0004-6361/200913076
  * Kanagawa et al. (2015) Kanagawa, K. D., Muto, T., Tanaka, H., et al. 2015, ApJ, 806, L15, doi: 10.1088/2041-8205/806/1/L15
  * Kenyon & Hartmann (1987) Kenyon, S. J., & Hartmann, L. 1987, ApJ, 323, 714, doi: 10.1086/165866
  * Law et al. (2021a) Law, C. J., Loomis, R. A., Teague, R., et al. 2021a, arXiv e-prints, arXiv:2109.06210.  https://arxiv.org/abs/2109.06210
  * Law et al. (2021b) Law, C. J., Teague, R., Loomis, R. A., et al. 2021b, arXiv e-prints, arXiv:2109.06217.  https://arxiv.org/abs/2109.06217
  * Li & Draine (2001) Li, A., & Draine, B. T. 2001, ApJ, 554, 778, doi: 10.1086/323147
  * Linsky (1998) Linsky, J. L. 1998, Space Sci. Rev., 84, 285
  * Lynden-Bell & Pringle (1974) Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603, doi: 10.1093/mnras/168.3.603
  * Mathis et al. (1977) Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425, doi: 10.1086/155591
  * McClure et al. (2016) McClure, M. K., Bergin, E. A., Cleeves, L. I., et al. 2016, ApJ, 831, 167, doi: 10.3847/0004-637X/831/2/167
  * McMullin et al. (2007) McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell, 127
  * Miotello et al. (2014) Miotello, A., Bruderer, S., & van Dishoeck, E. F. 2014, A&A, 572, A96, doi: 10.1051/0004-6361/201424712
  * Miotello et al. (2016) Miotello, A., van Dishoeck, E. F., Kama, M., & Bruderer, S. 2016, A&A, 594, A85, doi: 10.1051/0004-6361/201628159
  * Monnier et al. (2017) Monnier, J. D., Harries, T. J., Aarnio, A., et al. 2017, ApJ, 838, 20, doi: 10.3847/1538-4357/aa6248
  * Muro-Arena et al. (2018) Muro-Arena, G. A., Dominik, C., Waters, L. B. F. M., et al. 2018, A&A, 614, A24, doi: 10.1051/0004-6361/201732299
  * Najita & Ádámkovics (2017) Najita, J. R., & Ádámkovics, M. 2017, ApJ, 847, 6, doi: 10.3847/1538-4357/aa8632
  * Nealon et al. (2019) Nealon, R., Pinte, C., Alexander, R., Mentiplay, D., & Dipierro, G. 2019, MNRAS, 484, 4951, doi: 10.1093/mnras/stz346
  * Öberg & Bergin (2016) Öberg, K. I., & Bergin, E. A. 2016, ApJ, 831, L19, doi: 10.3847/2041-8205/831/2/L19
  * Öberg et al. (2011) Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011, ApJ, 743, L16, doi: 10.1088/2041-8205/743/1/L16
  * Oberg et al. (2021) Oberg, K. I., Guzman, V. V., Walsh, C., et al. 2021, arXiv e-prints, arXiv:2109.06268.  https://arxiv.org/abs/2109.06268
  * Okuzumi et al. (2016) Okuzumi, S., Momose, M., Sirono, S.-i., Kobayashi, H., & Tanaka, H. 2016, ApJ, 821, 82, doi: 10.3847/0004-637X/821/2/82
  * Penna et al. (2013) Penna, R. F., Sądowski, A., Kulkarni, A. K., & Narayan, R. 2013, MNRAS, 428, 2255, doi: 10.1093/mnras/sts185
  * Pilbratt et al. (2010) Pilbratt, G. L., Riedinger, J. R., Passvogel, T., et al. 2010, A&A, 518, L1, doi: 10.1051/0004-6361/201014759
  * Pinte et al. (2018a) Pinte, C., Price, D. J., Ménard, F., et al. 2018a, ApJ, 860, L13, doi: 10.3847/2041-8213/aac6dc
  * Pinte et al. (2018b) —. 2018b, ApJ, 860, L13, doi: 10.3847/2041-8213/aac6dc
  * Pinte et al. (2020) —. 2020, ApJ, 890, L9, doi: 10.3847/2041-8213/ab6dda
  * Poglitsch et al. (2010) Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A, 518, L2, doi: 10.1051/0004-6361/201014535
  * Powell et al. (2019) Powell, D., Murray-Clay, R., Pérez, L. M., Schlichting, H. E., & Rosenthal, M. 2019, ApJ, 878, 116, doi: 10.3847/1538-4357/ab20ce
  * Qi et al. (2011) Qi, C., D’Alessio, P., Öberg, K. I., et al. 2011, ApJ, 740, 84, doi: 10.1088/0004-637X/740/2/84
  * Qi et al. (2015) Qi, C., Öberg, K. I., Andrews, S. M., et al. 2015, ApJ, 813, 128, doi: 10.1088/0004-637X/813/2/128
  * Rab et al. (2020) Rab, C., Kamp, I., Dominik, C., et al. 2020, arXiv e-prints, arXiv:2008.05941.  https://arxiv.org/abs/2008.05941
  * Rich et al. (2020) Rich, E. A., Wisniewski, J. P., Sitko, M. L., et al. 2020, ApJ, 902, 4, doi: 10.3847/1538-4357/abb2a3
  * Rosenfeld et al. (2013) Rosenfeld, K. A., Andrews, S. M., Hughes, A. M., Wilner, D. J., & Qi, C. 2013, ApJ, 774, 16, doi: 10.1088/0004-637X/774/1/16
  * Schoonenberg & Ormel (2017) Schoonenberg, D., & Ormel, C. W. 2017, A&A, 602, A21, doi: 10.1051/0004-6361/201630013
  * Schwarz et al. (2018) Schwarz, K. R., Bergin, E. A., Cleeves, L. I., et al. 2018, ApJ, 856, 85, doi: 10.3847/1538-4357/aaae08
  * Seok & Li (2017) Seok, J. Y., & Li, A. 2017, ApJ, 835, 291, doi: 10.3847/1538-4357/835/2/291
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 500, 33
  * Sloan et al. (2003) Sloan, G. C., Kraemer, K. E., Price, S. D., & Shipman, R. F. 2003, ApJS, 147, 379, doi: 10.1086/375443
  * Stevenson & Lunine (1988) Stevenson, D. J., & Lunine, J. I. 1988, Icarus, 75, 146, doi: 10.1016/0019-1035(88)90133-9
  * Teague (2019) Teague, R. 2019, The Journal of Open Source Software, 4, 1632, doi: 10.21105/joss.01632
  * Teague et al. (2019) Teague, R., Bae, J., & Bergin, E. A. 2019, Nature, 574, 378, doi: 10.1038/s41586-019-1642-0
  * Teague & Foreman-Mackey (2018) Teague, R., & Foreman-Mackey, D. 2018, Research Notes of the American Astronomical Society, 2, 173, doi: 10.3847/2515-5172/aae265
  * Teague et al. (2021) Teague, R., Bae, J., Aikawa, Y., et al. 2021, arXiv e-prints, arXiv:2109.06218.  https://arxiv.org/abs/2109.06218
  * Tilling et al. (2012) Tilling, I., Woitke, P., Meeus, G., et al. 2012, A&A, 538, A20, doi: 10.1051/0004-6361/201116919
  * Trapman et al. (2017) Trapman, L., Miotello, A., Kama, M., van Dishoeck, E. F., & Bruderer, S. 2017, A&A, 605, A69, doi: 10.1051/0004-6361/201630308
  * van der Marel et al. (2016) van der Marel, N., van Dishoeck, E. F., Bruderer, S., et al. 2016, A&A, 585, A58, doi: 10.1051/0004-6361/201526988
  * van der Marel et al. (2018) van der Marel, N., Williams, J. P., & Bruderer, S. 2018, ApJ, 867, L14, doi: 10.3847/2041-8213/aae88e
  * van ’t Hoff et al. (2017) van ’t Hoff, M. L. R., Walsh, C., Kama, M., Facchini, S., & van Dishoeck, E. F. 2017, A&A, 599, A101, doi: 10.1051/0004-6361/201629452
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, doi: https://doi.org/10.1038/s41592-019-0686-2
  * Warren & Brandt (2008) Warren, S. G., & Brandt, R. E. 2008, Journal of Geophysical Research (Atmospheres), 113, D14220, doi: 10.1029/2007JD009744
  * Wichittanakom et al. (2020) Wichittanakom, C., Oudmaijer, R. D., Fairlamb, J. R., et al. 2020, MNRAS, 493, 234, doi: 10.1093/mnras/staa169
  * Williams & Best (2014) Williams, J. P., & Best, W. M. J. 2014, ApJ, 788, 59, doi: 10.1088/0004-637X/788/1/59
  * Williams & McPartland (2016) Williams, J. P., & McPartland, C. 2016, ApJ, 830, 32, doi: 10.3847/0004-637X/830/1/32
  * Wilson (1999) Wilson, T. L. 1999, Reports on Progress in Physics, 62, 143, doi: 10.1088/0034-4885/62/2/002
  * Wisniewski et al. (2008) Wisniewski, J. P., Clampin, M., Grady, C. A., et al. 2008, ApJ, 682, 548, doi: 10.1086/589629
  * Woitke et al. (2009) Woitke, P., Kamp, I., & Thi, W. F. 2009, A&A, 501, 383, doi: 10.1051/0004-6361/200911821
  * Woitke et al. (2016) Woitke, P., Min, M., Pinte, C., et al. 2016, A&A, 586, A103, doi: 10.1051/0004-6361/201526538
  * Woitke et al. (2019) Woitke, P., Kamp, I., Antonellini, S., et al. 2019, PASP, 131, 064301, doi: 10.1088/1538-3873/aaf4e5
  * Woodall et al. (2007) Woodall, J., Agúndez, M., Markwick-Kemper, A. J., & Millar, T. J. 2007, A&A, 466, 1197, doi: 10.1051/0004-6361:20064981
  * Zhang et al. (2015) Zhang, K., Blake, G. A., & Bergin, E. A. 2015, ApJ, 806, L7, doi: 10.1088/2041-8205/806/1/L7
  * Zhang et al. (2021) Zhang, K., Booth, A. S., Law, C. J., et al. 2021, arXiv e-prints, arXiv:2109.06233.  https://arxiv.org/abs/2109.06233
  * Zhang et al. (2018) Zhang, S., Zhu, Z., Huang, J., et al. 2018, ApJ, 869, L47, doi: 10.3847/2041-8213/aaf744