THE IONIZATION FRACTION IN THE DM Tau PROTOPLANETARY DISK

DM Tau was observed with SMA on top of Mauna Kea, Hawaii, on 2010 October 7. Six of the array antennas were available for these observations, arranged in a compact configuration spanning baselines of 16–77 m. Using the dual receiver mode, we simultaneously targeted the ortho-H₂D⁺ 1_{1, 0}–1_{1, 1} line at 372.4213 GHz and the N₂H⁺ J = 4–3 at 372.6725 GHz with the "high-frequency" receiver, and the 345.796 GHz CO J = 3–2 line with the "low-frequency" receiver. Table 1 summarizes the spectral line setup. At these frequencies, the synthesized beam sizes were 22 × 17 and 24 × 20, respectively.

The observing conditions were generally very good, with atmospheric opacity τ_225 GHz ∼ 0.04 measured at the nearby CSO and stable atmospheric phase. This opacity corresponds to a zenith transmission of ∼50% at the H₂D⁺ line frequency. Short observations of the two calibrators 3C120 and J0449+113 were interleaved with the observations of DM Tau for gain calibration. The bandpass response was calibrated with observations of Uranus and the available bright quasars 3C273, 3C454.3, and 3C84. Observations of Uranus and Callisto provided the absolute scale for flux densities, via comparison to theoretical models for their emission. The two calibrators agree within the expected uncertainties, but all reported data were calibrated with Callisto. The systematic uncertainty in the absolute flux scale is ∼10%. The data were edited and calibrated with the IDL-based MIR software package.² Continuum and spectral line images were generated and CLEANed using MIRIAD. To check that the flux calibration is accurate, we compared the continuum fluxes of 0.22 Jy at 0.87 mm and 0.30 Jy at 0.81 mm with previously published values (Andrews & Williams 2005, 2007) and found that they agree within the reported uncertainties.

The H¹³CO⁺ J = 3–2 line at 260.255 GHz was observed toward DM Tau on 2011 January 25 with the IRAM 30 m telescope. The observations were carried out with the EMIR 330 GHz receiver, with a beam size of ∼95 (FWHM). The receiver was connected to a unit of the autocorrelator with a spectral resolution of 320 kHz and a bandwidth of 240 MHz, equivalent to an unsmoothed velocity resolution of ∼0.4 km s⁻¹. Typical system temperatures were 300–400 K. The observations were carried out using wobbler switching with a 100'' throw. Pointing was checked every ∼2 hr on J0430+052 and J0316+413 with a typical accuracy of <2''. The main-beam brightness temperature was calculated from the antenna temperatures in pako using reported main beam and forward efficiencies (B_eff and F_eff) of 88% and 53%, respectively. The data were reduced with the CLASS program, part of the GILDAS software package.³ Linear baselines were determined from velocity ranges without emission features and then subtracted from the spectra.

To estimate the N₂H⁺ rotational temperature, we assume LTE and co-spatial emission from the N₂H⁺ 3–2 and 4–3 lines. This is a reasonable approximation because chemical models predict that N₂H⁺ should be abundant in a narrow temperature regime between where CO begins to freeze out and N₂ freezeout is complete (Bergin et al. 2006; Dutrey et al. 2007). Assuming optically thin emission, unity filling factor within the CO 3–2 disk, and molecules thermalized at a single rotational temperature,

$$\begin{equation} N = \frac{1.67\times 10^{14}}{\nu \mu ^2S} Q(T_{\rm rot})e^{E_{\rm u}/T_{{\rm rot}}}\int T_{\rm b}dv \end{equation} \tag{ 1 }$$

from Thi et al. (2004), where N is the total number of molecules, Q(T_rot) is the temperature-dependent partition function, E_u is the energy of the upper level in K, and T_rot is the rotational temperature. We obtain the partition function Q(T_rot) by second-order interpolation of the listed values in CDMS,⁴ μ²S is gathered directly, and E_u is calculated from the listed E_low values: E_u = E_low + 0.0334 × ν, where ν is the line frequency in GHz. Then an upper limit on the N₂H⁺ temperature can be calculated from

$$\begin{eqnarray} &&T_{{\rm rot}} = \frac{E_u^{4-3}-E_u^{3-2}}{ln\left(\frac{\nu _{4-3}\mu _{4-3}^2S_{4-3}}{\nu _{3-2}\mu _{3-2}^2S_{3-2}}\frac{\int T_{\rm b}^{3-2}dv}{\int T_{\rm b}^{4-3}dv}\right)}. \end{eqnarray} \tag{ 2 }$$

Using the integrated intensity of the J = 3–2 line and the 2σ upper limit on the J = 4–3 line, we obtain a 2σ upper limit on the N₂H⁺ rotational temperature of <26 K. This is consistent with the predictions that N₂H⁺ should become abundant below the CO freezeout temperature of 20 K, but does not provide any stronger constraints. It is also consistent with the N₂H⁺ J = 1–0 line reported by Dutrey et al. (2007). We estimate the line flux from their Figure 1 to be <0.1 Jy km s⁻¹. Using this limit and the integrated intensity of the J = 3–2 line results in a lower limit on the N₂H⁺ rotational temperature of 11 K. More sensitive observations of the N₂H⁺ lines would be valuable to constrain this excitation temperature further. In the absence of more detailed information, we calculate LTE column densities of midplane ions for a range of temperatures, from 10 to 20 K. The 20 K boundary is set by laboratory experiments on CO and N₂ freezeout—above this temperature a majority of CO and N₂ are maintained in the gas phase quickly destroying any formed H⁺₃ (Bisschop et al. 2006). The CO snowline was recently reported to correspond to a freezeout temperature of 19 K, consistent with predictions (Qi et al. 2011).

In general, column densities are calculated most accurately by constructing a self-consistent chemical-physical disk model and applying a radiative transfer code to predict line emission profiles that can be compared directly to observations. Without a detection of the H₂D⁺ line, however, it is difficult to constrain where the H₂D⁺ emission arises in the disk. We therefore opt for a classical LTE approach (e.g., van Dishoeck et al. 2003), calculated for the range of reasonable disk midplane temperatures (see Section 4.1) to put limits on the disk-averaged N(o-H₂D⁺).

A significant complication for the LTE calculation is that ortho (o) and para (p) H₂D⁺ can behave as separate species at low temperatures (Flower et al. 2004; Sipilä et al. 2010) resulting in a non-thermal partitioning between the o and p states. To calculate the partition function for o-H₂D⁺, we use a two-level approximation, which assumes no transitions between the o- and p-H₂D⁺, and that transitions between the ground state and first excited state dominate within each "species" of H₂D⁺. This approximation is expected to be valid at temperatures of 10–20 K, since E_u is 104 K for the first excited state of o-H₂D⁺. The o-H₂D⁺ partition function is calculated from

$$\begin{equation} Q=g_u\times e^{-E_u/T_{\rm ex}}+g_l\times e^{-E_l/T_{\rm ex}}, \end{equation} \tag{ 3 }$$

where the degeneracies, g_u and g_l, are equal to 3.

Assuming LTE (Equation (1)), excitation temperatures in the range 10–20 K and the partition function appropriate for these temperatures, N(o-H₂D⁺) limits are calculated from the observed upper limit on the H₂D⁺ J = 1_{1, 0}–1_{1, 1} line flux. As shown in Figure 3(a), the limits on the disk average N(o-H₂D⁺) range from 1.3 to 2.7 × 10¹² cm⁻² (with the highest value corresponding to the lowest temperature).

The limit on the total N(H₂D⁺), i.e., o + p, depends on the o-H₂D⁺ column density and the temperature-dependent o/p ratio. This ratio has been modeled by Sipilä et al. (2010) from 4 to 20 K for a molecular hydrogen density of 10⁶ cm⁻³ and an interstellar grain distribution. The derived o/p ratio is 0.5–0.1 for temperatures of 10–20 K, with the lowest ratio at 13–17 K. While the assumed grain size distribution and density are more appropriate for dense clouds than disk midplanes, the grain size distribution mainly affects the temperature profile, and Sipilä et al. (2010) find that the H₂D⁺ o/p is nearly constant with increasing density. Thus, we expect the values calculated with these assumptions should be valid for disk midplane conditions.

Figure 3(b) shows the result of applying the literature o/p ratios to the previously calculated o-H₂D⁺ column density limits to obtain limits on N(H₂D⁺) as a function of temperature. These limits range from 4 to 21 × 10¹² cm⁻². The presence of a peak value at 13 K is a direct consequence of the fact that the limit on N(o-H₂D⁺) decreases with temperature, while the o/p ratio has a minimum at 13–17 K.

The limit on the total column density of ions in the midplane, N(∑ H_{3 − x}D⁺_x), depends on N(H₂D⁺) and the ratio of N(H₂D⁺) to N(∑ H_{3 − x}D⁺_x). Like the o/p ratio, this ratio is expected to vary with density and temperature, as well as other environmental properties such as grain size (Flower et al. 2004; Ceccarelli & Dominik 2005; Sipilä et al. 2010). This ratio will also depend on whether or not CO and N₂ are frozen out completely, or if low abundances of these species can be maintained in the gas phase through non-thermal desorption (Asensio Ramos et al. 2007). Caselli et al. (2008) modeled the ratios of all deuterated H⁺₃ isotopologues for a molecular hydrogen density 10⁵ cm⁻³ and interstellar grains as a function of temperature. As mentioned previously, the disk midplane density is likely higher than assumed in these models. Sipilä et al. (2010) found an increasing importance of D⁺₃ with increasing density at 10 K, but it is unclear whether this effect is important also at higher temperatures. Another potential unknown is the abundance of H⁺, though Ceccarelli & Dominik (2005) found that its contribution is negligible at the high densities of the disk midplane.

With these caveats in mind, we use the ratios derived by Caselli et al. (2008) with maximum depletion, i.e., minimum amounts of CO and N₂ in the gas phase to calculate upper limits on N(∑ H_{3 − x}D⁺_x) in the disk midplane as a function of assumed temperature. Figure 3(c) shows the results, with limits that range from 1.5 to 8.5 × 10¹³ cm⁻², with the peak value at 13 K. The modification of the shape compared to the H₂D⁺ curve is due to the fact that the N(H₂D⁺)/N(∑ H_{3 − x}D⁺_x) ratio peaks near 13 K and then falls off quickly at higher temperatures.

To derive limits on the midplane ion abundances, we adopt the model of the DM Tau disk density and temperature structure described by Andrews et al. (2011) based on parametric fitting of the spectral energy distribution and resolved observations of millimeter-wave dust emission. This model does not include any chemistry, but rather provides an observationally constrained physical structure, which can be used to derive the gas masses in different temperature layers in the disk. This is needed to convert an ion column density into an ion abundance. We apply a standard interstellar gas-to-dust ratio of 100 to convert the dust densities in this model to gas densities. To determine a midplane mass, we define the disk midplane as all material ≲16 K, the temperature for complete N₂ freezeout in most astrophysical environments based on laboratory measurements (Öberg et al. 2005; Bisschop et al. 2006), i.e., the region where it is most likely that H⁺₃ and its isotopologues become abundant. In this disk layer, the density-weighted temperature is 11 K, and the disk-averaged gas column density N$_{\rm H_2}$ is 15 × 10²² cm⁻². Using this column density, and the upper limit on the ion column density in the midplane of 4 × 10¹³ cm⁻² from Table 3 and Figure 4, x_i < 2.7 × 10⁻¹⁰ n_H. Alternatively, if the onset of CO freezeout at 20 K defines the region where H₂D⁺ is abundant rather than N₂ freezeout at 16 K, then the density-averaged ion temperature is 13 K, and x_i < 4.5 × 10⁻¹⁰ n_H. Table 3 summarizes the midplane (<16 K) ionization results and also shows how the ionization fractions compare if H₂D⁺ is only assumed to be present in the midplane or in the entire T <20 K region. It is important to keep in mind that all of these calculations are based on model-dependent values, and better limits on the midplane ionization abundance require not only more sensitive o-H₂D⁺ observations, but also direct constraints on the o/p ratio and the relative abundances of different isotopologues. In particular, the latter estimate may be off by an order of magnitude from different model predictions of the importance of D⁺₃ compared to H₂D⁺.

Zoom In Zoom Out Reset image size

It is interesting to compare the estimates of the ionization fraction in different layers of the disk, since this may impact, e.g., MRI activation and accretion flow. The mass of each layer, defined by a temperature interval, can be estimated from the same DM Tau physical model as was used in Section 4.5. Between 20 K, where CO begins to freeze out, and the CO photodissociation region, HCO⁺ is predicted to be the most abundant ion (e.g., Willacy 2007). Below 20 K, as CO starts to freeze out, N₂H⁺ becomes more abundant, since reactions with CO are the main destruction pathway of N₂H⁺. Below 16 K, CO and N₂ are completely frozen out and H⁺₃ and its isotopologues become the main charge carriers.

The DCO⁺ abundance peaks between the region where HCO⁺ dominates and the midplane (Willacy 2007; Aikawa & Nomura 2006) and therefore is assumed to coexist with N₂H⁺ in the intermediate layer. In this layer, DCO⁺/HCO⁺ exceeds unity and the HCO⁺ contribution to the charge balance should be small. This model result is supported by comparing DCO⁺ line intensities: Guilloteau et al. (2006) detected DCO⁺ 2–1 at 0.29 ± 0.03 Jy km s⁻¹ and the intensity ratio of the 3–2/2–1 lines implies a rotational temperature of 11–23 K, assuming LTE.

For the midplane and the low-temperature intermediate layer, the emission regions of the ions are expected to be well confined, and the assumption of LTE at a single rotational temperature should be a reasonable first approximation. For the warm molecular layer where HCO⁺ dominates, the assumption of LTE is not a priori a good approximation since this layer spans a temperature range from 20 to 100 K or higher. However, the derived HCO⁺ column density is not very sensitive within this range; it varies by less than a factor of three when using Equation (1) over the full range of plausible temperatures and the LTE approach retains utility. As discussed below, the outcome of this LTE calculation compares well with a radiative transfer calculation of HCO⁺ emission in the DM Tau disk (Piétu et al. 2007).

We estimate the disk-averaged N(HCO⁺), N(N₂H⁺), and N(DCO⁺) using Equation (1), the observed intensities from Table 2, and the mass-weighted temperatures from the DM Tau structural model of 38 K for the CO-dominated layer (>20 K) and 18 K for the intermediate layer where N(N₂H⁺) and N(DCO⁺) are most abundant (16–20 K). Table 3 lists the results and Figure 5 shows the resulting column densities together with the H⁺₃ limit (LTE at 11 K), within their respective temperature regimes (<16 K, 16–20 K, and >20 K), corresponding to the different disk layers. Note that N(HCO⁺) is calculated from H¹³CO⁺ assuming a standard ¹²C/¹³C ratio of 70, since the high opacity of the main isotopologue underestimates N(HCO⁺) by a factor of seven.

Zoom In Zoom Out Reset image size

The calculation results in a disk-averaged DCO⁺/HCO⁺ ratio of ∼0.05. This is comparable to a previous estimate of 0.035 toward the disk around TW Hya (Thi et al. 2004). Both exceed the cosmic D/H ratio of ∼10⁻⁵ by three orders of magnitude, demonstrating an efficient deuterium fractionation in these disks. The derived ratio is an order of magnitude above the ratio found by Guilloteau et al. (2006), who used resolved HCO⁺ 1–0 and single-dish DCO⁺ 2–1 emission and a radiative transfer calculation rather than LTE. The derived disk-averaged N(HCO⁺) are almost identical over the inner 400 AU of the disk: 7 × 10¹² cm⁻² from their tabulated column density at 300 AU and power-law index, compared to 9 × 10¹² cm⁻² calculated here assuming LTE at 38 K. The different ratios must thus be due to different assumptions on the DCO⁺ emission conditions. The HCO⁺ column density compares well with model predictions of the outer disk (100–400 AU). It is for example consistent with the disk chemistry model by Willacy (2007) within a factor of three.

Dutrey et al. (2007) observed N₂H⁺ 1–0 with the Plateau de Bure interferometer and derived N₂H⁺/HCO⁺ of ∼0.02–0.03. This ratio agrees well with the results in Table 3, which implies a disk-averaged ratio of 0.03.

The disk-averaged N(H₂) within each temperature regime are 15, 2.4, and 2.6 × 10²² cm⁻², and these values are used with the derived ion column densities to calculate disk-averaged x_i in each layer. Figure 5 shows that the HCO⁺ abundance and the midplane ion limits are comparable at ∼3–4 × 10⁻¹⁰ n_H, while the abundance of DCO⁺ + N₂H⁺ is an order of magnitude lower at 3 × 10⁻¹¹ n_H.

Overall, the derived ion abundances suggest a decreasing x_i with depth in the disk, which might be expected from attenuation of ionizing radiation. The assumption that DCO⁺ and N₂H⁺ alone are representative for the degree of ionization in the intermediate layer may be oversimplified, however. The 16 < T < 20 layer is a transition region where some HCO⁺ may still be present and H⁺₃ may already have started to become important. A combination of better limits on the midplane ions as well as more detailed modeling of the intermediate layer with DCO⁺ and N₂H⁺ will be needed to confirm the depth-dependent trend on x_i.

The main sources of ionizing radiation in disks are ultraviolet (UV) photons, cosmic rays, radioactive decay, and X-rays (e.g., Willacy 2007), and the ionization level therefore depends on a combination of the incident fluxes and the attenuation of external photons with dust and gas column density. Willacy (2007) includes both the interstellar radiation field and a stellar UV field that is 500 times stronger at 100 AU. Because of efficient absorption by dust and H₂, the UV ionization is mainly important in the disk atmosphere and upper molecular layer. Cosmic rays are assumed to penetrate unperturbed to the disk midplane resulting in an ionization rate of 1.3 × 10⁻¹⁷ s⁻¹. Radioactive decay also contributes to the midplane ionization, at a rate of 0.6 × 10⁻¹⁷ s⁻¹. The effect of X-rays is not included in the model and the model may therefore underestimate x_i in lower disk layers.

In this model, the total ion column density (H⁺₃, H₂D⁺, D₂H⁺, D⁺₃, HCO⁺, DCO⁺, N₂H⁺) at 250 AU is 2–3 × 10¹³ cm⁻², corresponding to an average x_i of 2–3 × 10⁻¹⁰ n_H. This value is consistent with the observed abundances and limits. In the related disk model of Asensio Ramos et al. (2007), the ion abundance depends on the assumed cosmic ray flux, as expected.

In disks, including the model from Willacy (2007), the temperature increases both inward and upward in the disk, resulting in a vertically layered structure outside of an inner truncation radius. In the model, HCO⁺ extends through the disk, while DCO⁺ is only abundant outside of 50 AU and H₂D⁺ outside of 100 AU. As mentioned above, beyond 100 AU the different ion column densities are basically constant, justifying the comparison between disk-averaged observations and model predictions at 250 AU.

Because all UV photons and probably some of the X-rays and cosmic rays are absorbed in the outer disk layers, the ionization fraction should decrease between the outer disk layers and the disk midplane. At 250 AU, in the disk chemistry models (Figure 3 in Willacy 2007) the ion abundance in the upper molecular layer, where CO is in the gas phase and HCO⁺ is the dominant ion, is a factor of five higher than in the midplane (∼5 × 10⁻¹⁰ compared to ∼1 × 10⁻¹⁰ n_H). In the intermediate region, where DCO⁺ and N₂H⁺ abundances peak, the ion abundance is about a factor of two to three below the CO layer ion abundance. This is qualitatively consistent with the results of the analysis of the DM Tau observations (see Figure 4).

HCO⁺ and DCO⁺. The HCO⁺ column density predicted by Willacy (2007) is ∼3 × 10¹² cm⁻² at 250 AU and is found to not depend on the different desorption mechanisms considered. In contrast, the DCO⁺ column density predictions vary by an order of magnitude, resulting in a DCO⁺/HCO⁺ ratio of 0.1–1. The higher DCO⁺/HCO⁺ ratios appear with desorption through cosmic ray heating, which maintains CO in the gas phase in the disk midplane where the temperatures are low enough for efficient deuterium fractionation and enhanced DCO⁺ production. The observations indicate that the disk-averaged DCO⁺/HCO⁺ abundance ratio is below 0.1, which suggests that CO is confined to the warm >20 K molecular layer. The observational results are consistent with including UV photodesorption, which is efficient only in the warm molecular layer.

N₂H⁺. The N₂H⁺/HCO⁺ abundance ratio is predicted to be in the range from 2 to 9 × 10⁻³, but the models with the highest values also predict DCO⁺/HCO⁺ ratios of unity and can be excluded. Compared to the value of 0.03 derived from the observations, the models are an order of magnitude too low. The reason for this discrepancy is unknown, demonstrating that even the simplest nitrogen-based chemistry remains poorly constrained.

H₂D⁺. N(H₂D⁺) is predicted to be 10¹² cm⁻², with midplane abundance ∼10⁻¹¹ n_H. These values are well below the observational limits, and the reason is that D⁺₃ is the dominant carrier of charge in the midplane in these models, 20–50 times more abundant than H₂D⁺. In the disk model by Aikawa & Nomura (2006), the H₂D⁺ column density is predicted to be significantly higher, 10¹³ cm⁻², likely because multiple deuteration is not considered. Both models are thus consistent with the current limits, and better constraints on the chemical makeup of the midplane will require not only more sensitive observations of H₂D⁺, but also robust o/p determinations and additional constraints on the abundance ratios of the isotopolgues.