© ESO, 2013

1. Introduction

In the past decades, the abundance of deuterated molecules has been shown to be significantly enhanced in low-temperature environments (e.g. Snell & Wootten 1977; Wootten et al. 1982; Roueff et al. 2000; Loinard et al. 2001) compared to the D/H elemental ratio in the interstellar medium (~1.5 × 10^-5, Linsky 2003). The level of deuterium fractionation of molecules is a potentially good tracer of their formation process. This characteristic was demonstrated and used to infer chemical routes of formation for, e.g., methanol (Parise et al. 2004, 2006) and ammonia (Roueff et al. 2005).

Water is an important molecule, because it is the main constituent of grain mantles and an essential molecule in the oxygen chemistry in dense interstellar clouds. It can be formed by three mechanisms: dissociative recombination in gas phase by ion-molecule chemistry (Bergin et al. 2000), combination of atoms on the surface of cold dust grains at low temperature (Tielens & Hagen 1982; Cuppen et al. 2010), and gas phase reaction at high temperature (>250 K) (Wagner & Graff 1987; van der Tak et al. 2006). To constrain the formation process of water, studying the D/H ratio of water is a promising method, because the three scenarios would result in different fractionation.

The D/H ratio of water has been studied toward several hot-cores (Jacq et al. 1990; Pardo et al. 2001; van der Tak et al. 2006) and low-mass protostars (Parise et al. 2005a; Liu et al. 2011). The studies of the high-mass hot cores, Orion IRc2 and AFGL 2591, and the low-mass protostars, IRAS16293-2422 and NGC1333-IRAS2A, suggest that the HDO emission mainly originates in the envelope and with different HDO/H₂O ratios in the two different environments (~10^-2 for low-mass cases and  ~10^-4 for high-mass cases). In addition, the D/H ratios of water are significantly lower than the formaldehyde and methanol fractionation in the same sources (Parise et al. 2006; Turner 1990), which might be surprising if all species formed simultaneously on dust surfaces. This low deuterium enrichment of water, if confirmed, is a very valuable constraint for astrochemical models that strive to explain the chemical processes involved in the formation of water.

Previous studies have shown that the distribution of the HDO abundance can be constrained in low-mass protostellar envelopes with the radiative transfer analysis of several HDO lines, which span different energy conditions (Parise et al. 2005b; Liu et al. 2011).

While some detailed studies that also involve the ground transitions of HDO have recently been performed toward low-mass star-forming regions (Parise et al. 2005a; Liu et al. 2011; Coutens et al. 2012), most of the water fractionation studies in high-mass hot cores were performed more than ten years ago, when the submillimeter spectrum was still extremely difficult to observe. Thanks to the advance of sensitive submillimeter receivers, such as those on the APEX telescope, and to the launch of Herschel, it is now possible to carry out much more detailed radiative transfer analyses of HDO and H₂O. We aim here to carry out such a detailed study on one massive protostellar object, G34.26+0.15, to constrain the formation of water in a high-mass hot core.

G34.26+0.15 is a well-studied hot core associated with an ultra-compact HII region. The HDO/H₂O ratio has been derived to be 1.1 × 10^-4 toward this source from the analysis of one HDO and one H₂O transition line (Gensheimer et al. 1996). With the five transition lines including the ground-state lines in the submillimeter-wave range (241, 225, 266, 464, and 893 GHz with APEX) and with higher angular resolution observations (241 and 225 GHz with SMA), we can study the HDO abundance and distribution at the core scale. This source was one of the targets of the WISH Key Program (van Dishoeck et al. 2011) with the Herschel telescope, ensuring that information on the abundance of water is available (Wyrowski et al. 2010). This hot core has been the target of molecular line surveys and showed high levels of deuterium fractionation (Hatchell et al. 1998b; Jacq et al. 1990; Gensheimer et al. 1996). The D/H ratios of formaldehyde and HCN toward G34.26 are  ~1% and  ~0.1% (Roberts & Millar 2007; Hatchell et al. 1998a).

This paper is organized as follows. In Sect. 2 we present the observations. We present the observational results and the radiative transfer modeling in Sects. 3 and 4. The discussion and conclusion are given in Sects. 5 and 6.

2. Observations

2.1. Single-dish observations

HDO observations were carried out with the APEX telescope. We observed the 225, 241, 266, 464, and 893 GHz lines (Table 1) toward G34.26+0.15 at the hot core position α₂₀₀₀ = 18_h53^m18$\stackrel{\mathrm{s}}{.}$57 and δ₂₀₀₀ = 01^o14^m58$\stackrel{\mathrm{s}}{.}$3 in September 2010. The focus was checked on Saturn, and the local pointing on G34.26+0.15 itself. We used the wobbler-switching mode with a throw of 240′′. Table 1 lists the characteristics of the HDO observations. The temperature scale here was converted from ${\mathit{T}}_{\mathrm{A}}^{\mathrm{\ast }}$ to T_mb using the beam efficiencies indicated in Table 1, which we took from the APEX¹ website. The beam efficiency for CHAMP⁺ was measured in July 2010 on Mars.

2.2. Interferometer observations

The observations were carried out on April 26 and May 10, 2011 with the Submillimeter Array (SMA) on Mauna Kea, Hawaii. We used the 230 GHz receivers to observe the 225 and 241 GHz lines as well as the 1.3 mm continuum simultaneously. The primary-beam size (HPBW) of the 6-m-diameter antennas at 230 GHz was measured to be  ~54′′. The phase tracking center of this observation is the same as for the APEX observations. The spectral correlator covers 4 GHz bandwidth in each of the two sidebands separated by 10 GHz. The frequency coverages are 225.76−229.73 GHz in the lower sideband and 237.73−241.71 GHz in the upper sideband. Each band is divided into 48 chunks. A hybrid resolution mode was set with 512 channels per chunk (16.2 kHz resolution) for the HDO J_K − K+ = 3_1,2−2_2,1 (s47), and 256 channels per chunk (32.4 kHz) for the HDO J_K − K+ = 2_{1, 1}−2_1,2 (s46).

The visibility data were calibrated using the MIR software package, which was originally developed for the Owens Valley Radio Observatory. The absolute flux density scale was determined from observations of Neptune for all data. The pair of nearby compact radio sources 1743-038 and 1911-201 was used to calibrate relative amplitude and phase. We used 3C279 to calibrate the bandpass.

The calibrated visibility data were Fourier transformed and CLEANed using the MIRIAD package. The map was made using uniform weighting. The synthesized beam size was 3$\stackrel{\mathrm{″}}{.}$4   ×   2$\stackrel{\mathrm{″}}{.}$5 with a position angle of –64°. The continuum map was obtained by averaging all line-free channels in the two bands except chunks s44-s48 because they have different frequency resolution.

3. First results

3.1. Single-dish observations

The spectra observed toward G34.26+0.15 are presented in Fig. 1. All transitions are detected. They all appear in emission, except for the 893 GHz line, which shows a complex profile, dominated by an absorption component. This is because the 893 GHz line is a ground transition and has a high critical density. The line profiles of the three high-excitation transitions (225, 241, and 266 GHz) are close to Gaussian and their linewidths are  ~6.3−6.7 km s^-1, obtained from a Gaussian fit. Not only are the profiles of 225 and 241 GHz HDO lines very similar, but their intensities are very close to each other. The linewidth of the 464 GHz HDO line is slightly wider than the high-excitation lines (~7.2 km s^-1), although this does not appear to be significant, owing to the large error bars, and as shown in Fig. 1f, which displays the superposition of the 464 and 241 GHz line profiles. It is found that the difference between these two spectra is still within the 3σ noise level, and that there is no obvious evidence for the existence of a broad component in the 464 GHz line, as was found in some ground transition H₂O lines (Kristensen et al. 2010b; Chavarría et al. 2010).

Table 2 lists the properties of the observed HDO lines. Using the integrated line intensities, we performed a rotation diagram analysis assuming that the lines are optically thin and in local thermodynamic equilibrium (LTE) (Goldsmith & Langer 1999). We considered only the three high-excitation lines, because the two ground-state lines are suggested to be optically thick (see Sect. 4.3). We assumed a 1″ source size, as suggested from the interferometric observations (Table 3). Because this size is smaller than the three observational beams, the derived column density and source size are degenerate, and the temperature does not depend on the size assumption. The resulting rotation diagram is displayed in Fig. 2. The rotation temperature obtained is  ~$\mathrm{14}{\mathrm{1}}_{-62}^{\mathrm{+}\mathrm{489}}$ K, which is consistent with the estimates of the kinetic temperature (160  ±  30 K) from the analysis of OCS, CH₃OH, SO, and HC₃N (Mookerjea et al. 2007). The derived total molecular column density of the hot HDO is 4.1 × 10¹⁷ cm^-2 averaged on the assumed 1′′ extended source.

3.2. Interferometer observations

3.2.1. Continuum at λ = 1.3 mm

The total continuum flux retrieved within the SMA primary beam is 8.0 Jy. The continuum emission is found to be slightly resolved (Table 3). Figure 3 shows an overlay of the 1.3 mm continuum image with the free-free continuum emission observed by Avalos et al. (2006). The continuum peak detected at 1.3 mm does not coincide with HII regions A and B, but is roughly consistent with HII region C within the beam size (detail see Fig. 4). At 107 GHz, the integrated continuum intensity of 6.7 Jy is suggested to be dominated by free-free emission and only  ~10% can be contributed from dust (Mookerjea et al. 2007). To date, neither the turnover point of free-free emission nor the spectral index are clear. Thus, if the free-free emission is already optically thin at a wavelength longer than 3 mm, the free-free emission will contribute  ~6 Jy to the emission at 1.3 mm. In other words, the detected continuum flux with the SMA is the combination of free-free and dust emission. Figure 4 also supports this argument, because the peak of the 1.3 mm continuum is located between the peak of the free-free emission (from which it is offset by  ~0.3″) and the HDO peak. This point will be important when modeling the emission with RATRAN (which models only the dust continuum).

3.2.2. Spectral lines

The observational results for the two HDO lines are listed in Table 3. The integrated intensity maps of the two HDO lines (Fig. 4) show that their emission distributions are similar, suggesting that the emission comes from the same gas. In addition, the emission is not coincident with the A and B compact HII regions. The peak of the HDO emission is offset from the 1.3 mm dust continuum peak by approximately 0.3−0.4′′ (Table 3), under the theoretically positional error of approximately 0.04′′ with a signal-to-noise ratio (S/N) of 43. We believe that this offset is significant, because the continuum and the lines were observed simultaneously. This offset seems to be much smaller than the offset (2′′) observed for CH₃CN by Watt & Mundy (1999), and on the lower end of the different offsets observed for different molecules (CH₃OH, SO, and OCS) by Mookerjea et al. (2007). Those molecules have lower critical densities² compared to the two HDO lines, so excitation effects might partly explain the discrepancy between the distribution of the molecules.

Comparing the results observed with APEX, we find that the missing flux of the HDO 225 GHz and 241 GHz line emission are  ~37.1% and  ~36.4%, suggesting that most HDO emission is coming from the center of the molecular cloud core, which is consistent with the fact that the line shapes are very similar. Here, the SMA map was convolved with the beam size of APEX. The calibration uncertainties of APEX and SMA are usually considered to be  ~20% and 5−15%. It is difficult to believe that the apparent missing flux is real, because that would require that these high-excitation HDO lines show some extended emission, which is very unlikely in view of their critical density and energy. Indeed, the differences in fluxes between the SMA and APEX spectra can be also caused by the side lobes in the SMA images. Because the declination (δ₂₀₀₀) of our source is around 1 degree, the side lobes are quite strong and cannot be easily removed thoroughly. The negative side lobes will lower the obtained integrated flux (SMA) in this case. We therefore assume that the flux discrepancy is only originating in the uncertainty of the calibration and imaging of the data.

We resolved the CN 2–1 hyperfine structure of J = 3/2−1/2 and J = 5/2−3/2. All CN absorption lines show redshifted absorption; we discuss in more detail in Sect. 5.2.2.

4. Modeling

4.1. Physical profiles of the source

To analyze the HDO data, we used the 1 D Monte Carlo code RATRAN, developed by Hogerheijde & van der Tak (2000). The radiative transfer model takes as input the physical (n, T) profiles of the source. The heating source of the hot core and its density structure are not clearly understood toward G34.26+0.15, because there is no high angular resolution study (≤ 1″) in the submillimeter band. It has been argued that the hot core could be externally heated by this HII region (Mookerjea et al. 2007), because the hot core lies off the cometary-shaped HII region (Watt & Mundy 1999). This cometary-shape may be caused by the bow shock produced by the wind from a young star that moves at supersonic speed through an ambient molecular cloud (Wood & Churchwell 1989; van Buren et al. 1990; Watt & Mundy 1999) or by the champagne outflow (Andersson & Garay 1986; Fey et al. 1992; Carral & Welch 1992; Gaume et al. 1994). The separations between the HII region and the emission peaks of most N- and O-bearing molecules are  ≤1″ (Mookerjea et al. 2007, see Fig. 3). However, recent studies with HCN, HCO⁺, CO, and their isotopes observed with APEX show no evidence for this scenario and, moreover, the results show self-absorption and strong red asymmetries, which suggest an in-falling source with internal heating (Rolffs et al. 2011). In addition, the studies of the NH₃ (3₂−2₂) line observed with SOFIA also present a clear in-fall absorption in the spectra which are well reproduced also with Rolffs’ model (Wyrowski et al. 2012). This may be because they probe larger scales and the structures of the heating area are too small to be resolved in their observations (Rolffs et al. 2011). Therefore, the external heating structures could be internal in the model. Neither our APEX nor our SMA data have a sufficiently high resolution to probe the scale of the heating region; therefore, the two models we adopted here have the temperature structures derived by the inner heating source.

We use in the following the physical profiles of two different models (Nomura & Millar 2004; Rolffs et al. 2011, Fig. 5). Nomura’s model is constrained by the observed spectral energy distribution (SED) data at wavelengths from far-infrared to millimeter. Most data are from Chini et al. (1987) and were observed with the ESO 1 m telescope, except for an additional point at around 3 mm (Nomura & Millar 2004, see Fig. 1c). The 3-mm data point is taken from Watt & Mundy (1999) and were observed with BIMA (Nomura, private communication). As in Nomura’s study, the model with the column density of N = 10²⁵ cm^-2 is used here. The density and temperature profiles of the Rolffs model aim to reproduce the radial intensity profile from LABOCA data (345 GHz continuum). We adopted the opacity and velocity field setting of the Rolffs model for both models (opacity: ice mantle coagulation at a density of 10⁵ cm^-3, Ossenkopf & Henning 1994; velocity field: $\mathit{V}\mathrm{=}\mathrm{-}\mathrm{0.2}\mathrm{(}\frac{\mathrm{2}\mathit{G}{\mathit{M}}_{\mathrm{in}}}{\mathit{r}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}$, Rolffs et al. 2011). The turbulent linewidth, db, is set to be 3.0 km s^-1. The adopted distance is 3.7 kpc (Kuchar & Bania 1994). The HDO collisional rates used in this study are computed by Faure et al. (2012) for para-H₂ and ortho-H₂ separately. In the modeling, we assumed that the ortho-to-para ratio of H₂ is in LTE in each cell of the core.

4.2. Modeling procedure

The HDO emission was modeled with a jump model, where the fractional abundances of deuterated water, relative to H₂, in the inner part of the source (T    >    100 K, assumed evaporation temperature) and in the outer part (T    ≤    100 K) are two free parameters, X$\begin{array}{}\mathrm{HDO}\\ \mathrm{in}\end{array}$ and X$\begin{array}{}\mathrm{HDO}\\ \mathrm{out}\end{array}$. We made two sets of comparisons, one with the APEX data and the other with the SMA data. For comparison with the APEX data, the modeled maps were convolved with the APEX beam. For comparison with the SMA data, we used MIRIAD to create mock observations based on the modeled data. We first inserted the values of the observed parameters (epoch, rest frequency, and the coordinate of the observed reference position) into the blank header of the modeled data. Then, we multiplied the modeled image by the SMA primary beam. Afterward, we applied the (u, v) sampling retrieved from the SMA data to the modeled image and inverted the resulting (u, v) data to obtained the simulated maps/spectra. We then performed χ² analyses for X$\begin{array}{}\mathrm{HDO}\\ \mathrm{in}\end{array}$ ranging from 1  ×  10^-9 to 4  ×  10^-7 and for X$\begin{array}{}\mathrm{HDO}\\ \mathrm{out}\end{array}$ ranging from 1  ×  10^-15 to 2  ×  10^-11 for the Rolffs model and for X$\begin{array}{}\mathrm{HDO}\\ \mathrm{in}\end{array}$ ranging from 1  ×  10^-12 to 5  ×  10^-9 and for X$\begin{array}{}\mathrm{HDO}\\ \mathrm{out}\end{array}$ ranging from 1  ×  10^-17 to 5  ×  10^-10 for the Nomura model. We intend to also model the line velocity profiles and all spectra, including the spectra from the SMA data. Here the definition of χ² is $\mathrm{\Sigma }\frac{\mathrm{(}{\mathit{T}}_{\mathrm{mb,ob}}\mathrm{-}{\mathit{T}}_{\mathrm{mb,mod}}{\mathrm{)}}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}$, where the sum is over each channel in each spectrum. The σ within this χ² analysis includes the statistical errors and uncertainties in flux calibration (~10% for 225 and 241 GHz observed with SMA,  ~30% for 893 GHz and  ~20% for other lines observed with APEX), but does not include any uncertainty in the adopted collisional rate coefficients used in the excitation calculation (Faure et al. 2012).

4.3. Results

The black solid lines in Fig. 6 present the contours delimitating the 1σ (68.3%), 2σ (95.4%) and 3σ (99.7%) confidence intervals. These contours are derived with the method described in Lampton et al. (1976), who showed that the $\mathrm{\Delta }\mathit{S}\mathrm{=}{\mathit{\chi }}^{\mathrm{2}}\mathrm{-}{\mathit{\chi }}_{\min }^{\mathrm{2}}$ random variable follows a χ² distribution with p variables, p being the number of parameters. Here p = 2 (X_in and X_out), so that these contours correspond to χ² = ${\mathit{\chi }}_{\min }^{\mathrm{2}}$+2.3, ${\mathit{\chi }}_{\min }^{\mathrm{2}}$+6.17 and ${\mathit{\chi }}_{\min }^{\mathrm{2}}$+11.8. Here we used 1σ, 2σ, and 3σ in analogy with Gaussian-distributed random variables.

The fitting parameters and obtained results are listed in Table 4. The observed spectra are overlaid with the modeling results of the best fit in Fig. 7. Here the system velocity of 58 km s^-1 is assumed (V_LSR). The minimum ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}$ are not close to 1, which indicates that the spectra are not well reproduced by the two different physical models. This is also obvious in Fig. 7. When the Rolffs model is adopted (red lines in Fig. 7), the emission at 225, 241, and 266 GHz is not high enough at the large scale to fit the observed data (APEX data), but the most compact emission in the 225 and 241 GHz lines is well reproduced (SMA data). At the same time, the simulated spectra at 464 and 893 GHz are optically thick (τ = 2.1 and τ = 24.1 at lines and source center) and the line profile in 893 GHz spectrum does not fit the observed spectrum either. In Fig. 7, the blue and green lines indicate the results of the best-fit model using the Nomura physical profile, the green line was obtained excluding the spectrum at 893 GHz from the χ² analysis. In this case, the fractional abundances are 2 × 10^-9 and 5 × 10^-12 for the inner and outer region. This model fails to reproduce the central emission observed with the SMA.

Comparing the results reproduced by the two models (red and blue), we find that the best fit based on the Rolffs model is closer to the observed data (${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}\mathrm{=}\mathrm{6.0}$). The fitting problems in the ground transitions (464 and 893 GHz) are slightly less severe than for the Nomura model (see Fig.7). Moreover, from the comparison of the two fitting results of the Nomura model (blue and green), we find that the absorption in the 893 GHz spectra and emission in 464 GHz are too strong if we try to increase the emission in the high-excitation lines to fit the spectra (225, 241, 266 GHz). Indeed, the continuum flux reproduced with the Nomura model (~49 K) is much higher than the real flux (~13.0 K). This is also found in Appendices A.2. and A.3. Figure 8 shows the angular distribution of the 225 GHz line from the SMA observation and the two models. Here again the model based on the Rolffs profiles fits the data much better. This is similar for the 241 GHz line and is therefore not shown.

The model based on the Rolffs centrally heated core reproduces the SMA observed line distribution well. This does not exclude the possibility that the heating of the very central dust core is external, but it shows that even at the resolution of our SMA data, we do not see clear evidence of a non-spherically symmetric heating of the envelope, and that the large-scale properties of the envelope can therefore be derived using such a model. Distinguishing between internal and external heating at small scales would require observations with much higher angular resolution.

Based on the superior results obtained with the Rolffs profiles, we focus on the Rolffs model below and try to reproduce the lines with a two-jump model of the HDO abundance distribution. Finally, we modify Rolffs physical structure to try to fit the observed data. The results of the two-jump model and modified models are discussed in Sects. 5.1 and 5.2 and in Appendices A.2. and A.3.

5. Discussion

5.1. Two-jump model

To reproduce more emission in the high-excitation lines and decrease the emission in the 893 GHz line, we modeled the spectra with a two-jump model instead of a one-jump model. The additional jump was arbitrarily set to be located at T = 200 K. Figure 9 show the results of four testing models. The adopted HDO fractional abundances of the models are listed in Table 5. The red lines indicate the results of the original best-fit one-jump model, which underestimates the high-excitation lines. The two-jump models (2J-1 and 2J-3) are found to be somewhat better for these transitions (Fig. 9). This suggests that the two-jump model can produce higher emission in the high-excitation lines. This property has already been found for Sgr B2 (M) (Comito et al. 2010) and this is due to the very high critical density of these lines (see Table 1). Moreover, although we decreased the HDO fractional abundance in the second inner region of the core in model 2J-2 and 2J-4 by a factor of 2, the emission in the 893 GHz line is little affected (<1σ). In addition, the line shapes of different lines do not change from model to model. Therefore, the results suggest that adopting a two-jump model can help in increasing the emission in the high-excitation lines, but cannot improve the 893 and 464 GHz line shapes. We note, however, that because the inner part is very compact, the enhancement of the high-excitation lines in the APEX beam when increasing the HDO abundance in this inner part is linked to an even faster relative enhancement of the flux in the SMA beam. As a consequence, one cannot increase the inner HDO abundance too much, even if the effect on the APEX lines is almost unnoticeable.

5.2. Modified physical profiles

In the previous sections, we have seen that the observed spectra are not perfectly reproduced using the Rolffs physical profile, even with a two-jump model. This is also the case for some HCN lines (see Fig. 5 in Rolffs et al. 2011). We now investigate how the model can be improved by modifying the physical profiles from the Rolffs model. We first investigate where the HDO emission and absorption take place by means of a shell analysis and then test the effect of modfying the physical profiles of the Rolffs model.

5.2.1. Shell analysis and modified density/temperature profiles

To analyze the discrepancy between the best-fit model and the observed data, we performed a shell analysis by sequentially turning the abundance of HDO in the innermost shells to 0 and then simulating new spectra (Fig. A.1). The results of this study are shown in Appendix A.1. We find that the emission of all transitions (even the 893 GHz line) originates in the inner region, whereas the absorption originates in the outer region. This explains why it is impossible to improve the line profile of the 893 GHz within the frame of a jump model. Improvement can therefore only be achieved by modifying the physical profiles (density, temperature, velocity) of the source. In Appendices A.2 and A.3, we show that modifications of the density and temperature profiles are unable to improve the quality of the fit. In the next section, we investigate the effect of modifying the velocity field.

5.2.2. Modified velocity profile

– Modified infall velocity field So far, we used the velocity profile provided by Rolffs et al. (2011) for all modeling. This is a modified free-fall velocity field. However, the line shape in the 893 GHz spectra cannot be simply reproduced by our present model, which might imply that there is an additional velocity field in the core.

Figure 10 shows the moment-0 map overlaid with the dust continuum map (a) and the position-velocity diagram (b) of the CN 2–1 absorption line observed with the SMA simultaneously to the HDO lines. Clearly, the distribution of the CN absorption is consistent with that of the dust continuum, as in the case of HDO. Moreover, the central velocity of the position-velocity diagram of the CN absorption is about 61 km s^-1, not 58 km s^-1. This redshifted absorption has already been found for many molecules, for instance NH₃, CO, ¹³CO, and HCO⁺ (Heaton et al. 1989; Gómez et al. 2000; Wyrowski et al. 2012). Thus, it is possible that the gas that produces the absorption feature in the HDO 893 GHz spectra has a similar velocity distribution. Therefore we tried to modify the velocity field by adding a 3 km s^-1 difference from the outer region to the inner region of the core (Fig. 11). The results are shown in Fig. 12. The spectra are simulated with the fixed fractional abundances 6 × 10^-8 (X$\begin{array}{}\mathrm{HDO}\\ \mathrm{in}\end{array}$) and 5 × 10^-12 (X$\begin{array}{}\mathrm{HDO}\\ \mathrm{out}\end{array}$) (best-fit model). The additional velocity (3 km s^-1) widens the emission of the optically thin lines (225, 241, and 266 GHz) somewhat, which enables a better fit of the line-width of the SMA lines. The opacities of the HDO 464 and 893 GHz lines are reduced and the 464 GHz line becomes optically thinner in models V3 and V4 (τ = 1.3 and 0.8). Moreover, the peak positions of the absorption in the simulated 893 GHz spectra gradually approach those of the observed 893 GHz spectra. Therefore, carefully modifying the velocity field might be a good way to improve the fits. However, the detailed velocity profile around the boundary region of the core cannot be constrained from our data now. To obtain better modeling results, it is necessary to observing the HDO line map at 893 GHz with high angular resolution. In addition, this result also suggests that studying the HDO spectrum at 893 GHz can help us to understand the velocity field of the core.

– Modified turbulence velocity width Because the line profiles are also affected by gas motions, the turbulence is also responsible for the line emission/absorption widths. The assumption that turbulence varies over the core has already been used for analyzing the spectra lines toward several massive sources and brought nicely fitting results (Herpin et al. 2009). Our previous study shows that the width of all simulated 893 spectra is obviously much wider than that of the observed spectrum, while the widths of the high-excitation lines are well fitted. Thus, we modified the turbulent velocity field by changing the parameter db in the model. In the Fig. 13a, the distribution of the populations of the upper level of high-excitation lines (3_1,2, 2_2,0, and 2_2,1) is different from those of the ground-transition lines (1_1,1 and 1_0,1). The peaks of the high-excitation line populations are mainly located in the layer 15−16 (1.4    ×    10^-2 − 1.7    ×    10^-2 pc) and those of the ground-transition lines are placed in the layer 21−24 (6.4    ×    10^-2 − 1.0    ×    10^-1 pc). We used these considerations for elaborating our new db settings with the aim to narrow down the width of the simulated 893 spectra without significantly affecting the width of the high-excitation lines. Figure 13b illustrates our different db settings. The HDO fractional abundances in the model db1 – db4 are fixed with the values of the best-fit model: X$\begin{array}{}\mathrm{HDO}\\ \mathrm{in}\end{array}\mathrm{=}\mathrm{6}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathrm{10}}^{-8}$ and X$\begin{array}{}\mathrm{HDO}\\ \mathrm{out}\end{array}\mathrm{=}\mathrm{5}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathrm{10}}^{-12}$.

Comparing the spectra at 893 GHz with those of the previous best-fit model (Fig. 14), we find that the line profile has significantly improved much in db2 – db4. The width of the absorption part is narrower than before and more similar to the observed one. Moreover, we have better fits compared to previous spectra at 266 GHz and the opacity is reduced in the 464 GHz spectra. In addition, the ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}$ model db2 is 5.5, which is better than that of the Rolffs original model (6.0), showing that the fits can be improved by modifying the turbulent velocity field. We also performed χ² analyses for the same range of X$\begin{array}{}\mathrm{HDO}\\ \mathrm{in}\end{array}$ and X$\begin{array}{}\mathrm{HDO}\\ \mathrm{out}\end{array}$ in Sect. 4.2 based on model db2. The fractional abundances of the new best-fit in this test are still the same as in our previous results (X$\begin{array}{}\mathrm{HDO}\\ \mathrm{in}\end{array}\mathrm{=}\mathrm{6}\mathrm{\times }{\mathrm{10}}^{-8}$ and X$\begin{array}{}\mathrm{HDO}\\ \mathrm{out}\end{array}\mathrm{=}\mathrm{5}\mathrm{\times }{\mathrm{10}}^{-12}$). This implies that the turbulence velocity is not constant throughout the envelope of G34.26+0.15. However, the widths at 225 and 241 GHz are narrower than the previous and observed ones, suggesting that the db in the inner region of the core is not wide enough in the models. Here again, high angular resolution data at 893 GHz would be very helpful.

– Composition of the turbulent and infall velocity field effects Because modifying the infall velocity field can help to move the peak of the absorption at 893 GHz (model V1), and modifying the turbulence velocity field can narrow this absorption, we ran models with modified velocity (V1 and db2). The result is shown in Fig. 15. Obviously, we can now reproduce part of the emission and all the absorption with better width and proper position, which allows us to study the fractional abundance in the outer region of the core. On the other hand, as expected, the emission at 893 GHz is not reproduced properly, because it originates from the inner region of the core (T > 100 K) and the model V1 only changes the velocity field in the outer region of the core (T < 100 K). Additional constraints in this direction will hopefully come from interferometeric data.

In conclusion, we have shown that it is possible to improve the original best-fit model by modifying the velocity field of the source. This improvement is obtained while keeping the fractional abundance of HDO to the value we derived with the original velocity setup of Rolffs et al. (2011). This study confirms the validity of the fractional abundance of HDO derived from the best-fit model in Sect. 4.3. Further steps toward constraining the velocity field and therefore a better HDO model will need high angular resolution observations of, e.g., the HDO 893 GHz line.

5.3. HDO/H₂O ratio

${\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}$O is always used to trace H₂O, because H₂O is easily optically thick. Here we compare the emission and absorption of the two para-${\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}$O lines 3₁₃ − 2₂₀ and 1₁₁ − 0₀₀ (Gensheimer et al. 1996; Wyrowski et al. 2010) with our HDO results. Because the analysis of the HDO lines suggests that high-excitation lines can trace and constrain the inner fractional abundance quite well (Parise et al. 2005b; Liu et al. 2011), we can also analyze the high-excitation ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}$O 3₁₃ − 2₂₀ line (E_u = 204.7 K) (Gensheimer et al. 1996) to study the ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}$O fractional abundance in the inner region of the core. We reproduced the reported integrated flux and line width (20.5 K and 8.4 km s^-1 listed in Table 2 in Gensheimer et al. 1996) of this ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}$O line with RATRAN using a one-jump model based on the Rolffs original profiles. The resulting fractional abundance of the para-${\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}$O in the inner region of the core is (1.0 ± 0.7) × 10^-7 (3σ). Thus, the HDO/H₂O is 3.0 × 10^-4 with the assumption that the standard ortho/para ratio is 3 and the ¹⁸O/¹⁶O ratio is 500. This ratio is consistent with the estimate of 1.1 × 10^-4 with an uncertainty of a factor of 4 (Gensheimer et al. 1996). Here, the adopted db for this high-excitation ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}$O line is 4.5 km s^-1 instead of 3 km s^-1, because the width of this ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}$O line is wider than that of the high-excitation HDO lines (8.4 km s^-1 for the ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}$O line compared to  ~6−7 km s^-1 for the HDO lines). The different properties of the velocity field of the high-excitation ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}$O line may indicate that some part of the emission does not come from the quiescent envelope and the high-excitation HDO lines may not probe exactly the same gas. Therefore, the derived D/H ratio of water in the inner region here should be a lower limit.

Figure 16 displays the two spectra of the ground-state 1_1,1 − 0_0,0 transition of HDO and ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}$O. The line shapes are very similar, suggesting they may originate from the same gas. Under the two assumptions that the excitation temperature of the line is negligible with respect to the temperature of the background continuum source and that only the ground state of the molecule is populated, the total column density of the absorbing HDO/${\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}$O can then be computed from (Comito et al. 2003; Tarchi et al. 2004) ${\mathit{N}}_{\mathrm{tot}}\mathrm{=}\frac{\mathrm{8}\mathit{\pi }{\mathit{\nu }}^{\mathrm{3}}}{{\mathit{A}}_{\mathit{ul}}{\mathit{c}}^{\mathrm{3}}}\frac{{\mathit{g}}_{\mathit{l}}}{{\mathit{g}}_{\mathit{u}}}\mathit{\tau }\mathrm{\Delta }\mathit{v,}$(1)where Δv is the line width and g_l and g_u are the statistical weights of the lower and upper level. The optical depth of the absorbed line can be determined by $\mathit{\tau }\mathrm{=}\mathrm{-}\ln \left(\mathrm{1}\mathrm{-}\frac{{\mathit{T}}_{\mathrm{L}}}{{\mathit{T}}_{\mathrm{C}}}\right)\mathit{,}$(2)where T_L is the brightness temperature of the line and T_C is the brightness temperature of the continuum. The line parameters from the Gaussian fit and the estimated values are listed in Table 6. The estimated total column densities are 1.8 × 10¹² cm^-2 for HDO and 1.8 × 10¹² cm^-2 for para-${\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}$O. Because the profiles of these two lines are very similar to each other, we assumed that both of them are originating from the same gas. Therefore, the [HDO]/[H₂O] ratio is  ~4.9 × 10^-4, assuming that the standard ortho/para ratio is 3 and the ¹⁸O/¹⁶O ratio is 500. However, the abundance of para-${\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}$O might be underestimated here, because the emitted area of the continuum at 1101 GHz might be larger than that of the ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}$O due to the large Herschel beam. Assuming that the dust emissivity spectral index β is 2 as for the interstellar medium (Draine & Lee 1984; Schnee et al. 2010), the derived continuum flux at 1101 GHz of the core seen with APEX will be 2.8 K in the Herschel beam and the estimated total column ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}$O density becomes 4.7 × 10¹² cm^-2. Then the obtained D/H ratio of water is  ~1.9 × 10^-4, which is close to what we obtained for the inner region of the core. We note that this value is uncertain, as the absorption takes place on top of emission in both cases, which we have ignored here.

The comparison of the HDO fractional abundances between different sources is shown in Table 7. The water deuterium fractionation toward the two low-mass protostars, NGC1333 IRAS 2A and IRAS 16293-2422, have been studied with several HDO transitions (Liu et al. 2011; Parise et al. 2005a; Coutens et al. 2012). Most transitions were also used in this work.

The abundance of HDO jumps by more than four orders of magnitude at the 3σ confidence level in G34.26+0.15. Similarly to this work, jump models are widely used to analyze the HDO and H₂O spectra of other high- and low-mass sources in Table 7, except for Orion KL. Ices are therefore also shown to evaporate from the grains in the inner part of the envelope in protostars of different masses. However, the D/H ratios of water in G34.26+0.15 is much lower than in low-mass protostars and similar to that in AFGL 2591 and other high-mass sources, such as Ori-IRc2, G31.41+0.31, and G10.47+0.03A (van der Tak et al. 2006; Gensheimer et al. 1996). In addition, the D/H ratios of formaldehyde and HCN in G34.26+0.15 are also lower than those in low-mass protostars. This can be explained if the dense and cold pre-collapse phase only lasts a short time and less CO freezes onto grain mantles. This hypothesis has been confirmed in the high-mass source AFGL 2591 (van der Tak et al. 2006).

Figure 17 shows the diagram of HDO/H₂O vs. HDCO/H₂CO and the values are taken from Table 7. It is found that the fractionation of water is at least  ~2.4 times lower than that of formaldehyde. This under-fractionation of water is found to be the most extreme for G34.26+0.15. One possible explanation is that if two molecules are mainly formed on the grain surface, the different deuterium fractionation can be due to the different pathways of each ice formation in translucent clouds: using a very simple chemical network on the grains, Cazaux et al. (2011) showed that the deuteration of formaldehyde is sensitive to the gas D/H ratio during the cloud collapse, while the deuteration of water depends on the dust temperature when ices form. More realistic models including water formation routes involving O₂ and O₃ also predict a higher fractionation of formaldehyde than water (Du et al. in prep.). However, recent studies showed that formaldehyde can also be formed in the gas phase from photodissociation reactions of grain species (Roueff et al. 2006). Formaldehyde can also efficiently be fractionated in the gas phase up to warm temperatures (~70 K), through the CH₂D⁺ route (Turner 2001; Roueff et al. 2007; Parise et al. 2009). It is therefore likely that the HDCO/H₂CO ratio observed in the gas is not directly the result of grain chemistry. A better comparison would be with methanol fractionation, which is generally thought to form essentially on the grains.

Confirming the observed trend of under-fractionation of water would require observationally characterizing the deuterium fractionation (for both water and methanol) toward a bigger sample of sources, spanning a large mass interval. This under-deuteration of water with respect to other species is an important feature, which should be reproduced by astrochemical models aiming at explaining deuterium fractionation.

6. Conclusion

We presented five HDO transitions observed with APEX and two HDO high-excitation transitions observed with SMA toward the high-mass protostar, G34.26+0.15. With the 1D radiative transfer code RATRAN, and the physical profiles obtained from Rolffs et al. (2011), we derived the HDO fractional abundances relative to H₂ in the inner and outer region of the core to be X$\begin{array}{}\mathrm{HDO}\\ \mathrm{in}\end{array}\mathrm{=}\mathrm{6}\mathrm{\times }{\mathrm{10}}^{-8}$ (T    >   100 K) and X$\begin{array}{}\mathrm{HDO}\\ \mathrm{out}\end{array}\mathrm{=}\mathrm{5}\mathrm{\times }{\mathrm{10}}^{-12}$ (T    ≤    100 K). This result shows that the HDO abundance is enriched in the inner region because of the sublimation of the ice in the same way as for other studied low- and high-mass sources, such as NGC1333 IRAS 2A, IRAS 16293-2422, and AFGL 2591 (Liu et al. 2011; Parise et al. 2005a; Coutens et al. 2012; van der Tak et al. 2006).

Although the abundance can be constrained, the HDO line profiles are not well reproduced, especially for the ground-state lines. We found that the fitted emission can be improved with two-jump models, suggesting the enhancement of the fractional abundance in the hot core. Modifying the physical profile, we showed that the velocity field is responsible for the HDO line profile at 893 GHz. Higher angular resolution observations of the HDO 893 GHz line are needed to constrain the velocity profiles of deuterated water in G34.26+0.15.

The H₂O abundance is estimated from one high-excitation and one ground transition para-H¹⁸O line (Gensheimer et al. 1996; Wyrowski et al. 2010). The D/H ratios of water are 3.0    ×    10^-4 in the inner region and (1.9−4.9)    ×    10^-4 in the outer region of the core. The deduced HDO/H₂O ratios in G34.26+0.15 are similar to those in other high-mass sources, such as AFGL 2591 (van der Tak et al. 2006), and are much lower than in low-mass star-forming regions, suggesting the possibility that the dense and cold pre-collapse phase is short for high-mass stars. Water is also found to be less fractionated than other molecules such as formaldehyde and methanol. This confirms this characteristic, which was observed previously in other sources, and provides a challenge for chemical models of deuterium fractionation on grain surfaces.

Online material

Appendix A: Model analysis

A.1. Shell analysis

To analyze the discrepancy between the best-fit model based on the Rolffs physical profiles and the observed data for the two ground transitions (464 and 893 GHz), we first constrained the physical shells responsible for the emission and absorption features. To this end, we simulated new spectra by sequentially turning the abundance of HDO in the innermost shells to 0 (Fig. A.1). The jump around 0.02 pc is at the shell whose temperature is about 100 K (jump model) The green dashed line represents the model where the HDO abundance is zero in all inner regions of the core (model SH6). Figure A.2 presents the synthetic spectra obtained with these different models (SH1 – SH6). Obviously there is almost no emission in any spectrum with model SH6, suggesting that all emission, including the emission in the 893 GHz spectrum, is produced in the inner region of the core. Therefore, although loweringn the abundance in the inner region can help to reduce the emission in the 893 GHz line, it will also decrease the emission in other high-excitation lines (225, 241, and 266 GHz). On the other hand, the absorption in the 893 GHz spectrum is entirely produced in the outer region of the core and its line-width is wider than the observed one, implying that the assumed turbulent linewidth is wider than the observed one. Thus, we cannot adequately reproduce all spectra with the Rolffs original physical profiles.

A.2. Modified density profile

From previous test, it was found that the emission in the 893 GHz spectrum cannot be managed by the Rolffs physical profiles, but only by changing the HDO abundance. We therefore modified in the first place the density profile of the Rolffs model. Figure A.3 compares the density profiles of three different modified models and original the Rolffs model (red line). Here we modified the density power law index (p), keeping the density at 0.1 pc fixed. The indices are 1.1 for model D1 (green), 1.7 for model D2 (dark-blue), and 2 for model D3 (purple). Modifying the density profile has a direct impact on the distribution of the continuum emission. The comparison of the radial profiles of the continuum emission of different models at multi-wavelength bands (225, 241, 335, and 848 GHz) is shown in Fig. A.4. Ignoring the misfitting parts within the beam (orange) and considering that free-free emission most likely also significantly contributes to the continuum emission, model D1 and the Rolffs model better reproduce the continuum emission in the 225 and 241 GHz bands than other models. For the 335 and 848 GHz band, the Rolffs model and all modified models are superior to the Nomura model. The best-fit results of each model are shown in Fig. A.5 (model D1: 1.0 × 10^-7(X$\begin{array}{}\mathrm{HDO}\\ \mathrm{in}\end{array}\mathrm{)}$ and $\mathrm{1.0}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathrm{10}}^{-11}\mathrm{\left(}\mathit{X}\begin{array}{}\mathrm{HDO}\\ \mathrm{out}\end{array}\mathrm{\right)}$; model D2: $\mathrm{5.0}\mathrm{\times }{\mathrm{10}}^{-8}\mathrm{\left(}\mathit{X}\begin{array}{}\mathrm{HDO}\\ \mathrm{in}\end{array}\mathrm{\right)}$ and $\mathrm{5.0}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathrm{10}}^{-12}\mathrm{\left(}\mathit{X}\begin{array}{}\mathrm{HDO}\\ \mathrm{out}\end{array}\mathrm{\right)}$; model D3: $\mathrm{3.0}\mathrm{\times }{\mathrm{10}}^{-8}\mathrm{\left(}\mathit{X}\begin{array}{}\mathrm{HDO}\\ \mathrm{in}\end{array}\mathrm{\right)}$ and $\mathrm{5.0}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathrm{10}}^{-14}\mathrm{\left(}\mathit{X}\begin{array}{}\mathrm{HDO}\\ \mathrm{out}\end{array}\mathrm{\right)}$). Comparing all spectra, we find that model D3 produced intensity similar to the Rolffs model in high-excitation lines and less emission in the 893 GHz line. However, the 464 and 893 GHz lines are optically thicker (τ = 2.9 at 464 GHz and τ = 29.6 at 893 GHz) than the lines produced by the Rolffs model. The high opacities result in self-absorption in the 464 and 893 GHz spectra (less emission). Moreover, results (464 and 893 GHz spectra) simulated with models D1 and D2 are very similar to the results produced with the original Rolffs model. Thus, we find that these models are not significantly better than the original the Rolffs model. In other words, the fits cannot be improved by simply modifying the density profile.

A.3. Modified temperature profile

To decrease the emission in 893 GHz spectra and maintain the emission in high-excitation lines (225, 241, 266 GHz), modifying the temperature profile is a possible way. We increased the temperature in the very inner region (models T1, T3, and T4) of the core and decreased the fractional abundance in the inner region of the core. The modified and original temperature profiles are shown in Fig. A.6. Model T2 contrasts with other models.

Like for changing the density profile, modifying the temperature profile also influences the continuum emission. Figure A.7 compares the radial profiles of the continuum emission of different models in the four bands. Again, we find that the Rolffs model and all modified models are superior to the Nomura model in all bands. For the 225 and 241 GHz band, model T2 and the Rolffs model fit better, because the dust continuum is less than 25% of the total continuum flux, implying that increasing the temperature in the central region of the core is a poor solution. In the 335 and 848 GHz bands, T1, T2, and the Rolffs models are superior to others. The indistinguishable continuum at 848 GHz is due to a combination of optical depth and beam dilution effect: the temperature profile is changed only in the very central region of the core where the opacity is high and the beam size is much larger than the scale of the modified region (beam is  ~7.4″ at 848 GHz and 2 × 10^-2 pc corresponds to  ~1.1″). For 241 and 225 GHz, the angular resolution is sufficiently high and the opacity sufficiently low to detect the very inner region of the core.

To compare the effect on the spectra, Fig. A.8 shows the fitting results of the different modified models with a fixed HDO fractional abundance in the inner and outer region (X$\begin{array}{}\mathrm{HDO}\\ \mathrm{in}\end{array}\mathrm{=}\mathrm{6}\mathrm{\times }{\mathrm{10}}^{-8}$ and X$\begin{array}{}\mathrm{HDO}\\ \mathrm{out}\end{array}\mathrm{=}\mathrm{5}\mathrm{\times }{\mathrm{10}}^{-12}$). The fitting spectra of the T1 and T4 models are similar to spectra of the Rolffs model in the high-excitation spectra at 225 and 241 GHz, while the T2 model produces less emission in all spectra. However, models T1, T3, and T4 produce much poorer fits to the two ground-transition spectra, and T2 already becomes optically thick (τ = 3.4) at 464 GHz. Hence, simply modifying the temperature profiles does not improve the agreement between the model and the observations.

---

Acknowledgments

The authors warmly thank the anonymous referee and M. Walmsley for very constructive comments and suggestions, which greatly improved this paper. We are grateful to R. Rolffs for many fruitful discussions as well as to F. Herpin for some interesting ideas. We also would like to thank M. Hogerheijde and F. van der Tak for their help in dealing with RATRAN. The authors are grateful to the WISH team for providing access to the H₂O data. F.-C. Liu and B. Parise are funded by the German Deutsche Forschungsgemeinschaft, DFG Emmy Noether project number PA1692/1-1.

References

All Tables

All Figures