EXCLUSION OF COSMIC RAYS IN PROTOPLANETARY DISKS. II. CHEMICAL GRADIENTS AND OBSERVATIONAL SIGNATURES

T Tauri stars are known to be X-ray luminous with X-ray luminosities typically ranging L_XR ∼ 10²⁹–10³¹ 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 E_XR = 1 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 Section 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 Lyα line transfer, where photons scatter isotropically off atomic hydrogen present at the disk surface. Treatment of Lyα is important as this line has been observed to carry ∼80%–90% 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 Lyα 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, d = 55 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 G₀ (Fatuzzo & Adams 2008), where G₀ is the mean value for the ISM (Habing 1968; Draine 1978). A G₀ = 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 Figure 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 backward and forward directions). For the primary model, we assume a characteristic T Tauri star X-ray luminosity of L_XR = 10^29.5 erg s⁻¹ with the quiescent X-ray spectral template given in Paper I.

The most abundant low energy X-rays (E ∼ 1 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 5–7 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 R ∼ 50 AU is simply geometrical dilution. The X-ray ionization floor falls in the range ζ_XR ∼ (1–10) × 10⁻²¹ 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 ζ_XR depending upon the shape of the stellar X-ray spectrum (see also Section 6.1).

Zoom In Zoom Out Reset image size

The disk ionization contribution from the decay of SLRs has been studied extensively (Consolmagno & Jokipii 1978; Umebayashi & Nakano 1981, 2009; Finocchi & Gail 1997; 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 (t_half = 0.74 Myr; Schramm 1971; MacPherson et al. 1995; Umebayashi & Nakano 2009) and ⁶⁰Fe (t_half = 2.6 Myr³; Rugel et al. 2009) corresponds to approximately t_half ∼ 1.2 Myr (Appendix). 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 R > 50 AU (Figure 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 ∼ 1–10 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 Section 4.2 assume the SLR ionization rate is constant with time, which is acceptable for disks with t < 1 Myr or less. We relax this requirement in Section 6.2 where we consider SLR time evolution within the chemical calculation itself.

Ionization by CRs at the interstellar rate (1–5 × 10⁻¹⁷ s⁻¹) is a commonly assumed ingredient in models of disk chemistry and physics. However, detection limits of H₂D⁺ 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 analog to the modern day "Heliosphere." Within the Heliosphere, the solar wind strongly modulates CRs, especially those below E_CR ∼ 100 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 ζ_CR ∼ (3–100) × 10⁻²² 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 ζ_CR ≲ 10⁻¹⁸ 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 ζ_CR ∼ 10⁻¹⁹ 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 ζ_SLR ∼ (1–10) × 10⁻¹⁹ s⁻¹ for t < 1 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 (E > 1 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 ∼ 100 g cm⁻² (Umebayashi & Nakano 1981). Additional opacity can arise from magnetic mirroring and magnetic irregularities in the disk, which we discuss further in Section 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, ζ_CR ∼ 2 × 10⁻¹⁷ s⁻¹, while the M02 ionization rate is characteristic of CR ionization rates measured from H$_3^+$ 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.

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 Ohio State University (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 r_g = 0.1 μ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 1 mm 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 r_g 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., 7.5 × 10⁻³ μm² cm⁻³ at $n_{\rm H_2}=10^{10}$ 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 H₂D⁺ and N₂D⁺, which become enriched relative to the main isotopologue due to a chemical favorability toward the heavier isotopologue at low (T < 50 K) temperatures (Millar et al. 1989). Even though DCO⁺ is part of the network through the H₂D⁺+ 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 H₂D⁺ and closely related N₂D⁺ chemistry. An important facet of the H₂D⁺ chemistry is that the reaction rates depend strongly on the ortho-to-para ratio of the reactants, H₂ and H$_3^+$, 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 H₂D⁺ and include the ortho-to-para fraction as weights on the net reaction rate. In the following prescription we designate $f^1_o$ ($f^1_p$) as the fraction of H₂ in the ortho (para) state, and $f^2_o$ ($f^2_p$) as the fraction of H₂D⁺ in the ortho (para) state. The weighted reaction rate coefficients for (o)rtho and (p)ara H₂ and H₂D⁺ is

$$\begin{eqnarray} R(T) &=& R^{\rm oo}(T)f^1_{\rm o}(T)f^2_{\rm o}(T) + R^{\rm op}(T)f^1_{\rm o}(T)f^2_{\rm p}(T) +..., \nonumber\\ \end{eqnarray} \tag{ 1 }$$

where $f^{1}_o + f^{1}_p = 1$, $f^{2}_o + f^{2}_p = 1$, and R^oo(T) is the reaction rate coefficient for reactions between o-H₂ and o-H₂D⁺, for example. The o/p ratio of H₂ is given by

$$\begin{eqnarray} {\rm o/p}_{\rm H_2}(T_{\rm {gas}})=&5.3534\times 10^{-3} + \frac{(3.0346-5.3534 \times 10^{-3})}{[1+(96.5330/T_{\rm {gas}})^{3.5096}]}, \nonumber\\ \end{eqnarray} \tag{ 2 }$$

for gas at temperatures above T_g ⩾ 80 K and for temperatures below that value,

$$\begin{equation} {\rm o/p}_{\rm H_2}(T_{\rm {gas}})=9 \exp {(-170.5/T_{\rm {gas}})}. \end{equation} \tag{ 3 }$$

Flower et al. (2006) finds that at very low temperatures (T_g ≲ 10 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 Equation (3) from dropping below this value. The o/p ratio of H₂D⁺ is given by

$$\begin{eqnarray} &&{\rm o/p}_{\rm H_2D^+}(T_{\rm {gas}}) = {-}1.6977\times 10^{-2} + \! \frac{(3.0375+1.6977\times 10^{-2})}{ [1+(47.9640/T_{\rm {gas}})^{3.0692}]}, \nonumber\\ \end{eqnarray} \tag{ 4 }$$

where both o/p formalisms in Equations (2) and (4) are from Lee & Bergin (2014).

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 toward 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 H₂D⁺ electron recombination and reactions with atomic H (Roberts et al. 2004), H⁺ and D⁺ reactions with atomic/molecular hydrogen/deuterium, DCO⁺, N₂D⁺ 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₂ + H₂D⁺ reaction rates, the endothermicity of the back reactions for the overall less abundant heavier isotopologues, HD$_2^+$ and D$_3^+$ 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₂, H₂O/HDO, and H₂CO/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 ∼6200 reactions and ∼600 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 ∼10⁻⁸. The atomic (ionized) sulfur, S⁺, is reduced to a low value, where we assume 10⁻¹¹. The motivation behind the total reduced sulfur abundance, ∼10⁻⁸ rather than the diffuse ISM value of ∼10⁻⁶, 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 2 × 10⁻⁹ to 5 × 10⁻⁸ 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 (10⁻¹¹) 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 2.0 ± 0.35 × 10⁻⁵ (Geiss & Gloeckler 2003). The H₂D⁺, HD$_2^+$, and D$_3^+$ abundances are estimated from Walmsley et al. (2004) for the small (0.025 μm) grain case at high ($n_{\rm H_2}=10^7$ 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 N₂H⁺ abundance would increase but the shape of the column density profile remains largely unchanged (Schwarz & Bergin 2014; L. I. Cleeves et al. 2014, in preparation). Starting off with a larger fraction of ammonia, for example, tends to reduce the overall N₂H⁺ 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 N₂H⁺ is still sensitive to the CR ionization rate even without additional information on the disk nitrogen abundances.

Table 3. Initial Chemical Abundances, χ

Species	log	Species	log
	(χ)		(χ)
H2	5.00 × 10−1	H2O(gr)	2.50 × 10−4
He	1.40 × 10−1	N	2.25 × 10−5
CN	6.00 × 10−8	H$_3^+$	1.00 × 10−8
CS	4.00 × 10−9	SO	5.00 × 10−9
Si+	1.00 × 10−11	S+	1.00 × 10−11
Mg+	1.00 × 10−11	Fe+	1.00 × 10−11
C+	1.00 × 10−9	CO	1.00 × 10−4
N2	1.00 × 10−6	C	7.00 × 10−7
NH3	8.00 × 10−8	HCN	2.00 × 10−8
HCO+	9.00 × 10−9	H2CO	8.00 × 10−9
HD	2.00 × 10−5	H2D+	1.30 × 10−10
HD2+	1.00 × 10−10	D3+	2.00 × 10−10
C2H	8.00 × 10−9

Download table as:  ASCIITypeset image

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 t = 1 Myr of chemical evolution are shown in Figure 3. The column headings indicate the underlying CR ionization model as described in Section 3.3 (see also Table 2) with the X-ray luminosity fixed (L_XR = 10^29.5 erg s⁻¹). In addition to the abundances as a function of spatial position, we provide the vertically integrated column density, shown in Figure 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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, H₂D⁺ and N₂D⁺, are highly abundant (χ ∼ 10⁻¹¹ relative to H₂) in the midplane for CR ionization rates exceeding ζ_CR ≳ 10⁻¹⁷ 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 (1) a lack of abundant SLRs or (2) an old (>5 Myr) disk in which the initial reservoir of SLRs has decayed away (<3% remaining), the scattered intermediate energy (5–7 keV) X-ray photons set the absolute floor to the ionization rate, with magnitude ζ_XR ∼ (1–10) × 10⁻²¹ s⁻¹ for a stellar X-ray luminosity of L_XR = 10^29.5 erg s⁻¹. In the TTX model, the midplane abundances of molecular ions are at their minimum (see especially H₂D⁺, HD$_2^+$, and N₂D⁺in Figure 4). Correspondingly, H₂D⁺ and HD$_2^+$ 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 ζ_SLR ∼ (1–10) × 10⁻¹⁹ 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 R ∼ 200 AU. Consequently, distinguishing between a SLR and CR driven chemistry will be challenging. Outside of R ∼ 200 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 H₂D⁺.

Focusing now on the individual ions, the column density of HCO⁺ changes by an order of magnitude between the high (ζ_CR ≳ 10⁻¹⁷ 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 Section 6.1). A second layer of HCO⁺ forms at z/r ∼ 0.2, 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 ζ_CR ≳ 10⁻¹⁷ 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 R ∼ 200 AU not seen for the low CR models. The second peak is actually a consequence of a deficit of CO at R ∼ 50 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 N₂H⁺ 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 N₂H⁺ is that a tenuous (χ ∼ 10⁻¹¹) surface layer of N₂H⁺ is sustained from the high X-ray photon flux (see Figure 3) even in the presence of CO, which is the major reactant of N₂H⁺. For more X-ray luminous stars, the surface N₂H⁺ may contribute more to total column and mask the midplane abundances. Alternatively, the contribution from the surface N₂H⁺ 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 H₂O ice substrate (Collings et al. 2004), but this surface N₂H⁺ is never completely absent in any of our models.

The deuterated molecular ions considered here, H₂D⁺, N₂D⁺  and HD$_2^+$, 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 toward deuterated isotopologues being favored at low temperatures, resulting in their overabundance relative to the main isotopologue in cold (T < 50 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 R ∼ 100 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 (Figure 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.

For our molecular ion abundance results (Section 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 D = 100 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 002 (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 v_dop = 100 m s⁻¹, as observations indicate that the turbulent broadening in disks is small (Hughes et al. 2011).

For the HCO⁺ and N₂H⁺ 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 N₂D⁺ are unknown, we treat the N₂D⁺ level populations as in local thermodynamic equilibrium (LTE), which is a satisfactory approximation as the critical density of the (3–2) and (4–3) transitions are n_crit ∼ 2 × 10⁵ cm⁻³ and ∼8 × 10⁵ cm⁻³, respectively, as estimated from the Flower (1999) collision rates for HCO⁺ and line parameters (Einstein A-coefficients, frequencies, statistical weights) from the JPL Database⁴ (Pickett et al. 1998), originally measured in Anderson et al. (1977) and Sastry et al. (1981). We note that Hugo et al. (2009) provide inelastic collisional rates for the excitation of H₂D⁺; however, to simplify the calculations we assume the H₂D⁺ level populations are also in LTE. This assumption is appropriate as the H₂D⁺ emission originates from dense gas exceeding $n_{\rm {H_2}}>10^{8}$ cm⁻³, much larger than the critical density for o-H₂D⁺ (1₁₀–1₁₁), n_crit ∼ 10⁵ cm⁻³ (Hugo et al. 2009).

For the carbon and oxygen isotopes, we adopt ¹⁶O/¹⁸O = 500 (Kahane et al. 1992; Prantzos et al. 1996) and ¹²C/¹³C = 60 (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 D/H = 2 × 10⁻⁵ (see Section 4.1). Leiden LAMDA formatted H₂D⁺ and N₂D⁺ input files were compiled from the CDMS database⁵ (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), and Yonezu et al. (2009) for H₂D⁺. 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 (τ₀ < 1), while N₂H⁺ (4–3) and (3–2) reach τ₀ ∼ 5 just beyond the CO snowline for the W98 and M02 CR models. Observations of H¹³CO⁺ (3–2) and HCO⁺ (3–2) are reported in Öberg et al. (2011b) toward DM Tau. The disk integrated flux ratio of the (3–2) transition for the ¹²C/¹³C isotopologues is ∼7.6, and assuming a respective isotope ratio of 60 (Keene et al. 1998), corresponds to an H¹³CO⁺ (3–2) optical depth of τ = 0.14. 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 N₂D⁺ and H₂D⁺ lines are thin throughout the disk for all models. C¹⁸O is thick, with τ₀ ≳ 1 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 Section 6.3).

Zoom In Zoom Out Reset image size

The HCO⁺ isotopologues show a peak in their opacity offset from the radial center of the disk, where the column is highest (see Figure 4). The relative height between the central peak and the outer broad 100–300 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 (3–2) transition where depopulation becomes important toward hotter inner disk gas, whereas H¹³CO⁺ (5–4) appears more centrally peaked.

The N₂H⁺ 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 R ∼ 300 AU is fueled by primarily CR ionization. The H₂D⁺ optical depth also peaks in this same outer region in the models where CR ionization is active. The N₂D⁺ has a somewhat complicated τ₀ profile, which is reflected in the N₂D⁺ column densities (Figure 4). The high CR ionization models sustain N₂D⁺ co-spatial with N₂H⁺ and, when combined with warm temperatures close to the upper-state excitation temperatures (E_u = 22.2 K and 37.0 K for J = 3 and J = 4 rotational levels, respectively), leads to N₂D⁺ emission inside of R < 100 AU.

From the LIME emission models described in Section 5.1, we simulate full-science ALMA observations for the ground-based accessible transitions of HC¹⁸O⁺, H¹³CO⁺, N₂H⁺, N₂D⁺, and H₂D⁺. 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 N₂D⁺ (5–4) and HD$_2^+$ (1₁₀–1₀₁), 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.

Table 4. CASA Simulation Parameters

Emission	12 m Configurationa	Sensitivityb
Line	(12 m + 7 m Beam)	(mJy)
HC18O+ (3–2)	03 (128)	1.0
HC18O+ (4–3)	03 (096)	1.4
H13CO+ (3–2)	10 (033)	0.97
H13CO+ (4–3)	10 (025)	1.4
H13CO+ (5–4)	10 (020)	3.6
N2H+ (3–2)	10 (031)	1.3
N2H+ (4–3)	10 (024)	2.6
N2D+ (3–2)	01 (155)	1.3
N2D+ (4–3)	01 (116)	1.3
o-H2D+ (110–111)	01 (096)	2.5

Notes. ^ahttp://casaguides.nrao.edu/index.php?title=Antenna_List Beam averaged over major and minor axes from the combined 12 m and ACA observations. ^b6.1h in 0.2 km s⁻¹ bins from the ALMA Sensitivity Calculator.

Download table as:  ASCIITypeset image

5.2.1. Observational Parameters

The simulations presented here reflect 6h of total on-sky time with one of the 12 m array configurations and an ACA 7 m observation. While the choice of 6h is somewhat arbitrary, it represents a relatively "deep" observation, and ideally one would target more than one emission line toward 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 H₂D⁺, 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 N₂D⁺ 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 H¹³CO⁺, and N₂H⁺, 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 R = 400 AU disk at d = 100 pc on the sky is 8'', and for the higher frequency (ν > 300 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. Rather than showing images or channel maps for every emission line/model considered both in model emission (Section 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.

Zoom In Zoom Out Reset image size

5.2.2. ALMA Sensitivity Calculation

To calculate the uncertainty on the ALMA measurements, we use the uncertainties from the ALMA 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 t_int = 1000 s for the simulations), we scale the sensitivity calculator's uncertainty σ₀ by $\sigma _i=\sigma _0 \sqrt{N_A(N_A-1)} \sqrt{N_{\rm {time}}}\sqrt{N_{\rm {pol}}}$, where N_A is the number of antennas (50 for the 12 m array and 12 for the 7 m array), N_time is the number of time samplings during the course of the simulated observation's integration time (e.g., t_obs = 22, 000 s =6.1h on the 12 m array such that N_time = 22) and N_pol = 2 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, $\sigma _A=\sigma _i/\sqrt{N_{\rm bin}}$, where N_bin is typically ≳ 100.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 Section 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. 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 1'' = 100 AU. To create the profiles, we average 10 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.

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 (3–2) transition, less prominently seen in the (4–3). The reason that the (3–2) emission is brightest away from the column density peak is due to excitation of the line, which has an upper state energy of E_u ∼ 25 K. The abundant HCO⁺ gas inside of R ≲ 50 AU has temperatures typically exceeding 100 K. At this point, higher J-states of H¹³CO⁺ and HC¹⁸O⁺ 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 H¹³CO⁺, (4–3) and (5–4), there is an inner (θ < 05) 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 HC¹⁸O⁺ abundances, but the lower spatial resolution simulations adopted for HC¹⁸O⁺ (see Table 4) do not resolve the inner 50 AU. In summary, the HCO⁺ emission sensitively traces (1) stellar X-ray processes and (2) 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.

N₂H⁺ is expected to arise from cold gas where its destroyer, CO, is frozen out at temperatures below T < 20 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 N₂H⁺ 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 N₂H⁺ 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 N₂H⁺ 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 ∼seven in brightness (see the zoom-in in Figure 9), but detailed modeling may be required to estimate the potential X-ray contribution to the N₂H⁺ column density. An additional caveat in using N₂H⁺ as an ionization tracer is that its emission can be partially optically thick (if not completely, τ ∼ 5; Figure 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, N₂H⁺ 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 ζ_CR ∼ 10⁻¹⁹ 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 3σ and 6σ level for the (4–3) and (3–2) transitions of N₂H⁺, respectively. Given the N₂H⁺abundance's concurrent chemical dependence on temperature, a warmer disk where the N₂H⁺ abundance peaks further from the X-ray bright star may allow N₂H⁺ 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.

Zoom In Zoom Out Reset image size

N₂D⁺ shows interesting emission behavior such that the M02 models show less N₂D⁺ emission than the more weakly CR-ionized W98 models (Figure 8). This behavior is reflected in the column densities (Figure 4) and the abundances (Figure 3). More specifically, the M02 N₂D⁺ models show a deficit in abundance relative to the W98 models centered at R ∼ 200 AU. The inner N₂D⁺ gas is co-spatial with the high N₂H⁺ (fueled by X-ray flux) and with the boundary of the H₂D⁺ abundant region (T < 50 K), creating an N₂D⁺ layer. This layer is reflected in the W98 models as well. At intermediate distance radii R ∼ 150–300 AU, the M02 models show an N₂D⁺ 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 N₂H⁺ but instead sequestered into nitrogen bearing ices like NH₃, HCN and NO, and therefore is not available as N₂ in the gas to reform N₂H⁺ and N₂D⁺. This process is an analog 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 ∼10 brighter X-ray luminosity would have a similar effect (see Section 6.1).

The true utility of N₂D⁺ 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 (3–2) and (4–3) and in the column density plane. N₂D⁺ 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 ∼20, more than double that of N₂H⁺, allowing these CR models to be more easily disentangled than with N₂H⁺.

Finally, H₂D⁺ 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 ζ_CR ≲ 10⁻¹⁹ s⁻¹, the H₂D⁺ (1₁₀–1₁₁) 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 H₂D⁺ column toward 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 H₂D⁺ as a tracer. The utility of H₂D⁺ is highlighted when used in tandem with N₂D⁺. As mentioned above, N₂D⁺ decreases in brightness for both high ionization rates (when the precursors of N₂D⁺ are chemically destroyed) and low ionization rates (when N₂D⁺ is not produced). Observations of H₂D⁺ in conjunction with N₂D⁺ 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 N₂H⁺ and N₂D⁺ 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.

The results presented in Figures 7–9 consider a single X-ray luminosity, L_XR = 10^29.5 erg s⁻¹. To understand the sensitivity of the lines to X-ray ionization, we have computed an additional chemical model for a 10-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.

Zoom In Zoom Out Reset image size

The high X-ray ionization rate changes the HCO⁺ column density most significantly in the inner R < 100 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 three to four. This trend is explained by the balance between ionization and recombination, such that the steady state abundance of ions (or electrons) is proportional to $\propto \sqrt{\zeta /(\alpha n_{\rm {H_2}})}$, where α is the recombination rate and $n_{\rm {H_2}}$ is the number density of H₂. With everything else constant, a factor of 10 increase in ionization results in a factor of $\sqrt{10}\approx 3$ 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 N₂H⁺, N₂D⁺, and H₂D⁺ 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.

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 the Appendix), the ensemble of SLRs that dominate the ionization, namely ²⁶Al and ⁶⁰Fe, have an effective net half-life of approximately t_half ∼ 1.2 Myr (see Section 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 (>3 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), 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 t = 3 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 t = 0 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 N₂D⁺, H₂D⁺, and HD$_2^+$ 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.

Zoom In Zoom Out Reset image size

The results of our study are for a particular disk model (Section 2) and disk mass (M_g = 0.04 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_g = 0.02 M_☉) and double mass (M_g = 0.08 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 ∼two in density would change the disk opacity and would result in a correspondingly cooler (warmer) disk by ≲ 10% 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 Section 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 (∼1.2 times more ionization; Cleeves et al. 2013b).

In Figure 12, we show the column of the indicated molecular ion, N_i, normalized to the column of model CO, N_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 C¹⁸O, however we emphasize that even C¹⁸O has a minimum optical depth of τ ∼ 1 at line center in this particular model (see Figure 5). To determine the CO column observationally, one must consider either (1) more rarified CO isotopes than C¹⁸O or (2) 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.

Zoom In Zoom Out Reset image size

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 Section 6.1, where the ion abundance is proportional to $\sqrt{ \zeta / (\alpha n_{\rm {H_2}})}$. In terms of column densities, the ratio of the ion-species to CO is then approximately given by $N_i/N_{\rm {CO}}\sim \chi _i N_{\rm {H_2}}/\chi _{\rm {CO}}N_{\rm {H_2}} = \chi _i/\chi _{\rm {CO}} = \chi _{\rm {CO}}^{-1}\sqrt{ \zeta / (\alpha n_{\rm {H_2}})}$. For a more massive disk, ζ_XR decreases due to increased gas opacity, and $n_{\rm H_2}$ increases as the mass increases. Thus the quantity $\zeta / n_{\rm {H_2}}$ 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, ζ_XR, leaving $\sqrt{n_{\rm H_2}^{-1}}$ as the dominant term. M02 shows substantially more variation than can be attributed to changes in $n_{\rm H_2}$ or in X-ray flux, otherwise the variation would be seen in W98 as well. This behavior can be understood by looking toward the denominator of the abundance ratio: χ_CO. 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 H₂D⁺ and HD$_2^+$, 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 N_i/N_CO depends sensitively on both ζ_XR and $n_{\rm H_2}$. 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.

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 R ∼ 1 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_eff = 6000 K, where all other structural parameters are held fixed. At R = 15 AU, the disk surface (z/r = 0.3) temperature increases substantially from T_d ∼ 100 K to ∼140 K. The warm molecular layer (z/r = 0.1) increases from T_d = 40 K to 55 K, and the midplane increases by ΔT_d = 10 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. 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 Sections 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 R > 250 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.

Zoom In Zoom Out Reset image size

For all CR-ionization models, the column density of N₂H⁺ shows a decrease in the inner R < 100 AU and a surplus beyond R > 150 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 N₂H⁺ 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 N₂H⁺. The N₂D⁺ column density shows similar morphological changes as N₂H⁺ 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 H₂D⁺ relative to H$_3^+$. The H₂D⁺ 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 R > 50 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 (R > 50 AU). At R ∼ 175 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 ∼60% 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. N₂H⁺ and N₂D⁺ remain sensitive observational tracers of the CR ionization rate at both high and low CR fluxes beyond R > 100 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 N₂H⁺ 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 N₂H⁺ and the same behavior holds true for N₂D⁺. 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 C¹⁸O 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.

The CR contribution to disk ionization in both magnitude and scope is the least observationally constrained parameter in disks. A CR ionization rate of ζ_CR ≳ 3 × 10⁻¹⁷ s⁻¹, consistent with dense ISM values, was ruled out by Chapillon et al. (2011) using observations of H₂D⁺. 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 analog to the solar system's heliosphere. Within this paradigm we have ignored radial variations in the CR rate that may include (1) a gradual increase (1%/AU) in the CR rate with radius where modulation is weaker further out in the disk (Paper I) and (2) an edge to the region of CR-modulation, i.e., a heliopause. As compared to the negative ionization gradient expected for SLR ionization, (1) would create a positive ionization gradient, though the two effects are similar (∼one 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: (1) 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 (2) 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 20%–30% 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 toward the midplane (see Equation (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 ζ_CR ≳ 10⁻¹⁷ 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 (Sections 6.1 and 6.3, respectively). Alternatively, if neither of these effects are important, wind modulation nor magnetic effects, then the ionization should be ζ_CR ≳ 10⁻¹⁷ s⁻¹ in the outer disk and the more weakly emissive, difficult species to observe, like H₂D⁺, can help provide constraints on midplane ionization.