© ESO 2022

1 Introduction

Carbon and oxygen are two of the most abundant elements in interstellar clouds and play a crucial role in the chemistry of star forming regions and planetary atmospheres (Langer 2009; Madhusudhan 2019). If volatile carbon is carried mostly by dust grains, these grains can grow to produce planetesimals rich in complex organic molecule ices, the building blocks of life. CO isotopologue surveys have recently revealed that CO is frequently underabundant by orders of magnitude in the gas phase in the protoplanetary disks of T Tauri stars (van Zadelhoff et al. 2001; Dutrey et al. 2003; Chapillon et al. 2008; Favre et al. 2013; Schwarz et al. 2016; Ansdell et al. 2016, 2018; Long et al. 2017; Miotello et al. 2017). The amount of volatile depletion observed cannot be explained solely by freeze-out, photodissociation, or chemical processing in the warm molecular layer (van Zadelhoff et al. 2001; Williams & Best 2014; Miotello et al. 2017; Schwarz et al. 2018; Bosman et al. 2018; Zhang et al. 2019). Detailed analysis of the TW Hya disk, using [C I], HD, CO, and C₂H emission, shows that the depletion in CO is related to elemental carbon depletion (Bergin et al. 2013; Kama et al. 2016b; Calahan et al. 2021) and cannot be due to high dust-to-gas ratios only.

The outer disk gas-phase carbon depletion can be explained by the freeze-out of CO onto dust grains, coupled with a combination of chemical processing into less volatile carbon-carrying molecules and dust evolution locking the carbon in the ice on large dust grains in the midplane (Kama et al. 2016b; Krijt et al. 2020). The long timescales for carbon depletion processes are expected to result in stronger depletion, by several orders of magnitude, in older Class II disks (Krijt et al. 2020). Efficient vertical gas mixing is necessary to replenish the CO gas available for freeze-out near the disk midplane, which eventually results in lower carbon abundances in the atmospheric layers of the disk. Large dust grains eventually drift radially inward toward the central star due to the gas pressure (Takeuchi & Lin 2002; Trapman et al. 2021). If nothing stops the radial drift, the CO ice is released inside the CO snow line, resulting in super-solar C/H ratios, as observed in HD 163296 (Zhang et al. 2019, 2020a).

However, Zhang et al. (2019) show that not all protoplanetary disks have an increase in carbon content inside the CO snow line, a result of carbon locking into less volatile species in the outer disk. Thus, for example, if CO is processed into CO₂, this carbon will not return to the gas until the CO₂ snow line, which is unresolved in most data. However, the situation may be more complicated than this: a study of accreting T Tauri stars demonstrates that the general trend is that carbon is depleted in the inner disk, which we define as the gas accretion disk inside the dust sublimation radius, typically at 0.01–0.1 AU; McClure (2019) finds depletion factors ranging between 2 and 25 relative to the diffuse interstellar medium (ISM) value of 1.35 × 10⁻⁴. The volatile and refractory depletion trend across different atoms in the well-studied inner disk of TW Hya (Bosman & Banzatti 2019; McClure et al. 2020) is the inverse of the super-solar abundances observed in Solar System meteorites (Lodders 2010; Bergin et al. 2015; McClure et al. 2020), which demonstrates that the observed carbon depletion likely results from the trapping of ice-coated rocky pebbles in the outer disk, either in a pressure maximum or in bodies that grow quickly to kilometer-sized planetesimals, stopping the radial drift. Since some of the disks in McClure (2019) are compact in the dust and do not show any signs of strong pressure bumps, it has been proposed that some of these sources have produced planetesimals via the streaming instability, creating an initial dust trap before any planets have carved a gap. However, the carbon abundance in the outer regions of the disks in the McClure (2019) survey has never been unambiguously determined. Therefore, it is still unclear whether the carbon abundance in the outer disk is even lower than that of the inner disk, as expected based on models that include radial drift, or if the dust is efficiently trapped in dust traps, pebbles, and planetesimals.

Observational estimates of the elemental volatile carbon abundance in protoplanetary disks are difficult to constrain. The main gas-phase carbon carriers in the atmospheric layers of protoplanetary disks are C⁺, C⁰, CO, and CO₂, ordered by increasing shielding from ultraviolet photons (Bruderer et al. 2012; Bergin et al. 2014; Bosman et al. 2018; Krijt et al. 2020). CO is the dominant reservoir in molecular gas in the bulk of the disk, but it is not the most accurate probe of the C/H abundance as most isotopologues are optically thick (see, e.g., Zhang et al. 2020b; Booth & Ilee 2020) and can be affected by direct photodissociation and freeze-out (Kama et al. 2016a). In the outer and upper regions of the disk where CO and small dust grains are depleted and the effects of photodissociation are magnified, neutral and ionized carbon become the main carbon reservoirs in the gas. Atomic carbon is not often observed in protoplanetary disks and may have non-disk emission components; however, in combination with spatially resolved CO emission, it is a better tracer of the elemental volatile carbon abundance than CO alone (Kama et al. 2016a).

In this work, we present [C I] observations using the Atacama Compact Array (ACA) of the Atacama Large Millimeter/submillimeter Array (ALMA) for a sample of seven disks, four of which have a known carbon abundance in the inner disk as determined from accreting material emitting in the near-infrared (McClure 2019). Three of these four disks were identified as having a dust trap on the basis of their inner disk carbon abundances. The sample, observations, and data analysis are described in Sect. 2 and the analysis of the observed spectra in Sect. 3. Using the atomic carbon emission in combination with archival CO data, we determine a radial profile of the elemental carbon abundance in the three sources with known inner disk carbon abundance. The modeling methodology is described in Sect. 4, the results of the modeling in Sect. 5, and the implications for the disk chemistry and formation of planetesimals in Sect. 6.

2 Observations and Data Analysis

The full sample contains seven well-studied disks, four in Taurus, two in Chamaeleon I, and one in Ophiuchus. The sample was selected to have strong near-infrared permitted C⁰ lines tracing the inner disk in either the McClure (2019) data or archival Very Large Telescope (VLT) X-shooter spectra. Subsequent criteria were that the sources already have CO isotopologue 12m ALMA observations from other programs that could provide additional constraints on our modeling. Table 1 lists the properties of the stars in our ALMA sample. Most of the stellar properties are taken from McClure (2019). The rest velocities of the sources are taken from high resolution CO data.

2.1 Observations

Our ALMA ACA Band 8 observations were conducted in August and September of 2018 as program 2017.1.00857.S (PI: M.K. McClure). The continuum emission was recorded in two spectral windows (SPWs), centered on 479.59 GHz and 477.7 GHz, respectively, each with a bandwidth of 2.0 GHz. One of the line SPWs was centered on the [C I] J = 1−0 line at 492.165 GHz, and the other line SPW was centered on the CS J = 10−9 line at 489.755 GHz in an attempt to constrain the sulfur chemistry in the outer disk. The latter emission line is not detected in any of the sources, and hence it is not used in this analysis. The bandwidth of the line SPWs were 250 MHz each, resulting in a spectral resolution of 122 kHz or 0.086 km s⁻¹. All disks are unresolved in the beam of the ACA, which is typically 3.5″ × 2.7″.

The selected sample was split into three observing groups, based on their position on the sky and the integration time that was needed to get to the required sensitivity. The sources in Taurus were observed in seven different executions for ~ 1–2 h integration time per source, the sources in Chamaeleon in five different executions for ~2.5 h integration time per source and AS 205 A in two different executions for 37 min in total. Full details of the observations and used calibrators are summarized in Table A.1.

2.2 Data Analysis

The data were reduced with the standard ALMA pipeline in CASA (McMullin et al. 2007). The sources in Taurus and AS 205 A were calibrated using CASA v5.1.1, the sources in Chamaeleon were calibrated with CASA v5.4.0. Self-calibration was performed on one of the continuum spectral windows of the combined executions, except for the sources FZ Tau and FM Cha that have a very weak continuum emission (see Table 3). The self-calibration tables were then applied to the line SPWs. For the sources that have high enough S/N to apply self-calibration, this procedure provided an improvement in continuum image quality of 50–100% and improved the line peak S/N by ~10%. In order to gain S/N in the [C I] J = 1−0 line, especially for the sources with weak line emission, we re-binned the Measurement Set using the cvel2 command in CASA, reducing the velocity resolution by a factor of ~2 to 0.15 km s⁻¹.

The continuum is subtracted in the uv-plane after deselection of the channels with line emission based on the dirty image. All sources, except FZ Tau where no [C I] J = 1−0 emission was detected, were imaged using the tclean task in CASA with Briggs weighting and a robust factor of 0.5. The choice of this robust factor is a trade-off between high sensitivity for the weak lines and a small beam to maximize the ability to discriminate between disk and cloud emission. The emission and absorption features were masked during the imaging procedure using hand-drawn regions in each channel. The line spectra were extracted using a circular mask with a radius of 2.5″, centered on the continuum peak, and are plotted in Fig. 1. Continuum fluxes, continuum RMS noise and the final beam size of each source are listed in Table A.1. The disk component of the line flux of the [C I] J = 1−0 emission was selected based on literature V_lsr values taken from CO emission (Table 1) and the modeling.

2.3 Complementary Data

In addition to the [C I] mini-survey, we analyzed archival ALMA 12 m array data of ¹³CO and C¹⁸O emission for the three sources that we modeled and studied in more detail: DR Tau, DO Tau and DL Tau. The details of the different molecular transitions can be found in Table 2. The Band 6 CO J = 2–1 isotopologue data were observed as part of program 2016.1.01164.S (PI: G. J. Herczeg). Details of the data and initial calibration can be found in Long et al. (2019). The ¹³CO and C¹⁸O J = 2–1 is imaged using natural weighting, with a typical beam-size in the channel maps of 0.14″ × 0.11″ at a velocity resolution of 0.2 km s⁻¹. Spectra were extracted using a circular aperture with a radius of 1.5″ for DL Tau and DR Tau and 1″ for DO Tau.

For DR Tau we included additional Band 7 CO J = 3−2 isotopologue data, as part of program 2016.1.00158.S (PI: E.F. van Dishoeck). This data set was calibrated using CASA v4.7.3. After the standard pipeline calibration, we performed self-calibration on the continuum and applied the gain solutions to the line channels. The beam-size in the final images is 0.24″ × 0.20″ at a velocity resolution of 0.055 km s⁻¹. Self-calibration on the continuum reduced the noise by a factor of 2 to 0.13 mJy/beam. The total source integrated continuum flux is 333.7 ± 0.2 mJy. The spectrum is extracted using a circular aperture with a radius of 2″.

3 Observational Results

Six out of the seven sources are detected in [C I] with high confidence. All spectra are presented in Fig. 1. Emission is seen to be within −2 and 2 km s⁻¹ for the sources in Taurus and AS 205 A and within −4 and 4 km s⁻¹ in the sources in Chamaeleon I. Therefore, we use all channels with a flux density above zero within these limits to construct the moment zero maps as shown in Fig. 2. The integration limits for each source are given in Table 3. In most sources where [C I] emission is detected, the emission is either extended or suffers from foreground absorption. Integrated line fluxes, as well as the disk contribution, are listed in Table 3. In the forthcoming section, we describe the emission of each source separately.

DR Tau is a relatively young source (0.9 Myr; McClure 2019, see Table 1) with a compact, almost face-on disk (dust R_out = 54 AU, i = 5.4°, Long et al. 2019) that has no signs of structure in the millimeter continuum emission image (Long et al. 2019). In the [C I] J = 1−0 spectrum we can identify three different emission components. One emission peak, centered at the cloud velocity, is emission associated with the disk. This emission is single-peaked as expected based on the low inclination. The disk component is spaced over five channels, which means that the emission is spectrally resolved. The disk has an integrated line flux of 0.92 ± 0.16 Jy km s⁻¹.

On the redshifted side (centered at 0.7 km s⁻¹) of the [C I] J = 1−0 spectrum, we detect an additional emission component from a stream of material that is likely infalling. We discuss this extra component further in Sect. 5. Some additional emission is observed on the blueshifted side of the spectrum. The [C I] channel maps at this velocity show small extended emission and absorption features (see Fig. D.3), but the spatial resolution is insufficient to see clear structure in the emission.

DO Tau is a 2 Myr old source (see Table 1) harboring a compact dust disk with an outer radius of 37 AU and an inclination of 27.6°, without any substructure in the millimeter continuum emission image (Long et al. 2019). This source has a clear detection in [C I] J = 1−0, with two peaks separated by an absorption feature at cloud velocity. This absorption feature can be caused by the presence of extended emission in the data beyond the maximum resolvable scale of the interferometer, foreground cloud absorption or a combination of both. The redshifted non-Keplerian emission centered at 1.5 km s⁻¹ is connected with extended emission in the channel maps that appears to be colocated with the blueshifted part of the pole-on outflow observed in ¹²CO (Fernández-López et al. 2020), as illustrated in Appendix B. The moment zero map of DO Tau shows very little emission, mostly due to the foreground cloud. For DO Tau, we determined two different fluxes by treating the absorption feature in two different ways. The integrated emission of the disk contribution is 0.59 ± 0.10 Jy km s⁻¹, which masks the absorption feature. Assuming a constant flux of 0.6 Jy at velocities where we observe cloud absorption results in an integrated flux of 1.18 ± 0.10 Jy km s⁻¹. The latter approach is more useful in the modeling and is justified based on our modeling and the clean line profiles of the other two sources that have an inclination similar or higher than DO Tau (AS 205 A and DL Tau).

DL Tau is an old source (7.8 Myr, McClure 2019, we discuss its age more specifically in Sect. 6.2) and has a large continuum disk with a 95% effective radius of 147 AU, an inclination of 27.6° and three dust rings (Long et al. 2018). On-source [C I] emission is detected without any signs of extended emission in the clean channel maps (see Fig. D.1). The spectrum reveals a slightly asymmetric double-peaked Keplerian profile, with the blueshifted peak slightly brighter than the redshifted emission. The integrated line flux of DL Tau is 1.36 ± 0.10 Jy km s⁻¹. Even though the source has a large dust disk, no CO isotopologue emission is detected, which suggests that this source is highly depleted in CO.

FZ Tau is a 0.9 Myr old source (see Table 1) with a small (~10 AU Piétu et al. 2014) dust disk and an inclination of <70°. FZ Tau has no detected [C l] J = 1−0 emission, but does show a similar absorption feature as DO Tau. The angular distance between FZ Tau and DO Tau is small (2.3°), which means that the observed absorption is probably extended cloud absorption. The most likely reason for the non-detection is the small disk size (~10 AU based on the continuum Piétu et al. 2014). FZ Tau has strong near-infrared permitted C⁰ emission (McClure 2019) that requires a large enough gas-rich hot region interior to the dust sublimation radius and sufficient ultraviolet levels to photodissociate CO. Unfortunately, the beam dilution is too large to detect the [C I] emission tracing the outer disk using ACA.

WW Cha is a 6.4 Myr old source (see Table 1) with a 135 AU structured dust disk at an inclination of 37.2° (Kanagawa et al. 2021). WW Cha has a [C I] spectrum with a classic P-Cygni profile, indicative of an outflow, consisting of a broad intense emission feature on the redshifted side of the spectrum together with a narrower and a weaker absorption feature in the blueshifted emission. Besides the outflow we detect a disk component near the cloud velocity and a ridge of infalling material on the blueshifted side of the spectrum (see also the channel maps in Fig. D.5). The source is known to have a jet (Robberto et al. 2012) and lately confirmed, using the Very Large Telescope Interferometer (VLTI), to be a close binary with a separation of 1.01 AU (Anthonioz etal. 2015; GRAVITY Collaboration 2021). The P-Cygni profile could be the result of an outflow caused by the binary interactions with the disk.

FM Cha is 5.3 Myr old source (see Table 1) that harbors a 82 AU dust disk at an inclination of 51° (Hendler et al. 2020). FM Cha has a firm detection in [C I] J = 1−0 with two asymmetric peaks in the spectrum on either side of the systemic velocity. The [C I] line in FM Cha is much broader than in AS 205 A, which we add in gray in Fig. 1 for reference, even though both disks have a similar inclination. Most of the emission is extended, contaminated cloud emission at high velocities (>1 km s⁻¹) relative to the local standard of rest. Although we observe some hints of Keplerian emission close to the source velocity, we cannot constrain the disk contribution in the emission due to the low spatial resolution of the [C I] observations.

AS 205 A is a binary system with a large, 53 AU, dust disk to the north (AS 205 N) that is at an inclination of 20.1° relatively face-on, and a small, 24 AU, dust disk to the south (AS 205 S) that is inclined at 66.3° (Kurtovic et al. 2018). The system shows a lot of substructures in CO, due to gravitational interaction between the two systems. The distance between the two disks is 1.3″, which means that the binary system is unresolved in the beam of our observations. The neutral carbon emission is singly peaked, centered on the systemic velocity of the largest disk AS 205 N. The systemic velocity of the secondary star is offset by −2.5 km s⁻¹ in CO (Kurtovic et al. 2018), as indicated in Fig. 1. There is no significant emission detected at this relative velocity, which means that AS 205 S is not detected and that it is reasonable to assume that all the observed emission comes exclusively from AS 205 N.

Concluding, four disks out of the seven targets have been detected in [C I] J = 1−0 in this study. In the following we describe the modeling of the three Taurus disks in the sample for which the [C I] J = 1−0 emission can be easily decomposed. These three sources already have well-determined inner disk C abundances in the literature, which makes them more suitable for studying the radial profile of the carbon abundance.

Fig. 1 Spatially integrated emission spectra of the [C I] J = 1−0 line at 492.165 GHz for the seven sources in our sample. All the spectra are shifted to their respective local standard of rest velocity (Table 1). The shaded region marks the 1σ uncertainty for each source. The spectrum of AS 205 A is shown in gray with FM Cha as a reference of the typical linewidth for a modestly inclined disk. “N” and “S” depict the Vlsr of AS 205 Ν and AS 205 S, respectively.

Fig. 2 Moment zero maps of the [C I] J = 1−0 emission. Solid white contours show the 0.61 mm continuum at 10, 30, 100, and 300σ. The dashed white contour indicates the region used to extract a spectrum. Beam sizes are shown in the lower-left corner of each panel. The 3σ detection limit is illustrated as a white tick mark in the colorbar.

4 Description of the Modeling

The aim of our modeling was to develop a physical disk model tailored to each source that describes both the dust and the gas structure in each disk. Using this model, we could then determine the volatile carbon abundance in the outer disk by comparing the flux in the [C I] J = 1−0 observations with a grid of models, varying carbon and oxygen abundances. In order to determine the physical structure of the disks, we combined information from the spectral energy distribution (SED) and the outer radii of the disks in ALMA millimeter continuum and CO observations. Although the model may not be unique, it constrains a structure that agrees well with the observations. A detailed overview of the impact of changes in variables of the physical disk model on the line flux in [C I] J = 1−0 is available in Kama et al. (2016a).

The azimuthally symmetric physical-chemical Dust And LInes (DALI) model (Bruderer et al. 2012; Bruderer 2013) is used to model the sources using a smooth parameterized gas and dust disk density structure, neglecting all azimuthal substructures, gaps and rings seen in the dust. DL Tau is known to have multiple gaps in the millimeter continuum (Long et al. 2018), but the influence of the millimeter dust on the gas temperature in the outer disk is small (see, e.g., Facchini et al. 2018; van der Marel et al. 2018; Alarcon et al. 2020), so it is safe to use a smooth disk density profile to determine the gas temperature and chemistry in the outer disk. Both DO Tau’s and DR Tau’s disk are smooth in the millimeter dust. DALI first solves the radiative transfer of the continuum based on a dust density structure (Sect. 4.1) and an input external radiation field (Sect. 4.2) to determine the dust temperature and radiation field from ultraviolet to millimeter wavelengths at all locations in the disk. After that, the heating-cooling balance of the gas is solved simultaneously with steady state chemistry in an iterative sequence until convergence is reached. Once the solution is converged, we ray-trace the model in continuum and molecular lines using non-local thermodynamic equilibrium radiative transfer (Sect. 4.4).

4.1 Disk Parameters

The density structure of our model is fully parameterized and is based on a full disk, without inner cavity. We used the standard form of the dust surface density profile of a power law with an outer exponential taper: $${\displaystyle {\sum }_{\text{dust}}={{\displaystyle {\sum }_c\left(\frac{r}{R_c}\right)}}^{-\gamma }\exp \left[-{\left(\frac{r}{R_c}\right)}^{2-\gamma }\right],}$$(1)

defined by the surface density Σ_c at the characteristic radius R_c, and assuming an initial power-law index γ =1. The scale height h at distance r is given by $$h=h_c{\left(\frac{r}{R_c}\right)}^{\psi },$$(2)

defined by the scale height h_c at the characteristic radius with a flaring index ψ. The dust was split into two populations following the approach of D’Alessio et al. (2006), a small dust population that ranges from 0.005 to 1 micron and a large dust population that ranges from 0.005 to 1000 micron. The vertical structure of the small dust grains is $${\rho }_{\text{d,small}}=\frac{\left(1-f\right){\Sigma }_{\text{dust}}}{\sqrt{2\pi }rh}\exp \left[-\frac{1}{2}{\left(\frac{\pi /2-\theta }{h}\right)}^2\right],$$(3)

where f is the mass fraction of the large grain distribution, assumed to be 0.99, and θ is the opening angle from the mid-plane as seen from the central object. The large grains are settled with respect to the full disk height with a factor χ ∈ (0, 1), following the same distribution as the small grains (Eq. (3)) but with f replacing (1 − f) and χ × h replacing h. The radial distribution of the two dust grain populations is the same, excluding grain evolution processes. The gas density distribution is scaled up from the dust density distribution using the ISM gas-to-dust ratio of 100. The vertical gas structure follows the distribution of the small grains, preserving the global gas-to-dust ratio. We make an initial estimate of the sublimation radius in the disk using Eq. (14) in Dullemond et al. (2001). By assuming a dust sublimation temperature of 1500 K, negligible self-irradiation (H_rim/R_rim ≪ 1), and solar units for the luminosity, Eq. (14) reduces to r_subl [AU] = 0.07 (L[L_⊙])^1/2 (see also Monnier & Millan-Gabet 2002). For DO Tau we increased this result to 0.14 AU to match the mid-infrared emission from the inner disk, similar to the value found for the Brγ area for DO Tau as determined in Eisner et al. (2014).

4.2 Stellar Parameters

The stellar spectra are approximated by a perfect blackbody defined by the effective temperature (T_eff) of each star with an additional component due to accretion, assuming an accretion temperature of 10 000 K of the emitting radiation, using $$L_{\text{acc}}(v)=\pi B_v\left(T_{\text{acc}},v\right)\frac{G{M}_{*}M}{R_{*}}\frac{1}{\sigma T_{\text{acc}}^4}.$$(4)

The observed far infrared luminosity of the sources (Yang et al. 2012) agree within an order of magnitude with this first order approximation.

The properties of all modeled stars are taken from McClure (2019), which uses the procedures from McClure et al. (2013) to derive self-consistent stellar and accretion parameters. We included a X-ray luminosity of 10³⁰ erg s⁻¹ at a X-ray plasma temperature of 7 × 10⁷ K, and a cosmic ionization rate of 5 × 10⁻¹⁷s⁻¹ incident on the disk surface.

4.3 Chemical Network

The standard chemical network adopted in DALI is based on UMIST06 (Woodall et al. 2007). This chemical network includes neutral and charged polycyclic aromatic hydrocarbons, 109 molecular species and 1463 reactions. The code includes two-body interactions, freeze-out, photodissociation and – ionization, and thermal- and photo-desorption. There are no grain surface reactions considered, except for the hydrogenation process. It is reasonable to assume that this has no impact on the species that are used in this paper, which trace the atmospheric layers of the disk. For DR Tau and DO Tau we include the full treatment of CO freeze-out and isotope selective processes as described by Miotello et al. (2014, 2016). We adopt fiducial ISM gas-phase elemental carbon and oxygen abundances of C/H = 135 parts per million (ppm) and O/H = 288 ppm, respectively. These values are close to the median abundances observed in translucent interstellar clouds (Cardelli et al. 1996; Meyer et al. 1998; Parvathi et al. 2012) and consistent with the solar values within a factor of a few (Asplund et al. 2009). A full justification of these fiducial abundances is provided in Bruderer et al. (2012). In order to make the carbon abundances easy to interpret, we define a carbon depletion factor (f_depl) as (C/H)_ISM/(C/H). When changing the carbon abundance we vary the oxygen abundance such that the C/O ratio stays at the ISM level of 0.47, unless specified differently. Without a similar depletion of oxygen w.r.t. carbon, the neutral atomic carbon abundance would be more strongly dependent on the elemental carbon abundance, as there would be more oxygen available to create CO.

4.4 Radiative Transfer

After the chemical network has converged to a steady state equilibrium and gas temperature solution, atomic and molecular transitions can be imaged using the ray-tracer in DALI, which includes the radiative transfer. The excitation of the line is determined explicitly in a non-local thermodynamic equilibrium fashion, using the collision rate coefficients from the Leiden Atomic and Molecular Database (LAMDA) (Schöier et al. 2005; van der Tak et al. 2020). In this step, the sources are projected on the sky using the observed distance, inclination, and position angle (PA) given in Tables 1 and 4, respectively. For comparison with our observations, we convolved the model with a Gaussian beam with the same shape as the beam of the ALMA observations.

4.5 Fitting Procedure

In the fitting procedure, we followed a similar approach as Bruderer et al. (2014), creating a physical model that fits the SED through consecutive parameter value refinements, until they converge to a model that is best representative of the source. Using a χ² or Markov chain Monte Carlo method was not possible because of the long computation time of the continuum radiative transfer. We used the same characteristic radius (R_c), inclination (i) and PA for the gas as observed in high angular resolution dust continuum observations. The gas disk can be significantly larger than expected from millimeter continuum observations if radial drift and grain growth are important in the disk (Facchini et al. 2017; Trapman et al. 2019), but for the sake of simplicity we do not include dust evolution in our modeling. The high spatial resolution submillimeter continuum constrains the mass, R_c, i and the PA of the dust in the disks. Additionally, we look at the spectral line profiles of the CO isotopologue and [C I] emission to further constrain the density and temperature structure of the disks. These properties together constrain Σ_dust, γ, χ, f, h_c and ψ.

4.6 Main Sources of Uncertainty

We stress that the gas-phase carbon abundances found in the modeling analysis are subject to a few uncertainties. First, we assume a smooth disk with a system-averaged elemental volatile carbon depletion. It is thought that most CO processing preferentially occurs at the deepest, coldest layers (Bosnian et al. 2018), while the tracers that we analyze, such as CO and C⁰ are known to trace the atmospheric layers of the disk (z/r > 0.2 Bruderer et al. 2012). Therefore, both radial and vertical dependence in (C/H)_gas can cause discrepancies in the abundance found in the system-averaged CO isotopologue and [C I] emission. Based on the carbon abundance measured in the upper layers of the disk, we infer that the bulk of the disk can be depleted in carbon, which is a result from efficient mixing by turbulence. Additionally, abundances of carbon-carrying molecules depend on the elemental oxygen abundance in the disk, since an excess in carbon atoms can shift the balance between CO₂, CO, C, and hydrocarbon abundances. Elevated C/O ratios are not uncommon in disks and can arise by a stronger depletion of oxygen, because of freeze-out and processing of CO₂ and H₂O rather than CO (which has been inferred directly from water observations; see Bergin et al. 2010; Hogerheijde et al. 2011; Du et al. 2017) or direct photo-ablation of carbon rich grains (Bosnian et al. 2021). We assumed ISM C/O ratio where possible and discuss a potential elevated C/O ratio in Sect. 5.2.

Second, there is an uncertainty in the (C/H)_gas found in these three systems, resulting from the uncertainty in the modeling parameters. The biggest sources of uncertainty are the disk mass, R_c and the gas-to-dust ratio (Kama et al. 2016a). There are no HD observations of the sources in our sample, which means that there is high uncertainty on the gas disk mass and the gas-to-dust ratio. Based on the few sources that are observed in HD, we do not have reasons to assume that the gas-to-dust ratio in the bulk of the disk is much different from the ISM value of 100 (McClure et al. 2016; Bergin et al. 2013). Kama et al. (2016a) show that an underabundance of gas-phase carbon by an order of magnitude, such as found here, cannot be easily masked by a choice of parameters. Larger depletions will be even more reliably identified.

Fig. 3 SEDs of the three modeled sources with our model fit in orange. The black dots with error bars represent observations, the dashed black curve denotes the stellar component, and the solid black curve is a combination of the stellar component and the ultraviolet excess due to accreting material as used in the modeling.

5 Modeling Results

Our best-fit results to the SED for the dust density structure are presented in Fig. 3. All data points in the SEDs are de-reddened using the CCM89 extinction curve (Cardelli et al. 1989) for sources with A_V < 3 (DR Tau and DL Tau). For DO Tau we used the extinction curve for 0.3 ≤ A_k < 1 from McClure (2009) to correct the emission at wavelengths shorter than 25 micron. This source has A_V = 3.6 magnitudes (McClure 2019, and see our Table 1), which are higher than the range over which the CCM89 extinction curve is accurate. Reasonable fits are found for the SED, including spatial information from the submillimeter continuum. The inferred disk gas masses (Table 4) are within a factor of a few of the estimated disk mass based on multi-wavelength continuum observations and assuming a gas-to-dust ratio of 100 (see for example Kwon et al. 2015 for DL Tau (4.1 ± 0.5 × 10⁻² M_⊕) and DO Tau (1.4 + 0.1 × 10⁻² Μ_⊕) and Tazzari et al. 2016 for DR Tau (1.4 ± 0.4 × 10⁻² Μ_⊕)).

The mid-infrared part of the SED of DR Tau is underestimated by the best representative model. Varying h_c between 0.05 and 0.3, ψ between 0.05 and 0.2 and R_Sub1 between 0.05 and 0.3 AU, we find no plausible configuration of the parameters that agrees with the inner disk contribution of the SED < 10 micron. The strong mid-infrared excess in the SED could be caused by an inner disk gas component that is optically thick to its own radiation, producing continuum emission. Evidence for a gas contribution to the continuum interior to the dust sublimation radius in similar systems can be found in many studies (see, e.g., Kraus et al. 2008; Malygin et al. 2014; Lazareff et al. 2017). Since we focus on outer disk tracers, we do not expect the accreting gas column to have a strong influence on the results. It should be noted that our lack of a fit of the inner disk’s contribution to the continuum increases the model parameter space for fitting the SED at longer wavelengths.

Figure 4 presents the results of the total [C I] J = 1−0 line flux as a function of the elemental gas-phase carbon abundance for the three models. Observations are presented as black stars with their corresponding depletion factor and uncertainty in the figure legend. Model carbon abundances are presented as dots at depletion factors of [1, 3, 10, 30, 100, 300]. The derived depletion factor of the sources is the interpolation of the connecting curves with the error set by the measurement error on the [C I] J = 1−0 flux. A decrease by a factor of 10 in the elemental C abundance results typically in a factor of 3 in the total line flux. The CO emission and spectra will be discussed source individually in the forthcoming sections.

Fig. 4 Modeling results of the total flux as a function of the carbon depletion factor defined as fdepl = (C/H)/(C/H)ISM. Observations are plotted at the interpolated carbon depletion fraction. For DO Tau, we plot two values for the total flux: the star is the corrected flux for absorption with a constant value, and the diamond is only masked for absorption, as discussed in the main text. The derived carbon depletion fraction by interpolation with the model is presented in the top-right corner, with errors given by the measurement uncertainty of the [C I] observations.

Fig. 5 Emission line profiles of 13CO J = 2−1,13CO J = 3−2, C18O J = 2−1, C18O J = 3–2, and [C I] J = 1−0 for DR Tau, with the best fitted model overlaid in color. 1σ error bars are indicated in gray for each spectrum. The data and models of the weak CO isotopologue J = 2−1 transitions have been multiplied by a factor of 5 for illustrative purposes.

5.1 DR Tau: Moderate Depletion of Carbon and Oxygen

The results of the line emission fit of DR Tau are presented in Fig. 5 for the best-fit model with a carbon depletion factor of 5. A map of the surface brightness of the ¹³CO J = 3−2 transition is presented in Fig. 6. Using least squares fitting, a Gaussian disk component with an additional Gaussian component centered at 0.7 km s⁻¹ accounting for the non-disk emission is fitted to the [C I] and ¹³CO J = 3−2 spectra to determine the disk contribution.

The low ratio between the [C I] and CO isotopologue fluxes of the disk could only be explained by assuming γ < 1. Decreasing γ in the model results in a sharper exponential cut-off from the surface density, which affects the emitting region of the [C I] emission more than the CO emitting layer. Following the approach described in Miotello et al. (2018), we fitted the slope of the emission in the so-called pivot-region of the ¹³CO and C¹⁸O emission where the line is optically thin, but has enough column density to efficiently self-shield against isotope selective processes. We find γ = 0.5 and γ = 0.8 for the ¹³CO and C¹⁸O radial intensity profiles, respectively. In the modeling we used γ = 0.5, a value that is consistent with theoretical values for a viscously heated disk (Chambers 2009).

The C¹⁸O fluxes together with [C I] provide the strongest constraints on (C/H)_gas. Therefore, we include the isotopologueselective photodissociation of CO (Miotello et al. 2014, 2016), which prevents underestimating (C/H)_gas or the gas mass. We adopt ISM-like ¹³C and O¹⁸ abundance ratios of 1.753 ppm and 0.5143 ppm with respect to H (Wilson 1999). A model with (C/H) = 27 ppm, or a depletion factor of 5, reproduces the [C I] emission and most CO isotopologue lines. It overproduces the¹³CO. J = 2−1 line by a factor of 3. The modeling suggests similar radial emitting regions for both CO transitions. Using the peak ratio between the ¹³CO and C¹⁸O J = 2−1 spectra we determine the opacity of the line and the corresponding kinetic temperature using Eqs. (1) and (2) from Schwarz et al. (2016). The peak ratio between the two CO isotopologue transitions is only 1.45, while a ratio of 8 is expected based on the ISM abundances of ¹³C and ¹⁸O. This results in an opacity of 9.4 and a temperature of the emitting layer of 42 K. A similar analysis for the CO J = 3−2 isotopologue transition results in a peak ratio of 2.85, an opacity of 3.3 and a kinetic temperature of the emitting layer of 30 K. The reason why CO J = 2−1 emission seems to be much more optically thick, despite having a much lower Einstein A coefficient, could be the lack of short spacing in the uv-plane, which means that the data are not sensitive to large-scale structures that are expected to be more dominant in ¹³CO than in C¹⁸0 emission. In depth analysis of the temperature structure of the disk is beyond the scope of this paper. We discuss the CO isotopologue emission and the methods used to constrain γ in more detail in Appendix C.

The best-fit carbon and oxygen abundances are found to be 2.7 × 10⁻⁵ and 5.76 × 10⁻⁵, which corresponds to a carbon and oxygen deficiency of a factor of $5_{-1}^{+2}$. The analysis of the photospheric abundance and the corresponding depletion factor found in the inner disk of DR Tau is ${13}_{-1}^{+2}$ (McClure 2019). The low abundance in the inner disk in comparison with the outer disk is consistent with radial trapping of carbon-rich ices beyond the CO snow line. We return to this in Sect. 6.

The non-Keplerian, redshifted emission that is observed in both the [C I] J = 1−0 and ¹³CO J = 3−2 spectrum at a velocity of 0.7 km s⁻¹ originates likely from the circum-disk material that is also visible in the moment zero map of the ¹³CO emission to the northeast (indicated in Fig. 6). The high [C I]/¹³CO 3–2 flux ratio of this feature suggests that the emission is coming from material that is elevated above the disk surface and thus exposed to a high ultraviolet flux. The contaminating material is confined to a single velocity, is not kinematically connected to the disk and may be linked with the late infalling material seen in scattered light observations (Mesa et al. 2022). The C¹⁸O 3−2 line is slightly broader than expected, likely caused by contamination with the same accreting material, but here we do not try to make a two-component fit. Increasing the C/O ratio in a similar manner to our treatment of DO Tau in Sect. 5.2 could improve the fit to the ¹³CO J 2−1 spectrum. However, values of up to 3 are insufficient to fit the observed line flux. Furthermore, this would affect the higher S/N ¹³CO J = 3−2 emission as well, which is reasonably matched by models with the local interstellar standard C/O.

Fig. 6 Total intensity or moment zero map of the 13CO J = 3−2 emission in DR Tau. The black contours show the continuum emission at 0.8 mm. The gray contours show the residuals after subtracting an azimuthally averaged radial disk profile to highlight the infalling stream. The beam size is shown in the lower-left corner.

5.2 DO Tau: Moderate Depletion of Carbon and Oxygen

The data and best representative modeled spectra for DO Tau are illustrated in Fig. 7. We neglected the redshifted excess >1 km s⁻¹ that is associated with a pole-on outflow, as explained in the observation section and illustrated in Appendix B.

The ¹³CO 2−1 line is optically thick, with a ratio of 2 with respect to the total C¹⁸O 2−1 flux, while a ratio of 8 is expected based on the relative abundance of the two isotopes in the diffuse ISM. We present two models in Fig. 7, one with f_depl = 10 and an ISM-like C/O ratio and one with f_depl = 20 and an enhanced C/O ratio of 1. In the later model, oxygen is less abundant with respect to carbon than the first model. If oxygen is less abundant in the model but the carbon abundance remains constant at the ISM level, less CO reacts on the ice with OH to form CO₂ and H (e.g., Ioppolo et al. 2011) and thus the total CO abundance increases. By using an increased C/O ratio there will be more carbon in atomic state, because there will be more CO in the atmospheric layers to be photodissociated in C° and C⁺ and less oxygen to convert C back in CO. Since ¹³CO is optically thick, changing the oxygen abundance does not change the total line flux of the model significantly. The best fitted C/O value based on the [C I] and CO isotopologue emission is ~ 1, which is about twice the ISM value.

C₂H can provide much stronger constraints on the C/O ratio than [C I] and CO, as it is predominantly formed in regions where C/O ratios are high (Bergin et al. 2016; Miotello et al. 2019; Bosman et al. 2021). The observed upper limit of C₂H = 73 mJy km s⁻¹ found by Bergner et al. (2019) is low compared to the other sources in their sample, and more consistent with C/O close to the ISM value. To prove this we ran an extended chemical model based on the chemical model from Visser et al. (2018) that includes C₂H and C₂H₂ chemistry using the same temperature and density structure as for the atomic carbon and CO emission. We find that for a C/O ratio of 1 and carbon depleted by a factor of 20 with respect to the ISM value, the C₂H flux would be 550 mJy km s⁻¹. An ISM-like C/O ratio of 0.47 and a carbon depletion fraction of 10 results in a total flux of 64mJykm s⁻¹. C₂H traces regions in the disks that are closer to the star than [C I], which could mean that the discrepancy is caused by a radially varying C/O ratio in the disk. However, the systematic uncertainty of the models are too large to say anything conclusive about that, and full modeling of C₂H is beyond the scope of this paper. Deeper C₂H observations are needed to better constrain the C/O ratio. In the rest of this paper we use a C/H ratio that is ${17}_{-2}^{+3}$ times lower than the ISM value, set by the carbon abundance grid in Fig. 4.

Our carbon abundance in the outer disk of DO Tau is similar to, or higher, than the abundance in the inner disk, which is found to be depleted with gas-phase carbon by a factor of ${22}_{-11}^{+20}$ (McClure 2019). This radial pattern is consistent with retaining solids beyond the CO snow line. We return to this in Sect. 6.

Fig. 7 Observations and modeling of DO Tau. Two models are shown: the dashed line is a model with an ISM C/O ratio of 0.47 and a carbon depletion factor of 10, and the solid line is a model with an enhanced C/O ratio of 1 and a volatile carbon depletion fraction of 20. 1σ error bars are indicated in gray for each spectrum. The data and models of the 13CO and C18O J = 2−1 transitions have been multiplied by a factor of 5 and 20, respectively, for illustrative purposes.

5.3 DL Tau: Severe Depletion of Carbon and Oxygen

The data and best representative modeled spectra for DL Tau are illustrated in Fig. 8. The [C I] spectrum has a disk component that consists of a clean double-peaked Keplerian profile, without cloud contamination. The two CO isotopologue J = 2−1 lines are both upper limits, but the upper limits are low, so that they give additional constraints on the gas-phase carbon abundance.

The best representative model for DL Tau has R_c = 112 AU, similar to the value found for the dust disk (Kwon et al. 2015; Long et al. 2019). However, convolving the modeled [C I] emission with the same beam as the observations and comparing the azimuthally averaged radial profiles suggests that the radial extent is two times larger than the ACA beam of ~3″, while the observations appear as an unresolved point-source down to the same RMS level as the data. The most likely explanation to this dichotomy is that the emission in the outer parts of the disk is buried in the noise. To be consistent with the spatial information of the [C I] emission, we find that R_c has to be smaller than 65 AU. The SED constrains the total mass, but does not set limits to R_c. The characteristic radius of the dust is 112 AU (Kwon et al. 2015; Long et al. 2019), and the CO outer disk radius is 600 + 150 AU in Submillimeter Array (SMA) observations, which is consistent with R_c = 112 AU. On top of that, a model with R_c = 112 AU fits the [C I] J = 1−0 line profile better than a model with the same mass, but with R_c = 65 AU (see Fig. 8). The line profile is, besides the disk mass distribution, also dependent on the mass of the star. There is some uncertainty on the masses of stars, based on their evolutionary track, but the stellar mass would have to be off by a factor of 2–3 to fit the Keplerian profile, which is beyond the usually assumed uncertainties for the stellar mass (Baraffe et al. 1998).

For the R_c = 112 AU model, we find a best fitting depletion factor of ${157}_{-15}^{+17}$. If the disk were smaller, with R_c = 65, DL Tau would still be depleted in carbon by a factor of ${52}_{-8}^{+6}$. The high depletion value for DL Tau is consistent with the age of the source, as shown by chemical modeling of CO emission (Krijt et al. 2020). We return to this in Sect. 6.2. In order to put further constraints to the model for DL Tau and exclude that the disk is depleted in gas, resulting in a small outer gas radius, we need ALMA ¹²CO and CO isotopologue data and high resolution [C I] observations with more sensitivity in the outer parts of the disk.

Fig. 8 Observations and modeling of DL Tau. Two models are plotted: the solid line is a model with Rc = 112 AU as observed in the dust and the dashed line as a model with Rc = 65, which fits the spatially resolved [C I] emission best. Further details can be found in the main text. 1σ error bars are indicated in gray for each spectrum. The data and models of the weak CO isotopologue J = 2−1 transitions have been multiplied by a factor of 5 for illustrative purposes.

5.4 Systematic Uncertainty

Additional models are presented for the three sources in which the most relevant disk parameters are varied in order to obtain an indication of the systematic uncertainty on the carbon abundance in the modeling, as discussed in their respective subsection. These additional models show that the carbon depletion factor can be trusted within a factor of 3, with the assumption that the gas-to-dust ratio is equal to the canonical value of 100. Kama et al. (2016a) showed how the C/H abundance systematically depends on a large number of other disk parameters, including the dust-to-gas ratio. They found that in the worst case scenario, when taking into account a possibly varying dust-to-gas ratio, the [C I] flux in a disk is robustly sensitive to carbon depletions of at least a factor of 10 in the outer disk. Therefore, we can conclude that the significant carbon depletions in DO Tau and DL Tau are robust to the systematic uncertainties. In contrast, the order of magnitude carbon depletion in DR Tau, while significant under our current model assumption of a gas-to-dust ratio of 100, could be consistent with the ISM carbon abundance if the full range of systematic uncertainties in the DALI models is taken into account.

6 Discussion

We present the first successful ALMA mini-survey of neutral carbon to constrain the elemental carbon depletion in the outer disk of a sample of seven Τ Tauri disks. Two sources in our sample, DR Tau and DO Tau, have moderate depletion of volatile carbon and oxygen in the outer disk, and one source, DL Tau, has severe depletion. The observed sources show two distinct radial carbon abundance behaviors: the moderately depleted sources have carbon abundances in the outer disk that are consistent with the inner disk, while the severely depleted source has a much higher carbon abundance in the inner disk than in the outer disk. This work strengthens the case that some sources are depleted in gas-phase carbon in the outer disk in addition to direct freeze-out and photodissociation of CO. Figure 9 summarizes our results of the carbon abundance in the outer disk compared to the inner disk. We added TW Hya (Zhang et al. 2019; Kama et al. 2016b; McClure et al. 2020), HD 163296 (Kama et al. 2015; Zhang et al. 2019), and HD 100546 (Kama et al. 2015,2016b) and a combined point for ten Class 0/I sources (Harsono et al. 2014; Zhang et al. 2020b; Bergner et al. 2020) to the figure for reference. These three sources are the only other sources for which the carbon abundance is determined in both the inner disk and outer disk. For the Herbig Ae/Be stars we make the assumption that the composition of the stellar photo-sphere is representative of the accreting material in the inner disk, as in Kama et al. (2015).

6.1 Implications of the Radial Abundance Profile

Almost all sources with observed elemental carbon abundances have a considerable increase in carbon abundance in the inner disk with respect to the outer disk, except for the two young and compact disks in our sample, DR Tau and DO Tau (see Fig. 9). Large dust grains, the formation of which causes gaseous carbon depletion, will eventually drift inward toward the star, driven by the headwind of the gas that moves slightly sub-Keplerian due to the pressure support. With efficient drifting of dust grains and without any dust trapping, the carbon will be released inside the CO snow line. Any carbon that is processed into less volatile or refractory structures will ultimately sublimate at the dust sublimation rim. The inner disk gas (inside the dust sublimation rim) is constantly replenished with carbon-poor gas from the outer disk, which means that an enhancement there has to be a result of the sublimation of carbon-rich ices as a result of radial drift.

Fig. 9 Radial profile of the volatile carbon depletion fraction based on one point inside the dust sublimation rim and a disk-averaged carbon abundance in the outer disk. The outer depletion fractions of the sources in blue are from this work. Inner disk carbon abundances of these three sources are taken from McClure (2019). We added TW Hya in green (Zhang et al. 2019; Kama et al. 2016b; McClure et al. 2020) and one point in the outer disk for ten Class 0/I sources in yellow (Harsono et al. 2014; Zhang et al. 2020b; Bergner et al. 2020) to illustrate the increase in the level of carbon depletion during the system’s evolution. Two Herbig Ae/Be stars, HD 163296 (Kama et al. 2015; Zhang et al. 2019) and HD 100546 (Kama et al. 2015, 2016b), are shown in red to complete the picture of carbon depletion that we have. These sources do not follow the evolutionary trend, because they are much warmer.

6.1.1 Dust Traps

Strong pressure maxima and efficient planetesimal formation locations, incorporating the grains in larger bodies, can both stop the grains from radial drifting. DL Tau has three rings in the dust continuum observations that are located outside the CO snow line (Long et al. 2019). Sufficiently large dust grains will be trapped in the pressure maxima at the edge of these gaps, so that those grains remain in the outer disk, but trapped very locally, which means that the CO ice that has formed on these grains cannot efficiently return to the gas phase. The dust trapping in the pressure maxima that are observed in DL Tau are efficient enough to prevent a super-solar carbon abundance in the inner disk, contrary to what is observed in the radial CO abundance profiles of other protoplanetary disks such as HD 163296 (Zhang et al. 2019, 2020a). However, the fact that we see an enhancement in the inner disk with respect to the outer disk implies that the traps are “leaky” (i.e., not 100% efficient at trapping dust), that they formed in a later stage of the disk evolution (i.e., a shorter time ago than the viscous time needed to accrete the material from the dust traps to the inner disk), or that carbon depletion mechanisms started later than the viscous time. The last option is unlikely as recent studies pointed out that volatile CO rapidly evolves during the Class I and early Class II stage (Zhang et al. 2020b; Bergner et al. 2020).

The formation of a perfectly efficient dust trap at the CO snow line would prevent depleted carbon in the form of CO ice from entering the gas interior to the CO snow line. However, before it is observable via the near-infrared lines used by McClure (2019) this depleted carbon would still need to viscously accrete from the CO snow line to ~0.1 AU on a timescale of $$t_{\text{visc}}=\frac{r^2\Omega \text{K}}{\alpha c_{\text{s}}^2},$$(5)

with Ω_K the Keplerian angular velocity and c_s the sound speed in the midplane of the disk given by $\sqrt{k_{\text{B}}T/(2.3m_{\text{p}})}$ with k_B the Boltzmann constant, Τ the midplane temperature, m_p the proton mass, and α the viscosity parameter (Shakura & Sunyaev 1973). The innermost gap outside of the CO snow line in DL Tau is located at ~40 AU, and the temperature at the midplane is provided by DALI. However, α is not determined directly for DL Tau and is subject to large uncertainty. Observational measurements of a from turbulent broadening of CO J = 2−1 emission lines range from α = 0.08 in DM Tau to α < 0.006 in MWC 480, while theory predicts still lower a around 10⁻⁴ (Flaherty et al. 2020, and references therein). Using α = 0.08 and α = 0.006, we find that the formation of the dust trap at the CO snow line should be less than 0.16–2.2 Myr ago. This could not explain DL Tau’s carbon-depleted inner disk if the substructure were formed at the start of the Class II phase (~7 Myrs ago), as Krijt et al. (2020) showed that for a dust trap with >99% trapping efficiency, the inner disk would be depleted within 3 Myr by more than a factor of 3 with respect to the ISM abundance that is observed in DL Tau. However, if α were <0.003, which still agrees well with the observations, the viscous timescale would be long enough to cover even substructure that has formed right at the start of the Class II phase.

An alternative explanation for the large discrepancy between the carbon abundance in the inner and outer disk is that the dust traps in DL Tau are not 100% efficient in trapping dust. This idea is supported by modeling that shows that low drift efficiency will result in less enhancement in the inner disk (Krijt et al. 2020). McClure et al. (2020) devised a simple analytic model to calculate dust retention efficiencies. The total dust trap efficiency at different positions in the TW Hya disk, based on the distribution of carbon, is 99.8%, hence the strong depletion in the inner disk (Fig. 9). Following Eqs. (1)–(3) in Appendix C of McClure et al. (2020), we derived a simple expression for the dust trapping efficiency of our disks beyond the inner dust sublimation radius: $$f_{\text{trap}}=1-\frac{X_{\text{g,in}}-X_{\text{g,out}}}{X_{\text{g,ISM}}-X_{\text{g,out}}},$$(6)

where X_g,in is the carbon abundance in the inner disk, X_g,out the carbon abundance in the outer disk and X_g,ISM the carbon abundance in the ISM. Using this formula we find that the dust trap efficiency of the DL Tau disk is 61%, assuming an inner disk abundance that is lower by a factor of 2.5 with respect to the ISM value (McClure 2019). We note that this approximation is highly susceptible to variation in the carbon abundance close to the solar value and is hence an upper limit. If no carbon depletion is assumed in the inner disk, which is still within the systematic uncertainty of the measurement, then the locking fraction would be 0%.

Long et al. (2020) find low spectral index, α = 2−2.5, in large regions of the disk of DL Tau by comparing resolved 1.3 and 2.9 mm ALMA continuum data. They find similar disk outer radii at both wavelengths and show that this is consistent with a dust evolution model that includes a long lasting pressure bump. A low spectral index is in line with advanced grain evolution, which is also expected based on the evolutionary stage of the source. These authors show also that the spectral index in DL Tau varies only by ~10% near the pressure maxima, which is in line with the low trapping efficiencies of the dust traps we find in this work. We note, however, that this low variation is a lower limit, since the gaps are unresolved and therefore less pronounced due to beam smearing.

6.1.2 Dust Retention in Smooth, Compact Systems

The carbon abundances found for DR Tau and DO Tau are consistent within the observational uncertainties with what is observed in the inner disk, based on the accreting material onto the star (McClure 2019). This implies either that the carbon in ices on dust grains is prevented from drifting toward the star or that the chemistry outpaces the pebble drift by converting CO to less volatile species such as CO₂ and CH₃OH, preventing a resurgence of gas-phase CO just inside the snow line (Krijt et al. 2020). Based on the approximation of the locking fraction in Eq. (6), both DR Tau and DO Tau have a trapping efficiency close to 100% (97% and >100% for DO Tau and DR Tau, respectively). Both freeze-out and dust evolution scale with the orbital timescale, which means that carbon is more efficiently depleted closer in to the star, but outside the CO snow line. In extreme scenarios, this could result in a lower abundance in the inner disk than in the disk-averaged abundance if radial drift is stopped close to the CO snow line and the inner disk is accreting the gas that is locally highly depleted in carbon just outside the CO snow line. Also, additional dust processing could happen in the 20–40 AU between the CO snow line and the inner disk that could influence the carbon abundance. This could be tested by measuring the C/H ratio with James Webb Space Telescope (JWST) or Cryogenic High-Resolution Infrared Echelle Spectrograph (CRIRES+) data in the future. Better resolved radial carbon abundance profiles are necessary to put further constraints on the amount of material that is left behind in the disk to form planets.

There are no signs of dust trapping pressure maxima in the dust continuum images of both DO Tau and DR Tau, often assumed to be a result from planets (Long et al. 2019). However, small substructure will be unresolved in the ~15 AU beam of the continuum observations and the dust is optically thick up to ~25 AU, beyond the CO snow line, which could hide the presence of a gap or dust ring. Future high resolution studies, revolving around substructure near the habitable zone of compact disks, are necessary to determine whether substructure is exclusively present in extended systems or that compact disks host more substructure than we think. For DR Tau, R_gas/R_dust = 3.5 based on the 95% effective radius of the ¹³CO J = 3−2 emission and millimeter continuum, which is close to be a clear sign of dust evolution and radial drift (see Trapman et al. 2019). Considering that the method in Trapman et al. (2019) is based on ¹²CO observations, which in general has a more extended emitting area than the ¹³CO for opacity arguments, and the fact that the spatial extent of the CO isotopologue data is slightly larger than can be explained by the modeling (see Appendix C), these systems likely have advanced dust grain evolution and radial drift. This contrast with the radial carbon abundance profile could potentially be explained by dust traps relatively close to the star, not affecting the radial drift in the outer disk. It is also possible that these systems have had efficient radial drift in the past, which is stopped due to recent formation of one or more planetary cores, but not enough to show up in the submillimeter continuum.

6.2 Evolutionary Trend in Carbon Depletion

Combining all sources where the carbon abundance in the outer disk is determined using multiple carbon-carrying molecules (Fig. 9), we find a potential evolutionary trend in the level of elemental gas-phase carbon depletion, similar to those found in CO (Zhang et al. 2020b; Bergner et al. 2020). Recent measurements of C¹⁸O in ten Class 0/I sources in multiple star forming regions (shown as a yellow star in Fig. 9 Harsono et al. 2014; Zhang et al. 2020b; Bergner et al. 2020) indicate that the youngest sources are not yet depleted in CO, relative to the ISM value, therefore likely not in elemental carbon. As disks age the snow line moves inward and the vertical snow surface decreases, leading to the presence of a massive cold trap in the outer disk. Any CO that mixes through the dust growth layer or midplane will be lost to never return to the gas. This leads to gas-phase carbon depletion as seen in the outer disks of DO Tau and DR Tau by 1−2 Myrs (blue crosses and triangles, respectively, in Fig. 9). As the disks cool further, and CO-rich icy pebbles drift inward over time, the outer disk becomes increasingly depleted, leading to depletion by up to 2 orders of magnitude within 8 Myrs as in DL Tau and TW Hya (blue stars and green squares, respectively, in Fig. 9). DL Tau has a similar age and size as TW Hya, and has also an extended ring system, which explains the striking similarity between the two sources in terms of their outer disk carbon depletion. The only dissimilarity arises from the fact that DL Tau is less efficient in trapping dust, resulting in a much larger difference between the inner and outer disk carbon abundances, which we addressed in Sect. 6.1.1. TW Hya has an inner cavity (Andrews et al. 2016), which might act as an inner pressure bump that prevents the final drift of pebbles and reduces the carbon enhancement in the inner disk. Similarly old disks around hotter stars, for example the ~10 Myr Herbig Ae/Be stars HD 163296 and HD 100546 (red diamonds and triangles, respectively, in Fig. 9), remain only as depleted as the younger T Tauri disks, possibly due to increased heating and ultraviolet radiation from the central star preventing efficient CO ice formation over as large an area as for the T Tauri stars. Larger samples of disks are necessary in order to validate this evolutionary trend.

Calculations of the stellar age can vary considerably depending on the methods used to determine the effective temperature and luminosity to place the sources on evolutionary tracks. The ages for the sources in Taurus that we use here were calculated by McClure (2019) from the Siess et al. (2000) evolutionary tracks using stellar properties determined self-consistently from near-infrared spectra. The age estimates for DR Tau and DO Tau agree well with the young age of 1–2 Myr of the Taurus-Aurigae star forming region (Kenyon & Hartmann 1995), and are close to literature values, 0.37–2.8 Myr and 0.37–0.93 Myr, respectively (see, e.g., Andrews et al. 2013). The age of DL Tau is relatively high for the Taurus-Aurigae star forming region and compared to the literature values of 0.7–2.6 Myr (Andrews et al. 2013) and ${3.5}_{-1.6}^{+2.8}$ Myr (Long et al. 2019). These differences are a result of Andrews et al. (2013) and Long et al. (2019) not taking into account the high veiling of DL Tau (r_J ~2; McClure 2019) when determining the luminosity, namely using a direct fit to the de-reddened photometry, which results in overestimating of the stellar luminosity by a factor of 2 and therefore underestimating of the stellar age. Knowing the precise stellar properties is therefore crucial in our understanding of the evolutionary sequence for the outer regions of the disks.

Recent dynamical modeling of CO depletion processes shows that only a combination of chemical processing, freeze-out and sublimation, vapor diffusion, dust and ice dynamics, and pebble formation and dynamics can explain the level of cold, gaseous CO depletion observed in CO surveys (Krijt et al. 2020; Schwarz et al. 2018). The processes that are involved have in general long timescales, requiring several Myrs to achieve appreciable levels of carbon depletion in the outer regions of a disk. Krijt et al. (2020) find that multiple orders of volatile carbon depletion can only be explained by efficient depletion integrated over at least 3 Myr. Our observations are in line with this modeling, and show that the high level of outer disk carbon depletion of DL Tau and TW Hya is likely due to the fact that more volatile CO has had the chance to freeze out and be locked up in more complex molecules, in large dust grains or transported away from the outer disk by radial drift.

A secondary impacting mechanism on the level of depletion may be the size of the disks, either set by evolutionary or mass arguments. The carbon depletion mechanism is thought to be more efficient in disks with a larger freeze-out zone. HD 100546 and HD 163296 both have no strong carbon depletion, even though the properties of the disk are similar to TW Hya and both sources may be even older than TW Hya (Kama et al. 2016b). However, these two sources are much warmer because of the bright star in the center, and hence they have a small region where CO can freeze out. Similarly, DO Tau and DR Tau are both compact and thus have a relatively small volume beyond the CO snow line where CO freeze-out may occur. Future observations and modeling are needed to get a better understanding of the carbon depletion processes during the different stages of evolution of the disk.

7 Conclusions

In this paper, we present observations and modeling of the main gas-phase carbon carriers in protoplanetary disk atmospheres, C⁰ and CO isotopologues. We observed a sample of seven disks around T Tauri stars in [C I]. Six out of the seven sources have a firm detection in [C I], and in the spectra of four we can constrain the disk component. Three of the sources are modeled with the thermo-chemical disk model DALI to determine the C/H abundance in the outer disk. We can conclude the following:

• We find depletion factors in gas-phase elemental carbon abundance of $5_{-1}^{+2},{17}_{-7}^{+39}$, and ${157}_{-15}^{+17}$ for the outer disks of DR Tau, DO Tau, and DL Tau, respectively.

• DL Tau is much more depleted in carbon in the outer disk than in the inner disk, following a similar pattern observed in other disks. This is likely the result of advanced grain evolution, efficient radial drift, and only minor locking of carbon in the outer disk. We propose that either the substructure formed later than the viscous time or the dust rings observed in the millimeter continuum images of this disk are leaky (i.e., have <100% efficiency in trapping dust).

• DO Tau and DR Tau show, within the systematic error, similar depletion in volatile carbon in the outer disk as in the inner disk, which indicates very efficient dust trapping in the outer disk. Though neither disk shows any dust substructures from current high resolution millimeter continuum observations, our radial carbon abundance profiles suggest the presence of small-scale pressure maxima in the inner disk of these compact dust disks. In addition, these systems could have had efficient planetesimal formation that stopped most of the radially drifting dust grains.

• Combining all the literature volatile carbon depletion factors in the inner and outer disks of T Tauri stars, we find hints for an evolutionary trend in carbon depletion that is consistent with CO isotopologue studies and dynamical models of CO depletion processes that show increasing depletion for older sources.

Expanding the number of points in the radial profiles of carbon depletion will be crucial to understanding where the dust is locked up in the disk, potentially in the form of small planetesimals. Future higher spatial resolution observations of the [C I] J = 1−0 line, in combination with high resolution CO and continuum data, have the potential to be one of the few ways to directly probe the formation of small planetesimals in T Tauri disks.

Acknowledgements.

This paper makes use of the following ALMA data: ADS/JAO.ALMA#2017.1.00857.S, ADS/JAO.ALMA#2016.1.00158.S, ADS/JAO.ALMA#2016.1.01164.S, ALMA is a partner-ship of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. We thank the referee for comments that helped to improve the manuscript.

Appendix A Observing Details

In Table A.1 we present the observing log of the [C I] observations.

Appendix B Pole-on Outflow in [C I]

The redshifted excess emission in the spectrum of DO Tau (see Fig. 1) originates from an off-center, spatially resolved ring of emission, likely emitted from the pole-on outflow that is also observed in CO (Fernández-López et al. 2020). Some of the channels at this velocity are presented in Fig. B.1. The white ellipse in this figure marks the best-fit ellipse to the pole-on outflow observed in ¹²CO (B_maj = 3″, B_mn = 2″, PA = 167°) by Fernández-López et al. (2020). Studying the outflow in detail is beyond the scope of this paper; additional research to the nature of the wind and its impact on the disk dynamics is necessary.

Fig. B.1 Channel maps of the [C I] J = 1−0 emission of the non-Keplerian redshifted component in DO Tau. Channel velocities are relative to the cloud velocity. The white ellipse denotes the best-fit ellipse to the pole-on outflow observed in CO (Bmaj = 3″, Bmin = 2″, PA = 167°; Fernández-López et al. 2020). The arc-like emission is striking similar to the pole-on outflow observed in CO.

Appendix C Surface Density Profile DR Tau

DR Tau is a very complex source with a lot of substructure in the different gas tracers. We explored a range of modeling parameters that fit the SED reasonably well, but all of the models have their own caveats. Here we present some additional models to give a feeling for the uncertainty in carbon abundance of this source-specific model. Unfortunately, there is no archival high resolution ¹²CO data available and the isotopologues ¹³CO and C¹⁸O J = 2−1 are too weak to determine the spatial scale of the gas disk, as the outer regions of the disk are not detected with high confidence. This means that the spatial information we have of DR Tau is limited. In Fig. C.1 A. we present the azimuthally averaged radial profile of the normalized flux of the ¹³CO J = 3−2 emission, together with a small grid of models. In Fig. C.2 we present the corresponding emission spectra of the different transitions. We lowered the carbon abundance to match the modeled [C I] line flux with our observations. Our fiducial model has R_c = 52 AU, similar to the continuum observations, and γ = 1. This model fits the spatial extent of the ¹³CO emission (Fig. C.1 A) reasonably well, but under-predicts the line fluxes of most CO isotopologue transitions (Fig. C.2). Using this model we find volatile outer disks carbon depletion fractions up to 50 to match the [C I] J = 1−0 observations, while a depletion fraction of 5 suffices to fit the CO isotopologue transitions.

Following the procedure described in Miotello et al. (2018), we fitted a power-law surface density profile in the so called pivot-region of the radial profile of the CO isotopologue J = 3−2 emission, where the line is optically thin, but has enough column density that isotope selective processes are not a dominant factor. We present the results of this fit in Fig. C.1 B. For ¹³CO we performed a linear least-squares fit between 75 – 175 AU. For C¹⁸O the pivot region extends from 30 – 150 AU. We find γ = 0.5 and γ = 0.7 for ¹³CO and C¹⁸O, respectively, strengthening the case that γ < 1 for DR Tau.

The best representative model of the line fluxes, presented in the main text, has γ = 0.5 and R_c = 52 AU, with the same mass and geometry as the fiducial model. This model has a smaller spatial extent than the ¹³CO observations (Fig. C.1 A), but fits both the CO isotopologue data as well as the [C I] line fluxes (Fig. 5). We tested a model with γ = 1 and a smaller disk, R_c = 35 AU, that matches the spatial distribution of the CO isotopologue emission of the model with γ = 0.5 and R_c = 52 AU, but this model underproduces the CO isotopologue emission with respect to the [C I] emission, similar as the fiducial model. The fit could be improved in the future by using a different sized gas disk than the dust disk, as follows from observations, which is beyond the scope of this paper.

The ¹³CO J = 2−1 flux is over-predicted by a factor of a few in all models that we analyzed, which is probably due to the low UV-sampling at shorter baselines of the data. We neglect this line further in our modeling. In the main text we used the model that is closest to the continuum observations and reproduced the line emission spectra best. The effect of under-predicting the size of the disk, if any, could be that the disk is more depleted in volatile carbon than assumed in the text, but less than f_depl = 50.

Fig. C.1 Left: Radial profile of the 13CO data (in black) with three models that vary Rc and γ, that are convolved with a Gaussian beam similar to the observations. Parameters that are changed with respect to the fiducial model used in the main text are denoted in the legend. The shaded region marks the 1σ uncertainty on the azimuthal average of the data. Right: Radial profile of the 13CO (orange) and C18O emission (blue) in DR Tau. The dotted lines are a linear fit in log-log space. The fit is only performed in the pivot region (bold dotted line), where the line is optically thin but with enough column density to efficiently self-shield. The shaded regions mark the 1σ uncertainty on the azimuthal average of the data.

Fig. C.2 Line profiles for DR Tau, same as Fig. 7 but with changes in γ = 1 and Rc = 35, as indicated in the figure legend.

Appendix D Channel Maps

In Figs. D.1 – D.7 we present the channel maps of the [C I] observations for each source.

Fig. D.1 Primary-beam-corrected channel maps of the [C I] emission in DR Tau. Solid white contours show the millimeter continuum at 10, 50, and 200 σ.

Fig. D.2 Same as Fig. D.1 but for FZ Tau.

Fig. D.3 Same as Fig. D.1 but for DL Tau.

Fig. D.4 Same as Fig. D.1 but for DO Tau.

Fig. D.5 Same as Fig. D.1 but for WW Cha

Fig. D.6 Same as Fig. D.1 but for FM Cha.

Fig. D.7 Same as Fig. D.1 but for AS 205 A.

References

All Tables

All Figures

Fig. 1 Spatially integrated emission spectra of the [C I] J = 1−0 line at 492.165 GHz for the seven sources in our sample. All the spectra are shifted to their respective local standard of rest velocity (Table 1). The shaded region marks the 1σ uncertainty for each source. The spectrum of AS 205 A is shown in gray with FM Cha as a reference of the typical linewidth for a modestly inclined disk. “N” and “S” depict the Vlsr of AS 205 Ν and AS 205 S, respectively.

In the text

	Fig. 2 Moment zero maps of the [C I] J = 1−0 emission. Solid white contours show the 0.61 mm continuum at 10, 30, 100, and 300σ. The dashed white contour indicates the region used to extract a spectrum. Beam sizes are shown in the lower-left corner of each panel. The 3σ detection limit is illustrated as a white tick mark in the colorbar.
In the text

	Fig. 3 SEDs of the three modeled sources with our model fit in orange. The black dots with error bars represent observations, the dashed black curve denotes the stellar component, and the solid black curve is a combination of the stellar component and the ultraviolet excess due to accreting material as used in the modeling.
In the text

Fig. 4 Modeling results of the total flux as a function of the carbon depletion factor defined as fdepl = (C/H)/(C/H)ISM. Observations are plotted at the interpolated carbon depletion fraction. For DO Tau, we plot two values for the total flux: the star is the corrected flux for absorption with a constant value, and the diamond is only masked for absorption, as discussed in the main text. The derived carbon depletion fraction by interpolation with the model is presented in the top-right corner, with errors given by the measurement uncertainty of the [C I] observations.

In the text

	Fig. 5 Emission line profiles of 13CO J = 2−1,13CO J = 3−2, C18O J = 2−1, C18O J = 3–2, and [C I] J = 1−0 for DR Tau, with the best fitted model overlaid in color. 1σ error bars are indicated in gray for each spectrum. The data and models of the weak CO isotopologue J = 2−1 transitions have been multiplied by a factor of 5 for illustrative purposes.
In the text

	Fig. 6 Total intensity or moment zero map of the 13CO J = 3−2 emission in DR Tau. The black contours show the continuum emission at 0.8 mm. The gray contours show the residuals after subtracting an azimuthally averaged radial disk profile to highlight the infalling stream. The beam size is shown in the lower-left corner.
In the text

Fig. 7 Observations and modeling of DO Tau. Two models are shown: the dashed line is a model with an ISM C/O ratio of 0.47 and a carbon depletion factor of 10, and the solid line is a model with an enhanced C/O ratio of 1 and a volatile carbon depletion fraction of 20. 1σ error bars are indicated in gray for each spectrum. The data and models of the 13CO and C18O J = 2−1 transitions have been multiplied by a factor of 5 and 20, respectively, for illustrative purposes.

In the text

Fig. 8 Observations and modeling of DL Tau. Two models are plotted: the solid line is a model with Rc = 112 AU as observed in the dust and the dashed line as a model with Rc = 65, which fits the spatially resolved [C I] emission best. Further details can be found in the main text. 1σ error bars are indicated in gray for each spectrum. The data and models of the weak CO isotopologue J = 2−1 transitions have been multiplied by a factor of 5 for illustrative purposes.

In the text

Fig. 9 Radial profile of the volatile carbon depletion fraction based on one point inside the dust sublimation rim and a disk-averaged carbon abundance in the outer disk. The outer depletion fractions of the sources in blue are from this work. Inner disk carbon abundances of these three sources are taken from McClure (2019). We added TW Hya in green (Zhang et al. 2019; Kama et al. 2016b; McClure et al. 2020) and one point in the outer disk for ten Class 0/I sources in yellow (Harsono et al. 2014; Zhang et al. 2020b; Bergner et al. 2020) to illustrate the increase in the level of carbon depletion during the system’s evolution. Two Herbig Ae/Be stars, HD 163296 (Kama et al. 2015; Zhang et al. 2019) and HD 100546 (Kama et al. 2015, 2016b), are shown in red to complete the picture of carbon depletion that we have. These sources do not follow the evolutionary trend, because they are much warmer.

In the text

	Fig. B.1 Channel maps of the [C I] J = 1−0 emission of the non-Keplerian redshifted component in DO Tau. Channel velocities are relative to the cloud velocity. The white ellipse denotes the best-fit ellipse to the pole-on outflow observed in CO (Bmaj = 3″, Bmin = 2″, PA = 167°; Fernández-López et al. 2020). The arc-like emission is striking similar to the pole-on outflow observed in CO.
In the text

Fig. C.1 Left: Radial profile of the 13CO data (in black) with three models that vary Rc and γ, that are convolved with a Gaussian beam similar to the observations. Parameters that are changed with respect to the fiducial model used in the main text are denoted in the legend. The shaded region marks the 1σ uncertainty on the azimuthal average of the data. Right: Radial profile of the 13CO (orange) and C18O emission (blue) in DR Tau. The dotted lines are a linear fit in log-log space. The fit is only performed in the pivot region (bold dotted line), where the line is optically thin but with enough column density to efficiently self-shield. The shaded regions mark the 1σ uncertainty on the azimuthal average of the data.

In the text

	Fig. C.2 Line profiles for DR Tau, same as Fig. 7 but with changes in γ = 1 and Rc = 35, as indicated in the figure legend.
In the text

	Fig. D.1 Primary-beam-corrected channel maps of the [C I] emission in DR Tau. Solid white contours show the millimeter continuum at 10, 50, and 200 σ.
In the text

	Fig. D.2 Same as Fig. D.1 but for FZ Tau.
In the text

	Fig. D.3 Same as Fig. D.1 but for DL Tau.
In the text

	Fig. D.4 Same as Fig. D.1 but for DO Tau.
In the text

	Fig. D.5 Same as Fig. D.1 but for WW Cha
In the text

	Fig. D.6 Same as Fig. D.1 but for FM Cha.
In the text

	Fig. D.7 Same as Fig. D.1 but for AS 205 A.
In the text