© ESO 2019

1 Introduction

The distribution of carbon in protoplanetary disks is critical to understanding the organic content of exoplanets and therefore the origins of life. The bulk C/O ratio of planetary atmospheres may be linked back to the solid material accreted in the protoplanetary disk (Madhusudhan 2012). Recent observations of protoplanetary disks demonstrate that CO may be more depleted from the gas phase of the outer tens of AU than can be explained simply by photodissociation and freeze-out. In the three disks with gas masses derived using both CO and HD, the gas masses derived from CO are lower by factors of 5–100 (Bergin et al. 2013; McClure et al. 2016). Full chemical models of spatially resolved ALMA CO isotopolog ue observations in TW Hya confirm that CO is chemically depleted by a factor of 100 in the outer disk, even taking into account photodissociation and freeze-out (Schwarz et al. 2016). Additional chemical modeling of single-dish observations of the atomic [C⁰] 492.161 GHz line, which originates in the photodissociated layer of the disk atmosphere, demonstrates that the depletion seen in CO can be generalized to the bulk carbon budget as well (Kama et al. 2016). These authors find that the outer disk gas phase CO and carbon depletion can be explained by freeze-out of volatile carbon, in the form of CO ice, with subsequent grain surface reactions converting these ices into less volatile ices, like CO₂ or methanol, which keeps the gas phase CO abundance low interior to the CO snowline (Schwarz et al. 2016; van ’t Hoff et al. 2017).

Transformation of solid CO into a more refractory form of carbon on dust grains could enhance the amount of carbon-rich solids available to make planetary cores. If the coagulation of such grains into large bodies beyond the innermost carbon-bearing solid sublimation radius is efficient, then the carbon would be “locked” into the disk and the gas phase carbon abundance would remain low. In contrast, inefficient solid growth would lead to smaller bodies that radially drift towards the star, returning the “missing” carbon to the gas phase in the inner disk when the last icy grains have sublimated. Further modeling of HD, spatially resolved CO isotopolog ues, and resolved continuum observations in TW Hya suggests that up to 2.4 M_⊕ of dust may have grown to larger than centimeter sizes just outside of 5 AU (Zhang et al. 2017), which could permanently sequester volatile carbon. However, that analysis depends on the dust properties assumed for the continuum, and observations at sufficiently high spatial resolution are difficult to obtain at the distances of typical star forming regions, >100 pc. To measure the degree of volatile locking in a larger sample, a better solution is to measure the gas phase carbon abundances inside of the silicate sublimation radius, where the gas is thought to be optically thin for mass accretions rates of less than 10⁻⁷ M_⊙ yr⁻¹, based on models using molecular gas opacities (Muzerolle et al. 2004). For high-mass Herbig AeBe stars, these abundances can be inferred from stellar photosphere abundances because the photospheric convective layer is negligible, so recently accreted gas remains at the surface (Kama et al. 2016). However, the deeper convective layer in low-mass T Tauri stellar atmospheres prevents their stellar abundances to be used in this way to infer inner disk abundances. Disk abundance measurements of T Tauri stars must therefore be made through observations of the accretion columns (Ardila et al. 2013) or the inner disk gas directly.

In 2013, I identified neutral atomic carbon emission in a set of accreting T Tauri stars, as part of a pilot study of the inner dust rim (McClure et al. 2013a). These lines appeared in the strongest accretor in our sample, DR Tau, and grew weaker for the lower mass accretion rate objects. To investigate these lines and the inner dust rim, I observed an enlarged sample of T Tauri stars in Taurus from 0.8 to 2.5 μm. Here I present an initial analysis of the carbon line fluxes in this sample, focusing on the temperatures and densities required to excite them. The sample, extraction of the inner disk excess, and line fluxes are described in Sect. 2, the analysis of the temperatures, densities, carbon ionization fraction, and carbon abundance in Sect. 3, and the implications for the inner disk chemistry and formation of kilometer size bodies are discussed in Sect. 4.

2 Observations and data analysis

The full sample contains 26 well-studied accreting, “classical” T Tauri stars (CTTS) in the nearby Taurus-Auriga molecular cloud complex (~140 pc, Kenyon et al. 1994). The initial selection criteria were that the stars must be (1) single within a detection limit of δmK = 2 at separations greater than 20 mas, with (2) no evidence for gaps or radial structure based on their infrared SEDS, and (3) have Spitzer Infrared Spectrograph (IRS) spectra. Subsequent criteria were that the disks needed to cover a broad range of spectral types and mass accretion rates, based on previous optical publications, e.g. Kenyon & Hartmann (1995). Based on these optical spectral types, I also observed bright non-accreting, “weak-line” T Tauri stars (WTTS) for each target spectral type to act as photospheric templates for the CTTS. Data for ten stars were already published in McClure et al. (2013a) as part of a pilot study.

2.1 Observations

The spectra were obtained with SpeX (Rayner et al. 2003) at the NASA Infrared Telescope Facility (IRTF) on December 1–3, 2010 and January 8–12, 2013. I observed the targets with the short wavelength, cross-dispersed mode (SXD) with the 0′′.3 by 15′′.0 slit (R = 2000) from 0.8 to 2.5 μm with integration times that were selected to produce a signal-to-noise ratio S/N > 100 at H band. The data were obtained with the slit rotated to the parallactic angle and in an ABBA nod pattern. Full details of the data reduction are given in McClure et al. (2013a). In brief, the spectra were reduced with the standard Spextool package (Cushing et al. 2004; Vacca et al. 2003), resulting in sky subtracted, telluric corrected, and flux calibrated spectra.

2.2 Data analysis: inner disk excess extraction

The spectra of the WTTS and CTTS were corrected for radial velocities by shifting them to match the position of the strong 1.31 μm Al I doublet seen in the velocity-corrected IRTF spectral library G, K, and M dwarfs (Rayner et al. 2009). Then the WTTS wereextinction-corrected to match the continuum between the molecular bands shortward of 1.1 μm of the dwarf spectrum of the same spectral type from the IRTF standards library. Using the analysis routines from McClure et al. (2013a) and the additional WTTS spectra, I constructed equivalent width, W_λ, versus SpT trends for the WTTS in this sample and compared them with the IRTF library dwarf and giant equivalent widths. I then computed ratios between equivalent widths of nearby lines to use in determining the veiling-independent spectral type. However, some of the line ratios used in McClure et al. (2013a) showed a local maximum around K0; given the larger range of CTTS and WTTS spectral types (F8-M4), it was necessary to select additional ratios.

Having determined the spectral types of the CTTS, I used the W_λ measurements for individual lines in the CTTS and WTTS templates to determine the amount of veiling in these lines. The veiling measurements and continuum WTTS fluxes were used to compute the extinction for each CTTS as in McClure et al. (2013a). After extinction correction, the WTTS photospheric template is offset from the CTTS spectrum using the average measured veiling between 1.0 and 1.35 μm. The stellar luminosity, L_*, is calculatedfrom the correctly offset WTTS template, using the colors and bolometric corrections given in Kenyon & Hartmann (1995) and new stellar distances from Gaia DR2 (Bailer-Jones et al. 2018) where available or the standard value of 140 pc otherwise (Kenyon et al. 1994). Together with the effective temperatures given by the near-IR spectral types, these luminosities imply a stellar radius R_*. These parameters are compared with the Siess evolutionary tracks (Siess et al. 2000) to derive the stellar mass, M_*. From the excess emission spectrum, I measure the luminosity in the H⁰ Brγ line and use the trend from Muzerolle et al. (2003) to compute the accretion luminosity, L_acc, and mass accretion rate, Ṁ. The filling factor of the accretion shock on the stellar surface can be computed from L_acc, L_*, and T_eff assuming an average temperature of 8000 K in the post-shock region at the stellar surface: $f=\frac{L_{\text{acc}}/L_{*}}{L_{\text{acc}}/L_{*}+{(8000/T_{\text{eff}})}^4}$. The major caveat on this method for calculation of the accretion column filling factor is that the Brγ emission potentially originates in three regions: the accretion shock at the stellar surface, the accretion columns, and the inner disk. A more accurate estimate of the filling factors can be obtained by fitting accretion shock models to the UV excess (Ingleby et al. 2013); however I cannot make this calculation from the present data. The derived stellar and accretion parameters are given in Table 1.

Subtraction of the scaled WTTS photospheric template from the observed extinction-corrected CTTS spectrum reveals the combined excess emission from the accretion shock, inner gaseous disk, and dust sublimation rim.

2.3 Data analysis: carbon line identification and fluxes

The excess, which is a mostly smooth continuum overlaid with emission lines, is fit with a polynomial that is subsequently subtracted off. The continuum subtracted excesses for the five disks with the highest S/N spectra are plotted in Fig. 1. In the excess emission, I detect a series of emission lines at 5–20% of the continuum level. Using the NIST atomic line database (Kramida et al. 2018) and the list of lines in Escalante & Victor (1990) and Haris & Kramida (2017), I identify this emission as carbon recombination lines arising in a predominantly neutral carbon reservoir. The atomic data for all detected transitions are listed in Appendix A. The lines are permitted, with upper states energies of 8–9 eV. They are not seen in the WTTS stars used as photospheric templates, nor in similar IRTF SXD, Magellan FIRE, or VLT X-shooter spectra of transitional disks or pre-transitional disks, whichhave reduced gas and dust masses interior to some disk radius (Vacca & Sandell 2011; Espaillat et al. 2014).

For all of the disks, I calculate an integrated flux for the main line clusters (see Table 2). Since many of the C⁰ lines are blended, this flux is a manual integration of multiple lines rather than separate Gaussian fits. The strongest set of the series is a five-line complex at 1.069 μm; the integrated flux of this set of lines correlates with the flux of Brγ, an accretion indicator, suggesting that it is either associated with the accretion columns/shock or with the innermost gaseous disk inside the dust sublimation rim (Fig. 2). This region is predicted to be optically thin for higher mass Herbig AeBe stars with mass accretion rates of less than ~1 × 10⁻⁷ M_⊙ yr⁻¹ (Muzerolle et al. 2004).

3 Model analysis

To test whether the lines originate in the inner gaseous disk, one must determine how to produce a reservoir of neutral, atomic carbon and how to excite these particular lines. As a first order approximation to the geometrically thin midplane layer in the inner disk, Imodel a simple slab using a combination of a photoionization code and atomic database, as described below.

3.1 Carbon ionization fraction via Cloudy models

To test what fraction of the carbon in the inner disk should be neutral atomic carbon, C⁰, rather than ionized carbon or CO, I used version 17.02 of the photoionization code Cloudy (Ferland et al. 2017) to compute the expected carbon ionization fraction for two sets of slabs: one with a vertical power law in the hydrogen density as a function of distance from the star that roughly approximated the Gaussian density distribution in height expected from a disk in hydrostatic equilibrium, and one with a radially dependent hydrogen density with power −2. The midplane values in each of these slabs at 0.03 AU ranged in a grid from 10¹⁴ < n_H < 10¹⁷ cm⁻³, the range of midplane densities within 0.1 AU in the output of the full disk (D’Alessio et al. 2006) models used by McClure et al. (2016). For the radial slab, I set a stopping criterion of 0.15 AU, as an approximate location of the dust sublimation radius in typical TTS systems (McClure et al. 2013b), and for the vertical slab I set a stopping depth that would slightly overshoot the midplane location at 0.03 AU. Since the region of interest is inside the dust sublimation radius, I did not include a dust contributionin either slab. The physics and chemistry of molecular hydrogen is important for this region, so I enabled the larger H₂ molecule module in Cloudy (Shaw et al. 2005).

The radiation field was a combination of ATLAS9 stellar models with solar metallicity (Castelli & Kurucz 2004) at the stellar effective temperature, T_eff, and the stellar accretion shock temperature of 8000 K, taking the stellar and accretion luminosities found in Sect. 2.2. I also included a contribution to the radiation field from X-rays, assuming a 10 MK bremsstrahlung emission spectrum, with a luminosity of L_X = 10³⁰ erg s⁻¹ between 0.3 and 10 keV (Rab et al. 2018). For the vertical slab, I multiplied the incident flux at the top of the slab by the cosine of the incident angle to the slab normal. I included accretion heating by implementing Cloudy’s α disk option using the stellar mass found in Sect. 2.2.

The resulting physical conditions are shown in Fig. 3 for the stellar parameters of DR Tau. In the vertical slab (left panel), with a starting density of n_H = 10^10.5, the carbon content of the slab is predominantly neutral at densities less than n_H ~ 10¹⁵ cm⁻³. At the midplane, with hydrogen densities n_H ~ 10¹⁶ cm⁻³, the fractional abundance of carbon in CO increases to a limiting value of ~0.5; this limit results from collisional dissociation of CO at such high densities. So conservatively, C⁰ carries at least 50% of the bulk carbon at the midplane. The majority of hydrogen is also in neutral atomic form, with H⁺ making up the next largest contribution. Interestingly, at the midplane density, n_H ~ 10¹⁶ cm⁻³, the fractional molecular hydrogen plateaus around 10⁻⁶, while the ionized hydrogen fraction remains at 10^−3.5, resulting in an electron density larger than n_e ~ 10¹³. For all of the stars in my sample, the relationship between n_e and n_H was a power law form for log (n_H) ≥ 14.5, as seen in Fig. 4: log (n_e) = 1.03 log (n_H)−3.69.

In the radial slab model (Fig. 3, right panel), the temperature and electron density peak at a column density of N_H = 10¹⁹ cm⁻², where most of the stellar radiation is deposited. At midplane densities greater than n_H ~ 10¹⁵ cm⁻³ the slab quickly becomes optically thick to the stellar radiation in the radial direction, reaching an N_H~10²⁵ at R = 3.008R_*. At this column density, even the X-rays in the model are attenuated. Consequently the slab has a short, illuminated “face” that is both radiatively and collisionally excited as well as a larger segment for which the ionization of H⁰ is dominated by charge exchange collisions. The electron temperature exhibits an inversion, with a minimum near N_H~10²⁵ cm⁻² and a plateau of ~7000 K at N_H~10²⁸ cm⁻². As with the vertical slab, the C⁰ fraction is always at least 50% of the total carbon population.

The causes and implications of the high T_e and n_e are discussed in Sect. 4.2. I carry forward the minimum C⁰ fraction (>50%), range of atmosphere and midplane T_e and n_e, and the power law dependence of n_e on n_H for the midplane in order to model the C⁰ line emission in the following section.

Fig. 1 Continuum subtracted excess emission covering the C0 lines for the five disks with the best SN spectrum: FZ Tau, DR Tau, CW Tau, DO Tau, and DL Tau (from top to bottom). Grey vertical lines indicate positions of C0 lines of interest. Note the forbidden lines at 0.982 and 0.985 μm are not detected, in contrast with the permitted lines.

Fig. 2 Integrated luminosity of the three strongest lines in the main 1.07 μm C0 complex plotted against the mass accretion rate, Ṁ, calculated from the H0 Brγ line at 2.166 μm. Uncertainties in the C0 line luminosities are included, along with a linear regression fit to the data (dotted line). The Pearson correlation coefficient is 0.79, indicating a high likelihood of association with Ṁ.

Fig. 3 Two 1D Cloudy slab models of the inner disk around ~0.03 AU. Input stellar parameters are for DR Tau. Left: physical properties are plotted as a function of hydrogen column density for a vertical slab located at 0.03 AU with a power law density structure extending towards the disk midplane. These include: fractional abundance of dominant ionization stages for carbon (top panel) and hydrogen (second panel), hydrogen and electron densities (third panel), and electron temperature (bottom panel). The dashed vertical orange line indicates the location where the slab reaches nH = 1016 cm−3. Middle: analog ous plot, but for a radial slab of initial nH = 1016 cm−3 at the co-rotational radius and a power law dependence on disk radius extending away from the star. Right: same as the middle plot, but plotted against the log arithm of the disk radius instead of nH. Most of the very high density and temperature structure seen in the NH plot occurs at the corotational radius due to the high optical depth.

3.2 Electron densities and temperatures

Under the assumption that the emission arises in an optically thin region, which I verify a posteriori in Appendix C, the line emission depends on the electron density, n_e, and temperature, T_e, as well as the amount of emitting material, $n_{{\text{C}}^0}$. The latter depends on the abundance of carbon with respect to hydrogen, X_C, the fraction of carbon in C⁰, $n_{{\text{C}}^0}$/n_C, and the size of the emitting region. By taking flux ratios between different C⁰ emission lines, one can constrain n_e and T_e independent of X_C or $n_{{\text{C}}^0}$/n_C, assuming that both lines originate in the same region. Using the atomic data and routines provided by the version 8.0.2 of the CHIANTI atomic database (Dere et al. 1997; Landi et al. 2013, and references therein) to calculate the level populations, I determined the line intensity produced by a slab extending from the approximate co-rotational radius at 3 stellar radii, R_*, to an outer radius, R_out, defined below: $$I_{{\lambda }_{ij}}=\frac{h{\nu }_{ij}}{4\pi }{\displaystyle \int n_j}{A}_{ji}\text{d}r=\frac{h{\nu }_{ij}}{4\pi }\frac{A_{ji}}{n_{\text{e}}}{X}_{\text{C}}\frac{n_{{\text{C}}^0}}{n_{\text{C}}}\frac{n_{j,{\text{C}}^0}}{n_{{\text{C}}^0}}{\displaystyle {\int }_{3R_{*}}^{R_{\text{out}}}{n}_{\text{H}}}{n}_{\text{e}}\text{d}r$$(1)

The outer radius of the composite slab is taken from literature estimates based on kinematics of the CO fundamental or H⁰ Brγ interferometry, as described in Appendix B. The effect of optical depth to the stellar radiation is approximated by taking a two-component slab with the same electron density and temperature in each component, but with photoexcitation allowed only interior to 3.008 R_*. To find the best-fitting n_e and T_e independent of emitting area or hydrogen density, there are line pairs that are excellent measures of the density (1.1757/1.0694 μm and 1.0686/1.0694 μm) and temperature (0.9408/1.0694 μm). However, the strongest lines in my sample at ~1.07 and ~1.175 μm are three-line blends that are not spectrally resolved here. Therefore, rather than constraining n_e and T_e from integrated line fluxes, I normalized the continuum subtracted excess spectra to the flux at 1.0694 μm, and performed a least-squares fit of the data to a similarly normalized grid of models.

Figure 5 shows the 1σ contours for the best-fitting range of T_e and n_e for all lines in one disk, DO Tau. The main lines at 1.07 μm are well-fit by a wide range of parameter space, spanning high densities and low temperatures (~ 10¹³ cm⁻³ and 3000 K) to lower densities and higher temperatures (~ 10¹¹ cm⁻³ and 7000 K). However, the next strongest lines at 1.175 μm are only fit by high densities and lower temperatures, while the remaining lines are better fit by lower densities and higher temperatures. In all cases, it is possible to produce a “best fit” that does a reasonable job of fitting most of the lines at an intermediate temperature and density. Figure 6 shows the confidence intervals for the best-fitting slab. Generally, the electron densities are all high, indicating that the lines must arise in the disk midplane layer.

Fig. 4 Electron density versus hydrogen density for vertical slab models of stars in my sample. The orange dotted line represents the power law fit to the region with log nH > 14.5, and labeled equations gives the functional form of this fit.

3.3 Carbon abundances and uncertainties

Using the best-fit T_e and n_e, I fit the absolute flux of each line cluster individually to determine a carbon density. The mean and uncertainty of these carbon densities for the combined line set, $\overline{n_{\text{C}}}$, are listed inTable 3. It is then possible to use the relationship between the midplane n_H and n_e found in the Cloudy models of Sect. 3.1 to compute the n_H implied by these n_e values. I combine this with the observed n_c to obtain X_C, also listed in Table 3. The gas phase carbon abundances are all less than solar, by up to a depletion factor of ~42. I discuss the implications of this result in Sect. 4.1.1.

The abundances and depletion factors in Table 3 are reported with error bars propagated from the 1σ uncertainties on $\overline{n_{\text{C}}}$. Within the uncertainties two disks, DL Tau and FZ Tau, have carbon abundances that are consistent with solar values, while three others, DR Tau, DO Tau, and CW Tau, have robust carbon depletion.

In addition to the scatter in measurements of n_c for each disk, there are three sources of systematic uncertainty in the C/H values found with my method. The first is in the derivation of n_H from the relationship with n_e in the Cloudy slab models. At these high densities, the free electrons are generated mainly through charge exchange reactions, particularly with sulphur-bearing ions (Kingdon & Ferland 1999). The reaction coefficients are accurateto within a factor of 2 to 4 (Kingdon & Ferland 1996), an uncertainty which carries over into the C/H measurements. A second source of uncertainty is the C/O ratio of the gas. Solar C/O ratios produce CO/C_total ratios of at most 0.5, which is what I have taken into account in the final n_C ∕n_H calculation. This upper limit to CO/C holds down to at least C/O ~ 0.1; very low C/O values are unlikely, since the various mechanisms to deplete carbon would also deplete oxygen. One way to resolve this issue is to include fits to the flux of the CO 4.7 μm fundamental line in future numerical simulation efforts; this would directly measure CO/C⁰, particularly if combined with kinematic estimates of the emitting region size (see below).

A third source of uncertainty is the emitting volume, which is used to fit $\overline{n_{\text{C}}}$; the majority of the uncertainty is in the radial extension of the region. Although the vertical structure of the inner disk is important to understand the propagation of stellar radiation, for this sample of disks the C⁰ is well-constrained to originate in a narrow vertical region around the disk midplane, which is much smaller than the height of the dust sublimation rim (Muzerolle et al. 2004). This is because the n_e found from the C⁰ line ratios can only be explained by n_H values so large that the gas has become optically thick to the stellar radiation, which only happens close to the midplane (assuming that the disk is vertically in hydrostatic equilibrium). Therefore the uncertainty in emitting volume can be reduced to the radial extent of the slab. In this case, a smaller emitting region would drive $\overline{n_{\text{C}}}$ to higher values and therefore to a lower degree of carbon depletion. The radial location of the dust sublimation rim in the disk midplane is a hard upper limit for the outer edge of the emitting region, as it is only possible to prevent the reformation of H₂ from H⁰ in the absence of dust grains. To better constrain the outer radius of the emitting region, R_out, future studies should combine kinematically determined outer radii of the C⁰ emission itself with interferometrically determined outer radii of Br γ emission and kinematically determined inner radii of CO from the 4.7 μm fundamental line. According to the Cloudy slab models, these measurements should identify where the disk midplane transitions from mostly molecular to mostly atomic, and therefore where the C⁰ emission can originate.

The latter two uncertainties will be easily resolved with future observations using recently commissioned or upcoming instruments, such as the near- and mid-infrared interferometers VLT-PIONIER and -MATISSE, and high-resolution infrared spectrometers IRTF iSHELL and VLT CRIRES+ to kinematically constrain the emitting region. JWST NIRSPEC will be able to detect these lines at medium resolution even in samples of very faint disks. Additionally, both JWST NIRSPEC and MIRI observationsof the molecular gas outside of the atomic zone will provide robust column densities of all isotopolog ues of the main carbon-bearing molecular species, which will provide a lower-limit check on the total carbon column densities found in the atomic line region. Improvements to the current Cloudy simulations to make them 1+1D would demonstrate better the impact of the vertical disk structure and stellar radiation field, particularly for disks with lower values of n_H, in which photoprocesses play a larger role in the line emission. Systematic uncertainties on the n_H measurements would be improved by calculating or measuring more accurate charge exchange rates, and by separating out the H⁰ contributionfrom the accretion shocks from that of the disk. Finally, benchmarking my simulations with other disk chemical models (e.g. DIANA or DALI) using different reaction networks would also help to understand the impact of these systematic uncertainties on the accuracy of the C/H values calculated here.

Fig. 5 Left: electron temperatures and densities that fit the ratios between the 1.069 μm line and all other lines at a 1σ level (colored contours, wavelengths labeled). The lines are color coded to match the spectra in the right panel. The green fill in between the green contours is an example of how to read the contour plot; the filled region covers the best-fitting set of ne and Te for the 1.07 μm line ratio. Other contours are not filled, but the tickmarks point downhill. The gray contour indicates where the non-detection at 1.0544 μm is matched. The black contours are the 1, 2, and 3σ confidence intervals on the best-fit to the combined lines. Right: continuum subtracted spectra of the six main C0 transition clusters and the non-detections at 0.985 and 1.0544 μm, normalized to the flux at 1.0694 μm, in DO Tau. Model fits for Te = 5750 K and ne = 6.3 × 1012 cm−3 are overplotted. The spectra are color-coded to match the 1σ line ratios regions in the left panel.

4 Discussion

The detection of hot, dense gas inside the dust sublimation rim has important implications for several open topics in the disks and planet formation community, namely (1) the results of grain surface chemistry in the outer disk, (2) observational constraints on the efficiency and composition of planetesimals formed by ~2 Myr, and (3) the structure and physical properties of the little-explored inner gaseous disk.

4.1 Implications of measured carbon abundance

4.1.1 Confirmation of carbon depletion, rather than lower mass

This work is the first evidence for a gas phase carbon depletion, relative to the solar value of C/H = 2.69 × 10⁻⁴, inside the dust sublimation radius of protoplanetary disks. It is also noteworthy that this carbon abundance is inferred relative to hydrogen, via the electron densities, so there is no question of a major degeneracy between carbon abundance and hydrogen mass. Previous work using chemical models to analyze far-infrared and submillimeter HD, CO, and C⁰ line emission also found that carbon is depleted from the gas phase in the outer disk of TW Hya (Bergin et al. 2013; Schwarz et al. 2016; Kama et al. 2016), and potentially two younger disks, DM Tau and GM Aur (McClure et al. 2016). There has been discussion of whether these results are more widely applicable to other systems, given the degeneracy between C/H and the dust/gas ratio in interpreting ALMA observations using chemical models (Miotello et al. 2017). The detection of carbon depletion here in a different sample using independent wavelengths, techniques, and disk radii gives greater credence to the previous results and strengthens the case that, in some systems, a lower gas mass as derived from CO represents a change in the gas phase carbon abundance rather than a lower H₂ mass per se.

4.1.2 Causal mechanisms for C depletion

The outer disk carbon depletion found with ALMA could be accomplished by freeze-out of CO and hydrocarbons onto large dust grains (Yu et al. 2017; Favre et al. 2013). This “missing” carbon should be “released” to the gas phase when the dust grains interior to the snowline corresponding to the sublimation temperature for the solid carbon carrier. Depending on the size of the grains, they may either move approximately with the gas at the disk accretion rate or by radially drifting in faster than the gas accretes. In both cases, the gas measured in the inner disk would have returned to at least the initial C/H abundance (see Fig. 7, middle panel). For the gas in the inner disk to be depleted of carbon, one must avoid accreting the volatile-rich midplane gas resulting from the sublimation of the icy grains. This can be accomplished either by preventing the gas accretion or by preventing the grains from sublimating, e.g. stopping their radial drift (Fig. 7, bottom panel, options B and C).

Midplane volatile C gas accretion can be circumvented by the presence of a deadzone in the midplane, at radii and depths where the disk is not longer thermally or X-ray ionized (Gammie 1996). Volatile-poor gas from the outer disk could then be accreted through the ionized surface layer, while the volatile-rich gas is sequestered in the deadzone. However, more than 50% of the total ISM carbon budget is thought to be contained in refractory-rich volatile grains, which sublimate at temperatures where the MRI should already be active at the disk midplane. Therefore, unless radial drift is inefficient, refractory carbon-rich dust grains would drift through the deadzone and return the inner disk gas phase content to within a factor of 2 of the initial value. In the sample of five modeled disks, four have measured outer dust radii that are substantially less than their outer gas radii (Najita & Bergin 2018, and van Dishoeck, priv. comm., for DR Tau), with the final disk, FZ Tau, lacking the gas disk outer radius measurement. The discrepancy between the dust and gas radii suggests that these disks have experienced efficient radial drift. To deplete their inner disk gas substantially of carbon, some process must be preventing the radial drift of carbon-rich dust grains.

The formation of pressure bumps, which trap dust, or sufficiently large and fast growth can both stop grain from radially drifting. Dust traps typically produce a visible discontinuity in the millimeter grain distribution, e.g. rings or asymmetries. Of these disks, DL Tau shows clear gaps and rings in recent millimeter continuum observations (Long et al. 2018); however, it is also the disk with the least carbon depletion. In contrast, DR Tau has a smooth distribution, but an order of magnitude carbon depletion. In this disk, dust traps cannot be responsible for the observed carbon depletion, leaving grain growth to large bodies as the best explanation for carbon locking in at least a subset of disks.

Fig. 6 Grids of ne and Te for the C0 line ratios of each of the five stars. Black contours are the 1, 2, and 3σ confidence intervals on the best-fitting combination of ne and Te for the detected emission lines, given by the black dot. Gray contours are the 1, 2, and 3σ confidence intervals on the non-detected lines at 0.985 and 1.055 μm. The non-detections preclude the high density, high temperature solutions to the detected line fits seen in, e.g., CW Tau. The bottom right panel is a representative cartoon of the slab model, with the photoexcited cap and optically thick, collisionally excited core.

Fig. 7 Top: schematic of the physical locations in the disk corresponding to different carbon carriers (dust, molecular gas, and atomic gas), with labels for terms used in sketches below. Notional locations of midplane water ice snowline and N2 ice snowline are indicated by dashed blue lines that, for clarity, fade out with radius. Middle: illustration of the carbon release scenario, which returns the inner disk gas carbon abundance to the original value. Bottom: three options for the carbon locking scenario: (a) a deadzone blocking accretion through the midplane, and (b) dust traps and (c) planetesimal formation preventing the radial drift of volatile-rich dust grains. See text in Sect. 4.1.2 for more information.

4.1.3 Properties of kilometer size, carbon-rich planetesimals

I estimate the size of the formed bodies by equating the drift timescale (Birnstiel et al. 2012), $${\tau }_{\text{D}}=\frac{R_{\text{D}}}{u_{\text{D}}}=\frac{R_{\text{D}}{V}_{\text{K}}}{c_{\text{s}}^2\gamma }(\frac{1}{St}+\frac{1}{S{t}^{-1}})$$(2)

to the total age of the system found in Sect. 2.2. The quantities in the expression for the drift timescale are taken from the model of DO Tau published in McClure et al. (2015) and are defined as: R_D, the distance between the formation region to the CO or N₂ snowlines at ~30 AU, u_D, the drift velocity, V_K, the Keplerian velocity, c_s, the sound speed, γ, the absolute value of the power-law index of the gas pressure profile, and St, the Stokes number, which is dependent on the solid body size. If the large body forms at ~31 AU, then it must quickly grow to 3.2 km in order to avoid drifting past the CO snowline within 2 Myr. In contrast, if it forms near 100 AU, it must only be 12.6 m in radius to remain beyond the CO snowline in that time.

Direct comparison of the inner and outer radii carbon abundances in the same protoplantery disk would test the efficiency of the large body formation process. Unfortunately, given the high electron densities required to detect these C⁰ lines, they would not be expected to appear in transitional or pre-transitional systems like TW Hya, DM Tau, or GM Aur, for which the outer carbon depletion factor has been estimated, so it is not yet possible to directly compare the inner and outer disk carbon depletions in the same protoplanetary disk. Future observations of the [C⁰] line at 610 μm and CO ladder in the current sample of sources would make this comparison possible. However, with an assumption about the initial bulk composition of the gas, I can calculate a lower limit to the mass of “missing” carbon that has been locked into kilometer size bodies. Since the mass of hydrogen in the slab models used to fit the line fluxes, M_H, is simply V_slabn_Hm_p, then the mass of gaseous carbon is M_{C, gas} = X_{C, dep}M_H. If the combined gas and ice abundances were solar in the outer disks of the stars in my sample, then the minimum mass of carbon locked in planetesimals is M_{C, solids} = 12 M_H(X_C,solar − X_{C, dep}). These values are listed in Table 3 and are typically around 3 × 10⁻³ M_⊕. Since the disk is not static, this is obviously a gross lower limit.

4.1.4 Inferred origin beyond the N₂ snowline

The carbon masses found above would imply drastically different planetesimal masses, depending on the type of planetesimal into which the carbon was locked. For example, the bulk carbon abundance in planetesimals ranges from C = 31 ± 7 wt.% (Comet 67P, Fray et al. 2017) and C = 3.65 wt.% (CI: chondrites Orgueil/Ivuna, Alexander et al. 2012) to C = 0.0044 wt.% (Earth, Kargel & Lewis 1993). These bodies originate beyond the water ice snowline, between the snowline and the organic sootline, and between the sootline and the inner dust rim where even graphite dust will sublimate, respectively.

The question of which solid-state species is the dominant carrier of this carbon is of interest to discussions of optimizing “stickiness” during collisional grain growth (Musiolik et al. 2016) and exoplanet composition (Madhusudhan 2012). In order to probe the composition, ideally one would measure absolute abundances of silicon as well in the gas phase. Ardila et al. (2013) derive relative luminosities of silicon, nitrogen, and carbon in the pre- and post-shock regions of the accretion streams of T Tauri stars from UV emission from high ionization state lines. Although their technique is different from mine, with a different range of temperatures, ionizations, and physical structures, it is interesting to compare the relative abundances measured for the one star in common between our samples, DR Tau. If one assumes that the variations seen in their values of L_{Si IV}∕L_{C IV} between objects are due to differences in the abundance of Si/C, rather than variations in the local conditions, then one can calculate the gas phase C/Si and N/Si in the accretion stream. These values are a factor of 3 and 10.5 times less than the solar C/Si and N/Si ratios, respectively (see Table 4).

To determine whether this depletion arises predominately from volatile or rocky elements, I use my C/H value to calculate N/H and Si/H. Indeed, silicon is depleted from the accreting gas by a factor of 4.4, and nitrogen by a factor of 45.5, relative to the solar values.

Combining the absolute abundance of carbon, C/H, found here with Si/C and N/C in their paper yields estimates of Si/H_gas and N/H_gas. The amountof “missing” Si and N can be computed with the same method used above for C, resulting in masses on the order of 10⁻⁴ M_⊕. The number density ratios of C/Si and N/Si for the kilometer size planetesimals are then enhanced over solar values by 20 and 26%, respectively (Table 4). This result can directly tie the missing carbon and nitrogen to a particular formation region, i.e. the proto-Kuiper belt. In Pontoppidan et al. (2014), C/Si ratios for meteorites and Earth are on the order of 1 to 10⁻³, respectively, both of which represent a significant solid state depletion of carbon rather than an enhancement. In contrast, comets show C/Si ~10, an enhancement of carbon similar to what I find here.

The contrast is more striking for nitrogen, which is depleted in all bodies, including the two comets in Pontoppidan et al. (2014). However, it is worth noting that the nitrogen content of planetesimals would depend strongly on whether they formed outside of the N₂ snowline. Molecular nitrogen is expected to be the bulk carrier of nitrogen in molecular clouds and the outer regions of disks. If embedded in a water ice matrix, the sublimation temperature of N₂ shifts from 18 to 31 K, placing it at nearly the same location as the CO snowline (T_sub = 34 K in water ice, vs. 21 K otherwise) and the CH₄ snowline at 30 K, or around 30 AU, in a “typical” T Tauri disk (Piso et al. 2016). Icy planetesimals formed in this region are thought to be the parent bodies of Ultracarbonaceous Antarctic Micrometeorites, which contain large fractions of irradiated N- and C-rich polyaromatic organic compounds with N/C values ranging up to solar (Dartois et al. 2017; Engrand et al. 2018). The “locked” solid state values found for carbon and nitrogen in DR Tau are broadly consistent with the formation of kilometer size planetesimals in a proto-Kuiper belt region by ~1 Myr, the age found for the stellar parameters derived here.

It is worth noting that for another disk in which I measure a smaller carbon depletion, DO Tau, the presence of icy planetesimals in the midplane beyond 30 AU is inferred from a tentative detection of water ice at 63 μm with the Herschel Space Observatory (McClure et al. 2015). The crystalline nature of the ice, given the timescale for ice amorphitization at these radii, may arise through collisions between large icy bodies. Unfortunately PACS spectra for the remaining sources in my current sample are either non-existent or too noisy to be conclusive regarding the presence of ice. Further observations with SOFIA HIRMES or a future far-IR space telescope, e.g. OST or SPICA, would be illuminating. High sensitivity observations of ices in edge-on disks, through the upcoming JWST programs IceAge ERS (P.I.: McClure) and MIRI GTO (P.I.: van Dishoeck), will probe the nitrogen to water ratio of the ice inthe upper layers of the outer disk, which will help to determine the abundance and purity of N₂ and hence the probable location of such kilometer size bodies.

4.1.5 Comparison with higher mass stars

An analog ous analysis has been made for the composition of accreting gas in Herbig AeBe star systems. Due to their exterior radiative layer, which does not mix into their convective core, the photospheric abundances of these stars should reflect the composition of recently accreted gas. Kama et al. (2016) demonstrated that Herbig stars do not show any depletion in carbon or oxygen; instead they have on average solar abundances. In contrast, many of the Herbigs with transition disks show refractory depletion in their stellar photosphere, indicating that non-icy kilometer size bodies have formed. These patterns arealso broadly consistent with a scenario in which ices are the prime candidate for the carbon carrier. Since Herbig disks are warmer than T Tauri stars at large radii, due to the higher stellar radiation field, they will have a smaller region in which ices can form after infall from the molecular cloud. Additionally, their larger UV fluxes can photodesorb the ice mantles from their grains over a large radial area, leading to a lower total mass in ices. Since it is less efficient to form planetesimals via core accretion at large radii, and there may be less ice available in Herbig disks, it may be more difficult to form kilometer size planetesimals beyond the CO snowline in these disks than in T Tauri disks. Therefore, more of the carbon locked into the ices would be likely to return to the gas phase inside of the ice snowlines, and the Herbig photospheric C and O abundances would remain solar.

4.2 Inner disk physical conditions and chemistry

Another interesting result from this work is the unusual physical conditions within the dust sublimation radius from the Cloudy models. Specifically, the bulk hydrogen and carbon are simultaneously in the neutral atomic phase, contrary to the expectations of a “typical” photon-dominated or X-ray-dominated region, in which carbon photoionizes before H₂ photodissociates (Maloney et al. 1996). The result found here appears to be a real consequence of the high densities in this region, combined with the lack of solids in the region inside the dust sublimation rim. Specifically, the collisional dissociation rate coefficient for H₂ –H₂ collisions exceeds 10⁻¹⁰ cm³ s⁻¹ when T > ~8000 K and $n_{{\text{H}}_2}$ > 10⁷ cm⁻³ (Roberge & Dalgarno 1982). For higher densities, the temperature threshold is lower. Without dust grains, H₂ cannot reform efficiently and there is a large neutral H reservoir.

Then dissociative collisions between CO and H proceed for T > ~5000 and n_H > 10¹² cm⁻³. Charge transfer at these densities dominates the ionization of H⁺, producing a fraction of ~10⁻⁴ of H⁺ relative to the total hydrogen content. Given the greater abundance of H⁺ than C⁺, it is easy for ionized carbon to find an electron to recombine with, keeping the carbon predominantly neutral despite the large x_e. Recombination of the small amount of H⁺ also produces enough UV emission to electronically excite the smaller amount of H₂, which then heats the gas to ~7000 K through collisional de-excitation. The low fraction of singly ionized carbon is consistent with the non-detection of emission from its two lowest-lying excited states to its ground state (C⁺, 1334.532 and 2324.69 Å) in archival HST COS or STIS spectra of DR Tau. Molecular hydrogen emission has also been found to come from slightly larger disk radii than the CO emission in DR Tau (France et al. 2012).

The temperatures found from the slab models using CHIANTI are lower than those produced by the Cloudy slab models at the disk midplane. However, the Cloudy models are simple 1D models that only allow heating and cooling in a single dimension. At such close distances to the star, a 2D model becomes necessary to fully explain the propagation of stellar irradiation. If the model is allowed to cool radiatively and vertically, the midplane temperature may be lower than in the 1D versions. Alternatively, the charge exchange rate coefficients, which are ultimately responsible for the heating at the midplane, may have a higher degree of uncertainty at such high densities. I will take this into consideration for the next stage of modeling; the temperature and molecular content at the midplane is particularly important for interpreting the observational signatures of molecular emission coming from this region, which may be weaker than previously anticipated.

4.3 Shadowing of the inner dust sublimation rim

The Cloudy models indicate that for full disks with midplane densities greater than n_H ~10¹⁵ cm⁻³, there is a region near the midplane that is optically thick to the stellar radiation. The temperatures found by Cloudy for this hot region are similar to the midplane temperatures found by Muzerolle et al. (2004) for their optically thick inner disk, although their opacity sources mostly molecular. The optically thick region can have a two-fold effect on observations of the inner disk. First, it willcast a shadow on the inner silicate dust rim. Within this shadow, a thin extension of dust and molecules can potentially exist interior to the dust sublimation rim. A radial extension of midplane dust has the potential to produce an additional contribution to the interferometric signal from the inner rim, with a component originating inside of the minimum possible radius defined by gray-opacity grains at the most refractory dust sublimation temperatures (T_sub ~ 1800 K).

Similarly, if the inner disk gas is optically thick to its own radiation, then it could produce continuum emission, which would also produce a component that appears interior to the sublimation radius. Evidence for a contribution interior to the dust rim can be found in many studies, e.g. Eisner et al. (2007), Tannirkulam et al. (2008), Najita et al. (2009), Lazareff et al. (2017). Understanding whether this component is a radial extension of the inner rim is important for theories of planet formation in the inner disk. Boley et al. (2014) appeal to a pressure bump and semi-liquid rocky material at the inner rim in order to build planetary cores rapidly in situ; in their model these cores can eventually become systems of tightly packed inner planets (STIPs) or hot Jupiters without requiring substantial migration (Boley et al. 2016).

Further modeling with a 2D model that includes both the gas and a dust rim is needed to determine the spatial extent of an optically thick region, as well as its temperature, to see if it can explain those objects with interferometric signals inside the dust sublimation rim. In principle for my current model with T Tauri stars, the optically thick gas region would be difficult to distinguish from the star, as it would be unresolved and at a similar temperature.

5 Conclusions

I have analyzed a set of near-infrared spectra of single, accreting T Tauri stars in the Taurus star-forming region. After determining self-consistently the stellar and accretion parameters using updated distance measurements from Gaia, I extracted the excess from the inner disk. In 18 stars of the 26 star sample, I detect recombination lines of C⁰ at a ≥3σ level between 0.9 and 1.5 μm. For the five stars with a complete set of C⁰ lines in this region and a high degree of veiling at ~1.25 μm, I modeled the line emission using two slab models: a Cloudy photoionization model to determine the C⁰ abundance relative to C and an optically thin model with atomic data taken from the CHIANTI database to find the density, temperature, and carbon abundance of the C⁰ line emitting region. All five stars show evidence for the carbon emission arising in a very dense (n_H ~ 10¹⁶ cm⁻³), warm (T_e ~ 4500 K), and moderately ionized (log X_e ~ 3.3) region with a carbon depletion of up to a factor of 42 times with respect to the solar abundance. These parameters match well with predictions of the region interior to the dust sublimation radius, using Cloudy. Analysis of the gas phase carbon depletion at the inner edge of the gas disk implies the following:

• 1.

  Carbon has been depleted from the disk gas onto dust grains and locked into kilometer size bodies somewhere in the disk.

• 2.

  The most likely carrier for the locked carbon is C- and N-rich cometary ices located beyond the N₂ snowline, consistent with formation of a proto-Kuiper belt by 1 Myr.

• 3.

  The disk inside the dust sublimation radius should be largely atomic, with a midplane layer that is optically thick to stellar radiation and shadows the dust rim. This, combined with a potentially optically thick emission component from the4500 K gas, may explain detections of an interferometric component interior to the dust rim in Herbig AeBe stars.

Expanding this model of the inner disk to take into account the proximity to the star, include a 2D treatment of the stellar radiation, and disentangle the hydrogen emission from the inner disk will be crucial to determining the uncertainties on these abundance measurements and understanding the implications of the physical conditions at the dust rim for forming hot planets in situ. Further higher spectral resolution observations of these C lines and other atomic emission from inside the dust sublimation rim have the potential to be one of the few ways to probe the bulk composition of small planetesimals in the disk. I will explore these possibilities in future work (McClure & Dominik, in prep.).

Acknowledgements

This publication is part of a project that has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement ICED No 749864. It also benefited from suggestions by an anonymous referee. Data was taken as a Visiting Astronomer at the Infrared Telescope Facility, which is operated by the University of Hawaii under contract NNH14CK55B with the National Aeronautics and Space Administration.

Appendix A Atomic line data

Most of the C⁰ line complexes are composed of multiple, blended lines. For each individual line in these complexes, Table A.1 lists the wavelengths in vacuum, Einstein A coefficients, level energies, and lower and upper configuration terms distilled from the NIST Atomic Spectra Database (Kramida et al. 2018). The line complexes given around 0.96 and 1.26 μm are very weak and only visually suggested by the spectra upon comparison between targets. Their fluxes are not measured or used for any of the analysis here.

Appendix B Outer radii assumed for slab models

The radial 1D Cloudy model can only cool in one dimension, so it never reaches temperatures where CO would dominate over C⁰. Therefore, it is necessary to infer the outer radius of the C⁰ emission from literature measurements of H⁰, CO, and the dust sublimation rim, as shown in Table B.1.

Two of the disks, DR Tau and DO Tau, have all of these measurements. In H⁰, the outermost radius contributingto Brγ emission hasbeen measured interferometrically, at 0.05 and 0.1 AU (Eisner et al. 2014), respectively. A radius of the peak emission from a hot, broad component of CO has been measured kinematically from CRIRES spectra of the CO fundamental rovibrational band at ~4.7 μm (Banzatti & Pontoppidan 2015), at 0.06 and 0.12 AU, respectively. This peak emission comes from just outside of the Brγ outer radius. In the case of DO Tau, weaker CO emission may extend to an innermost CO radius of 0.03 AU (Salyk et al. 2011). For both DR Tau and DO Tau, the dust sublimation rim has been interferometrically measured at 0.12 and 0.18 AU, respectively, outside of the peak CO emission radius (Eisner et al. 2014). For DO Tau, this rim radius is also consistent with SED modeling of the inner rim excess (McClure et al. 2015). For these two disks, I assume the H⁰ Brγ outer radius as the outer radius for C⁰. Realistically, a large population of H⁰ is necessary in order to produce free electrons from the formation of H⁺, which can then recombine with C⁺ to maintain the neutral carbon population and produce these lines.

The other three disks in the sample lack the complete set of radii observations that DR Tau and DO Tau have. Since the carbon abundance required to fit the observed flux will increase if the emitting region is made smaller, I have chosen to use the smallest radii associated with CW Tau, FZ Tau, and DL Tau in order to present a conservative upper limit on the carbon abundance.

Appendix C Line optical depth verification

The optical depth of the C⁰ lines can be expressed as: $${\tau }_{\text{u},\text{l}}=\frac{n_{\text{u}}}{n_{\text{c}}}{X}_C{N}_{\text{H}}{A}_{\text{ul}}\frac{{\lambda }^3}{8\pi \Delta v}(e^{\frac{h\nu }{kT}}-1),$$(C.1)

where Δv is the width of the line in km s⁻¹, A_ul is the Einstein transition probability in s⁻¹, and N_u, the column density of C⁰ in the upperlevel of the transition, is the product of the first three terms: the fraction of C atoms in the upper level state at a given n_e and T_e, the carbon abundance relative to hydrogen, and the hydrogen column density. In Table C.1 I calculate τ_u,l using the n_e, T_e, n_H, and X_C parameters found for each disk in the main paper, with n_u∕n_c from calculations using the CHIANTI atomic line database, and assuming a Gaussian density profile to estimate N_H. Note that all of the lines are optically thin for each star except for CW Tau; for this star, the density is high enough that the main line in the 0.91 μm complex and the two brightest lines in the 1.07 μm complex have τ_line ~1.

References

All Tables

All Figures

	Fig. 1 Continuum subtracted excess emission covering the C0 lines for the five disks with the best SN spectrum: FZ Tau, DR Tau, CW Tau, DO Tau, and DL Tau (from top to bottom). Grey vertical lines indicate positions of C0 lines of interest. Note the forbidden lines at 0.982 and 0.985 μm are not detected, in contrast with the permitted lines.
In the text

	Fig. 2 Integrated luminosity of the three strongest lines in the main 1.07 μm C0 complex plotted against the mass accretion rate, Ṁ, calculated from the H0 Brγ line at 2.166 μm. Uncertainties in the C0 line luminosities are included, along with a linear regression fit to the data (dotted line). The Pearson correlation coefficient is 0.79, indicating a high likelihood of association with Ṁ.
In the text

Fig. 3 Two 1D Cloudy slab models of the inner disk around ~0.03 AU. Input stellar parameters are for DR Tau. Left: physical properties are plotted as a function of hydrogen column density for a vertical slab located at 0.03 AU with a power law density structure extending towards the disk midplane. These include: fractional abundance of dominant ionization stages for carbon (top panel) and hydrogen (second panel), hydrogen and electron densities (third panel), and electron temperature (bottom panel). The dashed vertical orange line indicates the location where the slab reaches nH = 1016 cm−3. Middle: analog ous plot, but for a radial slab of initial nH = 1016 cm−3 at the co-rotational radius and a power law dependence on disk radius extending away from the star. Right: same as the middle plot, but plotted against the log arithm of the disk radius instead of nH. Most of the very high density and temperature structure seen in the NH plot occurs at the corotational radius due to the high optical depth.

In the text

	Fig. 4 Electron density versus hydrogen density for vertical slab models of stars in my sample. The orange dotted line represents the power law fit to the region with log nH > 14.5, and labeled equations gives the functional form of this fit.
In the text

Fig. 5 Left: electron temperatures and densities that fit the ratios between the 1.069 μm line and all other lines at a 1σ level (colored contours, wavelengths labeled). The lines are color coded to match the spectra in the right panel. The green fill in between the green contours is an example of how to read the contour plot; the filled region covers the best-fitting set of ne and Te for the 1.07 μm line ratio. Other contours are not filled, but the tickmarks point downhill. The gray contour indicates where the non-detection at 1.0544 μm is matched. The black contours are the 1, 2, and 3σ confidence intervals on the best-fit to the combined lines. Right: continuum subtracted spectra of the six main C0 transition clusters and the non-detections at 0.985 and 1.0544 μm, normalized to the flux at 1.0694 μm, in DO Tau. Model fits for Te = 5750 K and ne = 6.3 × 1012 cm−3 are overplotted. The spectra are color-coded to match the 1σ line ratios regions in the left panel.

In the text

Fig. 6 Grids of ne and Te for the C0 line ratios of each of the five stars. Black contours are the 1, 2, and 3σ confidence intervals on the best-fitting combination of ne and Te for the detected emission lines, given by the black dot. Gray contours are the 1, 2, and 3σ confidence intervals on the non-detected lines at 0.985 and 1.055 μm. The non-detections preclude the high density, high temperature solutions to the detected line fits seen in, e.g., CW Tau. The bottom right panel is a representative cartoon of the slab model, with the photoexcited cap and optically thick, collisionally excited core.

In the text

Fig. 7 Top: schematic of the physical locations in the disk corresponding to different carbon carriers (dust, molecular gas, and atomic gas), with labels for terms used in sketches below. Notional locations of midplane water ice snowline and N2 ice snowline are indicated by dashed blue lines that, for clarity, fade out with radius. Middle: illustration of the carbon release scenario, which returns the inner disk gas carbon abundance to the original value. Bottom: three options for the carbon locking scenario: (a) a deadzone blocking accretion through the midplane, and (b) dust traps and (c) planetesimal formation preventing the radial drift of volatile-rich dust grains. See text in Sect. 4.1.2 for more information.

In the text