Constraining Gas-Phase Carbon, Oxygen, and Nitrogen in the IM Lup Protoplanetary Disk
We present new constraints on gas-phase C, N, and O abundances in the molecular layer of the IM Lup protoplanetary disk. Building on previous physical and chemical modeling of this disk, we use new ALMA observations of CH to constrain the C/O ratio in the molecular layer to be , i.e., higher than the solar value of . We use archival ALMA observations of HCN and HCN to show that no depletion of N is required (assuming an interstellar abundance of per H). These results suggest that an appreciable fraction of O is sequestered in water ice in large grains settled to the disk mid-plane. Similarly, a fraction of the available C is locked up in less volatile molecules. By contrast, N remains largely unprocessed, likely as N. This pattern of depletion suggests the presence of true abundance variations in this disk, and not a simple overall depletion of gas mass. If these results hold more generally, then combined CO, CH, and HCN observations of disks may provide a promising path for constraining gas-phase C/O and N/O during planet-formation. Together, these tracers offer the opportunity to link the volatile compositions of disks to the atmospheres of planets formed from them.
Jansky Fellow of the National Radio Astronomy Observatory
Gas-rich circumstellar disks around young stars provide a window to study the materials that are incorporated into forming planetary systems. Chemistry and dynamics in the disk can shift the balance of gas- versus ice-phase volatiles containing, e.g., carbon, nitrogen, and oxygen. In the core accretion paradigm (Pollack et al., 1996), those materials that end up as rocks or ices are incorporated into the “solid” planetesimals, while the remaining gas can be accreted into natal planets’ atmospheres. Correspondingly, it is essential to understand the form volatiles containing carbon, nitrogen, and oxygen take in the disk spatially and over time, to understand the initial elemental compositions of forming planets (e.g., Öberg et al., 2011a; Piso et al., 2016; Cridland et al., 2016, 2017).
Many processes can alter the gas versus solid abundances in the disk, including snow lines (e.g., Öberg et al., 2011a), chemistry (Bergin et al., 2014; Furuya & Aikawa, 2014; Eistrup et al., 2016; Schwarz et al., 2018), mixing and/or diffusion of gas (e.g., Semenov & Wiebe, 2011; Kama et al., 2016), and redistribution of ices as dust grows and evolves (e.g., Hogerheijde et al., 2011; Piso et al., 2016; Krijt et al., 2016; Öberg & Bergin, 2016). To date, observations of various carbon and oxygen carriers have suggested a substantial “missing” volatile mass within the disk molecular layer as traced by CO in the submillimeter (Favre et al., 2013; Cleeves et al., 2015; Zhang et al., 2017) and HO vapor in the far-infrared with Herschel (Bergin et al., 2010; Hogerheijde et al., 2011; Du et al., 2017), though see also Kamp et al. (2013) regarding model dependencies. For the few disks where estimates for both exist, more water in the disk surface is “missing” than CO when compared to interstellar abundances, and both have abundances lower than what simple desorption (thermal and non-thermal) models predict (e.g., in TW Hya’s disk; Hogerheijde et al., 2011; Favre et al., 2013; Schwarz et al., 2016; Kama et al., 2016; Cleeves et al., 2015).
Absolute elemental abundances are often difficult to estimate due to uncertain hydrogen disk masses. In this context, relative elemental abundances, such as C/O and N/O are promising avenues toward robustly characterizing disk volatile compositions. Recently, Bergin et al. (2016) reported bright hydrocarbon rings of CH and CH in the TW Hya (see also Kastner et al., 2015) and DM Tau disks. Using chemical models, Bergin et al. (2016) found the abundances of these small hydrocarbons were especially sensitive to the gas-phase C/O ratio of the disk. The observations required high C/O values, , to reproduce the observed line intensities (see also Du et al., 2015). Such prospects for measuring C/O in disks are now especially exciting as we enter an era where the elemental compositions, including C/O, of exoplanets’ atmospheres (e.g., Madhusudhan et al., 2011; Kreidberg et al., 2014; Macintosh et al., 2015; Bonnefoy et al., 2016; Lavie et al., 2017).
In contrast to carbon and oxygen abundance estimates, there are few constraints on total nitrogen abundances or N/O ratios in disks owing in large part to the difficulty of observing the likely primary nitrogen carrier, N, (e.g., Schwarz et al., 2016, and references therein). Abundant nitrogen bearing species, such as NH and HCN, are sensitive to other disk parameters than total N abundance, such as temperature structure, carbon abundance, and ionization rate. As such, interpreting these species in the context of bulk nitrogen abundance requires detailed knowledge of the source.
In this work, we constrain the carbon, nitrogen, and oxygen content of the warm molecular layer in the IM Lup protoplanetary disk using results from our previous study of CO and its isotopologues (Cleeves et al., 2016) and new and archival ALMA observations of CH and HCN and HCN. The solar mass star IM Lup harbors a massive gas rich disk, M based upon continuum (SED and resolved millimeter images) and CO multi-isotopologue multi-line data presented in (Cleeves et al., 2016). The source is relatively young at an age of Myr (Mawet et al., 2012).
In Cleeves et al. (2016), we found that IM Lup’s CO is under-abundant by a factor of compared to an interstellar CO abundance of per H based upon the dust-inferred disk mass. For comparison, this younger object appears to be “missing” less CO than the older TW Hya system, whose CO abundance is less abundant than IM Lup (e.g., Favre et al., 2013). However, based on this data alone it is difficult to tell whether the observed “missing” gas-phase CO is a result of missing carbon or missing oxygen or both; or, alternatively, missing total gas mass compared to dust. In the present paper, we explore what C/O and N/O abundance ratios are required in the disk’s warm molecular layer to reproduce the observed CH, HCN, and HCN line intensities.
2.1 ALMA Observations and Data Reduction
IM Lup was observed as part of a survey of nitrogen isotope chemistry and deuterium chemistry presented in Guzmán et al. (2017) and Huang et al. (2017), respectively. Despite this source being one of the most gas-rich disks and bright in HCN (see also Öberg et al., 2011b), the HCN line was not detected at a per channel RMS of 3.6 mJy per beam in 0.5 km s channels (Huang et al., 2017). The observations of HCN and HCN are presented in Huang et al. (2017) (see also Guzmán et al., 2017) and are not reproduced here.
The CH hyperfine complex was observed with ALMA as part of the Cycle 3 2015.1.00964.S program (PI: Öberg) on 01 May 2016 with 41 antennae for 12 minutes on source. The observations were calibrated by ALMA/NAASC staff using J1517-2422 for the bandpass, J1610-3958 for the phase and amplitude, and Titan for the flux calibration. We performed one additional round of phase self-calibration using the continuum within the spectral window containing the line. For the self-calibration, we use a solution interval of 30 seconds and average both polarizations. The continuum is estimated from line-free channels within the same spectral window, which is then subtracted in the -plane to obtain continuum-subtracted images. We detect two blended hyperfine pairs of CH , the and pair and the and pair, see Table 1. While the pairs are blended in velocity space, the inclined orientation of the Keplerian disk causes the emission from each of the pairs to be mostly spatially resolved on the sky. Images are generated using CASA 4.7 (McMullin et al., 2007) and the clean task. Figure 1 shows velocity-averaged channel maps and moment-0 maps tapered to 1” to improve signal to noise. The moment-0 maps show the integrated combined emission from each pair after clipping the individual channels below 1. The disk-integrated spectrum of the observed CH transitions is shown in Figure 2 for each of the hyperfine pairs integrated within a 4” radius circular aperture. Table 1 provides the CH disk integrated fluxes within this aperture, with errors estimated from a 4” circular aperture in line-free portions of the spectrum combined with 10% calibration uncertainty added in quadrature.
|Freq.||Flux Density (Blended)|
|[GHz]||[Jy km s]|
Note. – Uncertainties are quadrature combined RMS scatter and 10% calibration uncertainty.
The following sections describe our methods to constrain carbon, oxygen, and nitrogen abundances in IM Lup’s disk using chemical models to find what elemental abundances best reproduce the CH, HCN, and HCN observations, while also remaining consistent with the previous CO constraints (Cleeves et al., 2016).
3.1 Disk Physical Model
We use the underlying disk structure derived in Cleeves et al. (2016), which we summarize here. The disk surface density follows the self-similarity solutions of Lynden-Bell & Pringle (1974), with an inner power law and outer exponential taper in the disk surface density. The disk is vertically extended and flared, with a gas scale height of 12 AU at 100 AU. The millimeter grains form a thinner layer (25% of the gas scale height) and contain 99% of the total dust mass, M. We furthermore assume small grains follow the gas distribution and allow the millimeter grains to have a separate power law distribution with a cut off at the millimeter emission edge at AU. The small grains provide the greatest surface area per unit mass, and also fill most of the disk volume, and thus are the most important for the grain surface chemistry in the observable layers of the disk ( scale height). In the model, we track the total surface area per unit volume throughout the disk as an input to the chemical calculations.
The gas and dust temperature structures are fixed, where the dust temperature is calculated assuming radiative equilibrium with a stellar K and (Pinte et al., 2008; Cleeves et al., 2016). We note Alcalá et al. (2017) provides updated stellar parameters for IM Lup; however, to within the uncertainties the new values are similar, and so we have chosen to keep these values fixed to more readily compare to previous work. The gas temperature deviates from the dust temperature in the disk upper layers where there is both FUV heating from the central star and the external radiation field, where we fix the external radiation field to be Habing (Haworth et al., 2017; Pinte et al., 2018; Cleeves et al., 2016).
For the stellar high energy radiation field, we use the UV spectral template of TW Hya normalized to the observed Swift UVM2 flux density from IM Lup, and the “quiescent” X-ray template presented in Cleeves et al. (2013) normalized to the observed integrated X-ray luminosity provided in Günther et al. (2010) of erg s. The wavelength dependent radiation transport is calculated for both UV and X-rays using the code of Bethell & Bergin (2011).
3.2 Chemical Modeling Procedure
We calculate the 2D chemical abundances of CH and HCN using a time-evolving gas-grain model first presented in Fogel et al. (2011) and updated and expanded in Cleeves et al. (2014). The chemical model takes into account the spatial changes in dust surface area per unit volume as described in Section 3.2.1 and Appendix A. The simulations are run over 0.5 Myr, corresponding to the lower age limit of IM Lup. However, we reach steady state in the warm molecular layer (the region probed by the observations) well before this time and the results are not affected by this choice. Therefore, we find this assumption does not significantly impact our results. We use a non-deuterated chemical network with 5974 reactions and 638 species. We do not take into account carbon isotope chemistry. For the HCN line radiation transfer, we simply take the model HCN abundance and assume an isotope ratio of (Prantzos et al., 1996).
The two main variables considered in the modeling are:
the initial C, N, and O abundances (Section 3.2.2)
the cosmic ray ionization rate, .
Based on theory (Cleeves et al., 2013) and observations of the TW Hya disk (Cleeves et al., 2015), the typical cosmic ray ionization in disks may be low, s per H. These low rates may be related to magnetized wind modulation or deflection by tangled magnetic fields within the disk (Cleeves et al., 2013). We consider models that have a cosmic ray ionization rate similar to TW Hya and values approximately one order of magnitude higher and one order of magnitude lower. These models correspond to the “T Tauri minimum modulation” (ttm; s per H), “Solar System maximum” (ssx; s per H), and the “Solar System minimum” (ssm; s per H) from Cleeves et al. (2013).
3.2.1 Model updates
For the present study, the chemical code has been updated in three ways. First, we have improved the grain surface chemistry to provide more “realistic” reaction rates. In the past, we have allowed all ice at a given time to participate in grain-surface chemistry, since our code does not treat the multi-layered nature of the ice like the models of, for example, Hasegawa & Herbst (1993), Vasyunin & Herbst (2013), Garrod (2013), and Furuya et al. (2016). As such, without taking layered ices into account, we had implicitly enhanced the efficiency of grain surface chemistry, since more realistically we expect primarily the ice surface to be reactive. To approximate this behavior, we multiply the reaction rates by (number of monolayers), except for atomic H, which we expect to exist mainly as part of the reactive surface. Incorporating this effect is an especially important addition given that most disks appear to have a large amount of settled dust mass, i.e., have a deficit of small grains in their surface, including IM Lup. Essentially, at a fixed volatile abundance, the ice-coating on a given grain can become “thicker” with decreasing total grain-surface area per volume, making less of the ice mobile and reactive.
The second chemical code update is to include a temperature-dependent sticking coefficient as described in He et al. (2016), rather than assuming perfect sticking for all species. We adopt the generic fit from He et al. (2016) Table 1, except for gas-phase water, where we still assume perfect sticking due to its highly polar nature.
The third model update is the incorporation of N self-shielding using the Li et al. (2013) shielding functions and the vertically calculated H and N column densities. Without the latter, we over-predict the atomic N density and under-predict the N density.
3.2.2 Model abundances
Note. – Note: N abundance varied in Section 4.3.
The primary goal of this work is to constrain the total amount of volatile carbon, nitrogen, and oxygen in the upper layers of the disk, which may not have solar or interstellar bulk composition. Furthermore, we wish to remain agnostic to the specific process leading to deviations in the bulk abundances from interstellar values. To accomplish this, we take the simple approach of adjusting the initial chemical abundances in our models to explore a range of possible total C/H, O/H, and N/H abundances to jointly reproduce the ALMA observations of CH and HCN within our Cleeves et al. (2016) model framework.
We have examined the models to ensure that they reach a pseudo-steady state for the species of interest in the warm molecular layer between radii of 20 and 300 AU within a relatively short timescale ( years, or 0.2% of the simulation time). Essentially, the chemistry in this layer quickly re-adjusts to the local conditions and is most sensitive to the bulk C, N, and O content rather than the details of the initial abundance profile. This feature of the warm molecular layer chemistry allows us to constrain the abundances without needing the full (and uncertain) chemical and physical history of the gas. To confirm this behavior, we have tested models where we move the carbon, oxygen, and nitrogen into different carriers, e.g., N versus N versus NH and achieve similar output abundances in the warm molecular layer () to within a few percent.
The abundances of the species that are not varied are listed in Table 2, and are motivated by molecular cloud model abundances (Fogel et al., 2011; Cleeves et al., 2016). Species whose abundances are altered in the initial conditions are in Table 3. The baseline water abundance has been updated to reflect the high-end of interstellar water ice measurements, per H (Boogert et al., 2015), and thus relative depletion factors in water reported here may be higher if the intrinsic interstellar water content is higher, e.g., hidden in larger interstellar grains (i.e. van Dishoeck et al., 2014).
The starting point of our depletion models is our 2016 paper which revealed CO to be under-abundant by a factor of 19 relative to an interstellar CO abundance of per H (Ripple et al., 2013), updated from previously per H. For each abundance mixture, we require there to be sufficient elemental C and O to be able to produce a CO abundance of , but not excess. For example, the carbon abundance can be per H, and the oxygen abundance can be this value or greater (carbon-limited models). Similarly, oxygen can have the same per H abundance, but with equal or excess carbon (oxygen-limited models). For the nitrogen depletion factors, our primary carrier is N, and to simulate nitrogen “removal” we directly reduce the initial N abundance.
Note. – Models with C/O are carbon limited, while models with C/O are oxygen limited. In the latter case, the oxygen is split between CO and HO(gr).
Even though the chemical reprocessing timescales are short in the upper layers () of the disk, we have nonetheless attempted to create “realistic” initial compositions rather than purely atomic. As seen in Table 3, the oxygen and carbon are primarily in HO, CO and C when needed. This choice does not affect the main goals of this study (reproducing the observables), but results in more realistic midplane chemical abundances.
3.3 Model - Observation Comparison
We simulate the CH , HCN , and HCN observations using the non-LTE radiative transfer code LIME v1.8111https://github.com/lime-rt/lime (Brinch & Hogerheijde, 2010). We use the collisional rates input files provided by the Leiden LAMDA database (Schöier et al., 2005). All calculations take into account full non-LTE radiation transfer. The HCN and HCN collision rate data is assumed to be the same and sourced from Green & Thaddeus (1974) and does not include hyperfine structure. The CH collision rate data is from Spielfiedel et al. (2012). The non-LTE analysis is especially important given that the CH is abundant in our models across a wide range of densities, down to cm. We simulate each of the CH and line pairs together since their emission covers the same frequency space, even if the emission is not blended spatially. While the CO gas extends out as far as AU, the HCN and CH and HCN extend to about AU, and so we limit our emission models to this radius to only constrain the abundances where CH is well-detected, and discuss the implications of this in Section 5.2.
The line and millimeter continuum data are modeled in LIME simultaneously as the millimeter continuum opacity is known to be high in this source, especially in the inner disk (Cleeves et al., 2016). The continuum at these wavelengths is is entirely midplane dominated, in a thin vertical layer. We have not attempted to adjust the inner disk opacity as in Cleeves et al. (2016), and note that the main effect would be to reduce the observable line flux from the inner AU, where the CH indeed shows an inner depression (see Figure 1).
The input gas velocities include Keplerian rotation around a solar mass star, thermal broadening, and a fixed turbulent velocity of 100 m s. The final spectral resolution of the CH, HCN, and HCN simulations is 0.28, 0.28, and 0.35 km s, but we simulate the cubes at higher spectral resolution than observed and average down to take into account blurring due to channel averaging. The simulations assume a fixed distance of 161 pc (Gaia Collaboration et al., 2016, DR1), and a fixed inclination of the disk midplane of (Cleeves et al., 2016). The updated DR2 distance is pc, which is sufficiently consistent with the DR1 value, and so we have kept the DR1 distance for ease of comparison.
The models are compared to the ALMA data in the visibility plane, where the vis_sample package (Loomis et al., 2018)222https://github.com/AstroChem/vis_sample is used to sample the LIME channel maps at the same spatial frequencies as were observed to create the model visibilities. To assess the goodness-of-fit for models, we compare the between the observed and simulated continuum subtracted visibilities. The total model grid size is nine values of C/O and three values of for 27 models total.
4.1 Chemical Model Results
Figure 3 presents the 2D chemical abundances for a selection of C/O ratios and the intermediate CR rate value. It is clear that the CH is strongly sensitive to the C/O ratio in the gas, confirming early results of Bergin et al. (2016). The CH abundance is most sensitive for our models below C/O of 1.86, where there is a sharp column density transition straddling C/O of , clearly seen in the two orders of magnitude jump in CH column densities going from C/O of 1.86 to 0.95 in Figure 4. CH is essentially unaffected by the cosmic ray ionization rate for the three model values considered. This lack of dependence occurs because CH is abundant in a layer wrapping around the disk that is UV dominated rather than cosmic ray dominated (see Figure 3 and Bergin et al., 2016).
We can also see that the radial morphology of the column densities changes substantially with C/O (Figure 4). For high values , the CH column density forms a narrow ring or is centrally peaked. At low C/O ratios, the CH column density forms a wide ring that peaks outside of the millimeter dust disk ( AU). This transition occurs once the layer of CH interior to AU becomes thin, with little column density compared to the outer disk. Correspondingly, the radial morphology of CH can change from peaked to ringed even with a uniform C/O ratio, not necessarily requiring (but also not excluding) radial variations in C/O (Bergin et al., 2016).
HCN is sensitive to both the C/O ratio in the gas and the cosmic ray ionization rate, especially for models with C/O . Higher values of C/O and generally increase the HCN abundance. At a given radius, similar column densities are achieved in models with high CR rates and low C/O, and in models with low CR rates and high C/O, though the overall effect on the radial column density profile is different (Figure 4). The 2D abundance morphologies are similar between HCN and CH, though we find the HCN abundance is generally less than that of CH for a given C/O value. One morphological difference is that HCN extends vertically deeper, with a base of AU at AU compared to CH, which disappears below AU for C/O .
4.2 Fits to CH Observations
Based upon the chemical modeling results in Section 4.1, we can use the CH observations to constrain the C/O ratio in the IM Lup disk’s warm molecular layer using the grid of C/O models described in Table 3 and illustrated in Figures 3 and 4. The data and models are compared in the visibility plane with the procedure described in Section 3.3.
Figure 5 shows the - minimum between the data visibilities and our model visibilities varying initial gas composition (C/O) and , where smaller values indicate better fits. Note, adding 1 allows us to plot the difference on log scale. For the CH pair of lines a C/O ratio of is favored, while for the CH pair, C/O of is the best match. Combining both line pairs favors the model. We also find that the CH does not strongly distinguish between cosmic ray ionization rate values, as we expected from the models in Section 4.1. This C/O ratio corresponds to a water ice depletion factor of in the layer CH is present compared to an interstellar water ice abundance of per H (Boogert et al., 2015).
In the same figure, we also show the HCN and HCN results for comparison. C/O values of also fit the HCN data reasonably well, while C/O of predicts detectable HCN image-plane emission, regardless of cosmic ray ionization rate, inconsistent with observations.
4.3 Fits to HCN and HCN Observations
We now consider models with bulk nitrogen depletion to see whether we can arrive at a better fit for HCN at a fixed value for C/O. Taking bulk from Section 4.2, we reduce the nitrogen abundance from the fiducial value of N per H (i.e., N per H, see Table 2). Figure 6 shows the - values for this sub-grid of reduced nitrogen models.
The global best fit is the low value, s, and no nitrogen depletion. If instead we compare models within a fixed CR ionization rate, the intermediate CR ionization rate model favors no N-depletion, while the the high s model, favors anywhere between no and reduction in bulk nitrogen. However, even this factor is small compared to the depletion of CO, or that implied for water and oxygen not in CO based on the CH results, i.e., a depletion.
5.1 Potential mechanisms behind elevated C/O and N/O ratios
The relative differences in the abundances of the key volatile carriers are consistent with a picture of sequential loss of volatiles from the warm molecular layer. Water is the least volatile, freezing out at relatively high temperatures. Therefore it will be most impacted by the evolution of grains, through growth and settling, removing oxygen from the surface over time (e.g., Hogerheijde et al., 2011). CO plays an active role in gas and grain surface chemistry, and as such it is gradually converted to species like CO and CHOH, all of which can freeze out onto grains and then become depleted from the surface via dust settling (Bergin et al., 2014; Schwarz et al., 2018). N, however, is not as chemically active in the disk surface. At cooler temperatures, below the region of CO freeze-out, NH can form directly from N and survive. In the surface, in the presence of CO, this formation pathway is hindered, and most of the surface chemistry requires N be dissociated before forming other nitrogen-bearing species. As such, nitrogen will be less chemically “processed” than the key oxygen and carbon carriers and is expected to stay in its volatile N form given its relatively low desorption temperature (Öberg et al., 2005). Consequently, sequestration into ice and subsequent settling should theoretically be less effective for nitrogen-bearing species, consistent with our results from the ALMA observations.
5.2 Gas-phase C/O and N/O constraints
The model grid presented in Section 3.2.2 focuses on the total elemental abundances in both the gas and solid phases, as these are most relevant for the chemical modeling. The relevant values for comparing to observations of exoplanet atmospheres formed via core accretion are the C/O and N/O ratios specifically in the gas phase. Both CO and N are primarily in the gas phase where our observations are most sensitive (i.e., the region where CH emits, ), while HO is primarily ice with some gas-phase HO from UV photo-desorption (also producing OH).
Figure 7 plots the 2D distributions of the gas-phase C/O and N/O ratios for our best fit model with an intermediate CR ionization rate ( s) and no nitrogen depletion. Whether we use an intermediate or one order of magnitude lower CR ionization rate does not significantly impact these results. At the surface, all oxygen that starts in HO ice is quickly dissociated to gas-phase atomic oxygen, such that the gas C/O ratio is equal to the bulk ratio of . Where water is both frozen-out and shielded from UV deeper in the disk, CO becomes the primary carbon and oxygen carrier in the gas phase, resulting in C/O . The layered C/O structure is visible here, where C/O in the gas is above normalized heights of and below this layer until CO starts to freeze-out, at .
The N/O ratio on the other hand is throughout the disk atmosphere, and when most of the nitrogen is in the gas and CO is the primary oxygen carrier. In the layer where CO begins to freeze-out, N is still in the gas due to its slightly lower binding energy (Fayolle et al., 2016), and the gas-phase N/O ratio can be .
At high velocities we find a discrepancy between the data and modeled CH, which can be seen in Figure 8. The high velocity emission of CH is weaker in the models than in the data. Given the signal to noise of the CH data, we did not attempt to radially or vertically vary the carbon, oxygen, and nitrogen abundances spatially in our fitting; however, from other sources with varied ring-like morphologies in CH (Bergin et al., 2016), there is support for local variations in bulk C/O, N/O, etc. In this instance, the bright high velocity CH suggests that the inner disk has more excess carbon than the outer disk, which may point to slower sequestration of carbon-bearing volatiles in this region due to, e.g., reduced freeze-out in the warm inner disk.
We also examined models where we allowed the disk to extend out to 1000 AU, beyond the AU radial region where CH and HCN are observed. Since the model CH abundance continues to rise in the outer disk for the case (see Figure 3), we predict some emission beyond 500 AU, which was not seen. There may be a change in the gas surface density profile near this radius (Cleeves et al., 2016; Avenhaus et al., 2018), or the disk may have a lower C/O ratio in this region, i.e., may have less efficient volatile sequestration at these very low disk densities. Higher signal to noise, resolved observations of CH toward this source will help disentangle these scenarios.
Using detailed physical and chemical models constrained by CO and dust observations, along with CH and HCN data from ALMA, we constrain the C/O and N/O ratios in the molecular layer of the IM Lup disk. Our observations trace the properties of the disk where CH and HCN emit, primarily at normalized heights above . In these layers, the CH observations favor a super-solar elemental C/O ratio of . This high ratio is consistent with preferential loss of water ice from the surface, e.g., due to sequestration by an evolving and growing population of grains. We do not need to sequester any nitrogen for our best fit model, such that the N/O ratio is also super-solar, . The gas phase values of C/O and N/O vary spatially depending on the degree of UV shielding at a given location (Section 5.2).
While grain sequestration of ices tends to remove volatiles from the surface, it implicitly carries them to the midplane and into the inner disk through settling and radial drift (e.g. Öberg & Bergin, 2016; Piso et al., 2016; Krijt et al., 2016). If this interpretation is correct, this process would result in a large enhancement over interstellar values of water ice in the inner disk midplane, a mild enhancement of carbon-bearing ice, and relatively little nitrogen ice transport in the solid phase. Correspondingly, we expect the C/O and N/O ratios in midplane solids to both be lower than solar, with N/O much less than solar in the disk midplane.
Where settled ices eventually end up radially is still an open question. Radial drift of solids is thought to be quite efficient (see review of Testi et al., 2014), but this process may be slowed by the emergence of pressure variations in the disk that can effectively trap solids (e.g., Weidenschilling, 1980). Now with ALMA, ringed radial structures are being observed, and may even be common (e.g., ALMA Partnership et al., 2015; Andrews et al., 2016; Isella et al., 2016; Loomis et al., 2017; Huang et al., 2018; Fedele et al., 2018). In the absence of pressure traps, these grains should travel inward, thermally desorb, and enhance primarily oxygen, followed by carbon, and relatively little nitrogen. Salyk et al. (2011) indeed found very low N/O ratios, , in the inner disk with Spitzer; however, the nitrogen “correction factors” in the inner disk gas from HCN to total N are uncertain and require additional chemical modeling to constrain.
These results have interesting consequences for the debate regarding abundance measurements in disks (see summary of Bergin & Williams, 2017). At the low gas temperatures typical of disks, H does not emit appreciably (Bergin et al., 2013; McClure et al., 2016; Bergin & Williams, 2017). As a result, other uncertain mass tracers are typically used, such as the total dust mass multiplied by a conversion factor, or even optically thin CO emission itself. From our modeling, we find that we do not need to deplete nitrogen significantly, with a factor of difference between the CO depletion factor and that for nitrogen. These results would point to abundance variations between these species rather than an overall under-accounting of disk mass, which would impact all volatile abundances at similar if not equal levels. We of course cannot rule out with these data alone some missing gas mass, since for the higher CR models, a small amount of nitrogen depletion (a factor of a few) is allowed, even though these are not the global best fits. However, missing gas alone cannot explain the higher degree of depletion needed for CO and water ice. In the future, additional observations of nitrogen-bearing molecules like HCN may help to break this CO mass / gas mass degeneracy, where in larger disk surveys CO masses appear globally low relative to the dust masses scaled by the interstellar gas-to-dust ratio (Ansdell et al., 2016; Miotello et al., 2016; Long et al., 2017).
Going forward, these results show that readily observable molecular tracers like CO, CH, and HCN, and their isotopologues combined with astrochemical models can be used to constrain the gas-phase C/O and N/O ratios within planet-forming disks. Such measurements may furthermore shed light on the inferred volatile composition of the ices, once we better understand the mechanism(s) of volatile loss from the warm molecular layer (e.g., Furuya & Aikawa, 2014; Kama et al., 2016). Future high sensitivity, resolved observations of many sources may further help shed light on “typical” gas-phase C/O and N/O ratios, how they spatially vary, and how these ratios may vary with time as planets are forming in the disk. Together, these can one day be compared with C/O measurements of exoplanet atmospheres, to eventually help unravel their formation locations and histories, where observations have already shown a wide range of exoplanet C/O values, even within a single planetary system (Bonnefoy et al., 2016; Lavie et al., 2017).
- Alcalá et al. (2017) Alcalá, J. M., Manara, C. F., Natta, A., et al. 2017, A&A, 600, A20, doi: 10.1051/0004-6361/201629929
- ALMA Partnership et al. (2015) ALMA Partnership, Brogan, C. L., Pérez, L. M., et al. 2015, ApJ, 808, L3, doi: 10.1088/2041-8205/808/1/L3
- Andrews et al. (2016) Andrews, S. M., Wilner, D. J., Zhu, Z., et al. 2016, ApJ, 820, L40, doi: 10.3847/2041-8205/820/2/L40
- Ansdell et al. (2016) Ansdell, M., Williams, J. P., van der Marel, N., et al. 2016, ApJ, 828, 46, doi: 10.3847/0004-637X/828/1/46
- Avenhaus et al. (2018) Avenhaus, H., Quanz, S. P., Garufi, A., et al. 2018, ApJ, 863, 44, doi: 10.3847/1538-4357/aab846
- Bergin et al. (2014) Bergin, E. A., Cleeves, L. I., Crockett, N., & Blake, G. A. 2014, Faraday Discussions, 168, 61. https://arxiv.org/abs/1405.7394
- Bergin et al. (2016) Bergin, E. A., Du, F., Cleeves, L. I., et al. 2016, ApJ, 831, 101, doi: 10.3847/0004-637X/831/1/101
- Bergin & Williams (2017) Bergin, E. A., & Williams, J. P. 2017, in Astrophysics and Space Science Library, Vol. 445, Astrophysics and Space Science Library, ed. M. Pessah & O. Gressel, 1
- Bergin et al. (2010) Bergin, E. A., Hogerheijde, M. R., Brinch, C., et al. 2010, A&A, 521, L33, doi: 10.1051/0004-6361/201015104
- Bergin et al. (2013) Bergin, E. A., Cleeves, L. I., Gorti, U., et al. 2013, Nature, 493, 644, doi: 10.1038/nature11805
- Bethell & Bergin (2011) Bethell, T. J., & Bergin, E. A. 2011, ApJ, 739, 78, doi: 10.1088/0004-637X/739/2/78
- Bonnefoy et al. (2016) Bonnefoy, M., Zurlo, A., Baudino, J. L., et al. 2016, A&A, 587, A58, doi: 10.1051/0004-6361/201526906
- Boogert et al. (2015) Boogert, A. C. A., Gerakines, P. A., & Whittet, D. C. B. 2015, ARA&A, 53, 541, doi: 10.1146/annurev-astro-082214-122348
- Brinch & Hogerheijde (2010) Brinch, C., & Hogerheijde, M. R. 2010, A&A, 523, A25, doi: 10.1051/0004-6361/201015333
- Cleeves et al. (2013) Cleeves, L. I., Adams, F. C., & Bergin, E. A. 2013, ApJ, 772, 5, doi: 10.1088/0004-637X/772/1/5
- Cleeves et al. (2014) Cleeves, L. I., Bergin, E. A., & Adams, F. C. 2014, ApJ, 794, 123, doi: 10.1088/0004-637X/794/2/123
- Cleeves et al. (2015) Cleeves, L. I., Bergin, E. A., Qi, C., Adams, F. C., & Öberg, K. I. 2015, ApJ, 799, 204, doi: 10.1088/0004-637X/799/2/204
- Cleeves et al. (2016) Cleeves, L. I., Öberg, K. I., Wilner, D. J., et al. 2016, ApJ, 832, 110, doi: 10.3847/0004-637X/832/2/110
- Cridland et al. (2016) Cridland, A. J., Pudritz, R. E., & Alessi, M. 2016, MNRAS, 461, 3274, doi: 10.1093/mnras/stw1511
- Cridland et al. (2017) Cridland, A. J., Pudritz, R. E., Birnstiel, T., Cleeves, L. I., & Bergin, E. A. 2017, MNRAS, 469, 3910, doi: 10.1093/mnras/stx1069
- Du et al. (2015) Du, F., Bergin, E. A., & Hogerheijde, M. R. 2015, ApJ, 807, L32, doi: 10.1088/2041-8205/807/2/L32
- Du et al. (2017) Du, F., Bergin, E. A., Hogerheijde, M., et al. 2017, ApJ, 842, 98, doi: 10.3847/1538-4357/aa70ee
- Eistrup et al. (2016) Eistrup, C., Walsh, C., & van Dishoeck, E. F. 2016, A&A, 595, A83, doi: 10.1051/0004-6361/201628509
- Favre et al. (2013) Favre, C., Cleeves, L. I., Bergin, E. A., Qi, C., & Blake, G. A. 2013, ApJ, 776, L38, doi: 10.1088/2041-8205/776/2/L38
- Fayolle et al. (2016) Fayolle, E. C., Balfe, J., Loomis, R., et al. 2016, ApJ, 816, L28, doi: 10.3847/2041-8205/816/2/L28
- Fedele et al. (2018) Fedele, D., Tazzari, M., Booth, R., et al. 2018, A&A, 610, A24, doi: 10.1051/0004-6361/201731978
- Fogel et al. (2011) Fogel, J. K. J., Bethell, T. J., Bergin, E. A., Calvet, N., & Semenov, D. 2011, ApJ, 726, 29, doi: 10.1088/0004-637X/726/1/29
- Furuya & Aikawa (2014) Furuya, K., & Aikawa, Y. 2014, ApJ, 790, 97, doi: 10.1088/0004-637X/790/2/97
- Furuya et al. (2016) Furuya, K., van Dishoeck, E. F., & Aikawa, Y. 2016, A&A, 586, A127, doi: 10.1051/0004-6361/201527579
- Gaia Collaboration et al. (2016) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2016, A&A, 595, A2, doi: 10.1051/0004-6361/201629512
- Garrod (2013) Garrod, R. T. 2013, ApJ, 778, 158, doi: 10.1088/0004-637X/778/2/158
- Green & Thaddeus (1974) Green, S., & Thaddeus, P. 1974, ApJ, 191, 653, doi: 10.1086/153006
- Günther et al. (2010) Günther, H. M., Matt, S. P., Schmitt, J. H. M. M., et al. 2010, A&A, 519, A97, doi: 10.1051/0004-6361/201014386
- Guzmán et al. (2017) Guzmán, V. V., Öberg, K. I., Huang, J., Loomis, R., & Qi, C. 2017, ApJ, 836, 30, doi: 10.3847/1538-4357/836/1/30
- Hasegawa & Herbst (1993) Hasegawa, T. I., & Herbst, E. 1993, MNRAS, 263, 589, doi: 10.1093/mnras/263.3.589
- Hasegawa et al. (1992) Hasegawa, T. I., Herbst, E., & Leung, C. M. 1992, ApJS, 82, 167, doi: 10.1086/191713
- Haworth et al. (2017) Haworth, T. J., Facchini, S., Clarke, C. J., & Cleeves, L. I. 2017, MNRAS, 468, L108, doi: 10.1093/mnrasl/slx037
- He et al. (2016) He, J., Acharyya, K., & Vidali, G. 2016, ApJ, 823, 56, doi: 10.3847/0004-637X/823/1/56
- Hogerheijde et al. (2011) Hogerheijde, M. R., Bergin, E. A., Brinch, C., et al. 2011, Science, 334, 338, doi: 10.1126/science.1208931
- Huang et al. (2017) Huang, J., Öberg, K. I., Qi, C., et al. 2017, ApJ, 835, 231, doi: 10.3847/1538-4357/835/2/231
- Huang et al. (2018) Huang, J., Andrews, S. M., Cleeves, L. I., et al. 2018, ApJ, 852, 122, doi: 10.3847/1538-4357/aaa1e7
- Isella et al. (2016) Isella, A., Guidi, G., Testi, L., et al. 2016, Physical Review Letters, 117, 251101, doi: 10.1103/PhysRevLett.117.251101
- Kama et al. (2016) Kama, M., Bruderer, S., van Dishoeck, E. F., et al. 2016, A&A, 592, A83, doi: 10.1051/0004-6361/201526991
- Kamp et al. (2013) Kamp, I., Thi, W.-F., Meeus, G., et al. 2013, A&A, 559, A24, doi: 10.1051/0004-6361/201220621
- Kastner et al. (2015) Kastner, J. H., Qi, C., Gorti, U., et al. 2015, ApJ, 806, 75, doi: 10.1088/0004-637X/806/1/75
- Kreidberg et al. (2014) Kreidberg, L., Bean, J. L., Désert, J.-M., et al. 2014, ApJ, 793, L27, doi: 10.1088/2041-8205/793/2/L27
- Krijt et al. (2016) Krijt, S., Ciesla, F. J., & Bergin, E. A. 2016, ApJ, 833, 285, doi: 10.3847/1538-4357/833/2/285
- Lavie et al. (2017) Lavie, B., Mendonça, J. M., Mordasini, C., et al. 2017, AJ, 154, 91, doi: 10.3847/1538-3881/aa7ed8
- Leger et al. (1985) Leger, A., Jura, M., & Omont, A. 1985, A&A, 144, 147
- Li et al. (2013) Li, X., Heays, A. N., Visser, R., et al. 2013, A&A, 555, A14, doi: 10.1051/0004-6361/201220625
- Long et al. (2017) Long, F., Herczeg, G. J., Pascucci, I., et al. 2017, ApJ, 844, 99, doi: 10.3847/1538-4357/aa78fc
- Loomis et al. (2017) Loomis, R. A., Öberg, K. I., Andrews, S. M., & MacGregor, M. A. 2017, ApJ, 840, 23, doi: 10.3847/1538-4357/aa6c63
- Loomis et al. (2018) Loomis, R. A., Öberg, K. I., Andrews, S. M., et al. 2018, AJ, 155, 182, doi: 10.3847/1538-3881/aab604
- Lynden-Bell & Pringle (1974) Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603
- Macintosh et al. (2015) Macintosh, B., Graham, J. R., Barman, T., et al. 2015, Science, 350, 64, doi: 10.1126/science.aac5891
- Madhusudhan et al. (2011) Madhusudhan, N., Harrington, J., Stevenson, K. B., et al. 2011, Nature, 469, 64, doi: 10.1038/nature09602
- Mathis et al. (1977) Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425, doi: 10.1086/155591
- Mawet et al. (2012) Mawet, D., Absil, O., Montagnier, G., et al. 2012, A&A, 544, A131, doi: 10.1051/0004-6361/201219662
- McClure et al. (2016) McClure, M. K., Bergin, E. A., Cleeves, L. I., et al. 2016, ApJsubmitted
- McMullin et al. (2007) McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell, 127
- Miotello et al. (2016) Miotello, A., van Dishoeck, E. F., Kama, M., & Bruderer, S. 2016, A&A, 594, A85, doi: 10.1051/0004-6361/201628159
- Öberg & Bergin (2016) Öberg, K. I., & Bergin, E. A. 2016, ApJ, 831, L19, doi: 10.3847/2041-8205/831/2/L19
- Öberg et al. (2011a) Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011a, ApJ, 743, L16, doi: 10.1088/2041-8205/743/1/L16
- Öberg et al. (2005) Öberg, K. I., van Broekhuizen, F., Fraser, H. J., et al. 2005, ApJ, 621, L33, doi: 10.1086/428901
- Öberg et al. (2011b) Öberg, K. I., Qi, C., Fogel, J. K. J., et al. 2011b, ApJ, 734, 98, doi: 10.1088/0004-637X/734/2/98
- Pinte et al. (2008) Pinte, C., Padgett, D. L., Ménard, F., et al. 2008, A&A, 489, 633, doi: 10.1051/0004-6361:200810121
- Pinte et al. (2018) Pinte, C., Ménard, F., Duchêne, G., et al. 2018, A&A, 609, A47, doi: 10.1051/0004-6361/201731377
- Piso et al. (2016) Piso, A.-M. A., Pegues, J., & Öberg, K. I. 2016, ApJ, 833, 203, doi: 10.3847/1538-4357/833/2/203
- Pollack et al. (1996) Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62, doi: 10.1006/icar.1996.0190
- Prantzos et al. (1996) Prantzos, N., Aubert, O., & Audouze, J. 1996, A&A, 309, 760
- Ripple et al. (2013) Ripple, F., Heyer, M. H., Gutermuth, R., Snell, R. L., & Brunt, C. M. 2013, MNRAS, 431, 1296, doi: 10.1093/mnras/stt247
- Salyk et al. (2011) Salyk, C., Pontoppidan, K. M., Blake, G. A., Najita, J. R., & Carr, J. S. 2011, ApJ, 731, 130, doi: 10.1088/0004-637X/731/2/130
- Schöier et al. (2005) Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005, A&A, 432, 369, doi: 10.1051/0004-6361:20041729
- Schwarz et al. (2016) Schwarz, K. R., Bergin, E. A., Cleeves, L. I., et al. 2016, ApJ, 823, 91, doi: 10.3847/0004-637X/823/2/91
- Schwarz et al. (2018) —. 2018, ApJ, 856, 85, doi: 10.3847/1538-4357/aaae08
- Semenov & Wiebe (2011) Semenov, D., & Wiebe, D. 2011, ApJS, 196, 25, doi: 10.1088/0067-0049/196/2/25
- Spielfiedel et al. (2012) Spielfiedel, A., Feautrier, N., Najar, F., et al. 2012, MNRAS, 421, 1891, doi: 10.1111/j.1365-2966.2011.20225.x
- Testi et al. (2014) Testi, L., Birnstiel, T., Ricci, L., et al. 2014, Protostars and Planets VI, 339, doi: 10.2458/azu_uapress_9780816531240-ch015
- van Dishoeck et al. (2014) van Dishoeck, E. F., Bergin, E. A., Lis, D. C., & Lunine, J. I. 2014, Protostars and Planets VI
- Vasyunin & Herbst (2013) Vasyunin, A. I., & Herbst, E. 2013, ApJ, 762, 86, doi: 10.1088/0004-637X/762/2/86
- Weidenschilling (1980) Weidenschilling, S. J. 1980, Icarus, 44, 172, doi: 10.1016/0019-1035(80)90064-0
- Zhang et al. (2017) Zhang, K., Bergin, E. A., Blake, G. A., Cleeves, L. I., & Schwarz, K. R. 2017, Nature Astronomy, 1, 0130, doi: 10.1038/s41550-017-0130
Appendix A Updated grain surface chemistry
The amount of dust surface area impacts the rate of grain surface chemistry that can occur. Dust surface area can be altered by removing dust grains (but keeping their size spectrum constant), or changing the size spectrum. The default used in the chemical code is 0.1 micron-sized grains with an abundance of per H atom (Fogel et al., 2011). The reason for this default size, which is bigger than the minimum grain size used to compute the temperature structure, is because grains smaller than this can be heated by single-photon events important below m (Leger et al., 1985), and thus have stochastic chemical behavior not adequately described by the present treatment in the code.
In the sweeping approximation for surface chemistry (e.g. Hasegawa et al., 1992), the rates of surface reactions are given by:
where represents the total number of species on the average grain and is the number density of dust grains. The dimensions of are per volume per time. The term factors in for reactions with a barrier, where the probability for tunneling can be approximated by:
where is the area of a site taken to be 1 Å, and is the activation energy of a given reaction. The rate coefficient for the grain-surface reaction of species and is then:
For grains where there is a substantial ice mantle, only the surface monolayer(s) should participate in the chemistry, such that the and terms should be reduced by the ratio of the volume density of ice in the reactive surface area, , to the total volume density of ice, , such that the dilution is given by:
where is the number of surface sites per grain. The numerator, , is the number of sites per volume, which equals .
We define to represent the ratio of adjusted surface area per volume due to grain removal or growth compared to the standard case of 0.1 m sized grains:
If we assume an MRN grain size distribution (Mathis et al., 1977) where the number of grains is proportional to the size of the grains to , we can calculate the population integrated quantity to be:
where is 0.06 m, i.e., the single-photon heating limit of Leger et al. (1985) and varies. is the volumetric density of dust and is the density of silicates. Specifically, we assume two distinct dust populations, “big” grains that have an upper size limit of mm, and “small” grains that have an upper size limit of m as described in Cleeves et al. (2016). Taking into account both populations, and substituting this expression into gives:
From this, we can now express the grain surface area per unit volume by
allowing us to express the dilution term, , as:
In the case where both and are spread over the full ice mantle (i.e., are diluted), the rate decreases by :
The final piece we need is to estimate in terms of . We can do this approximately by rearranging our expression for :
and then calculate the average grain cross section again by computing the weighted average with the MRN grain size distribution:
Using this expression, we estimate the grain number density to be:
where the prefactor comes from , . We can now solve for the rate coefficient :
In this case, for decreasing surface area per unit volume , the reaction rate decreases since the amount of reactive ice (i.e., that on the surface and not buried in the mantle) goes down when there is less surface area for freeze out and individual mantles become thicker.
For the case of one diluted species and one not (like reactions with atomic hydrogen, where we assume that any atomic hydrogen is part of the reactive surface) we only multiply the rate by one dilution factor,
The dependence on drops out because the diluted species is penalized by more inert ice locked up in the non-reactive mantle, while the other more mobile species gains from because there are fewer possible grains to occupy, so any single grain harbors more of the undiluted hydrogen atoms. For the case where neither species is diluted, like H(gr) + H(gr) H,
the rate now increases with decreasing grain surface area per unit volume for the same reason as the case above, i.e., the probability of two H-atoms existing on the same grain and reacting to form H goes up when there are fewer possible grains for them to occupy. The main result of incorporating these reactions is an overall decrease in the efficiency grain-surface chemistry for molecules other than molecular hydrogen, where before under certain conditions CO would be quickly converted into CO, HCO, or CHOH ice. These species still form, but at a reduced abundance, a few percent up to tens of percent of the CO ice abundance just at or below the CO snow surface.