Mapping Protoplanetary Disk Vertical Structure with CO Isotopologue Line Emission

We searched the ALMA archive for observations of protoplanetary disks with inclinations of 30°–75° that covered several CO isotopologues in one or more lines, namely CO, ¹³CO, and C¹⁸O, and J = 2–1 and J = 3–2. We restricted our search to those sources observed at sufficiently high angular resolutions (≲03), line sensitivities (a few kelvins), and velocity resolutions (≲0.25 km s⁻¹) necessary to derive emission surfaces. After excluding previously published sources (Law et al. 2021a, 2022; Paneque-Carreño et al. 2021, 2022, 2023; Rich et al. 2021), we identified four such sources—the disks around DM Tau, Sz 91, LkCa 15, and HD 34282.

All data were obtained from the ALMA archive, except for the J = 3–2 lines of CO, ¹³CO, and C¹⁸O in the LkCa 15 disk (Jin et al. 2019) and for CO J = 3–2 in the HD 34282 disk (van der Plas et al. 2017), which were taken from previously published ALMA observations. To achieve the necessary data quality, we often combined multiple programs and executions. Each observational program is listed in Table 2 and described in detail in Appendix A.

Each archival project was initially calibrated by ALMA staff using the ALMA calibration pipeline and the required version of CASA (McMullin et al. 2007). Subsequent self-calibration was performed using CASA v5.4.0 for all sources, except for Sz 91, where we were unable to derive solutions that improved image quality.

Our self-calibration strategy closely followed that of the MAPS ALMA Large Program, which is described in detail by Öberg et al. (2021). We created pseudo-continuum visibilities by flagging line emission in each spectral window, which were then combined with continuum-only spectral windows, when available. We then averaged down these data into 125 MHz channels, imaged each execution block, and measured the phase centers with the imfit task in CASA. To account for any source proper motions or atmospheric/instrumental effects between observations, we aligned each execution block to a common phase center using the fixvis and fixplanets tasks to apply the necessary phase shifts and assign common labels to the phase centers, respectively.

We first self-calibrated the aligned, short-baseline data. When multiple short-baseline observations were available, all executions were initially concatenated together. Then, we concatenated the self-calibrated short-spacing data with the long-baseline data, and the combined visibilities were self-calibrated together. We considered phase solution intervals beginning at infinity and decreasing to 60 s, or stopping sooner if the signal-to-noise ratio (S/N) of the data did not improve. When it resulted in a S/N improvement, a single round of amplitude self-calibration was performed on the combined visibilities with an infinite solution interval. The self-calibration typically improved the continuum S/N by a factor of 2–3. We ultimately applied the resulting calibration solutions to the unflagged and spectrally nonaveraged visibilities, before subtracting the continuum with a first-order polynomial using the uvcontsub task.

We then switched to CASA v6.3.0 for all imaging. We used tclean to produce images of the J = 2–1 lines of CO, ¹³CO, and C¹⁸O for each source with Briggs weighting and Keplerian masks generated with the keplerian_mask (Teague 2020) code. Each mask was based on the stellar+disk parameters listed in Table 1 and was visually inspected to ensure that it contained all emission present in the channel maps. If required, manual adjustments to mask parameters were made, e.g., maximum radius and beam convolution size. Briggs robust parameters were chosen manually to prioritize spatially resolved and well-defined line emission in the channel maps (see Appendix B). Channel spacings ranged from 0.08 to 0.25 km s⁻¹ depending on the source and line. For several images, we also applied Gaussian uv-tapers of 0070–0150 to improve sensitivity to larger-scale emission. All images were made using the multiscale deconvolver with pixel scales of [0, 5, 15, 25], except for the CO J = 2–1 image in LkCa 15, which used [0, 8, 16, 32, 64]. All images were CLEANed down to a 4σ level, where σ was the rms measured in a line-free channel of the dirty image. Table 2 summarizes all image properties, including the ALMA project codes from which each image was generated.

Table 1. Stellar and Disk Characteristics

Source	Spectral	Distance a	Incl.	PA	M* b	L*	Age c	vsys b	Cloud
	Type	(pc)	(°)	(°)	(M⊙)	(L⊙)	(Myr)	(km s−1)	Contam.
DM Tau	M1 [1]	143	${36.0}_{-0.09}^{+0.12}$ [2]	154.5 ± 0.24 [3]	0.53 ± 0.02	0.24 [1] d	3–7 [1]	5.99 ± 0.04	...
Sz 91	M0 [4]	158	49.7 ± 0.2 [5]	17.0 ± 0.26 [3]	0.54 ± 0.07	0.26 ± 0.02 [5]	3–7 [5,6]	3.39 ± 0.04	Moderate
LkCa 15	K5 [7]	157	50.2 ± 0.03 [8]	61.9 ± 0.04 [8]	1.20 ± 0.07	${1.05}_{-0.21}^{+0.27}$ [7]	1–5 [1,7]	6.28 ± 0.04	...
HD 34282	A0-A1 [10]	306	59.3 ± 0.4 [9]	117.1 ± 0.3 [9]	1.69 ± 0.07	14.5 ± 0.7 [10]	4–7 [11]	−2.35 ± 0.01	...

Notes.

^a All distances are from Gaia Early Data Release 3 (Gaia Collaboration et al. 2021; Bailer-Jones et al. 2021). ^b Dynamical stellar masses and systemic velocities (in the LSR frame) are derived in this work and represent the mean values computed from all available CO isotopologues (see Appendix D). ^c Stellar ages are likely uncertain by at least a factor of 2. ^d Pegues et al. (2020) do not provide uncertainties on the stellar luminosity of DM Tau.

References: (1) Pegues et al. (2020), (2) Flaherty et al. (2020), (3) Law et al. (2022), (4) Romero et al. (2012), (5) Maucó et al. (2020), (6) Tsukagoshi et al. (2019), (7) Donati et al. (2019), (8) Facchini et al. (2020), (9) van der Plas et al. (2017), (10) Guzmán-Díaz et al. (2021), (11) Vioque et al. (2018).

Download table as:  ASCIITypeset image

Table 2. CO Isotopologue Image Cube Properties

Source	Beam	JvM a	robust	uv-taper	Chan. δv	rms b	ALMA
Transition	('' × '', deg)			('')	(km s−1)	(mJy beam−1, K)	Project Code(s)
DM Tau
CO J = 2–1	0.13 × 0.12, −37.9	0.37	0.5	0.075	0.08	1.0, 1.5	2013.1.00498.S, 2016.1.00724.S, 2017.1.01460.S
13CO J = 2–1	0.20 × 0.18, −46.1	0.48	0.5	0.150	0.08	1.5, 1.0	2016.1.00724.S, 2017.1.01460.S
C18O J = 2–1	0.20 × 0.18, −44.0	0.49	0.5	0.150	0.08	1.1, 0.8	2016.1.00724.S, 2017.1.01460.S
Sz 91
CO J = 2–1	0.23 × 0.16, 17.5	0.41	0.5	...	0.16	2.1, 1.3	2013.1.00663.S, 2013.1.01020.S, 2015.1.01301.S
13CO J = 2–1	0.24 × 0.17, 21.9	0.41	0.5	...	0.20	2.2, 1.4	2013.1.00663.S, 2013.1.01020.S, 2015.1.01301.S
C18O J = 2–1	0.41 × 0.38, 30.0	0.21	2.0	0.300	0.20	0.8, 0.1	2013.1.00663.S, 2013.1.01020.S, 2015.1.01301.S
LkCa 15
CO J = 2–1	0.37 × 0.27, −9.2	0.73	0.0	...	0.20	2.3, 0.4	2013.1.00226.S, 2018.1.01255.S
13CO J = 2–1	0.14 × 0.10, −18.8	0.38	0.0	0.070	0.25	0.7, 1.3	2018.1.00945.S
C18O J = 2–1	0.18 × 0.14, 10.0	0.14	2.0	...	0.25	0.2, 0.2	2018.1.00945.S
HD 34282
CO J = 2–1	0.11 × 0.09, −69.7	0.21	−2.0	...	0.08	0.4, 1.0	2015.1.00192.S, 2017.1.01578.S
13CO J = 2–1	0.12 × 0.10, −69.9	0.16	−2.0	...	0.20	0.2, 0.5	2015.1.00192.S, 2017.1.01578.S
C18O J = 2–1	0.12 × 0.10 −70.2	0.16	−2.0	...	0.20	0.2, 0.4	2015.1.00192.S, 2017.1.01578.S

Notes.

^a The ratio of the CLEAN beam and dirty beam effective area used to scale image residuals to account for the effects of non-Gaussian beams. See Section 2.2 and Jorsater & van Moorsel (1995) and Czekala et al. (2021) for further details. ^b rms values were calculated for the corresponding image cube channel (δv) and brightness temperatures were calculated assuming the Rayleigh–Jeans limit.

Download table as:  ASCIITypeset image

After generating the image cubes, we applied the "JvM" correction proposed in Jorsater & van Moorsel (1995) and described in more detail in Czekala et al. (2021). This correction scales the image residuals by a factor, , equal to the ratio of the effective areas of the CLEAN beam and dirty beam, to be in units consistent with the CLEAN model. Table 2 lists all values. While we used the JvM-corrected images here, we have also verified that line emission surfaces extracted from either the JvM- or non-JvM-corrected images yield consistent results.

We also imaged the non-continuum-subtracted line data and the continuum data by adopting the same imaging parameters as the line-only emission image cubes. The non-continuum-subtracted image cubes were required for the calculation of gas temperatures (Section 4.5) and provided an additional check that continuum subtraction did not influence the extracted emission surfaces, while the continuum images were used to define the disk centers. The line-only and line+continuum image cubes, as well as all zeroth moment maps, are publicly available. ²²

All four sources in our sample host transition disks and span a range in stellar properties, such as masses (0.53–1.69 M_⊙), spectral types (M0–A1) and bolometric luminosities (0.24–14.5 L_⊙), and disk characteristics, such as overall CO emission radial extents (≈400–1000 au). Table 1 shows a summary of source characteristics.

Figure 1 shows an overview of the disk sample in CO isotopologue velocity-integrated intensity, or "zeroth moment," maps, and in millimeter continuum emission. All continuum images were generated from the corresponding programs from which the line images were produced. We generated zeroth moment maps of line emission from the image cubes using bettermoments (Teague & Foreman-Mackey 2018) and closely followed the procedures outlined in Law et al. (2021b). We did not use a flux threshold for pixel inclusion, i.e., sigma clipping, to ensure accurate flux recovery and used the same Keplerian masks employed during CLEANing.

Zoom In Zoom Out Reset image size

The Sz 91 disk exhibits moderate cloud contamination in CO J = 2–1 between v_LSR ≈ 4 and 7 km s⁻¹, in which the ambient cloud significantly absorbs disk line emission with overlapping velocities. This is identified through visual inspection of channel maps (Appendix B) and manifests as a north–south spatial brightness asymmetry in images of the CO line emission (see Figure 1). Cloud obscuration was previously identified in this disk in a similar velocity range in the CO J = 3–2 line (Canovas et al. 2015; Tsukagoshi et al. 2019).

We derived vertical emission heights on a per-pixel basis directly from the line emission image cubes using the disksurf (Teague et al. 2021) Python code, based on the methodology presented in Pinte et al. (2018). We closely followed the methods outlined in Law et al. (2021a), which we briefly summarize below.

Before extraction, we first masked the image cubes with the same Keplerian masks used during CLEANing and manually excluded all channels where the front and back disk sides could not be clearly distinguished. We adopted an inclination and position angle (PA) for each disk (Table 1), which yields a deprojected radius r, emission height z, surface brightness I_ν , and channel velocity v for each pixel associated with the emitting surface.

We then filtered pixels based on priors of expected disk physical structure, i.e., removing those pixels with unphysically high z/r or large negative z values, as the emission must arise from at least the midplane. To avoid positively biasing our averages to nonzero z values, we did not remove those with small negative values, i.e., z/r > −0.1. For the HD 34282 disk, we instead only removed pixels with z/r < 0.05 to mitigate confusion due to the high inclination of this source and visually confirmed that the resulting surface was not artificially distorted by this cut. To avoid the misidentification of peaks due to noise, we also filtered pixels with low surface brightnesses, ranging from 1× rms (Sz 91) to 6× rms (HD 34282). The wide range in thresholds was a result of our heterogeneous sample with differing line sensitivities. Throughout this process, we prioritized obtaining the maximum number of robust emission surface pixels and at each step visually confirmed the fidelity of derived surfaces.

After completing these filtering steps, we binned the surfaces in two ways: (i) into radial bins equal to half of the FWHM of the beam major axis; and (ii) by computing a moving average with a minimum window size of 1/2× the beam major axis FWHM. While both methods are effective at reducing scatter in the extracted surfaces, the radially binned surfaces benefit from a uniform radial sampling, while the moving averages maintain a finer radial sampling that is sensitive to subtle vertical perturbations, e.g., from putative embedded planets. In some cases, we radially bin these surfaces further for visual clarity, but all quantitative analysis is performed with the original binning of each type of emission surface. All three types of line emission surfaces—individual measurements, radially binned, and moving averages—are made publicly available as data behind the figure.

We fitted parametric models to all line emission surfaces to facilitate comparisons with other observations and incorporation into models. For the majority of sources and lines, we use the form of an exponentially tapered power law, which describes the flared inner surfaces as well as the plateau and turnover region at larger radii. We adopted the same functional form as in Law et al. (2021a):

$$\begin{eqnarray}&&z(r)={z}_{0}\times {\left(\displaystyle \frac{r}{1^{\prime\prime} }\right)}^{\phi }\times \exp \left(-{\left[\displaystyle \frac{r}{{r}_{\mathrm{taper}}}\right]}^{\psi }\right).\end{eqnarray} \tag{ 1 }$$

However, in a few cases the derived surface was not well fitted by an exponentially tapered power-law profile, as no turnover of the emission surface was detectable (as discussed in detail in the following section). In those cases, we adopted a single power-law profile.

We used the Monte Carlo Markov Chain (MCMC) sampler implemented in emcee (Foreman-Mackey et al. 2013) to estimate the posterior probability distributions for z₀, ϕ, r_taper, and ψ. Each ensemble uses 64 walkers with 1000 burn-in steps and an additional 500 steps to sample the posterior distribution function. Table 3 shows the median values of the posterior distribution, with uncertainties given as the 16th and 84th percentiles. We stress that these represent the statistical uncertainties associated with fitting a functional form to the extracted data points and do not include the systematic uncertainties associated with extracting those data points (e.g., uncertainties in disk inclination and PA). We also list the maximum radius (r ${}_{\mathrm{fit},\ \max }$) that we considered for each surface fit.

Table 3. Parameters for CO Isotopologue Emission Surface Fits

Source	Line	Exponentially Tapered Power Law
		r ${}_{\mathrm{fit},\ \max }$ ('')	z0 ('')	ϕ	rtaper ('')	ψ
DM Tau	CO J = 2−1	2.70	${0.63}_{-0.01}^{+0.01}$	${1.81}_{-0.02}^{+0.02}$	[2.00] a	${1.51}_{-0.03}^{+0.03}$
	13CO J = 2−1 b	1.30	${0.22}_{-0.00}^{+0.00}$	${1.14}_{-0.02}^{+0.02}$	...	...
Sz 91	CO J = 2−1 b	2.50	${0.16}_{-0.01}^{+0.01}$	${0.72}_{-0.05}^{+0.05}$	...	...
LkCa 15	CO J = 2−1	4.75	${0.25}_{-0.00}^{+0.00}$	${1.42}_{-0.05}^{+0.06}$	${3.65}_{-0.08}^{+0.07}$	${3.39}_{-0.25}^{+0.26}$
	13CO J = 2−1	2.25	${0.16}_{-0.02}^{+0.07}$	${1.27}_{-0.29}^{+0.44}$	${1.93}_{-0.42}^{+0.19}$	${3.22}_{-1.26}^{+1.59}$
	CO J = 3−2	5.00	${0.20}_{-0.01}^{+0.02}$	${2.27}_{-0.17}^{+0.22}$	${2.64}_{-0.36}^{+0.29}$	${2.02}_{-0.27}^{+0.29}$
	13CO J = 3−2	2.75	${0.17}_{-0.01}^{+0.02}$	${1.14}_{-0.13}^{+0.24}$	${2.60}_{-0.27}^{+0.14}$	${3.76}_{-1.46}^{+0.90}$
HD 34282	CO J = 2−1	1.50	${0.83}_{-0.05}^{+0.02}$	${1.37}_{-0.04}^{+0.04}$	${1.03}_{-0.02}^{+0.05}$	${1.41}_{-0.05}^{+0.07}$
	CO J = 3−2	2.00	${0.38}_{-0.02}^{+0.02}$	${0.80}_{-0.08}^{+0.10}$	${1.65}_{-0.04}^{+0.03}$	${4.36}_{-0.59}^{+0.66}$

Notes.

^a Due to parameter degeneracies encountered when fitting the extended, elevated emission, we fixed r_taper. ^b Single power-law profiles were used.

Download table as:  ASCIITypeset image

Figures 2 and 3 show the CO and, when available, ¹³CO emission surfaces, respectively, derived for the disks in our sample. In all cases, the CO surfaces lie at higher altitudes than those of ¹³CO, consistent with line optical depths of ∼1 being reached at deeper layers for the less abundant ¹³CO isotopologue. The inability to derive a ¹³CO emission surface in the Sz 91 disk was driven by both low S/N and insufficient angular resolution, while for the HD 34282 disk, despite having the best angular resolution in our sample, the larger source distance prevented us from spatially resolving the ¹³CO emitting heights. We do not resolve the vertical structure of the C¹⁸O emission in any disk, which suggests that C¹⁸O arises at or near the disk midplanes.

Zoom In Zoom Out Reset image size

As shown in Figure 2, all disks show elevated CO emission with maximum absolute heights of ≈100 au, with the exception of the Sz 91 disk, which has CO heights of only ≈50 au. There is considerable variation in the typical z/r values spanned by the surfaces from ≈0.1 to ≳0.5. The DM Tau disk has by far the most elevated CO emission surface, while Sz 91 hosts the flattest disk. Figure 3 shows that for those sources where we could derive ¹³CO surfaces, the ¹³CO surface traces between z/r ≳ 0.2 (DM Tau) and z/r = 0.1–0.15 (LkCa 15).

Zoom In Zoom Out Reset image size

In all lines, the LkCa 15 disk shows the characteristic emission surface profile of a sharply rising surface in the inner disk, followed by a gradual flattening and then turnover at large radii which is, presumably, due to decreasing gas surface densities. The CO J = 2–1 surface in the DM Tau disk, and to a lesser degree the HD 34282 disk, does not clearly show a turnover at large radii and instead shows a prolonged plateau phase. The DM Tau disk is particularly reminiscent of the IM Lup disk, which showed similar diffuse, elevated CO emission at large radii (Pinte et al. 2018; Law et al. 2021a). The Sz 91 disk also does not show a clear turnover and instead appears to rise out to large radii, with a possible upturn at ≈400 au. This occurs at the radial location where we have insufficient S/N to be certain of this feature, but note that the more sensitive CO J = 3–2 data (Figure 2) show elevated, diffuse emission at similar radii (Tsukagoshi et al. 2019; Law et al. 2022).

To better illustrate the three-dimensional (3D) geometry of these line emission surfaces, we overlay the inferred CO emission surfaces on CO peak line intensity maps in Figure 4. All peak intensity maps were generated using the "quadratic" method of bettermoments with the full Planck function.

Zoom In Zoom Out Reset image size

In addition to the vertical structure, we also characterized the radial morphology of the CO isotopologue emission in each disk. We first generated a set of radial line intensity profiles (Section 4.2.1) and then used these profiles to catalog the properties of all substructures present in each source (Section 4.2.2) and to compute the radial size of the CO isotopologue line emission (Section 4.2.3).

4.2.1. Radial Line Intensity Profile Generation

We computed radial profiles using the radial_profile function in the GoFish python package (Teague 2019a) to deproject the zeroth moment maps. During deprojection, we incorporated the derived emission surfaces listed in Table 3. This is particularly important for highly elevated surfaces, e.g., CO, in order to derive accurate radial locations of line emission substructures (e.g., Law et al. 2021b; Rosotti et al. 2021). We generated a set of radial profiles from azimuthally averaged profiles to those only including 15°, 30°, and 45° wedges along the disk major axis. Azimuthally averaged profiles are most sensitive to weaker emission at large radii, while profiles extracted from narrow wedges possess a higher effective spatial resolution, given adequate S/N. We chose profiles by visual inspection by prioritizing both the sharpness of substructures, which may otherwise be smoothed away in azimuthally averaged profiles, while maintaining profile fidelity for sources and lines with lower S/N. For the Sz 91 disk, we extracted the CO J = 2–1 radial profile in a ±90° wedge in the northern part of the disk to avoid substantial cloud obscuration, following Law et al. (2022). Figure 5 shows the resultant radial profiles.

Zoom In Zoom Out Reset image size

4.2.2. Radial Substructures

4.2.3. Gas Disk Sizes

For several sources, we have CO isotopologue emission surfaces for both the J = 2–1 and J = 3–2 lines, either derived directly in this work as for the HD 34282 and LkCa 15 disks or by combining our results with those in the literature. CO rotational lines with sufficiently different excitation properties are expected to probe distinct regions within a disk, due to changing excitation conditions combined with strong radial and vertical temperature gradients across disks (e.g., Dartois et al. 2003; Bruderer et al. 2012; Fedele et al. 2016). We can now directly assess if any such excitation-related effects, i.e., differing emission heights between lines, are present in a small sample of disks.

Figure 6 shows those sources where we have emission surfaces of both the J = 2–1 and J = 3–2 lines in either CO and/or ¹³CO. In this sample, we also included the HD 163296 and GM Aur disks, which have CO and ¹³CO J = 2–1 surfaces, respectively, derived as part of the MAPS program (Law et al. 2021a), as well as ¹³CO J = 3–2 in GM Aur (Schwarz et al. 2021). We also extracted the CO J = 3–2 surface of HD 163296 from ALMA science verification data (e.g., Rosenfeld et al. 2013; de Gregorio-Monsalvo et al. 2013) using the same methods described in Section 3.1. While the ALMA science verification data of CO J = 3–2 in the HD 163296 disk and the ¹³CO J = 3–2 data in the GM Aur disk have considerably coarser beams than used in this work, we do not expect these larger beam sizes to significantly influence the derived surfaces for these disks (see Appendix D in Law et al. 2021a).

Zoom In Zoom Out Reset image size

In all disks and in either CO or ¹³CO, the J = 3–2 and J = 2–1 lines show consistent line emission heights. In a few cases, there are hints that the J = 3–2 surface is more elevated than that of the J = 2–1, particularly at large radii. This is most notable in the ¹³CO emission surfaces of the GM Aur disk and, to a lesser degree, in the LkCa 15 disk, but in neither case do we consider these putative differences conclusive. The overall similarity in line emission heights is perhaps not surprising, given that these surfaces are derived from consecutive low-J lines that span a relatively narrow range in excitation conditions, i.e., the upper-state energy of the CO J = 3–2 line (E_u ≈ 33.2 K) is only twice that of the J = 2–1 line (E_u ≈ 16.6 K). This is consistent with the traditional assumption that low (J_u < 6) CO rotational lines are tracing the cold gas in the outer disk. Thus, as the line emission is originating from similar heights, we expect both lines to be tracing gas at approximately the same temperature in each disk. We return to this in more detail in Section 4.5.

The vertical distribution of micron-sized dust grains in disks should be largely correlated with the vertical gas structure as small dust grains are expected to be strongly coupled to the gas. A detailed understanding of this relationship is critical, as dust evolution from micron-sized grains to pebbles is an important step in the planet formation process. Here, we compare observations of line emission surfaces with scattering heights measured in the near-IR (NIR; e.g., Ginski et al. 2016; Monnier et al. 2017; Avenhaus et al. 2018; Garufi et al. 2020).

Both the HD 34282 (de Boer et al. 2021; Quiroz et al. 2022) and LkCa 15 (Thalmann et al. 2015, 2016; Oh et al. 2016; Currie et al. 2019) disks have existing NIR polarimetric imaging that show well-defined rings, from which direct estimates of the small dust scattering heights can be inferred. ²³ The NIR heights of the rings in the HD 34282 disk were derived by de Boer et al. (2021), but, to our knowledge, no such heights have been reported for the LkCa 15 disk. Instead, we estimated the scattering height of the outer NIR ring in the LkCa 15 disk from the existing SPHERE images of Thalmann et al. (2016), following the procedures outline in Rich et al. (2021) and described in detail in Appendix C. We do not attempt to derive a scattering height for the inner disk (≲30 au) component.

Figure 7 shows these NIR heights compared to the CO isotopologue line emission surfaces derived in this work. In both the LkCa 15 and HD 34282 disks, the NIR ring heights lie at approximately the same height as the CO line emission surfaces. In general, but not exclusively, the scattering surfaces in disks have been found to lie lower than the CO surface and are instead often approximately vertically colocated with ¹³CO line emitting heights (Law et al. 2021a, 2022; Rich et al. 2021). NIR rings within 100 au, however, such as in the HD 163296 disk (Law et al. 2021a; Rich et al. 2021), are located at nearly the same height as CO, which is similar to what is observed here in both the LkCa 15 and HD 34282 disks.

Zoom In Zoom Out Reset image size

Overall, these results suggest an unexplored diversity in the relative distributions of small dust scattering surfaces and line emitting layers in disks. Comparisons between directly mapped line emission surfaces in disks with known scattering heights provide a powerful empirical way to probe multiple disk components at the same time, especially in light of the growing number of high-angular-resolution ALMA observations of CO isotopologue emission lines in those disks with known NIR features.

4.5.1. Calculating Gas Temperatures

4.5.2. Radial Temperature Profiles

Figure 8 shows the CO and, when available, ¹³CO radial temperature distributions along the emission surfaces. Derived brightness temperatures are generally consistent with expectations based on stellar luminosity and spectral classes, with the disk around Herbig Ae star HD 34282 showing warmer temperatures than those around the T Tauri stars.

Zoom In Zoom Out Reset image size

Many sources show drops in brightness temperatures toward their inner disks, as seen in Figure 8. At the smallest radii, this is primarily due to beam dilution as the emitting area becomes comparable to or smaller than the angular resolution of the observations. For instance, in the Sz 91 disk, the higher-angular-resolution CO J = 3–2 observations show an increasing temperature profile toward the inner disk, which indicates that the observed turnover at ∼100 au in the CO J = 2–1 temperature profile is primarily a beam effect. However, some sources show flattening or declining temperature profiles (toward the central star) at radii beyond that of the beam size. There are several possible explanations, e.g., CO gas depletion, dust absorption, or unresolved CO emission substructure. The dip in CO temperatures in the DM Tau disk is likely, at least in part, due to dust optical depth, since the ring at ∼20 au shows τ ≳ 1 (Hashimoto et al. 2021), while it is not clear what is driving the sharp change in the ¹³CO gas temperature. The inner disk of LkCa 15 shows marginally optically thick (τ ∼ 0.5) dust (Facchini et al. 2020) and decreased gas column density (Leemker et al. 2022), which both likely contribute to the observed temperature drops in Figure 8.

Beyond these inner dips, the radial temperature profiles are generally smooth, with a few exceptions. We identify a local temperature bump in the CO J = 2–1 temperatures of both the Sz 91 and HD 34282 disks. The feature in the HD 34282 disk is located at a radius of ≈280 au, while the enhancement in the Sz 91 disk occurs at approximately 310 au, which is coincident with a similar enhancement at ≈300 au seen in the CO J = 3–2 temperature profile (Law et al. 2022). Both features have relatively broad widths of ≈125–150 au and are listed in Table 6. A tentative bump is also seen at 300 au in the CO J = 3–2 line in the HD 34282 disk, which is consistent with the location of the enhancement in the J = 2–1 line. However, due to the larger beam size, it is difficult to confirm the reality of this feature.

Table 6. Radial Temperature Profile Fits

Source	Line	rfit,in (au)	rfit,out (au)	T100 (K)	q	Feat. a
DM Tau	CO J = 2−1	80	401	38 ± 0.5	0.56 ± 0.02	...
	13CO J = 2−1	140	393	26 ± 0.4	0.42 ± 0.02	...
Sz 91	CO J = 2−1 b	90	270	37 ± 1.2	0.76 ± 0.07	B310
	CO J = 3−2 b , c	69	250	47 ± 0.7	0.70 ± 0.03	B300
LkCa 15	CO J = 2−1	90	743	31 ± 1.0	0.46 ± 0.03	...
	13CO J = 2−1	98	352	23 ± 0.4	0.39 ± 0.02	...
	CO J = 3−2	140	781	33 ± 0.8	0.52 ± 0.02	...
	13CO J = 3−2	116	437	19 ± 0.7	0.25 ± 0.04	...
HD 34282	CO J = 2−1 b	105	238	52 ± 0.5	0.26 ± 0.01	B280
	CO J = 3−2 b	190	575	58 ± 6.6	0.34 ± 0.09	B300

Notes.

^a Local temperature bumps (B) are labeled according to their approximate radial location in astronomical units. ^b Temperature bumps in the Sz 91 and HD 34282 disks are excluded in the power-law fits. ^c Temperature fit is from Law et al. (2022).

Download table as:  ASCIITypeset image

In most outer disk regions, the gas temperatures drop below the CO freeze-out temperature, indicating that our approach breaks down as we approach the low-density, partially optically thin outer disk.

The temperatures that we derive agree well with previous estimates in the same disks. The CO gas temperatures of the DM Tau disk are consistent with those inferred using previous lower-angular-resolution (≳03) CO J = 2–1 observations (Flaherty et al. 2020; Law et al. 2022), with peak gas temperatures inferred here being only 2–3 K greater than those seen previously. The CO and ¹³CO gas temperatures in the inner 200 au of the LkCa 15 disk, which were estimated using observations at ≈03 (Leemker et al. 2022), are a few kelvins warmer, but are generally consistent, with our estimates. The Sz 91 disk has a steep temperature profile in the inner disk, which was derived in the same way using CO J = 3–2 observations at 014 (Law et al. 2022) and, for comparison, we also show this temperature in Figure 8. Although the angular resolution of the CO J = 2–1 data used here is only ≈1.5 times larger that of the CO J = 3–2 data, we are not sensitive to warm gas in the inner 100 au due to beam dilution. For instance, gas temperatures derived from CO J = 3–2 observations are nearly 70 K at ≈50 au (Law et al. 2022). We also note that, at any particular radius, the CO J = 2–1 temperature is 5–8 K lower than that of the J = 3–2 profile, which, unlike in the DM Tau disk, suggests that the temperatures derived here are being modestly lowered by nonunity beam filling factors. While it is possible that in the Sz 91 disk the lower CO J = 2–1 temperatures relative to the J = 3–2 line may instead reflect an origin of the J = 2–1 line in colder disk material, this is unlikely given the similar emitting heights of both lines (Figure 6).

We fitted all temperature profiles with power-law profiles, parameterized by slope q and T₁₀₀, the brightness temperature at 100 au, i.e.,

$$\begin{eqnarray}&&T={T}_{100}\times {\left(\displaystyle \frac{r}{100\,\mathrm{au}}\right)}^{-q}.\end{eqnarray} \tag{ 2 }$$

We excluded all temperatures within 2 times the FWHM of the major axis of the synthesized beam as these are likely affected by beam dilution, as discussed above. We also excluded the temperature bumps in both the Sz 91 (> 270 au) and HD 34282 (≈240–380 au) disks. Each profile was then fit using the Levenberg–Marquardt minimization implementation in scipy.optimize.curve_fit. Table 6 lists the fitting ranges and derived parameters. Most sources are fit well by power-law profiles with q ≈ 0.4–0.6, with CO J = 2–1, 3–2 in HD 34282 and ¹³CO J = 3–2 being considerably shallower (q ≈ 0.25). The Sz 91 disk shows the steepest profile with q = 0.76, which is consistent with the similarly steep CO J = 3–2 profile (q = 0.70) found in Law et al. (2022).

4.5.3. Two-dimensional Temperature Profiles

Combining all line data for each source, Figure 9 shows the thermal structure of the CO and, when available, ¹³CO emitting layers as a function of (r, z) for each source.

Zoom In Zoom Out Reset image size

For those sources where we have multiple CO isotopologue surfaces, which trace different heights in the same disk, we can construct a two-dimensional (2D) model of the temperature distribution. We adopt the same two-layer model used in Law et al. (2021a), which is similar to the one proposed by Dartois et al. (2003), and then subsequently modified with a different connecting term by Dullemond et al. (2020).

The midplane temperature, T_mid, and atmosphere temperature, T_atm, are assumed to have a power-law profile with slopes q_mid and q_atm, respectively:

$$\begin{eqnarray}&&{T}_{\mathrm{mid}}(r)={T}_{\mathrm{mid},0}{\left(r/100\,\mathrm{au}\right)}^{{q}_{\mathrm{mid}}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{T}_{\mathrm{atm}}(r)={T}_{\mathrm{atm},0}{\left(r/100\,\mathrm{au}\right)}^{{q}_{\mathrm{atm}}}.\end{eqnarray} \tag{ 4 }$$

Between the disk midplane and atmosphere, the temperature is smoothly connected using a tangent hyperbolic function:

$$\begin{eqnarray}&&{T}^{4}(r,z)={T}_{{\rm{mid}}}^{4}(r)+\displaystyle \frac{1}{2}\left[1+\tanh \left(\displaystyle \frac{\left.z-\alpha {z}_{q}(r\right)}{\left.{z}_{q}(r\right)}\right)\right]{T}_{{\rm{atm}}}^{4}(r),\end{eqnarray} \tag{ 5 }$$

where ${z}_{q}(r)={z}_{0}{\left(r/100\,\mathrm{au}\right)}^{\beta }$. The α parameter defines the height at which the transition in the tanh vertical temperature profile occurs, while β describes how the transition height varies with radius.

We fit the individual, per-pixel temperature measurements of CO and ¹³CO using emcee (Foreman-Mackey et al. 2013) to estimate the following seven parameters: T_atm,0, q_atm, T_mid,0, q_mid, α, z₀, and β. We used 256 walkers, which take 500 steps to burn in, and then an additional 5000 steps to sample the posterior distribution function. We took parameter values and uncertainties as the 50th, 16th, and 84th percentiles from the marginalized posterior distributions, respectively, and list them in Table 7.

Table 7. Summary of 2D Temperature Structure Fits

Source	Tatm,0 (K)	Tmid,0 (K)	qatm	qmid	z0 (au)	α	β
DM Tau	${38}_{-0.4}^{+0.5}$	${26}_{-0.4}^{+0.3}$	−0.74${}_{-0.02}^{+0.02}$	−0.39${}_{-0.01}^{+0.02}$	${18}_{-0.9}^{+1.0}$	${2.14}_{-0.09}^{+0.09}$	−0.07${}_{-0.05}^{+0.05}$
LkCa 15	${35}_{-0.9}^{+0.9}$	${21}_{-0.3}^{+0.3}$	−0.59${}_{-0.03}^{+0.03}$	−0.29${}_{-0.01}^{+0.01}$	${17}_{-1.5}^{+1.5}$	${2.57}_{-0.2}^{+0.22}$	${0.10}_{-0.04}^{+0.04}$

Download table as:  ASCIITypeset image

In the case of the LkCa 15 disk, since both J = 3–2 and J = 2–1 lines trace the same vertical layers, we chose to use the J = 2–1 line, which has a higher angular resolution. We only considered those regions with well-constrained temperature data in our fits, as shown in Figure 10. We also included all temperature data points, even those below T_B < 20 K in our fits. As gas temperatures close to or below the CO freeze-out temperature indicate that the line emission is likely partially optically thin, and thus a lower limit on the true gas temperature, we note that the true gas temperatures at large radii or closer to the disk midplane may be underestimated in our model.

Zoom In Zoom Out Reset image size

Figure 10 shows the 2D fitted models in comparison with the data. For both disks, the residuals between the fitted model and measured temperatures are typically no more than 15%, with the exception of ¹³CO in the DM Tau disk, where our model overpredicts gas temperatures by ≈25%–30%. As the emitting surfaces do not provide direct constraints in the disk midplanes, the empirically derived T_mid should be treated with caution.

The derived 2D temperature structure of the DM Tau disk is consistent with that of Flaherty et al. (2020), who used a similar two-layer model but a single power law for both T_mid and T_atm (i.e., q = q_mid = q_atm) and fixed values of α and β. Our model better captures the warmer temperatures along the inner (< 150 au) rising part of the CO emission surface. We also find a similar vertical temperature distribution in the LkCa 15 disk to that reported by Jin et al. (2019), who used an iterative, radiative-transfer model to derive the 3D disk temperature structure. ²⁴ The primary difference is that the model of Jin et al. (2019) predicts a cooler disk midplane but, given our observational constraints as mentioned above, combined with the difficulty associated with deriving accurate midplane temperatures, this is not surprising.

To enable homogeneous, source-to-source comparisons, we first needed to compute the characteristic height of the newly derived emission surfaces. To do so, we adopted the definition introduced by Law et al. (2022), i.e., the mean of all z/r values interior to a cutoff radius of r_cutoff = 0.8 × r_taper, where r_taper is the fitted parameter from the exponentially tapered power-law profiles from Table 3. This serves as a convenient definition, as averaging within 80% of the fitted r_taper ensures that we only include the rising portion of the emission surfaces, which we also visually confirmed for all surfaces. Since the CO J = 2–1 surface in the Sz 91 disk was fit with a single power-law profile, we manually fixed r_cutoff = 300 au to include only the rising portion of the surface. Since some disks are more flared than others, i.e., z/r changes rapidly with radius, we also computed the 16th to 84th percentile range within these same radii as a proxy of the overall flaring of each disk. Table 8 lists the characteristic z/r, flaring ranges, and r_cutoff values for all sources in our sample.

Figure 11 shows these representative z/r values for our disks and a large sample of literature sources as a function of stellar host mass, mean gas temperature, and CO gas disk size. All literature sources have directly mapped emission surfaces with source parameters derived in the same way as the disks in this work, i.e., dynamically derived masses from rotation maps, mean gas temperatures from directly mapped surfaces, and disk sizes computed as the radius enclosing 90% of the total flux (Law et al. 2021a, 2021b, 2022; Paneque-Carreño et al. 2021; Rich et al. 2021; Teague et al. 2021; Veronesi et al. 2021).

Zoom In Zoom Out Reset image size

Assuming that CO line emission surfaces scale with gas pressure scale heights, we expect that $z/r\sim {M}_{* }^{-1/8}$ and z/r ∼ T^−1/6 (Law et al. 2022), as shown in orange in Figure 11. The size of the gas disk, R_CO, and z/r were also shown to be tightly positively correlated (Law et al. 2022), as is evident in Figure 11. Here, we add two new sources (LkCa 15 and HD 34282) to these trends and refine the characteristic z/r value of the CO J = 2–1 surface of the DM Tau disk using new higher-angular-resolution observations.

The CO emission height of the LkCa 15 disk is consistent with all previous trends, while the HD 34282 disk appears to be an outlier with a CO emission surface with a characteristic z/r roughly a factor of 2 larger than expected for its stellar mass, disk gas temperature, and CO emission extent. Although we revise the characteristic z/r in the DM Tau disk down by ≈20% from the value derived in Law et al. (2022), who used lower-angular-resolution data, it still shows the second-highest z/r of all known CO surfaces (with only the HD 34282 disk being more elevated).

Both stellar mass and mean temperature trends remain highly scattered and the weak trend suggested by scaling laws means that, even with additional sources, it is difficult to assess the reality of such trends. Observations of disks with more diverse properties, such as substantially warmer gas temperatures and around either lower-mass (<0.5 M_⊙) or higher-mass (>3 M_⊙) stars, are required to provide meaningful anchor points for these potential trends.

The correlation between R_CO and z/r remains strong and the addition of several sources with large CO disks allows us to better quantify the R_CO–z/r correlation. To do so, we used the Bayesian linear regression method of Kelly (2007) using the linmix Python implementation. ²⁵ We find a best-fit relation of z/r = (2.3 ± 0.6 × 10⁻⁴)R_CO + (0.15 ± 0.04) with a 0.04 scatter of the correlation (taken as the standard deviation σ of an assumed Gaussian distribution around the mean relation). We find a correlation coefficient of $\hat{\rho }=0.77$ and associated confidence intervals of (0.27, 0.99), which represent the median and 99% confidence regions, respectively, of the 2.5 × 10⁶ posterior samples for the regression. Overall, this is consistent with the fit reported by Law et al. (2022) but with a slightly less steep (≈30%) slope.

Although not shown in Figure 11, the radially extended (R_CO ≳ 1400 au) but modestly elevated (z/r ∼ 0.3) edge-on CO disk around Gomez's Hamburger (Bujarrabal et al. 2009; Teague et al. 2020) may suggest a plateauing of z/r values in the largest of disks. Observations of more disks with extended CO gas disks (>800 au) are required to discern if—and at what gas disk size—line emitting heights begin to plateau.

We also calculated the characteristic z/r for each of the ¹³CO surfaces in our sample, as well as those in the literature from the MAPS disks (Law et al. 2021a) and the Elias 2–27 disk (Paneque-Carreño et al. 2021), which are listed in Table 8. For Elias 2–27, we only considered the uncontaminated western side of the disk to derive z/r and the radial size of the ¹³CO emission. Similar to the CO J = 2–1 surface in the Sz 91 disk, the majority of ¹³CO surfaces were either fit with a single power law or lacked a clear turnover in their surface (making previous fits of r_taper uncertain). Thus, for most disks we manually fixed the cutoff radii. In those cases where no plateauing of the surface was evident, we simply adopted the maximum radius for which a ¹³CO line emitting height was derived as r_cutoff and averaged over the entire radial range of the surfaces.

Figure 11 shows a potential trend between the characteristic z/r of the ¹³CO line emission surfaces with the stellar mass and mean ¹³CO gas temperatures, respectively. In fact, these trends are considerably steeper than those predicted by simple scaling laws. To illustrate this, we overlay the approximate best-fit power laws to each trend in blue in Figure 11 and find that the z/r of the ¹³CO surfaces traces $\sim {M}_{* }^{-3/4}$ and $\sim {T}_{{\rm{B}}}^{-1}$. We also find a tentative link between the size of the ¹³CO gas disk and the z/r emitting height of ¹³CO, which is approximately 4 times steeper in slope than that seen in CO. It is not immediately clear why the ¹³CO trends are so much steeper those of CO, but it may be related to the fact that these disks are not, in fact, vertically isothermal.

Given the small sample size and difficulty in measuring the lower heights associated with ¹³CO, it is possible that the apparent steepness of the ¹³CO trends relative to CO is, at least in part, due to a biased disk sample. To test this idea, we removed those sources without corresponding ¹³CO height measurements and reexamined the CO trends. We found that the mass and gas disk size trends are steeper than the full sample but are still not as steep as seen in ¹³CO, while the gas temperature trend showed no additional steepening. While this may suggest the presence of a possible bias in the ¹³CO disk sample, there still remains tentative evidence that the observed ¹³CO trends are tighter than those in CO. Firm conclusions are, however, limited by our small sample size, and additional high-spatial-resolution ¹³CO line observations of disks are needed. Moreover, higher-angular-resolution observations of the HD 34282 disk would be useful in determining if the ¹³CO surface in this source is also an outlier as in CO.