The effect of carbon grain destruction on the chemical structure of protoplanetary disks
The bulk composition of the Earth is dramatically carbon poor compared to that of the interstellar medium, and this phenomenon extends to the asteroid belt. To interpret this carbon deficit problem, the carbonaceous component in grains must have been converted into the gas-phase in the inner regions of protoplanetary disks prior to planetary formation. We examine the effect of carbon grain destruction on the chemical structure of disks by calculating the molecular abundances and distributions using a comprehensive chemical reaction network. When carbon grains are destroyed and the elemental abundance of the gas becomes carbon-rich, the abundances of carbon-bearing molecules, such as HCN and carbon-chain molecules, increase dramatically near the midplane, while oxygen-bearing molecules, such as \ceH2O and \ceCO2, are depleted. We compare the results of these model calculations with the solid carbon fraction in the solar system. Although we find a carbon depletion gradient, there are some quantitative discrepancies: the model shows a higher value at the position of asteroid belt and a lower value at the location of the Earth. In addition, using the obtained molecular abundances distributions, coupled with line radiative transfer calculations, we make predictions for ALMA to potentially observe the effect of carbon grain destruction in nearby protoplanetary disks. The results indicate that HCN, \ceH^13CN and c-\ceC3H2 may be good tracers.
The protoplanetary disks (PPD) around low-mass stars are analogs of the solar nebula. Disks are composed of dust and gas surrounding a central pre-main-sequence star. As dust grains settle down to the midplane, they collide, grow and form planetesimals. These are the ingredients which form rocky planets, the cores of gas-giant planets and small bodies in our solar system. Therefore, the composition of rocky bodies contain information on the refractory materials available during the formation of solar system and provide clues on both the local and large-scale chemical and physical processes in disks.
Figure 1 shows that the abundance of carbon relative to silicon has clear variances across the solar system bodies (e.g., Geiss 1987; Lee et al. 2010; Pontoppidan et al. 2014; Bergin et al. 2015; Anderson et al. 2017; Bardyn et al. 2017). The similar amount of carbon relative to silicon in the Sun and ISM suggests that this value is representative of the ratio at the beginning of the formation of the planets. Although the surface of our Earth is covered with organic compounds, including even live organisms, if one excludes the uncertain amount of carbon in the core, the Earth BSE value (Bulk Silicate Earth: includes only ocean, atmosphere and silicate mantle) is 4 orders of magnitude lower than the solar ratio. Even when taking into account that some light elements, including carbon, are possibly incorporated inside the core of Earth, the abundance ratio remains significantly lower than solar (Allègre et al. 2001). In the asteroid belt, the amount of carbon in meteorites relative to silicon is 1 to 2 orders of magnitude less than the solar value. Moving farther out, to the comet region, the carbon fractions become again comparable to that in the Sun. The C/Si ratio shows a clear correlation with heliocentric distance. Highly depleted elements in the solar system imply that they were in a phase too volatile to condense, and as such, they were not incorporated into planetesimals. Thus any gaseous carbon species which did not freeze onto the surface of grains nor become incorporated into refractory organic material will be either accreted onto the central star or dispersed out of the disk. In contrast, refractory carbon or carbon contained in ice mantles would have been incorporated into solid bodies producing, for example, the solar-like carbon abundance relative to silicon seen in comets.
The described trend of carbon depletion in the solar system is a long-standing problem in the community (Geiss 1987; Lee et al. 2010; Pontoppidan et al. 2014; Bergin et al. 2015; Anderson et al. 2017). This is the so-called carbon deficit problem. This issue exists not only in the solar system but also in white dwarf systems (Jura 2006). Therefore, the cause might be universal throughout the Galaxy. In the ISM about 60% of carbon is in the form of either graphite, amorphous carbon grains, or polycyclic aromatic hydrocarbons (PAHs) (Savage & Sembach 1996; Draine 2003). If the solar system inherited ISM-like material, then refractory carbon must have been destroyed prior to planetary formation to account for the observed carbon depletion. Some possible ways to destroy carbon grains have been discussed: (1) carbon can be released from the solid to the gas phase by hot oxygen atoms, if the gas temperature is higher than 500 K and the grain size is within 0.1-1 (Lee et al. 2010). (2) carbon aggregates such as graphite can react with oxygen-bearing species, e.g., OH (Finocchi et al. 1997; Gail 2001; Gail 2002). (3) PAH destruction by X-ray and extreme ultraviolet (EUV) photons (Siebenmorgen & Krügel 2010).
We start by considering the nature of a PPD, models of which are complex due to the interaction between physical and chemical processes of the dust and gas. The environment of the protoplanetary disk is diverse. The density varies by about 10 orders of magnitude and the temperature ranges from a few K to thousands of K. From the surface to the midplane, the main structure can be classified into three layers; a surface PDR (Photon Dominated Region)-like layer, a warm molecular layer and a cold midplane. The PDR-like layer mainly consists of atoms and ions due to the severe stellar UV radiation and interstellar radiation field the latter of which can influence the disk structure at larger radial distances from the star. In this layer, chemistry induced by photodissociation and photoionization dominates. Moving downward to the warm molecular layer, molecules here can survive due to partial shielding from the radiation field. In this layer, gas phase two-body reactions (neutral-neutral reactions and ion-neutral reactions) are active and set the molecular composition. Towards the midplane, excepting the very inner part, the temperature decreases due to the efficient shielding of radiation. Therefore, in the midplane, molecules can freeze on the surfaces of dust grains. Chemical modeling of PPDs has been widely described in, for example, Bergin et al. (2007), Henning & Semenov (2013) and Dutrey et al. (2014).
The distribution of molecules is important for planetary formation, especially near the midplane. The final composition of planets is related to the location of snowlines in comparison with the planet’s formation location. The balance between desorption and adsorption of gas-phase species onto grains determines the composition of the matter forming planetesimals and planets. From the observational point of view, several molecules, such as CO, \ceCO2, CN, HCN, \ceH2CO, \ceC2H2, OH (e.g., Dutrey et al. 1997; Thi et al. 2004; Pontoppidan et al. 2010) and some deuterated molecules, \ceDCO+, and DCN (e.g., van Dishoeck et al. 2003; Qi et al. 2008) have been detected in the disks so far. Recently, some relatively complex molecules in PPDs have been detected, such as \ceCH3CN, \ceCH3OH and \ceHCOOH (Öberg et al. 2015; Walsh et al. 2016; Favre et al. 2018; Loomis et al. 2018; Bergner et al. 2018). So far, the total molecular inventory in PPDs is 20. Thanks to the high spatial resolution and sensitivity of ALMA, we have opportunities to reveal further detail, even in the inner region of disks ( 10’s of au).
In this paper, we calculate the chemistry in protoplanetary disks taking into account carbon grain destruction in the inner region to study how it affects the molecular abundance distribution in the disk. We compare the results between cases with and without carbon grain destruction, finding the candidate species which have the largest variations in abundance between two cases. To investigate whether or not the effect of carbon grain destruction is observable with ALMA, we perform the radiative transfer calculations of molecular line emission from identified species in T Tauri and Herbig Ae disks.
We describe our disk model and the radiative transfer calculations in Section 2. The results and discussion are presented in Section 3: Section 3.1 presents the effect of carbon grain destruction on the molecular abundances and distributions (HCN, \ceCH4, \ceC2H2, \cec-C3H2 carbon-chain molecules, \ceH2O, OH, \ceO2 and \ceCO2). Section 3.2 discusses the radial distribution of the solid carbon fraction in the disk in our model in comparison with the solar system, and Section 3.3 presents the prediction for ALMA observation based on our calculations. We summarize and state the conclusions of this work in Section 4.
2 Protoplanetary disk model
2.1 Physical model
We use physical models of protoplanetary disks generated using the methodology described in Nomura & Millar (2005) with the addition of X-ray heating as described in Nomura et al. (2007). We model an axisymmetric Keplerian disk and two types of central stars are considered. The first is a typical T Tauri star with mass, = 0.5 , radius, = 2.0 , and effective temperature, = 4000 K (Kenyon & Hartmann 1995). The second is a Herbig Ae star with mass, = 2.5 , radius, = 2.0 , and effective temperature, = 10,000 K (e.g., Palla & Stahler 1993).
The gas density profiles are determined by assuming vertical hydrostatic equilibrium, the balance between gravity and the pressure gradient force. In order to obtain the gas surface density profiles, we adopt the viscous -disk model (Shakura & Sunyaev 1973), with a viscous parameter, and a mass accretion rate, . The temperature profile of the gas is calculated self-consistently with the gas density profiles by assuming local thermal balance between gas heating and cooling. The thermal processes of heating from X-ray ionization of hydrogen, grain photoelectric heating induced by far-ultraviolet photons and cooling via gas-grain collisions and line transitions are taken into account. The dust temperature profiles are obtained by assuming local radiative equilibrium between the blackbody emission of grains and the absorption of radiation from the central star as well as the surrounding dust grains.
Regarding dust, we assume that the dust and gas are well-coupled and that the dust-to-gas mass ratio is constant (0.01) throughout the disk. The dust size distribution model in Weingartner & Draine (2001), which reproduces the observational extinction curve of dense clouds, is adopted. The dust properties are important factors for both gas and dust temperatures. The adopted dust opacity is described in Appendix D of Nomura & Millar (2005).
For the UV radiation field two radiation sources are taken into account: radiation from the central star and the interstellar medium. In the case of the T Tauri star, we use a UV excess model that reproduces the observational data toward TW Hydrae. The components of the stellar UV radiation are blackbody emission, hydrogenic thermal Bremsstrahlung emission and Ly line emission (see Appendix C of Nomura & Millar 2005 for details). For the Herbig Ae star, the UV excess contributes a significantly smaller fraction to the total UV luminosity compared with that from the stellar blackbody radiation (e.g., Dent et al. 2013). Therefore, we use only the blackbody emission of the central star as the source of UV radiation field for the Herbig Ae model. For the X-ray radiation, we adopt a TW-Hydrae-like spectrum with a total luminosity, , for the T Tauri star (Preibisch et al. 2005). Due to the non-convective interiors of Herbig Ae stars, the typical X-ray luminosities are 10 times lower than those of T Tauri stars (Güdel & Nazé 2009). Therefore, we assume an X-ray spectrum with and for the Herbig Ae star based on observations (e.g., Zinnecker & Preibisch 1994; Hamaguchi et al. 2005).
For full details on the generation of the physical models, see Nomura & Millar (2005) and Nomura et al. (2007). The total grid numbers are 8699 and 12116 for the disks around T Tauri and Herbig Ae stars, respectively, where 129 radial steps are taken logarithmically for the disk radius from 0.04 to 305 au. Figure 2 shows the physical structure of the disk around the T Tauri star and the Herbig Ae star (see also Walsh et al. 2015; Notsu et al. 2016, 2017). From top, the gas number density decreases as the function of height and the radius. The densest region reaches about near the disk midplane, falling to in the surface layer, showing the diversity of the environment in the PPDs. The gas and dust temperature decouple in the diffuse region due to inefficient collisions between the gas and dust grains. As the density increases toward the disk midplane, the gas and dust temperature become well coupled. The UV flux from the central star is much stronger in the case of Herbig Ae star (right-hand-side) than in the T-Tauri star (left-hand-side). However, both UV and X-ray photons are absorbed by dust and gas and do not penetrate to the midplane. Since the Herbig Ae star has a higher luminosity, its enhanced gas and dust temperatures make line emission easier to detect than in the T Tauri case (see Section 3.3). Noted that different radial ranges are shown in Figure 2 for each type of disk.
2.2 Chemical reaction network
We use a chemical reaction network, containing both gas-phase reactions and gas-grain interactions to study the chemical evolution of the disks. In order to examine the effect of carbon grain destruction (Lee et al. 2010; Anderson et al. 2017), we focus on those small molecules that are not much affected by grain surface reactions and consider two initial values for the elemental abundance of carbon. Except the abundance of carbon, for both models without and with carbon grain destruction, we use the typical low-metal value often adopted to model the chemistry of dark clouds such as TMC-1 (Table 8 of Woodall et al. 2007; see also Table 1). In the model without carbon grain destruction, the elemental abundance of carbon in gas is the same as that in diffuse clouds, in which the gas-phase elemental abundance of carbon is lower than that of oxygen. For the model with carbon grain destruction, we consider an extreme case assuming that all the carbon in the grains is destroyed and released into the gas-phase, and set the abundance for ionized carbon to be the solar abundance of carbon, (e.g., Asplund et al. 2009; Draine 2011). In this case, the gas-phase elemental abundance of carbon is larger than that of oxygen. We note that the C/O ratio is larger than unity in the gas-phase in this extreme case, although the C/O ratio of the solar photosphere is around 0.5 (Allende Prieto et al., 2002). This is because a large fraction of available elemental oxygen is incorporated into silicate grains which are more difficult to destroy than carbon grains. According to Lee et al. (2010), if the grain size is small enough, carbon grains could be sputtered into the gas phase in the surface layer of the disks where the gas temperature is high. If turbulent mixing in the vertical direction of the disk is efficient enough, a significant amount of the small carbon grains could be sputtered inside a disk radius of around 10 au in the case of T Tauri disks. Therefore, in this paper we only consider chemistry in the very inner region and treat the T Tauri disk up to 10 au, and the Herbig Ae disk, which the UV irradiation is stronger, up to 50 au. We calculate the chemical evolution up to years, which is the typical life time for protoplanetary disks and at which point the chemical structure has almost reached steady state in the disk. We note that for a more realistic model, the evolution of physical conditions, such as density, temperature, and radiation field, should be treated together with the time evolution of molecular abundances. This is beyond the scope of this work, and for simplicity, we use the chemical composition of diffuse clouds as an initial condition, and use the fixed physical conditions of each protoplanetary disk for the calculation of chemical reactions. Also, and again for simplicity, we ignore the change in the dust mass and surface area caused by the carbon grain destruction and the freeze-out of gas-phase molecules on grains during the calculation of chemical reactions.
|Without C grain destruction||With C grain destruction|
2.2.1 Gas-phase reactions
We use the RATE06 release of the UMIST Database for Astrochemistry, for the gas-phase chemistry (Woodall et al. 2007). In this chemical network, there are 375 species and 4336 reactions, including 3957 two-body reactions, 214 photoreactions, 154 X-ray/cosmic-ray-induced photoreactions, and 11 reactions of direct X-ray/cosmic-ray ionization. Since the reactions including fluorine, F, and phosphorus, P, have a low impact on the chemistry (Walsh et al. 2010), we remove reactions containing these two elements to reduce the computation time. For the photoreactions, we calculate the UV radiation at each point in the disk as
The UV radiation is scaled by the interstellar UV flux, (van Dishoeck et al. 2006), and photoreaction rates at each point in the disk can be approximated as
where is the unshielded photoreaction rates due to the interstellar UV radiation field as compiled in RATE06 (see also Millar et al. 1997).
2.2.2 Gas-grain interaction
The gas-grain interactions included are the adsorption (freeze-out) and desorption of molecules on and from dust grain surfaces, respectively. When the gas temperature is low enough, the adsorption rate becomes larger than the desorption rate, and therefore, molecules freeze onto dust grains. Otherwise, molecules sublimate from dust grains and into the gas-phase, so-called “thermal desorption”. The rate of thermal desorption, and thus the temperature at which sublimation happens, will depend on the binding energy of each molecule to dust grains. The binding energies of several important molecules and adopted in their work are listed in Table 2. In addition to thermal desorption, we consider non-thermal desorption mechanisms: cosmic-ray-induced thermal desorption (Leger et al. 1985; Hasegawa & Herbst 1993) and UV photodesorption (Westley et al. 1995; Willacy & Langer 2000; Öberg et al. 2007). The details of reaction rates for these processes can be found in the Appendix A (see also Walsh et al. 2010; Notsu et al. 2016).
Considering the processes of freeze out, thermal desorption, cosmic-ray-induced desorption and photodesorption, the differential equation for the number density of species on the dust grain is written as (e.g., Hasegawa et al. 1992)
where is the adsorption rate, is the thermal desorption rate, represents the cosmic-ray-induced thermal desorption rate, and denotes the photodesorption rate of species . is the gas phase number density of species and is the number density of species located in uppermost active layers of the ice mantle. The value of is given by (Aikawa et al., 1996) as
where is the total ice number density of all species (see also Walsh et al. 2014) and = represents the number of active surface sites in the ice mantle per unit volume and is the number of surface layers considered as active, assumed to be two. is the dust grain radius, represents the surface density of sites and is the number density of dust grains. For more detail, please refer to Appendix A.
2.3 Radiative transfer and the line emission
In order to predict observational signatures of carbon grain destruction in the inner region of protoplanetary disks, we perform ray-tracing calculations to estimate the intensity of molecular line emission and obtain intensity maps for both models with and without carbon grain destruction. The line data are taken from the Leiden Atomic and Molecular Database (LAMDA, Schöier et al. 2005), the Cologne Database for Molecular Spectroscopy (CDMS, Müller et al. 2005), or the Jet Propulsion Laboratory (JPL) molecular spectroscopic database (Pickett et al. 1998). We modify the original 1D code, RATRAN (Hogerheijde & van der Tak 2000) to calculate the ray-tracing using an axisymmetric 2D disk structure under the assumption of local thermodynamic equilibrium (LTE).
The intensity of the line profile at a frequency , , is obtained by solving the radiative transfer equation along the line-of-sight in the disk
where is the total extinction coefficient of dust grains and molecular lines and is the source function given by
respectively, where and are the Einstein A and B coefficients for spontaneous and stimulated emission and is the Einstein B coefficient for absorption. and are the number densities of the molecules in the upper and lower levels, respectively. is the mass density of dust grains where we simply apply the dust to gas mass ratio of . is the dust absorption coefficient at a frequency and is Plank’s constant.
Here is the line profile function which is affected by thermal broadening and Keplerian rotation in the disk. As a result, , is given by
where is the Doppler width, is the speed of light, represents the temperature of gas and is the mass of the species, . denotes the Doppler shift due to the projected Keplerian velocity along the line-of-sight and is given by
where is the gravitational constant, is the mass of the central star, is the distance of the line emitting region from the central star, is the azimuthal angle between the semimajor axis and the line linking the point in the disk along the line of sight and the center of the disk, and is the inclination angle of the disk.
In order to predict the observable flux density, we integrate Equation 5 along the line of sight . The intensity at in the projected plane is given by
where is the emissivity at the position and frequency , given by
and is the optical depth from the line emitting point to the disk surface at the frequency , expressed by
The observable line flux integrated all over the disk is given by
where is the distance between the observer and the target object.
3 Result and Discussion
3.1 Effect of carbon grain destruction on molecular abundance distribution
In this subsection, we show effect of carbon grain destruction on carbon-bearing species (HCN, \ceCH4, \ceC2H2, carbon-chain and hydrocarbon molecules) and oxygen-bearing species (\ceH2O, OH, \ceO2 and \ceCO2), some of which have been detected in protoplanetary disks at infrared wavelengths (e.g., Carr et al. 2004; Lahuis et al. 2006; Carr & Najita 2008; Salyk et al. 2008; Pascucci et al. 2009; Gibb et al. 2007). Here we mainly focus on the T Tauri disk model as an analogue of our solar system.
First, in order to see how the physical properties affect the molecular abundance profiles, the top panels of Figure 3 show the number density of hydrogen nuclei, the gas temperature and dust temperature as a function of height divided by radius (Z/r) at a disk radius of 1 au (left), 3 au (middle) and 10 au (right). The bottom panels show similar plots for the UV flux, the X-ray ionization rate and the cosmic-ray ionization rate.
Since gas-phase molecular abundances change dramatically across their snowlines, we estimate the location of the snowlines of some molecules by equating the adsorption rate (Equation A1) with the desorption rate (Equation A2). Because thermal desorption rate is related to the binding energy, Table 2 presents the location of the snowline of each species together with binding energies adopted in the model.
|Molecule||Binding Energy (K)||Snowline in T-Tauri Disk (au)||Snowline in HAe Disk (au)|
Figures 4, 6, 8 and 10 show molecular abundances relative to total hydrogen nuclei in both the gas-phase and ice mantle at years as a function of (Z/r) at the disk radii of 10 au, 3 au and 1 au. The dashed and the solid lines represent the case without and with carbon grain destruction, respectively. Figures 5, 7, 9 and 11 display the 2-dimensional fractional abundances of gas-phase and ice mantle molecules as a function of radius and height, out to a maximum radius of 10 au since the destruction of carbon grains mainly occurs in the inner hot region only. In the latter set of figures, the left-hand-side plots represent the case without carbon grain destruction (CO1) and right-hand-side plots those with carbon grain destruction (CO1).
Since the abundance of gas-phase CO has a large effect on molecular abundances differences between the models with and without the carbon grain destruction, as is shown below, we present first the gas-phase and ice mantle abundances of CO (Figure 4). In the disk surface, CO gas is dissociated by the FUV photons so that oxygen and carbon are not locked in CO but are mostly in the form of atoms and ions since the synthesis of molecules is difficult due to the high flux of FUV photons. Figure 5 shows the global CO abundance distribution, and CO is abundant (x(CO)=) throughout the whole region near the midplane, since it is easily desorbed inside 10 au due to its low binding energy. The abundance of CO increases slightly when carbon grains are destroyed, due to the additional elemental carbon now available in the gas phase.
The abundances of other molecules, however, are significantly affected by the gas-phase elemental C/O abundance ratio, especially in the region where CO gas is abundant. In general, carbon-bearing molecules are not expected to be very abundant in an oxygen-rich environment (C/O1) since most of the carbon is incorporated into CO. However, in the carbon grain destruction regions, excess carbon exists (C/O1) and carbon-bearing molecular abundances can increase dramatically. This type of chemistry is similar to that observed and successfully modeled around carbon-rich AGB stars, such as IRC+10216 (e.g., Millar & Herbst 1994; Glassgold 1996).
Figure 6 shows 1D plots of the abundances of some carbon-bearing species (HCN, \ceCH4, \ceC2H2 and c-\ceC3H2). The left and right panels show the molecular abundances in gas and ice, respectively. Figure 7 shows the 2-D abundance distribution of the carbon-bearing species. The effects of carbon grain destruction, and the resulting chemistry, are different between the surface and the midplane. Near the midplane, carbon-bearing molecules can form efficiently via gas-phase reactions in the carbon-rich case. The most significant difference between two cases appears in HCN especially inside 2 au, near the HCN snowline. Due to its high binding energy (see Table 2), HCN can remain on dust grains beyond 2 au in the T Tauri disk. In the carbon-rich case, the peak gas-phase abundance of HCN is near the midplane of the inner region. Near the midplane inside 2 au, the maximum differences between two cases can reach 8 orders of magnitude. In the oxygen-rich case, the peak abundance is in the molecular layer. Therefore, HCN appears to be a good tracer of the effect of carbon grain destruction within its snowline, and we will focus on it in section 3.3. We note, however, that the abundance distribution of HCN could be affected by the initial abundance of species and the ingredients adopted in the chemical model (see e.g., Figure 14 of Walsh et al. 2015).
Similar behavior can be seen in \ceCH4 and \ceC2H2, but the differences are less significant. The snowlines of \ceCH4 and \ceC2H2 are located outside 10 au and around 9.3 au, respectively. The gas phase abundances of both molecules increase near the midplane inside the snowlines in the carbon-rich case. However, inside 2 au, most of the carbon is incorporated into HCN and long carbon-chain species in the gas-phase (see Figure 8 and Figure 9) and the abundances of \ceCH4 and \ceC2H2 gas decrease here. We note that \ceCH4 is more abundant than HCN and carbon-chain molecules in the gas-phase near the midplane inside 2 au in the oxygen-rich case. In the oxygen-rich case \ceCH4 gas has its peak abundance ( ) in the molecular layer. In the case of c-\ceC3H2, it shows a differences of about 7 orders of magnitude between two cases. Though it is less abundant than HCN, we treat c-\ceC3H2 as a possible tracer in section 3.3 as well. We note that though the reaction network used here treats c-\ceC3H2 and l-\ceC3H2 separately, we should be cautious that there are some uncertainties since it is rarely possible to identify the specific isomeric products in a laboratory reaction and we need to rely on some calculations about energetics of the product and so on. Nevertheless, the model calculations show quite good agreement with observations of c- and l-\ceC3H2 in IRC+10216 and TMC-1 (McElroy et al. 2013).
Figure 8 presents the 1-D plots of C(), CH() and CH() molecules, where we sum up the abundances of the molecules for or . Figure 9 shows 2-D gas and ice distributions of \ceC6H and \ceC9. We chose these two species as representative of carbon-chain molecules since they are abundant on grain surfaces (Table 4) and, indeed, most of the carbon-chain molecules have spatial distributions similar to these two species and to HCN. We note that we treat the carbon-chain molecules together since grain surface reactions, which may reduce the abundances of carbon-chain molecules, for example through hydrogenation, are not included here. From our calculations, we find that carbon-chain molecules with larger number of carbon atoms have larger abundances because the destruction of carbon-chain molecules by molecular ions is inefficient near the midplane of the disk where the density is high and the ionization degree is very low (). In the carbon-rich case, the gas-phase abundances increase inside the snowline near the midplane. Beyond the snowline, the ice-mantle abundances increase significantly. The binding energies of carbon chain molecules are relatively large and thus the molecules can remain on the dust grains until very close to the parent star ( 1.5 - 3 au).
We note that the gaseous carbon-chain molecules are abundant in the molecular layers beyond their snowlines in both cases. This region corresponds to the CO gas depletion layer. Figure 4 and Figure 5 show CO gas is depleted at au. Similar CO gas depletion can be seen in other disk chemical models (e.g., Walsh et al. 2010; Furuya & Aikawa 2014). Though the self and mutual shielding of CO photodissociation are not included in the model, the shielding factors are not so effective especially in the inner region of the disks (Walsh et al. 2012). In our model, the decrease is caused by photodissociation of CO due to FUV irradiation from the central star, and the subsequent formation and freeze-out of carbon-chain molecules onto grains. If the FUV field is strong enough, FUV photons dissociate CO into atomic species. Above the CO gas depletion area, the UV destruction is fast, and due to the low density, the adsorption rate on to grains is low. Therefore, any synthesized molecules are difficult to stick on the grains and will eventually convert to CO gas or, nearer the surface remain as atomic and ionized carbon. Moving deeper towards the disk midplane, where UV photons are extinguished by grains in the high density molecular layers, CO can remain in the gas-phase. A fraction of the carbon-chain molecules in ice mantle in the CO gas depleted layer is photodesorbed into gas, which produces the gaseous carbon-chain molecule rich layer seen in Figure 8 and Figure 9.
Figure 10 shows the 1-D plots of oxygen-bearing species (\ceH2O, \ceOH, \ceO2 and \ceCO2) in the gas phase and ice mantle. Figure 11 shows the 2-dimensional distribution of gas-phase oxygen-bearing species, \ceH2O and \ceCO2. The results are opposite in behavior to those of carbon-bearing species. They are less abundant when carbon grains are destroyed. In Figure 11, in the oxygen-rich case, the volatile water abundance is enhanced inside its snowline around 2 au. In the carbon-rich case, the gas-phase water abundance is very low inside 10 au and there is no clear water ice snowline present. The distributions of \ceCO2 gas show clearly the snowline around 5 au in both cases. In the oxygen-rich case, the abundance peak of appears near the midplane inside 2 au, whereas in the carbon-rich case, the peak value can only reach close to the central star. We note that the abundance of \ceO2 gas can be affected by the initial abundance setting, atomic or molecular (Eistrup et al. 2016; Eistrup & Walsh 2018)
Finally, we briefly summarize in which molecules the carbon is mainly incorporated in both cases. Excluding CO gas, in the oxygen-rich case, carbon mainly resides in \ceCO2 gas inside the HCN snowline, and as \ceCO2 gas and HCN ice beyond the HCN snowline. Beyond the \ceCO2 snowline, carbon mainly resides in \ceCO2 and HCN ice. A fraction of carbon forms \ceCH4 gas throughout the inner disk. Meanwhile, in the carbon-rich case, carbon mainly forms HCN gas and long, gas-phase, carbon-chain molecules inside their snowlines. Beyond it, carbon mainly ends up in long, carbon-chain ice and gaseous \ceCH4, while inside the \ceCO2 snowline, some carbon is converted to gas-phase \ceCO2.
3.2 Radial distribution of solid carbon fraction
To compare the observed carbon depletion gradient in the inner solar system (Figure 1) with the modeled values, we calculate the carbon fraction in grains relative to the solar abundance of silicon as a function of the disk radius,
We assume that the solar abundance of silicon is with respect to hydrogen nuclei (Asplund et al. 2009; Draine 2011). The value for the refractory carbon is assumed to be the difference between the solar abundance of carbon, with respect to hydrogen nuclei (Draine 2011), and the gas-phase carbon abundance in the local diffuse ISM, (Table 1). In the model with carbon grain destruction, all the carbon has been released from the refractory material and thus there is no carbon left in this form. The released refractory carbon will be incorporated into species either in the gas-phase or in the ice mantle and the abundance of ice mantle carbon is obtained from our model calculation described in Section 3.1. Table 3 shows the ratio of the carbon in the solid phase relative to the solar silicon at each disk radius.
|without destruction||with destruction||with des. and increase by 25%|
|without destruction||with destruction|
|10 au||CO HCN C||CH CH CH|
|5 au||HCN C HCOOH||CH CH CH|
|4 au||HCN HCOOH C||CH CH C|
|3 au||HCN HCOOH CO||CH CH C|
|2 au||HCOOH C CS||CH CH C|
|1 au||HCOOH CH CO||C CH C|
In the case without carbon grain destruction, the fraction does not change very much at different radii (Table 3), since it is imposed that the majority of the carbon is locked in refractory form. Figure 12 shows the percentages of the form of carbon for the models with and without carbon grain destruction at a disk radius of 3 au as an example. It shows that most of the carbon (75 %) is locked in refractory form in the case without carbon grain destruction. In this case the elemental abundance of oxygen is larger than that of carbon (C/O 1) in the gas-phase, and thus most of the remaining carbon (24%) is stored in the form of CO gas with less than 1% in the form of larger organic species. Since the refractory carbon and CO gas do not change significantly in abundance over 1 au r 10 au in this model, there is almost no fractional variation in this case. On the other hand, in the case with carbon grain destruction, the fraction varies slightly at 2 au r 10 au and suddenly drops at r = 1 au. In this case, all the carbon in refractory form released to the gas phase reacts to form a range of species. Since the elemental abundance of carbon is larger than that of oxygen in gas-phase (C/O 1), oxygen is mainly stored in CO in this case, capturing (see Table 1), the corresponding amount of carbon, that is, 60 % of the total carbon as CO gas (Figure 12). The remaining carbon is mainly in the form of carbon-chain molecules in the ice mantle at 2au, as is shown in Sect. 3.1 (see e.g., Figure 9 and Figure 10). Table 4 shows the most abundant carbon-bearing species in the ice mantle in the disk midplane at different radii. The carbon-bearing molecules in the ice mantle evaporate into the gas phase inside their snowlines and the locations of the snowlines depend on their binding energies (see Table 2). Therefore, the fraction of carbon in the solid-phase changes across these snowlines of carbon-bearing species. Since the most abundant carbon-bearing molecules in the ice mantle are carbon-chain molecules (see Table 4) whose snowlines are located around the disk radii of 2 au (see Table 2), the solid carbon fraction does not decrease until r 2 au. Noted here, as we mentioned before, there are some uncertainties in the abundances of carbon-chain molecules in the ice mantle because grain surface reactions are not included in our calculation. Once grain surface reactions are considered, hydrogenation reactions on the surface may lead to the formation of alkanes.
While the model results reproduce qualitatively the trend of carbon depletion observed in the solar system, the comparison between our result (Table 3) and the solar system data (Figure 1) shows some quantitative discrepancies. There is a relatively high carbon fraction in the asteroid belt (a few to an order of magnitude larger) and too much depletion at 1 au ( three orders of magnitudes smaller) compared to the observation values. This is possibly caused by the omission of, for example, grain surface reactions and/or turbulent mixing, which we have not yet taken into account in our chemical model and will alter the chemical structure of the disk and the partitioning of carbon between solid and gaseous forms. With grain surface reactions, a higher carbon fraction in the ice mantle might be produced by large carbon-bearing molecules formed from species stuck on the grain surface that react with other species (e.g., Hasegawa et al. 1992; Hasegawa & Herbst 1993; Garrod & Herbst 2006; Garrod et al. 2008; Walsh et al. 2014). On the other hand, turbulent mixing can bring the dust grains from the midplane up to the warm surface layer where the ice mantle can be photodissociated or thermally desorbed. This leads to a decrease in the solid carbon fraction, perhaps to a level consistent with that in the asteroid belt (e.g., Furuya et al. 2013; Furuya & Aikawa 2014; Semenov & Wiebe 2011). We note that the timescale of vertical mixing is yr, where represents the parameter for vertical mixing, similar to the parameter in the viscous disk (e.g., Furuya et al. 2013), and that the timescale is shorter in the inner region. In addition, uncertainties in the binding energies of molecules on grains will affect the solid carbon fraction in grains. Desorption rates depend exponentially on binding energies (see Equation A2). However, due to the lack of appropriate data for many binding energies, the values are still very uncertain. We therefore enlarged the value of binding energies by 25% to investigate how the solid carbon fraction changes as binding energies change. We found that the value at 1 au becomes comparable with the observed value on the Earth when binding energies increase about 25% (see Table 3). We note that the increased binding energies of carbon-chain molecules are close to the data in complied for the UMIST Database for Astrochemistry 111http://www.udfa.net/ (McElroy et al. 2013).
3.3 Prediction for ALMA observations
Our results in Section 3.1 suggest that if carbon grains are destroyed in the inner region of protoplanetary disks, it will affect the molecular abundance profiles, a process which could be tested by observing the molecular line emission. In this subsection, we perform ray tracing calculations, using the molecular abundance profiles obtained from our model calculation and determine whether or not carbon grain destruction in the inner disk can be tested observationally.
While molecular lines have been observed towards protoplanetary disks at infrared and (sub)millimeter wavelengths, we choose lines which are observable with ALMA. Since, as we have seen in Section 3.1, significant differences in molecular abundances between the models with and without carbon grain destruction appear only near the disk midplane where CO is dominant and not photodissociated. Therefore, (sub)millimeter lines are more suitable to trace potential differences because infrared lines are optically thick and trace only the disk surface layer. In addition, ALMA enables us to spatially resolve the inner region of the disk, which allows us to see the difference between the models more clearly. Among the molecular lines observable at (sub)millimeter wavelengths, we choose the HCN lines because (1) the HCN abundance distribution is very much affected by the carbon grain destruction (see Section 3.1), and (2) they are known to be strong towards protoplanetary disks (e.g., Dutrey et al. 1997; Qi et al. 2008; Öberg et al. 2010; Öberg et al. 2011; Chapillon et al. 2012; Huang et al. 2017). Three HCN isotopologues, DCN, \ceH^13CN and \ceHC^15N, have been detected towards protoplanetary disks (e.g., Qi et al. 2008; Huang et al. 2017; Guzmán et al. 2015; Guzmán et al. 2017). We concentrate here on \ceH^13CN because the DCN/HCN ratio is known to be sensitive to the temperature profile (e.g., Aikawa & Herbst 2001; Willacy 2007; Cleeves et al. 2016; Aikawa et al. 2018), while \ceH^13CN is less affected (e.g., Woods & Willacy 2009). Meanwhile, \ceH^13CN lines are stronger than \ceHC^15N lines (Guzmán et al. 2017) and it is easier to analyze its spatial distribution. The combination of HCN and its isotopologue lines is useful to test the carbon grain destruction model since HCN lines are optically thick and trace mostly surface layer of the disk where CO is photodissociated and the carbon grain destruction does not affect molecular abundances significantly. The HCN isotopologue lines are less optically thick and trace the lower layer of the disk where CO is not photodissociated and the effect of carbon grain destruction is significant.
Figure 13 displays the profiles of the HCN and \ceH^13CN J=4-3 lines at 354.5 and 345.3 GHz respectively, with the spectral resolution of 0.4 , for the T-Tauri disk, calculated by performing ray tracing calculations (see Section 3.2) using the obtained physical structure (Figure 2) and the HCN abundance profile (Figure 7). The HCN/\ceH^13CN abundance ratio is simply assumed to be the typical interstellar value of 70 (Qi et al. 2011). According to Woods & Willacy (2009), it could be larger by a factor of 2 in the disk surface where the effect of self-shielding is significant and different CO isotopologues are selectively photo-dissociated. But the assumption is reasonable in the region closer to the disk midplane. We assume the distance to the T-Tauri disk is 70 pc and its inclination angle is 30 degrees. The left- and right-hand sides of Figure 13 show the line profiles of HCN and \ceH^13CN, respectively. The line profiles are calculated only inside a radius of 5 au since the significant differences between two models appears mainly inside 2 au (Figure 7). Solid and dashed lines represent models with and without grain destruction, respectively. Figure 13 indicates that the line emission of \ceH^13CN shows more obvious differences between the two cases. However, the peak flux density is less than 1 mJy for the \ceH^13CN line for the model without the carbon grain destruction, and it is too weak to use it for testing the carbon grain destruction model even with ALMA.
For this reason, we have also modeled a Herbig Ae disk since the radiation from the central source is stronger and the hot region is larger (Figure 2). Thus, the carbon grain destruction region spreads to larger radii and the observational test will be easier. Figure 14 shows the distribution of HCN in the Herbig Ae disk as a function of radius up to 50 au and the disk height divided by the radius (Z/r). Gas-phase reactions produce HCN efficiently only in the very hot region inside a disk radius of r 5 au for the model without carbon grain destruction, while gas-phase HCN can be very abundant in the whole region inside the HCN snowline at 20 au for the model with carbon grain destruction. Therefore, a dramatic change appears appears at radii of 2 - 30 au between two models. In contrast, the significant differences in T-Tauri disk appear only very close to the central star (inside 2 au) and the flux density is too low to be observed.
Figure 15 is the zeroth moment map, that is, the integrated intensity map of the HCN and \ceH^13CN lines, calculated as
where is the intensity as a function of position, , and velocity, , given by Equation 10.
The physical structure (Figure 2) and the HCN abundance profile (Figure 14) of the Herbig Ae disk are used. We assume that the distance to the disk is 140 pc and its inclination angle is 30 degrees. The \ceH^13CN line intensity map shows the difference more clearly than the HCN line. This is because the \ceH^13CN line traces a relatively lower layer of the disk due to relatively lower optical depth than the HCN line. Figure 16 shows the optical depth of the HCN and \ceH^13CN lines along the vertical direction of the disk for the models without and with carbon grain destruction. The HCN line becomes optically thick up to a disk radius of 30 au even for the model without the carbon grain destruction, that is, the HCN line only traces the disk surface where the difference in the HCN abundance is not significant between the models (Section 3.1). On the other hand, the isotopologue \ceH^13CN line is optically less thick and can trace the HCN abundance difference near the midplane at a disk radii of r 30 au. Therefore, the difference in the emitting regions of intensity map between the HCN and \ceH^13CN lines is a good tracer for testing the carbon grain destruction.
Even in those cases for which spatially resolved observations are difficult, the profiles of the HCN and \ceH^13CN lines can be a good tracer as well. Figure 17 shows the line profiles of the HCN line (left-hand-side) and the \ceH^13CN line (right-hand-side). In both profiles, the solid line represents the model with carbon grain destruction and the dotted line is the model without. Like the T-Tauri disk model (Figure 13), the line emission of \ceH^13CN shows more substantial differences between the two models. The differences in the \ceH^13CN line range from high velocity (10 km/s, tracing the inner region) to low velocity (tracing the outer region). In contrast, the differences in HCN emission mainly exist at velocity of -5 to + 5 km/s.
Figure 18 is the normalized cumulative line flux as functions of disk radius (left-hand-side) and the velocity (right-hand-side). The former is calculated from the disk center to the disk radius of 25 au and the latter is calculated over the range 0 20 km/s. Each cumulative flux is normalized using the following equations,
where is defined as 20 km and as 25 au. The blue and green lines are models with and without carbon grain destruction, respectively, with HCN represented by the solid lines and \ceH^13CN the dashed lines. Comparing the results in the same color, the differences between two models can be easily distinguished by the ratio of HCN and \ceH^13CN. In left-hand-side figure, the green dashed line reaches 60 % of the cumulative flux inside 5 au and the green solid line reaches the same level inside 15 au. This can be explained from the intensity map. For the model without carbon grain destruction, the intensity of the \ceH^13CN line is strong only in very compact region ( 5 au) as shown in Figure 15. Meanwhile, the HCN line is strong out to 30 au and thus the normalized cumulative flux increases smoothly with radius. On the other hand, for the model with carbon grain destruction, both HCN and \ceH^13CN lines are strong out to 30 au and the cumulative fluxes increase smoothly, with no significant difference between them, so that differences between the two models can be tested by the cumulative flux as a function of radius if it can be spatially resolved. Even in the case for which the emission is not spatially resolved, we can see these differences in the line profiles, that is, cumulative flux as a function of velocity (e.g., Zhang et al. 2017). The green lines in Figure 18 (right-hand-side) show the differences between two models clearly. Comparing the normalized cumulative flux with the line profiles (Figure 17), the shape and the distribution of the model without carbon grain destruction are very different from HCN and \ceH^13CN. However, the line shape and the distribution in the model with carbon grain destruction are very similar for both HCN and \ceH^13CN. Both figures show significant differences between the HCN and \ceH^13CN lines for the model without the carbon grain destruction, while the difference is not large for the model with carbon grain destruction. Therefore, the differences between two models can be easily distinguished by the ratio of the HCN and \ceH^13CN line.
HCN and \ceH^13CN lines have been detected toward some Herbig Ae stars by ALMA (e.g., Guzmán et al. 2015, Huang et al. 2017). Guzmán et al. (2015) have reported the observations of the HCN J = 4 - 3 and \ceH^13CN J = 3 - 2 lines with relatively low spatial resolution of 0.7″. Mapping the inner region of the disk (around the disk radius of 10 - 20 au), taking advantage of ALMA’s high spatial resolution and high sensitivity, would be useful to diagnose carbon grain destruction in the inner disk.
Because the abundance distribution of molecules has some uncertainty depending on the chemical model, here we suggest c-\ceC_3H_2 - (217.882 GHz) as the other target line. This line has been detected in a Herbig Ae disk (Qi et al. 2013). The line is blended with the c-\ceC3H2 - line, but we treat only the \ceC_3H_2 - since it is slightly stronger. We note that other transition lines of c-\ceC3H2 have been detected toward T Tauri disks (Bergin et al. 2016). Figure 19 shows that the distribution of c-\ceC_3H_2 in the Herbig Ae disk as a function of radius up to 50 au and the disk height divided by the radius (Z/r), which indicates the significant difference appears at r 30 au between the models, similar to the case of HCN. The difference in the abundance distribution results in the clear differences in both the intensity map (Figure 20) and the line profile (Figure 21). To sum up briefly, we suggest \ceHCN, \ceH^13CN and c-\ceC_3H_2 as possible tracers of testing the carbon grain destruction effect in protoplanetary disks.
In this work, we focus on two issues: (1) The carbon depletion gradient in the inner solar system, and (2) searching for observational evidence of carbon grain destruction in the disks of T-Tauri and Herbig Ae stars. We conclude that the carbon grain destruction affects the abundances and distribution of various molecules in the protoplanetary disk, showing significant differences especially near the midplane in the inner region of the disk, where CO gas is abundant and is not photodissociated. For example, the gas-phase HCN abundance shows 8 orders of magnitudes difference near the midplane inside the radius of 2 au in the T-Tauri disk.
The distribution of molecules is determined by their volatility such that volatile species can evaporate into gas even in the outer region of the disk. Those molecules which remain in ice mantles can influence the composition of subsequently forming planets. Therefore, we present the location of the snowlines and calculate the solid carbon fraction relative to the solar abundance of carbon as a function of the radius in the inner region of a protoplanetary disk. Although the carbon depletion gradient is reproduced in the model with carbon grain destruction, the resulting solid carbon fractions show a quantitative discrepancy from those in the solar system. The solid carbon fractions in the asteroid belt are about an order of magnitude larger than the measured value of the meteorites. Meanwhile, the value at 1 au is 3 orders of magnitude smaller than the measured value. Including grain surface reactions in the model may help to better reproduce the carbon depletion gradients. For example, grain surface reactions are expected to make more complex, less volatile, organic molecules in the disk and could further enlarge the solid carbon fraction at 1 au. We also examined the effect of carbon grain destruction on predictions for ALMA observations. We find that lines in T Tauri are too weak to probe carbon grain destruction but that ALMA can probe this effect through the line ratio of HCN/\ceH^13CN as well as c-\ceC3H2 in Herbig Ae disks.
Appendix A The calculation of adsorption and desorption rates
The adsorption(freeze-out) rate of species i onto grain surface, , is written as (Hasegawa et al. 1992)
where is the sticking coefficient and we set it as 0.4 for all species (Veeraghattam et al. 2014), is the geometrical cross section of a dust grain, is the dust grain radius and is the number density of dust grains. We fix the number density ratio of dust grains to hydrogen nuclei () times the geometrical cross section of grain (), as , according to Rawlings et al. (1992), which is consistent with a gas-to-dust mass ratio of 100. represents the thermal velocity of species , where is the Boltzmann’s constant, is the temperature of gas and is the mass of species .
The thermal desorption rate, , with which species evaporate from the dust grain surface, is represented as (Hasegawa et al. 1992)
where is the binding energy of species to the dust grain surface in units of Kelvin, is the dust temperature and is the vibrational frequency of each adsorbed species in its potential well (Hasegawa et al. 1992).
where represents the surface density of sites, , and is the mass of the adsorbed species .
For cosmic-ray induced thermal desorption, we assume that dust grains with a radius of 0.1 are heated by the impact of relativistic Fe nuclei with energy from 20 to 70 MeV and deposit an energy of 0.4 MeV on average into each dust grain (Leger et al. 1985; Hasegawa & Herbst 1993). Assuming that the majority of molecules will desorb around 70 K, the cosmic ray induced thermal desorption rate, , is expressed as
where is the cosmic ray ionization rate of \ceH2 scaled by the interstellar value used in the UMIST database, . is the thermal desorption rate of species at dust temperature of 70 K. is the fraction of time that dust grain spends above 70 K and it is roughly calculated by the ratio of the desorption cooling time ( s) to the total interval time for the temperature of dust grain to become 70 K ( s, Leger et al. 1985). Therefore, f(70 K) . Although X-rays can penetrate the disk and induce desorption as well, we do not include it in this chemical network because of the remaining uncertainty in the desorption rate.
Photodesorption is independent of the surface binding energy. The photodesorption rate adopted is based on the experiments of Westley et al. (1995) and Öberg et al. (2007). Their results show that each photon absorbed by a dust grain will release a particular number of molecules into the gas phase and it is related to the fractional abundance on the dust grain surface. The photodesorption rate, , is expressed by the following equation (Willacy & Langer 2000; Willacy 2007)
where represents the UV radiative field at each position in the disk in units of photon . is the photodesorption yield determined from the experiments in the unit of molecules , and we adopt the value of , which is determined by experiments of pure water ice (Westley et al. 1995) and pure CO ice (Öberg et al. 2007). = represents the number of active surface sites in the ice mantle per unit volume and is the number of surface layers considered as active, assumed to be two.
- Aikawa et al. (2018) Aikawa, Y., Furuya, K., Hincelin, U., & Herbst, E. 2018, ApJ, 855, 119
- Aikawa & Herbst (2001) Aikawa, Y., & Herbst, E. 2001, A&A, 371, 1107
- Aikawa et al. (1996) Aikawa, Y., Miyama, S. M., Nakano, T., & Umebayashi, T. 1996, ApJ, 467, 684
- Allègre et al. (2001) Allègre, C., Manhès, G., & Lewin, É. 2001, E&PSL, 185, 49
- Allende Prieto et al. (2002) Allende Prieto, C., Lambert, D. L., & Asplund, M. 2002, ApJ, 573, L137
- Anderson et al. (2017) Anderson, D. E., Bergin, E. A., Blake, G. A., et al. 2017, ApJ, 845, 13
- Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481
- Bardyn et al. (2017) Bardyn, A., Baklouti, D., Cottin, H., et al. 2017, MNRAS, 469, S712
- Bergin et al. (2007) Bergin, E. A., Aikawa, Y., Blake, G. A., & van Dishoeck, E. F. 2007, 2007prpl conf, 751
- Bergin et al. (2015) Bergin, E. A., Blake, G. A., Ciesla, F., Hirschmann, M. M., & Li, J. 2015, PNAS, 112, 8965
- Bergin et al. (2016) Bergin, E. A., Du, F., Cleeves, L. I., et al. 2016, ApJ, 831, 101
- Bergner et al. (2018) Bergner, J. B., Guzmán, V. G., Öberg, K. I., Loomis, R. A., & Pegues, J. 2018, ApJ, 857, 69
- Carr & Najita (2008) Carr, J. S., & Najita, J. R. 2008, Sci, 319, 1504
- Carr et al. (2004) Carr, J. S., Tokunaga, A. T., & Najita, J. 2004, ApJ, 603, 213
- Chapillon et al. (2012) Chapillon, E., Guilloteau, S., Dutrey, A., Piétu, V., & Guélin, M. 2012, A&A, 537, A60
- Cleeves et al. (2016) Cleeves, L. I., Bergin, E. A., O’D. Alexander, C. M., et al. 2016, ApJ, 819, 13
- Dent et al. (2013) Dent, W. R. F., Thi, W. F., Kamp, I., et al. 2013, PASP, 125, 477
- Draine (2003) Draine, B. T. 2003, ARA&A, 41, 241
- Draine (2011) Draine, B. T. 2011, Physics of the interstellar and intergalactic medium (Princeton University Press)
- Dutrey et al. (1997) Dutrey, A., Guilloteau, S., & Guelin, M. 1997, A&A, 317, L55
- Dutrey et al. (2014) Dutrey, A., Semenov, D., Chapillon, E., et al. 2014, 2014prpl conf, 317
- Eistrup & Walsh (2018) Eistrup, C., & Walsh, C. 2018, A&A, in press (arXiv:1808.03329)
- Eistrup et al. (2016) Eistrup, C., Walsh, C., & van Dishoeck, E. F. 2016, A&A, 595, A83
- Favre et al. (2018) Favre, C., Fedele, D., Semenov, D., et al. 2018, ApJ, 862, L2
- Finocchi et al. (1997) Finocchi, F., Gail, H.-P., & Duschl, W. J. 1997, A&A, 325, 1264
- Furuya & Aikawa (2014) Furuya, K., & Aikawa, Y. 2014, ApJ, 790, 97
- Furuya et al. (2013) Furuya, K., Aikawa, Y., Nomura, H., Hersant, F., & Wakelam, V. 2013, ApJ, 779, 11
- Gail (2001) Gail, H.-P. 2001, A&A, 378, 192
- Gail (2002) —. 2002, A&A, 390, 253
- Garrod & Herbst (2006) Garrod, R. T., & Herbst, E. 2006, A&A, 457, 927
- Garrod et al. (2008) Garrod, R. T., Widicus Weaver, S. L., & Herbst, E. 2008, ApJ, 682, 283
- Geiss (1987) Geiss, J. 1987, A&A, 187, 859
- Gibb et al. (2007) Gibb, E. L., Van Brunt, K. A., Brittain, S. D., & Rettig, T. W. 2007, ApJ, 660, 1572
- Glassgold (1996) Glassgold, A. E. 1996, ARA&A, 34, 241
- Güdel & Nazé (2009) Güdel, M., & Nazé, Y. 2009, A&A Rev., 17, 309
- Guzmán et al. (2017) Guzmán, V. V., Öberg, K. I., Huang, J., Loomis, R., & Qi, C. 2017, ApJ, 836, 30
- Guzmán et al. (2015) Guzmán, V. V., Öberg, K. I., Loomis, R., & Qi, C. 2015, ApJ, 814, 53
- Hamaguchi et al. (2005) Hamaguchi, K., Yamauchi, S., & Koyama, K. 2005, ApJ, 618, 360
- Hasegawa & Herbst (1993) Hasegawa, T. I., & Herbst, E. 1993, MNRAS, 263, 589
- Hasegawa et al. (1992) Hasegawa, T. I., Herbst, E., & Leung, C. M. 1992, ApJS, 82, 167
- Henning & Semenov (2013) Henning, T., & Semenov, D. 2013, ChRv, 113, 9016
- Hogerheijde & van der Tak (2000) Hogerheijde, M. R., & van der Tak, F. F. S. 2000, A&A, 362, 697
- Huang et al. (2017) Huang, J., Öberg, K. I., Qi, C., et al. 2017, ApJ, 835, 231
- Jessberger et al. (1988) Jessberger, E. K., Christoforidis, A., & Kissel, J. 1988, Nature, 332, 691
- Jura (2006) Jura, M. 2006, ApJ, 653, 613
- Kenyon & Hartmann (1995) Kenyon, S. J., & Hartmann, L. 1995, ApJS, 101, 117
- Lahuis et al. (2006) Lahuis, F., van Dishoeck, E. F., Boogert, A. C. A., et al. 2006, ApJ, 636, L145
- Lee et al. (2010) Lee, J.-E., Bergin, E. A., & Nomura, H. 2010, ApJ, 710, L21
- Leger et al. (1985) Leger, A., Jura, M., & Omont, A. 1985, A&A, 144, 147
- Lodders (2010) Lodders, K. 2010, ASSP, 16, 379
- Loomis et al. (2018) Loomis, R. A., Cleeves, L. I., Öberg, K. I., et al. 2018, ApJ, 859, 131
- McElroy et al. (2013) McElroy, D., Walsh, C., Markwick, A. J., et al. 2013, A&A, 550, A36
- Millar et al. (1997) Millar, T. J., Farquhar, P. R. A., & Willacy, K. 1997, A&AS, 121, 139
- Millar & Herbst (1994) Millar, T. J., & Herbst, E. 1994, A&A, 288, 561
- Müller et al. (2005) Müller, H. S., Schlöder, F., Stutzki, J., & Winnewisser, G. 2005, JMoSt, 742, 215
- Nomura et al. (2007) Nomura, H., Aikawa, Y., Tsujimoto, M., Nakagawa, Y., & Millar, T. J. 2007, ApJ, 661, 334
- Nomura & Millar (2005) Nomura, H., & Millar, T. J. 2005, A&A, 438, 923
- Notsu et al. (2016) Notsu, S., Nomura, H., Ishimoto, D., et al. 2016, ApJ, 827, 113
- Notsu et al. (2017) —. 2017, ApJ, 836, 118
- Öberg et al. (2007) Öberg, K. I., Fuchs, G. W., Awad, Z., et al. 2007, ApJ, 662, L23
- Öberg et al. (2015) Öberg, K. I., Guzmán, V. V., Furuya, K., et al. 2015, Natur, 520, 198
- Öberg et al. (2010) Öberg, K. I., Qi, C., Fogel, J. K. J., et al. 2010, ApJ, 720, 480
- Öberg et al. (2011) —. 2011, ApJ, 734, 98
- Palla & Stahler (1993) Palla, F., & Stahler, S. W. 1993, ApJ, 418, 414
- Pascucci et al. (2009) Pascucci, I., Apai, D., Luhman, K., et al. 2009, ApJ, 696, 143
- Pickett et al. (1998) Pickett, H., POYNTER, R., COHEN, E., et al. 1998, JQSRT, 60, 883
- Pontoppidan et al. (2014) Pontoppidan, K. M., Salyk, C., Bergin, E. A., et al. 2014, 2014prpl conf, 363
- Pontoppidan et al. (2010) Pontoppidan, K. M., Salyk, C., Blake, G. A., et al. 2010, ApJ, 720, 887
- Preibisch et al. (2005) Preibisch, T., Kim, Y.-C., Favata, F., et al. 2005, ApJS, 160, 401
- Qi et al. (2011) Qi, C., D’Alessio, P., Öberg, K. I., et al. 2011, ApJ, 740, 84
- Qi et al. (2013) Qi, C., Öberg, K. I., Wilner, D. J., & Rosenfeld, K. A. 2013, ApJ, 765, L14
- Qi et al. (2008) Qi, C., Wilner, D. J., Aikawa, Y., Blake, G. A., & Hogerheijde, M. R. 2008, ApJ, 681, 1396
- Rawlings et al. (1992) Rawlings, J. M. C., Hartquist, T. W., Menten, K. M., & Williams, D. A. 1992, MNRAS, 255, 471
- Salyk et al. (2008) Salyk, C., Pontoppidan, K. M., Blake, G. A., et al. 2008, ChRv, 676, L49
- Savage & Sembach (1996) Savage, B. D., & Sembach, K. R. 1996, ARA&A, 34, 279
- 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
- Semenov & Wiebe (2011) Semenov, D., & Wiebe, D. 2011, ApJS, 196, 25
- Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
- Siebenmorgen & Krügel (2010) Siebenmorgen, R., & Krügel, E. 2010, A&A, 511, A6
- Thi et al. (2004) Thi, W.-F., van Zadelhoff, G.-J., & van Dishoeck, E. F. 2004, A&A, 425, 955
- van Dishoeck et al. (2006) van Dishoeck, E. F., Jonkheid, B., & van Hemert, M. C. 2006, FaDi, 133, 231
- van Dishoeck et al. (2003) van Dishoeck, E. F., Thi, W.-F., & van Zadelhoff, G.-J. 2003, A&A, 400, L1
- Veeraghattam et al. (2014) Veeraghattam, V. K., Manrodt, K., Lewis, S. P., & Stancil, P. C. 2014, ApJ, 790, 4
- Walsh et al. (2010) Walsh, C., Millar, T. J., & Nomura, H. 2010, ApJ, 722, 1607
- Walsh et al. (2014) Walsh, C., Millar, T. J., Nomura, H., et al. 2014, A&A, 563, A33
- Walsh et al. (2012) Walsh, C., Nomura, H., Millar, T. J., & Aikawa, Y. 2012, ApJ, 747, 114
- Walsh et al. (2015) Walsh, C., Nomura, H., & van Dishoeck, E. 2015, A&A, 582, A88
- Walsh et al. (2016) Walsh, C., Loomis, R. A., Öberg, K. I., et al. 2016, ApJ, 823, L10
- Weingartner & Draine (2001) Weingartner, J. C., & Draine, B. T. 2001, ApJ, 548, 296
- Westley et al. (1995) Westley, M. S., Baragiola, R. A., Johnson, R. E., & Baratta, G. A. 1995, P&SS, 43, 1311
- Willacy (2007) Willacy, K. 2007, ApJ, 660, 441
- Willacy & Langer (2000) Willacy, K., & Langer, W. D. 2000, ApJ, 544, 903
- Woodall et al. (2007) Woodall, J., Agúndez, M., Markwick-Kemper, A. J., & Millar, T. J. 2007, A&A, 466, 1197
- Woods & Willacy (2009) Woods, P. M., & Willacy, K. 2009, ApJ, 693, 1360
- Zhang et al. (2017) Zhang, K., Bergin, E. A., Blake, G. A., Cleeves, L. I., & Schwarz, K. R. 2017, NatAs, 1, 0130
- Zinnecker & Preibisch (1994) Zinnecker, H., & Preibisch, T. 1994, A&A, 292, 152