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

We use observations of ¹²CO J = 2 − 1, C¹⁸O J = 2 − 1, 1 − 0, ¹³CO J = 2 − 1, 1 − 0, and C¹⁷O J = 1 − 0 toward the HD 163296 disk. For this study, we use images which have an average spatial resolution of 26 au for the 1 − 0 transitions (025) and 13 au for the 2 − 1 transitions (Öberg et al. 2021, and Table 1). The final image cubes used for this study had a spectral resolution of 200 m s⁻¹, and used a robust = 0.5 weighting. The robust = 0.5 images were used as they provided the highest absolute resolution. These images were "JvM corrected" (Jorsater & van Moorsel 1995), which refers to the method of scaling the residuals in the image cube to have identical units to the CLEAN model. This ensures that the final image used for moment zero maps, and subsequent radial profiles, have the correct units (Jy CLEAN beam⁻¹); see Czekala et al. (2021) for a detailed explanation. The MAPS program summary, along with data reduction and calibration specifics can be found in Öberg et al. (2021), and the methods of the imaging process are thoroughly described in Czekala et al. (2021). Radial intensity profiles are created for each CO line for comparison with our model results. Radial profile derivations are presented in Law et al. (2021a). The package bettermoments ²⁵ (Teague & Foreman-Mackey 2018) is used to extract the moment zero map which is then used to calculate the radial intensity profile using 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 toward HD 163296. The HD observations were made using the Photoconductor Array Camera and Spectrometer (PACS) instrument and were spatially and spectral unresolved. Kama et al. (2020) used archival data from the DIGIT program (PI N. J. Evans) to determine the upper flux limits for 15 Herbig Ae/Be disks, including HD 163296. They found an upper limit for the flux of the J = 1 − 0 line of ≤0.9 × 10⁻¹⁷ W m⁻² and for the J = 2 − 1 they determined a flux of ≤2.4 × 10⁻¹⁷ W m⁻².

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

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

$$\begin{eqnarray}&&{\rm{\Sigma }}(r)={{\rm{\Sigma }}}_{c}{\left(\displaystyle \frac{r}{{r}_{c}}\right)}^{-\gamma }\exp \left[-{\left(\displaystyle \frac{r}{{r}_{c}}\right)}^{2-\gamma }\right],\end{eqnarray} \tag{ 1 }$$

r_c is the characteristic radius at which the surface density is Σ_c /e and γ 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:

$$\begin{eqnarray}&&\rho (r,z)=\displaystyle \frac{{\rm{\Sigma }}(r)}{\sqrt{2\pi }h(r)}\exp \left[-\displaystyle \frac{1}{2}{\left(\displaystyle \frac{z}{h(r)}\right)}^{2}\right],\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&h={h}_{c}{\left(\displaystyle \frac{r}{{r}_{c}}\right)}^{{\rm{\Psi }}},\end{eqnarray} \tag{ 3 }$$

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

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

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

We adopt the reaction rates from the UMIST 2006 database (Woodall et al. 2007) for the gas-phase chemistry with additional rates considering the self-shielding of CO, H₂, H₂O, and OH, dust-grain surface chemistry driven by temperature, UV, cosmic rays, and two-body chemical reactions on dust-grain surfaces (see references given by Du & Bergin 2014). Chemical processes also provide heating or cooling to the surrounding gas. These mechanisms, along with stellar and interstellar radiation, drive the thermal gas structure. Our study explores models that account for 0.01 Myr of chemical evolution. This is a relatively short period of time compared to disk lifetimes and is motivated by an average CO-processing timescale, which is found to be on the order of 1 Myr (Bergin et al. 2014; Bosman et al. 2018; Schwarz et al. 2018). The calculated depletion profile of CO is motivated by observations, and thus accounts for any long-timescale chemical effects that have taken place in the disk's history; we did not want to further process CO by running the disk model for longer than 0.01 Myr. It is also worth noting that the exact timescale will not affect our final temperature structure nor CO column density distribution; it will only affect the CO depletion profile. CO will be created or destroyed in our chemical network over time, and the CO depletion profile acts to counter any over- or underabundance of CO being produced. While localized variations in the CO abundance may affect the gas temperature, the global temperature structure is not strongly affected by the CO abundance. Finally, at the end of a given run, we extracted the dust and gas thermal profiles.

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

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

Zoom In Zoom Out Reset image size

In Z21, depletion of CO was calculated based on a model of the disk based on previous determinations of disk parameters and the C¹⁸O J = 2 − 1 radial intensity profile. The CO was depleted throughout the disk at the start of the ray-tracing method that creates the simulated observations, i.e., after the temperature and chemical evolution of the disk. This depletion means the CO abundance is reduced from its expected value (∼10⁻⁴ relative to H₂) in layers where the dust temperature is above the CO sublimation point and the gas is self-shielded from radiation. Since CO is a significant coolant in the disk-surface layers where the dust and gas are thermally decoupled, any depletion of CO may affect the thermal structure. Thus, for this study, we recalculated the CO depletion factors. We started by applying the Z21 CO depletion to the initial CO abundance. After running the model for 0.01 Myr, there were small inconsistencies between the simulated and observed C¹⁸O J = 2 − 1 radial emission profiles (see Figure A1). We then calculated a new CO depletion profile by determining what factor of increase or decrease in CO flux would be needed to reproduce the observation at each radii (using the same radial resolution as Z21). That factor was then applied to the original CO depletion profile at the corresponding radius. The model was then run again with the new CO depletion profile. We needed to iterate this process three times, and stopped iterating once the χ² value (using stat.chisquare function from the scipy package and comparing the simulated and observed C¹⁸O J = 2 − 1 radial profile intensity values) had dropped well below that of the first attempts. The first three attempts had a χ² value of 8.06–9.72, and the final model χ² = 2.47. The final depletion profile is shown in Figure 1, and on average the newly calculated CO depletion is 33% more depleted than that derived in Z21. Most of the difference is at large radii where the profiles vary by only 16% within 200 au.

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

While our initial model with the updated CO depletion profile reproduced most of the observed CO lines, ¹²CO J = 2 − 1 was underpredicted on average by a factor of 1.6 within 100 au, and a factor of 2.2 within 35 au. We completed a thorough exploration of the parameter space including gas mass, small-dust mass, surface density power-law index (γ), flaring index (Ψ), scale height (h_c), critical radius (r_c), and outer radius (see Appendix A). However, we found no combination of parameters that could improve the ¹²CO J = 2 − 1 flux agreement while leaving the other lines consistent with observations. This result, together with the high optical depth of the ¹²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); ¹²CO J = 3 − 2 was observed to be brighter than their model, which matches the outer disk. To solve this issue they increased the gas temperature beyond the dust temperature within 80 au in their HD 163296 disk model.

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

$$\begin{eqnarray*}&&{T}_{\mathrm{new}}={T}_{\mathrm{original}}\times {\left(\frac{r}{100\,\mathrm{au}}\right)}^{-0.4}\end{eqnarray*}$$

above a constant ratio of height over radius (z/r) = 0.21 reproduces the observed ¹²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 ¹²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 ¹²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 ¹²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σ uncertainly, and five of the transitions fit within 2σ or 3σ uncertainty (see Table B1). While the flux predictions from the model all fall within 3σ, all are on the fainter end of the observed flux. Additionally, the HD J = 2 − 1 flux does not change significantly, increasing only by 7%, remaining well below the observed upper limit. A model without the additional heat produces a CO flux from the J = 11 − 10 transition that is more than twice as small as our final model. On average, for the transitions between J = 11 − 10 and J = 23 − 22, a model without the additional heat produces a CO flux over four times lower than our current model.

Zoom In Zoom Out Reset image size

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

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

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is worth noting that the C¹⁷O and C¹⁸O 1 − 0 lines are the least well-fit to the observations. We also find that these transitions originate from deep within the disk, only slightly above the midplane, and CO snow surface. In Z21 it was noted that the C¹⁷O column densities appear to be higher than those estimated from C¹⁸O (correcting for optical depth). This is also the case for an earlier ¹³C¹⁸O detection, which appears to require a higher CO column than inferred from C¹⁸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 M_⊙ using ¹³C¹⁷O J = 3 − 2, a value 50% more massive than our final model. However, if CO is enhanced near the midplane, a larger mass may not be necessary. Such an effect is not accounted for in our analysis, but, at face value, would be consistent with our results in that a locally higher abundance of CO in the midplane, perhaps inside the CO snowline, would increase the emission of optically thin tracers (such as C¹⁷O) while having a reduced effect on the emission of the more abundant isotopologues. This process has a negligible effect on the thermal structure as the dust and gas are strongly coupled in these layers.

3.3.1. CO and CO₂ Snowline

An additional observational constraint on our model is the resolved emission heights of ¹²CO, ¹³CO, and C¹⁸O J = 2 − 1 as measured in Law et al. (2021b). We calculate the emitting surfaces of ¹²CO, ¹³CO, and C¹⁸O using simulated image cubes and the same methods as applied to the observations by Law et al. (2021b). This method depends on the ability to resolve the front and back side of a given disk over multiple channels of an image cube, and assumes azimuthal symmetry and gas rotating in circular orbits (Pinte et al. 2018). The emitting height extraction tools are found in the Python package diskprojection. ²⁸ A comparison of the emitting layers of our model and observations is shown in Figure 6.

Zoom In Zoom Out Reset image size

The calculated emitting layers for both simulated and observed data from Law et al. (2021b) exist at similar heights, and agree within 1σ at most radii. No additional changes to the model were necessary to arrive at this agreement between the model and observation. Our modeled emission height for ¹²CO J = 2 − 1 reproduces what was derived using the ALMA observations up until ∼125 au; beyond this radius and up to 250 au, the model's ¹²CO emits on average 11 au lower in the disk. When looking at the ¹³CO-emitting height, the model reproduces the observed heights beyond ∼125 au. This suggests that our model does not completely represent the CO vertical distribution, as it systematically produces a lower ¹²CO J = 2 − 1 emission height beyond ∼100au. However, this fit is quite good considering the detailed physics and chemistry included in the model. The C¹⁸O emitting heights agree well with what is observed, especially within ∼100 au, where both our model and the observations pick up on some substructure, a bump at 70 au, and subsequent dip at 85 au. Our model shows another increase in the C¹⁸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).

In Law et al. (2021b), radial brightness temperatures (T_B) are calculated for each of the CO J = 2 − 1 isotopologue lines, with a resolution of a quarter of the beam size. Here, we derive the brightness temperature of our model using an identical procedure to enable a consistent comparison. We compare the empirically derived temperatures for each of the J = 2 − 1 lines in Figure 7. At most radii, the model's derived brightness temperature agrees within 10% of the observed brightness temperature. Regions where the brightness model temperature is less than ∼20% of the observed T_B corresponds to regions where the emitting heights diverge. For example, in the ¹²CO comparison, the brightness temperature derived from the model is 20% less than the observed value starting at ∼125 au, precisely where the ¹²CO-emitting heights diverge and the model ¹²CO-emitting height is 11 au deeper in the disk, where temperatures are cooler. Brightness temperatures measured from ¹³CO are in good agreement, up until ∼250 au. Beyond ∼250 au the model's emitting height sharply decreases while the observed emitting heights stay relatively constant (although with a large uncertainty), which explains the large disconnect between the model and observed brightness temperatures. The average brightness temperature derived from ALMA observations of C¹⁸O is 24 K, thus at radii where the model differs from the observed brightness temperature the most (∼20% less than observed), it is by just 4–5 K. Based on these comparisons, we show that this temperature structure based on a model derived by matching radial intensity profiles of CO, the disk SED, and unresolved fluxes of HD and upper transitions of CO, reproduces the observed brightness temperatures relatively well.

Zoom In Zoom Out Reset image size

This model uses a disk mass of 0.14 M_⊙. However, there is a degeneracy between the CO and H₂ surface densities. Increasing one parameter while decreasing the other by the same factor produces a model that reproduces radial profiles in CO nearly identical to the original model (Calahan et al. 2021). Any change in heating based on the H₂ or CO mass (heating via H₂ formation, H₂ and CO self-shielding) will not significantly affect the regions where the dust and gas temperatures are coupled. This applies to the vast majority of the disk, and where the bulk of the HD flux originates. Observing the flux of the HD molecule is a way to break this degeneracy and directly probe mass, as its ratio to H₂ is well understood (Linsky 1998). Unfortunately, there was a nondetection of HD toward HD 163296 when observed with the PACS instrument on Herschel (Pilbratt et al. 2010; Poglitsch et al. 2010). Kama et al. (2020) provides 3σ line flux upper limits for the HD 1 − 0 and 2 − 1 lines, 6.0 × 10⁻¹⁸ W m⁻² and 3.0 × 10⁻¹⁸ W m⁻², respectively. Our model predicts an HD 1 − 0 flux of 3.3 × 10⁻¹⁸ W m⁻² and a 2 − 1 flux of 8.26 × 10⁻¹⁹ W m⁻². While it remains below the flux upper limit, this results in considerable uncertainty in the mass.

To determine an upper mass limit, we use our final model and increase the disk gas mass while decreasing the CO abundance by the same factor. Due to the degeneracy between H₂ and CO, we can continue to reproduce the observed CO radial profiles while increasing the mass of the disk. An HD 163296 disk model with a mass of 0.35 M_⊙ produces a predicted HD 1 − 0 flux of 5.9 × 10⁻¹⁸ W m⁻² and a 2 − 1 flux of 2.0 × 10⁻¹⁸ W m⁻². This is our upper mass limit, as constrained by the HD flux. This upper mass limit is higher than the majority of mass estimates for HD 163296, and lower than the largest estimate for HD 163296 which is 0.58 M_⊙ 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, ¹³C¹⁷O, and predicts a mass of 0.21 M_⊙. This is a higher mass than we use in our model, but is also consistent with the HD flux limits and our thermochemical model. Using their derived HD flux limits, Kama et al. (2020) determined an mass limit for the disk of ≤0.067 M_⊙. They created a disk model that reproduced dust observations and a number of unresolved rotational transitions of ¹²CO as well as a resolved observation of ¹²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.

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

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

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

We use our final HD 163296 disk model and rerun it with the new gas surface density dependent on gap location and planet mass. The two models with a gap from a 1 M_J planet produce a negligible change in the observed radial intensity profile, while a model with gap depths corresponding to a 4.45 M_J planet shows a significant decrease in flux, see Figure 8. We then calculate a new CO depletion profile in order for the models with a 4.45 M_J planet-induced gap to match the observed radial 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 gas and small dust surface densities with deep gaps, the CO abundance increases at the location of the gaps, bringing it to a level on par with the depletion factors just outside of the gaps. Our method of CO depletion makes it difficult to disentangle chemical/gap effects on the CO abundance, so it is impossible to discern if this relative increase in CO in the gaps is an enhancement or a leveling off to a more constant CO depletion across the radius. Studies such as Alarcón et al. (2020) and van der Marel et al. (2016) have predicted that CO enhancement in the gaps is needed to match observations. One possible mechanism to enhance the CO abundance is the presence of meridional flows found at the gaps in HD 163296 by Teague et al. (2019), which can transport gas and small grains from the upper atmosphere, bringing CO sublimating from grains and a rich chemistry to an otherwise chemically inert midplane. Whether or not CO is enhanced locally, is dependent on the depth of the gas gap. At present, the best method to determine gas depletion in gaps is to use kinematics to constrain the H₂ pressure gradients (Teague & Foreman-Mackey 2018) and thus the amount of CO chemical processing that might be happening (Alarcón et al. 2021)

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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

Zoom In Zoom Out Reset image size