HD far infrared emission as a measure of protoplanetary disk mass

TL;DR

This study evaluates hydrogen deuteride (HD) emission as a reliable tracer for measuring the gas mass in protoplanetary disks, demonstrating its potential accuracy and outlining observational requirements for future missions.

Contribution

It provides a detailed analysis of HD emission's robustness as a gas mass indicator, including the effects of disk structure and observational strategies, which was previously less understood.

Findings

HD 1-0 line intensity scales with gas mass with a slope of ~0.8.

Gas mass can be estimated within a factor of 2 for low-mass disks using HD flux alone.

Adding the HD 2-1 line or structural info reduces uncertainty to a factor of ~3.

Abstract

Protoplanetary disks around young stars are the sites of planet formation. While the dust mass can be estimated using standard methods, determining the gas mass - and thus the amount of material available to form giant planets - has proven to be very difficult. Hydrogen deuteride (HD) is a promising alternative to the commonly-used gas mass tracer, CO. We aim to examine the robustness of HD as tracer of the disk gas mass, specifically the effect of gas mass on the HD FIR emission and its sensitivity to the vertical structure. Deuterium chemistry reactions relevant for HD were implemented in the thermochemical code DALI and models were run for a range of disk masses and vertical structures. The HD J=1-0 line intensity depends directly on the gas mass through a sublinear power law relation with a slope of ~0.8. Assuming no prior knowledge about the vertical structure of a disk and using only the HD 1-0 flux, gas masses can be estimated to within a factor of 2 for low mass disks (M$_{\rm disk} < 10^{-3}$ M$_\odot$). For more massive disks, this uncertainty increases to more than an order of magnitude. Adding the HD 2-1 line or independent information about the vertical structure can reduce this uncertainty to a factor of ~3 for all disk masses. For TW Hya, using the radial and vertical structure from Kama et al. 2016b the observations constrain the gas mass to $6\cdot10^{-3}$ M$_\odot$ < M$_{\rm disk} < 9\cdot10^{-3}$ M$_\odot$. Future observations require a 5$\sigma$ sensitivity of $1.8\cdot10^{-20}$ W m$^{-2}$ ($2.5\cdot10^{-20}$ W m$^{-2}$) and a spectral resolving power R > 300 (1000) to detect HD 1-0 (HD 2-1) for all disk masses above $10^{-5}$ M$_\odot$ with a line-to-continuum ratio > 0.01. These results show that HD can be used as an independent gas mass tracer with a relatively low uncertainty and should be considered as an important science goal for future FIR missions.

Tables

Table 1. Table 1 : Dust mass opacities

	wavelength [$\mu$m]	${\kappa }_{\mathrm{abs}.}$ [cm2 g-1]
small grains	112	29.9
	56	154
large grains	112	30.0
	56	46.3

Table 2. Table 2 : DALI parameters of the physical model.

Parameter	Range
Chemistry
Chemical age	1 Myr
$\left[D\right]/\left[H\right]$	$2\cdot {10}^{-5}$
Physical structure
$\gamma$	1.0
$\psi$	[0.1,0.2,0.3,0.4,0.5]
$h_c$	[0.05,0.1,0.2,0.3] rad
$R_c$	35 AU
$R_{\mathrm{out}}$	200 AU
Mdisk	$2.6\cdot {10}^{-5},2.6\cdot {10}^{-4},7.7\cdot {10}^{-4},$
	$2.6\cdot {10}^{-3},7.7\cdot {10}^{-3},2.3\cdot {\mathrm{𝟏𝟎}}^{-\mathrm{𝟐}}$
	$7.7\cdot {10}^{-2},2.6\cdot {10}^{-1}$ M⊙
Dust properties
Gas-to-dust ratio	200
$f_{\mathrm{large}}$	0.99
$\chi$	0.2
$f_{\mathrm{PAH}}$	0.001
Stellar spectrum1
type	T Tauri
LX	1030 erg s-1
LFUV	$2.7\cdot {10}^{31}$ erg s-1
L∗	0.28 L⊙
${\zeta }_{\mathrm{cr}}$	${10}^{-19}{s}^{-1}$
Observational geometry
$i$	6∘
d	140 pc

Table 3. Table 3 : Fit coefficients of the HD flux-mass relation.

	HD 1-0	HD 2-1
$\alpha$	0.79	0.54
A	5.44$\cdot {10}^{-19}$	2.83$\cdot {10}^{-19}$
B	9.95$\cdot {10}^{-19}$	3.50$\cdot {10}^{-19}$
Mtr (M⊙)	1.32$\cdot {10}^{-3}$	5.43$\cdot {10}^{-4}$

Table 4. Table 4 : Disk parameters for the fixed dust mass models.

Parameter	Range
Dust masses
Mdust	${10}^{-5},{10}^{-4},{10}^{-3}$ M⊙
Gas-to-dust ratios
${\Delta }_{\mathrm{gd}}$	10, 100, 1000
Gas masses1
Mgas	${10}^{-5}-1$ M⊙
Vertical structure
$(\psi ,h_c)$	$(0.1,0.1),(0.3,0.1),(0.3,0.3)$

Table 5. Table 5 : Definition vertical structure classification

	$\psi$	$h_c$
“weakly flared”	0.1-0.2	0.05-0.1
“intermediately flared”	0.3-0.5	0.05-0.1
	or1
	0.1-0.2	0.2-0.3
“strongly flared”	0.3-0.5	0.2-0.3

Table 6. Table 6 : Simple HD chemical network reactions and adopted rate coefficient.

Reaction	$\alpha$	$\beta$	$\gamma$	Ref.
D+ + H $\to$ H+ + D	1.0(-9)	0.0(0)	0.0(0)	1,2
D + H+ $\to$ H + D+	1.0(-9)	0.0(0)	4.1(1)	1,2
D+ + H2 $\to$ H+ + HD	2.1(-9)	0.0(0)	0.0(0)	1,2
D+ + H2 $\to$ H + HD+	1.75(-9)	9.8(-2)	0.0(0)	1,2
HD + H+ $\to$ H2 + D+	1.0(-9)	0.0(0)	6.345(2)	1,2
HD+ + H $\to$ HD + H+	6.4(-10)	0.0(0)	0.0(0)	2
H${}_{}{}^{+}{}_2{}^{}$ + D $\to$ H2 + D+	6.4(-10)	0.0(0)	0.0(0)	2
D + cr $\to$ D+ + e-	4.6(-20)	0.0(0)	0.0(0)	2
HD + cr $\to$ H + D+ + e-	4.73(-20)	0.0(0)	0.0(0)	2
HD + cr $\to$ D + H+ + e-	4.73(-20)	0.0(0)	0.0(0)	2
HD + cr $\to$ D + H	5.99(-19)	0.0(0)	0.0(0)	2
HD + cr $\to$ D+ + H-	8.4(-22)	0.0(0)	0.0(0)	2
HD + cr $\to$ HD+ + e-	5.16(-18)	0.0(0)	0.0(0)	2
HD+ + H2 $\to$ H${}_{}{}^{+}{}_3{}^{}$ + D	5.25(-9)	0.0(0)	0.0(0)	2
H${}_{}{}^{+}{}_2{}^{}$ + HD $\to$ H${}_{}{}^{+}{}_3{}^{}$ + D	5.25(-9)	0.0(0)	0.0(0)	2
HD + He+ $\to$ He + H+ +D	5.5(-14)	-2.4(-1)	0.0(0)	2
HD + He+ $\to$ He + D+ +H	5.5(-14)	-2.4(-1)	0.0(0)	2
D+ + H $\to$ HD+ + $\gamma$	3.9(-19)	1.8(0)	0.0(0)	2
H+ + D $\to$ HD+ + $\gamma$	3.9(-19)	1.8(0)	0.0(0)	2
HD+ + H $\to$ H${}_{}{}^{+}{}_2{}^{}$ + D	1.0(-9)	0.0(0)	1.54(2)	2
H${}_{}{}^{+}{}_2{}^{}$ + D $\to$ HD+ + H	1.0(-9)	0.0(0)	0.0(0)	2
HD+ + D $\to$ HD+ + D	6.4(-10)	0.0(0)	0.0(0)	2
HD+ + e- $\to$ H + D	3.4(-9)	-4.0(-1)	0.0(0)	2
D+ + e- $\to$ D + $\gamma$	3.61(-12)	-7.5(-1)	0.0(0)	2
H + D $\to$ HD + $\gamma$	3.0(-18)	5.0(-1)	0.0(0)	∗
D+ + Mg $\to$ D + Mg+	1.1(-9)	0.0(0)	0.0(0)	∗
D+ + Si $\to$ D + Si+	9.9(-10)	0.0(0)	0.0(0)	∗
D+ + Fe $\to$ D + Fe+	7.4(-9)	0.0(0)	0.0(0)	∗
D+ + S $\to$ D + S+	1.3(-9)	0.0(0)	0.0(0)	∗
D+ + PAH $\to$ D + PAH+	3.1(-8)	0.0(0)	0.0(0)	∗
D+ + PAH- $\to$ D + PAH	8.3(-7)	-5.0(-1)	0.0(0)	∗
HD+ + PAH $\to$ HD + PAH+	3.1(-8)	0.0(0)	0.0(0)	∗
HD+ + PAH- $\to$ HD + PAH	8.3(-7)	-5.0(-1)	0.0(0)	∗

Equations

$${\Sigma }_{\mathrm{gas}}(R)={\Sigma }_c{\left(\frac{R}{R_c}\right)}^{\gamma }\exp \left[-{\left(\frac{R}{R_c}\right)}^{2-\gamma }\right].$$

$$n(R,z)$$

$$=\frac{1}{\sqrt{2\pi }}\frac{1}{Rh}\exp \left[-\frac{1}{2}{\left(\frac{z/R}{h}\right)}^2\right],$$

$$h=h_c{\left(\frac{R}{R_c}\right)}^{\psi }.$$

$$R_{\mathrm{diss},\mathrm{HD}}=1.5\times 10^{-9}I(\nu ){\mathrm{s}}^{-1}.$$

$$F_{\mathrm{HD}}=\{\begin{array}{ll}A{\left({\mathrm{M}}_{\mathrm{disk}}/{\mathrm{M}}_{\mathrm{tr}}\right)}^{\alpha },& \text{if }{\mathrm{M}}_{\mathrm{disk}}\le {\mathrm{M}}_{\mathrm{tr}}\\ A+B{\log }_{10}\left({\mathrm{M}}_{\mathrm{disk}}/{\mathrm{M}}_{\mathrm{tr}}\right),& \text{otherwise}\end{array}$$

$$L/C\equiv \frac{\int \mathrm{d}\nu F_{\mathrm{line}}}{\int \mathrm{d}\nu F_{\mathrm{cont}}}\simeq \frac{F_{\mathrm{peak}}}{F_{\mathrm{cont},\mathrm{bin}}}.$$

$${\log }_{10}\Delta \mathrm{M}(F)=\frac{1}{2}\left|{\log }_{10}{\mathrm{M}}_{\max }(F)-{\log }_{10}{\mathrm{M}}_{\min }(F)\right|.$$

$$\rho \frac{v_{\varphi }}{r}$$

$$0$$

$$\varphi =-\frac{G{M}_{\ast }}{{\left(r^2+z^2\right)}^{1/2}},$$

$$\frac{\mathrm{d}P}{\mathrm{d}z}$$

$$\frac{1}{P}\frac{\mathrm{d}P}{\mathrm{d}z}$$

$$P(z)=c_s(z{)}^2\rho (z)=C\exp \left(\int \frac{zG{M}_{\ast }}{c_s^2{\left(r^2+z^2\right)}^{3/2}}\mathrm{d}z\right).$$

Full text

*Aims. *To examine the robustness of HD as tracer of the disk gas mass, specifically the effect of gas mass on the HD far infrared emission and its sensitivity to the vertical structure. Also, to provide requirements for future far-infrared missions such as SPICA.

*Methods. *Deuterium chemistry reactions relevant for HD were implemented in the thermochemical code DALI and more than 160 disk models were run for a range of disk masses and vertical structures.

Results. The HD J=1-0 line intensity depends directly on the gas mass through a sublinear power law relation with a slope of $\sim$0.8. Assuming no prior knowledge about the vertical structure of a disk and using only the HD 1-0 flux, gas masses can be estimated to within a factor of two for low mass disks (M${}_{\rm disk}\leq 10^{-3}$ M⊙). For more massive disks, this uncertainty increases to more than an order of magnitude. Adding the HD 2-1 line or independent information about the vertical structure can reduce this uncertainty to a factor of $\sim 3$ for all disk masses. For TW Hya, using the radial and vertical structure from Kama et al. (2016b) the observations constrain the gas mass to $6\cdot 10^{-3}\ \mathrm{M}_{\odot}\leq\mathrm{M}_{\rm disk}\leq 9\cdot 10^{-3}$ M*⊙. Future observations require a $5\sigma$ sensitivity of $1.8\cdot 10^{-20}$ W m-2* ($2.5\cdot 10^{-20}$ W m*-2*) and a spectral resolving power $R\geq 300\ (1000)$ to detect HD 1-0 (HD 2-1) for all disk masses above $10^{-5}\ \mathrm{M}_{\odot}$ with a line-to-continuum ratio $\geq 0.01$.

*Conclusions. * These results show that HD can be used as an independent gas mass tracer with a relatively low uncertainty and should be considered as an important science goal for future far-infrared missions.

1 Introduction

One of the key properties for understanding the evolutionary path of disks and the formation of planets is the disk gas mass (e.g., Armitage 2015). Determining gas masses is no easy task. The main component of the gas, molecular hydrogen (H2), is difficult to observe. As a symmetric rotor, H2 has no electric dipole moment, limiting its rotational lines to quadrupole transitions ($\Delta J=2$). The $J=2$ level lies at $549.2$ K, so in order to produce any appreciable H2 emission gas temperatures have to be at least 100 K (Thi et al. 2001; Carmona et al. 2011) which is much higher than the temperature of the bulk of the gas in the disk. In addition, these lines lie in the mid- and near-infrared (IR) where the continuum is bright, which further hampers the observations due to the high continuum optical depth and low line-to-continuum ratios. Even if detected, the H2 emission does not trace the bulk of the mass (e.g., Pascucci et al. 2013). This means that indirect tracers of the gas mass have to be used.

Historically, most of the disk mass estimates have been based on observations of the (sub-)millimeter (mm) emission of the dust grains (e.g., Beckwith et al. 1990; Dutrey et al. 1996; Mannings & Sargent 1997; Williams et al. 2005). In order to relate the observed flux to the mass of the emitting material, a mass opacity $\kappa_{\nu}$ and a temperature $T_{\rm dust}$ have to be assumed (cf. Hildebrand 1983). To obtain the total (gas+dust) mass, the estimate of the dust mass has to be scaled up by the gas-to-dust mass ratio. Usually the interstellar medium (ISM) value of 100 is assumed, but it is unclear whether this is applicable to protoplanetary disks.

The gaseous component of the disk can be traced using carbon monoxide (CO), the second most abundant molecule in the disk, which is readily detected towards many protoplanetary disks (e.g., Dutrey et al. 1996; Thi et al. 2001; Dent et al. 2005; Panić et al. 2008; Williams & Best 2014). To derive a gas mass from CO observations, a relation has to be found between the amount of CO in the disk and the amount of H2. In the simplest case where most of the available carbon is contained in gas-phase CO, the fractional abundance of CO is $n({\rm CO})/n({\rm H}_{2})\sim 2\cdot 10^{-4}$, in line with observations of warm dense clouds (e.g., Lacy et al. 1994).

In protoplanetary disks there are two processes that decrease the CO abundance below its maximum value (van Zadelhoff et al. 2001). In the upper layer of the disk, CO is being destroyed through photodissociation by ultraviolet (UV) photons. In the midplane, where temperatures are below $\sim 20$ K, CO freezes out onto the dust grains. As a result of these non-linear processes, thermochemical models have to be used to relate CO observations to the total gas mass. These processes are well understood and they are implemented in current thermochemical models. However, studies that include these effects find that overall gas phase carbon must still be additionally depleted by a factor 10-100. This leads to low CO emission even for massive disks (e.g., Bruderer et al. 2012; Favre et al. 2013; Du et al. 2015; Kama et al. 2016b; Bergin et al. 2016).

Several processes have been proposed to explain the underabundance of gas-phase carbon in regions where freeze-out is not important. A possible explanation could be grain growth, where carbon has been locked up in large icy bodies that no longer participate in the gas-phase chemistry (see, e.g., Bergin et al. 2010, 2016; Du et al. 2015; Kama et al. 2016b). Another cause could be the conversion of CO into more complex species. This process can happen in the gas-phase through reactions between CO and He+, which extracts atomic carbon that can then be used to construct more complex molecules. These molecules have higher freeze-out temperatures than CO and freeze out onto the cold grains, thus lowering the apparent carbon abundance (Aikawa et al. 1997; Favre et al. 2013; Bergin et al. 2014). Similarly, ice chemistry can play an important role in converting CO into more complex organics, such as CH3OH or turning it into CO2 or CH4 ice (see, e.g., Figure 3c in Eistrup et al. 2016).

Due to its high abundance, 12CO becomes optically thick at the disk surface, meaning it does not trace the bulk of the mass. Instead, less abundant CO isotopologues such as 13CO, C18O and C17O have to be used, which remain optically thin throughout the disk (van Zadelhoff et al. 2001; Dartois et al. 2003). Visser et al. (2009) and Miotello et al. (2014) showed that C18O/C16O ratio is reduced beyond the 18O/16O isotope ratio in disks due to isotope-selective photodissociation. Not taking this effect into account can change the estimated mass using C18O up to one order of magnitude (Miotello et al. 2014).

Using the Herschel Space Observatory, hydrogen deuteride (HD) has been detected towards TW Hya, DM Tau and GM Aur (Bergin et al. 2013; McClure et al. 2016). HD has been suggested as an alternative tracer of the disk gas mass. Because of their chemical similarities, the distribution of HD is expected to closely follow that of H2. HD has a small dipole moment, so dipole transitions ($\Delta J=1$) are allowed. The difference between the energy needed to excite the first rotational level of HD ($E/k_{B}\simeq 128.5$ K) and the energy needed to excite the second rotational level of H2 ($E/k_{B}\simeq 549.2$ K), means that at a temperature of $T\sim 20$ K the expected emission of HD is much larger than that of H2.

The energetically lowest transition of HD, $J$=1-0, emits at 112 $\mu$m (Müller et al. 2005). This places the transition in a wavelength range where the dust continuum is bright, the atmospheric opacity is high and the production of low-noise detectors is expensive. All these reasons combined make detecting HD a challenge. Even with Herschel, deep integrations were required to observe HD in disks. After the end of the Herschel mission, there are currently no far-IR observatories capable of detecting HD in disks. However, this may change with the advent of future far-IR missions such as the proposed SPICA mission (Nakagawa et al. 2014) and with the HIRMES instrument on SOFIA.

In this work, a simple HD chemistry is incorporated into the thermochemical code Dust and Lines (DALI, Bruderer et al. 2012; Bruderer 2013). In Section 2 a description of the code is given, including the deuterium chemistry that was implemented. A description of the models that were run is also given here. In Section 3 the dependence of the HD far-infrared emission on the disk gas mass is examined. In addition, the effect of other disk properties such as the vertical structure and the dust distribution on the HD emission is investigated. Furthermore, the results are examined in the light of possible future far-IR missions such as SPICA. The effect of prior knowledge of the disk vertical structure on the uncertainty of the HD based mass estimates is discussed in Section 4. Additionally, the models are compared to the current observations. Of particular interest is the protoplanetary disk of TW Hya ($59.54\pm 1.45$ pc, Astraatmadja & Bailer-Jones 2016), where the gas mass estimated using HD is more than one order of magnitude higher than the gas mass estimated using C18O (Favre et al. 2013). Different studies have shown that the TW Hya disk seems genuinely depleted in volatile carbon and oxygen (e.g., Bergin et al. 2010, 2013, 2016; Hogerheijde et al. 2011; Favre et al. 2013; Du et al. 2015; Kama et al. 2016b). Recently, new observations of the CO isotopologues have become available (see, e.g., Schwarz et al. 2016). Therefore, the case of TW Hya is revisited using a combined HD and CO isotopologue analysis. The conclusions can be found in Section 5.

2 Model

The thermochemical code Dust and Lines (DALI) was used to calculate the line fluxes for the disk models (Bruderer et al. 2012; Bruderer 2013). DALI is a 2D physical-chemical code that for a given physical disk model calculates the thermal and chemical structure self-consistently. In addition to the disk density structure, a stellar spectrum has to be provided in order to determine the UV radiation field inside the disk. The computation proceeds through three steps. First the dust temperature structure and the internal radiation field (from UV- to mm-wavelengths) are determined by solving the continuum radiative transfer using a 2D Monte Carlo method. Next, the abundances of the molecular and atomic species at each point of the disk are obtained from the solution of the time-dependent chemistry. The atomic and molecular excitations are then computed using a non-LTE calculation. Given the excitation levels, the gas temperature is determined by balancing the heating and cooling processes. Due to the fact that both the chemistry and the excitations depend on temperature, the problem is solved iteratively until a self-consistent solution is found. Lastly, the model is ray-traced to construct spectral image cubes and line profiles. A more detailed description of the code can be found in Appendix A of Bruderer et al. (2012).

2.1 Density structure

The density structure used in the code follows the simple parametric model proposed by Andrews et al. (2011). This model is based on the assumption that the disk structure is determined by viscous accretion, where the viscosity goes as $\nu\propto R^{\gamma}$ (Lynden-Bell & Pringle 1974; Hartmann et al. 1998). The resulting surface density is given by

[TABLE]

Here $\Sigma_{c}$ is the surface density at the characteristic radius $R_{c}$.

The vertical structure is given by a Gaussian density distribution, as motivated by hydrostatic equilibrium (Kenyon & Hartmann 1987)

[TABLE]

where $H=Rh$ is the physical height of the disk and the scale height angle $h$ is parametrized by

[TABLE]

Here $h_{c}$ is the scale height at $R_{c}$ and $\psi$ is known as the flaring angle. Equation (4) is a representation for a physical disk, where the vertical structure is set by balancing gravity and the vertical pressure gradient at each point in the disk.

2.2 Dust settling

In order to include dust settling into the model, the dust grains are split into two populations:

• •

small grains (0.005-1 $\mu$m) are included with a (mass) fractional abundance $1-f_{\rm large}$. These are located throughout the vertical extent of the disk.

• •

large grains (1-$10^{3}\ \mu$m) are included with a fractional abundance $f_{\rm large}$. They are limited to a vertical region with scale height $\chi h;\ \chi<1$, simulating the effect of dust settling.

Note that in the current version of DALI the gas surface density profile is only coupled to the density profile of the small grains, in contrast to earlier versions of DALI. The resulting disk structure is summarized in Figure 1.

A standard ISM dust composition is assumed for the dust opacities following Weingartner & Draine (2001), in line with Bruderer (2013) (see e.g. his Figure 2). The mass opacities for absorption at 112 $\mu$m and 56 $\mu$m for both dust grain populations are given in Table 1.

2.3 Chemical network

For the models in this work the list of chemical species adopted in Miotello et al. (2014) was extended to include HD, D, HD+ and D+ as separate species. By adding the element D to H, He, 12C, 13C, N, 16O, 17O, 18O, Mg, Si, S and Fe, the total number of species was increased to 280. The reaction types included are H2 and HD formation on dust, freeze-out, thermal and non-thermal desorption, hydrogenation of simple species on ices, gas-phase reactions, photodissociation, X-ray and cosmic-ray induced reactions, PAH/small grain charge exchange/hydrogenation, and reactions with H${}_{2}^{*}$ (vibrationally excited H2). The details of the implementation of these reactions are described in Appendix A.3.1 of Bruderer et al. (2012).

To properly simulate HD integrated line fluxes, simple deuterium chemistry has been added to the chemical network. The only D-bearing species considered for the chemical calculation are D, D+, HD, and HD+. Accordingly, the important reactions regulating their abundances have been included, with rate coefficients from Roberts & Millar (2000); Walmsley et al. (2004). Photodissociation of HD was included following Glover & Jappsen (2007). The photodissociation rate for a radiation field with a flat spectrum is given by Eq. 51 in Glover & Jappsen (2007):

[TABLE]

Here $I(\nu)$ is the mean intensity of the radiation field. Note that the prefactor has units such that the photodissociation rate has units [s*-1*].

HD self-shielding was also included in the model using the shielding function of H2 but slightly shifted in wavelength. Finally HD formation onto grains and ion-exchange reactions with metals and PAHs have been included, with the same rate coefficients as that assumed for the analogue reactions with H instead of D. The reactions involving D-bearing species added in the chemical network are presented in Table 6 in Appendix D.

2.4 Grid of models

As a basis for the models, the physical structure of the TW Hya disk from the Kama et al. (2016b) model is used. A combination of spatially and spectrally resolved and unresolved line fluxes, in combination with observations of the spectral energy distribution (SED), was used to constrain the parameters of their model, which can be found in Table 2. Note that Kama et al. (2016b) find a strongly flared vertical structure for TW Hya. Also, their model has a radially steeper temperature profile and smaller outer radius than the models of Cleeves et al. (2015).

In total, over a 160 models were run using DALI. The employed chemical network is the expanded version of that used by Miotello et al. (2014), where simple deuterium chemistry has been considered (see Section 2.3). The input total disk masses are M${}_{\mathrm{disk}}=[2.6\cdot 10^{-5},\ 2.6\cdot 10^{-4},\ 7.7\cdot 10^{-4},\ 2.6\cdot 10^{-3},\ 7.7\cdot 10^{-3},\ 2.3\cdot 10^{-2},\ 7.7\cdot 10^{-2},\ 2.6\cdot 10^{-1}]$ M*⊙*. For each mass, 20 models with different vertical structures were ran: $h_{c}\in[0.05,0.1,0.2,0.3],\psi\in[0.1,0.2,0.3,0.4,0.5]$ (cf. equation 4). The other parameters have been kept fixed to the values assumed by Kama et al. (2016b) (Table 2). For the overall carbon and oxygen abundances, typical ISM values of [C]/[H]${}_{\rm ISM}=1.35\cdot 10^{-4}$, [O]/[H]${}_{\rm ISM}=2.88\cdot 10^{-4}$ are used (Bruderer et al. 2012).

The stellar spectrum for TW Hya from Cleeves et al. (2015) was used as the central source of radiation. It closely resembles a 4000 K blackbody with a 10000 K blackbody component to represent the observed UV excess (following the prescription in Kama et al. 2016a, see also Figure 3 and Table 3 in Kama et al. 2016b).

The time-dependent chemistry is run for 1 Myr, which is less than the expected age for TW Hya, but the HD emission is not sensitive to the chemical age.

The overall deuterium abundance [D]/[H] of the gas in the disk has a direct impact on the HD fluxes. Prodanović et al. (2010) used observations of the atomic D and H lines in the diffuse ISM to find $(\mathrm{[D]/[H]})_{\rm ISM}\geq(2.0\pm 0.1)\times 10^{-5}$ for the local ISM (i.e., within $\sim 1-2$ kpc of the Sun). The uncertainties due to the error on the assumed [D]/[H] are much smaller than the spread in HD fluxes due to the disk structure. Note that for older star forming regions such as Orion the deuterium abundance is lower, due to the fact that a fraction of the deuterium has been burned up in the cores of stars (Wright et al. 1999; Howat et al. 2002).

3 Results

3.1 HD flux vs. disk gas mass

The first step in investigating the viability of HD as a tracer of the disk gas mass is examining the behaviour of its emission under variation of the gas mass. This is presented in Figure 2, where the integrated HD 1-0 and 2-1 line fluxes are plotted as function of disk mass. In the left panel of Figure 2 the blue dots shows the median HD 1-0 flux for each disk mass, while the coloured region presents fluxes between the 16th and 84th percentiles in each mass bin. The HD 1-0 line flux increases with mass, following a linear dependence in log-log space for the low-mass disks (M${}_{\rm disk}\leq 10^{-3}\ \mathrm{M}_{\odot}$).

Towards the high mass regime, the relation flattens. There are two reasons for this change in behaviour. For disk masses above 0.001 M*⊙* the HD 1-0 line emission starts to become optically thick, which means that the HD emission no longer traces the full gas reservoir. In addition, high mass disks have larger densities, which decreases the gas temperature in the disk. Miotello et al. (2016) found a similar decrease in line luminosity for 13CO and C18O for high mass disks, which they attributed to optical depth for 13CO and increased freeze-out for C18O.

The right panel of Figure 2 shows that the flux-mass relation for the HD 2-1 line is flatter compared to HD 1-0. This is likely due to increased temperature dependence of the HD 2-1 line, compared to the HD 1-0 line. The increased temperature dependence is reflected in the increased sensitivity of HD 2-1 to the vertical structure of the disk, as indicated by the large flux variations in each mass bin.

The behaviour of the flux-mass relation can be expressed with an analytical relation, such as:

[TABLE]

where Mtr is defined as the mass where the flux-mass relation starts to flatten.

This model was fitted to the median fluxes of both HD lines. The resulting coefficients can be found in Table 3. Note that the powerlaw slopes for both HD lines are sublinear. This is likely a result from the fact that the HD emission is produced by warm gas only.

3.2 HD emitting layers

The 2D abundance of HD in a typical model is presented in Figure 3. This distribution is the outcome of the deuterium chemical calculation discussed in Section 2.3. Note that $n(\mathrm{HD})/n_{\rm gas}=2\cdot 10^{-5}$ throughout most of the disk, indicating that all of the available deuterium is locked up in HD. The exceptions to this rule are in the upper layer of the disk, where HD is photodissociated and in an intermediate layer where HD abundance is decreased by a factor of $\sim$2. This layer is due to a combination of chemical processes. It should be noted here that this layer is very thin and only covers a very small fraction of the HD emission area shown in Figure 3. Therefore the depletion in this layer has no relevant effect on the HD flux.

The solid blue lines in Figure 3 denote where 25% and 75% of the HD 1-0 emission is produced. This region is elevated above the midplane ($0.1\leq z/r\leq 0.25$), indicating that most of the HD 1-0 emission is produced by warm gas, with an average temperature of $30-70$ K (cf. Figure 15). The emitting region extends radially from $r_{\rm in}\approx 9$ AU up to $r_{\rm out}\approx 70$ AU for M${}_{\rm disk}=2.3\cdot 10^{-2}\ \mathrm{M}_{\odot}$. The inner and outer radii vary slightly vertical structure (cf. Figure 16). Figure 17 shows that when the disk gas mass is increased the emitting region shifts outward, from $(r_{\rm in},r_{\rm out})\approx(1\ \mathrm{AU},40\ \mathrm{AU})$ for M${}_{\rm disk}=2.3\cdot 10^{-5}\ \mathrm{M}_{\odot}$ to $(r_{\rm in},r_{\rm out})\approx(10\ \mathrm{AU},100\ \mathrm{AU})$ for M${}_{\rm disk}=2.3\cdot 10^{-1}\ \mathrm{M}_{\odot}$. This is a result of the optically thick HD zone extending to larger radii as the gas mass increases.

The $\tau=1$ surfaces for the line and continuum optical depth are shown in dashed and dashed-dotted lines, respectively. They show that continuum optical depth has no effect on the line emission. This is likely due to the fact that most of the dust mass is in large grains, which are settled to the midplane in contrast with the models by Bergin et al. (2014). The dashed line shows that for this model with M${}_{\rm disk}=2.3\cdot 10^{-2}$ M*⊙*, the HD emission starts to become optically thick, which is in line with the flattening of the slope seen in Figure 2.

For the HD 2-1 emission, shown in light blue in Figure 3, the emission originates even higher up in the disk, where gas temperatures are in the range of $70-100$ K (cf. Figure 15). The contour for the line optical depth show that the HD 2-1 emission starts to become optically thick in the inner part of the disk for M${}_{\rm disk}=2.3\cdot 10^{-2}$ M*⊙*.

3.3 Influence of the vertical structure

The strength of the HD emission depends on the amount of material and the temperature of the material. In protoplanetary disks the temperature structure is set by the vertical structure of the disk (Chiang & Goldreich 1997). The effect of the vertical structure on the HD 1-0 and 2-1 fluxes is presented in Figure 4. It shows the HD fluxes as function of the flaring angle $\psi$ for different scale heights $h_{c}$ (colours). This was done for three different mass bins ($2.3\cdot 10^{-5}\ \mathrm{M}_{\odot}$, $2.3\cdot 10^{-3}\ \mathrm{M}_{\odot}$ and $2.3\cdot 10^{-1}\ \mathrm{M}_{\odot}$). The fluxes in each mass bin are normalized with respect to the median flux in the same bin (cf. the blue dots in Figure 2).

The normalized HD 1-0 integrated line fluxes are presented in top panels of Figure 4. It shows that the HD 1-0 flux increases with increasing $\psi$. The relation is roughly linear, with a slope that is independent of $h_{c}$ and Mdisk.

The spacing between the different $h_{c}$ lines indicate that the HD 1-0 flux also increases systematically with $h_{c}$. The shape of the flux-$h_{c}$ relation depends only weakly on $\psi$. It does depend strongly on disk mass, with the influence of $h_{c}$ on the flux increasing as function of disk mass.

For HD 2-1, presented in the bottom panels of Figure 4, the dependencies of the flux on $h_{c}$ and $\psi$ are stronger than those found for HD 1-0. The bulk of the HD 2-1 emission is produced by gas at temperatures of $70-100$ K (cf. Figure 15). This gas is located higher up in the disk. As the optically thin HD 2-1 emission depends on both the temperature and the density, the location of the HD 2-1 emitting material varies strongly with vertical structure, as can be seen in Figure 16. This is in contrast to the HD 1-0 emitting material, which stays at roughly the same location, independent of the vertical structure.

Figure 4 shows that the maximum variation in flux due to the different vertical structure increases with disk mass, from $0.75\ \times$ the median flux for $\mathrm{M}_{\rm disk}\sim 10^{-5}\ \mathrm{M}_{\odot}$ up to $1.9\ \times$ the median flux for $\mathrm{M}_{\rm disk}\sim 10^{-1}\ \mathrm{M}_{\odot}$. For the HD 2-1 integrated line flux, shown in the bottom panel of Figure 4, the maximum variation in flux is both larger than for HD 1-0 and is less dependent on disk mass. The maximum variation of the flux is $1.9\ \times$ the median flux for $\mathrm{M}_{\rm disk}\sim 10^{-5}\ \mathrm{M}_{\odot}$, increasing to $2.8\ \times$ the median flux for $\mathrm{M}_{\rm disk}\sim 10^{-1}\ \mathrm{M}_{\odot}$. As mentioned in Section 3.2, the regions where the HD lines are optically thick have larger outer radii as the gas mass is increased. On the one hand this pushes the HD emitting regions upwards to parts of the disk that are more sensitive to the vertical structure. On the other hand, the increased size of the optically thick HD zone means that, for massive disks, a larger mass fraction of HD will be much warmer as the vertical structure changes. The uncertainty in derived gas masses from HD due to vertical structure is further quantified in Section 4.1.

In the above models, the vertical structure is parametrized according to Equation (4). Another approach would be to calculate the vertical structure using hydrostatic equilibrium from the computed gas temperature structure. The effect of this approach on the HD fluxes is examined in Appendix C and is found to be within the flux variations shown in Figure 4.

3.4 Influence of the large grains

For the models discussed in Section 2.4 only three parameters (Mdisk, $\psi$ and $h_{c}$) were varied. However, several other parameters may also affect the HD flux. Their role is investigated here. In particular, the large dust grains in the disk can affect the HD flux, both by affecting heating and by changing the optical depth of the dust at the wavelengths where HD emits. Two dust parameters are examined here: the mass fraction of the large grains $f_{\rm large}$ and the dust settling parameter of the large grains $\chi$ (cf. Section 2.2). In both cases, the mass and vertical structure of the disk were fixed to the values from Kama et al. (2016b). The resulting HD 1-0 fluxes, normalized to the median flux of the M${}_{\rm disk}=2.3\cdot 10^{-2}\ \mathrm{M}_{\odot}$ mass bin, are shown in Figure 5.

Figure 5 shows that the HD 1-0 flux increases with dust settling (i.e., when the large grains are more settled towards the midplane). This can be understood by the fact the dust optical depth at a fixed height in the disk increases when the large grain population (cf. Section 2.2) is less settled. This is especially the case in the model for TW Hya, where most of the dust mass is in these grains.

A similarly large effect is found when the mass fraction of the large grains is changed. When $f_{\rm large}$ decreases, the amount of small grains increases significantly, which increases the opacity at 112 $\mu$m.

Comparing these results with those from the previous section, it is clear that both the large grains and the vertical structure have similar effects on the HD 1-0 line flux. It should be noted however that both $\chi$ and $f_{\rm large}$ are constrained by the spectral energy distribution (e.g., D’Alessio et al. 2001; Chiang et al. 2001), so their values cannot be chosen arbitrarily.

3.5 Influence of the gas-to-dust ratio

Deriving accurate gas-to-dust mass ratios has been one of the motivations for measuring the gas mass directly (e.g., Ansdell et al. 2016; Miotello et al. 2017). However, the gas-to-dust ratio could also affect the HD emission directly, e.g., by changing the coupling between the dust and the gas. If the effect is strong it may prevent accurate measurements of the gas-to-dust ratio using gas masses derived from HD emission. For fixed gas mass of M${}_{\rm gas}=2.3\cdot 10^{-2}\ \mathrm{M}_{\odot}$, Figure 6 presents the HD 1-0 fluxes as functions of increasing gas-to-dust ratio (and therefore decreasing dust mass). The fluxes are normalized with respect to the median flux of the M${}_{\rm disk}=2.3\cdot 10^{-2}\ \mathrm{M}_{\odot}$ mass bin.

The figure shows that the HD 1-0 flux increases systematically with $\Delta_{\rm gd}$. The primary reason for this is the decreased coupling between the gas and the dust for higher gas-to-dust ratios. As the gas in the HD emitting region is less efficient in cooling than the dust, this leads to slightly higher gas temperatures. Additionally, the increased dust mass for the models with a lower $\Delta_{\rm gd}$ enhances the dust opacities at the wavelengths where HD emits, which lowers the amount of observed line flux. However, inspection of the $\tau_{\rm dust}=1$ surface at 112 $\mu$m shows that the increase in opacity is relatively small.

3.5.1 Determining gas-to-dust ratios using HD 1-0

Figure 6 shows that varying the gas-to-dust ratio by two orders of magnitude for a fixed gas mass only leads to a variation in the HD 1-0 line flux of $\sim 1.5\times$ the median flux. In this analysis the dust mass is assumed to be a free parameter. However, in practice the dust mass of the observed disks will be constrained, e.g., by observation of their millimeter-continuum fluxes. A comparison of HD 1-0 fluxes and millimeter continuum fluxes for models with different $\Delta_{\rm gd}$ is shown in Figure 7.

Here the HD emission is combined with the millimeter continuum emission of the dust to provide a method for estimating the gas-to-dust ratio from observations. For this purpose an additional set of 27 models was run, where the dust mass, gas-to-dust ratio and the vertical structure were varied (cf. Table 4). Figure 7 shows the HD 1-0 line flux, which traces the gas mass, plotted against the 870 $\mu$m continuum flux, which traces the dust mass. The lack of overlap between models with different gas-to-dust ratios indicates that combined observations of HD and the dust continuum can be used to determine $\Delta_{\rm gd}$. For the three disks for which HD 1-0 observations are available (i.e., TW Hya, DM Tau and GM Aur) gas-to-dust ratios are found to be $\Delta_{\rm gd}\sim 100$, close to the ISM value. This would suggest that indeed the lower values found from similar measurements made using 13CO are due to an underabundance of volatile carbon (see Miotello et al. 2017 and references therein).

3.6 Line-to-Continuum ratios

As mentioned in the introduction, observing lines in the far-infrared is inherently difficult because the emission of the dust in protoplanetary disks peaks at $\sim 100\ \mu\mathrm{m}$. Superimposed on this bright continuum are the intrinsically weak far-infrared lines. As a result, the line-to-continuum ratio ($L/C$) is a crucial ingredient for observing the HD levels.

The line-to-continuum ratio can be defined as:

[TABLE]

Here $F_{\rm line}$ and $F_{\rm cont}$ are the continuum subtracted line flux and continuum flux, respectively. The approximation holds if the line is narrow, such that the integrated line flux can be approximated by the total flux in the center-most frequency bin ($F_{\rm peak}$) times the bin width. In this case $F_{\rm cont,bin}$ is the continuum flux in the same bin.

Due to the difficulties in characterizing the noise properties of infrared detectors as well as cosmic ray impacts on the detector and electronics, the current best achievable signal-to-noise ratio (SNR) for these detectors is $\sim$300. In the most optimal case, these detectors are able to distinguish variations at the level of $\frac{1}{300}\equiv 0.0033$ of the source continuum flux. By requiring that the HD lines are detected at the $3\sigma$ level, the SNR = 300 translates into a detection limit of $L/C$ $\geq 3\cdot 0.0033=0.01$, which can not be improved by longer integration times.

The spectral resolving power $R=\lambda/\Delta\lambda$ of the instrument plays a crucial role in determining the $L/C$. As long as the line is unresolved, a higher $R$ corresponds to a larger $L/C$, independent of the disk properties. This assumption is valid for the HD lines in protoplanetary disks, which are intrinsically narrow (FWHM $\sim 5$ km s*-1* in the outer disk). As explained in the introduction there is currently no facility capable of observing HD rotational transitions in disks. However, this work and previous studies (e.g., Bergin et al. 2013; McClure et al. 2016) have shown that HD is a unique tracer of the total gas mass. This means that the possibility of observing HD should be included in the consideration for instruments on future FIR observatories, such as the proposed SPICA mission. To be able to detect HD, the detection limit of $L/C\geq 0.01$ places a constraint on the requirements of the spectral resolving power of these future FIR missions.

The top panel of Figure 8 shows $L/C$ for HD 1-0 as function of disk mass for three different values of $R$. It is clear that the lower spectral resolving power mainly lowers the $L/C$. The current HD 1-0 detections with the Herschel PACS instrument are for $R=1000$ at 112 $\mu$m. This figure demonstrates that future instruments, such as the SPICA SAFARI instrument (Roelfsema et al. 2014), need a resolving power of $R\geq 300$ to detect HD 1-0 over the full mass range examined here, if a detection limit of $L/C\geq 0.01$ (so S/N on the continuum of 300) is assumed. In the case of the HD 2-1 transition, shown in the bottom panel of Figure 8 a spectral resolving power of $R\geq 1000$ is required to detect HD 2-1 towards all disk masses above $10^{-5}\ \mathrm{M}_{\odot}$.

3.7 Sensitivities of future FIR missions

In addition to good line-to-continuum ratios, future FIR missions also require good sensitivities to be able to detect the HD rotational lines in protoplanetary disks. In Figure 9 the models presented in Figure 2 are compared to the sensitivities of two FIR instruments: the proposed SPICA SAFARI instrument (Roelfsema et al. 2014) and the recently approved HIRMES instrument for SOFIA (E. Bergin, priv. comm.). The fluxes in the figure have been calculated for a distance of 140 pc, in line with typical distances to the closest star forming regions. The sensitivity of the Herschel PACS instrument is also shown here. The $5\sigma$ sensitivities after 1 hr of integration for SAFARI and Herschel PACS are shown in dark red and black respectively. The purple and orange lines represent the $5\sigma$ sensitivities at 112 $\mu$m after 10 hrs of integration for the SAFARI and HIRMES instruments, respectively. Figure 9 shows that the proposed SPICA SAFARI instrument will be able to detect HD 1-0 in disks with gas masses down to $5\cdot 10^{-4}$ M${}_{\odot}\approx 0.5\ \mathrm{M}_{\rm jup}$. If the integration time is increased to 10 hrs the mass sensitivity goes down to $8\cdot 10^{-5}$ M${}_{\odot}\approx 0.08\ \mathrm{M}_{\rm jup}$.

4 Discussion

4.1 Determining the disk gas mass

In Section 3.1 the behaviour of the HD emission under variation of the disk gas mass was discussed. Figure 2 shows that it is possible to constrain the gas mass using the HD 1-0 flux. However, the variations in the flux due to the uncertainties on vertical structure translate into uncertainties in the determined gas mass. In this section these uncertainties are quantified and the impact of complementary observations on these uncertainties is investigated.

Figure 2 shows that models with different combinations of disk gas masses and vertical structures can produce the same HD 1-0 flux. The difference between the minimum and the maximum disk mass, interpolated using the model fluxes shown in Figure 2, can be used to quantify the uncertainty on the gas mass estimates using the flux (cf. the small panel in Figure 10). Here the mass derivation uncertainty is defined as a factor of $10^{\Delta\mathrm{M}}$, where:

[TABLE]

Here $\mathrm{M}_{\rm min}(F)$, $\mathrm{M}_{\rm max}(F)$ are the minimum and maximum disk masses that can produce a HD 1-0 flux $F$, given the range of vertical structures examined.

The calculated mass uncertainties as function of HD 1-0 flux are presented in Figure 10. For the full set of models introduced in Section 2.4 (cf. Table 2), shown in black, $\Delta$M increases with HD flux. The increase is steep, with $\Delta$M $=1$ dex at a HD flux of $\sim 5\cdot 10^{-19}$ W m*-2*, indicating that this flux corresponds to an estimated mass between $10^{-3}\ \mathrm{M}_{\odot}$ and $10^{-1}\ \mathrm{M}_{\odot}$.

This estimate relies only on the HD 1-0 flux, without making any assumptions on the vertical structure of the disk. In practice, most disks will also have been observed at other wavelengths and with other chemical tracers. These can be used to put constraints on the vertical structure of the disk. This would in turn lead to lower $\Delta$M. If complementary observations can distinguish between “weakly flared” disks, “intermediately flared” disks and “strongly flared” disks (see Table 5 for the definitions), the uncertainties on the mass estimates can be significantly improved.

The top left panel of Figure 10 shows the three subsets with respect to the full set of models together with the mass uncertainties for each subset. For “intermediately flared” disks, shown in orange, $\Delta$M is significantly lower compared to the full set. Here the maximum mass uncertainty is $\sim 0.65$ dex, meaning that the mass can be estimated within a factor of $\sim 5$.

For the “strongly flared” disks, shown in light green, the uncertainty is minimized, with $\Delta$M never becoming larger than 0.2 dex. As a result, masses for this subset can be estimated to within a factor of two.

The mass uncertainty for the “weakly flared” disks follows the uncertainty of the full sample, with $\Delta$M being $\sim$0.1 dex lower than for the full sample. For these ranges of $h_{c}$ and $\psi$ the HD flux is quite sensitive to the vertical structure, which results in variations in the flux comparable to those in the full set. However, a difference of 0.1 dex is still a factor of two improvement on the mass uncertainty.

4.2 HD 1-0 and HD 2-1 line fluxes

An interesting complementary observable is the HD 2-1 line emission. The main disk parameter determining the difference between the $J=1$ and $J=2$ level populations of HD is the gas temperature, which is largely set by the vertical structure. By combining observations of HD 2-1 and HD 1-0 this difference can be exploited to determine the disk mass with low uncertainty.

As shown in Figure 11, the different dependencies of HD 1-0 and HD 2-1 line fluxes on disk mass and vertical structure result in a clear mass segregation, independent of vertical structure. For disks with M${}_{\rm disk}\geq 7.7\cdot 10^{-3}\ \mathrm{M}_{\odot}$ the mass bins start to overlap. As the bins differ by a factor of 3 in mass, the combination of two HD lines allow the gas mass to be determined to within a factor of $\sim$3, provided that the observational uncertainties are small enough. To be able to determine the mass accurately across the whole mass range, the observational uncertainties have to decrease with increasing disk mass. This is not a too strict requirement as the observed line flux increases with disk mass, allowing for a more precise flux determination. The observations of HD 1-0 and HD 2-1 for TW Hya shown in Figure 11 will be discussed in Section 4.4.

4.3 Comparing models to observations

The HD 1-0 transition has been observed for three disks: TW Hya (Bergin et al. 2013), DM Tau and GM Aur (McClure et al. 2016). These observations are shown in Figure 12, which represents the high mass part of the left panel of Figure 2. The observed flux for TW Hya was scaled to a distance of 140 pc, to allow for comparison with the same models.

For DM Tau the lower limit on disk gas mass based on the HD 1-0 observation is M${}_{\rm disk}\approx 4\cdot 10^{-3}\ \mathrm{M}_{\odot}$. This correspond to a high, strongly flared vertical structure. If instead the median vertical structure is assumed, the estimated gas mass becomes M${}_{\rm disk}\approx 1.5\cdot 10^{-2}\ \mathrm{M}_{\odot}$. For GM Aur, a similar analysis gives a lower limit on the disk gas mass of M${}_{\rm disk}\approx 1\cdot 10^{-2}\ \mathrm{M}_{\odot}$ for a high,strongly flared disk. Comparing the observed HD 1-0 flux to the medians, a gas mass of M${}_{\rm disk}\approx 8\cdot 10^{-2}\ \mathrm{M}_{\odot}$ is found.

Both estimates based on the medians are in line with the results of McClure et al. (2016), who found gas masses of $(1.0-4.5)\cdot 10^{-2}$ and $(2.5-19.5)\cdot 10^{-2}\ \mathrm{M}_{\odot}$ for DM Tau and GM Aur by modeling both the HD flux and the continuum SED.

The observations for TW Hya shown in Figure 12 will be discussed in Section 4.4.

4.4 Case study: TW Hya

The protoplanetary disk of TW Hya ($59.54\pm 1.45$ pc, Astraatmadja & Bailer-Jones 2016), one of the best-studied disks. Its disk gas mass has been estimated to be $\geq 0.05\ \mathrm{M}_{\odot}$ using HD (Bergin et al. 2013). This estimate is an order of magnitude higher than a similar estimate of the gas mass using C18O (Favre et al. 2013).

The observations of the HD 1-0 and 2-1 line fluxes for the TW Hya disk are presented in Figures 11 and 12. Comparing these observations to the models shows that the observed HD 1-0 flux for TW Hya places a lower limit of M${}_{\rm disk}\approx 3\cdot 10^{-3}\ \mathrm{M}_{\odot}$ if high, strongly flared vertical structure is assumed. Using the radial and vertical structure from Kama et al. (2016b), the observations constrain the disk gas mass to $6\cdot 10^{-3}\ \mathrm{M}_{\odot}\leq\mathrm{M}_{\rm disk}\leq 9\cdot 10^{-3}\ \mathrm{M}_{\odot}$. A similar gas mass is found from combining the observations of HD 1-0 and HD 2-1 (cf. Figure 11), which favor a gas mass between $7.7\cdot 10^{-3}\ \mathrm{M}_{\odot}\leq\mathrm{M}_{\rm disk}\leq 2.3\cdot 10^{-2}$ M*⊙. Both estimates are lower than previous estimates (e.g., Bergin et al. 2013; Kama et al. 2016b). Note that the estimates lie closer to the disk mass estimated from C18*O observations (Favre et al. 2013), suggesting a smaller tension between the two.

Using models that include detailed modeling of the CO and H2O chemistry, studies found that the disk surrounding TW Hya is genuinely depleted in volatile carbon and oxygen (e.g., Bergin et al. 2010, 2013; Hogerheijde et al. 2011; Favre et al. 2013; Du et al. 2015). Recent observations of [CI], [CII], [OI], C2H and C3H2 have reinforced this conclusion (Kama et al. 2016b; Bergin et al. 2016).

In their model, Kama et al. (2016b) incorporate the isotopes 13C, 17O, 18O and D parametrically. In this section, their best fitting model (M${}_{\rm disk}=2.3\cdot 10^{-2}\ \mathrm{M}_{\odot}$, $\psi=0.3$ and $h_{c}=0.1$) is rerun using the extended chemical network that includes the deuterium and CO isotopologue selective chemistry following Miotello et al. (2014). For these models, the carbon abundance was varied between the ISM abundance and two consecutive orders of underabundance, i.e., X(C) $\equiv$ [C]/[H] $\in[1.35\cdot 10^{-4},\ 1.35\cdot 10^{-5},\ 1.35\cdot 10^{-6}]$. The same was done for the oxygen abundance, X(O) $\in[2.88\cdot 10^{-4},\ 2.88\cdot 10^{-5},\ 2.88\cdot 10^{-6}]$. These underabundances are denoted using $\delta_{i}\equiv$ X($i$)/X($i$)ISM. Note that in order to match the hydrocarbon abundances Kama et al. (2016b) fine-tuned their final model to $\delta_{\rm C}=0.01$ and C/O = 1.5, while the disk structure was unchanged from the best-fit fiducial model.

Figure 13 shows the resulting line fluxes for several rotational transitions of HD, 13CO and C18O, similar to Figure 6 in Kama et al. (2016b). In the figure only the models with $\mathrm{X(C,O)}\sim 10^{-6},10^{-5}$ are shown. The full figure can be found in Appendix E. It should be noted that HD abundance used in Kama et al. (2016b) ($n(\mathrm{HD})/n(\mathrm{H}_{2})=1\cdot 10^{-5}$) is a factor four lower than the maximal HD abundance of the models in this work ($n(\mathrm{HD})/n(\mathrm{H}_{2})=4\cdot 10^{-5}$). In order to compare the output of the chemical calculations to their parametric method, the HD line fluxes of their model were recalculated using the higher HD abundance. This accounts for the differences in the HD line fluxes seen here and those in Kama et al. (2016b). Similar steps where taken so that the $\delta_{\rm C}=\delta_{O}=0.01$ model and the Kama et al. (2016b) model have the same oxygen and carbon abundances.

Starting with HD, one can see that the line fluxes from the models run here, which include the deuterium chemistry network, match the line fluxes from the model with the parametric HD within a few percent (8% and 15% for HD 1-0 and HD 2-1, respectively). This can be understood from the results in Section 3.2, which show that all the available deuterium is locked up in HD, thus mimicking the parametric approach used by Kama et al. (2016b). Note that all models overproduce the observed HD fluxes by a factor $\sim 2$, suggesting a slightly lower disk gas mass for TW Hya than used by Kama et al. (2016b).

For 13C, there is again good agreement between the fluxes from Kama et al. (2016b) and the $\delta_{\rm C}=\delta_{O}=0.01$ model (its CO isotopologue counterpart). This can be understood by reactions such as the ion-molecule reaction 13C*+* + 12CO $\rightleftharpoons$ 13CO + 12C*+* + 35 K. The forward reaction direction of this reaction is favoured at low temperatures, leading to an increased 13CO abundance which balances out the additional decrease of 13CO due to isotope-selective photodissociation (Miotello et al. 2014).

The line fluxes of C18O of the $\delta_{\rm C}=\delta_{O}=0.01$ model are systematically lower than the same lines from Kama et al. (2016b). This qualitatively agrees with previous results (e.g., Visser et al. 2009; Miotello et al. 2014), but the effect is smaller than has been suggested. A possible explanation can be found in the gas mass of the disk. Miotello et al. (2014, 2016) show that for disks with M${}_{\rm disk}\geq 7\cdot 10^{-3}\ \mathrm{M}_{\odot}$ freeze-out is the dominant process affecting the C18O flux, whereas isotope-selective photodissociation is more important for low mass disks (see, e.g., Fig. 3 in Miotello et al. 2016).

When the models are compared to the observations to determine the carbon and oxygen underabundance the results vary based on the tracer used. The 13CO observations point to $\delta_{\rm C}$ = 0.01, in line with the value found by Kama et al. (2016b). The observed 13CO 10-9 flux and the observations for C18O instead favour a lower underabundance of carbon, between $\delta_{\rm C}$ = 0.01 and $\delta_{\rm C}$ = 0.1. For all CO isotopologue lines $\delta_{\rm O}$ = 0.01 and $\delta_{\rm O}$ = 0.1 fit the observations equally well. Combined with the lower gas mass estimates found for TW Hya, these results show that the inclusion of isotope-selective processes decreases underabundance of carbon needed to explain the observations to $\delta_{\rm C}\approx 0.1$. It should be noted however that in this analysis only the CO isotopologues are considered, whereas the analysis in Kama et al. (2016b) also included observations of atomic carbon and atomic oxygen.

5 Conclusions

Accurately determining the amount of material in protoplanetary disks is crucial for understanding both the evolution of disks and the formation of planets in these disks. Recent observations of the rotational transitions of HD in protoplanetary disks (Bergin et al. 2013; McClure et al. 2016) have added a new possibility of tracing the disk gas mass. In this work, the robustness of HD as a tracer of the disk gas mass and its dependence on the vertical structure are investigated. The thermochemical code DALI (Bruderer et al. 2012; Bruderer 2013) was used to calculate the line fluxes for the disk models. The normal chemical network of DALI was expanded to include CO isotopologue selective chemistry, following Miotello et al. (2014). A simple D chemistry network was also implemented in DALI. A series of models was run spreading a range of disk masses (M${}_{\rm disk}\sim 10^{-5}-10^{-1}\ \mathrm{M}_{\odot}$) and vertical structures ($h_{c}\sim 0.05-0.3,\ \psi\sim 0.1-0.5$), with the large grains settled close to the midplane. From these models, the following observations can be made:

• •

The HD flux increases as a powerlaw for low mass disks (M${}_{\rm disk}\leq 10^{-3}$ M*⊙). The fitted powerlaw indices are 0.8 for HD 1-0 (112 $\mu$m) and 0.5 for HD 2-1 (56 $\mu$m). These indices are less than 1.0 due to the dependence of the emission on the gas temperature. For high mass disks (M${}_{\rm disk}>10^{-3}$ M⊙*) the HD flux scales with $\log_{10}\ \mathrm{M}_{\rm disk}$, which is a result of the increased line optical depth and decreased overall temperature of these disks.

• •

The maximum variation in HD 1-0 flux due to the different vertical structure increases with disk mass, from $0.75\ \times$ the median flux for $\mathrm{M}_{\rm disk}\sim 10^{-5}\ \mathrm{M}_{\odot}$ up to $1.9\ \times$ the median flux for $\mathrm{M}_{\rm disk}\sim 10^{-1}\ \mathrm{M}_{\odot}$. The variation for HD 2-1 is both larger than for HD 1-0 and is less dependent on disk mass. The maximum variation of this flux is $1.9\ \times$ the median flux for $\mathrm{M}_{\rm disk}\sim 10^{-5}\ \mathrm{M}_{\odot}$, increasing to $2.8\ \times$ the median flux for $\mathrm{M}_{\rm disk}\sim 10^{-1}\ \mathrm{M}_{\odot}$.

• •

The influence of the large grain population on the HD flux is less than that of the vertical structure, approximately 0.7$\times$ the median flux at M${}_{\rm disk}=2.3\cdot 10^{-2}\ \mathrm{M}_{\odot}$.

• •

Without making assumptions on the vertical structure of the disk, the HD 1-0 flux can be used to estimate the disk gas mass to within a factor of $\sim$3 for low mass disks (M${}_{\rm disk}\leq 10^{-3}$ M*⊙*). For more massive disks the uncertainty in the estimated mass increases to more than an order of magnitude.

• •

If complementary observations are employed to constrain the disk’s flaring, the uncertainty on the gas mass can be reduced down to factor of 2, even for massive disks. Moreover, a combination of HD 1-0 and the HD 2-1 fluxes can be used to a determine the disk gas mass to within a factor of 3, without making assumptions on the vertical structure.

• •

For DM Tau and GM Aur, gas mass estimates found by comparing the observed fluxes to the models agree with the results from McClure et al. (2016), confirming the high gas masses of these disks.

• •

For TW Hya, a gas mass between $6\cdot 10^{-3}\ \mathrm{M}_{\odot}\leq\mathrm{M}_{\rm disk}\leq 9\cdot 10^{-3}\ \mathrm{M}_{\odot}$ is found if the best fit vertical structure from (Kama et al. 2016b) is assumed. This estimate agrees with the combination of HD 1-0 and HD 2-1 line fluxes, which favour a gas mass of $7.7\cdot 10^{-3}\ \mathrm{M}_{\odot}$.

• •

Detailed modeling of the TW Hya disk shows that the difference between HD- and CO-based gas masses is mitigated by including CO isotopologue selective effects. A carbon underabundance of $\sim 10$ with respect to the ISM can also explain the CO isotopologue observations.

In the interest of future far-infrared observatories, the line-to-continuum ratios for HD 1-0 and HD 2-1 where calculated for different spectral resolving power. If a maximum signal-to-noise on the continuum of 300 is assumed, it was shown that future far-IR missions need a spectral resolving power $R\geq 300$ (equivalent to $\Delta v\leq 500$ km s*-1*) to detect HD 1-0 for all disk masses. To detect HD 2-1 towards all models, a spectral resolving power of $R\geq 1000$ (equivalent to $\Delta v\leq 150$ km s*-1*) is required. Furthermore, a $5\sigma$ sensitivity of $1.8\cdot 10^{-20}$ W m*-2* ($2.5\cdot 10^{-20}$ W m*-2*) is required to detect HD 1-0 (2-1) in disks with gas mass $\geq 3\cdot 10^{-5}\ \mathrm{M}_{\odot}\approx 0.03\ \mathrm{M}_{\rm jup}$. Both requirements can be met by proposed future missions such as SPICA.

Acknowledgements.

The authors would like to thank Inga Kamp and Ted Bergin for the useful discussions and for providing the details for the SPICA SAFARI and SOFIA HIRMES instruments. They also acknowledge Arthur Bosman for his support with the DALI code and Davide Fedele for providing the HD 2-1 line flux that was observed for TW Hya. Astrochemistry in Leiden is supported by the Netherlands Research School for Astronomy, by a Royal Netherlands Acadamy of Arts and Sciences (KNAW) professor prize, and by the European Union A-ERC grant 291141 CHEMPLAN.

Appendix A Abundance and temperature maps of TW Hya

In Figure 14 the abundance maps of 13CO and C18O are shown for the TW Hya model (M${}_{\rm disk}=2.3\cdot 10^{-2}\ \mathrm{M}_{\odot}$, $h_{c}=0.1$, $\psi=0.3$ and $\delta_{\rm C}=\delta_{\rm O}=0.01$.) The gas and dust temperature structures for the same model are shown in the panels of Figure 15.

Appendix B Abundance and emission maps of grid models

The HD abundance maps for several models with different vertical structures and disk masses are shown here. Figure 16 presents four models with different vertical structures: $(h_{c},\psi)=[0.05,0.1],[0.05,0.5],[0.3,0.1],[0.3,0.5]$. These represent the most extreme cases in the grid. Note that the top and bottom rows in the figure have different vertical extents. Also note that the HD 1-0 emission (blue solid contours) originates from approximately the same horizontal region in all four models, between a few AU and several tens of AU.

Figure 17 presents the abundance maps for three disks with different disk masses: M${}_{\rm disk}=[10^{-5},10^{-3},10^{-1}]\ \mathrm{M}_{\odot}$. Note that the height $\tau=1$ surface for the line optical depth (blue dashed line) increases with disk mass.

Appendix C Effects of including hydrostatic equilibrium

As shown in Figure 3, the HD emission originates from regions close to the disk surface. Accordingly, the HD emission will depend on the assumed disk vertical structure. In the models presented here the vertical structure is parametrized with a Gaussian (cf. Equation 2), to mimic a balance between gravitational force and pressure gradient, which in reality would depend on the temperature structure of the disk. Another possible approach is to the hydrostatic equilibrium to derive the vertical structure (e.g., Woitke et al. 2009). This results in a puffed-up inner rim, which shadows the outer disk from direct irradiation by the star. However, this may be a temporary feature that could be dispersed by e.g. stellar winds.

For completeness, we have investigated how the inclusion of hydrostatic equilibrium would affect the HD emission. The fiducial model described in Section 3.2 was rerun using eight iterations of the hydrostatic solver of DALI, where the temperature structure of the $(n-1)^{\rm th}$ model is used to calculate the vertical structure of the $n^{\rm th}$ model. The resulting vertical structure is then used to recalculate the dust and gas temperatures, the chemistry and the HD emissions.

C.1 The hydrostatic solver

One way to obtain the vertical structure of the disk is to solve the equations of hydrostatic equilibrium. The hydrostatic equations in cylindrical geometry read (e.g. Hartmann 2000)

[TABLE]

where the velocities in the $z$ and $r$ direction are assumed to be 0. Without self-gravity, the gravitational potential by the star is

[TABLE]

where $M_{*}$ is the mass of the star.

Assuming that the pressure gradient in the $r$ direction is negligible, the first hydrostatic equation yields the Keplerian velocity profile and the second hydrostatic equation can be solved independently for each radius. Instead of directly solving the equation for the density $\rho$, the equation can also be solved for pressure. The pressure and the density are related through $P=c_{s}^{2}\rho$, where $c_{s}$ is the sound speed which depends on the gas temperature $T_{\rm gas}$ and therefore varies with height $z$. The differential equation for $P$ reads

[TABLE]

which is solved by

[TABLE]

In the case of an isothermal disk, the above equation reduces to a Gaussian vertical density profile, as given in equation (2). The integration constant $C$ in the above equation results from the normalization conditions for the vertical column density and should thus adjusted such that $\int\rho(z)\mathrm{d}z$ is preserved.

C.2 Comparing with parametrized vertical structure

Figure 18 shows the effect of eight iterations of the hydrostatic solver on the density and gas temperature structure of the fiducial model (M${}_{\rm disk}=2.3\cdot 10^{-2}\ \mathrm{M}_{\odot}$, $\psi=0.3$ and $h_{c}=0.1$, cf. Section 3.2). It also shows the regions where the HD 1-0 emission is produced before and after applying the hydrostatic solver. The top panel shows that the HD 1-0 emission has moved slightly towards the midplane for the hydrostatic model. This is likely the result of the decrease in height of the line optical depth surface. The bottom panel shows that the gas temperature in the HD 1-0 emitting region has changed by less than 5%.

The effect of the hydrostatic solver on the HD emitting regions is shown in Figure 19. While disk models obtained with the hydrostatic calculations are vertically more extended, the HD emitting regions have not extremely changed. This is also reflected in the integrated fluxes, which for the hydrostatic model are $\sim 50\%$ (HD 1-0) and $\sim 66\%$ (HD 2-1) lower compared to the model with the parametrized vertical structure. While this affects the inferred disk masses by a similar absolute factor, all the trends presented here should be robust. Also note that varying the parametric vertical structure changes the HD fluxes by a similar or larger factor (cf. Figure 4). The mass uncertainties calculated in Section 4.1 therefore remain unchanged.

Appendix D Deuterium chemistry

Table 6 list the reactions that were included in the deuterium chemical network. The rate coefficients were adapted from Roberts & Millar (2000); Walmsley et al. (2004) and Glover & Jappsen (2007). Where the rate coefficients were not known, the coefficients for the analogue reactions with H instead of D were used.

Appendix E Line fluxes of the TW Hya model

The extended version of Figure 13 is shown here. Its setup is similar to Figure 6 in Kama et al. (2016b). The fluxes and upper limits for the observations used can be found in Table B.1 in Kama et al. (2016b). The line fluxes of 13CO 6-5, C18O 3-2 and C18O 6-5 were calculated from the data presented in Schwarz et al. (2016).