Constraining Gas-phase Carbon, Oxygen, and Nitrogen in the IM Lup Protoplanetary Disk

L. Ilsedore Cleeves; Karin I. Öberg; David J. Wilner; Jane Huang; Ryan A. Loomis; Sean M. Andrews; V. V. Guzman

doi:10.3847/1538-4357/aade96

1. Introduction

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; Cridland et al. 2016, 2017; Piso et al. 2016).

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; Krijt et al. 2016; Öberg & Bergin 2016; Piso et al. 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 H₂O 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 nonthermal) models predict (e.g., in TW Hya's disk; Hogerheijde et al. 2011; Favre et al. 2013; Cleeves et al. 2015; Kama et al. 2016; Schwarz et al. 2016).

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 C₂H and c−C₃H₂ in the TW Hya (see also Kastner et al. 2015) and DM Tau disks. Using chemical models, Bergin et al. (2016) found that 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, > 1, 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 N₂H⁺ and HCN, are more 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 Atacama Large Millimeter/Submillimeter Array (ALMA) observations of C₂H and HCN and H¹³CN. The solar mass star IM Lup harbors a massive gas-rich disk, M_disk ∼ 0.1–0.2 M_⊙ based upon continuum (spectral energy distribution 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 0.5–1 Myr (Mawet et al. 2012).

In Cleeves et al. (2016), we found that IM Lup's CO is under-abundant by a factor of ∼20 compared to an interstellar CO abundance of 1.4 × 10⁻⁴ 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 ∼3–5× 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 C₂H, HCN, and HC¹³N line intensities.

2. Observations

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 J = 3–2 (see also Öberg et al. 2011b), the H¹³CN J = 3–2 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 J = 3–2 and H¹³CN J = 3–2 are presented in Huang et al. (2017) (see also Guzmán et al. 2017) and are not reproduced here.

The C₂H N = 3–2 hyperfine complex was observed with ALMA as part of the Cycle 3 2015.1.00964.S program (PI: Öberg) on 2016 May 1, 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 s and average both polarizations. The continuum is estimated from line-free channels within the same spectral window, which is then subtracted in the uv-plane to obtain continuum-subtracted images. We detect two blended hyperfine pairs of C₂H N =3–2, the $J=\tfrac{5}{2}-\tfrac{3}{2}$ F = 2–1 and F = 3–2 pair and the $J=\tfrac{7}{2}-\tfrac{5}{2}$ F = 3–2 and F = 4–3 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 the 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 C₂H transitions is shown in Figure 2 for each of the hyperfine pairs integrated within a 4'' radius circular aperture. Table 1 provides the C₂H N = 3–2 disk-integrated fluxes within this aperture, with errors estimated from a 4'' circular aperture in line-free portions of the spectrum combined with a 10% calibration uncertainty added in quadrature.

**Figure 1.** Channel maps of C₂H N = 3–2 for the (a) $J=\tfrac{5}{2}-\tfrac{3}{2}$ pair and (b) $J=\tfrac{7}{2}-\tfrac{5}{2}$ pair, channel averaged by a factor of two for clarity. The solid contour line indicates 3σ. The beam is in the lower left corner. The velocity relative to the local standard of rest in km s⁻¹ is indicated in the bottom right corner of each panel. The purple ellipse on the right panels indicates the scale of the millimeter thermal dust emission, R = 313 au (Cleeves et al. 2016).
Download figure:
Standard image High-resolution image

**Figure 1.** Channel maps of C₂H N = 3–2 for the (a) $J=\tfrac{5}{2}-\tfrac{3}{2}$ pair and (b) $J=\tfrac{7}{2}-\tfrac{5}{2}$ pair, channel averaged by a factor of two for clarity. The solid contour line indicates 3σ. The beam is in the lower left corner. The velocity relative to the local standard of rest in km s⁻¹ is indicated in the bottom right corner of each panel. The purple ellipse on the right panels indicates the scale of the millimeter thermal dust emission, R = 313 au (Cleeves et al. 2016).
Download figure:
Standard image High-resolution image

**Figure 2.** C₂H N = 3–2 spectra extracted from a circular 4'' mask. The dotted line indicates the line centers for each pair in the rest frame of the F = 2–1 (top) and F = 3–2 (bottom) transitions.
Download figure:
Standard image High-resolution image

Table 1. C₂H N = 3–2 Line Observations

Transition	Rest	Disk-integrated
	Freq.	Flux Density (Blended)
	[GHz]	[Jy km s⁻¹]
$J=\tfrac{5}{2}-\tfrac{3}{2}$ F = 2–1	262.067	0.56 ± 0.08
$J=\tfrac{5}{2}-\tfrac{3}{2}$ F = 3–2	262.065
$J=\tfrac{7}{2}-\tfrac{5}{2}$ F = 3–2	262.006	0.76 ± 0.11
$J=\tfrac{7}{2}-\tfrac{5}{2}$ F = 4–3	262.004

Note. Uncertainties are quadrature combined rms scatter and 10% calibration uncertainty.

Download table as: ASCII Typeset image

3. Methods

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 C₂H, HCN, and H¹³CN 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_dust = 0.0017 M_⊙. We furthermore assume small grains follow the gas distribution and allow the millimeter grains to have a separate power-law distribution with a cutoff at the millimeter emission edge at 313 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 (≳1 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 T_eff = 3900 K and R_* = 2.5 R_⊙ (Pinte et al. 2008; Cleeves et al. 2016). We note Alcalá et al. (2017) provide updated stellar parameters for IM Lup; however, 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 G₀ = 4 Habing (Cleeves et al. 2016; Haworth et al. 2017; Pinte et al. 2018).

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 of 4.3 × 10³⁰ erg s⁻¹ provided in Günther et al. (2010). 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 C₂H and HCN using the 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 the Appendix. 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 H¹³CN J = 3–2 line radiation transfer, we simply take the model HCN abundance and assume an isotope ratio of ¹²C/¹³C = 70 (Prantzos et al. 1996).

The two main variables considered in the modeling are:

1.
The initial C, N, and O abundances (Section 3.2.2).
2.
The cosmic ray ionization rate, ζ_CR.

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, ζ_CR ≤ 2 × 10⁻¹⁹ 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; ζ_CR = 1.0 × 10⁻²⁰ s⁻¹ per H₂), "solar system maximum" (ssx; ζ_CR = 2.0 × 10⁻¹⁹ s⁻¹ per H₂), and the "solar system minimum" (ssm; ζ_CR = 1.3 × 10⁻¹⁸ 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 the grain-surface chemistry, since our code does not treat the multilayered 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 overpredict the atomic N density and underpredict the N₂ density.

3.2.2. Model Abundances

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 C₂H 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 (∼10³ 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 (z/r ≳ 0.2) 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, χ(H₂O_ice) = 8 × 10⁻⁵ per H (Boogert et al. 2015), and thus the 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).

Table 2. Fixed Chemical Abundances Besides CO/C⁺/H₂O Ice Relative to Total H Atoms

Molecule	Abundance	Molecule	Abundance
H₂	5.00 × 10⁻¹	He	1.40 × 10⁻¹
N₂	3.75 × 10⁻⁵	CS	4.00 × 10⁻⁹
SO	5.00 × 10⁻⁹	HCO⁺	9.00 × 10⁻⁹
${{\rm{H}}}_{3}^{+}$	1.00 × 10⁻⁸	C₂H	8.00 × 10⁻⁹
Si⁺	1.00 × 10⁻⁹	Mg⁺	1.00 × 10⁻⁹
Fe⁺	1.00 × 10⁻⁹

Note. N₂ abundance varied in Section 4.3.

Download table as: ASCII Typeset image

Table 3. Abundances of Key C and O Species

C/O	C⁺	CO	H₂O(gr)
0.08	0.0	7 × 10⁻⁶	8 × 10⁻⁵
0.47	0.0	7 × 10⁻⁶	8 × 10⁻⁶
0.64	0.0	7 × 10⁻⁶	4 × 10⁻⁶
0.81	0.0	7 × 10⁻⁶	1.6 × 10⁻⁶
0.9	0.0	7 × 10⁻⁶	8 × 10⁻⁷
0.95	0.0	7 × 10⁻⁶	4 × 10⁻⁷
1.0	3.5 × 10⁻⁶	3.5 × 10⁻⁶	3.5 × 10⁻⁶
1.86	9.5 × 10⁻⁶	3.5 × 10⁻⁶	3.5 × 10⁻⁶
3.71	2.25 × 10⁻⁵	3.5 × 10⁻⁶	3.5 × 10⁻⁶

Note. Models with C/O ≤ 1 are carbon limited, while models with C/O ≥ 1 are oxygen limited. In the latter case, the oxygen is split between CO and H₂O(gr).

Download table as: ASCII Typeset image

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 1.3 × 10⁻⁴ per H (Ripple et al. 2013), updated from previously 1.4 × 10⁻⁴ 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 7 × 10⁻⁶, but not excess. For example, the carbon abundance can be 7 × 10⁻⁶ per H, and the oxygen abundance can be this value or greater (carbon-limited models). Similarly, oxygen can have the same 7 × 10⁻⁶ 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.

Even though the chemical reprocessing timescales are short in the upper layers (z/r ≥ 0.25) 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 H₂O, 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 C₂H N = 3–2, HCN J = 3–2, and H¹³CN J = 3–2 observations using the non-LTE radiative transfer code LIME v1.8⁷ (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 H¹³CN collision rate data is assumed to be the same and sourced from Green & Thaddeus (1974) and does not include hyperfine structure. The C₂H collision rate data is from Spielfiedel et al. (2012). The non-LTE analysis is especially important given that the C₂H is abundant in our models across a wide range of densities, down to n_H ∼ 10⁵ cm⁻³. We simulate each of the C₂H $J=\tfrac{5}{2}-\tfrac{3}{2}$ and $J=\tfrac{7}{2}-\tfrac{5}{2}$ 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 ∼1000 au, the HCN and C₂H and HCN extend to about ∼500 au, and so we limit our emission models to this radius to only constrain the abundances where C₂H 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 R ≲ 40 au, where the C₂H 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 C₂H, HCN, and H¹³CN simulations is 0.28, 0.28, and 0.35 km s⁻¹, but we simulate the cubes at 40× 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 48° (Cleeves et al. 2016). The updated DR2 distance is 158 ± 3 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)⁸ 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 the 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 ζ_CR for 27 models total.

4. Results

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 C₂H is strongly sensitive to the C/O ratio in the gas, confirming the early results of Bergin et al. (2016). The C₂H abundance is most sensitive for our models below a C/O of 1.86, where there is a sharp column density transition straddling C/O of ∼1, clearly seen in the two orders of magnitude jump in C₂H column densities going from a C/O of 1.86–0.95 in Figure 4. C₂H is essentially unaffected by the cosmic ray ionization rate for the three model values considered. This lack of dependence occurs because C₂H 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 ≥1, the C₂H column density forms a narrow ring or is centrally peaked. At low C/O ratios, the C₂H column density forms a wide ring that peaks outside of the millimeter dust disk (R_mm = 313 au). This transition occurs once the layer of C₂H interior to R ≲ 300 au becomes thin, with little column density compared to the outer disk. Correspondingly, the radial morphology of C₂H 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 ≥ 1. Higher values of C/O and ζ_CR 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 C₂H, though we find the HCN abundance is generally less than that of C₂H for a given C/O value. One morphological difference is that HCN extends vertically deeper, with a base of z ∼ 60 au at R = 300 au compared to C₂H, which disappears below z ∼ 80 au for C/O < 1.

4.2. Fits to C₂H Observations

Based upon the chemical modeling results in Section 4.1, we can use the C₂H 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 χ²+1) between the data visibilities and our model visibilities varying initial gas composition (C/O) and ζ_CR, where smaller values indicate better fits. Note, adding 1 allows us to plot the difference on log scale. For the C₂H N = 3–2 $J=\tfrac{5}{2}-\tfrac{3}{2}$ pair of lines a C/O ratio of ∼0.81 is favored, while for the C₂H N = 3–2 $J=\tfrac{7}{2}-\tfrac{5}{2}$ pair, a C/O of ∼0.7 is the best match. Combining both line pairs favors the C/O = 0.81 model. We also find that the C₂H 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 50× in the layer C₂H is present compared to an interstellar water ice abundance of 8 × 10⁻⁵ per H (Boogert et al. 2015).

**Figure 5.** χ² comparison between the observed visibilities and models for each of the lines or pairs (see top of each panel). Blended hyperfine components for the C₂H lines are modeled simultaneously. From left to right, ${\chi }_{\min }^{2}$ is 4.58 × 10⁶, 4.39 × 10⁶, 6.43 × 10⁶, and 6.70 × 10⁶, so in all cases the reduced χ² is close to 1, with the minimum being the closet model-data fit. For H¹³CN the models are consistent with nondetection in the image plane for C/O < 1. Note: the C/O variations reflect the total elemental ratio in the simulations (gas and ice), see Section 5.2 for the gas-phase C/O values. These values should be interpreted as surface (z/r ≳ 0.2) constraints, i.e., the region of abundant C₂H (Section 4.1).
Download figure:
Standard image High-resolution image

In the same figure, we also show the HCN J = 3–2 and H¹³CN J = 3–2 results for comparison. C/O values of 0.7–0.9 also fit the HCN data reasonably well, while a C/O of ≳1 predicts detectable H¹³CN J = 3–2 image-plane emission, regardless of the cosmic ray ionization rate, inconsistent with observations.

4.3. Fits to HCN and H¹³CN Observations

We now consider models with bulk nitrogen depletion to see whether we can arrive at a better fit for HCN J = 3–2 at a fixed value for C/O. Taking bulk C/O = 0.81 from Section 4.2, we reduce the nitrogen abundance from the fiducial value of 7.5 × 10⁻⁵ N per H (i.e., 3.75 × 10⁻⁵ N₂ per H, see Table 2). Figure 6 shows the (χ²–χ_min²+1) values for this sub-grid of reduced nitrogen models.

**Figure 6.** χ² comparison between the ALMA visibilities and models for HCN J = 3–2 for differing amounts of initial nitrogen depletion. The colors are the same as those in Figure 5 and indicate the CR ionization rate. ${\chi }_{\min }^{2}$ is 6.43 × 10⁶, so again all models have a Δχ² close to 1, with smaller values being better fits to the data. For intermediate to low CR rates, no nitrogen depletion is favored. For the higher rate, the models are insensitive to nitrogen depletion factors of ≤4.
Download figure:
Standard image High-resolution image

The global best fit is the low ζ_CR value, 1 × 10⁻²⁰ 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 ζ_CR = 1.3 × 10⁻¹⁸ s⁻¹ model, favors anywhere between no and a 4× reduction in bulk nitrogen. However, even this factor is small compared to the ∼20× depletion of CO, or that implied for water and oxygen not in CO based on the C₂H results, i.e., a ∼50× depletion.

5. Discussion

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 CH₃OH, 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, N₂H⁺ 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 comparison 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 C₂H emits, z/r ≳ 0.2), while H₂O is primarily ice with some gas-phase H₂O 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 (ζ_CR ∼ 10⁻¹⁹ 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 H₂O ice is quickly dissociated to gas-phase atomic oxygen, such that the gas C/O ratio is equal to the bulk ratio of ∼0.8. 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 ∼1. The layered C/O structure is visible here, where C/O in the gas is 0.8 above normalized heights of z/r > 0.3 and ∼1 below this layer until CO starts to freeze out, at z/r ≲ 0.1.

The N/O ratio on the other hand is ≫1 throughout the disk atmosphere, and ∼10 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 >100.

At high velocities we find a discrepancy between the data and modeled C₂H, which can be seen in Figure 8. The high velocity emission of C₂H is weaker in the models than in the data. Given the signal to noise ratio of the C₂H 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 C₂H (Bergin et al. 2016), there is support for local variations in bulk C/O, N/O, etc. In this instance, the bright high velocity C₂H 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.

**Figure 8.** Sample model channel maps at the line center and at the line wings for HCN J = 3–2 (left) and C₂H N = 3–2 $J=\tfrac{7}{2}-\tfrac{5}{2}$ (right) lines, along with the data (top, respectively) for comparison. The rows show varying C/O ratios while the columns show the low end and high end of the cosmic ray ionization rate considered.
Download figure:
Standard image High-resolution image

We also examined models where we allowed the disk to extend out to 1000 au, beyond the ∼500 au radial region where C₂H and HCN are observed. Since the model C₂H abundance continues to rise in the outer disk for the C/O = 0.81 case (see Figure 3), we predicted 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 C₂H toward this source will help disentangle these scenarios.

6. Conclusions

Using detailed physical and chemical models constrained by CO and dust observations, along with C₂H 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 C₂H and HCN emit, primarily at normalized heights above z/r ≳ 0.2. In these layers, the C₂H observations favor a super-solar elemental C/O ratio of ∼0.8. 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, ∼10. 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., Krijt et al. 2016; Öberg & Bergin 2016; Piso 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; Fedele et al. 2018; Huang 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, 5 × 10⁻⁴, 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 4–20× 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, C₂H, 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).

This paper makes use of the following ALMA data: ADS/JAO.ALMA#2015.1.00964.S and ADS/JAO.ALMA#2013.1.00226.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan) and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. L.I.C. acknowledges the support of NASA through Hubble Fellowship grant HST-HF2-51356.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555. J.H. acknowledges support from the National Science Foundation Graduate Research Fellowship under Grant No. DGE-1144152. All simulations were carried out using the Smithsonian Institution High Performance Cluster (SI/HPC).

Appendix: 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 χ_gr⁰ = 6 × 10⁻¹² 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 0.06 μ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:

$\begin{eqnarray}&&{R}_{i,j}={b}_{i,j}({R}_{\mathrm{diff},{\rm{i}}}+{R}_{\mathrm{diff},{\rm{j}}}){N}_{i}{N}_{j}{n}_{\mathrm{gr}},\end{eqnarray} \tag{ 1 }$

where N_x represents the total number of species x on the average grain and n_gr is the number density of dust grains. The dimensions of R_i,j are per volume per time. The b_i,j term factors in for reactions with a barrier, where the probability for tunneling can be approximated by:

$\begin{eqnarray}&&{b}_{i,j}=\exp [-2(a/\hslash )\sqrt{2\mu {E}_{a}}],\end{eqnarray} \tag{ 2 }$

where a is the area of a site taken to be 1 Å, and E_a is the activation energy of a given reaction. The rate coefficient κ_i,j for the grain-surface reaction of species i and j is then:

$\begin{eqnarray}&&{\kappa }_{i,j}=\displaystyle \frac{{R}_{i,j}}{{n}_{i}{n}_{j}}\end{eqnarray} \tag{ 3 }$

$\begin{eqnarray}&&=\,\displaystyle \frac{{b}_{i,j}({R}_{\mathrm{diff},{\rm{i}}}+{R}_{\mathrm{diff},{\rm{j}}}){N}_{i}{N}_{j}{n}_{\mathrm{gr}}}{{n}_{i}{n}_{j}}\end{eqnarray} \tag{ 4 }$

$\begin{eqnarray}&&=\,\displaystyle \frac{{b}_{i,j}({R}_{\mathrm{diff},{\rm{i}}}+{R}_{\mathrm{diff},{\rm{j}}})\tfrac{{n}_{i}}{{n}_{\mathrm{gr}}}\tfrac{{n}_{j}}{{n}_{\mathrm{gr}}}{n}_{\mathrm{gr}}}{{n}_{i}{n}_{j}}\end{eqnarray} \tag{ 5 }$

$\begin{eqnarray}&&=\,{b}_{i,j}({R}_{\mathrm{diff},{\rm{i}}}+{R}_{\mathrm{diff},{\rm{j}}})\displaystyle \frac{1}{{n}_{\mathrm{gr}}}.\end{eqnarray} \tag{ 6 }$

For grains where there is a substantial ice mantle, only the surface monolayer(s) should participate in the chemistry, such that the N_i and N_j terms should be reduced by the ratio of the volume density of ice in the reactive surface area, N_sites n_ice, to the total volume density of ice, n_ice, such that the dilution d is given by:

$\begin{eqnarray}&&d=\displaystyle \frac{{n}_{\mathrm{gr}}{N}_{\mathrm{sites}}}{{n}_{\mathrm{ice}}},\end{eqnarray} \tag{ 7 }$

where N_sites is the number of surface sites per grain. The numerator, n_grN_sites, is the number of sites per volume, which equals n_grσ_gr/σ_site.

We define Π_adj 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:

$\begin{eqnarray}&&{{\rm{\Pi }}}_{\mathrm{adj}}=\displaystyle \frac{{\sigma }_{\mathrm{gr}}{n}_{\mathrm{gr}}}{{\sigma }_{\mathrm{gr}}^{0}{n}_{\mathrm{gr}}^{0}}.\end{eqnarray} \tag{ 8 }$

If we assume an MRN grain size distribution (Mathis et al. 1977) where the number of grains n_gr is proportional to the size of the grains to r^−3.5, we can calculate the population integrated quantity n_grσ_gr to be:

$\begin{eqnarray}&&{n}_{\mathrm{gr}}{\sigma }_{\mathrm{gr}}=-3{\rho }_{d}/{\rho }_{{si}}\left[\displaystyle \frac{{r}_{\max }^{-0.5}-{r}_{\min }^{-0.5}}{{r}_{\max }^{0.5}-{r}_{\min }^{0.5}}\right],\end{eqnarray} \tag{ 9 }$

where r_min is 0.06 μm, i.e., the single-photon heating limit of Leger et al. (1985) and r_max varies. ρ_d is the volumetric density of dust and ρ_si is the density of silicates. Specifically, we assume two distinct dust populations, "big" grains that have an upper size limit of r_max = 1 mm, and "small" grains that have an upper size limit of r_max = 1 μm as described in Cleeves et al. (2016). Taking into account both populations, and substituting this expression into Π_adj gives:

$\begin{eqnarray}&&{{\rm{\Pi }}}_{\mathrm{adj}}=100\displaystyle \frac{{\rho }_{d}}{{\rho }_{\mathrm{gas}}}\displaystyle \frac{{f}_{\mathrm{big}}\left[\tfrac{{r}_{\max ,\mathrm{big}}^{-0.5}-{r}_{\min }^{-0.5}}{{r}_{\max ,\mathrm{big}}^{0.5}-{r}_{\min }^{0.5}}\right]+{f}_{\mathrm{sm}}\left[\tfrac{{r}_{\max ,\mathrm{sm}}^{-0.5}-{r}_{\min }^{-0.5}}{{r}_{\max ,\mathrm{sm}}^{0.5}-{r}_{\min }^{0.5})}\right]}{\left[\tfrac{{r}_{\max ,\mathrm{sm}}^{-0.5}-{r}_{\min }^{-0.5}}{{r}_{\max ,\mathrm{sm}}^{0.5}-{r}_{\min }^{0.5}}\right]}.\end{eqnarray} \tag{ 10 }$

From this, we can now express the grain-surface area per unit volume by

$\begin{eqnarray}&&{\sigma }_{\mathrm{gr}}{n}_{\mathrm{gr}}={{\rm{\Pi }}}_{\mathrm{adj}}{\sigma }_{\mathrm{gr}}^{0}{n}_{\mathrm{gr}}^{0}={{\rm{\Pi }}}_{\mathrm{adj}}{\sigma }_{\mathrm{gr}}^{0}{\chi }_{\mathrm{gr}}^{0}{n}_{{\rm{H}}},\end{eqnarray} \tag{ 11 }$

allowing us to express the dilution term, d, as:

$\begin{eqnarray}&&d=\displaystyle \frac{{{\rm{\Pi }}}_{\mathrm{adj}}{\sigma }_{\mathrm{gr}}^{0}{\chi }_{\mathrm{gr}}^{0}{n}_{{\rm{H}}}/{\sigma }_{\mathrm{site}}}{{n}_{\mathrm{ice}}}.\end{eqnarray} \tag{ 12 }$

In the case where both i and j are spread over the full ice mantle (i.e., are diluted), the rate decreases by d²:

$\begin{eqnarray}&&{\kappa }_{i,j}={d}^{2}\displaystyle \frac{{b}_{i,j}({R}_{\mathrm{diff},{\rm{i}}}+{R}_{\mathrm{diff},{\rm{j}}})}{{n}_{\mathrm{gr}}}\end{eqnarray} \tag{ 13 }$

$\begin{eqnarray}&&=\,{\left(\displaystyle \frac{{{\rm{\Pi }}}_{\mathrm{adj}}{\sigma }_{\mathrm{gr}}^{0}{\chi }_{\mathrm{gr}}^{0}{n}_{{\rm{H}}}}{{\sigma }_{\mathrm{site}}{n}_{\mathrm{ice}}}\right)}^{2}\displaystyle \frac{{b}_{i,j}({R}_{\mathrm{diff},{\rm{i}}}+{R}_{\mathrm{diff},{\rm{j}}})}{{n}_{\mathrm{gr}}}.\end{eqnarray} \tag{ 14 }$

The final piece we need is to estimate n_gr in terms of Π_adj. We can do this approximately by rearranging our expression for σ_grn_gr:

$\begin{eqnarray}&&{n}_{\mathrm{gr}}=\displaystyle \frac{{{\rm{\Pi }}}_{\mathrm{adj}}{\sigma }_{\mathrm{gr}}^{0}{n}_{\mathrm{gr}}^{0}}{\langle {\sigma }_{\mathrm{gr}}\rangle },\end{eqnarray} \tag{ 15 }$

and then calculate the average grain cross section again by computing the weighted average σ with the MRN grain size distribution:

$\begin{eqnarray}&&\langle {\sigma }_{\mathrm{gr}}\rangle =5\pi \left[\displaystyle \frac{{r}_{\max }^{-0.5}-{r}_{\min }^{-0.5}}{{r}_{\max }^{-2.5}-{r}_{\min }^{-2.5}}\right]\end{eqnarray} \tag{ 16 }$

$\begin{eqnarray}&&\approx \,5\pi \left[\displaystyle \frac{-{r}_{\min }^{-0.5}}{-{r}_{\min }^{-2.5}}\right]=5\pi {r}_{\min }^{2}.\end{eqnarray} \tag{ 17 }$

Using this expression, we estimate the grain number density to be:

$\begin{eqnarray}&&{n}_{\mathrm{gr}}=\displaystyle \frac{{{\rm{\Pi }}}_{\mathrm{adj}}{\sigma }_{\mathrm{gr}}^{0}{\chi }_{\mathrm{gr}}^{0}{n}_{{\rm{H}}}}{5\pi {r}_{\min }^{2}}\end{eqnarray} \tag{ 18 }$

$\begin{eqnarray}&&\approx \,\displaystyle \frac{1}{1.8}{{\rm{\Pi }}}_{\mathrm{adj}}{\chi }_{\mathrm{gr}}^{0}{n}_{{\rm{H}}},\end{eqnarray} \tag{ 19 }$

where the prefactor comes from ${\sigma }_{\mathrm{gr}}^{0}=\pi {(0.1\mu m)}^{2}$ , r_min =0.06 μm. We can now solve for the rate coefficient κ:

$\begin{eqnarray}&&\begin{array}{l}{\kappa }_{i,j}={\left(\displaystyle \frac{{{\rm{\Pi }}}_{\mathrm{adj}}{\sigma }_{\mathrm{gr}}^{0}{\chi }_{\mathrm{gr}}^{0}{n}_{{\rm{H}}}}{{\sigma }_{\mathrm{site}}{n}_{\mathrm{ice}}}\right)}^{2}{b}_{i,j}({R}_{\mathrm{diff},{\rm{i}}}+{R}_{\mathrm{diff},{\rm{j}}})\displaystyle \frac{5\pi {r}_{\min }^{2}}{{{\rm{\Pi }}}_{\mathrm{adj}}{\sigma }_{\mathrm{gr}}^{0}{\chi }_{\mathrm{gr}}^{0}{n}_{{\rm{H}}}}\end{array}\end{eqnarray} \tag{ 20 }$

$\begin{eqnarray}&&=\,\displaystyle \frac{{N}_{\mathrm{sites}}^{2}}{{n}_{\mathrm{ice}}^{2}}{{\rm{\Pi }}}_{\mathrm{adj}}{\chi }_{\mathrm{gr}}^{0}{n}_{{\rm{H}}}{b}_{i,j}({R}_{\mathrm{diff},{\rm{i}}}+{R}_{\mathrm{diff},{\rm{j}}})\displaystyle \frac{5\pi {r}_{\min }^{2}}{{\sigma }_{\mathrm{gr}}^{0}}.\end{eqnarray} \tag{ 21 }$

In this case, for decreasing surface area per unit volume Π_adj, 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 d factor,

$\begin{eqnarray}&&\begin{array}{l}{\kappa }_{i,j}=\left(\displaystyle \frac{{{\rm{\Pi }}}_{\mathrm{adj}}{\sigma }_{\mathrm{gr}}^{0}{\chi }_{\mathrm{gr}}^{0}{n}_{{\rm{H}}}}{{\sigma }_{\mathrm{site}}{n}_{\mathrm{ice}}}\right){b}_{i,j}({R}_{\mathrm{diff},{\rm{i}}}+{R}_{\mathrm{diff},{\rm{j}}})\displaystyle \frac{5\pi {r}_{\min }^{2}}{{{\rm{\Pi }}}_{\mathrm{adj}}{\sigma }_{\mathrm{gr}}^{0}{\chi }_{\mathrm{gr}}^{0}{n}_{{\rm{H}}}}\end{array}\end{eqnarray} \tag{ 22 }$

$\begin{eqnarray}&&=\,\displaystyle \frac{{N}_{\mathrm{sites}}}{{n}_{\mathrm{ice}}}{b}_{i,j}({R}_{\mathrm{diff},{\rm{i}}}+{R}_{\mathrm{diff},{\rm{j}}})\displaystyle \frac{5\pi {r}_{\min }^{2}}{{\sigma }_{\mathrm{gr}}^{0}}.\end{eqnarray} \tag{ 23 }$

The dependence on Π_adj drops out because the diluted species is penalized by more inert ice locked up in the nonreactive mantle, while the other more mobile species gains from Π_adj 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) $\to$ H₂,

$\begin{eqnarray}&&\begin{array}{l}{\kappa }_{i,j}={b}_{i,j}({R}_{\mathrm{diff},{\rm{i}}}+{R}_{\mathrm{diff},{\rm{j}}})\displaystyle \frac{5\pi {r}_{\min }^{2}}{{{\rm{\Pi }}}_{\mathrm{adj}}{\sigma }_{\mathrm{gr}}^{0}{\chi }_{\mathrm{gr}}^{0}{n}_{{\rm{H}}}},\end{array}\end{eqnarray} \tag{ 24 }$

the rate now increases with decreasing grain-surface area per unit volume Π_adj 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₂, H₂CO, or CH₃OH 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.

Constraining Gas-phase Carbon, Oxygen, and Nitrogen in the IM Lup Protoplanetary Disk

Article metrics

Permissions

Author affiliations

Author notes

ORCID iDs

Dates

Abstract

1. Introduction

2. Observations

2.1. ALMA Observations and Data Reduction