Temperature Structures of Embedded Disks: Young Disks in Taurus Are Warm

Figure 1 shows the 1.3 mm continuum images and integrated intensity (zeroth moment) maps for C¹⁷O 2 − 1 and H₂CO 3_1,2 − 2_1,1 toward the five sources in our sample. The molecular emission toward IRAS 04302 is highlighted at slightly higher spatial resolution in Figure 2. Radial cuts along the major axis are presented in Figure 3. The continuum emission is elongated perpendicular to the outflow direction for all sources, consistent with a disk as observed before. For the first time, TMC1 is resolved into a close binary (∼85 au separation). We will refer to the two sources as TMC1-E (east) and TMC1-W (west).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Both C¹⁷O and H₂CO are clearly detected toward all sources with a velocity gradient along the continuum structures (see Figure A1). The velocity gradient suggests that the material in TMC1 is located in a circumbinary disk, but a detailed analysis is beyond the scope of this paper. For both molecules, integrated fluxes are similar (within a factor of 2–3) in all sources (Table 2), and both lines have a comparable (factor of 2–4) strength toward each source, with H₂CO brighter than C¹⁷O, except for TMC1A. The H₂CO emission is generally more extended than the C¹⁷O emission, both radially and vertically, except toward TMC1 and TMC1A, where both molecules have the same spatial extent. This is not a signal-to-noise issue, as can be seen from the radial cuts along the major axis (Figure 3).

The most striking feature in the integrated intensity maps is the V-shaped emission pattern of the H₂CO in the edge-on disk IRAS 04302 (see Figure 2), suggesting that the emission arises from the disk surface layers and not the midplane, in contrast to the C¹⁷O emission. The H₂CO emission displays a ringlike structure toward L1527. Given that this disk is also viewed edge-on, this can be explained by emission originating in the disk surface layers, with the outer component along the midplane arising from the envelope. As we will show later in this section, the emission toward IRAS 04302 shows very little envelope contribution, which can explain the difference in morphology between these two sources. The C¹⁷O emission peaks slightly offset (∼60 au) from the L1527 continuum peak, probably due to the dust becoming optically thick in the inner ∼10 au, as seen before for ¹³CO and C¹⁸O (van 't Hoff et al. 2018a). The current resolution does not resolve the inner 10 au; hence, the reduction in CO emission is more extended. In IRAS 04302, a similar offset of ∼60 au is found for both C¹⁷O and H₂CO, suggesting that there may be an unresolved optically thick dust component as well.

Toward L1489, C¹⁷O has a bright inner component (∼200 au) and a weaker outer component that extends roughly as far as the H₂CO emission (∼600 au). A similar structure was observed in C¹⁸O by Sai et al. (2020). The slight rise seen in C¹⁸O emission around ∼300 au to the southwest of the continuum peak is also visible in the C¹⁷O radial cut. Imaging the C¹⁷O data at lower resolution makes this feature clearer in the integrated intensity map. In contrast, the H₂CO emission decreases in the inner ∼75 au, but beyond that, it extends smoothly out to ∼600 au. The off-axis protrusions at the outer edge of the disk pointing to the northeast and southwest were also observed in C¹⁸O and explained as streams of infalling material (Yen et al. 2014).

The C¹⁷O emission peaks slightly (∼40–50 au) off source toward TMC1A. Harsono et al. (2018) showed that ¹³CO and C¹⁸O emission is absent in the inner ∼15 au due to the dust being optically thick. The resolution of the C¹⁷O observations is not high enough to resolve this region, resulting in only a central decrease in emission instead of a gap. A clear gap is visible for H₂CO with the emission peaking ∼100–115 au off source. The central absorption falling below zero is an effect of resolved-out large-scale emission.

Finally, toward TMC1, H₂CO shows a dip at both continuum peaks, while the C¹⁷O emission is not affected by the eastern continuum peak. As discussed for the other sources, this may be the result of optically thick dust in the inner disk. The protrusions seen on the west side in both C¹⁷O and H₂CO are part of a larger arc-like structure that extends toward the southwest beyond the scale shown in the image.

While it is tempting to ascribe all of the compact emission to the young disk, some of it may also come from the envelope and obscure the disk emission. To get a first impression as to whether the observed emission originates in the disk or envelope, position–velocity (pv) diagrams are constructed along the disk major axis for the four single sources (Figure 4). In these diagrams, disk emission is located at small angular offsets and high velocities, while envelope emission extends to larger offsets but has lower velocities. In all sources, C¹⁷O predominantly traces the disk, with some envelope contribution, especially in L1527 and L1489. The H₂CO emission also originates in the disk but has a larger envelope component. An exception is IRAS 04302, which shows hardly any envelope contribution. These results for L1527 are in agreement with previous observations (Sakai et al. 2014a). In L1489, a bright linear feature is present for H₂CO extending from a velocity and angular offset of −2 km s⁻¹ and −2'', respectively, to offsets of 2 km s⁻¹ and 2''. This feature matches the shape of the SO pv diagram (Yen et al. 2014), which was interpreted by the authors as a ring between ∼250 and 390 au. While a brightness enhancement was also identified by Yen et al. (2014) in the C¹⁸O emission (similar to that seen here for H₂CO), the C¹⁷O emission does not display such a feature.

Zoom In Zoom Out Reset image size

Another way to determine the envelope contribution is from the visibility amplitudes. Although a quantitative limit on the envelope contribution to the line emission requires detailed modeling for the individual sources, which will be done in a subsequent paper, a first assessment can be made with more generic models containing either only a Keplerian disk or a disk embedded in an envelope (see Appendix B). For IRAS 04302, both the C¹⁷O and H₂CO visibility amplitude profiles can be reproduced without an envelope. This suggests that there is very little envelope contribution for this source, consistent with the pv diagrams. A disk is also sufficient to reproduce the visibility amplitudes at velocities $\gt | 1|$ km s⁻¹ from the systemic velocity toward L1489, L1527, and TMC1A. For the low velocities, a small envelope contribution is required. The line emission presented here is thus dominated by the disk.

Although both the C¹⁷O and H₂CO emission originates predominantly from the disk, the C¹⁷O emission extends to higher velocities than the H₂CO emission in IRAS 04302, L1527, and TMC1A. This is more easily visualized in the spectra presented in Figure A2. These spectra are extracted in a 6'' circular aperture and only include pixels with >3σ emission. While H₂CO is brighter at intermediate velocities than C¹⁷O (even when correcting for differences in emitting area), it is not present at the highest velocities. Thus, H₂CO emission seems to be absent in the inner disk in these sources, which for TMC1A is also visible in the moment zero map (Figure 1). However, in L1489, both molecules have similar maximum velocities. Toward TMC1, they extend to the same redshifted velocity, while C¹⁷O emission is strongly decreased at blueshifted velocities as compared to the redshifted velocities.

To compare the C¹⁷O and H₂CO observations between the different sources more quantitatively, we calculate disk-averaged total column densities, N_T, assuming optically thin emission in local thermodynamic equilibrium (LTE) using

$$\begin{eqnarray}&&\displaystyle \frac{4\pi F{\rm{\Delta }}v}{{A}_{{ul}}{\rm{\Omega }}{{hcg}}_{\mathrm{up}}}=\displaystyle \frac{{N}_{T}}{Q({T}_{\mathrm{rot}})}{e}^{-{E}_{\mathrm{up}}/{{kT}}_{\mathrm{rot}}},\end{eqnarray} \tag{ 1 }$$

where FΔv is the integrated flux density; A_ul is the Einstein A coefficient; Ω is the solid angle subtended by the source; E_up and g_up are the upper-level energy and degeneracy, respectively; and T_rot is the rotational temperature.

The integrated fluxes are measured over the dust-emitting area (Table 3). We note that this does not necessarily encompass the total line flux, but it will allow for an abundance estimate as described below. A temperature of 30 K is adopted for C¹⁷O and 100 K for H₂CO, as these are slightly above their freeze-out temperatures. The C¹⁷O column density ranges between ∼2 and 20 × 10¹⁵ cm⁻², with the lowest value toward L1489 and the highest value toward TMC1A. The H₂CO column density is about an order of magnitude lower, with values between ∼4 and 18 × 10¹⁴ cm⁻². The lowest value is found toward TMC1A, and the highest value is found toward L1527. Changing the temperature for H₂CO to 30 K decreases the column densities by only a factor of ≲3.

The H₂CO column density toward L1527 is a factor of 3–6 higher than previously derived by Sakai et al. (2014a), possibly because they integrated over different areas and velocity ranges for the envelope, disk, and envelope–disk interface. Integrating the H₂CO emission over a circular aperture of 05 and excluding the central $| {\rm{\Delta }}v| \leqslant 1.0$ km s⁻¹ channels to limit the contribution from the envelope and resolved-out emission results in an H₂CO column density of 9.7 × 10¹³ cm⁻², only a factor of 2–3 higher than that found by Sakai et al. (2014a). Pegues et al. (2020) found H₂CO column densities spanning 3 orders of magnitude (∼5 × 10¹¹–5 × 10¹⁴ cm⁻²) for a sample of 13 Class II disks. The values derived here for Class I disks are thus similar to the high end (≲4 times higher) of the values for Class II disks.

An assessment of the molecular abundances can be made by estimating the H₂ column density from the continuum flux. First, we calculate the disk dust masses, M_dust, from the integrated continuum fluxes, F_ν, using

$$\begin{eqnarray}&&{M}_{\mathrm{dust}}=\displaystyle \frac{{D}^{2}{F}_{\nu }}{{\kappa }_{\nu }{B}_{\nu }({T}_{\mathrm{dust}})},\end{eqnarray} \tag{ 2 }$$

where D is the distance to the source, κ_ν is the dust opacity with the assumption of optically thin emission, and B_ν is the Planck function for a temperature T_dust (Hildebrand 1983). Adopting a dust opacity of ${\kappa }_{1.3\mathrm{mm}}=2.25$ cm² g⁻¹, as used for Class II disks by, e.g., Ansdell et al. (2016), and a dust temperature of 30 K similar to, e.g., Tobin et al. (2015) for embedded disks results in disk dust masses between 3.7 M_E for TMC1-E and 75 M_E for TMC1A. Using the same dust opacity as for Class II disks is probably reasonable if grain growth starts early on in the disk formation process. However, adopting ${\kappa }_{1.3\mathrm{mm}}=0.899$ cm² g⁻¹, as is often done for protostellar disks and envelopes (e.g., Jørgensen et al. 2007; Andersen et al. 2019; Tobin et al. 2020), only affects the molecular abundances by a factor of ∼2. Assuming a gas-to-dust ratio of 100 and using the size of the emitting region, these dust masses result in H₂ column densities of 2–90 × 10²³ cm⁻².

The resulting C¹⁷O and H₂CO abundances are listed in Table 3. For C¹⁷O, the abundances range between 1.2 × 10⁻⁸ and 1.2 × 10⁻⁷. Assuming a C¹⁶O/C¹⁷O ratio of 1792 (as in the ISM; Wilson & Rood 1994), a CO ISM abundance of 10⁻⁴ with respect to H₂ corresponds to a C¹⁷O abundance of 5.6 × 10⁻⁸. The derived C¹⁷O abundances are thus within a factor of 5 of the ISM abundance, suggesting that no substantial processing has happened, as observed for Class II disks where the CO abundance can be 2 orders of magnitude below the ISM value (e.g., Favre et al. 2013). These results are consistent with the results from Zhang et al. (2020) for three Class I disks in Taurus (including TMC1A) but not with the order-of-magnitude depletion found by Bergner et al. (2020) for two Class I disks in Serpens. For H₂CO, the abundance ranges between ∼3 × 10⁻¹⁰ and ∼8 × 10⁻⁹ in the different sources, except for TMC1A, where the abundance is ∼5 × 10⁻¹¹, probably due to the absence of emission in the inner region. Abundances around 10⁻¹⁰–10⁻⁹ are consistent with chemical models for protoplanetary disks (e.g., Willacy & Woods 2009; Walsh et al. 2014). However, H₂CO abundances derived for TW Hya and HD 163296 are 2–3 orders of magnitude lower, 8.9 × 10⁻¹³ and 6.3 × 10⁻¹², respectively (Carney et al. 2019).

A caveat in determining these abundances is the assumption that the continuum and line emission are optically thin. As discussed in Section 3.1, there is likely an optically thick dust component that would result in underestimates of the dust masses and overestimates of the abundances. On the other hand, optically thick dust hides molecular line emission originating below its τ = 1 surface, which leads to underestimates of the abundances. Based on the results from Zhang et al. (2020), C¹⁷O may be optically thick in Class I disks. This would also result in underestimating the abundances. Scaling the dust temperature used in Equation (2) with luminosity, as done by Tobin et al. (2020) for embedded disks in Orion, results in dust masses lower by a factor of ∼2 and therefore slightly higher abundances. Moreover, the integrated line flux is assumed to originate solely in the disk, but as shown in Figure 4, there can be envelope emission present. Finally, the H₂CO emission originates in the disk surface layers, which means the abundances are higher than derived here assuming emission originating throughout the disk. To take all these effects into account, source-specific models are required.

Water and methanol form on ice-covered dust grains and thermally desorb into the gas phase at temperatures ∼100–150 K. These molecules are thus expected to trace the hot region inside the water snowline. The observations cover one HDO (deuterated water) transition (3_1,2 − 2_2,1) with an upper-level energy of 168 K and 16 transitions in the CH₃OH J = 5_k − 4_k branch with upper-level energies ranging between 34 and 131 K. None of these lines are detected in any of the disks.

To compare these nondetections to observations in other systems, a 3σ upper limit is calculated for the disk-averaged total column density by substituting

$$\begin{eqnarray}&&3\sigma =3\times 1.1\sqrt{\delta v{\rm{\Delta }}V}\times \mathrm{rms}\end{eqnarray} \tag{ 3 }$$

for the integrated flux density, FΔv, in Equation (1). Here δv is the velocity resolution, and ΔV is the line width expected based on other line detections. The factor of 1.1 takes a 10% calibration uncertainty into account. Assuming the water and methanol emission arises from the innermost part of the disk, the rms is calculated from the baseline of the spectrum integrated over a central 05 diameter aperture (∼one beam) and amounts to ∼2.7 mJy for HDO and ∼3.0 mJy for CH₃OH. A line width of 4 km s⁻¹ and a rotational temperature of 100 K are adopted.

A 3σ column density upper limit of ∼8 × 10¹³ cm⁻² is then found for HDO. This is 1–2 orders of magnitude below the column densities derived for the Class 0 sources NGC1333 IRAS2A, NGC1333 IRAS4A-NW, and NGC1333 IRAS4B (∼10¹⁵–10¹⁶ cm⁻²; Persson et al. 2014) and more than 3 orders of magnitude lower than toward the Class 0 source IRAS 16293A (∼5 × 10¹⁷ cm⁻²; Persson et al. 2013). Taking into account the larger beam size of the earlier observations (∼1'') lowers the column density derived here by only a factor of ∼4. Furthermore, Taquet et al. (2013) showed that the HDO observations toward NGC1333 IRAS2A and NGC1333 IRAS4A are consistent with column densities up to 10¹⁹ and 10¹⁸ cm⁻², respectively, using a grid of non-LTE large velocity gradient (LVG) radiative transfer models.

For CH₃OH, the 5_0,5 − 4_0,4 (A) transition provides the most stringent upper limit of ∼8 × 10¹⁴ cm⁻². This upper limit is orders of magnitude lower than the column density toward the Class 0 source IRAS 16293 (2 × 10¹⁹ cm⁻² within a 70 au beam; Jørgensen et al. 2016) and the young disk around the outbursting star V883 Ori (disk-averaged column density of ∼1.0 × 10¹⁷ cm⁻²; van 't Hoff et al. 2018b). A similarly low upper limit (5 × 10¹⁴ cm⁻²) was found for a sample of 12 Class I disks in Ophiuchus (Artur de la Villarmois et al. 2019). However, this upper limit is not stringent enough to constrain the column down to the value observed in the TW Hya protoplanetary disk (peak column density of 3–6 × 10¹² cm⁻²; Walsh et al. 2016) or the upper limit in the Herbig Ae disk HD 163296 (disk-averaged upper limit of 5 × 10¹¹ cm⁻²; Carney et al. 2019).

For a better comparison with other sources, column density ratios are calculated with respect to H₂ and H₂CO and reported in Table 3. Using the H₂ column density derived from the continuum flux, upper limits of ∼1–40 × 10⁻¹⁰ are found for the HDO abundance. The CH₃OH upper limits range between 1 and 40 × 10⁻⁹. This is orders of magnitude lower than what is expected from ice observations (10⁻⁶ − 10⁻⁵; Boogert et al. 2015), and thus from thermal desorption, as observed in IRAS 16293 (≲3 × 10⁻⁶; Jørgensen et al. 2016) and V883 Ori (∼4 × 10⁻⁷; van 't Hoff et al. 2018b). Abundances for nonthermally desorbed CH₃OH in TW Hya are estimated to be ∼10⁻¹²−10⁻¹¹ (Walsh et al. 2016). Sakai et al. (2014a) detected faint CH₃OH emission (from different transitions than targeted here) toward L1527, with a CH₃OH/H₂CO ratio between 0.6 and 5.1. Our upper limit of 4.6 for L1527 is consistent with these values. CH₃OH/H₂CO ratios of 1.3 and <0.2 were derived for TW Hya and HD 163296, respectively, but our CH₃OH upper limit is not stringent enough to make a meaningful comparison. An assumption here is that the emitting regions of CH₃OH and H₂CO are cospatial. As noted in Section 3.1, H₂CO seems absent in the inner disk where CH₂OH is expected.

For (near) edge-on disks, CO freeze-out should be readily observable, as CO emission will be missing from the outer disk midplane (Dutrey et al. 2017; van 't Hoff et al. 2018a). Van 't Hoff et al. (2018a) studied the effect of CO freeze-out on the optically thick ¹³CO and C¹⁸O emission in L1527. The less abundant C¹⁷O is expected to be optically thin and mainly traces the disk. Here we employ the models from van 't Hoff et al. (2018a) to predict the C¹⁷O emission pattern for varying degrees of CO freeze-out (see Figure C1): a "warm" model (no CO freeze-out), an "intermediate" model (CO freeze-out in the outer disk midplane), and a "cold" model (CO freeze-out in most of the disk, except the inner part and surface layers). Briefly, in these models, gaseous CO is present at a constant abundance of 10⁻⁴ with respect to H₂ in the regions in the disk where T > 20 K and in the envelope. For the warm model, the L1527 temperature structure from Tobin et al. (2013) is adopted, and for the intermediate and cold models, the temperature is reduced by 40% and 60%, respectively. There is no CO freeze-out in the 125 au disk in the warm model, while the intermediate and cold models have the CO snowline at 71 and 23 au, respectively. Synthetic image cubes are generated using the radiative transfer code LIME (Brinch & Hogerheijde 2010), making use of the C¹⁷O LAMDA file (Schöier et al. 2005) for the LTE calculation, and are convolved with the observed beam size.

Figure 5 shows moment zero maps integrated over the low, intermediate, and high velocities for the warm and cold edge-on disk model. Models with and without an envelope are presented. The difference between the warm and cold model is most clearly distinguishable at intermediate velocities (Figure 5, middle row). In the absence of an envelope, the emission becomes V-shaped in the cold model, tracing the warm surface layers where CO is not frozen out. This V shape is not visible when there is a significant envelope contribution. The cold model differs from the warm model in that the envelope emission becomes comparable in strength to the disk emission when CO is frozen out in most of the disk. In the warm case, the disk emission dominates over the envelope emission. At low velocities (Figure 5, top row), the difference between a warm and cold disk can be distinguished as well when an envelope is present, although in practice, this will be much harder due to resolved-out emission at these central velocities. Without an envelope, the low-velocity emission originates near the source center due to the rotation, and the models are indistinguishable, except for differences in the flux. Due to the rotation, the emission at these velocities gets projected along the minor axis of the disk (that is, east–west). At the highest velocities (Figure 5, top row), the emission originates in the inner disk, north and south of the source. If CO is absent in the midplane, very high angular resolution is required to observe this directly through a V-shaped pattern.

Zoom In Zoom Out Reset image size

The C¹⁷O moment zero maps integrated over different velocity intervals for IRAS 04302 and L1527 are presented in Figure 5. The observations show no sign of CO freeze-out in L1527 and resemble the warm model (most clearly seen at intermediate velocities), consistent with previous results for C¹⁸O and ¹³CO (van 't Hoff et al. 2018a). On the other hand, IRAS 04302 displays a distinct V-shaped pattern at intermediate velocities, suggesting that CO is frozen out in the outer part of this much larger disk (∼250 au, compared to 75–125 au for L1527; Tobin et al. 2013; Aso et al. 2017; Sheehan & Eisner 2017).

The vertical distribution of the emission in both disks is highlighted in Figure 6 with vertical cuts at different radii. In L1527, the C¹⁷O emission peaks at the midplane throughout the disk, while for IRAS 04302, the peaks shift to layers higher up in the disk for radii ≳110 au. A first estimate of the CO snowline location can be made based on the location of the V shape. In the cold model, the CO snowline is located at 23 au, but due to the size of the beam, the base of the V shape and the first occurrence of a double peak in the vertical cuts are at ∼55 au. In IRAS 04302, the V shape begins at a radius of ∼130 au, so the CO snowline location is then estimated to be around ∼100 au.

Zoom In Zoom Out Reset image size

A clear V-shaped pattern is also visible in the H₂CO integrated emission map for IRAS 04302 (Figure 1). The V shape starts at around 55 au (∼1 beam offset from the continuum peak). If the reduction of H₂CO in the midplane is fully due to freeze-out, the snowline is then located around (or inward of) ∼25 au. In L1527, H₂CO emission also appears to come from surface layers, except in the outer disk (see Figures 1 and 6). The cold models show that CO emission from the envelope becomes comparable in strength to emission from the disk if CO is frozen out in a large part of the disk. Given that the envelope contribution is much larger in L1527 than in IRAS 04302, the emission peaking in the outer disk midplane is likely originating in the envelope. Instead of a clear V shape, the emission in the inner region forms two bright lanes along the continuum position. A similar pattern is seen in the individual channels. This suggests that the H₂CO snowline is unresolved at the current resolution and closer in than in IRAS 04302 (≲25 au).

A zeroth-order estimate of the midplane temperature profile for IRAS 04302 can be made from these two snowline estimates using a radial power law, T ∝ R^−q. For disks, often a power-law exponent q of 0.5 is assumed, but q can range between 0.33 and 0.75 (see, e.g., Adams & Shu 1986; Kenyon et al. 1993; Chiang & Goldreich 1997). A power law with q = 0.75 matches the two temperature estimates reasonably well (see Figure 7). This temperature profile is quite similar to the profile constructed for L1527 based on ¹³CO and C¹⁸O temperature measurements (van 't Hoff et al. 2018a). The L1527 temperature profile predicts an H₂CO snowline radius of ≲10 au, consistent with the results derived above. Thus, IRAS 04302 is warm like L1527, with freeze-out occurring only in the outermost part of this large disk.

Zoom In Zoom Out Reset image size

For less inclined disks, observing freeze-out directly is much harder; the projected area between the top and bottom layer becomes smaller (that is, the V shape becomes more narrow), therefore requiring higher spatial resolution to observe it. In addition, because now both the near and the far sides of the disk become visible, emission from the far side's surface layers can appear to come from the near side's midplane (see Figure C2 and Pinte et al. 2018), which makes a V shape due to emission originating only in the surface layers that are harder to observe. For the L1527 disk model, the intermediate and warm models become quite similar for an inclination of 60° at this angular resolution, and only a cold disk shows a clear V-shaped pattern (Figure 8).

Zoom In Zoom Out Reset image size

Figure 9 shows the C¹⁷O moment zero maps for the intermediate inclined disks TMC1A and L1489. The disk size, stellar mass, and stellar luminosity of TMC1A are comparable to L1527. At intermediate velocities, there is no sign of a V-shaped pattern, so these observations do not suggest substantial freeze-out of CO in TMC1A. In order to constrain the CO snowline a little better, models were run with snowline locations of 31, 42, and 56 au (that is, in between the cold and intermediate models). All three models show a V shape, suggesting that the CO snowline is at radii ≳70 au in TMC1A. This is consistent with the results from Aso et al. (2015), who found a temperature of 38 K at 100 au from fitting a disk model to ALMA C¹⁸O observations, and Harsono et al. (submitted), who found a temperature of 20 K at 115 au. There is no sign of a V-shaped pattern in the H₂CO emission.

Zoom In Zoom Out Reset image size

For L1489, the intermediate velocities show a more complex pattern, with CO peaking close to the source and at larger offsets (≳2''). A similar structure was seen in C¹⁸O (Sai et al. 2020). This could be the result of nonthermal desorption of CO ice in the outer disk if the dust column is low enough for UV photons to penetrate (Cleeves 2016) or due to a radial temperature inversion resulting from radial drift and dust settling (Facchini et al. 2017). Such a double CO snowline has been observed for the protoplanetary disk IM Lup (Öberg et al. 2015; Cleeves 2016). The structure of the continuum emission, a bright central part and a fainter outer part, makes these plausible ideas. Another possibility is that the extended emission is due to a warm inner envelope component. The UV irradiated mass of L1489 derived from ¹³CO 6–5 emission is similar to that of L1527 and higher than for TMC1A and TMC1 (Yıldız et al. 2015). This may provide a sufficient column along the outflow cavity wall for C¹⁷O emission to be observed. A high level of UV radiation is supported by O and H₂O line fluxes (Karska et al. 2018).

If the edge of the compact CO emission is due to freeze-out, the CO snowline is located at roughly 200 au. Models based on the continuum emission have temperatures of ∼30 or ∼20–30 K at 200 au (Brinch et al. 2007 and Sai et al. 2020, respectively), so CO could indeed be frozen out in this region. The H₂CO emission does not show a gap at 200 au, which could mean that the emission is coming from the surface layers. The C¹⁷O (and C¹⁸O) abundance in these warmer surface layers may then be too low to be detected at the sensitivity of these observations.

We have used observations of C¹⁷O and H₂CO toward five class I disks in Taurus to address whether embedded disks are warmer than more evolved Class II disks. While the C¹⁷O observations can indicate the presence or absence of ≲20 K gas, the addition of H₂CO observations allows one to further constrain the temperature profile. The picture that is emerging suggests that these young disks have midplanes with temperatures between ∼20 and ∼70 K: cold enough for H₂CO to freeze out but warm enough to retain CO in the gas phase (Figure 10). This suggests that, for example, the elemental C/O ratio in both the gas and ice could be different from that in protoplanetary disks. If planet formation starts during the embedded phase, the conditions for the first steps of grain growth are then different than generally assumed.

Zoom In Zoom Out Reset image size

Young disks being warmer than protoplanetary disks can also have consequences for the derived disk masses from continuum fluxes. This has been taken into consideration in recent literature by adopting a dust temperature of 30 K for solar-luminosity protostars (Tobin et al. 2015, 2016b; Tychoniec et al. 2018; Tychoniec et al. 2020), although not uniformly (e.g., Andersen et al. 2019; Williams et al. 2019), while 20 K is generally assumed for protoplanetary disks (e.g., Ansdell et al. 2016). In their study of Orion protostars, Tobin et al. (2020) took this one step further by scaling the temperature by luminosity based on a grid of radiative transfer models resulting in an average temperature of 43 K for a 1 L_⊙ protostar. Since higher temperatures will result in lower masses for a certain continuum flux, detailed knowledge of the average disk temperature is crucial to determine the mass reservoir available for planet formation. While the current study shows that embedded disks are warmer than protoplanetary disks, and the radial temperature profiles for L1527 and IRAS 04302 hint that 30 K may be to low for the average disk temperature, source-specific modeling of the continuum and molecular line emission is required to address what would be an appropriate temperature to adopt for the mass derivation. However, an increase in temperature by a factor of 2 will lower the mass by only a factor  of 2 (see Equation (2)), and Tobin et al. (2020) still found embedded disks to be more massive than protoplanetary disks by a factor >4. Differences in temperature can thus not account for the mass difference observed between embedded and protoplanetary disks.

5.1.1. The Textbook Example of IRAS 04302

5.1.2. Comparison with Protostellar Envelopes and Protoplanetary Disks

No sign of CO freeze-out is detected in the C¹⁷O observations of L1527, and while freeze-out is much more difficult to see in non-edge-on disks, TMC1A does not show hints of freeze-out at radii smaller than ∼70 au. A first estimate puts the CO snowline at ∼100 au in IRAS 04302, and the CO snowline may be located around ∼200 au in L1489. These young disks are thus warmer than T Tauri disks, where the snowline is typically at a few tens of au, as can be seen in Figure 11. We only include class II disks for which a CO snowline location has been reported based on molecular line observations, either ¹³C¹⁸O (for TW Hya; Zhang et al. 2017) or N₂H⁺ (Qi et al. 2019). There is no clear trend between CO snowline location and bolometric luminosity for either Class, but the Class I disks have CO snow lines at larger radii compared to Class II disks with similar bolometric luminosities.

Zoom In Zoom Out Reset image size

In protostellar envelopes, snowline radii larger than expected based on the luminosity have been interpreted as a sign of a recent accretion burst (Jørgensen et al. 2015; Frimann et al. 2017; Hsieh et al. 2019). During such a time period of increased accretion, the circumstellar material heats up, shifting the snow lines outward. Once the protostar returns to its quiescent stage, the temperature adopts almost instantaneously, while the chemistry takes longer to react. During this phase, the snow lines are at larger radii than expected from the luminosity. The results in Figure 11 could thus indicate that small accretion bursts have occurred in the Class I systems and that the CO snow lines have not yet shifted back to their quiescent location. When such a burst should have happened depends on the freeze-out timescale, τ_fr,

$$\begin{eqnarray}&&{\tau }_{\mathrm{fr}}=1\times {10}^{4}\,\mathrm{yr}\sqrt{\displaystyle \frac{10\,{\rm{K}}}{{T}_{\mathrm{fr}}}}\displaystyle \frac{{10}^{6}\,{\mathrm{cm}}^{-3}}{{n}_{{{\rm{H}}}_{2}}},\end{eqnarray} \tag{ 4 }$$

where T_fr is the freeze-out temperature and ${n}_{{{\rm{H}}}_{2}}$ is the gas density (Visser et al. 2012). For densities ≳10⁸ cm⁻³, the CO freeze-out timescale is ≲100 yr. This could suggest that Class I protostars frequently undergo small accretion bursts. Alternatively, these young disks may have lower densities than more evolved disks. As shown by the model results from N. M. Murillo et al. (2020, in preparation), decreasing the density while keeping the luminosity constant shifts the snow lines outward. If this is what is causing the results in Figure 11, this means that embedded disks not only have different temperature structures from protoplanetary disks but also have different density structures. However, the larger disk masses derived for embedded disks compared to protoplanetary disks for similar disk radii make this unlikely (Tobin et al. 2020).

Another comparison is made in Figure 7, where the radial temperature profiles inferred for L1527 and IRAS 04302 are shown together with those for the younger Class 0 disklike structure around IRAS 16293A (van 't Hoff et al. 2020) and the Class II disk TW Hya (Schwarz et al. 2016). The young disks are warmer than the more evolved Class II disk but much colder than the Class 0 system IRAS 16293A. When making this comparison, one should keep in mind that IRAS 16293A reflects an envelope where the temperature will be larger at larger scales because of the spherical rather than disk structure. In a disk, the temperature will drop more rapidly in the radial direction due to the higher extinction compared to an envelope. Nevertheless, such an evolutionary trend is expected because the accretion rate decreases as the envelope and disk dissipate. As a consequence, heating due to viscous accretion diminishes, and hence the temperature drops, as shown by two-dimensional physical and radiative transfer models for embedded protostars (D'Alessio et al. 1997; Harsono et al. 2015). In addition, the blanketing effect of the envelope decreases as the envelope dissipates (Whitney et al. 2003).

As a first comparison between the observations and model predictions, models from Harsono et al. (2015) are overlaid on the observationally inferred temperature profiles in Figure 7 (right panel). In these models, the dust temperature is determined based on stellar irradiation and viscous accretion. Models are shown for a stellar luminosity of 1 L_⊙, an envelope mass of 1 M_⊙, a disk mass of 0.05 M_⊙, a disk radius of 200 au, and different accretion rates. The disk mass has a negligible effect on the temperature profiles (see Harsono et al. 2015 for details). The observations for IRAS 16293A match reasonably well with the temperature profile for a heavily accreting system (10⁻⁴ M_⊙ yr⁻¹), consistent with estimates of the accretion rate (e.g., ∼5 × 10⁻⁵ M_⊙ yr⁻¹; Schöier et al. 2002). However, because in these models, the total luminosity is based on the stellar and accretion luminosity (and a contribution from the disk), the match for IRAS 16239A with a strong accretion model may just reflect the system's bolometric luminosity of 20 L_⊙. In contrast, the temperature profiles for L1527 and IRAS 04302 are comparable to the colder 10⁻⁷ M_⊙ yr⁻¹ model, consistent with the accretion rate of ∼3 × 10⁻⁷ M_⊙ yr⁻¹ for L1527 (see van 't Hoff et al. 2018a). Similar accretion rates on the order of 10⁻⁷ M_⊙ yr⁻¹ have been reported for L1489, TMC1A, and TMC1 (e.g., Mottram et al. 2017; Yen et al. 2017) based on the bolometric luminosities (see, e.g., Stahler et al. 1980; Palla & Stahler 1993). We are not aware of a measurement toward IRAS 04302, but our very preliminary modeling results (M. L. R. van 't Hoff et al. 2020, in preparation) are consistent with an accretion rate on the order of 10⁻⁷ M_⊙ yr⁻¹. Measured accretion rates for TW Hya range between ∼2 × 10⁻¹⁰ and 2 × 10⁻⁹ M_⊙ yr⁻¹ (e.g., Herczeg & Hillenbrand 2008; Curran et al. 2011; Ingleby et al. 2013), and accretion rates of ∼10⁻¹⁰−10⁻⁸ M_⊙ yr⁻¹ are typically measured for protoplanetary disks around T Tauri stars (see Hartmann et al. 2016 for a review).

The results presented here thus provide observational evidence for cooling of the circumstellar material during evolution. More sources need to be observed to confirm this trend and answer more detailed questions, such as, when has a disk cooled down sufficiently for large-scale CO freeze-out? Does this already happen before the envelope dissipates? The object IRAS 04302 is a borderline Class I/Class II object embedded in the last remnants of its envelope, but it still has a temperature profile more similar to L1527 than TW Hya. Although a caveat here may be the old age of TW Hya (∼10 Myr), this hints that disks may stay warm until the envelope has fully dissipated.

5.1.3. TMC1

One of the major questions regarding the chemical composition of planetary material is whether it contains complex organic molecules (COMs). Due to the low temperatures in protoplanetary disks, observations of COMs are very challenging because these molecules thermally desorb at temperatures ≳100–150 K, that is, in the inner few au. In contrast, COMs are readily detected on disk scales in protostellar envelopes (e.g., IRAS 16293, NGC1333 IRAS2A, NGC1333 IRAS4A, and B1-c; Taquet et al. 2015; Jørgensen et al. 2016; van Gelder et al. 2020) and in the young disk V883 Ori, where a luminosity outburst has heated the disk and liberated the COMs from the ice mantles (van 't Hoff et al. 2018b; Lee et al. 2019).

Although young disks seem warmer than protoplanetary disks, the CH₃OH and HDO nondetections with upper limits orders of magnitude below the column densities observed toward Class 0 protostellar envelopes suggest that they are not warm enough to have a hot core–like region with a large gas reservoir of COMs. This is consistent with recent findings by Artur de la Villarmois et al. (2019) for a sample of Class I protostars in Ophiuchus. More stringent upper limits are required for comparison with the Class II disks TW Hya and HD 163296. However, the detection of HDO and CH₃OH may have been hindered by optically thick dust in the inner region or the high inclinations of these sources. Modeling by N. M. Murillo et al. (2020, in preparation) shows that the water snowline is very hard to detect in near-edge-on disks. These nondetections thus do not rule out the presence of HDO and CH₃OH; in fact, if the region where HDO and CH₃OH are present is much smaller than the beam, they may have higher columns than the upper limits derived here. This is corroborated by the weak detection of CH₃OH in L1527 (Sakai et al. 2014a). These results thus merely show that Class I disks do not have an extended hot core–like region, making the detection of COMs just as challenging as in Class II disks.

A question related to the chemical composition is whether the disk material is directly inherited from the cloud, processed en route to the disk, or even fully reset upon entering the disk. Young disks like L1527, where no CO freeze-out is observed, suggest that no full inheritance takes place, at least not for the most volatile species like CO. Ice in the outer disk of IRAS 04302 could be inherited. However, the freeze-out timescale for densities >10⁶ cm⁻³ is <10⁴ yr, so this CO could have sublimated upon entering the disk and frozen out as the disk cooled (see, e.g., Visser et al. 2009). Without CO ice, additional grain-surface formation of COMs will be limited in the young disks. So if COMs are present in more evolved disks, as, for example, shown for V883 Ori, they must have been inherited from a colder precollapse phase. Physicochemical models show that prestellar methanol can indeed be incorporated into the disk (Drozdovskaya et al. 2014).

While the H₂CO emission is brighter than the C¹⁷O emission at intermediate velocities, no H₂CO emission is detected at the highest velocities in IRAS 04302, L1527, and TMC1A, suggesting a reduction in H₂CO flux in the inner ≲20–30 au in these disks. This is not just a sensitivity issue, as, for example, C¹⁷O and H₂CO have similar strengths and emitting areas in channels around +1.9 km s⁻¹ with respect to the source velocity in L1527, while 3.05 km s⁻¹ is the highest velocity observed for C¹⁷O and 2.60 km s⁻¹ the highest velocity for H₂CO. The decrease in H₂CO emission is also unlikely to be due to the continuum being optically thick because this would affect the C¹⁷O emission as well, unless there is significantly more C¹⁷O emission coming from layers above the dust millimeter τ = 1 surface than H₂CO emission. Given the observed distributions, with H₂CO being vertically more extended than C¹⁷O, this seems not to be the case. Moreover, the drop in H₂CO in TMC1A occurs much further out than where the dust becomes optically thick.

Formaldehyde rings have also been observed in the protoplanetary disks around TW Hya (Öberg et al. 2017), HD 163296 (Qi et al. 2013a; Carney et al. 2017), DM Tau (Henning & Semenov 2008; Loomis et al. 2015), and DG Tau (Podio et al. 2019). Interestingly, a ring is only observed for the 3₀₃ − 2₀₂ and 3₁₂ − 2₁₁ transitions and not for the 5₁₅ − 4₁₄ transition. Öberg et al. (2017) argued that the dust opacity cannot be the major contributor in TW Hya because the dust opacity should be higher at higher frequencies, thus for the 5₁₅ − 4₁₄ transition. Instead, they suggested a warm inner component that is visible in the 5₁₅ − 4₁₄ transition (E_up = 63 K) and not in the 3₁₂ − 2₁₁ transition (E_up = 33 K). For L1527, we observe the 3₁₂ − 2₁₁ transition, and radiative transfer modeling for the L1527 warm disk model shows that both the C¹⁷O (E_up = 33 K) and H₂CO emission go down by a factor of ∼2 if the temperature is increased by 80%. An excitation effect thus seems unlikely, unless the C¹⁷O emission is optically thick. The latter is not expected, given that the C¹⁸O in L1527 is only marginally optically thick (van 't Hoff et al. 2018a). The absence of H₂CO emission in the inner disk thus points to a reduced H₂CO abundance. A lower total (gas + ice) H₂CO abundance (more than an order of magnitude) in the inner 30 au is seen in models by Visser et al. (2011), who studied the chemical evolution from prestellar core into disk, but these authors did not discuss the H₂CO chemistry.

The H₂CO abundance in the inner disk can be low if its formation is inefficient. It can form in both the gas and ice (e.g., Willacy & Woods 2009; Walsh et al. 2014; Loomis et al. 2015). On the grain surfaces, the dominant formation route is through hydrogenation of CO (Watanabe & Kouchi 2002; Cuppen et al. 2009; Fuchs et al. 2009). Since there seems to be no CO freeze-out in these young disks, or only at radii ≳100 au, H₂CO is expected to form predominantly in the gas. Ring-shaped H₂CO emission due to increased ice formation outside the CO snowline, as used to explain the ring observed in HD 163296 (Qi et al. 2013a), is thus not applicable to the disks in this sample.

In the gas, the reaction between CH₃ and O is the most efficient way to form H₂CO (e.g., Loomis et al. 2015). Therefore, a decrease in gas-phase H₂CO formation would require a low abundance of either CH₃ or O. CH₃ is efficiently produced by photodissociation of CH₄ or through ion–molecule reactions. A low CH₃ abundance thus necessitates the majority of carbon to be present in CO, in combination with a low X-ray flux, as carbon can only be liberated from CO by X-ray-generated He⁺. Atomic oxygen is formed through photodissociation of H₂O and CO₂ or dissociation of CO via X-ray-generated He⁺. A low atomic oxygen abundance would thus require a low UV and X-ray flux.

Besides a low formation rate, a high destruction rate would also decrease the amount of H₂CO. However, the destruction products have a limited chemistry, and re-creation of H₂CO is the most likely outcome. Willacy & Woods (2009) showed that a third of the ions formed by H₂CO destruction through HCO⁺ and DCO⁺ form CO instead of reforming H₂CO, leading to a depletion between 7 and 20 au for their disk model. However, this only reduces H₂CO in the midplane, not in the surface layers. In addition, Henning & Semenov (2008) suggested the conversion of CO into CO₂-containing molecules and hydrocarbons that freeze out onto dust grains (see also Aikawa et al. 1999). However, the C¹⁷O observations do not suggest heavy CO depletion.

Another effect that could contribute is photodesorption of methanol ice that is inherited from earlier phases. Laboratory experiments have shown that methanol does not desorb intact upon vacuum ultraviolet (VUV) irradiation but rather leads to the release of smaller photofragments including H₂CO (Bertin et al. 2016; Cruz-Diaz et al. 2016). This could lead to an increase of H₂CO outside the region where CH₃OH ice thermally desorbs (∼100–150 K). Finally, turbulence may play a role, as models by Furuya & Aikawa (2014) show the formation of H₂CO rings when mixing is included. However, these rings are due to a decrease of H₂CO inside the CO snowline and an increase outside this snowline, and these results may not be applicable to embedded disks without CO freeze-out. Observations of higher-excitation H₂CO lines and chemical modeling with source-specific structures may provide further insights.

It is worth noting that Pegues et al. (2020) found both centrally peaked and centrally depressed H₂CO emission profiles for a sample of 15 protoplanetary disks. A reduction of H₂CO emission toward three out of the five disks in our sample could mean that the H₂CO distribution is set during the embedded stage.