ApJ accepted: March 18, 2008
Preprint typeset using LATEX style emulateapj v. 03/07/07
RESOLVING THE CHEMISTRY IN THE DISK OF TW HYDRAE
I. DEUTERATED SPECIES
Chunhua Qi1 , David J. Wilner1 , Yuri Aikawa2 , Geoffrey A. Blake3 , Michiel R. Hogerheijde4
arXiv:0803.2753v1 [astro-ph] 19 Mar 2008
ApJ accepted: March 18, 2008
ABSTRACT
We present Submillimeter Array (SMA) observations of several deuterated species in the disk around
the classical T Tauri star TW Hydrae at arcsecond scales, including detections of the DCN J=3–2
and DCO+ J=3–2 lines, and upper limits to the HDO 31,2 –22,1 , ortho-H2 D+ 11,0 –11,1 and para-D2 H+
11,0 –10,1 transitions. We also present observations of the HCN J=3–2, HCO+ J=3–2 and H13 CO+
J=4–3 lines for comparison with their deuterated isotopologues. We constrain the radial and vertical
distributions of various species in the disk by fitting the data using a model where the molecular
emission from an irradiated accretion disk is sampled with a 2D Monte Carlo radiative transfer code.
We find that the distribution of DCO+ differs markedly from that of HCO+ . The D/H ratios inferred
change by at least one order of magnitude (0.01 to 0.1) for radii <30 AU to ≥70 AU and there is a
rapid falloff of the abundance of DCO+ at radii larger than 90 AU. Using a simple analytical chemical
model, we constrain the degree of ionization, x(e-)=n(e-)/n(H2 ), to be ∼ 10−7 in the disk layer(s)
where these molecules are present. Provided the distribution of DCN follows that of HCN, the ratio of
DCN to HCN is determined to be 1.7±0.5 × 10−2 ; however, this ratio is very sensitive to the poorly
constrained vertical distribution of HCN. The resolved radial distribution of DCO+ indicates that in
situ deuterium fractionation remains active within the TW Hydrae disk and must be considered in
the molecular evolution of circumstellar accretion disks.
Subject headings: circumstellar matter —comets: general —ISM: molecules —planetary systems: protoplanetary disks —stars: individual (TW Hydrae) —stars: pre-main-sequence
1. INTRODUCTION
Millimeter-wave interferometers have imaged the gas
and dust surrounding over a dozen T Tauri and Herbig
Ae stars (see Dutrey et al. 2007 for a review). These
studies demonstrate the potential to dramatically improve our understanding of disk physical and chemical
structure, providing insights that will ultimately enable
a more comprehensive understanding of star and planet
formation. As analogs to the Solar Nebula, these circumstellar disks offer a unique opportunity to study the conditions during the planet formation process, especially
the complex chemical evolution that must occur. In the
outer parts of the disk directly accessible to millimeterwave interferometry, observations of deuterated species
are particularly important because they can constrain
the origin of primitive solar system bodies such as comets
and other icy planetesimals.
Deuterated molecule chemistry is sensitive to the temperature history of interstellar and circumstellar gas, as
well as to density-sensitive processes such as molecular
depletion in the cold (< 20 K) and dense (> 106 cm−3 )
disk midplane. The similarity of molecular D/H ratios
between comets and high mass hot cores has been used
to argue for an “interstellar” origin of cometary matter, but ambiguity remains and this argument is not se1 Harvard–Smithsonian Center for Astrophysics, 60 Garden
Street, MS 42, Cambridge, MA 02138, USA; cqi@cfa.harvard.edu,
dwilner@cfa.harvard.edu.
2 Department of Earth and Planetary Sciences, Kobe University,
Kobe 657-8501, Japan; aikawa@kobe-u.ac.jp.
3 Divisions of Geological & Planetary Sciences and Chemistry
& Chemical Engineering, California Institute of Technology, MS
150–21, Pasadena, CA 91125, USA; gab@gps.caltech.edu.
4 Leiden Observatory, Leiden University, PO Box 9513, 2300 RA,
Leiden, The Netherlands; michiel@strw.leidenuniv.nl.
cure (Bergin et al. 2007). Aikawa & Herbst (1999) investigate the chemistry of deuterium-bearing molecules
in outer regions of protoplanetary disks and find that
molecules formed in the disk have similar D/H ratios
as those in comets. To date, the determination of D/H
ratios in disks has been limited to the measurements of
DCO+ /HCO+ in two classical T Tauri Stars (cTTs) with
single dish telescopes (TW Hydrae, van Dishoeck et al.
2003; DM Tauri, Guilloteau et al. 2006), but DCO+ has
not been observed in comets. Although H2 D+ and HDO
lines have been detected in disks (Ceccarelli et al. 2004,
2005), there are no corresponding H+
3 and H2 O observations that allow derivations of the D/H ratio. Therefore, direct measurements in disks of species observed in
comets, such as DCN and HCN, will be important to help
relate deuterium fractionation in disks to that in comets.
Deuterated molecules in cold pre-stellar and protostellar cores are found to be enhanced by orders of magnitude
over the elemental D/H abundance ratio of 1.5×10−5
through fractionation in the gas phase at low temperatures, an effect driven primarily by deuterated H+
3 , which
is formed by exchange reaction between H+
3 and HD,
then transfers its deuteron to neutral species. Theoretical models of disks (Aikawa et al. 2002; Willacy 2007)
with realistic temperature and density structure show
that the DCO+ /HCO+ ratio increases with radius due
to the decreasing temperature moving out from the star.
Spatially resolved data of deuterated molecules will help
test disk physical models through these chemical consequences.
The current capabilities of millimeter-wave observatories are limited both by sensitivity and by the small angular size of circumstellar disks, which makes spatially
resolving the deuterium fractionation difficult. In this
2
work we take advantage of the proximity of the most
nearby classical T Tauri star, TW Hydrae, which is surrounded by a disk of radius 3.′′ 5, or 200 AU at a distance
of 56 pc (Qi et al. 2004), and a possible orbiting planet
of mass (9.8±3.3) MJupiter at 0.04 AU (Setiawan et al.
2008), to study the physical and chemical structure of a
protoplanetary environment at high spatial and spectral
resolution using interferometry. The TW Hydrae disk is
viewed nearly face-on, and so is well posed to investigate
the radial distribution of molecular emission. Although
the current angular resolution of ∼1–2′′ places only a
few independent beams ( or ‘pixels’) across the disk, effectively improved resolution of the disk chemistry can
be achieved thanks to the fact that the gas kinematics are essentially Keplerian. Beckwith & Sargent (1993)
show that molecular line emission from Keplerian disks
displays a “dipole” pattern near the systematic velocity due to the shear created by the orbital motion. The
emission near the systemic velocity is dominated by the
outer disk regions, and the separation of the emission
peaks depends sensitively on the abundance gradient of
the emitting species with radius. This important feature,
commonly seen in velocity channel maps of disks with a
range of inclination angles (even as small as 6–7 degrees
in the case of TW Hydrae, Qi et al. 2004), may be used
to study the radial distribution of molecular emission far
more precisely than is afforded by the limited number of
“pixels” provided by the available resolution.
Here we report on Submillimeter Array (SMA)5 (Ho
et al. 2004) observations of deuterated species in the disk
around TW Hydrae, including the first detection and images of DCN and DCO+ . In §2 we describe the observations, while in §3 we introduce the analysis method used
and the molecular distribution parameters derived. We
describe the model fitting results and discuss the implications in §4, and we present a summary and conclusions
in §5.
2. OBSERVATIONS
All of the observations of TW Hydrae (R.A.:
11h 01m 51.s 875; Dec: -34◦ 42′ 17.′′ 155; J2000.0) were
made between 2005 February and 2006 December using
the SMA 8 antenna interferometer located atop Mauna
Kea, Hawaii. Table 1 and Table 2 summarize the observational parameters for the detected and undetected
species, respectively. The SMA receivers operate in a
double-sideband (DSB) mode with an intermediate frequency band of 4–6 GHz which is sent over fiber optic transmission lines to 24 “overlapping” digital correlator “chunks” covering a 2 GHz spectral window. Two
settings were used for observing the HCN/DCN and
HCO+ /DCO+ lines: at 265.7–267.7 GHz (lower sideband) the tuning was centered on the HCN J=3–2 line at
265.8862 GHz in chunk S22 while the HCO+ J=3–2 line
at 267.5576 GHz was simultaneously observed in chunk
S02; at 215.4–217.4 GHz (lower sideband) the tuning was
centered on the DCO+ J=3–2 line at 216.1126 GHz in
chunk S16 while the DCN J=3–2 line at 217.2386 GHz
was simultaneously observed in chunk S02 (the HDO
31,2 –22,1 line at 225.90 GHz was also covered in the USB
5 The Submillimeter Array is a joint project between the Smithsonian Astrophysical Observatory and the Academia Sinica Institute of Astronomy and Astrophysics, and is funded by the Smithsonian Institution and the Academia Sinica.
in chunk S06). Combinations of two array configurations
(compact and extended) were used to obtain projected
baselines ranging from 6 to 180 meters. Cryogenic SIS
receivers on each antenna produced system temperatures
(DSB) of 200-1400 K at 260 GHz and 100–200 K at 210
GHz. Observations of ortho-H2 D+ and H13 CO+ were
simultaneously carried out (in the dual receiver observational mode) on 28 December 2006 using only the compact array configuration. Three out of eight antennas
were equipped with 400 GHz receivers (DSB system temperature between 400 and 800 K) and were available for
observations of the 372.4213 GHz ortho-H2 D+ 11,0 –11,1
line, and all eight antennas were operating with 345 GHz
receivers (DSB system temperature between 150 and 300
K) for observations of the 346.9985 GHz H13 CO+ J=4–3
and 345.796 GHz CO J=3–2 lines. The correlator was
configured with a narrow band of 512 channels over the
104 MHz chunk, which provided 0.2 MHz frequency resolution, or 0.28 km s−1 at 217 GHz, 0.23 km s−1 at 267
GHz and 0.16 km s−1 at 372 GHz. Observations of the
691.6604 GHz para-D2 H+ 11,0 –10,1 line were shared with
observations of the CO J=6–5 line made on 17 February
2005 (details provided in Qi et al. 2006). Calibrations
of the visibility phases and amplitudes were achieved
with interleaved observations of the quasars J1037-295
and J1147-382, typically at intervals of 20-30 minutes.
Observations of Callisto provided the absolute scale for
the calibration of flux densities. The uncertainties in the
flux scale are estimated to be 15%. All of the calibrations were done using the MIR software package 6 , while
continuum and spectral line images were generated and
CLEANed using MIRIAD.
3. SPECTRAL LINE MODELING
Several model approaches have been developed to interpret interferometric molecular line and continuum observations from disks (Dutrey et al. 2007). In brief, the
approach of Qi et al. (2004, 2006) works as follows: the
kinetic temperature and density structure of the disk is
determined by modeling the spectral energy distribution
(SED) assuming well mixed gas and dust, with the results confirmed by the resolved (sub)mm continuum images. Then, a grid of models with a range of disk parameters including the outer radius Rout , the disk inclination i, position angle P.A. and the turbulent line width
vturb are produced and a 2D accelerated Monte Carlo
model (Hogerheijde & van der Tak 2000) is used to calculate the radiative transfer and molecular excitation.
The collisional rates are taken from the Leiden Atomic
and Molecular Database 7 (Schoier et al. 2005) for nonLTE line radiative transfer calculations. Specifically in
the models used in this paper, the rate coefficients for
HCO+ and DCO+ in collisions with H2 have been calculated by Flower (1999); the rate coefficients for HCN
and DCN in collisions with H2 have been scaled from the
rate of HCN-He calculated by Green & Thaddeus (1974).
The model parameters are fitted using a χ2 analysis in
the (u, v) plane.
In the studies by Qi et al. (2004, 2006), the distribution
of CO molecules was assumed to follow the H nuclei as
derived from the dust density structure and a gas of solar
6
7
http://www.cfa.harvard.edu/∼cqi/mircook.html
http://www.strw.leidenuniv.nl/∼moldata
3
TABLE 1
Observational Parameters for SMA TW Hydrae (detected species)
Rest Frequency (GHz):
Observations
Antennas used
Synthesized beam:
Channel spacing (km s−1 ):
R.M.S.a (Jy beam−1 ):
Peak intensityb (Jy)
a
b
c
HCN 3–2
265.886
2005 Mar 04
2005 Apr 21
2005 April 26
7
1.′′ 6 × 1.′′ 1 PA -0.5◦
0.23
0.35
4.7
HCO+ 3–2
267.558
2005 Mar 04
2005 Apr 21
2005 Apr 26
7
1.′′ 6 × 1.′′ 1 PA -6.3◦
0.23 km
0.29
2.0
H13 CO+ 4- -3
346.999
2006 Dec 28
DCN 3–2c
217.239
2006 Apr 28
DCO+ 3–2
216.113
2006 Apr 28
2006 Feb 03
8
4.′′ 1 × 1.′′ 8 PA 3.3◦
0.70
0.16
0.70
8
5.′′ 9 × 3.′′ 2 PA -1.5◦
0.56
0.10
0.31
8
2.′′ 6 × 1.′′ 6 PA 2.8◦
0.28
0.10
0.56
SNR limited by the dynamic range.
Intensity averaged over the corresponding beam.
Only compact configuration data used for DCN observation.
TABLE 2
Observational Parameters for SMA TW Hydrae (undetected species)
HDO 31,2 –22,1
225.897
2006 April 28
8
5.′′ 7 × 3.′′ 1 PA -0.6◦
0.54
0.11
<2.0 × 1014
composition. However, molecules in disks do not necessarily share the same distribution as molecular hydrogen.
Theoretical models (e.g. Aikawa et al. 1996; Aikawa &
Nomura 2006) predict so-called three-layered structure;
most molecules are photo-dissociated in the surface layer
of the disk and frozen out in the mid-plane where most
of hydrogen resides, with an abundance that peaks in the
warm molecular layer at intermediate scale heights. To
approximate this more complex behavior, we introduce
new molecular distribution parameters for use in spectral
line modeling.
For first-order analysis in the radial molecular distributions, the column densities are assumed to vary as a
power law as a function of radius:
r
)pi
Σi (r) = Σi (10AU)(
10AU
Here Σi and pi describe column density distribution of
a specific species (i) rather than disk (hydrogen) column
density. In most disk models to date, a single radial
power-law is adopted for fitting the hydrogen or dust surface density. By assuming CO follows the distribution of
hydrogen column density, Qi et al. (2004, 2006) also find
satisfactory power-law fits for multiple CO transtions.
When a single power law is insufficient for fitting the
radial distribution, then a broken power law with two
different indices will be used.
To calculate the surface density at different radii, the
vertical molecular distribution is needed, but this distribution may well vary with distance from the star.
However, theoretical models (Aikawa & Nomura 2006)
indicate that the vertical distribution of molecules at
different radii is similar as a function of the hydrogen
column density measured from the disk surface, Σ21 ≡
ΣH /(1.59 × 1021 cm−2 ), where the denominator is the
conversion factor of the hydrogen column density to Av
for the case of interstellar dust. As indicated by Figure
o-H2 D+ 11,0 –11,1
372.421
2006 Dec 28
3
4.′′ 7 × 3.′′ 8 PA 14.5◦
0.65
1.17
<1.7 × 1012
p-D2 H+ 11,0 –10,1
346.999
2005 Feb 17
4
3.′′ 3 × 1.′′ 3 PA 7.5◦
0.35
7.39
<9.0 × 1014
102
Normalized Fraction
Rest Frequency (GHz):
Observations:
Antenna used:
Synthesized beam:
Channel spacing (km s−1 ):
R.M.S. (Jy beam−1 ):
Nmolecule (1 σ) (cm−2 ):
100
10-2
10-4
10-6
10-8
0.001
σs
0.010
0.100
Σ21
σm
1.000
10.000 100.000
Fig. 1.— This schematic diagram shows an arbitrary molecular
vertical distribution as a function of Σ21 measured from the disk
surface at a certain radius. σm and σs are the midplane and surface
boundary parameters for model fitting (see §3 of the text).
5 of Aikawa & Nomura (2006), the vertical distribution
of molecular abundances shows a good correlation with
Σ21 . Figure 1 shows a schematic diagram of an arbitrary distribution of a molecule in disk where the x-axis
shows Σ21 and the y-axis shows the normalized molecular fraction. Smaller Σ21 values (to the left of the plot)
approach the surface of the disk, while larger Σ21 values
(to the right) approach the midplane of the disk. For
modeling, we make the further simplifying assumption
that gaseous molecules exist with a constant abundance
in layers between σs and σm , the surface and midplane
boundaries of Σ21 , respectively, as indicated in Figure 1
by the shaded area. While the vertical distribution of
molecules in disks may have a more complex distribution than assumed here, the adopted parameters provide
a gross approximation to the vertical location where the
species is most abundant. This is adequate for a first
description, given the quality of the data available. For
4
4
0.6
HCO+ 3-2
DCO+ 3-2
0.4
Flux (Jy)
Flux (Jy)
3
2
1
0.2
0.0
0
-0.2
-2
0
2
4
Velocity (km/s)
6
6
-2
8
HCN 3-2
2
4
Velocity (km/s)
6
8
DCN 3-2
0.3
4
0.2
Flux (Jy)
Flux (Jy)
0
2
0.1
0.0
0
-0.1
-2
-2
0
2
4
Velocity (km/s)
6
8
-2
0
2
4
Velocity (km/s)
6
8
Fig. 2.— Top: The HCO+ , DCO+ , HCN and DCN J=3–2 spectra at the peak continuum (stellar) position. The fluxes of HCO+ and
DCO+ are averaged over the beam of DCO+ J=3–2 (2.′′ 6 × 1.′′ 6 PA 2.8◦ ). The fluxes of HCN and DCN are averaged over the beam of
DCN J=3–2 (5.′′ 9 × 3.′′ 2 PA -1.5◦ ). The vertical dotted lines indicate the positions of the fitted VLSR for each molecular transition except
for DCN 3-2 where we adopt the VLSR from that of HCN J=3–2. Bottom: Velocity channel maps of the HCO+ , DCO+ , HCN and DCN
J=3–2 emissions toward TW Hydrae. The angular resolutions are 1.′′ 6 × 1.′′ 1 at PA -6.3◦ for HCO+ J=3–2 and 1.′′ 6 × 1.′′ 1 at PA -0.5◦ for
HCN J=3–2. The cross indicates the continuum (stellar) position. The axes are offsets from the pointing center in arcseconds. The 1σ
contour steps are 0.4, 0.12, 0.35 and 0.09 Jy beam−1 for HCO+ , DCO+ , HCN and DCN J=3–2 respectively and the contours start at 2σ.
example, the model of Aikawa & Nomura (2006) shows
two vertical peaks of HCO+ , but the secondary peak is
an order of magnitude smaller, and provides effectively
negligible emission. Also, at least for the nearly face-on
disk of TW Hydrae, the uncertainties in vertical distributions do not affect the derived radial distribution.
This simple model captures the basic characteristics of
three-layered structure predicted by theoretical models.
Using the new distribution parameters: Σi (10AU) and
pi (radial), σs and σm (vertical), the radial and vertical
distributions of the molecules in disks can be constrained
by observation. The best fit model is obtained by mini-
mizing the χ2 , the weighted difference between the real
and imaginary part of the complex visibility measured at
the selected points of the (u, v) plane. The χ2 values are
computed by the simultaneous fitting of channels covering LSR velocities from 1 to 4 km s−1 . Rout and pi are
calculated on grids in steps of 5 (AU) and 0.2, respectively. log (σs ) and log (σm ) are calculated on grids in
steps of 0.2 within the range from −2 to 2. Σi (10AU) and
i are found with each pair of radial and vertical distribution parameters by minimizing χ2 . The systemic velocity
VLSR is fit separately since it is not correlated with those
distribution parameters. For each fit parameter, the 1σ
5
√
uncertainties are estimated as χ21σ = χ2m + 2n, where
n is the number of degrees of freedom and χ2m is the χ2
value of best-fit model, as in Isella et al. (2007). We
tested values for the turbulent velocity in the range from
0.0 to 0.15 km s−1 and found that exact value does not
have a significant impact, in part because of the coarse
spectral resolution of the data. Therefore, we fixed the
turbulent velocity at an intermediate value, 0.08 km s−1 .
4. RESULTS AND DISCUSSION
40
140
120
30
80
6
10
60
7
0
10 3
40
8
10
20
20
10
109
6
0
20
40
50
Z(AU)
100
7
10 10
8
0
4190
10
30
50
100
150
R(AU)
200
Fig. 3.— The solid and dotted contours show the temperature
and density profiles from the TW Hya model of Calvet et al. (2002).
The blue and red lines confine the locations of HCO+ and HCN in
the model (see text).
Figure 2 shows the spectra of HCO+ , DCO+ , HCN
and DCN at the stellar (peak continuum) position and
their velocity channel maps. Table 3 summarizes the
power-law fitting results on HCO+ , DCO+ and HCN.
The DCN and H13 CO+ lines are too weak for a χ2 analysis of their distribution parameters, and so only the
ratios of DCN/HCN and HCO+ /H13 CO+ are fit. Figure 3 shows the density and temperature contours of the
disk model (adopted from Calvet et al. 2002) and the
locations of HCO+ and HCN derived from the model fitting procedure. Figure 4 shows the radial distribution of
the molecular column densities of the best fit models for
HCO+ , DCO+ and HCN. Not surprisingly, we found that
the radial distributions are better constrained than the
vertical ones: the local minimum of χ2 for different grids
of vertical parameters are within the noise limit, i.e. the
vertical parameters are the least well determined, due in
part to the face-on nature of the TW Hydrae disk. We
therefore treat the best-fit vertical results as fixed, and
investigate the uncertainty of other parameters.
4.1. HCO+ and DCO+
The first detection of DCO+ in the disk around
TW Hydrae was obtained by van Dishoeck et al. (2003)
with the James Clerk Maxwell Telescope (JCMT), who
reported a beam-averaged DCO+ /HCO+ abundance ratio of 0.035. Our spatially resolved observations of the
DCO+ and HCO+ J=3–2 emission from the disk suggest
a more complex chemical situation.
Figure 5 shows the channel maps of HCO+ J=3–2, together with the best-fit model and the data-model residuals. Table 3 lists the best-fit model parameters. The χ2
surface for the radial power index pHCO+ and the outer
radius Rout is shown in the top panel of Figure 6. The
1σ contour confines Rout to 200±10 AU and pHCO+ to
−2.9±0.3. The H13 CO+ J=4–3 line was also detected,
but the emission is not strong enough to constrain its distribution. Assuming that H13 CO+ has the same distribution as HCO+ , then fitting the ratio of HCO+ /H13 CO+
to match the intensity of the H13 CO+ emission indicates
an HCO+ /H13 CO+ ratio of 100±25. This value is consistent with the nominal solar system value of 89. Figure 7
presents the best-fit model spectra of H13 CO+ J=4–3
overlayed on the SMA data.
The vertical distribution of DCO+ is even less well constrained than that of HCO+ , and we choose to adopt the
same values of σs and σm for both HCO+ and DCO+
(as indicated in Table 3). Because of the nearly face-on
viewing geometry for the disk, the fitting of radial distribution power-law index pi is not affected significantly by
the uncertainties in vertical structure, but only the value
of Σi (10AU) needs to be adjusted.
Examination of the channel maps shows upon inspection that the radial distributions of DCO+ and HCO+
are different. To demonstrate how differences in the radial distribution affect the resulting images, Figure 8
presents three models of DCO+ radial distribution and
Figure 9 shows the corresponding simulated channel
maps of the DCO+ J=3–2 line.
Model 1 assumes that the DCO+ distribution follows the best-fit model of HCO+ .
The minimum
χ2 is determined to be 575396 and the corresponding
DCO+ /HCO+ is found to be 4.7 × 10−2 . Comparing the
simulated maps from Model 1 with the data in Figure 9
shows a distinct difference in that the double peaked nature of the central channel in the data is more obvious
than in this model. The contrast of the contour levels between the central channel and the adjacent channels are
also smaller in the data than in this model. Because the
emission of the central channel mostly originates at large
disk radii, these differences suggest that the DCO+ emission arising from the outer regions of the disk is stronger
than predicted by this model. The slope of radial DCO+
distribution does not decrease as steeply as does HCO+ .
In other words, the D/H ratio must increase with increasing radius.
Model 2 shows the best-fit result for the radial distribution of DCO+ assuming a single power index. As
shown in Figure 4, the radial distribution of DCO+ is
strikingly different from that of HCO+ , with a positive
power index of 2.4 and a smaller but better constrained
Rout 90±5 AU. The contours of the iso-χ2 surface for
DCO+ in the middle panel of Figure 6 indicate that the
uncertainty of the radial power index is very large but
that the index is still larger than 1.6 within the 1σ error, much larger than −2.9 found for HCO+ . The simulated DCO+ channel maps for Model 2 shown in Figure 9
are an improvement over Model 1 in matching the data.
However in this model, Rout (∼ 90 AU) is much smaller
than the disk radius observed with HCO+ and CO (∼
200 AU) and DCO+ increases with radius but then disappears sharply at Rout ∼90 AU, as a step function, which
is hard to understand. Comparison with the data shows
that there are fewer complete contours around the double peaks in the central channel, which suggests that the
DCO+ emission is maximum at an intermediate radius
6
1014
1.000
DCN/HCN=1.7±0.5×10-2
[DCO+]/[HCO+]
Column Density [cm-2]
HCO+/H13CO+=100±25
1013
0.100
[HCO+]
[HCN]/10
1012
0.010
[DCO+]
[DCO+]/[HCO+]
1011
10
0.001
100
R [AU]
Fig. 4.— Radial distribution of molecular column densities and DCO+ /HCO+ ratio for the best-fit models for TW Hydrae. Solid lines
depict single power-law fits, while the dashed lines are for DCO+ Model 3 where two power-law indices are used.
TABLE 3
Fitting results
Parameters
Stellar Mass: M∗ (M⊙ )
Inclination: i(deg)
Systemic velocity: VLSR (km s−1 )
Position angle: P.A.(deg)
Turbulent line width: vturb (km s−1 )
Outer radius: Rout (AU)
Column Density at 10AU: Σ(10AU) (cm−2 )
Radial power index: p
Vertical Parameters: σs ,σm
Minimum χ2 :
Reduced χ2 : χ2r
a DCO+ Model 2.
b DCO+ Model 3, no error estimation.
c DCO+ turning point at radius 70 AU.
HCO+
0.6
6.8±0.3
2.88±0.05
-27.4
0.08
200 ± 10
3.8 ±0.5 × 1015
-2.9 ± 0.3
0.1, 10
329588
1.21
rather than at the edge.
For Model 3, instead of using a single power index (p)
to fit for the radial distribution of DCO+ column density
N(DCO+ ), the fit uses two power indices (p1 and p2) and
a turning point (Rt ) where the power law index changes
from p1 to p2, i.e. the location of the peak of N(DCO+ ).
The parameters p1, p2 and Rt are searched within limited grids to minimize χ2 . In this model, N(DCO+ ) is
found to increase with radius out to ∼ 70 AU, and then
to decrease. Table 3 presents the best-fit parameters;
no error estimation is provided due to the computation
difficulties. Figure 9 shows the best-fit model images.
Comparison of Model 3 with the data shows a visual improvement over Model 2 in the central channel, although
the χ2 value of the two models are not distinguishable.
We thus cannot clearly discriminate with the χ2 statistic
if there is indeed a peak with radius for N(DCO+ ) as
DCO+a
0.6
7.4±0.6
2.94±0.06
-27.4
0.08
90 ± 5
1.9±0.2 × 1010
2.4 ± 0.8
0.1, 10
575347
1.95
DCO+b
0.6
7.4
2.9 4
-27.4
0.08
160
4.8 × 106
7,-6 c
0.1,10
575347
1.95
HCN
0.6
6.6±0.8
2.73±0.06
-27.4
0.08
100 ± 10
2.4±0.4 × 1014
-1.0 ± 1.2
0.3, 30
297198
1.35
in Model 3, or if DCO+ increases with radius and disappears suddenly around 90 AU as in Model 2. Even
with this ambiguity, however, both models imply that
the D/H ratios change by at least an order of magnitude (0.01 to 0.1) from radii <30 AU to >70 AU and
that there is a rapid falloff of N(DCO+ ) at radii larger
than 90 AU. Because the emission in the central channel
comes from the outer part of the disk (Keplerian rotation) projected along the line-of-sight with a very small
inclination of around 7 degrees for TW Hydrae, the central velocity channel is most important for constraining
the radial distribution. Based on the difference in the
central velocity channel map between the models, we believe Model 3 is the more plausible description of the
radial distribution of DCO+ . Figure 10 shows the channel maps of the DCO+ J=3–2 emission, together with
those of Model 3 and the data-model residuals.
7
Fig. 5.— Top: Velocity channel maps of the HCO+ J=3–2 emission toward TW Hydrae. The angular resolution is 1.′′ 6 × 1.′′ 1 at PA
-6.3◦ . The cross indicates the continuum (stellar) position. The axes are offsets from the pointing center in arcseconds. The 1σ contour
step is 0.4 Jy beam−1 and the contours start at 2σ. Middle: channel map of the best-fit model with the same contour levels. Bottom:
difference between the best-fit model and data on the same contour scale.
Observations of deuterated molecular ions and the level
of deuterium fractionation have been used to estimate the
ionization degree in molecular clouds, and a similar analysis can be applied to circumstellar disks. If we consider
only the ionization balance determined by HCO+ , H+
3,
DCO+ and electrons in steady state as shown in Equation 14 of Caselli (2002), the electron fractional abundance can be derived to be around 10−7 . Of course,
this value is only valid in the intermediate layer where
HCO+ is abundant and multiply deuterated H+
3 is less
abundant than HCO+ . Several important complications
are also neglected in this analysis, including the presence
of other atomic and molecular ions, neutral species besides CO which destroy H2 D+ , and negatively charged
dust grains and refractory metals. Still, accurate measurements of DCO+ and HCO+ are the first steps toward
an understanding of the ionization fraction in the disk.
The increase of D/H ratios from inner to outer disk is
generally consistent with the current theoretical models
of the gas-deuterium fractionation processes that consider the effect of cold temperature. But the quick disappearance of DCO+ at radii beyond 90 AU (comparing
with Rout around 200 AU for CO and HCO+ ) is puzzling,
since DCO+ is expected to be abundant and observable
in the cold outer region of the disk where HCO+ is still
available. More theoretical work is needed to explain the
disappearance of DCO+ in the outer part of the disk.
4.2. HCN and DCN
Figure 11 shows the HCN J=3–2 channel maps, together with the best fit model and residuals. Table 3
lists the best-fit model parameters, and Figure 4 shows
the radial distribution of column density derived. The
χ2 surface shown in the bottom panel of Figure 6 indicates Rout 100±10 AU and pHCN −1.0±1.2. The radial
power index is poorly constrained probably due to more
complex distributions for HCN. A detailed comparison
of the molecular distributions of HCN and CN will be
presented elsewhere.
Although the best fit vertical parameters seem to indicate that HCN (σm,HCN = 30) is found much deeper
toward the midplane than is HCO+ (σm,HCO+ = 10), we
emphasize that we are not able to accurately constrain
the vertical distributions from the present data. This ambiguity strongly affects the column density of HCN (i.e.
ΣHCN (10AU)) needed to fit the data (not the power index pHCN of radial distribution). A worse fit to the data
(χ2 larger by 3σ over the best-fit model) can be obtained
by assuming that the vertical distributions of HCN and
HCO+ are the same, but the HCN column density is 1.5
times larger than that needed for the best fit model due
to higher density near the midplane.
The DCN J=3–2 transition is detected at a signal-tonoise ratio of 3 near the fitted HCN VLSR of 2.73 km s−1
(Figure 2). While this signal-to-noise ratio is not high,
the significance of the detection is further supported by
the channel maps (Figure 12 upper panel) where the velocity gradient along the disk position angle of ∼ −30◦
is consistent with that seen in CO J=2–1 and J=3–2 (Qi
et al. 2004) and the other molecular lines presented in
this paper. Since the DCN 3–2 emission is weak, we
are not able to fit for the molecular distribution and so
make the simplifying assumption that the distribution of
DCN follows that of HCN. As with H13 CO+ , we fit the
DCN/HCN ratio over the whole disk and determine the
DCN/HCN ratio to be 1.7±0.5 × 10−2 . To again empha-
8
-2.0
H13CO+ 4-3
0.6
HCO
-2.5
0.4
Flux (Jy)
Radial Power Index
+
-3.0
-3.5
0.0
-4.0
-0.2
150
200
250
Outer Radius (AU)
300
-5
DCO+
2
5
Velocity (km/s)
10
1013
1
0
-1
70
80
90
100
110
Outer Radius (AU)
120
Model 2
1012
1011
Model 1
Model 3
1010
109
2
HCN
1
100
R [AU]
Fig. 8.— Models of DCO+ with different distributions of radial
column densities.
0
-1
-2
-3
70
Column Density [cm-2]
Radial Power Index
0
Fig. 7.— The beam-averaged H13 CO+ J=4–3 spectra at the
continuum (ste llar) position. The SMA data are presented by the
solid histogram, and the simulated model by the dashed histogram.
The vertical dotted line indicates the position of the fitted VLSR
for HCO+ J=3–2.
3
Radial Power Index
0.2
80
90 100 110 120 130 140
Outer Radius (AU)
Fig. 6.— Iso-χ2 surfaces of (Rout ,pi ) for HCO+ , DCO+ (as in
Model 2) and HCN. Contours correspond to the 1 to 6 σ errors. For
DCO+ , the index values larger than 3 at around 90 AU indicate
the ratios of DCO+ /HCO+ larger than 1, so the χ2 surface is not
calculated beyond that.
size the impact of the assumed vertical distribution on
the derived column densities, the DCN/HCN ratio could
be as high as 5 × 10−2 if HCN and DCN are distributed
vertically over the same region as is HCO+ .
Highly fractionated DCN/HCN ratios have been measured in comets. In the coma of comet Hale-Bopp, for example, Meier et al. (1998) reported the ratio to be around
2.3 × 10−3 , but higher DCN/HCN ratios – (D/H)HCN,jet
≈0.025 are detected from the pristine material sublimed
from icy grains ejected in jets from the nucleus which
may present a more representative sampling of cometary
ices that have not experienced significant thermal processing (Blake et al. 1999). Such ratios are consistent
with those found here in the TW Hydrae disk, indicating that high D/H ratios in comets could originate from
material in the outer regions of disks where in situ deuterium fractionation is ongoing, rather than requiring an
inheritance from interstellar material.
4.3. Upper limits for H2 D+ , D2 H+ and HDO
In the disk midplane, H2 is expected to be gaseous
and the molecular ion formed by the cosmic ray ionization of H2 , H+
3 , is expected to be the most abundant
ion. Unfortunately, H+
3 is only detectable in cold gas
via infrared absorption. In its deuterated forms, however, H2 D+ and even D2 H+ are expected to be abundant in the cold, dense gas (Ceccarelli & Dominik 2005).
The ground-state transition of ortho-H2 D+ was first detected in a young stellar object (NGC 1333 IRAS 4A)
9
Fig. 9.— DCO+ J=3–2 channel maps toward TW Hydrae and the simulated model distributions depicted in Figure 8.
by Stark et al. (1999) and in a prestellar core (L1544)
by Caselli et al. (2003). Both H2 D+ and D2 H+ have
been detected toward another pre-stellar core 16293E via
their ground-state submm rotational lines (Vastel et al.
2004). The inclusion of multiply deuterated H+
3 in chemical models leads to predictions of higher values of the
D/H ratio in cold, high-density regions of the interstellar
medium. Similarly, in the dense, cold disk midplane, CO
is depleted, and high abundances of H2 D+ and D2 H+
are expected. For this reason, Ceccarelli et al. (2004)
searched for the ground-state transition of ortho-H2 D+
with the Caltech Submillimeter Observatory (CSO), and
reported 3.2σ and 4.7σ detections toward TW Hydrae
and DM Tauri, respectively. With the 400 GHz and 690
GHz receiver-equipped SMA antennas, we searched for
the 372 GHz ortho-H2 D+ 11,1 –11,0 and 690 GHz paraD2 H+ 110 –101 lines toward TW Hydrae. No significant
emission signals were detected. Here we discuss the upper limits and their implications.
Our 3-antenna SMA observations give a 1σ upper
limit for the ortho-H2 D+ 11,1 –11,0 line emission of 1.2
Jy beam−1 km s−1 with a 4.′′ 7 × 3.′′ 8 synthesized beam.
To compare the result with single dish data, the extent
of the source emission must be known. Since H2 D+ was
observed along with the H13 CO+ 4–3 line in a dualreceiver observation on 28 December 2006 and H13 CO+
4–3 has been clearly detected at JCMT (van Dishoeck
et al. 2003), we can compare the intensities of these two
lines between the SMA data and the single dish observations to constrain the emitting areas. The H13 CO+ 4–3
line was detected at JCMT with an integrated intensity
of 0.07 K km s−1 in a 13′′ beam (van Dishoeck et al.
2003); while for the SMA the integrated intensity of this
line is determined to be 0.61 Jy beam−1 km s−1 with
a beam of 4.′′ 1 × 1.′′ 8. If the extent of the emission of
H13 CO+ is similar to that of H2 D+ , our interferometric
upper limit of Jy beam−1 km s−1 for H2 D+ (considering the small change of the beam sizes) becomes 0.09
K km s−1 (1σ) or 0.27 K km s−1 (3 σ) upper limits,
which is consistent with 2σ upper limits of about 0.2
K km s−1 by the JCMT (Thi et al. 2004) but slightly
lower than the 3.2σ detection of 0.39 K km s−1 by Ceccarelli et al. (2004). We estimate the 1σ upper limit of
the ortho-H2 D+ column density to be 1.7 × 1012 cm−2
according to the Equation 4 of Vastel et al. (2004), which
is slightly less than the 2σ upper limit estimate of 4.4 ×
1012 cm−2 by Thi et al. (2004) since the SMA data had
a smaller noise level. Of course this analysis assumes the
extent of H2 D+ is similar to that of H13 CO+ . With the
deployment of more 400 GHz receivers on the SMA and
further H2 D+ observations, it should be possible to provide rather better constraints on the H2 D+ abundance
10
Fig. 10.— Top: Velocity channel map of the DCO+ J=3–2 toward TW Hydrae. The angular resolution is 2.′′ 6 × 1.′′ 6 at PA 2.8◦ . The
cross indicates the continuum (stellar) position. The axes are offsets from the pointing center in arcseconds. The 1σ contour step is 0.12 Jy
beam−1 and the contours start at 2σ. Middle: channel map of Model 3 with the same contour levels. Bottom: difference between Model
3 and data on the same contour scale.
Fig. 11.— Top: Velocity channel map of the HCN J=3–2 toward TW Hydrae. The angular resolution is 1.′′ 6 × 1.′′ 1 at PA -0.5◦ . The
cross indicates the continuum (stellar) position. The axes are offsets from the pointing center in arcseconds. The 1σ contour step is 0.35
Jy beam−1 and the contours start at 2σ. Middle: channel map of the best-fit model with the same contour levels. Bottom: difference
between the best-fit model and data on the same contour scale.
11
0.3
DCN 3-2
Flux (Jy)
0.2
0.1
0.0
-0.1
-0.2
-5
0
5
Velocity (km/s)
10
Fig. 12.— Top left: DCN J=3–2 velocity channel maps (red: 2.84 km s−1 , blue: 2.28 km s−1 ) from TW Hydrae, overlaid on the 217
GHz dust continuum map (gray scale). The cross indicates the position of the continuum peak. Top right: the simulated model for DCN
J=3–2. The 1σ contour step is 0.09 Jy beam−1 and the contours start with 2σ. Bottom: the beam-averaged DCN 3–2 spectra at the
continuum (stellar) position. The SMA data are presented by the solid histogram, and the simulated model by the dashed histogram. The
vertical dotted line indicates the position of the fitted VLSR for HCN J=3–2.
12
in the disk.
We estimate the 1σ upper limit for para-D2 H+ 110 –101
to be 5.35 Jy beam−1 km s−1 with a beam of 3.′′ 3 × 1.′′ 3.
The 1σ upper limit to the para-D2 H+ column density is
estimated to be 9.0 × 1014 cm−2 . This is less constrained
than H2 D+ due to the relatively poor system sensitivity
at 690 GHz.
For HDO 31,2 –22,1 , the 1σ upper limit is 0.10
Jy beam−1 km s−1 . Assuming an excitation temperature
of 30 K, the upper limit for the HDO column density is
2.0 × 1014 cm−2 . Since the lower state energy level of
this line is nearly 160 K, it must originate from warm
regions of the disk which are quite distinct from the cold
layers where HDO ground state transition absorption, as
found in DM Tauri (Ceccarelli et al. 2005), must arise –
although we note that the detection of the HDO absorption line in DM Tauri has been disputed by Guilloteau
et al. (2006).
5. SUMMARY AND CONCLUSIONS
Observations of deuterated species in circumstellar
disks are important to understand the origin of primitive
solar system bodies in that they can directly constrain
the deuterium fractionation in the outer regions where
cometary ices are likely formed. Spatially resolved observations of the D/H ratios in disks enable the comparison
of the fractionation measured in comets such as HaleBopp (Blake et al. 1999) with the specific conditions at
each disk radius. We have presented the first images
of DCO+ and DCN emission from the disk around a
classical T Tauri star, TW Hydrae, along with images
of the HCN and HCO+ J=3–2 lines. The observations
of deuterium fractionation serve as a clear measure of
the importance of low-temperature gas-phase deuterium
fractionation processes. These observations strongly support the proposed link among high gas densities, cold
temperature and enhanced deuterium fractionation. De-
tailed chemical models are still needed to explain how
DCO+ disappears from the outer part of the disk.
The similarity of the D/H ratios in cold clouds, disks
and pristine cometary material has been used to argue
that the gas spends most of its lifetime at low temperatures and is incorporated into the disks before the envelope is heated, i.e. before the Class I stage. By combining
self-consistent physical models and 2D radiative transfer
codes to interpret high spatial resolution millimeter-wave
molecular images, we are only now beginning to investigate the radial and vertical distributions of molecules in
disks. The radial distribution of DCO+ in the disk of
TW Hydrae indicates that in situ deuterium fractionation is ongoing. The molecular evolution within disks
must therefore be considered in the investigation of the
origin of primitive solar system bodies.
We have obtained less stringent constraints on the vertical distributions of molecules in the disk of TW Hydrae.
To address the ambiguity present in the analysis of single
objects, data from a robust sample of disks is needed, in
particular one that covers a range of disk inclinations.
More sensitive observations are also needed for the rare
isotopologues H13 CN, H13 CO+ and, of course, DCN, to
understand the radial and vertical gradient of deuterated
species in these disks. In the future, observations of DCN
and other species with the Atacama Large Millimeter Array will provide further insight into the chemical state of
protoplanetary disks.
Partial support for this work comes from NASA Origins of Solar Systems Grant NNG05GI81G. M.R.H is
supported by a VIDI grant from the Netherlands Organization for Scientific Research. C.Q. acknowledges Paola
Caselli for her help and useful suggestions. We thank the
referee for very useful comments.
REFERENCES
Aikawa, Y. & Herbst, E. 1999, ApJ, 526, 314
Aikawa, Y., Miyama, S. M., Nakano, T., & Umebayashi, T. 1996,
ApJ, 467, 684
Aikawa, Y. & Nomura, H. 2006, ApJ, 642, 1152
Aikawa, Y., van Zadelhoff, G. J., van Dishoeck, E. F., & Herbst,
E. 2002, A&A, 386, 622
Beckwith, S. V. W. & Sargent, A. I. 1993, ApJ, 402, 280
Bergin, E. A., Aikawa, Y., Blake, G. A., & van Dishoeck, E. F.
2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt,
& K. Keil, 751–766
Blake, G. A., Qi, C., Hogerheijde, M. R., Gurwell, M. A., &
Muhleman, D. O. 1999, Nature, 398, 213
Calvet, N., D’Alessio, P., Hartmann, L., Wilner, D., Walsh, A., &
Sitko, M. 2002, ApJ, 568, 1008
Caselli, P. 2002, Planet. Space Sci., 50, 1133
Caselli, P., van der Tak, F. F. S., Ceccarelli, C., & Bacmann, A.
2003, A&A, 403, L37
Ceccarelli, C. & Dominik, C. 2005, A&A, 440, 583
Ceccarelli, C., Dominik, C., Caux, E., Lefloch, B., & Caselli, P.
2005, ApJ, 631, L81
Ceccarelli, C., Dominik, C., Lefloch, B., Caselli, P., & Caux, E.
2004, ApJ, 607, L51
Dutrey, A., Guilloteau, S., & Ho, P. 2007, in Protostars and Planets
V, ed. B. Reipurth, D. Jewitt, & K. Keil, 495–506
Flower, D. R. 1999, MNRAS, 305, 651
Green, S., & Thaddeus, P. 1974, ApJ, 191, 653
Guilloteau, S., Piétu, V., Dutrey, A., & Guélin, M. 2006, A&A,
448, L5
Ho, P. T. P., Moran, J. M., & Lo, K. Y. 2004, ApJ, 616, L1
Hogerheijde, M. R. & van der Tak, F. F. S. 2000, A&A, 362, 697
Isella, A., Testi, L., Natta, A., Neri, R., Wilner, D., & Qi, C. 2007,
A&A, 469, 213
Meier, R., Owen, T. C., Jewitt, D. C., Matthews, H. E., Senay,
M., Biver, N., Bockelee-Morvan, D., Crovisier, J., & Gautier, D.
1998, Science, 279, 1707
Qi, C., Ho, P. T. P., Wilner, D. J., Takakuwa, S., Hirano, N.,
Ohashi, N., Bourke, T. L., Zhang, Q., Blake, G. A., Hogerheijde,
M., Saito, M., Choi, M., & Yang, J. 2004, ApJ, 616, L11
Qi, C., Wilner, D. J., Calvet, N., Bourke, T. L., Blake, G. A.,
Hogerheijde, M. R., Ho, P. T. P., & Bergin, E. 2006, ApJ, 636,
L157
Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black,
J. H. 2005, A&A, 432, 369
Setiawan, J., Henning, T., Launhardt, R., Müller, A., Weise, P., &
M., K. 2008, Nature, 451, 38
Stark, R., van der Tak, F. F. S., & van Dishoeck, E. F. 1999, ApJ,
521, L67
Thi, W.-F., van Zadelhoff, G.-J., & van Dishoeck, E. F. 2004, A&A,
425, 955
van Dishoeck, E. F., Thi, W.-F., & van Zadelhoff, G.-J. 2003, A&A,
400, L1
Vastel, C., Phillips, T. G., & Yoshida, H. 2004, ApJ, 606, L127
Willacy, K. 2007, ApJ, 660, 441