# The Anatomy of an Unusual Edge-on Protoplanetary Disk. II. Gas Temperature and a Warm Outer Region

C. Flores, G. Duchêne, S. Wolff, M. Villenave, K. Stapelfeldt, J. P. Williams, C. Pinte, D. Padgett, M. S. Connelley, G. van der Plas, et al.

► To cite this version:

C. Flores, G. Duchêne, S. Wolff, M. Villenave, K. Stapelfeldt, et al. The Anatomy of an Unusual Edge-on Protoplanetary Disk. II. Gas Temperature and a Warm Outer Region. The Astronomical Journal, 2021, 161, 10.3847/1538-3881/abeb1e. insu-03705312

HAL Id: insu-03705312

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

Submitted on 2 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

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

# The Anatomy of an Unusual Edge-on Protoplanetary Disk. II. Gas Temperature and a Warm Outer Region

C. Flores<fnref n="1" />, G. Duchêne<fnref n="2" />,<fnref n="3" />, S. Wolff<fnref n="4" />, M. Villenave<fnref n="3" />, K. Stapelfeldt<fnref n="5" />, J. P. Williams<fnref n="6" />, C. Pinte<fnref n="3" />,<fnref n="7" />, D. Padgett<fnref n="5" />, M. S. Connelley<fnref n="1" />, G. van der Plas<fnref n="3" />, F. Ménard<fnref n="3" />, and M. D. Perrin<fnref n="8" />

Received 2020 October 15; revised 2021 March 1; accepted 2021 March 1; published 2021 April 26

## Abstract

We present high-resolution ¹²CO and ¹³CO 2–1 ALMA observations, as well as optical and near-infrared spectroscopy, of the highly inclined protoplanetary disk around SSTC2D J163131.2–242627. The spectral type we derive for the source is consistent with a $1.2 M_{\odot}$ star inferred from the ALMA observations. Despite its massive circumstellar disk, we find little to no evidence for ongoing accretion on the star. The CO maps reveal a disk that is unusually compact along the vertical direction, consistent with its appearance in scattered light images. The gas disk extends about twice as far away as both the submillimeter continuum and the optical scattered light. CO is detected from two surface layers separated by a midplane region in which CO emission is suppressed, as expected from freeze-out in the cold midplane. We apply a modified version of the tomographically reconstructed distribution method presented by Dutrey et al. to derive the temperature structure of the disk. We find a temperature in the CO-emitting layers and the midplane of ~33 K and ~20 K at $R < 200$ au, respectively. Outside of $R > 200$ au, the disk’s midplane temperature increases to ~30 K, with a nearly vertically isothermal profile. The transition in CO temperature coincides with a dramatic reduction in the submicron and submillimeter emission from the disk. We interpret this as interstellar UV radiation providing an additional source of heating to the outer part of the disk.

*Unified Astronomy Thesaurus concepts: Protoplanetary disks (1300); T Tauri stars (1681); Circumstellar gas (238)*

## 1. Introduction

Protoplanetary disks are the birth site of planetary systems, and their structure can provide key insights on the planet formation process itself (Andrews 2020, and references therein). The radial structure of protoplanetary disks has been thoroughly studied in both gas and dust components thanks in large part to the Atacama Large Millimeter/submillimeter Array (ALMA; e.g., Andrews et al. 2018). These observations have revealed the prevalent occurrence of small-scale structures such as rings and spiral arms, which have been associated with planets and other physical mechanisms acting in the disks (Christiaens et al. 2014; Zhang et al. 2015; Dong et al. 2016; Gonzalez et al. 2017; Pinte et al. 2019). For disks that have been mapped in both continuum and line emission, it has been revealed that the submillimeter-emitting dust component is much less extended radially than the gas component (Mannings & Sargent 1997; Ansdell et al. 2018; Boyden & Eisner 2020). This is generally interpreted as the result of radial drift induced by the drag forces that the gas induces on large dust grains (e.g., Birnstiel & Andrews 2014; Cleeves et al. 2016), although differences in integrated optical depth could also be at play (Facchini et al. 2017). To fully characterize a protoplanetary disk, it is therefore important to use multiple tracers. In addition to submillimeter continuum and line emission maps, scattered light can be used to trace micron-sized grains, which are expected to be well coupled to the gas (e.g., Barrière-Fouchet et al. 2005) but whose surface brightness is independent of the local temperature.

The vertical structure of protoplanetary disks, on the other hand, has not been studied in as much detail due to a combination of optical depth and projection effects. Disks are expected to be in vertical hydrostatic equilibrium, whereby the vertical component of stellar gravity is balanced by the thermal pressure of the gas. Scattered light images confirm that small particles of dust are lifted up to a flared upper surface, (e.g., Avenhaus et al. 2018) but the height of this $\tau = 1$ surface above the midplane can be any number of scale heights depending on the local surface density. While the mass of the central star, which determines the first force, can be derived kinematically from the gas’ orbital motion (e.g., Huélamo et al. 2015; Simon et al. 2019), the pressure gradient is set by the vertical density and temperature profiles in the disk, which is not as readily accessible to observations. Interestingly, recent studies have provided model-independent methods to tomographically determine the gas temperature structure of protoplanetary disks (Dutrey et al. 2017; Pinte et al. 2018; Dullemond et al. 2020; Teague et al. 2020). These methods are applicable to spatially and kinematically resolved disks where the gas component is sufficiently optically thick so that the gas intensity can be transformed into brightness temperature via the Plank equation. These methods open new avenues to address questions related to vertical snow lines (Pinte et al. 2018), the extent of gas depletion in frozen-out regions (Dullemond et al. 2020), and the role of external sources of irradiation such as the interstellar radiation field (Dutrey et al. 2017).

While these temperature retrieval methods can be used for disks at moderate inclination (Pinte et al. 2018; Dullemond et al. 2020),

<footnote n="1">¹Institute for Astronomy, University of Hawaii at Manoa, 640 N. Aohoku Place, Hilo, HI 96720, USA; [caflores@hawaii.edu](mailto:caflores@hawaii.edu)</footnote>

<footnote n="2">²Astronomy Department, University of California, Berkeley, CA 94720, USA</footnote>

<footnote n="3">³Université Grenoble Alpes, CNRS, Institut de Planétologie et d'Astrophysique, IPAG, F-38000 Grenoble, France</footnote>

<footnote n="4">⁴Steward Observatory and the Department of Astronomy, The University of Arizona, 933 N. Cherry Avenue, Tucson, AZ 85721, USA</footnote>

<footnote n="5">⁵Jet Propulsion Laboratory, California Institute of Technology, Mail Stop 321-100, 4800 Oak Grove Drive Pasadena, CA 91109, USA</footnote>

<footnote n="6">⁶Institute for Astronomy, University of Hawaii at Manoa, 2680 Woodlawn Drive, Honolulu, HI 96822, USA</footnote>

<footnote n="7">⁷School of Physics and Astronomy, Monash University, Clayton Vic 3800, Australia</footnote>

<footnote n="8">⁸Space Telescope Science Institute, Baltimore, MD 21218, USA</footnote>
the spatial confusion introduced by projection effects can only be alleviated by kinematics under favorable circumstances. Disks at high inclinations (i ≳ 75°), often dubbed edge-on disks, offer a distinct advantage by preventing the two disk surfaces from being projected on top of one another, making the tomographic temperature reconstruction much easier (Dutrey et al. 2017). Our team has designed a method to identify new edge-on protoplanetary disks using Spitzer mapping of star-forming regions and obtaining confirmation images with the Hubble Space Telescope (HST), identifying a dozen new such systems to complement the existing sample (Stapelfeldt et al. 2014).

In this study, we focus on SSTC2D J163131.2-242627, hereafter Oph163131, a member of the Ophiuchus star-forming region (Evans et al. 2009; Hsieh & Lai 2013) located at a distance of *d* = 147 ± 3 pc (Ortiz-León et al. 2017). The central source of Oph163131 does not have a published spectrum, and therefore the stellar parameters of this source were highly uncertain prior to this study. The disk surrounding Oph163131 was unambiguously revealed to be edge-on with a radius of about 190 au in HST imaging (Stapelfeldt et al. 2014), as well as in subsequent ALMA continuum survey observations (Cox et al. 2017; Cieza et al. 2019). We have obtained dedicated high-resolution ALMA observations, which are presented in Villenave et al. (2020). In a companion paper (Wolff et al. 2021, hereafter Paper I), we simultaneously modeled the HST scattered light image and the ALMA continuum map. The latter is vertically much thinner, and the radiative transfer model demonstrates that the disk has undergone significant settling.

In this paper, we present and analyze new ALMA CO observations of the Oph13131 disk, as well as the first optical and infrared spectra of this star. In particular, we develop a slightly modified version of the tomographically reconstructed distribution (TRD) method (Dutrey et al. 2017) to account for non-perfectly edge-on disk, and we apply it to the disk around Oph163131 in both ¹²CO J = 2–1 and ¹³CO J = 2–1 isotopologues. Interferometric and spectroscopic observations of Oph163131 are presented in Section 2. Observational results such as spectral typing of the central source and dynamical mass estimates are presented in Section 3. In Section 4, we describe the method used to reconstruct the temperature structure of Oph163131 and present its temperature distribution. In Section 5, we put Oph163131 in context by comparing it with other disk sources with independently measured temperature profiles. We contrast the inner disk region’s temperature against models that reproduce the ALMA dust continuum and scattered light images. In addition, we present models that include photodissociation, photodesorption, and enhanced UV irradiation as possible mechanisms to explain the observations. In Section 6, we present our conclusions.

## 2. Observations and Data Reduction

### 2.1. ALMA Continuum and CO Maps

We observed Oph163131 with ALMA in Band 6 as part of project 2016.1.00771.S (PI: Duchêne). The spectral setup was divided into two 1.9 GHz continuum spectral windows, of rest frequency 216.5 GHz and 219.6 GHz, and three spectral windows that were set up to include the ¹²CO, ¹³CO, and C¹⁸O J = 2-1 transitions at 230.538 GHz, 220.399 GHz, and 219.56 GHz, respectively. The line spectral windows have native velocity resolutions of 0.16 km s⁻¹ for ¹³CO, and C¹⁸O

and 0.085 km s⁻¹ for ¹²CO. The source was observed using two different array configurations, specifically a compact and an extended configuration on 2017 April 25 and on 2017 July 7, respectively. The baselines ranged from 15 m to 2650 m between the two configurations, and the total observing time on source was 40 minutes. The raw data were calibrated using the Common Astronomy Software Applications (CASA) pipeline version 4.7.2.

To maximize the dynamical range, we performed phase self-calibration from the continuum data from both configurations. We extracted the continuum images using the CASA `tclean` task on the combined data, with a Briggs parameter of 0.5. The resulting beam size is 0.′23 × 0.′13. The continuum image is further analyzed in Paper I. After applying the continuum self-calibration solutions to all spectral windows, we extracted the emission lines from the calibrated visibilities by subtracting continuum emission using the `uvcontsub` function in CASA. In addition to the continuum subtracted data, we produced a non-continuum subtracted data set to test whether our modeling method is affected by possible over-subtraction of the continuum in optically thick regions of the disk. We binned the native spectral resolution to 0.25 km s⁻¹ and used a Briggs parameter of 0.5 to create higher signal-to-noise ratio (S/N) line images. The resulting beams were 0.′25 × 0.′20 for ¹²CO and 0.′26 × 0.′20 for ¹³CO, and C¹⁸O. The S/N of the C¹⁸O isotopologue was not high enough to produce a data set we could analyze.

### 2.2. JCMT CO Map

On 2011 April 16, we obtained ¹²CO J = 3–2 data from the Heterodyne Array Receiver Program (HARP) at the 15 m James Clerk Maxwell telescope (JCMT). An 8' × 8' map around Oph 163131 was acquired at a rest frequency of 345 GHz using 8192 Auto Correlation Spectral Imaging System channels with a bandwidth coverage of 250 MHz at a spectral resolution of 30 kHz (0.025 km s⁻¹). The observations were performed in good weather conditions with measured zenith opacities at 225 GHz (τ<sub>225 GHz</sub>) varying from 0.09 to 0.1.

The HARP data was reduced using the `SMURF` code, and a data cube was created with the `MAKECUBE` task. The pixel size of the data cube was set to 6", which oversamples the 14" JCMT main beam size (at 345 GHz). In order to increase the S/N of the data, the native velocity resolution was binned to 0.25 km s⁻¹. At this velocity resolution, the measured antenna temperature noise is σ<sub>TA*</sub> = 0.3 K per 0.25 km s⁻¹, or equivalently σ<sub>Jy</sub> = 9 Jy per 0.25 km s⁻¹. This conversion was performed using S(Jy) = 15.6 T<sub>A</sub>(K*)/η<sub>a</sub>, where η<sub>a</sub> = 0.52 is the aperture efficiency of HARP.<sup>9</sup>

### 2.3. SOAR Optical Spectroscopy

We obtained an optical spectrum of Oph163131 on 2013 May 19, with the Goodman spectrograph installed on the Southern Astrophysical Research (SOAR) 4.1 m telescope at Cerro Tololo. Two 900 second spectra were acquired while dithering along the 0.′46 slit. The object was observed at an airmass of 1.01. We used the 400 lines/millimeter grating and the red camera, which yields a spectral resolution of ~1850 over the full 3400–7500 Å wavelength range. Spectral images were reduced using the IRAF `ccdproc` package. A narrow

<footnote n="9">⁹ http://www.eaobservatory.org/jcmt/instrumentation/heterodyne/calibration/</footnote>
feature in the quartz flat-field lamp was removed using procedures within the IRAF `twodspec` package, as customized by SOAR observatory staff (“Removing arc from GOODMAN 600b flatfield”; T. Ribeiro 2021, private communication). Spectral extraction and wavelength calibration referencing SOAR mercury/argon arcs were performed using IRAF `doslit` procedures.

## 2.4. IRTF Near-infrared Spectroscopy

Our near-infrared (NIR) observations of Oph163131 were carried out with the 3.0 m NASA Infrared Telescope Facility (IRTF) on Maunakea, Hawaii on 2019 May 16. We used the upgraded SpeX instrument in the short crossdispersed (SXD) mode, which covers the 0.7–2.5 µm region in each exposure (Rayner et al. 2003), and the 0.′′5 wide slit, which gives a resolving power of $R = 1200$. Our target was nodded along the slit, with two exposures taken at each nod position. The total exposure time was driven by our goal to get a spectrum with S/N > 30 in the *H* band. The individual exposure times were limited to three minutes to ensure that the telluric emission lines would cancel when consecutive images taken at alternate nod positions were subtracted. After observing Oph163131, we acquired an A0 telluric standard star within 0.1 air masses for telluric correction purposes.

A thorium-argon lamp was observed for wavelength calibration, and a quartz lamp spectrum was acquired for flat-fielding. The SpeX data were flat-fielded, extracted, and wavelength-calibrated using Spextool (Cushing et al. 2004). After extraction and wavelength calibration, the individual extracted spectra were co-added using `xcombspec`. `Xtell-cor` was then used to construct a telluric correction model using the observed A0 standard, after which the observed spectrum of the target was divided by the telluric model. Finally, `xmergeorders` was used to combine the spectra in the separate orders into one continuous spectrum.

Although the SXD mode of SpeX provides a wavelength coverage down to 0.7 µm, the low S/N at wavelengths shortward of about 1 µm precludes connecting the NIR spectrum of Oph163131 with its optical counterpart.

# 3. Observational Results

## 3.1. Spectral Characterization

To the best of our knowledge, the optical and NIR spectra of Oph163131 presented here are the first of this source. In edge-on disks, the smooth scattering properties of astronomical dust preserve the spectral features associated with the central source (Appenzeller et al. 2005). Therefore, we use our SOAR and SpeX spectra to (1) estimate the spectral type, and hence the physical properties of the central star, and (2) evaluate the accretion status of the system based on the strength of its emission lines.

In the absence of spectrophotometric calibration and since scattering off the disk can alter the apparent colors of the central source for our optical spectrum, we focus our spectral analysis on narrow spectral features rather than on the global shape of the continuum. After extracting the 4700–7500 Å spectral range, where the signal-to-noise is sufficient, we fit the spectrum with a sixth-order polynomial and divide it by the spectrum. Polynomial orders lower than four or higher than seven leave small unwanted residuals, which do not strongly impact the stellar classification but look artificial in the stellar spectrum. The resulting continuum-normalized spectrum is

shown in the top left panel of Figure 1. The most prominent feature in the spectrum is the sharp Na I D absorption line at 5890 Å and the complex of lines surrounding the 5150 Å Mg line. The main TiO bands (at around 5900, 6200, 6800, and 7100 Å), which are typically seen in K5 and later-type stars, are weak in the spectrum of Oph 163131. The combination of these observations suggests that, in the optical, Oph163131 is an early- to mid-K spectral-type source (Jacoby et al. 1984; Allen & Strom 1995). In the NIR, the spectrum of Oph 163131 (also shown in Figure 1) shows a number of atomic features across the *J*, *H*, and *K* bands, as well as modest CO bandheads beyond 2.29 µm. The relative strengths of these features, when compared to the spectral library of Rayner et al. (2009), as well as the overall shape of the continuum, are consistent with a mid- to late-K-type star.

We further this analysis by comparing the spectrum of Oph 163131 to that of templates from the Medium resolution Isaac Newton Telescope (INT) Library of Empirical Spectra (MILES) and Stellar Library (STELIB; in the optical; Le Borgne et al. 2003; Sánchez-Blázquez et al. 2006; Falcón-Barroso et al. 2011) and IRTF/SpeX (in the NIR; Rayner et al. 2009) libraries. We perform the comparison with both dwarfs and giants.<fnref n="10" /> To circumvent issues such as stellar rotation and line veiling induced by accretion, we limit our analysis to a visual inspection that focuses equally on the strength of individual absorption features and on small-scale continuum structure related to molecular bands. This is sufficient for the level of precision needed for our analysis. The closest match for the optical spectrum of Oph 163131 is a K3–K4 giant, while in the infrared, the spectrum resembles a K4–K5 dwarf. Such discrepant spectral-type measurement at different wavelengths could be simply due to the visual characterization of the stellar spectra. However, it might also be a real phenomenon, as spectroscopic studies at high resolution have shown that starspots on the surfaces of young stars produce an observed temperature dependency with wavelength (see, e.g., Sokal et al. 2018; Flores et al. 2020). The stellar luminosity (and hence the luminosity classification) of Oph163131 is uncertain because of dilution by an unknown amount through scattering off the disk. Therefore, it cannot be calculated using the distance to the star and photometry. In any case, since our best-match luminosity type lies between dwarfs (V) and giants (III), we adopted a subgiant (IV) luminosity classification for Oph 163131. This is consistent with previous work where subgiant surface gravities have been measured for young stars (Doppmann et al. 2005; Flores et al. 2019). In summary, and considering the optical and IR spectra of the star, we classify Oph163131 as a K4 IV star, with a spectral-type uncertainty of about one subtype. This implies a stellar mass of ~1 $M_\odot$ based on recent dynamical mass measurements for stars in Taurus and Ophiuchus (Simon et al. 2019) and also from stellar evolutionary models (e.g., Somers & Pinsonneault 2015; Feiden 2016).

The optical spectrum of Oph 163131 stands out from that of other edge-on disks systems for the absence of strong emission lines (e.g., Appenzeller et al. 2005). Indeed, we only detect a very weak $H\alpha$ emission component, with an equivalent width of 0.6 ± 0.1 Å, superimposed on a significant photospheric absorption feature. Such a weak $H\alpha$ line, especially in the absence of other accretion- or outflow-related emission lines, generally indicates that the central source is currently not accreting at an

<footnote n="10">¹⁰ Unfortunately, neither library contains spectra of subgiants.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Spectrum of Oph 163131 in the optical (top left panel) and NIR (other panels) compared to dwarfs and giants from the MILES, STELIB, and SpeX libraries. In the top left panel, a known telluric feature, weak TiO bandheads, and Mg and Na features are labeled under the top, middle, and bottom spectra, respectively. The spectra are offset vertically for clarity. The optical spectrum is normalized to a low-order polynomial fit to the continuum whereas the NIR spectra are spectrophotometrically calibrated.</figcaption></figure>

appreciable level (e.g., Barradoy Navascués & Martín 2003; Duchêne et al. 2017), i.e., Oph163131 can be classified as a weak-line T Tauri star. This assessment is confirmed at first glance by our NIR spectrum, in which no significant Pa γ, Pa β, or Br γ emission is detected (see Figure 2); although, it is worth pointing out that these lines are less sensitive to low levels of accretion onto T Tauri stars (e.g., Folha & Emerson 2001; Alcalá et al. 2014).

The only other significant emission line detected in our spectra is the 1.08 µm He I line, with an equivalent width of 2.0 ± 0.2 Å. This line has a complex origin, with possible contributions from chromospheric activity, disk and/or stellar wind, and accretion onto the central star (Edwards et al. 2006). Surprisingly, all T Tauri stars with an He I emission line as strong as Oph163131 are accretors (see Figure 3), leading to an apparent paradox. The ultimate diagnostic of the origin of the He I line lies in the line profile (Fischer et al. 2008; Thanathibodee et al. 2018), which our low-resolution data do not allow us to analyze. We are therefore left with an ambiguous classification for Oph163131. We note that the edge-on nature of the object could lead to biased estimates of line equivalent width if the emission is spatially distinct from the star itself. However, if the emission is chromospheric in nature, its origin is at the stellar surface, and thus our equivalent width estimates would be unbiased.

In summary, while Oph163131 is unlikely to be accreting at a high rate, existing data are consistent with both: (1) an object

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Main emission lines in the spectrum of Oph 163131. From left to right, the lines are Hα, He I 1.08 µm and H Pa γ, H Pa β, and H Br γ. The wavelength of each feature is indicated by a vertical bar above the spectrum.</figcaption></figure>

with a modest amount of accretion (and a remarkably weak Hα line), and (2) a non-accreting, chromosphere-dominated emission line spectrum (albeit with an unusually strong He I 1.08 µm line). Only higher spectral resolution observations can solve this ambiguity.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Comparison between the equivalent width of the Hα and He I 1.08 µm lines among accreting (red diamonds) and non-accreting (blue triangles) T Tauri stars. Data are from Edwards et al. (2006), Manara et al. (2013), and Alcalá et al. (2014). Arrows mark the upper limits on the He I line. The black star indicates the results for Oph163131.</figcaption></figure>

## 3.2. Dynamical Mass Measurement

In the previous section, we used the spectral type derived for Oph163131 and results from stellar evolutionary models to infer the mass of the star. However, it is possible to directly measure the stellar mass using the disk's CO data and the "position-velocity (PV) diagram" technique (see, for example, Huélamo et al. 2015).

In general, dynamical mass measurements using this technique work better when the emission from a disk is optically thin (e.g., $^{13}$CO or $^{18}$O). This is because optically thin molecules trace the emission closer to the disk's midplane, or equivalently because there is less confusion with emission arising high up in the disk's surface. In the case of highly inclined disks $i \gtrsim 85^{\circ}$, however, the position of the disk's major axis approximately coincides with the disk midplane. Therefore low-elevation emission can be observed even for optically thick tracers such as $^{12}$CO. Since Oph163131 is known to host an almost perfectly edge-on disk, we can confidently calculate the dynamical mass of the star from either of the CO tracers.

The computed PV diagrams in $^{12}$CO and $^{13}$CO along the major axis of the disk are shown in Figure 4. In each case, Keplerian velocity profiles of an inclined disk with $i = 87^{\circ}$ (the best-fit inclination to the ALMA continuum map from Paper I) with a stellar mass of $1.2 M_{\odot}$ were overlaid. The systemic velocity of the disk was estimated to be $v_{\text{sys}} = -4.4 \text{ km s}^{-1}$, which is further confirmed in Section 3.4 using $^{12}$CO integrated spectrum and channel maps. Although the Keplerian profiles were matched to the data by eye, we are confident in this technique, as we have tested it with synthetic disk models computed with MCFOST (Pinte et al. 2006, 2009). We assign conservative upper and lower mass limits to Oph163131 of $1.4 M_{\odot}$ and $1.0 M_{\odot}$, respectively. We conclude that the dynamical mass measured for Oph163131 is in agreement with the spectral type derived in Section 3.1.

## 3.3. Cloud Contamination

A crucial part of the temperature reconstruction method explained in Section 4 requires us to determine whether there is emission on a scale much larger than the disk, for example, caused by a molecular cloud at the position of Oph163131. To

investigate whether this large-scale emission is present or not, we use $^{12}$CO $J=3-2$ single-dish observations acquired with HARP/JCMT (see Section 2.2).

A $\sim 6' \times 6'$ JCMT CO map centered at the position of Oph163131 is shown in Figure 5. Part of the large L1689 molecular cloud is detected E of Oph163131 at velocities ranging from 3.3 to $4.6 \text{ km s}^{-1}$. The spectrum measured at the location of Oph163131 shows no significant emission above $2\sigma$ in any of the velocity channels. Conversely, the three spectra measured in different locations of the molecular cloud show prominent emission at $\sim 4 \text{ km s}^{-1}$. This velocity range does not overlap with the CO emission from the Oph163131 disk in our ALMA data (see Figure 6). Therefore, we conclude from the JCMT data that the cloud emission does not coincide in the spatial or velocity dimensions with that of the disk.

To further investigate whether low levels of cloud emission could affect the disk emission while being undetected in our JCMT data, we consider the $^{12}$CO and $^{13}$CO ALMA spectra of Oph163131 (Figure 6). The spectra were computed by summing the flux in a rectangular region of $5''4$ along the major axis and $1''5$ along the minor axis of the disk. The $^{12}$CO emission is detected above a $3\sigma$ level between $-9.75$ and $+1.0 \text{ km s}^{-1}$, while the $^{13}$CO emission is detected above a $3\sigma$ level in the slightly narrower range from $-8.0$ to $-0.75 \text{ km s}^{-1}$. This difference is likely due to the lower S/N detection of the latter isotopologue.

Both spectra show the characteristic double peak profile seen in other disks (e.g., Beckwith & Sargent 1993; Mannings & Sargent 1997). The nearly perfectly symmetric profile displayed in $^{12}$CO is a strong indicator that no emission is added or subtracted to the disk. In cases where the line emission from a disk is contaminated by either a background or foreground cloud, the observed spectra always display a strongly asymmetric double profile in optically thick tracers such as $^{12}$CO (see for example, Perez et al. 2015; Guilloteau et al. 2016).

## 3.4. Disk Morphology

We now describe the morphology of the Oph163131 disk as seen by ALMA. Figure 7 shows the $^{12}$CO $J=2-1$ and $^{13}$CO $J=2-1$ channel maps of the disk. In both isotopologues, Oph163131 exhibits the typical Keplerian velocity pattern of a highly inclined disk (see, for example, Dutrey et al. 2017), with $^{12}$CO having a stronger detection due to the higher peak S/N of $\sim 24$ per $0.25 \text{ km s}^{-1}$ compared to the peak S/N for $^{13}$CO of only S/N $\sim 12$ per $0.25 \text{ km s}^{-1}$. Oph163131 is spatially resolved in both isotopologues with the $^{12}$CO extending $4''6$ along the major axis and $1''2$ in the minor axis direction, and $^{13}$CO extending $3''6$ and $0''7$ along the major and minor axes, respectively.

In $^{12}$CO, the disk displays two clear bright layers along the minor axis (in the SE and NW directions) separated by a darker middle lane. The projected distance between the two bright layers is $\sim 0''45$, and the SE layer is $10-20 \text{ mJy}$ brighter than the NW layer. Other ALMA results where CO maps show this type of brightness asymmetry are the Flying Saucer (2MASS J16281370-2431391; Dutrey et al. 2017) and IM Lup (Pinte et al. 2018). The physical explanation for this phenomenon is the viewing angle of the observer with respect to the disk. Oph163131 is tilted slightly away from a perfectly edge-on configuration, with the SE side being the "top" side as seen from the observer. Geometrical and optical depth effects, in
<figure><img src="image_6.png" /><figcaption>**Figure 4.** PV diagrams along the major axis of Oph163131. The ¹²CO(2–1) isotopologue is shown on the left panel with the LSR velocity of the system in the left y-axis. The ¹³CO(2–1) is shown on the right panel using the velocity offset coordinates and assuming a systemic velocity of $v_{\text{sys}} = -4.4 \text{ km s}^{-1}$. The white dashed lines correspond to a Keplerian velocity profile for a star of 1.2 $M_{\odot}$. The white filled region is bracketed by upper and lower limits of 1.4 $M_{\odot}$ and 1.0 $M_{\odot}$, respectively.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** The left panel shows the JCMT ¹²CO *J* = 3–2 intensity map of the large-scale structure around Oph163131. In color, we displayed four different velocities channels where the emission of the cloud is the strongest (with $v_{\text{LSR}}$ ranging from 3.6 to 4.4 km s⁻¹). Four 14″ diameter (i.e., beam-sized) circular apertures are overlaid at locations within the cloud (light gray, gray, and black) and at the location of Oph163131 (pink). The right panel shows the JCMT spectra obtained at the locations mentioned above, and uses the same colors as the apertures in the left panel. The dashed blue vertical lines and the horizontal arrow depicts the velocity range over which the Oph163131 disk is detected above 3σ in the ¹²CO *J* = 2–1 ALMA data. There is no spatial or spectral overlap between the Oph163131 disk and the cloud.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** ALMA ¹²CO (in solid blue) and ¹³CO (in dashed red) spectra of the disk around Oph163131. The ¹³CO emission was scaled up by a factor of three for illustrative purposes.</figcaption></figure>

particular the thickness of the CO-emitting layer and absorption by the cooler midplane layer, result in the “bottom” layer being partially hidden whereas we can directly see a much larger surface area of the “top” layer. Interestingly, scattered light images of Oph163131 also show the SE layer being brighter than the NW layer, as can be seen from Figure 8, where the ~0.6 µm HST F606W scattered light image of Oph163131 is overlaid on the ¹²CO integrated intensity (moment zero) map (see also Paper I).

Although the ¹³CO emission is resolved along the minor axis of the disk (over-sampled 2.8 times), it is not possible to distinguish the two emitting layers seen in ¹²CO. The combination of a lower S/N and smaller vertical size of the disk in ¹³CO prevents us from distinguishing any additional structure in this isotopologue.

In Figure 9, we present a comparison between ¹²CO, ¹³CO, and 1.3 mm dust continuum brightness profiles along the minor and major axes of the disk. The left panel of Figure 9 shows
<figure><img src="image_9.png" /><figcaption>**Figure 7.** <sup>12</sup>CO (top) and <sup>13</sup>CO (bottom) ALMA channel maps of the Oph163131 disk. Each channel shows a square region of 4" in size centered on the predicted stellar position marked with a yellow star. The LSR velocity (in kilometers per second) is indicated in the top right of each channel panel. The synthesized beam is shown as a gray ellipse in the lower-leftmost channel panel.</figcaption></figure>

that the dust disk is smaller radially by a factor of almost two
than the <sup>12</sup>CO gas. Similar results where the gas component
extends farther in the radial direction than the dust component
have been reported in the literature and are interpreted as

inward radial drift of dust grains (Ansdell et al. 2018). We note
that the double peak component in the integrated gas profile is
the natural result of the high inclination of the disk, with some
self-absorption occurring along the line of sight to the central
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Left panel: overlay of the ¹²CO zero moment map shown in color and the HST F606W image described in Paper I in contours of 3σ, 10σ, 20σ, 40σ, and 60σ. The effective ALMA beam is indicated in the bottom left corner whereas two unrelated stars in the HST image are representative of the point-spread function in the optical. Right panel: overlay of the ¹²CO first moment map and the 1.3 mm continuum map described in Paper I. Contours for the continuum map are shown at the 5σ, 50σ, 100σ, 150σ, 200σ, 250σ, and 300σ level; the first moment map is only shown where the moment zero exceeds the 5σ level.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Brightness profiles along the major (left panel) and minor axis (right panel), integrated over the other dimension in both cases, of the dust millimeter continuum, ¹²CO, and ¹³CO data of Oph163131 (dashed blue, solid red, and dotted–dashed green, respectively). The horizontal ¹³CO profile has been convolved by a Gaussian with a width equal to the minor axis of the beam to improve signal-to-noise. All brightness profiles are displayed in arbitrary units, after normalization to their respective maxima. The horizontal bar in the top right represents the FWHM of the beam along its major axis.</figcaption></figure>

star. Villenave et al. (2020) measured a full radial extent of 2.50 ± 0.01 for the disk based on the continuum map (see also Paper I). The right panel of Figure 9 demonstrates that the dusty disk is vertically thinner than both gas components, and possibly much thinner, as the continuum is unresolved in this direction. From this panel, one can also see an asymmetric profile in the ¹²CO, caused by the brightness differences between the SE and NW layers, but a much more symmetric profile for the ¹³CO component, supporting the fact that no significant structure is detected in ¹³CO.

# 4. Tomographic Gas Temperature Retrieval

## 4.1. Method

Temperature estimates of protoplanetary disks are often obtained through detailed modeling of (sub)millimeter

observations or from long wavelength coverage spectral energy distribution fitting (e.g., Woitke et al. 2019). In both of these approaches, it is common to parameterize the temperature and the density structure of a disk and then vary many of these parameters until a good match between models and observations is reached. These many parameter fitting techniques are often under-constrained given the real complexity found in high-angular resolution disk observations.

Direct measurement of the temperature structure of disks is therefore not only important in our understanding of the radiative processes of disks, but it can, in principle, also be used to constrain current modeling work. Here we describe a model-independent technique developed to measure the temperature structure of protoplanetary disks. This technique, referred to as TRD, was first implemented by Dutrey et al. (2017) and was applied to the edge-on protoplanetary disk “the
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Demonstration of the TRD method assuming three different disk's inclination angles: $i = 90°$ (left), $i = 87°$ (middle), and $i = 84°$ (right). All three figures show the same PV diagram along the disk's major axis. The white dashed lines correspond to isoradial lines ranging from 40 to 300 au; these lines remain the same with inclination. The cyan solid lines represent is altitude lines (regions at the same vertical distance above the midplane). Assuming a perfectly edge-on disk (left panel), the flux from different radii emerge at a fixed altitude ($z = 0$ au in this case). As the disk is declined out of the edge-on position, the observed flux displays a radial and vertical dependency within the PV cut. At an inclination of $i = 87°$ (middle panel), the emission emerges from altitudes ranging between $1 \lesssim z \lesssim 10$ au, while at an inclination of $i = 84°$ (right panel), the emission arises from even higher altitudes, $2 \lesssim z \lesssim 20$ au.</figcaption></figure>

Flying Saucer.” A similar method was recently implemented by Teague et al. (2020) for Gomez’s Hamburger (2MASS J18091339–3210500).

The TRD technique is a geometrical approach designed to identify parcels of gas on a spatially and spectrally resolved protoplanetary disk and map these parcels of gas onto a two-dimensional radial and vertical structure. This technique assumes that the disk’s inclination angle is close to edge-on and that the disk follows a Keplerian velocity field. However, the latter point is not exactly correct, as gaseous disks are pressure supported systems, and therefore they do not rotate exactly at Keplerian velocities. It has been demonstrated that omitting the pressure term ($-\frac{1}{\rho}\frac{dP}{dr}$) only causes a small overestimation of the gas velocity at large disk radii. Pinte et al. (2018) analyzed the velocity pattern of the disk around IM Lupi (whose mass is similar to that of Oph163131) and found that the pressure gradient term only increases the rotational velocity of the disk by $0.1 \text{ km s}^{-1}$ at distances larger than 300 au. Since our disk does not extend beyond 300 au and our velocity resolution is already larger than this effect, we are safe to assume pure Keplerian rotation.

In order to convert the incident flux emitted from the gaseous component of the disk into brightness temperature, we assume that the gas is in local thermodynamic equilibrium (LTE). The LTE condition is valid in high-density regions, for example, toward the midplane of most protoplanetary disks, but not necessarily toward the disk’s photosphere. Besides, when the emission is optically thick and fills the beam, the brightness temperature becomes comparable to the gas’s excitation temperature. This is typically the case for $^{12}$CO, but it can also be true for the rarer isotopologue $^{13}$CO (see, for example, Pinte et al. 2018). Furthermore, the excitation temperature of low-$J$ CO lines, as is our case here, follows very closely the kinematic temperature of the gas (Pavlyuchenkov et al. 2007).

In the following paragraphs, we briefly describe the steps involved in the reconstruction of the temperature distribution of a protoplanetary disk. We refer the reader to Appendix A for the derivation of the equations we used to reconstruct the disk emission.

Our modeling technique depends on the stellar mass ($M_*$), the distance to the source ($d$), and the inclination angle of the disk relative to the observer ($i$). These values are input into Equations (A5) and (A8), which parameterize the radial position ($\rho$), and the altitude above the midplane ($z$) of small parcels of gas in a protoplanetary disk as a function of the observed line-of-sight velocity ($v_y$), the impact parameter ($x$), and the projected altitude above the midplane ($z'$). A PV diagram is then selected at a projected altitude $z'$, and the code calculates the radial $\rho$ and altitude $z$ coordinates, in addition to the intensity value $I$ of the disk at such locations. This process is repeated by selecting a new PV cut at a different projected elevation ($z' + \Delta z'$) to obtain a new set of intensity values for all of the radial positions $\rho$ and elevations ($z + \Delta z$). This iteration continues until all of the projected altitudes above and below the projected midplane are analyzed.

It is worth noting that in a perfectly edge-on case ($i = 90°$), the emission obtained from a PV diagram emerges from a fixed elevation in the disk, allowing us to construct a one-to-one mapping between projected and true midplane positions. In cases where $i \neq 90°$, however, the emission from a PV diagram obtained at a given projected elevation $z'$ comes from a range of elevations $\Delta z$ above and below such projected elevation. This is exemplified in Figure 10 where we show the same PV diagram, measured at $z' = 0$, but we use three different disk inclination angles. In the left panel, where $i = 90°$, the emission from different radial positions is expected to arise from a unique vertical elevation $z = 0$ au. When the inclination is assumed to be $i = 87°$ (middle panel), the emission from different radial positions emerges from a variety of elevations above the midplane $1 \lesssim z \lesssim 10$ au. Lastly, if the inclination is assumed to be $i = 84°$, then the range of elevations where the emission is assumed to be emitted is even larger than before $2 \lesssim z \lesssim 20$ au. Therefore, we emphasize the importance of considering small deviations out of edge-on inclinations when using the TRD technique.

After all of the intensity values $I$ and spatial coordinates $\rho$ and $z$ from the PV cuts are collected, we construct a two-dimensional array of radial ($R$) and vertical ($Z$) positions with steps of 10 au in both directions. This spacing size was chosen
<figure><img src="image_13.png" /><figcaption>**Figure 11.** TRD of the ¹²CO (top panel) and the ¹³CO (bottom panel) lines. The reconstructed temperature maps are shown in physical units of radius and altitude, both in astronomical unit. On the ¹²CO reconstructed temperature, we distinguish the near and the far emitting layers of the disk. The ¹³CO reconstructed temperature, on the other hand, does not show any significant structure. The gray ellipses shown on the top left part of the bottom panel illustrate how the Keplerian velocity shear affects the resolution as a function of radius.</figcaption></figure>

to oversample by factors of ~3 and ~4 the resolution of our observations in the radial and vertical directions, respectively. Intensity values are then assigned to each of these R-Z cells, and median values are calculated within each cell. To calculate the statistical uncertainty in our method, we additionally measure the 16th and 84th percentiles of the intensity values for each of the cells, and find a 2-4 mJy (1-2 K temperature) uncertainty in most parts of the disk. Before the median intensity is converted into brightness temperature, each of the R-Z cells with 10 or fewer individual intensity values or cells with median intensity below 9 mJy (or 4.5 K) were masked. This ensures that the statistical computations of the median, the standard deviation, and the percentiles are meaningful. At last, the average intensity values are transformed into brightness temperature using Planck's Equation (1), where $T_{bg}$ is the microwave background emission of 2.73 K.

$$T_b = \frac{h\nu}{k} \left[ \ln \left( \frac{2h\nu^3}{I_\nu c^2} + 1 \right) \right]^{-1} + T_{bg}. \quad (1)$$

In Appendix B, we show a variety of tests performed using synthetic disk models computed with MCFOST. These tests, which include the disk's inclination and spatial resolution variations, confirm that the temperature retrieval method works well in the interior regions of the disk where the density is large enough to meet our optical depth criteria.

## 4.2. The Temperature Structure of Oph163131

We now apply the TRD method to the ¹²CO J = 2-1 and ¹³CO J = 2-1 data cubes of Oph163131¹¹. We fixed the stellar mass to 1.2 M☉, the disk inclination angle to i = 87°, and the

distance to the source to d = 147 pc. The top panel of Figure 11 shows the result of the reconstructed temperature distribution for the ¹²CO line. The reconstruction technique was applied to both sides of the disk (positive and negative velocities with respect to vsys), both of which show a maximum radial extent of Rmax = 350 au and altitude above the midplane of Zmax = 80 au. From this plot, we can clearly recognize the “top” and “bottom” layers of the disk, which are separated by a colder midplane that extends from about 100 to 200 au in radius. The peak brightness temperature on both sides of the disk reaches about the same temperature of 35 ± 2 K at a distance of 135 ± 20 au from the center of the disk. The displacement of the peak brightness temperature outside of the inner part of the disk is mostly caused by the beam convolution effects, which are described in Appendix B. However, toward the innermost regions of the disk R ≲ 1 au (which cannot be accessed at our spatial resolution), dust absorption could also play a role in decreasing the apparent brightness temperature of the gas and therefore shift the peak temperature position.

In the disk, emission detected at radial distances smaller than ~130 au is more heavily affected by beam dilution due to the smaller scale height at the inner parts of a flared disk. As demonstrated in Appendix B.3, beam convolution results in underestimating the maximum temperature of the disk and overestimating the gas temperature close to the midplane. The temperature gradient from the ¹²CO reconstructed temperature is therefore only a lower limit to the true temperature gradient in Oph163131.

The bottom panel of Figure 11 shows the reconstructed temperature distribution for the ¹³CO isotopologue. This reconstructed distribution is radially and vertically more compact than the ¹²CO temperature distribution, and we only detect emission at vertical distances smaller than 45 au and radial distances between 60 au and 250 au. A similar reduction in the relative height of the CO isotopologue emission layers

<footnote n="11">¹¹ We also applied this method to the non-continuum subtracted data cubes, but no significant differences were found.</footnote>
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Panels (a) and (b) show the left-to-right averaged temperature distribution from Figure 11. The colored lines correspond to a visual trace over the near, far, and cold midplane layers for the ¹²CO, and just center for the ¹³CO. The wider transparent region next to the white traces corresponds to the three pixel (five pixel) vertical average used to create panel (c) for ¹²CO and ¹³CO, respectively. Panel (c) shows the traced temperatures as a function of radius for the cold midplane, bottom, and top layers of the ¹²CO map, and the center temperature for the ¹³CO map. The shaded regions around the solid lines correspond to the temperature uncertainties for each layer. The horizontal bars in panel (c), drawn at distances of 100 and 250 au, depict the changing radial resolution of this method (see Section 4.2).</figcaption></figure>

was identified by Teague et al. (2020) for Gomez’s Hamburger. A combination of lower S/N, lower optical depth, and beam smearing prevents us from identifying clear structures in this isotopologue. The detected ¹³CO emission is unlikely to arise from the midplane due to CO freeze-out; instead, the emission we detect is likely emitted by thin layers just above and below the freeze-out zone, which are completely unresolved by our limited ~0.2′′2 spatial resolution.

In the lower left panel of Figure 11, two gray ellipses represent the spatial resolution variations at different radial distances. These radial resolution differences are due to the Keplerian shear velocity of the disk, calculated as $dr = 2r dv/v$ ($r$), where $v(r)$ is the Keplerian velocity of the disk, and $dv$ is the local line width. For distances smaller than 160 au, the beam size dominates the radial resolution limit of the map; however, at distances larger than 160 au, the radial resolution is limited by the Keplerian velocity shear. The two ellipses shown at 160 au and 250 au represent then radial resolution variations from 30 au and 60 au. We note that the local line width in our data set is dominated by our velocity resolution of 0.25 km s⁻¹, and not by the combined thermal and turbulent line broadening, as was the case for the analysis performed on the Flying Saucer (Dutrey et al. 2017).

### 4.3. Radial and Vertical Temperature Profiles

To better understand the radial temperature structure of the disk, we show the temperature variation as a function of radius for ¹²CO and ¹³CO in Figure 12. First, we averaged both sides (left and right in Figure 11) of the temperatures maps, and then visually identified three distinct layers for the ¹²CO (panel (a)) and a single layer going through the center of the ¹³CO (panel (b)). A three pixel average in the vertical direction (24 au) for the ¹²CO, and a five pixel average (40 au) for the ¹³CO were used to compute the temperature profiles shown in panel (c).

Panel (c) of Figure 12, shows that the temperature of the top layer (in blue) increases from 27 K at 60 au to a peak temperature of 35 K at 150 au, then it declines to a temperature

of roughly 31 K at a distance of 210 au and remains constant toward the outer part of the disk. The initial increase in temperature as a function of radius is not physical, but mostly caused by beam dilution (as discussed in Section 4.2). On the other hand, the subsequent temperature decline is meaningful and allows us to understand how rapidly the temperature changes as a function of radius. The bottom (purple) and cold (green) layers share very similar temperature profiles. The bottom layer reaches a minimum temperature of 23 K at 150 au, and it increases toward the outer edge of the disk to a steady temperature of 25 K beyond 210 au. The cold midplane, on the other hand, reaches a minimum temperature of 21 K at 150 au and then increases to a temperature of 28 K beyond 210 au. The temperature increase of the cold midplane region toward the outer parts of the disk is an interesting feature and is further discussed in Section 5.3.

The ¹³CO temperature profile (in red) reaches a maximum temperature of 19 K at distances of 150–200 au, and due to beam dilution, it decreases toward shorter and longer radial distances. The lower temperature found for ¹³CO compared to ¹²CO likely results from the combination of (1) a lower optical depth (similar to, or lower than, unity), (2) a smaller vertical extent, whereby the two emitting layers of ¹³CO (above and below the presumed frozen-out midplane) are only partially resolved by our observations, and (3) a lower temperature in these more embedded layers.

Although it seems clear from Figure 11 that in ¹²CO there are a top and a bottom emitting layers separated by a colder midplane, the three pixel averaged profiles of Figure 12 suggest that the midplane is not significantly cooler than the bottom layer. We further investigate this ambiguity by looking at vertical temperature profiles of the ¹²CO map. As shown in Figure 13, we chose two regions to perform the vertical cuts, one region at a radial distance of 135–175 au, and a second region at radial a distance of 245–285 au. The first region was chosen to go through the cold midplane region in the disk (marked as a light blue region), while the second one was chosen to go through the outer disk (marked as a light green
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Top panel: left-to-right average of the ¹²CO temperature map. The two shaded regions at distances of 135–175 (blue) and 245–285 (green) mark the distances where radial averages are performed. Bottom panel: temperature as a function of the disk's altitude taken at 135–175 (blue) and 245–285 (green) as shown in the top panel. The shaded regions around the solid lines correspond to the temperature uncertainties calculated from the maps. The horizontal bar shown in the bottom panel corresponds to the spatial resolution of 36 au in this direction.</figcaption></figure>

region). In the bottom panel of Figure 13, we show the temperature of the disk as a function of vertical distance above and below the midplane in the two regions marked above. The blue line clearly depicts the top and bottom emitting regions and the colder midplane zone mentioned before, with the midplane region being $\delta T = 3 \pm 2$ K cooler than the bottom region and $13 \pm 2$ K cooler than the top region. This is qualitatively consistent with TRD maps of the Flying Saucer and Gomez's Hamburger (Dutrey et al. 2017; Teague et al. 2020). The temperature structure changes significantly when the cut is performed toward the outer region. The green curve does not show a valley structure, which shows that the midplane temperature has increased at radial distances larger than R ~ 200 au, but there is still a marginal level of asymmetry showing that the top side of the disk still appears to be slightly warmer than the bottom side at larger distances from the center.

The temperature asymmetry between the top and bottom sides of the disk is an optical depth effect caused by a slight deviation out of the edge-on inclination of Oph163131 coupled with the significant line optical depth and vertical thickness of the warm molecular layer. At the top of the disk ($z > 0$), the optical depth $\tau = 1$ is reached toward the disk's surface, therefore providing the temperature of the disk's high-elevation layers, which are warmer as they are directly exposed to the star. This effect is further amplified by the fact that, even at

small deviations out of the edge-on position, the farside of the disk, on the top side, becomes visible to the observer (even if the disk is not resolved), causing a larger area of the disk to be exposed to the observer and consequently artificially increasing the measured temperature (see Appendix B.2). On the bottom side ($z < 0$ au) of the disk, however, the $\tau = 1$ surface is reached closer to the midplane from the observer's perspective, providing a temperature measurement of the disk's lower-elevation regions. Additionally, it is expected that at the locations where the dust continuum is optically thick, the gas emission could be absorbed, further lowering the brightness temperature. Given the geometry of the disk, this effect only plays a role in regions below the midplane ($z < 0$). This is why even after the TRD technique transforms the disk into R–Z coordinates, there is still a temperature asymmetry between the two sides of the disk.

The temperature asymmetry, combined with the angular resolution of the observations and the intrinsic angular scales of the Oph163131 disk, also explains the off-centered cold disk midplane (by approximately −20 au) shown in Figure 13. The more exposed CO upper emitting layer appears broader and hotter than the less-exposed CO bottom emitting layer, which causes the whole disk midplane to be shifted toward the bottom side of the disk (see also Figure 17 of the Appendix). All of the above is then amplified by the flux smearing induced by the spatial resolution of the observations. Finally, the close match between the temperature of the bottom ¹²CO emission and of the (symmetric) ¹³CO is consistent with observations of IM Lupi (Pinte et al. 2018) and confirms a freeze-out temperature of ~20 K.

# 5. Discussion

## 5.1. Oph163131 Compared to Other Sources

We now put Oph163131 in context by comparing its properties against other sources with direct disk temperature measurements. We focus on four sources with spectrally and spatially resolved ¹²CO observations, where a model-independent temperature measurement was performed: IM Lupi (Pinte et al. 2018), the Flying Saucer (Dutrey et al. 2017), Gomez's Hamburger (Teague et al. 2020), and HD 163296 (Isella et al. 2018; Dullemond et al. 2020). Two of the sources (Gomez's Hamburger and the Flying Saucer) harbor edge-on disks and have been analyzed using the TRD method in a very similar fashion as we have done here for Oph163131, therefore a fully reconstructed temperature map is available for comparison. The other two sources (IM Lupi and HD 163296) have disks with intermediate inclination angles ($i \sim 45^\circ$) and report direct temperature measurements on specific parts of the disk (i.e., at the surface level or toward the midplane).

The central stars of these sources span a relatively large range of masses; Gomez's Hamburger and HD 163296 are regarded as Herbig stars with stellar masses of $2.5 \pm 0.5$ M$_\odot$ (Teague et al. 2020) and $2.04_{-0.14}^{+0.25}$ M$_\odot$ (Andrews et al. 2018), respectively. IM Lupi and the Flying Saucer, on the other hand, have stellar masses of $1.0 \pm 0.1$ M$_\odot$ and $0.58 \pm 0.01$ M$_\odot$, respectively, both consistent with being T Tauri stars, which makes them better candidates to be compared against the disk of Oph163131. The much higher luminosity of Herbig stars should result in substantially higher temperatures throughout the disk (e.g., Kamp et al. 2011).
The maximum temperatures measured from ¹²CO observations in the disks of the Herbig stars Gomez's Hamburger and HD 163296 reach values of ~60 K and ~90 K, respectively, which are significantly hotter than the temperatures measured in the T Tauri stars' disks of ~30 K for IM Lupi and Oph163131, and of ~20 K in the disk of the less-massive star, the Flying Saucer. Minimum temperatures of the disks are harder to compare, as they strongly depend on where the temperatures are measured, how well resolved the disks are, and the sensitivity of the observations.

Pinte et al. (2018) produced radial temperature profiles computed from the ¹²CO and the ¹³CO lines for IM Lupi (their Figure 4, right panel), which can be easily compared to our Figure 12. In ¹²CO, the maximum temperatures of the disks are 35 and 32 K, which are reached at distances of ~150 au and ~190 au for Oph163131 and IM Lupi, respectively. The temperature declines toward larger radii, however, and is markedly shallower in Oph163131 than in IM Lupi: ~2.5 K per 100 au and ~6 K per 100 au, respectively. In ¹³CO, the radial profiles of both sources reach a plateau with peak values of 19 and 21 K for Oph163131 and IM Lupi, respectively, with similar decreases toward the outermost regions of the disks. We note, however, that contrary to the IM Lupi study, our ¹³CO observations do not fully resolve the layers on either side of the midplane, and thus, the ¹³CO temperature profile reported in our study is only a vertical average, making the comparison less straightforward.

For most of the sources mentioned above, we have estimated the ratio between the disk's full vertical extent ($Z_{max}$) and the total disk's radial size ($R_{max}$) in ¹²CO. The thicker disk, using this metric, is the Flying Saucer with a ratio of $Z_{max}/R_{max} \sim 0.5$. IM Lupi and the Gomez's Hamburger come next with a similar vertical-to-radial ratio of $Z_{max}/R_{max} \sim 0.33$, and finally the thinnest disk is Oph163131 with a ratio of $Z_{max}/R_{max} \sim 0.23$. In this comparison, we have excluded HD 163296, as the vertical size measurement performed on this disk is strongly model dependent and therefore not a fair comparison against the other sources (Isella et al. 2018).

Despite contrasting disks sizes of sources obtained at different spatial resolutions, these comparisons can provide some insights into how the stellar masses, evolution, and temperature of the disks relate to the aspect ratio of the sources. For T Tauri stars, it appears that more massive stars correlate with thinner disks; although with only three sources in the sample, this result lacks statistical significance. An alternative view is that each disk could be in a different evolutionary state where vertical settling may have progressed by different amounts. In the case of Oph163131, the absence of strong accretion (see Section 3.1) means that there is much less UV radiation available to heat up the gas, which ultimately causes the flat aspect of this disk. Alcalá et al. (2017) measured an accretion luminosity of $0.08 L_\odot$ for IM Lupi by modeling the excess emission from the UV to the NIR of X-shooter data. The excess of accretion radiation might explain the larger $Z_{max}/R_{max}$ value for IM Lupi when compared to Oph163131. As for the Flying Saucer, we have not found accretion measurements in the literature. HST images show no signs of any jet or molecular outflows in the vicinity of this source, yet this does not rule out stellar accretion.

The accretion status of other edge-on protoplanetary disks that have a low height-to-radius ratio and display flat parallel

nebulae in scattered light images (as seen in Oph163131) is discussed next. HK Tau B, a well-studied edge-on disk (Koresko 1998; Stapelfeldt et al. 1998), was classified as a classical T Tauri star (CTTS) by Monin et al. (1998), albeit with an optical spectrum that displays a modest $H_\alpha$ emission and no other emission line. The optical spectrum of LkHα 263C, a similarly flat disk (Chauvin et al. 2002; Jayawardhana et al. 2002), reveals that it has a strong $H_\alpha$ emission along with prominent forbidden emission lines. Finally, ESO Hα 574, another flat-shaped disk imaged with HST (Stapelfeldt et al. 2014), has a prominent optical jet, which suggests ongoing stellar accretion. In summary, flat-shaped protoplanetary disks are found around actively accreting and non-substantially accreting stars. We encourage a study that would quantify the accretion radiation and aspect ratio of disks to test whether accretion-induced UV radiation could significantly affect the vertical thickness of protoplanetary disks.

## 5.2. The Inner Disk: Inside 200 au

The reconstructed ¹²CO temperature map of Oph163131 shows two distinct regions with a demarcation at a radius of ~200 au. Inside of this radius, the CO emission arises from two elevated layers separated by a midplane, whereas the outer region consists of a much more vertically uniform emission. Interestingly, this radial break coincides with the outer radius at which both the submillimeter continuum emission and the scattered light are detected (see Figure 8 and also Villenave et al. 2020). We, therefore, address the “inner” and “outer” regions separately. Considering first the inner disk, we note that the observed behavior of the gas emission is similar to that observed in other disks (e.g., IM Lupi and the Flying Saucer discussed above). We now wish to go beyond qualitative comparisons and evaluate whether detailed radiative transfer models of the dust components of Oph163131 can reproduce the temperature structure we obtain from the gas.

In Paper I, ALMA 1.3 mm dust continuum observations and 0.8 µm HST scattered light images of Oph163131 were modeled. To reproduce the observations, Paper I uses the MCFOST radiative transfer code and a Markov Chain Monte Carlo approach to fit the millimeter and submicron data of Oph163131 independently. The individual fits, which assumed a tapered-edge surface density profile, provided a set of density and geometrical parameters of the disk such as inclination angle, scale height, surface density, and flaring exponent (see their Equation (1) and Table 4). In this section, we compare the models from Paper I to our ¹²CO reconstructed temperature map, limiting the analysis to the inner 200 au.

To perform a meaningful comparison, we took the following approach. We used Paper I dust disk density and geometrical parameters, a gas-to-dust ratio of 100, and a CO/H₂ abundance of $10^{-4}$, and we imposed CO freeze-out (CO/H₂ = 0) in regions where the temperature of the disk falls below 20 K. We then computed the synthetic ¹²CO J = 2–1 observations of such a disk with MCFOST. From the resulting data cubes, we analyzed these models in the exact same manner as we did with gas observations in Section 4, i.e., using the same parameters for Oph163131, such as stellar mass, spatial and spectral resolution, and distance to the star. The two models we consider are: (a) a model that fits the HST 0.814 µm image (Paper I's model A), which has a disk inclination angle of i = 84°5, a scale height of $H_0 = 7.2$ au (at 100 au), and a critical
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Comparison between the two temperature maps produced using the MCFOST synthetic models with parameters adopted from Paper I. The top panel corresponds to the temperature map obtained using the 0.8 µm HST best-fit parameters (model A), and the bottom panel from ALMA 1.3 mm continuum parameters (model B).</figcaption></figure>

radius of $R_C = 89$ au, where $R_C$ has the usual meaning of the distance at which the disk surface density changes from a power law to an exponential taper; and (b) a model that reproduces the 1.3 millimeter continuum ALMA map (Paper I's model B), with an inclination angle of $i = 87^\circ 2$, a scale height of $H_0 = 9.5$ au (at 100 au), and a critical radius of $R_C = 74.5$ au. The model retrieved temperature maps are shown in Figure 14.

Both models produce a disk with clearly defined top and bottom CO emission layers separated by cold and frozen midplane. The maximum temperature of the gas, however, is almost twice as hot as the data, and the elevation of the CO layers in the 150–200 au radial range is also significantly higher (50 au) than the observations (~35 au). This means that the parameters that reproduce submicron- (model A) and millimeter-sized (model B) dust grains do not necessarily match the gas observations when adopting standard gas-to-dust and CO-to-H₂ abundance ratios.

One possible path to reconcile the dust and gas observations is to reduce the overall CO gas mass in the disk. This would decrease the altitude at which the CO becomes optically thick and also the temperature of the emitting layer, as the gas is found at deeper layers in the disk. In practice, this means either adopting a lower gas-to-dust ratio (<100) or a lower CO-to-H₂ abundance ratio (<10⁻⁴) in the disk. Ansdell et al. (2016) performed a large study of young stars in Lupus and found that most disks have gas-to-dust ratios of less than 100 for an assumed standard CO-to-H₂ abundance (10⁻⁴). Long et al. (2017) analyzed a large number of disks bearing sources in Chameleon I and found median gas-to-dust ratios of ~4–15, depending on the assumptions made. Recent observations of individual disks have shown that CO could be depleted by several orders of magnitude with respect to H₂ (Favre et al. 2013; Kama et al. 2016; Trapman et al. 2017; Zhang et al. 2017, 2020). These anomalous abundances have been explored using physical and chemical models that deplete CO in disks (e.g., Bosman et al. 2018). In particular, Krijt et al. (2020) used

models that combine chemical processing of CO with ice sequestration in the disk midplane to obtain strong CO depletion in the outer disk warm molecular layers.

Although the possibility of a lower gas-to-dust ratio and CO-to-H₂ abundance in disks is extensively suggested in the literature, such ad hoc adjustment in our models of Oph163131 does not necessarily provide a satisfactory match to the gas observations. One example of this is that using a lower gas-to-dust ratio would immediately affect the dust settling mechanism, and thus, other back-reactions could occur if lower CO-to-H₂ abundances were adopted. As discussed in the next subsection, only detailed modeling of the CO observations would provide the appropriate physical parameters of the gas component of the disk around Oph163131.

### 5.3. The Outer Disk: Outside 200 au

In Section 4.3, we pointed out that the outer disk ($R > 200$ au) of Oph163131 has an almost constant temperature of ~30 K at all elevations and radii. This is readily seen from panel (c) of Figure 12 where all three layers converge to a similar temperature, and in Figure 13 where the outer vertical temperature profile is almost flat-topped. This result is somewhat surprising in multiple ways. First of all, we did not observe this behavior in the tests performed on the synthetic disk models. In these tests, we varied the beam size and the disk inclination and yet always obtained a top and a bottom layer separated by a colder midplane (Appendix B), indicating that this is not a feature of the TRD technique. Therefore there appears to be an additional process that leads to a much more vertically uniform temperature than expected. Second, a sharp jump in temperature is inconsistent with typical trends as a function of distance from the central star. This indicates that there is an important transition in the disk properties at this location. Finally, the 30 K that is observed in the outer disk is much higher than the expected blackbody equilibrium temperature, which is about 20 K and 15 K at 200 and 300 au, respectively, assuming low-albedo dust grains. Since this
temperature has a very shallow dependency on stellar luminosity, it is highly unlikely that this could be due to the star being much more luminous than we assumed in our models. This points to the gas being overheated compared to our expectations. Here we speculate on a possible explanation for these behaviors.

The radial position where the midplane temperature of Oph163131 increases from 21 K to 30 K coincides with the position where the millimeter continuum emission is no longer detected (R ~ 200 au; see Figure 8 and Paper I). A similar behavior is seen in the Flying Saucer (Dutrey et al. 2017), in which a sharp (albeit smaller) temperature jump is observed at the outer radius of the millimeter continuum map. However, one significant difference between the two disks is that submicron dust is traced in scattered light out to much larger radii in the Flying Saucer, whereas in Oph163131, the radial extent of the millimeter continuum and scattered light images are very similar to one another (see Paper I). We speculate that the dust grain difference relates to the observed midplane temperature profile outside of the millimeter continuum: it is uniform for Oph163131 but gradually declining for the Flying Saucer.

Dutrey et al. (2017) proposed that this sharp temperature increase is related to a higher influx of UV radiation penetrating the disk from the outside, leading to enhanced thermal or direct CO desorption. This is supported by our observations of Oph163131. Broadly speaking, the depletion of dust in the outer disk leads to a lower optical depth and thus higher penetration of short-wavelength radiation. This process is even more dramatic in the case of Oph163131 where even small dust grains, which have high UV opacity, are absent or rare in the outer disk. In this case, we expect interstellar photons to penetrate deeper into the disk and to deposit more heat, leading to an even more uniform temperature in that region.

The much reduced dust density could also play a role in the heating of the outer parts as the gas could be more decoupled from the dust, allowing it to remain at a higher temperature. Large dust grains are known to be the most efficient emitters of far-IR and submillimeter radiation, and without them, the outer disk cannot radiatively cool as efficiently as in the inner parts. In this regard, we note that Oph163131 is positioned between the L1689 and L1688 Ophiuchus molecular clouds, in a region of substantially lower dust and gas column density than the Flying Saucer (Lombardi et al. 2008; Ladjelate et al. 2020). The remarkably uniform and warm temperature of the outer region of Oph163113 could thus be the consequence of enhanced interstellar UV radiation coupled with much lower dust opacity.

To go beyond pure speculation, we computed MCFOST synthetic disk models to test this hypothesis. We used the stellar parameter derived for Oph163131 and used disk geometrical and density parameters that visually matched the full extent of the disk around Oph163131. None of these models are meant to reproduce the observations of the gas component of Oph163131 but to demonstrate the effects of a higher external UV flux on the disks. We first constructed a “standard” model that consists of an inner region (R< 170 au) with a dust mass of $10^{-3} M_\odot$ and an outer region with a dust mass of $10^{-6} M_\odot$ that extends from 170 au to 300 au. Both regions have a gas-to-dust ratio of 100, but the outer region contains only small grains (the maximum grain size was set to 10 µm). This last requirement comes from the lack of submillimeter dust emission detected for Oph163131 at distances >200 au (see Paper I and Villenave et al. 2020).

Finally, to perform a meaningful comparison to our results, we applied the TRD technique to this synthetic data, just as we did in Section 5.2.

Figure 15 shows the TRD maps for above mentioned models. The first model (panel (a)) corresponds to the “standard model,” in which we did not include photodesorption or photodissociation of the CO molecule, as a purely illustrative exercise. This model shows two well separated top and bottom layers extending from the inner part of the disk and all the way to the outer disk. In this model, the temperature of the midplane in the outer parts is of the order of 17 K, while the top layer reaches a maximum temperature of 26 K, indicating a strong vertical gradient that is not seen in the data. In the second model (panel (b)), we show the same “standard” model, but now we have activated photodesorption and photodissociation in the model, and external illumination is implemented with a UV flux set to the UV field in the solar neighborhood value (Draine & Bertoldi 1996). We see from this model that the interstellar photons have eroded the upper and bottom layers of the disk, creating a vertically flatter outer disk. Although the temperature of the midplane in the outer regions is still of the order of 17 K, we notice that the temperature has an almost uniform value at all vertical elevations. The third model (panel (c)) is a variant of the same model where the UV irradiation field is 100 times stronger than the value in the solar neighborhood. The increased irradiation field does not alter the shape of the disk any further, but it homogeneously increases the temperature of the outer disk to 22 K. In this model, we see the distinctive feature of an outer disk’s midplane that is hotter than the inner disk’s midplane, i.e., a radial temperature increase in the midplane’s temperature at radii larger than ~100 au. Finally, in a fourth model (panel (d)), we used the same parameters as in the previous model, i.e., a UV flux 100 times the solar neighborhood value, but we also increased the mass of the outer disk by a factor of 10, i.e., to $10^{-5} M_\odot$ in dust ($10^{-3} M_\odot$ in gas). As can be seen in Figure 15, the increased outer disk mass leads to self-shielding of the CO gas in the outer disk, producing the formation of distinct upper and lower layers, a colder midplane, and an overall flared appearance of the outer disk. The temperature of the top layer increased to 28 K, but the midplane temperature decreased to 20 K. Although it is not shown in the figure, we tested models with an even larger UV irradiation field. By the time the outer disk reached a temperature of 30 K, the radiation in the models was so high (1000 times the solar neighborhood value) that it causes the whole disk to increase in temperature.

In conclusion, we showed that photodissociation and photodesorption are processes that can explain the vertically flat outer part of a disk, due to the erosion of the low-density top and bottom layers, as long as the overall gas column density (equivalently, mass) in the outer layer remains quite low, requiring a sharp decrease in surface density at 170–200 au. The stripping of these layers produces in turn a vertically isothermal outer disk. In addition, when the UV irradiation field is increased (compared to solar neighborhood values), it causes a rise in the temperature of the outer parts, allowing the possibility of having an outer disk midplane that is hotter than the inner disk midplane. Finally, we point out that to obtain optimal disk parameters for this source, it would be necessary to perform a detailed radiative transfer modeling of the CO, varying density and geometrical parameters, and
<figure><img src="image_17.png" /><figcaption>**Figure 15.** MCFOST models with external UV irradiation might explain the flat-topped outer regions of Oph163131 as well as the increased temperature toward the outer parts of the disk. Panel (a)) shows the generic model without external irradiation. Panel (b)) shows the generic model, but we turn on photodissociation and photodesorption with a UV irradiation environment equal to solar. Panel (c)) shows a 100 times increased UV irradiation environment. Panel (d)) shows the generic model with 100 times the solar UV irradiation (same as (c)), but the mass of the outer disk is increased by a factor of 10.</figcaption></figure>

including external UV irradiation along with photodesorption and photodissociation processes.

## 6. Conclusion

In this study, we present the first ALMA CO observations of the disk around Oph163131, as well as the first optical and NIR spectroscopic observations of the central star. A summary of our results is as follows:

1. We obtained optical and NIR low-resolution spectra ($R \sim 1000-2000$) of Oph163131 with the Goodman and SpeX spectrographs, respectively. By comparing the strengths of several absorption lines in both wavelength regions against appropriate stellar templates, we derived a spectral type of K4±1 IV for the star.

2. From the optical spectrum of Oph163131, we measured a $0.6 \pm 0.1$ Å equivalent width of the $H\alpha$ line. This measurement in addition to a null detection of the $Br\gamma$ or $Pa\beta$ emission lines indicates that Oph163131 is currently not accreting at an appreciable level. A detection of the He I line at 1.08 µm, on the other hand, suggests that there might be either outflow activity or a very low accretion.

3. We used the velocity pattern imprinted in the disk by the star to infer a dynamical mass for Oph163131. Assuming a distance of 147 pc and an inclination angle of $i = 87^\circ$, we derived a stellar mass of $1.2 M_\odot$ for the central source with conservative upper and lower mass limits of $1.4 M_\odot$ and $1.0 M_\odot$, respectively.

4. We have generalized the TRD method developed by Dutrey et al. (2017) to include non-perfectly edge-on disks. This generalization increases the number of disks that can be analyzed using this technique. In Appendix B, we have characterized some of the limitations of this method that can produce, for example, an apparent temperature decrease in the low-density regions of the disk and a temperature increase of the top side of a disk (closest to the observer) in the presence of small inclinations.

5. Using the TRD method, we obtained a model-independent measurement of the temperature structure of the disk surrounding Oph163131. Toward the inner $R > 200$ au regions of the disk, we detect two emission layers separated by a dark lane, which we interpret as the top and bottom parts of a highly inclined disk separated by a cold and probably frozen-out midplane.

6. The colder midplane temperature of 21 K in the disk around Oph163131 increases to ~30 K for distances $R > 200$ au. This temperature increment coincides with a decrease in the micron- and millimeter-sized dust detection in the disk. MCFOST models with increased UV irradiation qualitatively match the temperature increase as well as the vertically thinner structure of the outer regions. We propose that external UV irradiation could heat up the gas component in the outer disk, as the interstellar photons can more easily penetrate regions where dust particles are not present.

Overall, the disk surrounding Oph163131 is markedly more compact along the vertical direction than other well-studied
disks, in both its dust and gaseous components. Combined with strong evidence for dust settling, low-to-absent accretion on the central star, and a low-density, overheated outer region, this suggests that the system is in a relatively advanced state of evolution that warrants further analysis. In particular, detailed line radiative transfer modeling of the CO maps presented here and new observations that are focused on different molecular tracers could provide a thorough picture of the vertical temperature profile, from the midplane to the upper surface of the disk.

We thank the referee for an insightful review, which helped to improve this paper. C.F. and M.S.C. acknowledge support from the National Aeronautics and Space Administration (NASA) Infrared Telescope Facility, which is operated by the University of Hawaii under contract 80HQTR19D0030 with NASA. C.F. also acknowledges support from the Graduate Student Organization at the University of Hawaii at Manoa. G.D. acknowledges support from NASA grants NNX15AC89G and NNX15AD95G/NExSS as well as 80NSSC18K0442. K.R.S. and D.L.P. acknowledge support from HST GO grant 12514 from the Space Telescope Science Institute. F.Me., M.V., and G.vd.P. acknowledge funding from ANR of France under contract No. ANR-16-CE31-0013. C.F. is grateful to Nienke van der Marel for an early mentoring that later developed into this publication.

*Facilities*: ALMA, JCMT, SOAR, IRTF.

*Software*: Astropy (Astropy Collaboration et al. 2018), CASA (McMullin et al. 2007), MCFOST (Pinte et al. 2006, 2009), Spextools (Cushing et al. 2004), xtellcor (Vacca et al. 2003).

# Appendix A

## Equations to Perform the Temperature Reconstruction of a Non-perfectly Edge-on Disk

In this section, we provide the equations that relate a protoplanetary disk's velocity and geometrical coordinates as seen from an observer's point of view to the velocity and geometrical coordinates in the disk's coordinate system.

We define the disk's spatial coordinates as {$x, y, z$} and its velocity coordinates as {$v_x, v_y, v_z$}. An observer at a finite distance $d$ and whose line of sight is inclined at an angle $\alpha$ with respect to the disk's midplane defines its own coordinate system {$x', y', z'$} and velocity coordinates {$v_x', v_y', v_z'$}. Both coordinate systems are related by the following matrix transformation:

$$ \begin{pmatrix} x' \\ y' \\ z' \end{pmatrix} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & \cos \alpha & -\sin \alpha \\ 0 & \sin \alpha & \cos \alpha \end{pmatrix} \cdot \begin{pmatrix} x \\ y \\ z \end{pmatrix}. \quad (A1) $$

We further assume that the disk is rotating at Keplerian speed $v = \sqrt{GM/r}$ and that the velocity component in the y-direction in the disk's reference frame can be described as:

$$ v_y = x \sqrt{\frac{GM}{(\rho^2 + z^2)^{3/2}}}; \quad \rho^2 = x^2 + y^2 \quad (A2) $$

in the second part of Equation (A2), $\rho$ represents the disk's radial distance in cylindrical coordinates as viewed from the disk's reference frame. We now combine Equations (A1) and (A2) to derive the line-of-sight velocity of the disk from the

observer's perspective:

$$ \begin{align} v'_y &= v_y \cos \alpha - v_z \sin \alpha; \quad v_z &= 0 \nonumber \\ v'_y &= x \sqrt{\frac{GM}{(\rho^2 + z^2)^{3/2}}} \cos \alpha. \tag{A3} \end{align} $$

In $v_z = 0$, we assumed that the disk is not moving in the z-direction. We now would like to express the radial distance $\rho$ in the disk's rest frame as a function of the observer's coordinates $v'_y, x',$ and $z'$. From Equations (A3) and (A1), we have

$$ \begin{align} \rho^2 &= \left( \frac{x'^2 GM \cos^2 \alpha}{v_y'^2} \right)^{2/3} - z^2 \nonumber \\ \rho^2 &= \left( \frac{x'^2 GM \cos^2 \alpha}{v_y'^2} \right)^{2/3} - (z' \cos \alpha - y' \sin \alpha)^2, \tag{A4} \end{align} $$

which, after some algebra, we can write as:

$$ \rho^2 = x'^2 + (z' \sin \alpha \pm \cos \alpha \sqrt{gm - (x'^2 + z'^2)})^2 \quad (A5) $$

with

$$ gm = \left( \frac{x'^2 GM \cos^2 \alpha}{v_y'^2} \right)^{2/3}. \quad (A6) $$

We take the real and positive solution of $\rho$, appropriate for the radius in cylindrical coordinate.

We additionally solve for the vertical geometrical component $z$ from the disk's coordinate system as a function of the observer's coordinates. Using Equation (A1), we have

$$ \begin{align} z &= z' \cos \alpha - y' \sin \alpha \nonumber \\ y' &= \frac{\sqrt{\rho^2 - x'^2} - z' \sin \alpha}{\cos \alpha}, \tag{A7} \end{align} $$

which we combine to obtain

$$ z = z' \cos \alpha - (\sqrt{\rho^2 - x'^2} - z' \sin \alpha) \tan \alpha. \quad (A8) $$

Therefore, for a given point in the PV diagram ($x', v'_y$ space) and for a given $\alpha$ angle and $z'$ vertical position, we can use Equations (A5) and (A8) to find the radial and vertical position of an emitting parcel of gas in the disk.

# Appendix B

## TRD Model Testing

In order to understand the limitations of the TRD method, we performed several tests on synthetic disk models. MCFOST is a line and continuum radiative transfer code that can produce atomic and molecular line maps of complex structures such as protoplanetary disks. To create a synthetic disk model, MCFOST first computes the temperature structure of a model assuming one or several energy sources. Under the assumption that dust and gas are thermally coupled (LTE), MCFOST then calculates the intensity of a given atomic or molecular species, creating a line map that can be compared to observations. In our models, we assumed that the dominant source of energy is a single star with stellar parameters consistent with our results in Sections 3.1 and 3.2. We assumed a tapered-edge disk
<table><caption>**Table 1**
Disk and Stellar Parameters of the Test Model</caption><thead><tr><td>Model parameter</td><td>Value</td></tr></thead><tbody><tr><td>a<sub>max</sub> (μm)</td><td>1000</td></tr><tr><td>M<sub>dust</sub> (M<sub>☉</sub>)</td><td>5.5 10<sup>-4</sup></td></tr><tr><td>gas/dust</td><td>100</td></tr><tr><td>CO/H<sub>2</sub></td><td>10<sup>-4</sup></td></tr><tr><td>β</td><td>1.25</td></tr><tr><td>γ</td><td>-1.0</td></tr><tr><td>R<sub>c</sub> (au)</td><td>40</td></tr><tr><td>h (au)</td><td>15</td></tr><tr><td>M (M<sub>☉</sub>)</td><td>1.2</td></tr><tr><td>L (L<sub>☉</sub>)</td><td>0.96</td></tr><tr><td>T<sub>eff</sub> (K)</td><td>4500</td></tr></tbody></table>

Note. The scale height *h* is measured at 100 au.

geometry following that of Panić et al. (2009) and Pinte et al. (2018) with density and geometrical parameters qualitatively consistent with the observations. We emphasize that none of the disk models presented in this section are meant to be good a fit to the ALMA observations but to provide an idea of how different parameters such as the beam size of the observations and the inclination angle of the disk affect the TRD method. In these test models, we do not include photodesorption nor photodissociation effects. In Table 1, we summarize the stellar and disk parameters we used to create the synthetic disk models.

## B.1. The Ideal Model

The first synthetic model we created is a perfectly edge-on disk ($i = 90°$) with an extremely high spatial resolution ($0.′′01$ or $1.4$ au) and a spectral resolution that matches our observations ($<0.25$ kms$^{-1}$). Although we tried models with higher spectral resolution ($0.01$ kms$^{-1}$), we found no significant differences between these two. We applied the TRD technique to this synthetic disk and compared the reconstructed temperature map to the input temperature from MCFOST in Figure 16. The input MCFOST temperature and the reconstructed temperature maps are shown as black lines and color maps, respectively, with uneven steps of 5, 10, and 20 K to better display both temperatures. From the figure, we see that the temperature input and the retrieved temperature follow each other very closely at almost all radii and especially toward low elevations. The 15 K to 35 K temperature regions, for example, match at almost all radii and elevations. The hotter temperature of 55 K, on the other hand, starts deviating from the retrieved temperatures for radii larger than 200 au. The difference between the computed and the input temperatures is due to a lower density (lower optical depth, ultimately) of the gas, which happens at both very high elevations and larges distances from the center of the disk. In both cases, the reconstructed temperature is lower than the actual temperature of the disk. This behavior is expected, and it follows from the definition of brightness temperature, which is known to coincide with the kinetic temperature of the gas when the optical depth of a line is sufficiently large. However, it is important to note that the models presented here do not have any included noise; therefore, the optically thinner gas, often

<figure><img src="image_18.png" /><figcaption>**Figure 16.** Comparison between the TRD method applied to the synthetic disk model and the input disk temperature parameters. The model (shown in color) was created with the parameters defined in Table 1 and convolved with a $0.′′01$ beam. The MCFOST input temperature, shown by the black lines, ranges from 15 to 55 K.</figcaption></figure>

missed in the observations due to S/N detection limits, might not present a problem when using the TRD technique.

## B.2. TRD with Different Disk Inclinations

The probability of encountering a protoplanetary disk that is perfectly edge-on is very small. It is therefore important to consider the cases where there are small deviations in the inclination of the disk toward the observers. To explore the effects of inclining the disk’s out-of-edge-on position, we have computed two new MCFOST models with the same disk and stellar parameters as defined in Table 1, but this time we adopted inclinations of $i = 85°$ and $i = 80°$.

In Figure 17, we plotted these two new models with inclinations of $i = 85°$ (middle) and $i = 80°$ (right panel), in addition to the default model with $i = 90°$ (left panel). This time, both sides of the disk were plotted so that the inclination effects can be visualized for the regions above and below the midplane. For the regions above the midplane, we see that as the disk is declined out of the edge-on position, the farside of the disk becomes visible and lies above the photosphere of the disk overlapping with the nearside. When the disk is inclined at $i = 85°$, it is not possible to easily distinguish between the far- and the nearside of the disk in the reconstructed temperature map; thus, the emission from the farside is mixed with the emission from the photosphere on the nearside, resulting in an increase on the measured temperature at higher elevations and also an increase in the scale height at all radii. When the disk is further declined to $i = 80°$, the near- and farside of the disk are clearly separated, and the measured peak temperature returns to its original value. The layers below the midplane are also affected by changes in inclination. As the disk goes from $i = 90°$ to $i = 80°$, the nearside gets hidden behind the disk’s midplane, causing the light emitted from the bottom layer to arise from a place closer to the midplane (because of optical depth effects), and therefore we see a cooler temperature.

Although the TRD technique takes into consideration the inclination angle of the disk, the temperature reconstruction is affected by the optical depth. This means, for example, that for a disk that is inclined at $i = 80°$, the layer farther from the
<figure><img src="image_19.png" /><figcaption>**Figure 17.** Demonstration of the disk inclination variation when using the TRD method. From left to right, we plotted the synthetic disk at inclination angles of $i = 90^{\circ}$, $i = 85^{\circ}$, and $i = 80^{\circ}$. The disk's apparent temperature changes with changing inclination in both the layers above and below the midplane. The black contour lines correspond to the input MCFOST temperature.</figcaption></figure>

observer (bottom layer in our Figure 17) would not appear to be as hot as the closer layer (top layer), and this is not because the reconstruction is imperfect. Instead, this relates to the position in the disk where the optical depth of the line becomes much greater than unity. The same effect also causes one side of the disk (top and bottom layers) to appear larger than the other side. Although this may seem to be the case from Figure 17, the reality is that we are comparing different temperature contours, which are arbitrarily selected. If one allows the temperatures to range from much lower to much higher values, then the disks would have the same sizes, but there would be temperature differences between them.

### B.3. TRD with Different Spatial Resolution

Now we turn to the problem of how the spatial resolution affects the reconstructed temperature distribution of a protoplanetary disk. We use the same disk parameter as in the previous section (see Table 1) and keep the inclination to $i = 90^{\circ}$. To simulate observations at different spatial resolutions, we convolved the synthetic disk model with a two-dimensional circular Gaussian profile with an FWHM that varies from 0.01 to 0.3.

The spatial resolution effect is the most important limitation in our study of the temperature structure of a disk. It smears out the flux of the disk, limiting our ability to detecting high-temperature regions confined in small-scale structures. In

Figure 18, we illustrate this effect by convolving the same model shown in Figure 16 with different synthetic beam sizes. As we convolve the models with larger beam sizes, the retrieved peak temperature decreases, and the small-scale emission layers disappear. When the spatial resolution of the observations is degraded, the retrieved peak temperature shifts toward the center of the disk. Additionally, the disk appears to be thicker in the vertical direction extending beyond the nominal disk's photosphere and also it becomes thicker toward the bottom part of the disk, where most of the CO has frozen-out. When the spatial resolution of the observations (for a disk of this size) is better than 0.1, the convolution effects are much milder, and the retrieved temperature is comparable to the true temperature in the disk. It must be kept in mind, however, that even for very high spatial resolution observations, there are practical limits on the resolution of the reconstruction temperature maps. As the TRD method depends on the average of multiple pixels at a given radial position, on the shear Keplerian velocity, and on the turbulent velocity of the disk (as explained in Section 4).

To conclude, we would like to emphasize that if the spatial resolution of the observations was poor enough, the beam smearing effect would largely exceed the inclination effects. With the data we have presented, we have to acknowledge that inside of 125–150 au, we are likely to be affected by this effect; however, the outer disk is much less affected.
<figure><img src="image_20.png" /><figcaption>**Figure 18.** Comparison between the temperature retrieved model convolved with different beam sizes. The top left panel shows a model convolved with a 0.′′01 beam (same as Figure 16). The top right, bottom left, and bottom right panels show the same model convolved with 0.′′1, 0.′′2, and 0.′′3 beam sizes, respectively. The contours correspond to the input MCFOST temperature.</figcaption></figure>

## ORCID iDs

C. Flores https://orcid.org/0000-0002-8591-472X

G. Duchêne https://orcid.org/0000-0002-5092-6464

S. Wolff https://orcid.org/0000-0002-9977-8255

M. Villenave https://orcid.org/0000-0002-8962-448X

K. Stapelfeldt https://orcid.org/0000-0002-2805-7338

J. P. Williams https://orcid.org/0000-0001-5058-695X

C. Pinte https://orcid.org/0000-0001-5907-5179

D. Padgett https://orcid.org/0000-0001-5334-5107

M. S. Connelley https://orcid.org/0000-0002-8293-1428

G. van der Plas https://orcid.org/0000-0001-5688-187X

F. Ménard https://orcid.org/0000-0002-1637-7393

M. D. Perrin https://orcid.org/0000-0002-3191-8151

## References

Alcalá, J. M., Manara, C. F., Natta, A., et al. 2017, A&A, 600, A20

Alcalá, J. M., Natta, A., Manara, C. F., et al. 2014, A&A, 561, A2

Allen, L. E., & Strom, K. M. 1995, AJ, 109, 1379

Andrews, S. M. 2020, ARA&A, 58, 483

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

Ansdell, M., Williams, J. P., Trapman, L., et al. 2018, ApJ, 859, 21

Ansdell, M., Williams, J. P., van der Marel, N., et al. 2016, ApJ, 828, 46

Appenzeller, I., Bertout, C., & Stahl, O. 2005, A&A, 434, 1005

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

Avenhaus, H., Quanz, S. P., Garufi, A., et al. 2018, ApJ, 863, 44

Barradoy Navascués, D., & Martín, E. L. 2003, AJ, 126, 2997

Barrière-Fouchet, L., Gonzalez, J. F., Murray, J. R., Humble, R. J., & Maddison, S. T. 2005, A&A, 443, 185

Beckwith, S. V. W., & Sargent, A. I. 1993, ApJ, 402, 280

Birnstiel, T., & Andrews, S. M. 2014, ApJ, 780, 153

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

Boyden, R. D., & Eisner, J. A. 2020, ApJ, 894, 74

Chauvin, G., Ménard, F., Fusco, T., et al. 2002, A&A, 394, 949

Christiaens, V., Casassus, S., Perez, S., van der Plas, G., & Ménard, F. 2014, ApJL, 785, L12

Cieza, L. A., Ruíz-Rodríguez, D., Hales, A., et al. 2019, MNRAS, 482, 698

Cleeves, L. I., Öberg, K. I., Wilner, D. J., et al. 2016, ApJ, 832, 110

Cox, E. G., Harris, R. J., Looney, L. W., et al. 2017, ApJ, 851, 83

Cushing, M. C., Vacca, W. D., & Rayner, J. T. 2004, PASP, 116, 362

Dong, R., Fung, J., & Chiang, E. 2016, ApJ, 826, 75

Doppmann, G. W., Greene, T. P., Covey, K. R., & Lada, C. J. 2005, AJ, 130, 1145

Draine, B. T., & Bertoldi, F. 1996, ApJ, 468, 269

Duchêne, G., Becker, A., Yang, Y., et al. 2017, MNRAS, 469, 1783

Dullemond, C. P., Isella, A., Andrews, S. M., Skobleva, I., & Dzyurkevich, N. 2020, A&A, 633, A137

Dutrey, A., Guilloteau, S., Piétu, V., et al. 2017, A&A, 607, A130

Edwards, S., Fischer, W., Hillenbrand, L., & Kwan, J. 2006, ApJ, 646, 319

Evans, N. J., II, Dunham, M. M., Jørgensen, J. K., et al. 2009, ApJS, 181, 321

Facchini, S., Birnstiel, T., Bruderer, S., & van Dishoeck, E. F. 2017, A&A, 605, A16

Falcón-Barroso, J., Sánchez-Blázquez, P., Vazdekis, A., et al. 2011, A&A, 532, A95

Favre, C., Cleeves, L. I., Bergin, E. A., Qi, C., & Blake, G. A. 2013, ApJL, 776, L38

Feiden, G. A. 2016, A&A, 593, A99

Fischer, W., Kwan, J., Edwards, S., & Hillenbrand, L. 2008, ApJ, 687, 1117

Flores, C., Connelley, M. S., Reipurth, B., & Boogert, A. 2019, ApJ, 882, A75

Flores, C., Reipurth, B., & Connelley, M. S. 2020, ApJ, 898, 109

Folha, D. F. M., & Emerson, J. P. 2001, A&A, 365, 90

Gonzalez, J. F., Laibe, G., & Maddison, S. T. 2017, MNRAS, 467, 1984

Guilloteau, S., Piétu, V., Chapillon, E., et al. 2016, A&A, 586, L1
Hsieh, T.-H., & Lai, S.-P. 2013, ApJS, 205, 5

Huélamo, N., de Gregorio-Monsalvo, I., Macías, E., et al. 2015, A&A, 575, L5

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

Jacoby, G. H., Hunter, D. A., & Christian, C. A. 1984, ApJS, 56, 257

Jayawardhana, R., Luhman, K. L., D'Alessio, P., & Stauffer, J. R. 2002, ApJL, 571, L51

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

Kamp, I., Woitke, P., Pinte, C., et al. 2011, A&A, 532, A85

Koresko, C. D. 1998, ApJL, 507, L145

Krijt, S., Bosman, A. D., Zhang, K., et al. 2020, ApJ, 899, 134

Ladjelate, B., André, P., Könyves, V., et al. 2020, A&A, 638, A74

Le Borgne, J. F., Bruzual, G., Pelló, R., et al. 2003, A&A, 402, 433

Lombardi, M., Lada, C. J., & Alves, J. 2008, A&A, 489, 143

Long, F., Herczeg, G. J., Pascucci, I., et al. 2017, ApJ, 844, 99

Manara, C. F., Testi, L., Rigliaco, E., et al. 2013, A&A, 551, A107

Mannings, V., & Sargent, A. I. 1997, ApJ, 490, 792

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

Monin, J. L., Menard, F., & Duchene, G. 1998, A&A, 339, 113

Ortiz-León, G. N., Loinard, L., Kounkel, M. A., et al. 2017, ApJ, 834, 141

Panić, O., Hogerheijde, M. R., Wilner, D., & Qi, C. 2009, A&A, 501, 269

Pavlyuchenkov, Y., Semenov, D., Henning, T., et al. 2007, ApJ, 669, 1262

Perez, S., Casassus, S., Ménard, F., et al. 2015, ApJ, 798, 85

Pinte, C., Harries, T. J., Min, M., et al. 2009, A&A, 498, 967

Pinte, C., Ménard, F., Duchêne, G., et al. 2018, A&A, 609, A47

Pinte, C., Ménard, F., Duchêne, G., & Bastien, P. 2006, A&A, 459, 797

Pinte, C., van der Plas, G., Ménard, F., et al. 2019, NatAs, 3, 1109

Rayner, J. T., Cushing, M. C., & Vacca, W. D. 2009, ApJS, 185, 289

Rayner, J. T., Toomey, D. W., Onaka, P. M., et al. 2003, PASP, 115, 362

Sánchez-Blázquez, P., Peletier, R. F., Jiménez-Vicente, J., et al. 2006, MNRAS, 371, 703

Simon, M., Guilloteau, S., Beck, T. L., et al. 2019, ApJ, 884, 42

Sokal, K. R., Deen, C. P., Mace, G. N., et al. 2018, ApJ, 853, 120

Somers, G., & Pinsonneault, M. H. 2015, ApJ, 807, 174

Stapelfeldt, K. R., Duchêne, G., Perrin, M., et al. 2014, in IAU Symp. 299, Exploring the Formation and Evolution of Planetary Systems, ed. M. Booth, B. C. Matthews, & J. R. Graham (Cambridge: Cambridge Univ. Press), 99

Stapelfeldt, K. R., Krist, J. E., Ménard, F., et al. 1998, ApJL, 502, L65

Teague, R., Jankovic, M. R., Haworth, T. J., Qi, C., & Ilee, J. D. 2020, MNRAS, 495, 451

Thanathibodee, T., Calvet, N., Herczeg, G., et al. 2018, ApJ, 861, 73

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

Vacca, W. D., Cushing, M. C., & Rayner, J. T. 2003, PASP, 115, 389

Villenave, M., Ménard, F., Dent, W. R. F., et al. 2020, A&A, 642, A164

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

Wolff, S. G., Duchêne, G., & Stapelfeldt, K. R. 2021, AJ, 161, 238

Zhang, K., Bergin, E. A., Blake, G. A., Cleeves, L. I., & Schwarz, K. R. 2017, NatAs, 1, 0130

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

Zhang, K., Schwarz, K. R., & Bergin, E. A. 2020, ApJL, 891, L17