Exclusion of Cosmic Rays in Disks. II.

Exclusion of Cosmic Rays in Protoplanetary Disks. II. Chemical Gradients and Observational Signatures

L. Ilsedore Cleeves1 , Edwin A. Bergin1 and Fred C. Adams1 2
1affiliation: Department of Astronomy, University of Michigan, 1085 S. University Ave, Ann Arbor, MI 48109
2affiliation: Department of Physics, University of Michigan, 450 Church St, Ann Arbor, MI 48109

The chemical properties of protoplanetary disks are especially sensitive to their ionization environment. Sources of molecular gas ionization include cosmic rays, stellar X-rays and short-lived radionuclides, each of which varies with location in the disk. This behavior leads to a significant amount of chemical structure, especially in molecular ion abundances, which is imprinted in their submillimeter rotational line emission. Using an observationally motivated disk model, we make predictions for the dependence of chemical abundances on the assumed properties of the ionizing field. We calculate the emergent line intensity for abundant molecular ions and simulate sensitive observations with the Atacama Large Millimeter/Sub-millimeter Array (ALMA) for a disk at  pc. The models readily distinguish between high ionization rates ( s per H) and below, but it becomes difficult to distinguish between low ionization models when  s. We find that HD emission is not detectable for sub-interstellar CR rates with ALMA (6h integration), and that ND emission may be a more sensitive tracer of midplane ionization. HCO traces X-rays and high CR rates ( s), and provides a handle on the warm molecular ionization properties where CO is present in the gas. Furthermore, species like HCO, which emits from a wide radial region and samples a large gradient in temperature, can exhibit ring-like emission as a consequence of low-lying rotational level de-excitation near the star. This finding highlights a scenario where rings are not necessarily structural or chemical in nature, but simply a result of the underlying line excitation properties.

Subject headings:
astrochemistry — circumstellar matter — protoplanetary disks — radiative transfer — stars: pre-main sequence

1. Introduction

Ionization plays an important role in setting thermal (e.g., Glassgold et al., 2004), dynamical (Balbus & Hawley, 1991; Gammie, 1996; Matsumura & Pudritz, 2003), and chemical (e.g., Semenov et al., 2004; Öberg et al., 2011b) properties of protoplanetary disks. The chemistry occurring in the bulk mass of disks is especially sensitive to ionization for two reasons: (i) in the cold gas, ion-neutral reactions are the most efficient means towards chemical complexity, and (ii) in the ices, the crucial hydrogenation reactions (Tielens & Hagen, 1982; Hasegawa et al., 1992) are fueled by hydrogen atoms that are extracted by ionization of molecular H. The dominant ionization processes in disks are photoionization from stellar and interstellar UV and X-ray radiation, thermal ionization, ionization by the decay products of short-lived radionuclides (SLRs), and cosmic ray (CR) ionization. Additionally, the cluster environment can provide a source of ionization on the outer surface of the disk from interstellar FUV, which can exceed that of the galactic average interstellar FUV by a factor of (Fatuzzo & Adams, 2008). Of these sources, CRs and potentially SLR decay provide ionization in the densest and coldest layers of the disk, where UV and X-rays are highly attenuated. However, it is unclear as to whether or not CRs are actually present at rates derived for the interstellar medium (ISM), i.e.,  s in dense gas and  s in the diffuse ISM (Dalgarno, 2006, and references therein). Modulation of the CR flux can occur as a result of exclusion by natal stellar winds as explored in detail in Cleeves et al. (2013a, hereafter Paper I) and discussed in Glassgold (1999); Aikawa & Herbst (1999a) and Fromang et al. (2002), in addition to exclusion by magnetic fields (Dolginov & Stepinski, 1994; Padovani & Galli, 2011; Fatuzzo & Adams, 2014). At these reduced levels, the ionization from SLR decay products becomes as important, and perhaps even more so, than that of CRs.

In the present work, we set out to determine the chemical imprint of individual ionization processes in a generic protoplanetary disk model. We outline how observations of molecular species can be used as a blueprint to constrain the underlying ionization properties of the disk’s dense molecular gas. More specifically, we focus on the impact of the assumed CR flux on molecular abundances in tandem with a detailed treatment of ionization by SLRs and stellar X-rays, including a Monte Carlo treatment of the scattered X-ray radiation field. We extend these results to make testable predictions for future sensitive observations with the Atacama Large Millimeter/Submillimeter Array (ALMA) of molecular ion emission in protoplanetary disks. Such predictions will help more accurately determine not only the ionization fraction in disks, which is important for constraining models of turbulence and chemistry, but also the source of ionization traced by a given molecular ion and transition.

Current detections of molecular ions in disks include NH, HCO, and DCO (e.g., Dutrey et al., 1997; van Dishoeck et al., 2003; Thi et al., 2004; Öberg et al., 2010, 2011a). The elusive HD has yet to be detected in a disk, with reported detections of the () line towards DM Tau and TW Hya (Ceccarelli et al., 2004) not confirmed by reanalysis or more sensitive observations (Guilloteau et al., 2006; Qi et al., 2008; Chapillon et al., 2011). The HD upper limits nonetheless provide some constraints on the midplane ionization rate. Chapillon et al. (2011) investigated chemical models incorporating deuterium chemistry and find models with midplane ionization rates below  s are required to explain the line’s non-detection. Strictly speaking, CRs are the primary midplane ionizing agent considered by Chapillon et al. (2011), though we note this value more generally provides a limit on the total ionization rate, which is expected to be primarily due to SLRs and CRs in the HD emitting gas, discussed in more detail below. We note that the deuterium chemistry of HD, specifically the reactions leading back to H, are sensitive to the ortho-to-para ratio of both H and H (Chapillon et al., 2011; Hugo et al., 2009; Albertsson et al., 2014), which we include in this work (see §4.1). Furthermore, HD can undergo deuterium substitution towards the fully substituted and unobservable D (Walmsley et al., 2004), which complicates the interpretation of HD measurements in determining ionization from HD. Such progressive D-substitution may point towards future difficulty in detection experiments of what would otherwise be a useful tracer of ionization and cold chemistry. This point further emphasizes the utility of chemical models in the interpretation of molecular ion emission as a proxy for measuring disk ionization rates.

The paper is structured as follows: In Sections §2 and §3, we review the physical model and the sources of ionization considered. In Section §4.1 we describe the updated chemical reaction network used to predict the molecular abundances, and Section §4.2 examines the chemical abundance variations between the different ionization models. Section §5 makes predictions for sensitive, high-spatial resolution submillimeter observations that can be used to determine and to distinguish between the important ionization mechanisms. In Section §6, we discuss the effects of a more X-ray luminous source, time-decay of radionuclides, disk gas mass, and the impact of disk magnetic fields. In Section §7, we summarize our results and discuss their implications.

Figure 1.— Disk model used for the chemical calculations. Top row, left to right: i) gas density, ii) settled dust density, iii) dust temperature. Bottom row, left to right: i) gas temperature, ii) integrated X-ray radiation field and iii) FUV radiation field from the central star. The FUV flux is integrated between  Å using the observed TW Hydra FUV spectrum (Herczeg et al., 2002, 2004), including Ly-. The X-ray luminosity is  erg s integrated between  keV.

2. The Physical Model

We refer the reader to a detailed description of the disk model parameters in Paper I, and review only the main aspects of the model below. The fiducial disk mass is  M within a  AU radius, where the vertical structure and geometrical parameters (i.e. flaring) are typical of observed protoplanetary disks (Andrews et al., 2011). The gas and dust surface densities follow a power law with an exponential taper at the outer edge (Hughes et al., 2008; Andrews et al., 2011) and the dust is heated via passive irradiation from the central K5V star with T = 4300 K with the code TORUS (Harries, 2000; Harries et al., 2004; Kurosawa et al., 2004; Pinte et al., 2009). The gas temperature is calculated from a fitting function calibrated to detailed thermochemical models of disks heated by FUV photons from the central star (Simon Bruderer in private communication, 2013). In the heating calculation, we consider only the central star’s FUV field (described in §3.1). The vertical distribution of dust is broken up into two populations to simulate grain growth, with a population of large millimeter grains concentrated at the geometrical center of the disk (the midplane) and a population of “atmosphere” micron-sized grains distributed over larger scale heights. The gas and small micron-sized grains follow the same scale height (see Paper I for details, as well as the dust optical properties). Figure 1 shows the disk density, thermal and stellar radiation fields, and Table 1 outlines the main physical parameters of our model.

Stellar Model
Stellar Mass 1.06 M
Stellar Radius 1.83 R
Stellar T 4300 K
L 2.9 erg/s aaAs computed from the observed FUV spectrum of TW Hya integrated between 930 and 2000 Å (Herczeg et al., 2002, 2004).
L 10 erg/s
Disk Model
R 0.1 AU
R 400 AU
M 3.9 M
M 0.039 M
Table 1Stellar and disk model parameters.

3. Ionization Sources

Star/disk systems are subject to a variety of ionization sources, including FUV and X-ray radiation from the central stars, short-lived radionuclides, and CRs. These ionization sources are not constant in time or in space. Instead, they are expected to vary with the local environment, from system to system, and display time dependence. Some environments include ionizing radiation from nearby (more massive) stars (see the review of Adams, 2010), although this complication is left for future work. This section outlines the physical mechanisms that contribute to the disk ionization rate and the ranges of values considered at present.

3.1. Stellar Photoionization

T Tauri stars are known to be X-ray luminous with X-ray luminosities typically ranging  erg s (e.g., Feigelson & Decampli, 1981; Feigelson et al., 1993; Telleschi et al., 2007). Consequently, the disk’s X-ray ionization properties are perhaps the best (observationally) constrained ionization agent present in the disk molecular gas. That said, the permeability of the disk gas to X-rays is unknown owing to uncertainties in the disk mass (the column density) and composition (opacity) of the absorbing materials. The biggest uncertainty in the X-ray opacity is in the time-dependent effects of dust settling, which redistributes and removes the absorbing material from the upper layers. More specifically, by removing dust mass from the upper layers the opacity shifts from a well-mixed gas and dust regime to a gas-dominated opacity. For Asplund et al. (2009) abundances, the change in opacity corresponds to a factor of two decrease in absorption cross section between the well-mixed and fully settled cases at  keV (Bethell & Bergin, 2011a). However, with knowledge of the star’s X-ray luminosity, it may be possible to shed light on the permeability of the disk atmosphere to X-rays with the proposed molecular tracers in this work (see §4.2).

We calculate the stellar FUV and X-ray radiation fields as a function of position and wavelength within the disk using a Monte Carlo treatment described in Bethell & Bergin (2011b). For FUV transfer, we consider the dust model’s position-dependent opacities and compute the absorption and scattering on dust grains. In addition to the radiative transfer through the dust, we calculate the Lyman- line transfer, where photons scatter isotropically off atomic hydrogen present at the disk surface. Treatment of Lyman- is important as this line has been observed to carry of the total FUV flux (Herczeg et al., 2004; Schindhelm et al., 2012; France et al., 2014); as a consequence of its scattering properties, the Lyman- photons are able to penetrate deeper into the disk gas than the rest of the primarily dust-scattered FUV photons (Bethell & Bergin, 2011b). Because the present paper is mainly focused on understanding the ionization properties of the dense molecular gas, we hold the FUV radiation field constant in the models presented here since FUV mainly contributes to the ionization in the lower density surface layers where, for example, C emission originates. The assumed incident FUV spectrum is that of TW Hydra (Herczeg et al., 2002, 2004), the closest and least extincted T Tauri star along the line of sight,  pc (Perryman et al., 1997). In this work, we do not include the interstellar FUV radiation field as a source of ionization in the dense gas, but it is expected to play a role in the cluster environment where the external FUV field can be many thousands of (Fatuzzo & Adams, 2008), where is the mean value for the ISM (Habing, 1968; Draine, 1978). A =1 has a less significant effect on the chemistry since a small amount of UV opacity from small grains restricts the UV penetration to an outer shell on the disk surface, which is dwarfed by scattered UV from star itself (see Fig. 3 of Cleeves et al., 2013a). A similar ionization study should be carried out for more extreme cluster environments, where the high UV field will be accompanied by a higher SLR abundance and potentially CR abundance in the vicinity of massive stars.

The specific (energy-dependent) X-ray fluxes are calculated using the same code as was used for the FUV (Bethell & Bergin, 2011b) where we instead adopt the updated X-ray absorption gas and dust opacities of Bethell & Bergin (2011a). We note that the model has been updated since Paper I, where the dust-reduced opacity in the upper layers (90% reduction in dust, or “e0p1”) was originally assumed to be uniform throughout the disk. We have since updated this calculation to take into account our knowledge of the local gas-to-dust mass ratio. We determine the absorption opacity due to gas and dust individually at each location in the disk from the Bethell & Bergin (2011a) cross sections. The X-ray scattering is dominated by Thomson scattering and approximately isotropic (i.e., photons are scattered equally in the backwards and forwards directions). For the primary model, we assume a characteristic T Tauri star X-ray luminosity of = 10 erg s with the quiescent X-ray spectral template given in Paper I.

Figure 2.— Sources of H ionization present in the midplane as a function of disk radius. The orange points are the result of our Monte Carlo calculation and are dominated by intermediate energy X-ray photons, typically  keV. The solid blue-green line is the cosmic ray ionization rate including the modulation of stellar winds for the present-day Sun at solar maximum. The CR value typically assumed for disk chemical models is shown in dashed blue-green for comparison. The yellow line shows the initial contribution from short-lived radionuclide ionization, predominantly Al. The radial decrease in SLR ionization is due to losses of SLR decay products, which escape the disk. The effective half-life of the SLR curve is  Myr.

The most abundant low energy X-rays ( keV) are also the most easily attenuated, and consequently do not contribute to the dense midplane ionization. The most important X-rays in the densest gas are typically of intermediate energies, between  keV, which are less numerous but can more readily penetrate the gas and dust. We emphasize that a detailed treatment of X-ray scattering is crucial in understanding disk ionization, otherwise there would be no X-ray photons in the midplane (Igea & Glassgold, 1999; Ercolano & Glassgold, 2013, Paper I). Correspondingly, the scattered X-ray radiation field provides the absolute floor to the midplane ionization rate in the absence of CRs and SLRs (see Figure 2). In the Figure, the initial rise in midplane X-ray intensity occurs as the disk surface density drops (becomes more optically thin), while the fall-off beyond  AU is simply geometrical dilution. The X-ray ionization floor falls in the range  s. For a more massive or denser disk than the one considered here, the role of X-rays can be diminished in the midplane by more extreme gas extinction. On the other hand, a star with higher X-ray luminosity would have a proportionally higher depending upon the shape of the stellar X-ray spectrum (see also §6.1).

3.2. Short-Lived Radionuclide Ionization

The disk ionization contribution from the decay of SLRs has been studied extensively (Consolmagno & Jokipii, 1978; Umebayashi & Nakano, 1981; Finocchi & Gail, 1997; Umebayashi & Nakano, 2009; Umebayashi et al., 2013; Cleeves et al., 2013b). The initial abundances of radioactive species within a typical protoplanetary disk are uncertain but can be estimated for the case of the protosolar disk from isotopic abundance measurements in meteorites (e.g., Wasserburg et al., 2006, and references therein). Whether these values are representative of all disks is unknown; however, the frequency of differentiated extrasolar asteroids seem to indicate that the young Solar System was at least not atypical in its SLR content (Jura et al., 2013), a hypothesis which models successfully reproduce (Vasileiadis et al., 2013). Moreover, recent work indicates that both direct (disk) SLR injection in clusters and distributed SLR enrichment in molecular clouds can produce abundances comparable to those inferred for the early Solar Nebula (Adams et al., 2014).

Time evolution of abundances adds further uncertainty in estimating the SLR contribution to the ionization budget in disks. The characteristic half-life of the ensemble of important SLRs, primarily Al ( Myr; Schramm, 1971; MacPherson et al., 1995; Umebayashi & Nakano, 2009) and Fe ( Myr111We note that in Cleeves et al. (2013b), the value assumed for the Fe half-life predates Rugel et al. (2009), where the original value of  Myr used in Paper I originates from Kutschera et al. (1984). Once the Al has been depleted after approximately  Myr, Fe is the primary contributor to SLR ionization. However at this epoch, the ionization attributed to scattered stellar X-rays begins to exceed that of the SLR contribution (see Fig 2), and correspondingly the change in chemical properties or ionization fractions are not expected to be large for the updated half-life. Nonetheless, we have included an updated Figure 4 from Cleeves et al. (2013b) and new fits to the midplane ionization rate in Appendix A of this manuscript.; Rugel et al., 2009) corresponds to approximately  Myr (Appendix A). This implies that for disks aged 5 Myr, the contribution from SLRs is reduced by nearly 80%, and consequently scattered X-rays and SLRs would contribute to the midplane ionization at a similar magnitude outside of  AU (Fig. 2). Inside this radius, the SLR contribution exceeds that of X-rays due to the high attenuation of X-ray photons in the inner disk midplane.

Another complication is that a substantial fraction of the SLR decay products, e.g., beta particles and -rays, escape the disk prior to ionizing the gas when surface densities drop below  g cm so that not all of the available energy is deposited locally. The present model incorporates the effects of -ray, and loss in the outer disk (described in Cleeves et al., 2013b) for the settled disk slab model. We note that the main chemical model results presented in §4.2 assume the SLR ionization rate is constant with time, which is acceptable for disks with  Myr or less. We relax this requirement in §6.2 where we consider SLR time evolution within the chemical calculation itself.

3.3. Cosmic Ray Ionization

Ionization by cosmic rays at the interstellar rate ( s) is a commonly assumed ingredient in models of disk chemistry and physics. However, detection limits of HD emission suggest that the ionization rate is actually lower than expected for the dense ISM, pointing to some variety of additional attenuation. One possible explanation is CR exclusion by natal stellar winds, i.e, an analogue to the modern day “Heliosphere.” Within the Heliosphere, the solar wind strongly modulates CRs, especially those below  MeV. Pre-main-sequence stars such as T Tauri stars are significantly more magnetically active, have high rotation rates, and high mass loss rates, all of which may drive high levels of CR exclusion at all particle energies. Models of the Gyr-old Sun show significant modulation even at late times (Svensmark, 2006; Cohen et al., 2012), corresponding to incident CR ionization rates in the range  s. In Paper I, we present a simple model of a scaled-up Heliosphere, i.e. a “T-Tauriosphere,” to estimate the degree to which a T Tauri star could potentially exclude CRs with winds, and how this exclusion is imprinted on the ionization state of the disk. In the instance of relatively mild modern-day solar winds at solar minimum, the ionization rate by CRs at the disk surface is reduced to just  s. This value is already over an order of magnitude below the typically assumed (dense gas) interstellar rate. At solar maximum, the ionization rate at the disk surface is reduced to  s (see Figure 2). If young, energetic stellar winds are even more efficient at sweeping away CRs, the CR ionization rate will drop below the solar maximum value at which point SLRs now become the primary midplane ionization contributor, where  s for  Myr.

While the “T-Tauriosphere” models explored in Paper I were simple scaled up versions of the modern day solar wind, the spectra reflect the general behavior we expect under more extreme stellar wind conditions such as those explored in Cohen et al. (2012), where low energy CRs are excluded most severely, and high energy CRs ( GeV) are only weakly modulated. Here we explore two wind-modulated CR ionization models, including the present-day Sun (Solar System Maximum: SSX) and a scaled-up exclusion model (T Tauri Maximum: TTX) as presented in Paper I. These models provide a realistic framework that allows us to quantify how the chemistry is affected by modulated incident CR fluxes. We note that the wind modulation affects the incident spectrum and consequently the ionization rate at the disk surface, in addition to attenuation by the disk gas, where the typical attenuating surface density is  g cm (Umebayashi & Nakano, 1981). Additional opacity can arise from magnetic mirroring and magnetic irregularities in the disk, which we discuss further in §7 (see also Dolginov & Stepinski, 1994; Padovani & Galli, 2011; Fatuzzo & Adams, 2014).

For the scaled up (TTX) CR model discussed in Paper I, the ionization rate due to CR drops below that of X-rays at all radii, and the TTX models correspond to a purely X-ray dominated disk for models without SLR ionization. Table 2 provides characteristic rates in the disk surface layers for the different CR wind modulation models considered in this work. In addition to the wind modulated SSX and TTX spectra, we examine two “ISM-like” models, M02 (Moskalenko et al., 2002) and W98 (Webber, 1998). W98 is the closest model to what is typically assumed in models of disk chemistry,  s, while the M02 ionization rate is characteristic of CR ionization rates measured from H in the diffuse ISM (Indriolo & McCall, 2012), which are unattenuated. For disks without a significant amount of surrounding nebulous gas, such as TW Hya, the M02 model may be more representative of the interstellar (completely unmodulated) CR rate.

Model ID
Moskalenko et al. (2002) M02 s
Webber (1998) W98 s
Solar System Maximum SSX s
T Tauri Maximum TTX s
Table 2CR model ionization rates for cm.

4. The Chemical Network

The chemical nature of a parcel of interstellar gas and ice is highly sensitive to the properties of its environment. These properties include the radiation field intensity, grain size and volume density (regulating freeze-out and desorption), gas/dust temperatures and densities. In cold gas ( K), the gas-phase chemistry is particularly sensitive to the gas ionization fraction, as ion-neutral reactions are the most important (i.e., quickest) gas phase reactions that take place (Herbst & Klemperer, 1973). The geometry of disks, ranging between lower density, warm ionized surfaces to predominantly neutral, cold midplanes provides a diversity of physical conditions, which are directly translated into a rich chemical environment. Conversely, observations of the chemical composition of disks provide clues into the underlying physical environment and therefore are a powerful observational tool to help understand the important physics governing these systems. In the following section we examine how the chemical properties of disks, particularly the molecular ions, react to different assumptions regarding ionization processes, and how they may be used as diagnostic tools.

4.1. Chemical Model

The physical model and stellar UV and X-ray radiation fields provide the backbone on which we solve for the time-dependent chemical abundances with the Fogel et al. (2011) disk chemistry code. This chemical reaction network is based upon the OSU gas-phase network (Smith et al., 2004), where Fogel et al. (2011) substantially expanded the network to include important processes such as thermal and non-thermal sublimation, photodissociation, freeze-out onto grains, CO and H (and isotopologues) self-shielding, and stellar and non-stellar ionization of H and helium. The code is calculated as 1+1D, where different disk radii are treated independently and self-shielding is considered in the vertical direction. The code is run in parallel with the publicly available GNU Parallel software (Tange, 2011). We note that in the calculation of the temperature structure and UV transfer, we consider the spatial dependence of the grain size populations in detail. The grain-surface chemistry, however, is fixed to a single “typical” grain-size of m, where the underlying assumption is that the small grains dominate the surface area and are most important for the chemistry. Because the gas and small grains are uniformly distributed, this approximation is justified in the present work; however, we would nonetheless be “missing” surface area in the midplane where the 1mm grains are concentrated. We note, though, that this correction should be small considering that small grains are present throughout the disk at all scale-heights. We furthermore note that the physical size of is not the important quantity, but rather this number translates into an effective grain surface area per unit volume in the chemical code (e.g.,   cm at  cm). Nonetheless, future work should explore more than two size populations, including their vertical distribution in the disk (for example, in the formalism of Dullemond & Dominik, 2004) and how these will affect the grain-chemistry.

In order to make predictions for deuterated molecular ions, we have added to the Fogel et al. (2011) network a simple deuterium chemistry to predict the abundances of HD and ND, which become enriched relative to the main isotopologue due to a chemical favorability towards the heavier isotopologue at low ( K) temperatures (Millar et al., 1989). Even though DCO is part of the network through the HD+ CO formation pathway, we do not make predictions for the DCO abundances in the present work; the chemistry of DCO depends sensitively on the deuterated carbon chemistry, for which we have not included a complete network. Instead, the network is designed to reliably predict the relatively simpler HD and closely related ND chemistry. An important facet of the HD chemistry is that the reaction rates depend strongly on the ortho-to-para ratio of the reactants, H and H, and their isotopologues (Pagani et al., 2009; Hugo et al., 2009; Chapillon et al., 2011). We approximate this behavior by assuming that the ortho-to-para ratio is locally thermal for HD and include the ortho-to-para fraction as weights on the net reaction rate. In the following prescription we designate () as the fraction of H in the ortho (para) state, and () as the fraction of HD in the ortho (para) state. The weighted reaction rate coefficients for (o)rtho and (p)ara H and HD is:


where , , and is the reaction rate coefficient for reactions between o-H and o-HD, for example. The o/p ratio of H is given by


for gas at temperatures above  K and for temperatures below that value,


Flower et al. (2006) finds that at very low temperatures ( K) the ortho-to-para ratio of H exceeds the Boltzmann value, and to approximate this behavior we set a floor to the o/p ratio of H at 10 that limits Eq. (3) from dropping below this value. The o/p ratio of HD is given by


where both o/p formalisms in Eqs. (4.1) and (4.1) are from Lee & Bergin (2014), submitted.

In the chemical network, the deuterium extension to the reaction set – for the most part – mirrors the main isotopologues, and we assume statistical branching ratios where necessary. While the ionization of H and helium are the first step towards the molecular ion chemistry, we do allow for HD and D to be directly ionized by X-rays, CRs and SLRs, assuming the same cross sections as for H. The self-shielding functions for HD and D are from Wolcott-Green & Haiman (2011). For the reactions where rates are known for the deuterated isotopologue, we incorporate values from the literature for HD electron recombination and reactions with atomic H (Roberts et al., 2004), H and D reactions with atomic/molecular hydrogen/deuterium, DCO, ND reactions with H (Roberts & Millar, 2000), and neutral-neutral warm water-formation deuterium reactions including their significant barriers (Bergin et al., 1999). While we include the spin information for the H + HD reaction rates, the endothermicity of the back reactions for the overall less abundant heavier isotopologues, HD and D are taken to be a single value in the network, 187 K and 341 K, respectively (Roberts et al., 2004). In addition to the standard set of gas-phase reactions in the original network, we include a simple grain-surface chemistry allowing for H/HD/D, HO/HDO, and HCO/HDCO formation on grain surfaces through hydrogenation (Hasegawa et al., 1992) and CO formation through reactions of CO with OH and O with barriers (Garrod & Pauly, 2011). In total, the full network spans over reactions and species.

The initial chemical abundances are listed in Table 3 and are for the most part adopted from Fogel et al. (2011) based upon the dark cloud model of Aikawa & Herbst (1999b). The abundances of CS and SO have been adjusted to the observationally derived abundances of the TMC-1 dark cloud from Bachiller et al. (1990) originally compiled in Guélin (1988) and Rydbeck et al. (1980). Together, CS and SO constitute the main sulfur reservoir with an abundance relative to hydrogen number totaling . The atomic (ionized) sulfur, S, is reduced to a low value, where we assume . The motivation behind the total reduced sulfur abundance, rather than the diffuse ISM value of , comes from the observation that the volatile sulfur abundance in dense clouds is far lower than that of the diffuse ISM (Joseph et al., 1986; Tieftrunk et al., 1994). These findings were confirmed by observations of molecular sulfur-bearing species including CS and SO (e.g., Langer et al., 1996; Wakelam et al., 2004) where the molecular abundances measured in dense clouds are far less than those predicted by chemical calculations for a diffuse ISM sulfur abundance (Hatchell et al., 1998; Wakelam et al., 2004). In particular, Wakelam et al. (2011) find that based upon chemical modeling of dense (high mass) cores, the volatile sulfur abundance is found to be even further substantially reduced, where the observations account for an abundance of sulfur totaling relative to hydrogen across four sources. The interpretation is that the sulfur has been converted into a more refractory form, i.e., a “sulfur-rich residuum” (Wakelam et al., 2004), and is not chemically available to produce volatile sulfur-bearing molecules. Evidence of high levels of atomic sulfur in shocked gas near Class 0 protostars further supports this scenario, where in one explanation the observed sulfur atoms are sputtered from the sulfur-rich residuum (Anderson et al., 2013). Likewise, the ionized metal abundances of Si, Mg, and Fe have been reduced to low levels () based upon the chemical modeling results of Maret & Bergin (2007) for the B68 pre-stellar core. The HD abundance is based upon the protosolar bulk value, where D/H in HD is (Geiss & Gloeckler, 2003). The HD, HD, and D abundances are estimated from Walmsley et al. (2004) for the small (0.025 m) grain case at high ( cm) density. Regarding the nitrogen abundances, we note that in these models, the initial abundance of nitrogen is primary atomic. Had we begun with a larger fraction molecular nitrogen, the overall NH abundance would increase but the shape of the column density profile remains largely unchanged (Schwarz & Bergin 2014, submitted; Cleeves et al. 2014, in prep.). Starting off with a larger fraction of ammonia, for example, tends to reduce the overall NH abundance across the bulk disk (outside of the ammonia snowline) by removing nitrogen from the gas-phase. Consequently, the shape of the column density profile of NH is still sensitive to the CR ionization rate even without additional information on the disk nitrogen abundances.

Species log() Species log()
H HO(gr)
He N
Si S
Mg Fe
Table 3Initial chemical abundances, .

4.2. Chemical Abundance Results

Figure 3.— Disk chemical abundances relative to H as a function of radius and height, where = (0,0) AU is the location of the central star. Abundances are shown at time  Myr of chemical evolution. Columns correspond to different CR/SLR ionization models for a fixed X-ray luminosity, and rows are different molecular abundances as labeled on the leftmost column. Specific CR and/or SLR ionization model is indicated at the top of each column. The two rightmost columns include contribution from the decay of short-lived radionuclides adding at most  s to the total H ionization (see Figure 2). Their inclusion provides a floor for the molecular ion abundance (see especially NH).

In this model framework, we compute the time-dependent chemistry as a function of radial and vertical position within the disk. The molecular abundances for a select set of ALMA observable molecular ions and of CO after  Myr of chemical evolution are shown in Figure 3. The column headings indicate the underlying CR ionization model as described in §3.3 (see also Table 2) with the X-ray luminosity fixed ( erg s). In addition to the abundances as a function of spatial position, we provide the vertically integrated column density, shown in Figure 4.

In general, for decreasing CR ionization, the midplane ion abundances drop precipitously while the X-ray dominated surface layers are unchanged across CR models. For example, the important dense ionization tracers, HD and ND, are highly abundant ( relative to H) in the midplane for CR ionization rates exceeding  s, but drop in abundance by more than three orders of magnitude for SSX modulated rates and below.

It is important to point out that the TTX model provides a unique chemical picture of a purely X-ray dominated disk. In the instance of a modulated CR rate and either (i) a lack of abundant SLRs or (ii) an old ( Myr) disk in which the initial reservoir of SLRs has decayed away ( remaining), the scattered intermediate energy ( keV) X-ray photons set the absolute floor to the ionization rate, with magnitude  s for a stellar X-ray luminosity of  erg s. In the TTX model, the midplane abundances of molecular ions are at their minimum (see especially HD, HD, and NDin Figure 4). Correspondingly, HD and HD have column densities typically three orders of magnitude less than if CRs are present at ISM levels and will unfortunately be difficult to detect observationally.

In the rightmost two columns of Figure 3, we include (static) contribution from SLR ionization, which provides an ionization rate floor of magnitude  s for fixed initial solar nebula-like SLR abundances (Finocchi & Gail, 1997) of Al, Cl, and Fe. The vertical and radial profile of the SLR ionization rate is taken from the settled disk model of Cleeves et al. (2013b) (hybrid dust slab, see their Figure 2 as well as Figure 2 in the present work). For the case of SSX+SLR, the two sources contribute similarly to the total H and helium ionization, so chemically there is little change with the inclusion of SLRs, especially in the inner disk. The TTX+SLR model, however, increases the total ionization such that the molecular ion abundances resemble the SSX runs inside of  AU. Consequently, distinguishing between a SLR and CR driven chemistry will be challenging. Outside of  AU, however, the SLR losses become important and the TTX+SLR model drops off more steeply than the SSX model, such as can be seen in HD.

Figure 4.— Vertically integrated column densities (cm) of the indicated species as dependent upon the magnitude of non-stellar ionization sources. Line color indicates CR/SLR ionization model with the approximate  s value for the incident CR ionization rates (plus SLR, if applicable), where the X-ray luminosity is fixed as  erg s. HCO and NH reach a “floor” in their column density from a stellar X-ray ionization baseline. X-rays contribute less to ND, HD, and HD and so their column density is more sensitive to both the high and low ( s) ionization models. CO column densities are provided in the lower right, and at high cosmic ray ionization rates the CO abundance is eroded by reactions with He.

Focusing now on the individual ions, the column density of HCO changes by an order of magnitude between the high ( s) and low CR models. The low CR models pile up at a constant HCO column, which physically corresponds to the minimum amount of HCO provided by X-ray ionization in this model. In the absence of interstellar CR, the bulk of the HCO column is therefore sensitive to only the stellar X-ray ionizing radiation (see §6.1). A second layer of HCO forms at , above the vertical CO freeze-out region, when CRs are present (models W98 and M02). This layer all but disappears for the SSX and TTX models, pointing to the potential utility of disk vertical structure observations in understanding the underlying ionization environment. It is important to note that even though we do not vary the stellar UV field in this work, for very high UV fields the CO that would otherwise form HCO can be photo-dissociated and photo-ionized, thus forming C in abundance. Thus the HCO will be indirectly UV sensitive, particularly for UV luminous Herbig Ae/Be stars (Jonkheid et al., 2007), which can be orders of magnitude brighter than T Tauri stars, or in the case of extremely externally irradiated disks. Nonetheless, the general trends over different ionization models should still hold, albeit with a lower overall HCO column. Moreover, while the HCO column can provide constraints on the ionization rates exceeding  s, the HCO is not very sensitive to lower CR rates due to its precursor, CO, freezing out in the CR dominated region. Moreover, the HCO column densities for the W98 and M02 models have a second peak in column density at  AU not seen for the low CR models. The second peak is actually a consequence of a deficit of CO at  AU, seen in the same plot. Because HCO forms directly from CO, it is sensitive to the CO chemistry, both freeze-out and CO chemical processing as is the case for the W98 and M02 models.

The NH column shows a three order of magnitude spread between a high CR rate (M02) and an X-ray only disk (TTX), with one order of magnitude variation between the different low CR models (SSX and TTX) and is a potential candidate as a midplane ionization tracer. One caveat in using NH is that a tenuous () surface layer of NH is sustained from the high X-ray photon flux (see Figure 3) even in the presence of CO, which is the major reactant of NH. For more X-ray luminous stars, the surface NH may contribute more to total column and mask the midplane abundances. Alternatively, the contribution from the surface NH can be reduced if the N binding energy is higher than 855 K, i.e., the N value when CO and N ice are well-mixed where the CO binding energy is the same, 855 K (Öberg et al., 2005). Such is the case for an N ice layered on an HO ice substrate (Collings et al., 2004), but this surface NH is never completely absent in any of our models.

The deuterated molecular ions considered here, HD, ND  and HD, show a promising sensitivity to the low ionization models, spanning many orders of magnitude in their column densities. Observations of these species will be essential for measuring midplane ionization, especially in the absence of CRs. This result is a natural consequence of the pathways towards deuterated isotopologues being favored at low temperatures, resulting in their overabundance relative to the main isotopologue in cold ( K) gas. Naturally, the same gas is also the least affected by X-ray ionization, and therefore these species allow us to peer through the X-ray dominated upper layers directly to the midplane.

In addition to the ions, we plot the CO abundance and column density (Figures 3 and 4). CO is a commonly used tracer of gas mass assuming an ISM conversion factor of CO to hydrogen mass of 10 and is frequently used to determine chemical abundances per H from an observed column of optically thin gas. There is evidence, however, that the CO abundance may be lower, and hence gas masses from CO may be underestimated (Favre et al., 2013). A possible reason for the low observed abundance is that CO is chemically eroded over time by reactions with He. The majority of the carbon returns back to CO eventually, but some non-zero fraction of the carbon is recycled into other carbon-bearing molecules, reducing the CO abundance over time (Bergin et al., 2014; Furuya & Aikawa, 2014). For the highly ionized CR models W98 and M02, the abundance of CO is visibly reduced near the inner disk midplane at around  AU due to an increased abundance of He and consequently a speed-up of the CO erosion process. The column density plot (Figure 4) also reflects this behavior in the M02 and W98 models and is an example of how ions can have a long-lasting effect even on abundant neutral species. We emphasize that the layered structure of CO induced by the ionization-chemistry in the abundance plot (Fig. 3) may be smeared out in the presence of turbulent motions of the gas (Semenov & Wiebe, 2011; Furuya & Aikawa, 2014), which is beyond the scope of the current paper but should be explored in future work.

5. Line Emission Modeling

The chemical models demonstrate a sensitive link between abundances (and column densities) and ionization properties of the disk. For example, the HCO and NH column densities typically are sensitive to stellar X-rays, though not exclusively for high interstellar CR ionization rates. The deuterated ions trace deeper gas, and probe ionization properties near the midplane, tracing ionization due to CRs and SLRs. In this section we show how these effects are imprinted on the molecular emission and how excitation effects and opacity can “mask” the molecular column densities that allow us to discriminate between models. From these emission models, we simulate realistic ALMA observations to determine the utility of emission line tracers as probes of individual ionizing agents.

5.1. Line Radiative Transfer

For our molecular ion abundance results (§4.1), we have simulated observations of the ground-accessible submillimeter transitions of each species considered here. The strength of the line emission depends on both the total column of material as well as the temperature of the emissive gas within the column. To simplify the problem, we have simulated the emergent line intensity assuming the disk is observed face-on at a distance of  pc. The line radiation transfer is carried out with LIME (Brinch & Hogerheijde, 2010) where we create simulated “perfect” images of the line emission in 100 m s channels in pixels of size 0.02 (2 AU), which is much smaller than our desired resolution in the ALMA simulations. For gas motions, we assume the gas is in Keplerian rotation about the star, and we include a turbulent doppler B-parameter of  m s, as observations indicate that the turbulent broadening in disks is small (Hughes et al., 2011).

For the HCO and NH emission models, we calculate the level populations in non-LTE with the collision rates of Flower (1999) as compiled in the Leiden LAMDA database (Schöier et al., 2005). For the latter we do not separate out the hyperfine structure lines, and hence refer to the transitions by their rotational-J states. Given that the collisional rates for ND are unknown, we treat the ND level populations as in local thermodynamic equilibrium (LTE), which is a satisfactory approximation as the critical density of the () and () transitions are  cm and  cm, respectively, as estimated from the Flower (1999) collision rates for HCO and line parameters (Einstein A-coefficients, frequencies, statistical weights) from the JPL Database222http://spec.jpl.nasa.gov/ (Pickett et al., 1998), originally measured in Anderson et al. (1977); Sastry et al. (1981). We note that Hugo et al. (2009) provide inelastic collisional rates for the excitation of HD; however, to simplify the calculations we assume the HD level populations are also in LTE. This assumption is appropriate as the HD emission originates from dense gas exceeding n cm, much larger than the critical density for o-HD (),  cm (Hugo et al., 2009).

For the carbon and oxygen isotopes, we adopt (Kahane et al., 1992; Prantzos et al., 1996) and (Keene et al., 1998). The deuterium isotopologues of the molecular ions are calculated within the code where the main gas reservoir (molecular hydrogen) has initially (see §4.1). Leiden LAMDA formatted HD and ND input files were compiled from the CDMS database333http://www.astro.uni-koeln.de/cdms/catalog (Müller et al., 2001, 2005) and from the JPL database, where the primary literature regarding the line parameters can be found in Saito et al. (1985); Amano & Hirao (2005); Yonezu et al. (2009) for HD. All simulations are carried out using the same dust distribution and opacities from the physical structure model, where the continuum is subsequently subtracted from the resulting line emission profiles.

5.1.1 Line Opacity

The LIME code is capable of providing physical line intensities, e.g., Jansky/pixel and Kelvin, as well as the line optical depth, . In Figure 5, we show the line-center, vertical optical depth through the disk of the simulated transitions (direct from the emission model, i.e., no beam convolution) as a function of cylindrical disk radius. The HCO isotopologues are for the most part thin (), while NH () and () reach just beyond the CO snowline for the W98 and M02 CR models. Observations of HCO (3-2) and HCO (3-2) are reported in Öberg et al. (2011b) towards DM Tau. The disk integrated flux ratio of the (3-2) transition for the C/C isotopologues is , and assuming a respective isotope ratio of 60 (Keene et al., 1998), corresponds to an HCO (3-2) optical depth of . Though this value is optically thin, it is nonetheless higher than the values from our models. Potential explanations include a potentially higher X-ray luminosity from the star, which is unknown, or that the 19 AU inner cavity (Calvet et al., 2005; Andrews et al., 2011) may permit the X-rays to more directly ionize the outer disk gas. The ND and HD lines are thin throughout the disk for all models. CO is thick, with at all radii, though the CO opacity drops slightly near the star where high CR+X-ray ionization chemically destroys CO with He more quickly than it is replenished. We emphasize that the molecular ion opacities may be higher for higher ionization rates or if the column of gas is larger for a more massive disk, and vice versa for lower mass models (see discussion in §6.3).

The HCO isotopologues show a peak in their opacity offset from the radial center of the disk, where the column is highest (see Fig. 4). The relative height between the central peak and the outer broad  AU peak is due to a combination of column density and excitation. The ring-like feature is exaggerated for low rotational J lines such as the () transition where depopulation becomes important towards hotter inner disk gas, whereas HCO () appears more centrally peaked.

The NH lines have two main emission features, where the inner ring is caused by a contribution of stellar X-rays and CR ionization in tandem with the CO condensation front, while the outer tail at AU is fueled by primarily CR ionization. The HD optical depth also peaks in this same outer region in the models where CR ionization is active. The ND has a somewhat complicated profile, which is reflected in the ND column densities (Fig 4). The high CR ionization models sustain ND co-spatial with NH and, when combined with warm temperatures close to the upper-state excitation temperatures ( K and 37.0 K for and rotational levels, respectively), leads to ND emission inside of  AU.

Figure 5.— Face-on line center optical depth () of the labeled emission lines as a function of radial distance from the star. The y-axis values are to be multiplied by the value in the lower left corner of each plot. The emission lines for the most part are optically thin except in the case of NH () and (), which reaches of a few inside of  AU. CO, a commonly used tracer of gas-mass, is optically thick for most of the disk, and therefore either a more optically thin isotopologue is needed or sensitive observations of the line wings where the gas may become thin to use CO as a potential mass normalizer. Colors correspond to M02: purple, W98: yellow, SSX: blue, and TTX: green.

5.2. Simulated ALMA Observations

From the LIME emission models described in §5.1, we simulate full-science ALMA observations for the ground-based accessible transitions of HCO, HCO, NH, ND, and HD. The full set of emission lines considered at present are listed in Table 4. We note there are additional lines that are accessible by ALMA that we chose not to include, such as ND (5-4) and HD (), because these lines were unobservable for all emission models considered. The full-science array is comprised of fifty 12-m antennas, twelve 7-m antennas, and four 12-m antennas, i.e., the total power (TP) array, where the 7-m and TP arrays are referred to as the Atacama Compact Array (ACA). For each of the line models we assume the same set of simulated observing parameters, and note that in some instances there are more efficient (i.e., less ALMA time) means to achieve the same goals. The specific settings for actual observations should be tailored to the specific target and specific line being studied.

Emission 12-m Configuration444 Sensitivity555
Line (12-m + 7-m Beam) (mJy)
HCO () 03 ()
HCO () 03 ()
HCO () 10 ()
HCO () 10 ()
HCO () 10 ()
NH () 10 ()
NH () 10 ()
ND () 01 ()
ND () 01 ()
o-HD 01 ()
aafootnotetext: http://casaguides.nrao.edu/index.php?titleAntenna_List Beam averaged over major and minor axes from the combined 12-m and ACA observations.

. bbfootnotetext: 6.1 in 0.2 km/s bins from the ALMA Sensitivity Calculator.

Table 4CASA Simulation Parameters.

5.2.1 Observational Parameters

The simulations presented here reflect 6 of total on-sky time with one of the 12-m array configurations and an ACA 7-m observation. While the choice of 6 is somewhat arbitrary, it represents a relatively “deep” observation, and ideally one would target more than one emission line towards a given source in the same proposal. Thus this length of time may be somewhat optimistic, but is designed to detect the lines with high signal-to-noise. Some of the lines are undetectable for low ionization models as is the case for HD, and thus deeper observations would be needed to detect the line if that specific model reflects the true properties of the disk. In that case, going to a different tracer such as ND may be a better choice than a deeper observation of a weak line. Furthermore, these results can be approximately scaled for sources at different distances and the estimated signal-to-noise can be adjusted for different observation times.

Regarding antenna configurations, for stronger emission lines, such as HCO, and NH, we choose somewhat more extended 12-m configurations (see Table 4), while for weaker lines, we simulate the most extended configuration that still provides high signal-to-noise if possible, which in some instances is the most compact, least-resolved full-operations configuration, denoted 01. We include ACA 7-m observations because the maximum scale of the  AU disk at  pc on the sky is 8, and for the higher frequency ( GHz) lines considered here, the 7-m array recovers the flux at all scales without the need for the TP array. At lower frequencies, a second 12-m configuration can be used instead to recover all of the flux in significantly less time than for the 7-m observations, which is clearly the better – and default – option for on-sky ALMA observations. Additionally, the long duration of the simulated observations fills out the UV coverage and, due to sky-projection effects, can decrease the minimum baseline, allowing for slightly larger maximum recoverable scales.

In summary, the simulations presented here are designed to show the sensitivity of the emission line observations to chemical signatures of ionization processes in disks. These results provide a primer to aid in quantifying the contribution of different physical processes to disk ionization, a crucial ingredient in models of disk chemistry. We note that the specific choice of observing parameters will depend upon properties of the source, the ALMA cycle, which determines the available antenna configurations, and the lines targeted. Figure 6 illustrates the full process beginning with LIME emission models to creating ALMA on-sky simulations for two example emission lines.

Figure 6.— Schematic illustrating the calculation of simulated observations. Top: The face-on line emission models calculated from the chemical abundance models (Fig. 3) using the radiative transfer code LIME (Brinch & Hogerheijde, 2010) in non-LTE where available. Line emission models are in units of Jy/pixel and are observed at a distance of  pc. Middle: Using CASA’s simobserve task, we simulate observations with the full 12m ALMA array plus the 7m ACA array. UV-binned visibility amplitudes integrated over velocity and on-sky emission shown,  AU. Bottom: Cuts across the on-sky emission from the cleaned simulated observations, see text for details. Error bars (sensitivity) are estimated from the ALMA sensitivity calculator values assuming a 6.1 observation in 100 km/s channels. Average beam size is and for the left and right plots, respectively.

Rather than showing images or channel maps for every emission line/model considered both in model emission (§5.1) and CASA simulations, we instead present our results as circularly averaged visibility amplitudes and circularly averaged integrated line intensity cuts in sky coordinates, which facilitates easy comparison between observations of different ionization models for a given emission line.

5.2.2 ALMA Sensitivity Calculation

To calculate the uncertainty on the ALMA measurements, we use the uncertainties from the ALMA sensitivity calculator666https://almascience.nrao.edu/proposing/sensitivity-calculator, which set the width of the shaded regions in the image plots (rightmost column of Figures 7 and 8). We do not simulate thermal noise within the CASA simulations themselves to save computational time. To estimate the error on individual baselines at each time sampling (defined by integration time, which was s for the simulations), we scale the sensitivity calculator’s uncertainty by , where is the number of antennas (50 for the 12-m array and 12 for the 7-m array), is the number of time samplings during the course of the simulated observation’s integration time (e.g., on the 12-m array such that =22) and is the number of polarizations, where the emission is unpolarized. We note that if we had chosen a smaller (larger) sampling time within a factor of a few, we would have correspondingly more (less) points to bin over, and the final error on the binned measurement is not very sensitive to the choice of sampling time. From the full measurement set, we bin the velocity-integrated visibility amplitudes by projected UV distance. The uncertainty on the binned value (the mean) is taken to be the error on the individual measurements divided by the number of points within the UV distance bin, , where is typically .

5.2.3 Emission Model Results

We present the primary results of the ALMA simulated observations in Figures 7 and 8. The CR ionization models described in §3.3 are indicated by line color, where purple represents the highest, diffuse ISM CR rate, while the blue curve is representative of an X-ray dominated chemistry.

Figure 7.— ALMA simulations for the observable transitions of the indicated molecular ions. Line color represents different CR ionization models: M02 (purple), W98 (yellow), SSX (cyan), and TTX (blue). The simulated antenna configuration used in CASA’s simobserve task is indicated in parenthesis (## + 7m). Left: Velocity-integrated visibility amplitudes. The dashed vertical green line corresponds to the minimum baseline from the 12-m observations, which can be less than the physical separation of the antennas due to sky projection. The red dashed line is the maximum baseline sampled by the 7-m observations. Right: reconstructed image profiles in sky-coordinates integrated over velocity. Synthesized beam indicated by the black bar in the bottom left corner of the image plots. The width of the individual lines corresponds to the sensitivity of the simulated observations, see §5.2.2 for details.

The identifier of the full-operations array used to simulate each set of data is indicated in parenthesis next to the name of the species (see Table 4 for equivalent angular resolution). The left and right columns are different representations of the same simulated data. The left column shows the visibility amplitudes integrated over the line (Jy km s) where the quantity “integrated amplitude” corresponds to the fact that we have integrated over the channels that contain line flux rather than averaging, which is often done for continuum amplitudes. The visibility amplitude plots are labeled with a dashed green vertical line indicating the maximum recoverable scale for the 12-m array observations alone. The dashed red vertical line indicates the minimum scale sampled by the 7-m array. We note that the minimum baseline indicated by the green line may be less than what would be predicted for the same array based upon the physical baselines due to sky-projection effects mentioned previously. The right column shows the brightness profile of the reconstructed (cleaned) images on the sky, where is relative to the position of the central star and  AU. To create the profiles, we average ten radial slices across the face of the disk. The noise reflects the sensitivity on a single cut, but in practice averaging over slices can further reduce the noise on the profile.

Figure 8.— The same as Fig. 7, for “cold gas” ions typically tracing gas below  K. HD is unobservable for our models with ionization rates of  s. Consequently, ND is the only molecular ion which allows for some differentiation between ionization rates below  s, as seen clearly in the visibilities, especially for the ND (4-3) transition. The high CR ionization rate models overlap for NH because the line becomes very optically thick. We note, that while weak, the SSX model is detectable in both the () and () transitions of NH at the 3 and 6 levels, respectively. See also Figure 9 for a zoom-in on the low ionization models.

The HCO isotopologues’ line emission (Figure 7) is sensitive to both the stellar X-rays and the high CR ionization models (M02 and W98 in purple and yellow respectively). This behavior is reflected in the column densities (Figure 4) and in the simulated observations, both in the visibilities and reconstructed images. The line intensity is indistinguishable for lower CR rates for models TTX and SSX, green and blue, where stellar X-rays – and not CRs – set the minimum HCO abundance in this particular model. The high M02 CR rates form a clear emission ring (i.e. the HCO emission peaks offset from the star) for the () transition, less prominently seen in the (). The reason that the () emission is brightest away from the column density peak is due to excitation of the line, which has an upper state energy of  K. The abundant HCO gas inside of  AU has temperatures typically exceeding 100 K. At this point, higher J-states of HCO and HCO become more populated, thereby de-populating the J = 3 levels and decreasing the column of emissive gas, behavior reflected in the line optical depths (Figure 5) and emission plots (Figure 7). For higher rotational transitions of HCO, () and (), there is an inner () peak in the simulated observations (Figure 7) due to a confluence of high X-ray rates and abundant gas-phase CO. This feature is also present in the underlying HCO abundances, but the lower spatial resolution simulations adopted for HCO (see Table 4) do not resolve the inner 50 AU. In summary, the HCO emission sensitively traces i) stellar X-ray processes and ii) high CR fluxes, but is limited in its capacity as a midplane tracer due to CO freeze-out in the outer disk. Furthermore, the low-J HCO behavior highlights a situation in which observations of emission rings are excitation effects, rather than chemical or physical structure, such as a planet or a snowline.

NH is expected to arise from cold gas where its destroyer, CO, is frozen out at temperatures below  K (Aikawa et al., 2001; Bergin et al., 2002; Jørgensen et al., 2004). Nevertheless, in the strongly irradiated X-ray layers we find that some NH is sustained even in the presence of CO. Consequently, its strength as a diagnostic tool of midplane ionization is somewhat decreased by potential confusion with surface NH emission. As shown in Figure 8, the high ionization and low ionization models are clearly discernible, while the SSX and TTX models are – in a relative sense – far more difficult to tell apart. We note NH may nonetheless have utility as a midplane ionization tracer, and that the SSX and TTX models are both observable and differ by a factor of in brightness (see the zoom-in in Figure 9), but detailed modeling may be required to estimate the potential X-ray contribution to the NH column density. An additional caveat in using NH as an ionization tracer is that its emission can be partially optically thick (if not completely, ; Fig. 5). Indeed, the high interstellar CR models, M02 and W98, have similar emission line strengths because the column densities are quite large, and the correspondingly thick emission lines are no longer sensitive to column density. Therefore, in this model, NH observations would help distinguish between disks where CRs are present or where they are excluded, e.g., by winds, however it is perhaps not a precise tool for determining CR rates below  s for this particular disk structure. In Figure 8, we note that even though the SSX model emission is substantially lower than that of the W98 and M02 models, it nonetheless is still detectable. The width of the radial profile corresponds to the ALMA sensitivity, and thus disks with SSX-level of CR ionization are detectable at the and level for the () and () transitions of NH, respectively. Given the NHabundance’s concurrent chemical dependence on temperature, a warmer disk where the NH abundance peaks further from the X-ray bright star may allow NH to be a more sensitive tracer of non-stellar ionizing processes and be less optically thick at the emission peak due to lower outer disk gas densities.

ND shows interesting emission behavior such that the M02 models show less ND emission than the more weakly CR-ionized W98 models (Figure 8). This behavior is reflected in the column densities (Fig. 4) and the abundances (Fig. 3). More specifically, the M02 ND models show a deficit in abundance relative to the W98 models centered at  AU. The inner ND gas is co-spatial with the high NH (fueled by X-ray flux) and with the boundary of the HD abundant region ( K), creating an ND layer. This layer is reflected in the W98 models as well. At intermediate distance radii  AU, the M02 models show an ND deficit as a direct consequence of the high degree of CR ionization in the presence of cold gas where freeze-out becomes important. The combination of ionization and cold gas drives important non-equilibrium chemistry. More specifically, the nitrogen is not being recycled back into NH but instead sequestered into nitrogen bearing ices like NH, HCN and NO, and therefore is not available as N in the gas to reform NH and ND. This process is an analogue to a similar sequence that may occur for CO (Bergin et al., 2014). In the M02 model, these pathways are specifically triggered by M02’s high CR flux, but a factor of brighter X-ray luminosity would have a similar effect (see §6.1).

The true utility of ND is highlighted at low CR ionization rates. The SSX models are clearly distinguishable from the TTX models, which is reflected in the emission plane for both () and () and in the column density plane. ND emission is always very optically thin, making it a direct tracer of column regardless of CR rate. The difference between the SSX and TTX line intensities is a factor of about , more than double that of NH, allowing these CR models to be more easily disentangled than with NH.

Figure 9.— ALMA simulations of the low CR ionization rate models with and without a constant SLR ionization rate. Only the low ionization models are shown for clarity. Figure quantities are the same as Figure 7. The CR-only models are shown as TTX: blue and SSX: cyan. The same CR models now including SLR ionization are shown in green and purple, respectively. As can be seen from both visibility curves and sky-emission, the inclusion of SLR ionization does not significantly change the profile or intensity of the lines shown. A TTX model with SLR ionization (purple) would be difficult to distinguish from an SSX model without SLRs (cyan).

Finally, HD is a commonly used cold ionization tracer. However, as can be seen in Figure 8, its emission is only detectable for interstellar CR rates or higher, owing to its weak line strength. For ionization rates at or below  s, the HD () is undetectable even for full ALMA operations, and therefore it is only a useful tracer of interstellar CR rates, if they are present. These results are consistent with existing limits on the observed HD column towards the TW Hya protoplanetary disk (Chapillon et al., 2011; Qi et al., 2008); however, such limits are much higher than all of the line strengths predicted here and thus more sensitive observations are required to determine if CRs are present with HD as a tracer. The utility of HD is highlighted when used in tandem with ND. As mentioned above, ND decreases in brightness for both high ionization rates (when the precursors of ND are chemically destroyed) and low ionization rates (when ND is not produced). Observations of HD in conjunction with ND would allow one to break the degeneracy between these scenarios.

In Figure 9, we show simulated ALMA observations for the models including SLRs (no time decay) for the SSX and TTX models. The NH and ND lines are the only tracers for which there may be a measurable difference for the low ionization models. From these plots it is apparent that distinguishing between CR fueled chemistry and SLR chemistry is extremely difficult. Additional factors such as – highly uncertain – disk ages will be necessary to determine fractional contributions. Alternatively, if an otherwise unexpected “jump” in ionization is seen, e.g., at the boundary of a T-Tauriosphere, then the contribution from each component can be determined unambiguously. Without this additional information, however, measurements of dense gas ionization using emission line tracers will most likely reflect a combination of both CR and SLR effects.

6. Further Considerations

In this section, we relax certain assumptions of our model and explore how our model results depend upon these additional parameters, including X-ray luminosity, temporal decay of SLRs, and the assumed mass of the disk.

6.1. Higher X-ray Luminosity

Figure 10.— Chemical models for an X-ray rate of  erg s (gray solid lines) and an enhanced rate of  erg s (black lines). Rows correspond to the indicated emission lines as labeled in the left column, and columns correspond to the CR models as labeled at the top of the figure. For the enhanced models, the three line styles correspond to different X-ray spectral templates. The solid black line holds the spectral shape fixed and increases the overall luminosity by a factor of 10. The dashed black line corresponds to a harder X-ray spectrum with the same luminosity as the solid line, normalized such that there are fewer photons with  keV and more with . The dotted line is a model with the same 1 keV flux as the solid black line, but with a hard spectral template such that the X-ray luminosity is also slightly higher, a factor of 1.7, due to the hard X-ray tail).

The results presented in Figures 79 consider a single X-ray luminosity,  erg s. To understand the sensitivity of the lines to X-ray ionization, we have computed an additional chemical model for a ten-fold increase in X-ray luminosity for three spectral shapes: (1) the same “quiescent” spectral shape for the baseline model, solid line; (2) a harder X-ray spectrum as was used in Cleeves et al. (2013a), with the same normalized luminosity as model (1), dashed line; and (3) the hard X-ray spectrum normalized to have the same X-ray flux at 1 keV as model (1), dotted line. The results of these higher X-ray luminosity models are shown in Figure 10, where the standard model is shown as the gray solid line and the elevated models are shown in black.

The high X-ray ionization rate changes the HCO column density most significantly in the inner  AU for the M02 and W98 models; for the outer disk there is a reduced effect as the contribution by CRs is more important. For the reduced CR models, the HCO column density is enhanced throughout the disk by a factor of 3-4. This trend is explained by the balance between ionization and recombination, such that the steady state abundance of ions (or electrons) is proportional to , where is the recombination rate and is the number density of H. With everything else constant, a factor of ten increase in ionization results in a factor of increase in the abundance and consequently column density. Consequently, owing to the sensitivity of HCO to the X-ray ionization field, observations of optically thin isotopologues of HCO may help put constraints on the permeability (optical depth) of the disk gas and dust to X-ray photoionization if the stellar X-ray luminosity is known. Furthermore, additional constraints on the disk gas mass combined with X-ray measurements would provide an approximate measure of the opacity due to dust and gas in the X-ray irradiated layers. Variations induced by the difference in spectral templates (black lines) are smaller, typically a factor of 1.2-1.8 in the TTX models. The hard X-ray spectrum considered in (2) (dashed line) even has a slightly lower molecular ion column density compared to the softer X-ray spectrum of (1) owing to the initial destruction of CO and N by hard X-ray generated He. Model (3), i.e., the spectrum normalized at 1 keV with model (1), has a slightly higher energy integrated luminosity, 1.7 times that of models (1) and (2), at which point production overtakes the destruction of precursors. The NH, ND, and HD column densities do not become sensitive to variations in the X-ray luminosity and spectrum until very low CR rates, primarily TTX, where most of the ionization in these cases comes from stellar X-rays. The flat increase in ionization for the TTX models is a consequence of high X-ray ionization rates near the star and slow recombination at large radii, where the drop in density (recombination) balances the decrease in X-ray flux with distance from the star.

6.2. Short-Lived Radionuclide Time Decay

Another simplifying assumption made in the chemical abundance calculations in Figures 3 and 4 was the use of a constant, non-decaying SLR ionization rate for the models including radionuclide ionization sources. From the results of Cleeves et al. (2013b) (see also Appendix A), the ensemble of short-lived radionuclides that dominate the ionization, namely Al and Fe, have an effective net half-life of approximately  Myr (see §3.2 for further details). For disk lifetimes up to 1 Myr, the change in ionization rate is correspondingly less than a factor of two, but can become significant for older disks ( Myr). To simulate the decay of ionizing SLRs with time within the disk chemistry code, we have created a model where the SLR ionization rate is now time-dependent internal to the chemical calculations and the total (ensemble) rate decays with a 1.2 Myr half-life (Appendix A), as a simple first order approximation. In reality, the specific radionuclide decay products (i.e., -rays, particles, etc.) will evolve as each parent nuclide decays, i.e., Al versus Fe and will result in different amounts of energy deposited or lost depending upon the ionization cross sections and the decay rate of the parent nuclide. The changeover from a Al dominated SLR rate to Fe rate happens at around 5 Myr, and so strictly speaking a more detailed treatment that is beyond the scope of this paper would separate the individual contributions. We note that in addition to time variation, there will be variation between the starting abundances of SLRs between disks. In this analysis we assume solar nebula-like abundances from Cleeves et al. (2013b).

Figure 11 shows the effect of time-decaying SLRs on molecular column densities for the two lowest CR ionization models, SSX and TTX, where the change is most significant. The Figure 11 column densities are shown at  Myr of time evolution instead of 1 Myr, which corresponds to a decrease of nearly an order of magnitude in the SLR ionization rate. We show the model for fixed initial SLR abundances (blue), for CR only (orange), and a time-decaying SLR ionization rate (magenta). The inclusion of time-dependence for the SLRs does not sensitively affect the SSX model, which carries similar contribution from both CR and SLR ionization (at  Myr). In other words, in the absence of SLRs, the CRs provide similar levels of ionization in the SSX case. The CR “absent” TTX models (Figure 11, right) are far more sensitive to the change in SLR rate, where the time decay changes the column densities of ND, HD, and HD by at least an order of magnitude. Note the HCO column is insensitive to differences in the SLR model assumptions, regardless of CR rate, as is expected by an X-ray dominated chemistry.

Figure 11.— Column density dependence on radionuclide ionization now including time dependence. All curves are shown at  Myr of chemical evolution. The blue line shows a model with CRs and a fixed SLR ionization rate set by the initial value, no time decay included. Magenta curves include a decay on the SLR ionization rate with a half-life of  Myr (Cleeves et al., 2013b). The yellow line shows the column density in the absence of SLR ionization. The SSX model (as well as the higher CR ionization models, M02 and W98) do not depend on the time decay of SLRs as expected. The TTX model is not buffered by CR ionization and is hence far more sensitive to the time decay of SLRs, except in the case of HCO.

6.3. Dependence on Disk Mass

The results of our study are for a particular disk model (§2) and disk mass ( M). Observed protoplanetary disks show significant diversity in mass, diameter, stellar properties and environment (Williams & Cieza, 2011). It is therefore interesting to quantify the sensitivity of these results to the specific choice of disk mass. To test this scenario, we take our standard model and find the chemical signatures for a disk of half mass ( M) and double mass ( M). To facilitate a simple comparison, we do not change the UV radiation field or temperature structure, but we do recompute the X-ray and CR ionization field. In reality, an increase or decrease by a factor of in density would change the disk opacity and would result in a correspondingly cooler (warmer) disk by in dust temperature. In the present section we focus on the abundances determined from the column densities; however, we note that the emission line ratios for a particular species will also reflect the change in local temperature due to the mass change. Furthermore, a larger (smaller) mass also would make for a more (less) UV/X-ray shielded disk. However, by changing a single parameter we may investigate, in this case, the role of more or less efficient ion-recombination and how this effect plays into the measured abundances (see §6.4 for temperature dependence). We leave the geometrical parameters of the disk unchanged such that the disk density scales with the change in mass. CRs are slightly more excluded by the higher gas column (not a large overall effect) and we do not consider SLRs here. A higher disk mass will trap more SLR decay products prior to loss but it is only a small effect ( more ionization; Cleeves et al., 2013b).

In Figure 12, we show the column of the indicated molecular ion, , normalized to the column of model CO, , which acts as our observable mass-reference. Observations of optically thin lines provide their respective total (line-of-sight integrated) column densities. However, when one wants to determine relative abundances, typically to H, a normalization quantity is required. The total gas mass can be inferred from millimeter dust emission or from molecular gas traced by CO, in both cases requiring a calibrated conversion factor. Both methods have substantial caveats. The millimeter-wave dust emission evolves as the underlying population of dust evolves via grain growth, making the dust-to-gas conversion factor a time-dependent quantity. The CO abundance relative to H likewise can be unreliable, both owing to freeze-out at low temperatures and chemical processing initiated by ionization (Bergin et al., 2014), though this process may be slowed by vertical mixing (Furuya & Aikawa, 2014) and by grain growth through decreased surface area for freeze-out (Bergin et al., 2014). Ideally a less chemically reactive mass tracer like HD should be used for normalization, if available (Bergin et al., 2013). Nonetheless, given limited data on HD in protoplanetary disks, we provide column densities relative to the CO traced gas column, a far more widely available gas mass probe. In practice, CO column densities are extracted from optically thin isotopologues such as CO, however we emphasize that even CO has a minimum optical depth of at line center in this particular model (see Fig. 5). To determine the CO column observationally, one must consider either i) more rarified CO isotopes than CO or ii) emission restricted to the line-wings to ensure the line is thin. In this analysis, the column densities of both the ions and of CO are taken directly from the chemical models.

Figure 12.— Column-derived abundances of the indicated molecular ions, normalized to CO as a function of disk radius. Colors represent disk mass, where yellow, blue and magenta correspond to , , and disk gas mass models. CR models are as indicated in the column headings. Changes in the normalized column densities of a given molecular ion for different disk masses are, in general, far smaller than changes across different CR ionization rates. Ionization decreases as one goes from left to right. See text for details and discussion.

In general, for the intermediate ionization rate models, W98 and SSX, the molecular ion abundance is not very sensitive to changes in the disk mass (formally density). Most importantly, the abundances are far more sensitive to the CR ionization rate than they are to the disk mass. Nonetheless, some variation does exist, where in Figure 12 the line thickness increases for increasing disk mass. HCO in particular shows a decrease in abundance with disk mass across all CR models. These results can be understood by the same relation discussed in §6.1, where the ion abundance is proportional to . In terms of column densities, the ratio of the ion-species to CO is then approximately given by . For a more massive disk, decreases due to increased gas opacity, and increases as the mass increases. Thus the quantity will decrease for increasing mass, reducing the ion abundance, explaining the general trend for all of the molecular ions considered here.

Two special cases are the high M02 CR rate (leftmost column) and the low TTX CR rate (rightmost column). Generally speaking, the M02 CR models tend to damp out the change in X-ray ionization, , leaving as the dominant term. M02 shows substantially more variation than can be attributed to changes in or in X-ray flux, otherwise the variation would be seen in W98 as well. This behavior can be understood by looking towards the denominator of the abundance ratio: . As seen in the column density plots in Figure 4, the CO abundance is decreased due to chemical processing, and this deficit is reflected in the column integrated abundance, especially for HD and HD, which are far more spread out in abundance than the W98 and SSX models in Figure 12.

TTX is also a special case because the molecular ion abundances are very sensitive to the X-ray ionization rate, its only ionizing source. Thus TTX is far more sensitive to the mass of the disk because depends sensitively on both and . Regardless, we note that the abundance spread over different mass models – even for the TTX case – is far less than the spread resulting from different ionization rates. We thus conclude that ionization (X-ray, CR and SLRs) is the more important quantity regulating ion abundances, not disk gas mass.

6.4. Dependence on Stellar Spectral Type

The chemical properties of disks are especially sensitive to the dust and gas temperatures, which are set by irradiation from the central star at radii beyond the inner  AU, i.e., the region where accretion heating becomes significant. The temperature in the warm molecular layer is especially important in regulating the observable gas-phase chemistry (Aikawa et al., 2002). To explore the chemical dependence on the temperature of the central star, we recompute our baseline chemical model for a central star with T K, where all other structural parameters are held fixed. At  AU, the disk surface () temperature increases substantially from  K to  K. The warm molecular layer () increases from  K to 55 K, and the midplane increases by  K from 30 K to 40 K. This substantial increase in temperature has a sizable effect on the abundance profiles as shown in Figure 13.

Figure 13.— Vertical column densities of select species as a function of stellar spectral type. Solid lines correspond to the fiducial disk model with stellar effective temperature of  K. Dashed lines show column densities for the warmer (more luminous) central star,  K. See discussion in §6.4 for details.

Specifically, the increase in disk temperature pushes out the CO snow line, resulting in a larger abundance of CO in the inner disk, and a thicker layer of CO radially across the entire disk. The HCO column densities in Figure 13 for the M02 and W98 (high CR flux) cases reflect the enhanced HCO production from the additional CO. However, the effect is mediated by the same chemical processing of CO as discussed in §4.2 and §6.3, where for the warmer star, CO processing occurs deeper in the dense layers of the disk, closer to the midplane, and thereby has a larger effect on the total column of CO than for the cooler star. HCO for the low CR-ionization cases, SSX and TTX, shows a decrease in the outer disk, beyond  AU, due to the combination of CO reprocessing happening deeper in the warmer disk and the simultaneous loss of non-thermal CR desorption, enabling rapid carbon sequestration from CO into other ices.

For all CR-ionization models, the column density of NH shows a decrease in the inner  AU and a surplus beyond  AU in the cold versus warm disk model. The inner deficit is a direct result of the enhanced inner disk CO abundance and increased destruction rate. The outer enhancement follows from the increased gas-phase N abundance in the warmer disk model where it would otherwise freeze out, thereby enabling NH formation. This process operates in tandem with the loss of CO in the dense gas due to reprocessing, reducing one of the primary destruction agents of NH. The ND column density shows similar morphological changes as NH where it is pushed further out radially; however, its abundance is simultaneously affected by the net reduction in deuterium-bearing species in the warm gas and the immediate loss of HD relative to H. The HD column density has the same drop in the inner disk due to the increased CO abundance and overall warmer gas temperatures, though the change is less severe in the outer disk mainly because its precursor, H, does not freeze-out in either model.

The warmer disk model not only changes the chemical properties of the disk, but also the observational strategy necessary to experimentally determine the disk ionization processes, especially that of CRs. In the cold disk case, all molecular ions considered here were able to distinguish between high (M02 and W98) and low (SSX and below) CR fluxes outside of  AU; however, the column density of these species are less sensitive to changes at lower ionization rates due to freeze-out of CO and N. For the warm disk, the CO freeze-out region is pushed further out radially and the HCO column density becomes more sensitive to the CR ionization rate outside of the X-ray dominated region ( AU). At  AU, the SSX and TTX CR ionization rates are indistinguishable for the cold disk model because of CO freeze-out, but in the warm disk case, there is a difference between the SSX and TTX models. While this difference is not sufficient to measure the ionization rate to great accuracy for this particular model, a warm disk around a more X-ray faint star (or a less X-ray permeable disk) may engender conditions favorable for optically thin HCO isotopologue emission tracing CR-dominated layers down to low CR flux-levels. NH and ND remain sensitive observational tracers of the CR ionization rate at both high and low CR fluxes beyond  AU, where the warm disk acts to increase the “dynamic range” between the model column densities due to the reduced N freeze-out. For example, the variation in NH column density between SSX and TTX models was approximately one order of magnitude or less in the fiducial cold disk (see Figure 4). The warm disk increases the fractional difference to over two orders of magnitude in NH and the same behavior holds true for ND. Thus warmer disks may allow us to estimate more precisely the CR ionization rate by mediating the effects of freeze-out and increasing the overall column. One potential caveat is that the radially larger “snowline” makes for a higher disk-averaged CO column and disk-integrated CO opacity compared to a cooler disk. In this case, even rarer CO isotopes than CO may be necessary to interpret the HCO column densities or to use as a proxy for gas mass. It is also important to point out that the emission is not only sensitive to the column but also the emitting temperature, and therefore the particular transition or transitions targeted should also take into account the inherent temperature of the disk as estimated by the stellar luminosity.

6.5. Disk Magnetic “Opacity” to Cosmic Ray Ionization

The CR contribution to disk ionization in both magnitude and scope is the least observationally constrained parameter in disks. A cosmic ray ionization rate of  s, consistent with dense ISM values, was ruled out by Chapillon et al. (2011) using observations of HD. There are a few possible explanations for the low CR rates inferred. The present work focuses on the possibility of wind-modulated incident CR ionization rates due to the presence of an analogue to the Solar System’s heliosphere. Within this paradigm we have ignored radial variations in the GCR rate that may include (i) a gradual increase (1%/AU) in the CR rate with radius where modulation is weaker further out in the disk (Paper I) and (ii) an edge to the region of CR-modulation, i.e, a Heliopause. As compared to the negative ionization gradient expected for SLR ionization, (i) would create a positive ionization gradient, though the two effects are similar ( order of magnitude) and may conspire to cancel each other out. The other important radial effect is the extent of the wind-modulation zone, the “T-Tauriosphere.” We have assumed the disk is fully enclosed, but if winds can only punch out the inner tens of AU, perhaps by magnetic trapping of the winds (Turner et al., 2014), the region of exclusion may be much smaller. Alternatively, if the winds dominate, they could perhaps encircle a much larger region, hundreds to thousands of AU in size due to the high gas densities (and correspondingly higher ram pressures) of early stellar winds.

Magnetic fields provide an additional source of “opacity” to CRs in two ways: (i) for funnel-shaped magnetic field configurations seen in protostars (e.g., Girart et al., 2006), the magnetic field can mirror at most 50% of the CRs away (see also Padovani & Galli, 2011) and (ii) the presence of magnetic irregularities with size scales near the CR gyro-radius can scatter CRs (Cesarsky & Volk, 1978). Both effects can act in tandem, where magnetic irregularities on an hourglass magnetic field configuration can further enhance the fraction of mirrored particles (Fatuzzo & Adams, 2014). The role of CR exclusion by magnetic irregularities induced by disk gas turbulence has been explored in Dolginov & Stepinski (1994). It was found that modest magnetic field strengths with imposed turbulent irregularities significantly impeded the CR rate in the disk, such that only of CRs reach depths corresponding to the disk scale height at all radii beyond the inner few AU. Including irregularities causes the CR ionization rate to decay very quickly with vertical depth towards the midplane (see Eq. (9) of Dolginov & Stepinski, 1994), thus removing CRs entirely by the midplane. Plugging in the numbers typical of our own disk model, we find that magnetic irregularities act to reduce the CR rate at the midplane by six orders of magnitude at 100 AU compared to the CR rate at the disk surface. Furthermore, both winds and magnetic effects can operate simultaneously, such that the winds reduce the incident CR ionization rate and disk magnetic irregularities substantially curtail the CR propagation internal to the disk, much faster the classically assumed penetration depth of 100 g cm (Umebayashi & Nakano, 1981).

The two scenarios, winds and “magnetic opacity,” act similarly to reduce the CR flux in the disk’s midplane but where they differ is in the surface. While we have touted the HCO isotopologue emission as an excellent X-ray tracer, it is still sensitive to the higher CR ionization rate models considered here (with  s, i.e., M02 and W98). If winds are modulating the CR ionization, we would expect the HCO emission to reflect a uniformly low CR rate regardless of height. If magnetic irregularities dominate the attenuation of CRs, then the CR ionization rate should be normal in the HCO traced upper layers and absent in the midplane. To conduct such an experiment however, requires reasonably good constraints on the stellar X-ray luminosity and distributions of density and temperature in the disk, as both of these are expected to affect the HCO emission (§6.1 and §6.3, respectively). Alternatively, if neither of these effects are important, wind modulation nor magnetic effects, then the ionization should be  s in the outer disk and the more weakly emissive, difficult species to observe, like HD, can help provide constraints on midplane ionization.

7. Discussion and Conclusions

Using a generic, observationally motivated model of a T Tauri protoplanetary disk, we present chemical abundance models and simulated submillimeter emission line observations that can be used as a blueprint to constrain the detailed ionization environment within the gas disk. In particular, sensitive () ALMA observations of nearby protoplanetary disks ( pc) will readily be able to distinguish between systems with high, interstellar CR ionization levels ( s) and those with low (sub-interstellar) ionization levels, though determining very low CR rates ( s) will be made difficult by weak emission from molecular tracers and X-rays providing a lower ionization limit that may hide the effects of CRs. We emphasize that the chemical results presented here demonstrate relative trends across different ionization models and that the physical structure of the disk and properties of the star will determine how the chemical abundances measured directly map to ionization properties of the disk in an absolute sense. Consequently, creating detailed models of particular sources with as many observational constraints as possible will be crucial to mapping out the ionization in detail for any particular system.

We highlight the molecular ions that are useful tracers of specific ionizing agents in dense gas, e.g., the warm, X-ray irradiated molecular surface through HCO or the cold, dense SLR and/or CR dominated midplane through ND and HD. Moreover, by isolating individual ionization sources, whether by central stellar processes or otherwise, we can use these results to observationally quantify the relative importance of each ionizing agent to the total ionization state of the disk gas. Better constraints on the underlying ionization environment inform models of turbulence via the ionization-dependent magneto-rotational instability (Balbus & Hawley, 1991) and models of disk chemistry via ion-neutral and grain surface reactions. In other words, dense gas ionization drives both the fundamental processes that govern both the formation (by potentially regulating turbulence-free “dead zones”) and chemical make-up of planetary systems.

In particular, the least observationally constrained ionizing agents present in disks are those which dominate the midplane ionization – CRs and SLRs. More specifically, the CR rate incident on a protoplanetary disk is unknown and may be strongly modulated by winds and/or magnetic fields. Without CRs, the total midplane ionization rate is reduced by at least one to two orders of magnitude depending upon radial location and the degree to which SLR ionization contributes, which becomes the primary midplane ionization source in the absence of CRs (see also Paper I). High spatial resolution observations may furthermore reveal radial structure (e.g., gradients) in the CR rate with distance from the star, potentially belying the presence of an analogue heliosphere to that of the Solar System. The molecular tracers outlined in this paper thus provide signposts of the presence or absence of important physical processes related to disk ionization. In summary, the main results of this work are:

  1. Chemical abundances of different molecular ions trace different ionizing agents (and regions) in our fiducial disk model. The abundance variation is born out in their submillimeter emission, and thus the lines discussed here can be used as diagnostic observational tools of disk ionization.

  2. Optically thin isotopes of HCO trace primarily the incident X-ray ionizing flux, in particular X-rays with energies  keV. For a generally warmer disk, the CO snow line may exist further out, and in this scenario, radially extended (resolved) HCO emission may be sensitive to midplane ionization sources, i.e., CRs and SLRs.

  3. NH is sensitive to both cold ionization processes (CRs and SLRs) and warm ( K) gas ionization via stellar X-rays and is likely moderately optically thick inside  AU ( of a few). This behavior makes for somewhat difficult direct interpretation, requiring detailed chemical/line emission models and multiple transitions.

  4. We find that the cold gas ionization tracer HD becomes undetectable in a full ALMA observation at  pc for CR ionization rates below solar maximum (SSX models). To detect the amount of HD present in our models at the level of  s, at least an order of magnitude more sensitive observations would be required, or a less distant disk such as TW Hya.

  5. ND, rather than HD is a more sensitive tracer of midplane ionization due to CRs and SLRs and the least sensitive to stellar X-rays. ND can be used to measure ionization levels even at low CR rates ( s) and should be detectable with ALMA in a reasonable amount of time ().

  6. Cosmic rays may be excluded by wind processes or internal/external magnetic processes (or both), but combined observations of surface ionization tracers such as HCO isotopologues and midplane tracers like ND will help illuminate the CR-exclusion mechanism that dominates (see §6.5).

  7. Short-lived radionuclides are an important contributor to the ionization in disks, exceeding that of scattered stellar X-rays within the first few Myr. If CRs are not present, ND should still be detectable for young disks due to the SLR contribution, but its emission will fade over the disk lifetime.

  8. It will be very difficult to tell apart SLR ionization from CR ionization for measured midplane ionization rates of  s. A negative ionization gradient with radius would indicate a SLR dominated chemistry (SLRs can escape the tenuous outer disk radii), while a flat or positive gradient (if the inner disk is very dense with  g cm) may point to a CR dominated chemistry.

  9. The mass-normalized column density of ions (with respect to HD, CO, or dust) is far more sensitive to the CR ionization rate than to the mass of the disk itself, where we have considered disk masses spanning the range  M.

  10. Not all emission “rings” trace physical deficits in abundance or structure. Large temperature gradients present in disks can result in low-J rotational lines peaking offset from the star due to de-excitation in warm gas, such as can happen with HCO for high CR rates.

This paper provides a viable starting point to study the ionization sources acting in circumstellar disks and their corresponding chemical signatures. However, this work must be carried forward in several ways. First, the ionization models used here are preliminary. We need to construct more detailed models for how CRs are suppressed by both T Tauri winds and magnetic field fluctuations; we also need improved assessments of the ionization rates provided by background cluster environments. On another front, both existing and upcoming submillimeter facilities will provide important constraints on the actual ionization levels realized in these systems and will determine which molecules provide the most information. With improved theoretical and observational input, the chemical signatures considered here can then be revisited. In the end, we will thus obtain a good working understanding of both the relevant ionization processes and the chemical structure of planet forming disks. This information, in turn, can then be used to constrain disk evolution. More specifically, the chemical structure of the disk determines the locations of both the dead zones (where MRI cannot operate) and chemical gradients in the gas and ice, and these structures greatly influence the accompanying processes of disk accretion and planet formation.

Acknowledgements: The authors thank the anonymous referee and editor for their comments and suggestions, which have greatly improved the manuscript. This work was supported by NSF grant AST-1008800.


  • Adams (2010) Adams, F. C. 2010, ARA&A, 48, 47
  • Adams et al. (2014) Adams, F. C., Fatuzzo, M., & Holden, L. 2014, ApJ, 789, 86
  • Aikawa & Herbst (1999a) Aikawa, Y. & Herbst, E. 1999a, ApJ, 526, 314
  • Aikawa & Herbst (1999b) —. 1999b, A&A, 351, 233
  • Aikawa et al. (2001) Aikawa, Y., Ohashi, N., Inutsuka, S.-i., Herbst, E., & Takakuwa, S. 2001, ApJ, 552, 639
  • Aikawa et al. (2002) Aikawa, Y., van Zadelhoff, G. J., van Dishoeck, E. F., & Herbst, E. 2002, A&A, 386, 622
  • Albertsson et al. (2014) Albertsson, T., Semenov, D., & Henning, T. 2014, ApJ, 784, 39
  • Amano & Hirao (2005) Amano, T. & Hirao, T. 2005, Journal of Molecular Spectroscopy, 233, 7
  • Anderson et al. (2013) Anderson, D. E., Bergin, E. A., Maret, S., & Wakelam, V. 2013, ApJ, 779, 141
  • Anderson et al. (1977) Anderson, T. G., Dixon, T. A., Piltch, N. D., Saykally, R. J., Szanto, P. G., & Woods, R. C. 1977, ApJ, 216, L85
  • Andrews et al. (2011) Andrews, S. M., Wilner, D. J., Espaillat, C., Hughes, A. M., Dullemond, C. P., McClure, M. K., Qi, C., & Brown, J. M. 2011, ApJ, 732, 42
  • Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481
  • Bachiller et al. (1990) Bachiller, R., Menten, K. M., & del Rio Alvarez, S. 1990, A&A, 236, 461
  • Balbus & Hawley (1991) Balbus, S. A. & Hawley, J. F. 1991, ApJ, 376, 214
  • Bergin et al. (2002) Bergin, E. A., Alves, J., Huard, T., & Lada, C. J. 2002, ApJ, 570, L101
  • Bergin et al. (2014) Bergin, E. A., Cleeves, L. I., Crockett, N., & Blake, G. A. 2014, Faraday Discussions, Accepted
  • Bergin et al. (2013) Bergin, E. A., Cleeves, L. I., Gorti, U., Zhang, K., Blake, G. A., Green, J. D., Andrews, S. M., Evans, II, N. J., Henning, T., Öberg, K., Pontoppidan, K., Qi, C., Salyk, C., & van Dishoeck, E. F. 2013, Nature, 493, 644
  • Bergin et al. (1999) Bergin, E. A., Neufeld, D. A., & Melnick, G. J. 1999, ApJ, 510, L145
  • Bethell & Bergin (2011a) Bethell, T. J. & Bergin, E. A. 2011a, ApJ, 740, 7
  • Bethell & Bergin (2011b) —. 2011b, ApJ, 739, 78
  • Brinch & Hogerheijde (2010) Brinch, C. & Hogerheijde, M. R. 2010, A&A, 523, A25
  • Calvet et al. (2005) Calvet, N., D’Alessio, P., Watson, D. M., Franco-Hernández, R., Furlan, E., Green, J., Sutter, P. M., Forrest, W. J., Hartmann, L., Uchida, K. I., Keller, L. D., Sargent, B., Najita, J., Herter, T. L., Barry, D. J., & Hall, P. 2005, ApJ, 630, L185
  • Ceccarelli et al. (2004) Ceccarelli, C., Dominik, C., Lefloch, B., Caselli, P., & Caux, E. 2004, ApJ, 607, L51
  • Cesarsky & Volk (1978) Cesarsky, C. J. & Volk, H. J. 1978, A&A, 70, 367
  • Chapillon et al. (2011) Chapillon, E., Parise, B., Guilloteau, S., & Du, F. 2011, A&A, 533, A143
  • Cleeves et al. (2013a) Cleeves, L. I., Adams, F. C., & Bergin, E. A. 2013a, ApJ, 772, 5
  • Cleeves et al. (2013b) Cleeves, L. I., Adams, F. C., Bergin, E. A., & Visser, R. 2013b, ApJ, 777, 28
  • Cohen et al. (2012) Cohen, O., Drake, J. J., & Kóta, J. 2012, ApJ, 760, 85
  • Collings et al. (2004) Collings, M. P., Anderson, M. A., Chen, R., Dever, J. W., Viti, S., Williams, D. A., & McCoustra, M. R. S. 2004, MNRAS, 354, 1133
  • Consolmagno & Jokipii (1978) Consolmagno, G. J. & Jokipii, J. R. 1978, The moon and the planets, 19, 253
  • Dalgarno (2006) Dalgarno, A. 2006, Proceedings of the National Academy of Science, 103, 12269
  • Dolginov & Stepinski (1994) Dolginov, A. Z. & Stepinski, T. F. 1994, ApJ, 427, 377
  • Draine (1978) Draine, B. T. 1978, ApJS, 36, 595
  • Dullemond & Dominik (2004) Dullemond, C. P. & Dominik, C. 2004, A&A, 421, 1075
  • Dutrey et al. (1997) Dutrey, A., Guilloteau, S., & Guelin, M. 1997, A&A, 317, L55
  • Ercolano & Glassgold (2013) Ercolano, B. & Glassgold, A. E. 2013, MNRAS, 436, 3446
  • Fatuzzo & Adams (2008) Fatuzzo, M. & Adams, F. C. 2008, ApJ, 675, 1361
  • Fatuzzo & Adams (2014) —. 2014, ApJ, 787, 26
  • Favre et al. (2013) Favre, C., Cleeves, L. I., Bergin, E. A., Qi, C., & Blake, G. A. 2013, ApJ, 776, L38
  • Feigelson et al. (1993) Feigelson, E. D., Casanova, S., Montmerle, T., & Guibert, J. 1993, ApJ, 416, 623
  • Feigelson & Decampli (1981) Feigelson, E. D. & Decampli, W. M. 1981, ApJ, 243, L89
  • Finocchi & Gail (1997) Finocchi, F. & Gail, H.-P. 1997, A&A, 327, 825
  • Flower (1999) Flower, D. R. 1999, MNRAS, 305, 651
  • Flower et al. (2006) Flower, D. R., Pineau Des Forêts, G., & Walmsley, C. M. 2006, A&A, 449, 621
  • Fogel et al. (2011) Fogel, J. K. J., Bethell, T. J., Bergin, E. A., Calvet, N., & Semenov, D. 2011, ApJ, 726, 29
  • France et al. (2014) France, K., Schindhelm, E., Bergin, E. A., Roueff, E., & Abgrall, H. 2014, ApJ, 784, 127
  • Fromang et al. (2002) Fromang, S., Terquem, C., & Balbus, S. A. 2002, MNRAS, 329, 18
  • Furuya & Aikawa (2014) Furuya, K. & Aikawa, Y. 2014, ArXiv e-prints
  • Gammie (1996) Gammie, C. F. 1996, ApJ, 457, 355
  • Garrod & Pauly (2011) Garrod, R. T. & Pauly, T. 2011, ApJ, 735, 15
  • Geiss & Gloeckler (2003) Geiss, J. & Gloeckler, G. 2003, Space Sci. Rev., 106, 3
  • Girart et al. (2006) Girart, J. M., Rao, R., & Marrone, D. P. 2006, Science, 313, 812
  • Glassgold (1999) Glassgold, A. E. 1999, in IAU Symposium, Vol. 191, Asymptotic Giant Branch Stars, ed. T. Le Bertre, A. Lebre, & C. Waelkens, 337
  • Glassgold et al. (2004) Glassgold, A. E., Najita, J., & Igea, J. 2004, ApJ, 615, 972
  • Guélin (1988) Guélin, M. 1988, in Topics in Molecular Organization and Engineering, Vol. 2, Molecules in Physics, Chemistry, and Biology, ed. J. Maruani (Springer Netherlands), 175–187
  • Guilloteau et al. (2006) Guilloteau, S., Piétu, V., Dutrey, A., & Guélin, M. 2006, A&A, 448, L5
  • Habing (1968) Habing, H. J. 1968, Bull. Astron. Inst. Netherlands, 19, 421
  • Harries (2000) Harries, T. J. 2000, MNRAS, 315, 722
  • Harries et al. (2004) Harries, T. J., Monnier, J. D., Symington, N. H., & Kurosawa, R. 2004, MNRAS, 350, 565
  • Hasegawa et al. (1992) Hasegawa, T. I., Herbst, E., & Leung, C. M. 1992, ApJS, 82, 167
  • Hatchell et al. (1998) Hatchell, J., Thompson, M. A., Millar, T. J., & MacDonald, G. H. 1998, A&A, 338, 713
  • Herbst & Klemperer (1973) Herbst, E. & Klemperer, W. 1973, ApJ, 185, 505
  • Herczeg et al. (2002) Herczeg, G. J., Linsky, J. L., Valenti, J. A., Johns-Krull, C. M., & Wood, B. E. 2002, ApJ, 572, 310
  • Herczeg et al. (2004) Herczeg, G. J., Wood, B. E., Linsky, J. L., Valenti, J. A., & Johns-Krull, C. M. 2004, ApJ, 607, 369
  • Hughes et al. (2011) Hughes, A. M., Wilner, D. J., Andrews, S. M., Qi, C., & Hogerheijde, M. R. 2011, ApJ, 727, 85
  • Hughes et al. (2008) Hughes, A. M., Wilner, D. J., Qi, C., & Hogerheijde, M. R. 2008, ApJ, 678, 1119
  • Hugo et al. (2009) Hugo, E., Asvany, O., & Schlemmer, S. 2009, J. Chem. Phys., 130, 164302
  • Igea & Glassgold (1999) Igea, J. & Glassgold, A. E. 1999, ApJ, 518, 848
  • Indriolo & McCall (2012) Indriolo, N. & McCall, B. J. 2012, ApJ, 745, 91
  • Jonkheid et al. (2007) Jonkheid, B., Dullemond, C. P., Hogerheijde, M. R., & van Dishoeck, E. F. 2007, A&A, 463, 203
  • Jørgensen et al. (2004) Jørgensen, J. K., Schöier, F. L., & van Dishoeck, E. F. 2004, A&A, 416, 603
  • Joseph et al. (1986) Joseph, C. L., Snow, Jr., T. P., Seab, C. G., & Crutcher, R. M. 1986, ApJ, 309, 771
  • Jura et al. (2013) Jura, M., Xu, S., & Young, E. D. 2013, ApJ, 775, L41
  • Kahane et al. (1992) Kahane, C., Cernicharo, J., Gomez-Gonzalez, J., & Guelin, M. 1992, A&A, 256, 235
  • Keene et al. (1998) Keene, J., Schilke, P., Kooi, J., Lis, D. C., Mehringer, D. M., & Phillips, T. G. 1998, ApJ, 494, L107
  • Kurosawa et al. (2004) Kurosawa, R., Harries, T. J., Bate, M. R., & Symington, N. H. 2004, MNRAS, 351, 1134
  • Kutschera et al. (1984) Kutschera, W., Billquist, P. J., Frekers, D., Henning, W., Jensen, K. J., Xiuzeng, M., Pardo, R., Paul, M., Rehm, K. E., Smither, R. K., Yntema, J. L., & Mausner, L. F. 1984, Nuclear Instruments and Methods in Physics Research B, 5, 430
  • Langer et al. (1996) Langer, W. D., Castets, A., & Lefloch, B. 1996, ApJ, 471, L111
  • Lee & Bergin (2014) Lee, J.-E. & Bergin, E. A. 2014, ApJ, submitted.
  • MacPherson et al. (1995) MacPherson, G. J., Davis, A. M., & Zinner, E. K. 1995, Meteoritics, 30, 365
  • Maret & Bergin (2007) Maret, S. & Bergin, E. A. 2007, ApJ, 664, 956
  • Matsumura & Pudritz (2003) Matsumura, S. & Pudritz, R. E. 2003, ApJ, 598, 645
  • Millar et al. (1989) Millar, T. J., Bennett, A., & Herbst, E. 1989, ApJ, 340, 906
  • Moskalenko et al. (2002) Moskalenko, I. V., Strong, A. W., Ormes, J. F., & Potgieter, M. S. 2002, ApJ, 565, 280
  • Müller et al. (2005) Müller, H. S. P., Schlöder, F., Stutzki, S., & Winnewisser, G. 2005, Journal of Molecular Structure, 742, 215
  • Müller et al. (2001) Müller, H. S. P., Thorwirth, S., Roth, D. A., & Winnewisser, G. 2001, A&A, 370, L49
  • Öberg et al. (2010) Öberg, K. I., Qi, C., Fogel, J. K. J., Bergin, E. A., Andrews, S. M., Espaillat, C., van Kempen, T. A., Wilner, D. J., & Pascucci, I. 2010, ApJ, 720, 480
  • Öberg et al. (2011a) Öberg, K. I., Qi, C., Fogel, J. K. J., Bergin, E. A., Andrews, S. M., Espaillat, C., Wilner, D. J., Pascucci, I., & Kastner, J. H. 2011a, ApJ, 734, 98
  • Öberg et al. (2011b) Öberg, K. I., Qi, C., Wilner, D. J., & Andrews, S. M. 2011b, ApJ, 743, 152
  • Öberg et al. (2005) Öberg, K. I., van Broekhuizen, F., Fraser, H. J., Bisschop, S. E., van Dishoeck, E. F., & Schlemmer, S. 2005, ApJ, 621, L33
  • Padovani & Galli (2011) Padovani, M. & Galli, D. 2011, A&A, 530, A109
  • Pagani et al. (2009) Pagani, L., Vastel, C., Hugo, E., Kokoouline, V., Greene, C. H., Bacmann, A., Bayet, E., Ceccarelli, C., Peng, R., & Schlemmer, S. 2009, A&A, 494, 623
  • Perryman et al. (1997) Perryman, M. A. C., Lindegren, L., Kovalevsky, J., Hoeg, E., Bastian, U., Bernacca, P. L., Crézé, M., Donati, F., Grenon, M., Grewing, M., van Leeuwen, F., van der Marel, H., Mignard, F., Murray, C. A., Le Poole, R. S., Schrijver, H., Turon, C., Arenou, F., Froeschlé, M., & Petersen, C. S. 1997, A&A, 323, L49
  • Pickett et al. (1998) Pickett, H. M., Poynter, R. L., Cohen, E. A., Delitsky, M. L., Pearson, J. C., & Müller, H. S. P. 1998, J. Quant. Spec. Radiat. Transf., 60, 883
  • Pinte et al. (2009) Pinte, C., Harries, T. J., Min, M., Watson, A. M., Dullemond, C. P., Woitke, P., Ménard, F., & Durán-Rojas, M. C. 2009, A&A, 498, 967
  • Prantzos et al. (1996) Prantzos, N., Aubert, O., & Audouze, J. 1996, A&A, 309, 760
  • Qi et al. (2008) Qi, C., Wilner, D. J., Aikawa, Y., Blake, G. A., & Hogerheijde, M. R. 2008, ApJ, 681, 1396
  • Roberts et al. (2004) Roberts, H., Herbst, E., & Millar, T. J. 2004, A&A, 424, 905
  • Roberts & Millar (2000) Roberts, H. & Millar, T. J. 2000, A&A, 361, 388
  • Rugel et al. (2009) Rugel, G., Faestermann, T., Knie, K., Korschinek, G., Poutivtsev, M., Schumann, D., Kivel, N., Günther-Leopold, I., Weinreich, R., & Wohlmuther, M. 2009, Physical Review Letters, 103, 072502
  • Rydbeck et al. (1980) Rydbeck, O. E. H., Hjalmarson, A., Rydbeck, G., Ellder, J., Kollberg, E., & Irvine, W. M. 1980, ApJ,