© ESO 2020

1. Introduction

Relative abundances of common elements, such as, C, O, and Si, are used to characterize both planetary atmospheric compositions (Madhusudhan 2012) and those of solar system solid bodies (Bergin et al. 2015). Despite their high solar abundances and enhancement in gas giant atmospheres, C, O, and N are depleted relative to Si on Earth and all classes of meteorites (Bergin et al. 2015; Lodders 2010). The fraction of each element in rocks (“refractories”) versus ices (“volatiles”) may influence these trends. Ice sublimation depletes solids, and additional mechanisms may exist to destroy C-rich rocky grains in the inner regions of protoplanetary disks (PPDs) while retaining those that are rich in Si (Anderson et al. 2017). The radial transport of fresh solids from the cold outer disk could renew the C in the terrestrial planet forming region of PPDs, and it may be that the low C abundances in Solar System bodies are a consequence of Jupiter segregating the inner solar PPD, rather than a result of destruction by a universal mechanism (Klarmann et al. 2018). However, the midplane distribution of volatile-rich solids in protoplanetary disks has not been observationally well-determined yet due to high optical depths. Moreover, it is unclear where or when these solids start being retained in the disk or “locked” out of the gas accreting onto the central star, which is a key step in the formation of planets (Drake et al. 2005).

The presence of midplane ice reservoirs in the nearby disk TW Hya (Fig. 1) has been inferred from cold gas depletions of CO and water by factors of ∼100 and ∼800, respectively, relative to the values typically measured for the interstellar medium (ISM; Schwarz et al. 2016; Hogerheijde et al. 2011), depending on the disk chemistry (Kamp et al. 2013). These missing gas species are thought to be located beyond “snowlines”, that is to say the locations in PPDs where the midplane temperature is low enough for gas of a given species to freeze out onto dust grains, for example, 4 AU for water and 20 AU for CO in TW Hya (Zhang et al. 2013, 2017). If icy grains grow to at least millimeter sizes, then they no longer vertically circulate into the region where the ice could be photodesorbed by stellar radiation. Dust evolution models predict that these large icy grains will drift radially toward the star faster than the PPD gas, enhancing the gaseous carbon and oxygen abundances inside the CO and water snowlines (Du et al. 2015; Krijt et al. 2018), as is seen for HD 163296 (Zhang et al. 2020). In contrast, TW Hya does not show a strong gaseous volatile molecule enhancement inside its CO snowline at 20 AU (Zhang et al. 2017). Either CO ice is chemically processed into a more rocky carbon carrier or the radial drift of solids is inefficient due to grains collecting in traps (Bosman & Banzatti 2019), which may be the first step to forming planets or planetesimals themselves.

Fig. 1. Schematic of emitting regions for previous TW Hya observations, compared with the region probed in this work. For comparison purposes, observed disk structures in millimeter grains (ALMA, Andrews et al. 2016, VLTI, Menu et al. 2014), micron-sized grains (SPHERE, van Boekel et al. 2017), and CO gas (Banzatti & Pontoppidan 2015; Huang et al. 2018) are shown. The vertical scale is not meaningful. References: (a) Dupree et al. (2012), (b) Brickhouse et al. (2010), (c) Herczeg et al. (2002), (d) France et al. (2012), (e) Bergin et al. (2013), (f) Zhang et al. (2013), Bosman & Banzatti (2019), (g) Banzatti & Pontoppidan (2015), Bosman & Banzatti (2019), and (h) Zhang et al. (2017). The “zone” labels correspond to the regions described in the analytic calculation of Sect. 4.1.

Bulk gas elemental abundances measured inside the dust sublimation rim can can determine whether the icy C and O have returned to the gas phase or are still missing. The former would imply chemical processing of the ice into other ice species with higher sublimation temperatures, while the latter would suggest that the C and O have been locked into midplane ice reservoirs. For stars > 1.4 M_⊙, the inner disk elemental abundances can be inferred from the stellar photospheric abundances, which are due to the formation of a radiative outer envelope that isolates recently accreted material at the stellar surface (Kama et al. 2015; Jermyn & Kama 2018). However, this technique does not work with PPDs around lower mass T Tauri stars, which are the precursors of solar type stars as their outer layers are convective. Interior circulation thus dilutes the compositional signatures of their accreted gas. Recently, we directly measured the inner PPD gas around a set of these lower mass stars through the use of near-infrared atomic C emission lines (McClure 2019). That work inferred carbon locking in these PPDs, but it did not observationally measure abundances of hydrogen or other critical elements, such as O or Si, which could identify whether these elements were locked in rocky versus icy materials.

In this Letter, we model the hydrogen density and atomic abundances for the inner disk of TW Hya from near-infrared atomic emission lines. We then use these abundances to identify the locations where rock- and ice-forming elements are being locked into midplane dust grain reservoirs in this system.

2. Data analysis and physical scenario

To measure the disk abundances inside the dust sublimation radius, we analyzed spectra of TW Hya spanning 0.8−2.5 μm taken in three epochs: 2010a and 2010b (VLT X-shooter) as well as 2013 (Magellan FIRE). The observations and data reduction details are given in Appendix A; in brief, the reduced spectra were corrected for telluric absorption, heliocentric velocities, and convolved down to a uniform spectral resolution of R ∼ 6000. All three epochs show H⁰, He⁰, Ca⁺¹, and O⁰ emission lines (Figs. A.1 and A.2). In the 2013 spectrum, these lines were brighter and showed additional emission from C⁰ and S⁰ lines at 1.07 and 1.045 μm, respectively, as is shown in Figs. A.1 and A.2. We subtracted the continuum under the lines of interest in each spectrum using the automated continuum fitter from McClure et al. (2013a). For the detected Ca⁺¹, O⁰, C⁰, and S⁰ lines, we integrated the lines and list the fluxes in Table A.1. For the nondetected strong lines, we used the continuum root mean square (rms) value to obtain the 3σ upper limits.

The strengthening of all emission lines in 2013 suggests an increase in the density of the accreting, inner disk material for this epoch. We used the H recombination line profiles to confirm this scenario. As is shown in Fig. 2, the continuum-subtracted Paschen γ lines from 2010b and 2013 are both centrally peaked at the stellar rest velocity, with inverse P-Cygni profiles in the wings. However, the maxima and minima of the 2013 profile are stronger. By subtracting the 2010 line profiles from the 2013 spectrum, we removed the contributions of the star and bulk disk to the emission lines. We modeled the residual profile with a four-component Gaussian decomposition. A centrally peaked component is found within the observed velocity resolution uncertainty. Since any stellar chromospheric contribution was removed during the subtraction of the weaker accretion epoch, our fitting results suggest that the centrally peaked component could only come from the inner accretion disk for which the Keplerian motion is in the plane of the sky. In contrast, there is additional inverse P-Cygni emission peaking at −152 km s⁻¹ and absorbing at +114 km s⁻¹. Since TW Hya is face-on (Andrews et al. 2016), we interpret the blue-shifted emission component to be the contribution from magnetospheric accretion columns launching out of the disk plane and the red-shifted absorption component to be the post-accretion shock braking region at the stellar surface. A fourth Gaussian component absorbs at +389 km s⁻¹, which is close to the predicted free-fall velocity at the accretion shock interface and ∼4 times the post-shock velocity, as is expected theoretically (Hartmann et al. 2016). Therefore the shape of the H line profiles support our hypothesis of an increase in the density of the inner disk gas. The physical conditions and elemental abundances in this gas are modeled in the section below.

Fig. 2. Isolation of the hydrogen contribution from the inner disk edge. Panel A: comparison of the Paγ line emission for the 2010b epoch (gray) and 2013 epoch (black) after continuum subtraction. The dashed vertical line indicates the local stellar velocity to which the spectra were corrected. Panel B: Gaussian decomposition fit (gray) of the residual spectrum (black) obtained by subtracting the 2010b epoch spectrum from that of 2013. The different Gaussian components are indicated by colors and include blue-shifted emission (blue: −152 km s−1, FWHM ∼ 172 km s−1), a centrally peaked emission (green: 7 km s−1, FWHM ∼ 126 km s−1), and two redshifted absorption components (dark red: 114 km s−1, FWHM ∼ 137 km s−1; bright red: 389 km s−1, FWHM ∼ 124 km s−1). Panel C: cartoon of the physical scenario revealed by the variability and Gaussian fitting.

3. Modeling

The Gaussian decomposition fit to the H line profiles kinematically isolates the 2013 inner disk hydrogen emission. Using an existing 1D model of the inner disks of T Tauri stars (McClure 2019) driven by the photoionization and chemistry code Cloudy (Ferland et al. 2017), we compared the observed flux ratio between the H Paβ and Paγ emission lines to measure, observationally, the H number density, n_H, in the inner disk. Then we modeled the observed flux ratios between these elemental lines and Paγ, together with the derived n_H, in order to measure total elemental abundances relative to the total amount of H for the inner disk gas.

Our 1D model assumes the inner disk extends from the corotational radius at 0.015 AU out to 0.024 AU. This width was determined from the radial distance travelled on a viscous timescale of ∼2.5 years for the difference between epochs 2013 and 2010b, assuming a “typical” α = 0.01. We took a radially dependent power-law hydrogen density of power −2. We created a composite radiation field to represent the joint contributions from the star and accretion shock. For the star we used ATLAS9 stellar models with a solar metallicity (Castelli & Kurucz 2003) at the stellar effective temperature and luminosity of T_eff = 3850 K with L_* = 0.3 L_⊙. We treated the accretion shock as a combination of ultraviolet (UV) and X-ray emission. The UV emission was approximated by a stellar emission spectrum with a temperature and luminosity of T_shock = 8000 K with L_shock = 0.01 L_⊙. We included a hard X-ray component, assuming a 10 MK bremsstrahlung emission spectrum, with a luminosity of L_X = 10³⁰ erg s⁻¹ between 0.3 keV and 10 keV (Rab et al. 2018). Accretion heating enabled us to use the Hextra keyword, as described in Sect. 11.9.4 of the Cloudy user manual, assuming α = 0.01, the stellar mass of M_* = 0.6 M_⊙ (Sokal et al. 2018), and the radius set to 0.024 AU.

Cloudy includes a chemical network with atomic and molecular opacities that are appropriate for the hot, dense, dust-free region inside the inner dust rim, which enables a radiative transfer calculation in the absence of dust opacity. It is important to note that due to the high densities at this location, the chemistry does not behave as in a typical PDR, that is, the ionization of H is provided by charge exchange collisions, rather than photo- or collisional-ionization. Therefore, C⁰ and H⁰ coexist and they are the main C and H carriers despite the exposure to UV and X-rays or the expected self-shielding by H₂. The inner disk chemistry is discussed extensively in the text of Sect. 3 and Figs. 3 and 4 of McClure (2019), but we also briefly summarize the salient mechanisms, the sources of systematic uncertainty in the final abundances, and a complete list of the model parameters in Appendix B.

A value of log(n_H) = 14.40 cm⁻³ at the inner edge of the disk model is required to fit the observed Paβ/Paγ = 2.3 flux ratio in the inner disk kinematic component from 2010b–2010a. By 2013, the inner disk gas density increased to log (n_H) = 14.74 cm⁻³. This is consistent with the physical scenario derived from the hydrogen line profiles. We then fixed log(n_H) to the 2013 value of 14.74 cm⁻³ and varied the atomic abundances of C, O, S, and Ca. Plots of the model fluxes with the elemental abundance for each element are shown in Fig. A.3. After confirming that the Paβ/Pγ ratio was still valid, we determined the 3σ upper limits to the abundances of N, K, Na, Mg, Fe, Al, and Si, which have undetected strong lines in this region. The resulting abundances are given in Table A.1.

4. Results and discussion

The results of this comprehensive elemental abundance measurement for the inner disk of a T Tauri star contradict the prediction of dust evolution models. All measured elements are depleted with respect to the Sun, with volatiles depleted by up to two orders of magnitude and refractories by up to five orders of magnitude. We note that our C/H and O/H values of 3.39 × 10⁻⁶ and 1.60 × 10⁻⁵, respectively, are consistent within our systematic uncertainties with the following C/H and O/H values for TW Hya that were recently found by Bosman & Banzatti (2019): 3 × 10⁻⁶ and 6 × 10⁻⁶, respectively. Likewise, our N/H upper limit of < 1.6 × 10⁻⁵ is depleted by a factor of ∼4. This value encompasses the depleted N/H of ∼3 × 10⁻⁶ inferred for TW Hya by van ’t Hoff et al. (2017) from ALMA observations of N₂H⁺. It is also consistent with an earlier finding of N depletion on the stellar surface of TW Hya, ranging between factors of 1.5 and 5 (Brickhouse et al. 2010). As is shown in Fig. 3, the accretion of primary atmospheres from this gas would not explain the enhanced giant planet atmospheric volatile abundances seen in our Solar System. However, the abundance pattern seen between the elements in TW Hya’s gas is a close inverse of the CI chondrite abundance pattern, suggesting that the depletion could be caused by the retention of chondrite-like solids in the disk.

Fig. 3. Top: comparison of bulk elemental gas abundances inside TW Hya’s dust sublimation rim (black stars, gray stars with arrows for upper limits), with atmospheric abundances of (Owen & Encrenaz 2003) Jupiter (red circles), Saturn (orange circles), Uranus (light blue circle), and Neptune (purple circle), and CI chondrite solids (Lodders 2010, salmon triangles). All data are plotted with error bars. Abundances are with respect to hydrogen and are normalized to the solar abundances (Asplund et al. 2009), with the solar value indicated by the dotted horizontal line. Bottom: comparison of bulk elemental abundances of solids for: TW Hya (black stars), CI chondrites (Lodders 2010, salmon triangles), comet 67P (bulk, dust/ice 1, teal triangles; Bardyn et al. 2017; Rubin et al. 2019), and UCAMMs (violet triangles; Dartois et al. 2017; Mathurin et al. 2019).

4.1. Outer disk dust traps for C, N, and S

Bosman & Banzatti (2019) suggest that the C depletion in TW Hya’s inner disk could result from either a dust trap at the CO snowline (∼20 AU) or a CO₂ trap close to the H₂O snowline (4 AU). Due to the fact that our C/H measurement is entirely within the dust sublimation rim, we can use previous measurements of CO and a simple analytic treatment (Appendix C) to estimate the location where C is locked out.

First, we divided the disk into three zones (see Fig. 1), one exterior to the CO snowline (Zone 1), a second between the CO snowline and inner rim of the dust (Zone 2), and a third inside the inner rim of the dust (Zone 3). This division closely matches the dust traps found by Bosman & Banzatti (2019), with the dust traps at the end of Zone 1 and in the middle of Zone 2. For each zone, we took a representative value of C/H. In Zone 3, a C/H depletion factor of 79 was taken from the present work. For Zones 1 and 2, we calculated a C/H ratio from the C¹⁸O measurements of Zhang et al. (2019). That work finds a minimum depletion factor for CO/H₂ of ∼78 outside of the CO snowline, while the average depletion factor within the CO snowline is ∼20, which is considered relative to the interstellar medium CO/H₂ value of 10⁻⁴. These values convert into C/H values of 6.4 × 10⁻⁷ and 2.5 × 10⁻⁶, making the gross approximation that all gas phase C is in CO.

Then, we consider that an infinitely thin parcel of gas and dust is transported from the outer edge of Zone 1 through each zone sequentially. The dust is composed of a CO ice component and a “higher temperature” component with carbon carriers that sublimate at higher temperatures than CO, including CO₂, organic residues, or graphite, for example. The dust is either locked into the disk in a given zone or allowed to pass through into the next zone, where a portion of it sublimates. We assume a higher temperature component fraction of 0.83, which is the refractory C fraction observed in comet 67P (Rubin et al. 2019), and we used Eqs. (C.5) and (C.10) to calculate the fraction of the initially available carbon dust that is locked into the disk that the parcel moves through in each zone. We find that 95.9% of the initial carbon dust mass is locked into the disk outside of the CO snowline (Zone 1). Therefore, the combined efficiency of any dust traps outside of that point is also ∼96%. Only 4.1% of the initial dust is transported into Zone 2, with 0.7% of that dust in volatiles that sublimate. Between the CO snowline and the inner dust sublimation rim (Zone 2), an additional 3.06% of the initial dust mass is locked. Only 0.34% of the initial dust crosses the inner dust rim into Zone 3, where it sublimates. In total, only 1.04% of the initial Zone 1 carbon dust mass is returned to the gas phase that is to be accreted onto the central star. The agreement between our C/H abundance in Zone 3 and that of Bosman & Banzatti (2019), which is representative of the gap at 1 AU in Zone 2 and between the H₂O and CO₂ ice snowlines and the silicate dust trap, suggests that 3% of the initial dust locked in Zone 2 is located beyond the 1 AU gap. The efficiency of locking in traps between 1.5 and 20 AU is then 90%.

Most of the uncertainty in the above estimate depends on the higher temperature carbon fraction. Reasonable values range from 0.5, which would lock 99% of the initial carbon mass within Zone 1, at the expense of Zone 2. By decreasing the CO/H₂ depletion factor from 20 to 15, the value observed by Zhang et al. (2019) at 4 AU would reduce the locked dust mass beyond the CO snowline to 94% of the initial dust mass, suggesting that even if CO is converted into CO₂ ice, most of that material is already efficiently trapped in the outer disk rings. Locking solids chemically isolates the inner disk, thus preventing the replenishment of higher temperature C and potentially leading to Earth-like C abundances (Klarmann et al. 2018; Anderson et al. 2017) for forming terrestrial planets.

A quantative analysis of the N depletion validates our C transport model and choice of f_r. The N/H found by van ’t Hoff et al. (2017) is only ∼4% of the solar abundance, that is, the same fraction of N is missing as we find for C beyond the CO snowline in our model using f_r = 0.83. Since N₂ is likely the dominant solid N-carrier in the outer disk and as its midplane snowline at ∼18 K (Piso et al. 2016) should be just outside of the CO snowline, then the N₂ must also be trapped in the three submillimeter rings beyond 30 AU.

In contrast with TW Hya, several other disks show no evidence of C or N gas depletion. HD 163296 has an enhanced C/H in the gas interior to the CO snowline (Zhang et al. 2020) due to radial drift, while IM Lup shows depletions of a factor of 4 at most (Cleeves et al. 2018). For IM Lup, the gas temperature structure puts the N₂ snowline at > 250 AU (Cleeves et al. 2016, Fig. 9), which is well outside of the outermost dust ring in IM Lup at 134 AU (Huang et al. 2018). Therefore, only a little bit of N₂ ice is trapped in this disk and gaseous N is not as depleted. A similar lack of dust traps beyond the N₂ snowline could account for the ISM N abundance found by Anderson et al. (2019) for a set of ∼5−11 Myr disks in upper Sco; however, these disks lack the same high spatial resolution continuum observations that enabled the ring detections in IM Lup and TW Hya. For HD 163296, a ring with a deep gap does sit at 100 AU, just outside the CO snowline at 70 AU, which apparently does not trap dust efficiently (Isella et al. 2018). However, while TW Hya’s outer disk gaps are not as deep as in HD 163296, it does have three of them, which may explain their large cumulative trapping efficiency. Alternatively, if the dust in the rings of TW Hya has grown larger than in HD 163296, it may be filtered out more efficiently (Rice et al. 2006) or planetesimals may have formed, meaning that they would no longer drift efficiently regardless of trapping.

As a further check as to the locking location, in the bottom panel of Fig. 3, we compare the bulk atomic composition of the solids left behind in TW Hya’s disk with the abundances of some primitive Solar System bodies, which roughly originate in each of the three zones, for example, CI chondrites (Lodders 2010, Zone 3), comet 67P (Bardyn et al. 2017; Rubin et al. 2019, Zone 2), and ultra-carbonaceous Antarctic micro-meteorites (Dartois et al. 2017; Mathurin et al. 2019, UCAMMs, Zone 1). The TW Hya solids are all nearly solar; and within the uncertainties, chondritic abundances are generally a good match for the refractory and semi-volatile elements as well as S. However, chondrites underproduce TW Hya’s C, N, and O abundances by an order of magnitude in the case of the former two. Comet 67P and the UCAMMs are more abundant in all three of these elements, but with opposite abundance patterns, that is, there is low N/O in comets and high N/O in UCAMMs. A crude check, fitting a blend of abundances from CI chondrites and either Comet 67P or UCAMMs, or a combination of the two, to TW Hya’s solid abundances, confirms that neither 67P nor UCAMMs on their own provide a good fit to C, N, or O. The former overproduces O and underproduces N, while the latter underproduces N and O, and both fits are poor overall (χ² ∼ 5000−8000). Interestingly, a combination of all three bodies succeeds in fitting C and O, with a χ² ∼ 3, but not N.

There are several potential explanations as to the apparent “excess” of nitrogen in TW Hya’s solids. The UCAMMs are thought to originate from differentiated parent bodies beyond the N₂ snowline in which N₂ was converted to higher temperature residues (Dartois et al. 2017), but they are ultimately recovered in Antarctica. These residues are more stable against sublimation than the original N₂ ice, but is is possible that they encountered some loss on atmospheric entry, artificially lowering the observed abundances. Alternatively, there may be differences between the dust trapping patterns in the Solar System and TW Hya that produce different types of primitive bodies; that is to say, CI chondrites, comet 67P, and UCAMMs may not be an appropriate basis set for TW Hya’s solids. This is generally true to some extent, particularly as we have only used one comet, 67P, to represent the entire cometary class due to a lack of comprehensive elemental abundances for additional comets. Finally, the primitive bodies represent the cumulative elements locked into solids at a particular disk radius over the lifetime of the solar nebula, while TW Hya’s solid abundances represent a snapshot of the material that has been locked out of a moving parcel of gas over multiple radii and a later range of times, from ∼2−8 Myr based on viscous accretion timescales. If the N-rich dust trap formed later in TW Hya’s lifetime than it did in the solar nebula, then the net effect on the cumulative N abundances in the primitive bodies would be diluted.

Unlike N, S can be found in both volatile molecules, for example, H₂S, or minerals that contain more refractory elements, that is, FeS. An analysis of the atmospheres of accreting Herbig AeBe stars with disks suggests that of the total S content, < 1% is in gaseous form, with 89 ± 8% of S is in refractory minerals (Kama et al. 2019) and the remainder in ices. The depletion of S here is consistent with an 89% refractory fraction; if the trap at 2 AU lets through 0.3% of the total initial dust mass, as for the C transport model, then the sublimation of volatile ice with an initial fraction of 11% in that dust can account for all of the gas phase S that we see in the inner disk. The mineral sulphides pass through the 1 AU gap into the inner, refractory disk, which is discussed below.

4.2. Inner disk refractory dust trap

One of the refractory elements, Ca, was firmly detected in the gas phase with a depletion of nearly five orders of magnitude, relative to the solar value. While more depletion of refractories than volatiles is to be expected due to rocky cores in all icy dust grains, five orders of magnitude seems extreme at first glance. However, two of the upper limits to the abundances of the other refractory elements, Mg and Si, are depleted by approximately 10⁻⁴, which roughly supports these types of low values. The latter is qualitatively consistent with previous findings of weaker Si relative to C or N in TW Hya’s accretion streams (Ardila et al. 2013). Our extreme depletions suggest that there is an additional dust trap interior to the sublimation point of most species considered for our higher temperature carbon component, where the bare silicate grains are prevented from sublimating.

An obvious choice for this trap is the dust ring inside of 0.5 AU, detected with ALMA (Andrews et al. 2016). The inner disk temperature and density structures calculated self-consistently in hydrostatic equilibrium normally reach the silicate sublimation temperature of 1400 K at ∼0.1 AU in M0 T Tauri stars with mass accretion rates > 10⁻⁹ M_⊙ yr⁻¹ (e.g., GO Tau in McClure et al. 2013b, Fig. 10). However, using near- and mid-infrared (IR) interferometry, Menu et al. (2014) find the inner rim for TW Hya at 0.3 AU instead, suggesting that the inner edge is associated with a pressure maximum rather than the silicate sublimation front. Blocking accretion or drift of large silicate grains at 0.3 AU would retain the elements Mg, Si, and Ca in the disk. The cumulative amount of the initial dust mass passing through the three submillimeter rings outside of the CO and N₂ snowlines is ∼4%, and the trap near the water snowline only allows 10% of that first 4% of the dust to continue through. To reproduce the fraction of initial Si, Mg, and Ca that returned to the inner disk gas of ∼10⁻⁴ − 2 × 10⁻⁵, only 4.6−0.6% of the dust in the inner 0.3−0.5 AU can pass through the 0.3 AU dust trap and sublimate. This suggests that the trapping efficiency is inversely proportion to the disk radius, with a 96% trapping efficiency cumulatively over the outer three submillimeter rings, 90% at the 2 AU ring, and ∼95−99% at the 0.3 AU ring.

If the refractory C within the 0.3−0.5 AU ring is efficiently removed by photolysis or oxidation, as is suggested by Anderson et al. (2017), then the high trapping efficiency of the 2 AU and outer submillimeter dust traps cuts off C replenishment by drifting solids, as is proposed by Klarmann et al. (2018). Should the innermost, refractory trap of TW Hya eventually form terrestrial planets, they would form from very dry and poor in C and N dust, possibly resulting in Earth-like depletions of these elements.

5. Conclusions

To summarize, we piloted a novel method to determine the bulk atomic composition of unseen solids locked at the midplanes of protoplanetary disks, using the TW Hya disk. The volatile elements in this system are depleted by a factor of ∼100 with respect to the solar value, while the refractories are depleted by up to five orders of magnitude. The high degree of depletion and the difference between the depletion factors of the volatile and refractory elements can be explained if all of the submillimeter rings host dust traps, including the truncated ring from 0.3–0.5 AU. A simple model for the efficiency of dust retention in these traps finds that the inner two dust traps at 0.3 and 2 AU have a high efficiency (≥90% each), while the three submillimeter rings beyond the midplane CO and N₂ snowlines are collectively ∼96% and less individually. Such a high retention efficiency, relative to other disks, may be due to the positioning of the rings relative to individual molecular snowlines as well as the age and maximum grain size in the system. In any event, only a small fraction of the elements in dust end up accreted onto the central star. These results, combined with theoretical modeling of refractory C destruction, suggest that any terrestrial planets forming in the trap within TW Hya’s inner disk from the remaining solids will be relatively “dry” and C-poor, similar to those in our Solar System.

Acknowledgments

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 Skłodowska-Curie grant agreement ICED No. 749864 (M.M.). M.K. was supported by the University of Tartu ASTRA project 2014-2020.4.01.16-0029 KOMEET, financed by the EU European Regional Development Fund. We thank the anonymous referee for comments that significantly improved the manuscript. C. D. acknowledges support by Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO), project no. ALW-GO/15-01.

References

Appendix A: Observation and data reduction details

One spectrum of TW Hya was obtained with the FIRE spectrograph (Simcoe et al. 2013) at Las Campanas (R ∼ 6000, 0.8 < λ < 2.5 μm)(30) on January 2, 2013 (PI McClure). Two 6.0 s exposures in the Fowler 1 read mode were obtained with the 0.6″ slit under 0.45″ seeing at an airmass of 1.009 in excellent conditions. We obtained the usual suite of arc lamp and flat-field calibrations and extracted the FIRE spectrum using the standard FIREhose pipeline, with MCN 7202 as the telluric calibrator star. The other two spectra were obtained with VLT X-shooter (Vernet et al. 2011) (visible: 0.8−1 μm, R ∼ 18 000; near-IR arm: 1−2.5 μm, R ∼ 11 000) on April 7, 2010 and May 3, 2010. We corrected ESO’s Phase 3 reduced data for telluric absorption using the MolecFit software package. Spectra from all three epochs were corrected for their heliocentric velocities, and the X-shooter spectra were convolved down to the resolution of the FIRE spectrum. The processed spectra, before and after local continuum subtraction (as described in Sect. 2), are shown below. A comparable spectrum of hydrogen is shown in the main text (Fig. 2). The integrated line fluxes are given in Table A.1.

Fig. A.1. Top: comparison of the three observed epochs of TW Hya: 2010a (cyan, thin line), 2010b (gray, medium line), and 2013 (black, thick line). The oxygen line (left) is a triple line blend, and the calcium line (right) is blended with an H0 Paschen series line and an absorption artifact from the telluric correction of that line, which we fit out before determining the integrated flux. Bottom: residual spectra after subtraction of epoch 2 from epoch 3.

Fig. A.2. Same as Fig. A.1, but for S (left) and C (right). The integrated flux of the S line is includes the three-line blend, while the C flux includes only the first three lines. The line styles are the same as in Fig. A.1.

Fig. A.3. Upper left corner: model flux ratios between Paβ and Paγ (y-axis) versus nH (x-axis), compared with the observed flux ratio and error bars (horizontal lines). All other panels: model flux ratios of studied atomic lines relative to Paγ (y-axis) versus elemental abundance relative to H in solar units (x-axis) for all elements in this study. The diagonal dashed lines indicate the systematic error in the models due to the charge exchange rate coefficients, as is discussed in Appendix B. Observational flux ratios are indicated (horizontal lines) with error bars for the four elements with > 3σ detections (C, O, S, and Ca), while the 3σ upper limits are indicated for the other elements.

Appendix B: Inner disk chemical model mechanisms and uncertainties

Our 1D inner disk model uses version 17.02 of Cloudy (Ferland et al. 2017). As discussed in Sect. 3.3 of McClure (2019), the high midplane gas densities at 0.024 AU (n_H ∼ 10¹⁵ [cm⁻³], N_H ∼ 10²⁶ [cm⁻²]) lead to a chain of reactions that ultimately produce a mostly neutral gas population, but with a large fraction of free electrons donated mainly by hydrogen itself.

Since the physics and chemistry of molecular hydrogen is critical to understanding the atomic H emission, we enabled the larger H₂ molecule module in Cloudy (Shaw et al. 2005). The collisional dissociation of H₂ is efficient at these high densities; however, the formation of H₂ is hampered due to the absence of dust inside the silicate sublimation rim (at 0.3 AU). Therefore, a large population of H⁰ exists. In turn, at these densities, H⁺¹ is efficiently created through charge exchange reactions, particularly with sulfur as the collision partner. The fraction of hydrogen in H⁺¹ is ∼10^−3.5 and, therefore, H is the major contributor to the electron density, unlike a PDR where C⁺¹ dominates. The recombination to H⁰ produces Lyman photons. However, at these large column densities, the gas is optically thick to UV radiation, so the Lyman photons are absorbed locally by the remaining population of H₂, causing it to enter an excited electronic state. Collisional de-excitation of the H₂ then heats the surrounding gas to ∼6000 K. The net result is that at this location in the disk midplane, there is no dust, little molecular emission, and mostly neutral atomic gas for which emission line ratios directly probe the bulk elemental abundances. We visually summarize this series of processes in Fig. B.1.

The main source of uncertainty in these calculations comes from the charge exchange rate coefficients, which are accurate to within a factor of 2−4 (Kingdon & Ferland 1996, 1999). The coefficients determine the amount of H⁺¹ and, therefore, the electron density of the gas. In turn, the electron density impacts the level populations in the upper state of the transitions giving rise to the observed emission lines. The uncertainty in the rate coefficients leads to a systematic uncertainty of typically 20% on the synthetic line fluxes, as is discussed in Sect. 3.3.4 of Hazy, the Cloudy users’ manual. We have included that uncertainty in the abundances listed in Table A.1. A secondary source of uncertainty in these calculations is the geometry of the emitting region. However, because we are fitting line ratios, we are not overly sensitive to the geometry of the emitting area, assuming that all of the lines originate in the same region.

An additional complication in the N abundance determination in less dense regions would normally be self-shielding of N₂ against photodissociation, which is known to limit the conversion of N₂ into atomic N or other N-bearing molecules more efficiently than for CO. However, as the excitation temperature increases, the self-shielding relative to other modes of shielding (e.g., by H, H₂, or CO) becomes much less efficient (Heays et al. 2014), and a larger column of N₂ would be required in order to self-shield. At the high temperatures and densities in our model (T_e ∼ 6500 K, log(n_H) ∼ 14.75 cm⁻³), Cloudy indeed predicts that atomic N represents at least 99.8% of the total N budget. Critically, even if shielding were efficient at these temperatures, the fact that the UV field is generated locally by an H recombination would effectively prevent the formation of a sufficiently large column to shield, as the H+ and N₂ are intimately mixed. There is observational evidence, in any event, for N₂ photodissociation already occurring further out in the disk. Hily-Blant et al. (2019) detected selective N₂ photodissociation indirectly where the N₂ snowline starts in the upper layers of the disk around 40 AU, through a variation in the ¹⁵N/¹⁴N ratio in ALMA observations of HCN isotopologs. And in our own data, we see the 1.04 μm [N I] quadruplet at all three epochs; however, it is subtracted out of the epoch residual spectrum, suggesting that it originates further out in the disk rather than in the ring of denser material that we modeled moving toward the star in 2013.

As an independent check on the resulting hydrogen density, we ran CHIANTI models of line ratios taken between the C⁰ lines, as in McClure (2019), with additional line ratios between the S⁰ triplet. The range of electron densities and electron temperatures produced by these abundance- and geometry-independent measures for both C and S, combined, are n_e ∼ 2 × 10¹⁰ − 10¹² cm⁻³ at T_e ∼ 7500−5500 K. The best-fitting Cloudy model has a value of n_e ∼ 10¹¹ cm⁻³ at T_e ∼ 6500 K, which is consistent with CHIANTI’s predicted range.

Fig. B.1. Flowchart of the inner disk physical processes leading to the high temperatures and domination of H over H2 and C over CO at the disk midplane, in contrast to the typical photon-dominated region physics seen in molecular clouds and the upper layers of the outer disk.

Appendix C: Carbon locking fractions

As an independent check on the location where carbon is locked out of the accreting gas, we formulated a toy model, pictured in Fig. 1 of the main text. First, we divided the disk into three zones: Zone 1 (between the outer disk edge and the CO snowline), Zone 2 (between the CO snowline and inner rim of the dust), and Zone 3 (interior to the inner rim of the dust). Then we define the model as a parcel of constant density gas and dust that initially sits in the outer disk but is radially transported via accretion processes through the disk to the launching region for the stellar magnetospheric accretion columns. The parcel starts as carbon gas, M_g, 0, with overall solar carbon abundances of X_g, 0 = (C/H)_⊙ = 2.69 × 10⁻⁴ (Asplund et al. 2009), from which in Zone 1 solid carbon condenses out with a mass of

$$\begin{array}{c}\hfill M_{\mathrm{d},1}=M_{\mathrm{g},0}-M_{\mathrm{g},1},\end{array}$$(C.1)

leaving a depleted carbon gas mass of M_g, 1 and a corresponding abundance of X_g, 1 with respect to hydrogen. In the main text, we assume X_g, 1 = 6.4 × 10⁻⁷ (Zhang et al. 2019).

The solid carbon component can be split into volatile ice, as CO, and dust that sublimates at higher temperatures than CO, with the fraction f_r. The higher temperature carbon could be in any of the following materials: CO₂ or other C-rich ices, organic residues, amorphous carbon, or graphite. The grains are assumed to move radially along with the gas, unless they become decoupled from the accretion stream and locked in place in one of the zones. This locking could either be caused by the dust decoupling from the gas and entrapment of larger, millimeter-sized dust grains at pressure maxima or by the formation of bodies that are neither large enough to accrete with the gas nor radially drift. The model is agnostic to the exact locking mechanism. We do not consider complications from the decoupling of the gas and nonlocked dust due to radial drift; we assume that this dust moves with the gas.

We define a fraction, f_L, 1, of the initial dust, M_d, 1, that is locked in some manner in Zone 1, outside of the CO snowline, while the remaining dust exists in small, coupled grains that accrete with the gas across the CO snowline. The successfully transported dust is defined as:

$$\begin{array}{c}\hfill M_{\mathrm{d},\mathrm{d},1}=(1-f_{\mathrm{L},1})M_{\mathrm{d},1}.\end{array}$$(C.2)

At the CO snowline, the volatile dust sublimates, releasing carbon back into the gas phase. Therefore, the new carbon gas mass is

$$\begin{array}{c}\hfill M_{\mathrm{g},2}=M_{\mathrm{g},1}+(1-f_{\mathrm{r}})M_{\mathrm{d},\mathrm{d},1},\end{array}$$(C.3)

which can be rewritten as

$$\begin{array}{c}\hfill M_{\mathrm{g},2}=M_{\mathrm{g},1}+(1-f_{\mathrm{r}})(1-f_{\mathrm{L},1})(M_{\mathrm{g},0}-M_{\mathrm{g},1}),\end{array}$$(C.4)

and it corresponds to a carbon abundance of X_g, 2 with respect to hydrogen. In the main text, we assume X_g, 2 = 2.5 × 10⁻⁶ (Zhang et al. 2019). Solving Eq. (C.4) for the fraction of dust locked in Zone 1, f_L, 1, and rewriting the gas masses of each zone in terms of their observed C/H abundances, X_g, 1 and X_g, 2, yields

$$\begin{array}{c}\hfill f_{\mathrm{L},1}=1-\frac{1}{1-f_{\mathrm{r}}}·\frac{X_{\mathrm{g},2}-X_{\mathrm{g},1}}{X_{\mathrm{g},0}-X_{\mathrm{g},1}},\end{array}$$(C.5)

assuming that no appreciable hydrogen mass is locked into the grains. The interior of the dust mass to the CO snowline is then the fraction of the transported solids with a higher sublimation temperature than CO:

$$\begin{array}{c}\hfill M_{\mathrm{d},2}=f_{\mathrm{r}}{M}_{\mathrm{d},\mathrm{d},1}.\end{array}$$(C.6)

A fraction of the dust mass in Zone 2, f_L, 2, is also eventually locked, either through the process of growth to large grains or to planetesimal sizes. Therefore, the mass of dust that is transported in the parcel across the inner dust rim is the dust mass in Zone 2 less the fraction of dust mass M_d, 2 locked in Zone 2:

$$\begin{array}{c}\hfill M_{\mathrm{d},\mathrm{d},2}=(1-f_{\mathrm{L},2})M_{\mathrm{d},2}.\end{array}$$(C.7)

The mass of the gas phase carbon interior to the innermost dust rim is then the carbon gas mass of Zone 2, plus the additional mass of higher temperature transported dust, all of which now sublimates:

$$\begin{array}{c}\hfill M_{\mathrm{g},3}=M_{\mathrm{g},2}+M_{\mathrm{d},\mathrm{d},2}.\end{array}$$(C.8)

The latter can be rewritten as

$$\begin{array}{c}\hfill M_{\mathrm{g},3}=M_{\mathrm{g},2}+(1-f_{\mathrm{L},2})f_{\mathrm{r}}(1-f_{\mathrm{L},1})(M_{\mathrm{g},0}-M_{\mathrm{g},1}).\end{array}$$(C.9)

Zone 3 corresponds to the carbon abundance found in this work, X_g, 3 = 3.4 × 10⁻⁶. Solving Eq. (C.9) for the fraction of dust mass M_d, 2 locked into Zone 2, f_L, 2, yields

$$\begin{array}{c}\hfill f_{\mathrm{L},2}=1-\frac{1}{(1-f_{\mathrm{L},1})f_{\mathrm{r}}}·\frac{X_{\mathrm{g},3}-X_{\mathrm{g},2}}{X_{\mathrm{g},0}-X_{\mathrm{g},1}}.\end{array}$$(C.10)

This locked fraction, which is in terms of the Zone 2 dust mass, can be converted via the definition of M_d, 2 into a locking fraction in terms of the initial dust mass, M_d, 1. We cannot simultaneously constrain both the locking fraction and the higher temperature carbon fraction, but if we assume that the latter is equivalent to the refractory carbon fraction, then a reasonable value range is from f_r = 0.5 in the ISM to f_r = 0.83 in comet 67P, as is discussed in the main text. Using this range of values and the three observed carbon abundances discussed above produces locking fractions outside (Zone 1) and inside (Zone 2) the CO snowline ranging from 90−99%.

All Tables

All Figures

Fig. 1. Schematic of emitting regions for previous TW Hya observations, compared with the region probed in this work. For comparison purposes, observed disk structures in millimeter grains (ALMA, Andrews et al. 2016, VLTI, Menu et al. 2014), micron-sized grains (SPHERE, van Boekel et al. 2017), and CO gas (Banzatti & Pontoppidan 2015; Huang et al. 2018) are shown. The vertical scale is not meaningful. References: (a) Dupree et al. (2012), (b) Brickhouse et al. (2010), (c) Herczeg et al. (2002), (d) France et al. (2012), (e) Bergin et al. (2013), (f) Zhang et al. (2013), Bosman & Banzatti (2019), (g) Banzatti & Pontoppidan (2015), Bosman & Banzatti (2019), and (h) Zhang et al. (2017). The “zone” labels correspond to the regions described in the analytic calculation of Sect. 4.1.

In the text

Fig. 2. Isolation of the hydrogen contribution from the inner disk edge. Panel A: comparison of the Paγ line emission for the 2010b epoch (gray) and 2013 epoch (black) after continuum subtraction. The dashed vertical line indicates the local stellar velocity to which the spectra were corrected. Panel B: Gaussian decomposition fit (gray) of the residual spectrum (black) obtained by subtracting the 2010b epoch spectrum from that of 2013. The different Gaussian components are indicated by colors and include blue-shifted emission (blue: −152 km s−1, FWHM ∼ 172 km s−1), a centrally peaked emission (green: 7 km s−1, FWHM ∼ 126 km s−1), and two redshifted absorption components (dark red: 114 km s−1, FWHM ∼ 137 km s−1; bright red: 389 km s−1, FWHM ∼ 124 km s−1). Panel C: cartoon of the physical scenario revealed by the variability and Gaussian fitting.

In the text

Fig. 3. Top: comparison of bulk elemental gas abundances inside TW Hya’s dust sublimation rim (black stars, gray stars with arrows for upper limits), with atmospheric abundances of (Owen & Encrenaz 2003) Jupiter (red circles), Saturn (orange circles), Uranus (light blue circle), and Neptune (purple circle), and CI chondrite solids (Lodders 2010, salmon triangles). All data are plotted with error bars. Abundances are with respect to hydrogen and are normalized to the solar abundances (Asplund et al. 2009), with the solar value indicated by the dotted horizontal line. Bottom: comparison of bulk elemental abundances of solids for: TW Hya (black stars), CI chondrites (Lodders 2010, salmon triangles), comet 67P (bulk, dust/ice 1, teal triangles; Bardyn et al. 2017; Rubin et al. 2019), and UCAMMs (violet triangles; Dartois et al. 2017; Mathurin et al. 2019).

In the text

Fig. A.1. Top: comparison of the three observed epochs of TW Hya: 2010a (cyan, thin line), 2010b (gray, medium line), and 2013 (black, thick line). The oxygen line (left) is a triple line blend, and the calcium line (right) is blended with an H0 Paschen series line and an absorption artifact from the telluric correction of that line, which we fit out before determining the integrated flux. Bottom: residual spectra after subtraction of epoch 2 from epoch 3.

In the text

	Fig. A.2. Same as Fig. A.1, but for S (left) and C (right). The integrated flux of the S line is includes the three-line blend, while the C flux includes only the first three lines. The line styles are the same as in Fig. A.1.
In the text

Fig. A.3. Upper left corner: model flux ratios between Paβ and Paγ (y-axis) versus nH (x-axis), compared with the observed flux ratio and error bars (horizontal lines). All other panels: model flux ratios of studied atomic lines relative to Paγ (y-axis) versus elemental abundance relative to H in solar units (x-axis) for all elements in this study. The diagonal dashed lines indicate the systematic error in the models due to the charge exchange rate coefficients, as is discussed in Appendix B. Observational flux ratios are indicated (horizontal lines) with error bars for the four elements with > 3σ detections (C, O, S, and Ca), while the 3σ upper limits are indicated for the other elements.

In the text

	Fig. B.1. Flowchart of the inner disk physical processes leading to the high temperatures and domination of H over H2 and C over CO at the disk midplane, in contrast to the typical photon-dominated region physics seen in molecular clouds and the upper layers of the outer disk.
In the text